{
/***************************************************************************************************************/
/*                               Description of the program                                                    */
/*                                                                                                             */
/*     This program look for events with a good electorn in DIS (Q2>2, W>2) as the first particle and having   */
/* a positive pion as well.                                                                                    */
/* It also takes the leading pion of each events and plots the transverse momentum squared of these particles. */
/* It seperates the results for the solid target and the deuterium one.                                        */
/* It also give a plot of the total energy of the good event, and the perpendicular momentum as well as the    */
/* variation of number of protons with the transverse momentum of the positive pions.                          */
/*                                                                                                             */
/*                                                                                                             */
/***************************************************************************************************************/

gROOT->Reset();
gSystem->Load("/home/silvestr/eg2/porc1/slib/particle_identification.so");
ifstream in;
char x[200];
Int_t i=0,l=0,maxl[10],final;
Int_t numberEventsConElec=0;
Double_t visibleEnergyEvent;
Double_t pperpx,pperpy, pperp;
Int_t nentries=0; 
TDSTReader *R1 =new TDSTReader();
part_ident partident;
categorize p_categorize;
part_ident_pt p_pt;
TString presult;
TString *presult_ptr;
TFile *plots = new TFile("Q_sup_1.root","RECREATE");
Int_t maxl[10];
Int_t m ;

// histograms for electrons
//TH2F *betta_momentum_electrons = new TH2F("betta_momentum_electrons","betta vs. momentum for electrons",300,0,6,110,0.,1.1);
TH1F *NPHE_electrons = new TH1F("NPHE_electrons","NPHE for electrons:   ",300,0,300);
TH1F *z_coordinate_electrons = new TH1F("z_coordinate_electrons","z coordinate for electrons",300,-50,25);
TH1F *Q2_electrons = new TH1F("Q2_electrons","Q2 for electrons:   ",300,0.,5.);
TH1F *W_electrons = new TH1F("W_electrons","W for electrons:   ",300,0.,5.);
//TH2F *triangle_electrons = new TH2F("triangle_electrons","Q2 vs Nu for all eventes with a good electron",200,0,10,500,0,10);

// histograms for protons, pions-, kaons, positrons, unrecognised particules
TH2F *betta_momentum_protons = new TH2F("betta_momentum_protons","betta vs. momentum for protons",10,0,6,10,0.,1.1);
TH2F *betta_momentum_kaons = new TH2F("betta_momentum_kaons","betta vs. momentum for kaons",10,0,6,10,0.,1.1);
TH2F *betta_momentum_pions_minus = new TH2F("betta_momentum_pions_minus","betta vs. momentum for - pions",10,0,6,10,0.,1.1);
TH2F *betta_momentum_positrons = new TH2F("betta_momentum_positrons","betta vs. momentum for positrons",10,0,6,10,0.,1.1);
//TH2F *betta_momentum_unrecognised = new TH2F("betta_momentum_unrecognised","betta vs. momentum for unrecognised particules",10,0,6,10,0.,1.1);

// histograms for positivs pions
TH2F *betta_momentum_pions_plus = new TH2F("betta_momentum_pions_plus","betta vs. momentum for pion+",10,0,6,10,0.,1.1);
TH1F *NPHE_pions_plus = new TH1F("NPHE_pions_plus","NPHE for pion +",200,0,200);
TH1F *z_coordinate_pions_plus = new TH1F("z_coordinate_pions_plus","z coordinate for pion+",300,-50,25);
TH1F *z_pions_plus = new TH1F("z_pions_plus","z for pion+",200,0,5);
TH1F *Chi2cc_pions_plus = new TH1F("Chi2cc_pions_plus","Chi2cc for pion+",200,-5,5);
TH2F *triangle_leading_pions_plus = new TH2F("triangle_leading_pions_plus","Q2 vs Nu when leanding pions+",10,0,10,10,0,10);

// histrograms to see the depence of the # of protons with pt
TH1F *varproton_pt2 = new TH1F("varproton_pt2","Np vs pt squared",300,-1,3);
TH1F *varproton_low_energy = new TH1F("varproton2","Np vs pt squared if energy of proton < position pion+",300,-1,3);
TH1F *varproton_pt4 = new TH1F("varproton_pt4","Np vs pt^4",300,-1,3);

// histograms for each good events
TH1F *totalEnergy = new TH1F("totalEnergy","total visible energy",300,-1,6);
TH1F *pperpeny = new TH1F("pperpeny","total pperpen y",300,-50,50);
TH1F *pperpenx = new TH1F("pperpenx","total pperpen x",300,-50,50);
TH1F *pperpen = new TH1F("pperpen","total pperpen value",300,-10,100);

// histograms for leading pions+
TH2F *dependence_z_pt_leading = new TH2F("dependence_z_pt_leading","z vs pt^2 for pions+, z>0.5",20,0,5,20,0,1);
TH2F *dependence_z_nu_leading = new TH2F("dependence_z_nu_leading","z vs.nu for pions+, z>0.5",20,0,5,20,0,1);
//TH1F *z_leading = new TH1F("z_leading","z of pion+, z>0.5:  ",50,0,2);
TH1F *t_moment_leading = new TH1F("t_moment_leading","transverse momentum squared for pions+, z>0.5",300,-1,3);
TH2F *triangle_leading = new TH2F("triangle_leading","Q2 vs Nu for pions+, z>0.5",20,0,10,20,0,10);
TH2F *betta_momentum_leading_pperpnull = new TH2F("betta_momentum_leading_pperpnull","positive pion betta vs. momentum leading pions+ when pperp=0",10,0,6,10,0,1.1);

// histograms for leading pions+, solid target
TH1F *t_moment_solid = new TH1F("t_moment_solid","transverse momentum squared on solid target",300,-1,3);
TH1F *z_pions_plus_solid = new TH1F("z_pions_plus_solid","z for pion+ on solid target",50,0,5);
//TH2F *betta_momentum_solid = new TH2F("betta_momentum_solid","positive pion betta vs. momentum, higher momentum on solid target",300,0,6,300,0,1.1);
TH1F *moment_solid = new TH1F("moment_solid","momentum on solid target",300,0,5);
TH1F *pperpen_leading_solid = new TH1F("pperpen_leading_solid","total pperpen_leading_solid value",300,-50,50);
TH1F *totalEnergy_leading_solid = new TH1F("totalEnergy_leading_solid","total visible energy_leading_solid",300,-1,6);

// histograms for leading pions+, solid target, z>0.5 
TH1F *t_moment_solid_energetic = new TH1F("t_moment_solid_energetic","transverse momentum squared on solid target, z>0.5",300,-1,3);
TH1F *z_pions_plus_solid_energetic = new TH1F("z_pions_plus_solid_energetic","z for pion+ on solid target, z>0.5",50,0,5);
TH1F *moment_solid_energetic = new TH1F("moment_solid_energetic","momentum on solid target, z>0.5",300,0,5);
TH1F *totalEnergy_leading_solid_energetic = new TH1F("totalEnergy_leading_solid_energetic","total visible energy_leading_solid, z>0.5",300,-1,6);
TH1F *pperpen_leading_solid_energetic = new TH1F("pperpen_leading_solid_energetic","total pperpen_leading_solid value, z>0.5",300,-50,50);
TH2F *betta_momentum_solid_energetic = new TH2F("betta_momentum_solid_energetic","positive pion betta vs. momentum, leading on solid target, z>0.5",10,0,6,10,0,1.1);

// histograms for leading pions+, deuterium target
TH1F *t_moment_deuterium = new TH1F("t_moment_deuterium","transverse momentum squared on deuterium target",300,-1,3);
TH1F *z_pions_plus_deuterium = new TH1F("z_pions_plus_deuterium","z for pion+ on deuterium target",50,0,5);
//TH2F *betta_momentum_deuterium = new TH2F("betta_momentum_deuterium","positive pion betta vs. momentum, higher momentum on deuterium target",300,0,6,300,0,1.1);
TH1F *moment_deuterium = new TH1F("moment_deuterium","momentum on deuterium target",300,0,5);
TH1F *pperpen_leading_deuterium = new TH1F("pperpen_leading_deuterium","total pperpen_leading_deuterium value",300,-50,50);
TH1F *totalEnergy_leading_deuterium = new TH1F("totalEnergy_leading_deuterium","total visible energy_leading_deuterium",300,-1,6);

// histograms for leading pions+, deuterium target, z>0.5 
TH1F *t_moment_deuterium_energetic = new TH1F("t_moment_deuterium_energetic","transverse momentum squared on deuteriumd target, z>0.5",300,-1,3);
TH1F *z_pions_plus_deuterium_energetic = new TH1F("z_pions_plus_deuterium_energetic","z for pion+ on deuterium target, z>0.5",50,0,5);
TH1F *moment_deuterium_energetic = new TH1F("moment_deuterium_energetic","momentum on deuterium target, z>0.5",300,0,5);
TH1F *totalEnergy_leading_deuterium_energetic = new TH1F("totalEnergy_leading_deuterium_energetic","total visible energy_leading_deuterium, z>0.5",300,-1,6);
TH1F *pperpen_leading_deuterium_energetic = new TH1F("pperpen_leading_deuterium_energetic","total pperpen_leading_deuterium value, z>0.5",300,-50,50);
//TH2F *betta_momentum_deuterium_energetic = new TH2F("betta_momentum_deuterium_energetic","pion plus betta vs.momentum, leading on deuterium target, z>0.5",300,0,6,300,0,1.1);

//in.open("infile");
in.open("/home/silvestr/tutorials/mytest/filer.tex"); 
//in.open("/mss/clas/eg2a/production/Pass0/v0/ROOT/r041163_00.root");
gBenchmark->Start("while");
while(1)
{
in >> x;
if (!in.good()) break;
R1->AddFile(x);
}
gBenchmark->Show("while");

nentries = (Int_t)R1->GetEntries();
cout << nentries << endl;
    gBenchmark->Start("one_events");
    gBenchmark->Start("leading_if");

gBenchmark->Start("all_events");
  gBenchmark->Start("one_particule");
    TStopwatch timer;
    timer.Start();
gBenchmark->Start("protonloop");
          gBenchmark->Start("deuteriumtarget");
          gBenchmark->Start("solidtarget");
 gBenchmark->Start("msup0");
// loop through all the events
for(i=1;i<nentries;i++) {
  R1->Next();
  presult_ptr=(p_categorize->categorization(R1));
  
  m= 0;
  maxl[0]=-400;
  visibleEnergyEvent=0;
  pperpx=0;
  pperpy=0;
  pperp=0;
  
  // if good event with particule 0 = electron in DIS
  if((*presult_ptr).CompareTo("electron")==0 && partident.Q2(R1,0)>=1/* && partident.Q2(R1,0)<=1.3*/ && partident.W(R1,0)>2){
    
    numberEventsConElec++;


    // loop through all the particuls of each event
       for(l=0;l<10;l++) {
  
    presult=*(presult_ptr+l);
     

      // sum the visible energy
      if(presult.CompareTo("not recognized")!=0){//omiting the neutral particules and the ones giving a wrong beta
	visibleEnergyEvent+= p_pt.Etotal(R1,l);
      }
      // sum the component of the perpendicular momentum
      pperpx+=p_pt.pperpx(R1,l);
      pperpy+=p_pt.pperpy(R1,l);
      
      if(presult.CompareTo("electron")==0) {
	//betta_momentum_electrons->Fill(partident.moment(R1,l),partident.betta(R1,l));
	NPHE_electrons->Fill(partident.nphe(R1,l));
	z_coordinate_electrons->Fill(partident.z(R1,l));
	Q2_electrons->Fill(partident.Q2(R1,l));
	W_electrons->Fill(partident.W(R1,l));
	//triangle_electrons->Fill(partident.niu(R1,l),partident.Q2(R1,l));
      }
      if(presult.CompareTo("pion +")==0) {
	betta_momentum_pions_plus->Fill(partident.moment(R1,l),partident.betta(R1,l));
	NPHE_pions_plus->Fill(partident.nphe(R1,l));
	z_coordinate_pions_plus->Fill(partident.z(R1,l));
	z_pions_plus->Fill(partident.Z(R1,l));
	Chi2cc_pions_plus->Fill(partident.chi2cc(R1,l));
	
	// looking for the highest energy pion +
	if(m==0) {
	  maxl[m]=l;
	  m=1;
	  triangle_leading_pions_plus->Fill(partident.niu(R1,0),partident.Q2(R1,0)); 
	}
	if(m>0 && p_pt.Etotal(R1,l)>p_pt.Etotal(R1,maxl[m-1])){
	  maxl[m]=l;
	  m++;
	  triangle_leading_pions_plus->Fill(partident.niu(R1,0),partident.Q2(R1,0)); 
	}
      }		                     
      if(presult.CompareTo("pion -")==0) {
	betta_momentum_pions_minus->Fill(partident.moment(R1,l),partident.betta(R1,l));
      }
      if(presult.CompareTo("positron")==0){
	betta_momentum_positrons->Fill(partident.moment(R1,l),partident.betta(R1,l));
      }
      if(presult.CompareTo("proton")==0){
	betta_momentum_protons->Fill(partident.moment(R1,l),partident.betta(R1,l));
      }
      if(presult.CompareTo("kaon")==0){
	betta_momentum_kaons->Fill(partident.moment(R1,l),partident.betta(R1,l));
	}
      /*if(presult.CompareTo("not recognized")==0){
	betta_momentum_unrecognised->Fill(partident.moment(R1,l),partident.betta(R1,l));
	}*/
      gBenchmark->Stop("one_particule");
      timer.Stop();
    }   
    gBenchmark->Stop("one_events");
  
    totalEnergy->Fill(visibleEnergyEvent);
    pperpenx->Fill(pperpx);
    pperpeny->Fill(pperpy);
    pperp=sqrt(pperpx*pperpx+pperpy*pperpy);
    pperpen->Fill(pperp);


    // we fill the histograms with the leading pions in each event 
       if(m>0){
	 if(partident.Z(R1,maxl[m-1])>0.5){
	   dependence_z_pt_leading->Fill(pow(p_pt.pt(R1,maxl[m-1]),2),partident.Z(R1,maxl[m-1]));
	   t_moment_leading->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	   triangle_leading->Fill(partident.niu(R1,0),partident.Q2(R1,0)); 
	   dependence_z_nu_leading->Fill(partident.niu(R1,maxl[m-1]),partident.Z(R1,maxl[m-1]));
	   //z_leading->Fill(partident.Z(R1,maxl[m-1]));
	    if(pperp==0) {betta_momentum_leading_pperpnull->Fill(partident.moment(R1,maxl[m-l]),partident.betta(R1,maxl[m-l]));}
	  }    gBenchmark->Stop("msup0");
	     
	 
      // for the solid target
	  if(partident.z(R1,0)>-28){
	    z_pions_plus_solid->Fill(partident.Z(R1,maxl[m-1]));
	    t_moment_solid->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	    moment_solid->Fill(partident.moment(R1,maxl[m-1]));
	    pperpen_leading_solid->Fill(pperp);
	    totalEnergy_leading_solid->Fill(visibleEnergyEvent);
	if(partident.z(R1,maxl[m-1])>0.5){
	  z_pions_plus_solid_energetic->Fill(partident.Z(R1,maxl[m-1]));
	  t_moment_solid_energetic->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	  moment_solid_energetic->Fill(partident.moment(R1,maxl[m-1]));
	  pperpen_leading_deuterium_energetic->Fill(pperp);
	  totalEnergy_leading_solid_energetic->Fill(visibleEnergyEvent);
	  if(partident.betta(R1,maxl[m-1])>0.825){
	    betta_momentum_solid_energetic->Fill(partident.moment(R1,maxl[m-1]),partident.betta(R1,maxl[m-1]));	
	    }
	}
	/*if(partident.betta(R1,maxl[m-1])>0.825){
	  betta_momentum_solid->Fill(partident.moment(R1,maxl[m-1]),partident.betta(R1,maxl[m-1]));	
	  }*/
				    }
          gBenchmark->Stop("solidtarget");
	  
      // for the deuterium target
      if(partident.z(R1,0)<-28){
	z_pions_plus_deuterium->Fill(partident.Z(R1,maxl[m-1]));
	t_moment_deuterium->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	moment_deuterium->Fill(partident.moment(R1,maxl[m-1]));
	pperpen_leading_deuterium->Fill(pperp);
	totalEnergy_leading_deuterium->Fill(visibleEnergyEvent);
	if(partident.z(R1,maxl[m-1])>0.5){
	  z_pions_plus_deuterium_energetic->Fill(partident.Z(R1,maxl[m-1]));
	  t_moment_deuterium_energetic->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	  moment_deuterium_energetic->Fill(partident.moment(R1,maxl[m-1]));
	  pperpen_leading_deuterium_energetic->Fill(pperp);
	  totalEnergy_leading_deuterium_energetic->Fill(visibleEnergyEvent);
	  /*if(partident.betta(R1,maxl[m-1])>0.825){
	    betta_momentum_deuterium_energetic->Fill(partident.moment(R1,maxl[m-1]),partident.betta(R1,maxl[m-1]));	
	    }*/
	  }
	/*if(partident.betta(R1,maxl[m-1])>0.825){
	  betta_momentum_deuterium->Fill(partident.moment(R1,maxl[m-1]),partident.betta(R1,maxl[m-1]));
	  }*/
	}
          gBenchmark->Stop("deuteriumtarget");
	  
      // histograms of the dependance of the apparition of proton with transverse momentum
	  for(l=0;l<10;l++){
	presult=*(presult_ptr+l);
	if(presult.CompareTo("proton")==0){
	  varproton_pt2->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	  varproton_pt4->Fill(pow(p_pt.pt(R1,maxl[m-1]),4));
	  if(partident.moment(R1,l)<partident.moment(R1,maxl[m-1])){
	    varproton_low_energy->Fill(pow(p_pt.pt(R1,maxl[m-1]),2));
	    }
	}     gBenchmark->Stop("protonloop");

			    }
   gBenchmark->Stop("leading_if");
	 }
}
}
    gBenchmark->Print("one_events");
    gBenchmark->Print("one_particule");
 gBenchmark->Print("protonloop");
          gBenchmark->Print("deuteriumtarget");
          gBenchmark->Print("solidtarget");
 gBenchmark->Print("msup0");
     gBenchmark->Print("leading_if");
timer.Print();
      gBenchmark->Show("all_events");

cout << "number of events with a good electron as first particle= " <<numberEventsConElec <<  endl;

plots->Write();
in.close();
return 0;

}