#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <sstream>
#include <cmath>
#include <TRandom.h>
#include "TTree.h"
#include "TColor.h"
#include "TStyle.h"
#include <TFile.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TH3D.h> 
#include <TDirectory.h> 
#include <TSystem.h>
//#include "TPalette.h"
#include "TPaletteAxis.h"
#include <TRandom3.h>
#include <TString.h>
using namespace std;
#define HORIZSHIFT 0.4 
#define L 18374.0 //cm
// Routine to add the diffusion to FLUKA simulations
// Defining structure to hold data for the TTree

/*struct data {
  Int_t idelevo;
  Double_t neuteneo;
  Double_t bin_xo;
  Double_t bin_yo;
  Int_t num_elo;
//  Double_t zo;
  //Int_t index;
//  Double_t chargeo;
};
data data_tree;
*/


void run()
{
 TH1::AddDirectory(kFALSE);  // <- This avoids the messages "memory leak" when plotting more than one histogram
 TH2::AddDirectory(kFALSE);  // <- This avoids the messages "memory leak" when plotting more than one histogram
TTree::SetMaxTreeSize(500000);

// Output file
 TFile*fout_YES= new TFile("PROVAResults_YES.root","recreate");
 TTree*tout_YES = new TTree("tout_YES","TTree containing data in arrays");
 Double_t values_YES[6] = {0,0,0,0,0,0};
 TBranch *b2 = tout_YES->Branch("b2",values_YES, "a/D:b:c:d:e:f");

 TFile*fout_NO= new TFile("PROVAResults_NO.root","recreate");
 TTree*tout_NO = new TTree("tout_NO","TTree containing data in arrays");
 Double_t values_NO[6] = {0,0,0,0,0,0};
 TBranch *b3 = tout_NO->Branch("b3",values_NO, "a/D:b:c:d:e:f");

// Input file
 //TFile*f; 
 int nbins=1000;
 int energia_neut=5.e-09;
//*********** Some parameters and variables ******************
 double GeV_to_MeV=1000000000.;// 1e+9;
 double ioniz_potent=26.14; //w_Ar=26 ; w_CF4=28; w_C4H10=23
 double D=253.467 ;//MICRONS/CENTIMETER**0.5  +-    3.54%
 double D_CM=253.467e-4;//CM/CM**0.5
 // VARIABLES TO READ FORTRAN DATA
 int nb_story;
 int parID;
 int nn;
 float x0, y0, z0, x1, y1, z1;
 float endep;
 FILE *infile;
 char Formatted[200];
 float start_x;
 float start_y;
 float start_z;
 float end_x=0;
 float end_y=0;
 int el_nb=0;


 
//========================== LOOP ON OUTPUT FILES ========================
 for(int j=1;j<=1;j++) {// Loop over j
 
  if(j<=9) sprintf(Formatted,TString::Format("PXMG00%d_TRACK.txt",j));
  if(j>9) sprintf(Formatted,TString::Format("PXMG0%d_TRACK.txt",j));

  cout << Formatted << endl;
  infile = fopen (Formatted,"r");
  if (infile == NULL) {
                       printf("Unable to open file \n");
                      }
  int ncase = -1;
  int cont_stories=0;
  int count =0;

  while(1)  {
    fscanf(infile, "%d %d %d %f %f %f %f %f %f %f",&nb_story,&parID, &nn, &x0, &y0, &z0, &x1, &y1, &z1, &endep);
    if(nn!=1) cout << "Error at " << nb_story << "n = " << nn<<endl;
    if(nb_story!=ncase) {ncase=nb_story;
                        cont_stories++;}
    count++;
    if (feof(infile)!=NULL) break;}

   cout << " The file " << Formatted << " has "<< count<< " lines "<< endl;
   cout << " You have info about " << cont_stories << " events "<<endl;


   fseek(infile, 0, SEEK_SET); // Rewind file to beginning
   int *Anb = new int[count];
   int *AparID = new int[count];   
   int *Ann = new int[count];
   float *Ax0 = new float[count];   
   float *Ay0 = new float[count];
   float *Az0 = new float[count];   
   float *Ax1 = new float[count];
   float *Ay1 = new float[count];   
   float *Az1 = new float[count];
   float *Aendep = new float[count]; // In GeV  

   for(int i=0; i<count; i++) {
      fscanf(infile, "%d %d %d %f %f %f %f %f %f %f",&nb_story,&parID, &nn, &x0, &y0, &z0, &x1, &y1, &z1, &endep);
      Anb[i]=nb_story;
      AparID[i]= parID; 
      Ann[i]=nn;
      Ax0[i]=x0;
      Ay0[i]=y0;
      Az0[i]=z0;
      Ax1[i]=x1;
      Ay1[i]=y1;
      Az1[i]=z1;
      Aendep[i]=endep;
      if (feof(infile)!=NULL) break; }


  ncase = -1;
  int here=0;
  float en_neut=0;
  Int_t eve=0;

  TH2F*h_YES= new TH2F("h_YES","Histo d'appoggio ",nbins,-2.5,2.5,nbins,-2.5,2.5);
  h_YES->SetDirectory(0);
  TH2F*h_NO= new TH2F("h_NO","Histo d'appoggio ",nbins,-2.5,2.5,nbins,-2.5,2.5);
  h_NO->SetDirectory(0);
  for(int ii=1; ii<=nbins;ii++) {
   for(int kk=1; kk<=nbins;kk++) {
     h_YES->SetBinContent(ii,kk,0);
     h_NO->SetBinContent(ii,kk,0);
                                }
                               }


  for(int i=0; i<count; i++){//Here
     if (Anb[i]!=ncase) {
         ncase=Anb[i];
         for(int ii=1; ii<=nbins;ii++) {
           for(int kk=1; kk<=nbins;kk++) { 
             if(h_YES->GetBinContent(ii,kk)!=0) {
                 fout_YES->cd();
                 values_YES[0]=j;                  
                 values_YES[1]=energia_neut;//5e-09;//energy;
                 values_YES[2]=ncase;
                 values_YES[3]=ii;
                 values_YES[4]=kk;
                 values_YES[5]=h_YES->GetBinContent(ii,kk);
                 tout_YES->Fill();
                 h_YES->SetBinContent(ii,kk,0);  }

              if(h_NO->GetBinContent(ii,kk)!=0){
                fout_NO->cd();
                values_NO[0]=j;                  
                values_NO[1]=energia_neut;//energy;
                values_NO[2]=ncase;
                values_NO[3]=ii;
                values_NO[4]=kk;
                values_NO[5]=h_NO->GetBinContent(ii,kk);
                tout_NO->Fill();
                h_NO->SetBinContent(ii,kk,0);  }
                                         }
                                      }
                   }

    if (Anb[i]==ncase) {
       start_x= (Ax0[i]+Ax1[i])/2.;
       start_y= (Ay0[i]+Ay1[i])/2.;
       start_z= (Az0[i]+Az1[i])/2.;
       el_nb=Aendep[i]*GeV_to_MeV/ioniz_potent;
       h_NO->Fill(start_x,start_y, el_nb);
       for(int ii=0;ii<el_nb;ii++) {
            end_x = gRandom->Gaus(start_x,D_CM*sqrt(0.4-start_z));
            end_y = gRandom->Gaus(start_y,D_CM*sqrt(0.4-start_z));
            h_YES->Fill(end_x,end_y);}
                     }
//..........
        } //Here
   delete Anb;
   delete AparID;   
   delete Ann;
   delete Ax0;   
   delete Ay0;
   delete Az0;   
   delete Ax1;
   delete Ay1;   
   delete Az1;
   delete Aendep;
     } // Loop over j

#if 0
  fout_YES->cd();//tout->Write();
  tout_YES->Write();
  tout_YES->Delete();
  fout_YES->Close();
  fout_NO->cd();//tout->Write();
  tout_NO->Write();
  tout_NO->Delete();
  fout_NO->Close();

#endif

  fout_YES = tout_YES->GetCurrentFile();
  fout_YES->Write();
  delete fout_YES;
  fout_NO = tout_NO->GetCurrentFile();
  fout_NO->Write();
  delete fout_NO;
/*fout_YES->cd();//tout->Write();
fout_YES->Write();
fout_NO->cd();//tout->Write();
fout_NO->Write();*/

}