//********************************************************************
//     Macro to look for protons
//     coming from labdas
//********************************************************************

#if !defined( __CINT__) || defined(__MAKECINT__)
#include <TMath.h>
#include <TError.h>
#include <Riostream.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TParticle.h>
#include <TCanvas.h>
#include <TBenchmark.h>
#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>

#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliStack.h"
#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliRunLoader.h"
#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliRun.h"
#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliESD.h"
#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliHeader.h"
#include "/exports1/aliprod1/zampolli/zampolliSoft/AliRoot/AliRoot/STEER/AliGenEventHeader.h"

#endif

extern AliRun *gAlice;
extern TBenchmark *gBenchmark;
extern TROOT *gROOT;

void momrec2events(){
  Double_t pl=0.,pu=6.;
  Int_t nbins = 40;
  Double_t minres=-0.2,maxres=0.2;
  Int_t nbinsres = 400;
  Double_t delta = (pu-pl)/nbins;
  cout << " delta = " << delta << endl; 
  Char_t hnamepi[8];
  Char_t hnameka[8];
  Char_t hnamepr[8];
  Char_t fname[100];
  TH1F * piP = new TH1F("piP","Transverse Momentum distr for pions", nbins, pl, pu);
  TH1F * kaP = new TH1F("kaP","Transverse Momentum distr for kaons", nbins, pl, pu);
  TH1F * prP = new TH1F("prP","Transverse Momentum distr for protons", nbins, pl, pu);
  TH1F * piPcheck = new TH1F("piPcheck","Transverse Momentum distr for pions", nbins, pl, pu);
  TH1F * kaPcheck = new TH1F("kaPcheck","Transverse Momentum distr for kaons", nbins, pl, pu);
  TH1F * prPcheck = new TH1F("prPcheck","Transverse Momentum distr for protons", nbins, pl, pu);
  piP->Sumw2();
  kaP->Sumw2();
  prP->Sumw2();
  TH1F * piR = new TH1F("piR","Transverse Momentum distr for reco pions", nbins, pl, pu);
  TH1F * kaR = new TH1F("kaR","Transverse Momentum distr for reco kaons", nbins, pl, pu);
  TH1F * prR = new TH1F("prR","Transverse Momentum distr for reco protons", nbins, pl, pu);
  piR->Sumw2();
  kaR->Sumw2();
  prR->Sumw2();
  TH1F * piR2 = new TH1F("piR2","Transverse Momentum distr for reco pions", nbins, pl, pu);
  TH1F * kaR2 = new TH1F("kaR2","Transverse Momentum distr for reco kaons", nbins, pl, pu);
  TH1F * prR2 = new TH1F("prR2","Transverse Momentum distr for reco protons", nbins, pl, pu);
  piR2->Sumw2();
  kaR2->Sumw2();
  prR2->Sumw2();
  TH1F * piPgen = new TH1F("piPgen","Reconstructed Transverse Momentum distr for pions", nbins, pl, pu);
  TH1F * kaPgen = new TH1F("kaPgen","Reconstructed Transverse Momentum distr for kaons", nbins, pl, pu);
  TH1F * prPgen = new TH1F("prPgen","Reconstructed Transverse Momentum distr for protons", nbins, pl, pu);
  TH1F * piRpt = new TH1F("piRpt","Reconstructed Transverse Momentum distr for pions", nbins, pl, pu);
  TH1F * kaRpt = new TH1F("kaRpt","Reconstructed Transverse Momentum distr for kaons", nbins, pl, pu);
  TH1F * prRpt = new TH1F("prRpt","Reconstructed Transverse Momentum distr for protons", nbins, pl, pu);
  piRpt->Sumw2();
  kaRpt->Sumw2();
  prRpt->Sumw2();
  TH1F * ptrespi = new TH1F("ptrespi","Transverse Momentum Resolution distr for pions", nbinsres, minres, maxres);
  TH1F * ptreska = new TH1F("ptreska","Transverse Momentum Resolution distr for kaons", nbinsres, minres, maxres);
  TH1F * ptrespr = new TH1F("ptrespr","Transverse Momentum Resolution distr for protons", nbinsres, minres, maxres);

  TH1F * hcorrpi = new TH1F("hcorrpi","Correction for pt efficiency, pi", nbins, pl, pu);
  TH1F * hcorrka = new TH1F("hcorrka","Correction for pt efficiency, ka", nbins, pl, pu);
  TH1F * hcorrpr = new TH1F("hcorrpr","Correction for pt efficiency, pr", nbins, pl, pu);

  TH1F * hcorr2pi = new TH1F("hcorr2pi","Correction for pt efficiency, pi", nbins, pl, pu);
  TH1F * hcorr2ka = new TH1F("hcorr2ka","Correction for pt efficiency, ka", nbins, pl, pu);
  TH1F * hcorr2pr = new TH1F("hcorr2pr","Correction for pt efficiency, pr", nbins, pl, pu);

  TH1F **histospi = new TH1F*[nbins+1];
  TH1F **histoska = new TH1F*[nbins+1];
  TH1F **histospr = new TH1F*[nbins+1];
  for (Int_t i=1;i<=nbins;i++){
    if (i<10) {
      sprintf(hnamepi,"histopi0%i",i);
      sprintf(hnameka,"histoka0%i",i);
      sprintf(hnamepr,"histopr0%i",i);
    }
    else {
      sprintf(hnamepi,"histopi%i",i);
      sprintf(hnameka,"histoka%i",i);
      sprintf(hnamepr,"histopr%i",i);
    }
    histospi[i] = new TH1F(hnamepi,"",nbins,pl,pu);
    histoska[i] = new TH1F(hnameka,"",nbins,pl,pu);
    histospr[i] = new TH1F(hnamepr,"",nbins,pl,pu);
  }
  
  const char *workingdir = gSystem->WorkingDirectory();
  cout << " workingdir = " << workingdir << endl;
  char namedir[600];
  char namesc[600];

  for (Int_t ifile = 0;ifile<50;ifile++){
    sprintf(namesc,"/exports1/aliprod1/zampolli/ProdNewHEAD/Hij_5fm.Run1/Events/Hij_5fm.Evt%i/Fluct80psNew/PidComb_%i.root",ifile,ifile);

    TFile file(namesc,"READ");
    if(file.IsZombie()) continue;
    if (ifile==19) {
      piPcheck -> Add((TH1F*)file.Get("piP"));
      kaPcheck -> Add((TH1F*)file.Get("kaP"));
      prPcheck -> Add((TH1F*)file.Get("prP"));
    }
    piP -> Add((TH1F*)file.Get("piP"));
    kaP -> Add((TH1F*)file.Get("kaP"));
    prP -> Add((TH1F*)file.Get("prP"));
    piR -> Add((TH1F*)file.Get("piR"));
    kaR -> Add((TH1F*)file.Get("kaR"));
    prR -> Add((TH1F*)file.Get("prR"));
  }

  hcorr2pi->Divide(piP,piR,1,1,"b");
  hcorr2ka->Divide(kaP,kaR,1,1,"b");
  hcorr2pr->Divide(prP,prR,1,1,"b");

  for (Int_t ifile = 0;ifile<50;ifile++){
   if (gAlice) {
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice=0;
   }    
    sprintf(namedir,"/exports1/aliprod1/zampolli/ProdNewHEAD/Hij_5fm.Run1/Events/Hij_5fm.Evt%i",ifile);    
    cout << " ------------------------------------ " << endl;
    cout << " ----------FILE = " << namedir << "---------------------- " << endl;
    cout << " ------------------------------------ " << endl;
    sprintf(fname,"%s/galice.root",namedir);
    TArrayF vertex(3);
    AliRunLoader *rl = AliRunLoader::Open(fname);
    if (rl == 0x0) {
      cerr<<"Can not open session"<<endl;
      continue;
    }    
    cout << " rl = " << rl << endl;
    if (rl->LoadHeader()) {
      cerr<<"LoadHeader returned error"<<endl;
      delete rl;
      continue;
    }
    rl->LoadKinematics();
    
    sprintf(fname,"%s/AliESDs.root",namedir);
    TFile *ef=TFile::Open(fname);
    if (!ef || !ef->IsOpen()) {
      ::Error("chisquare","Can't AliESDs.root !"); 
      delete rl;
      continue;
    }
    AliESD* event = new AliESD;
    TTree* tree = (TTree*) ef->Get("esdTree");
    if (!tree) {
      ::Error("AliESDComparison.C", "no ESD tree found");
      delete rl;
      continue;
    }
    tree->SetBranchAddress("ESD", &event);
    Int_t e=0;
    while (tree->GetEvent(e)) {
      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
      
      rl->GetEvent(e);
      
      e++;
      
      Int_t ntrk=event->GetNumberOfTracks();
      AliStack *stack = rl->Stack();
      
      TArrayF vertex(3);
      rl->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
      Int_t npart = stack->GetNtrack();
      cout << " npart = " << npart << endl;
      while (ntrk--) {
	AliESDtrack *t=event->GetTrack(ntrk);
	Int_t label=TMath::Abs(t->GetLabel());
	UInt_t status=AliESDtrack::kESDpid;
	Double_t p[3]={0,0,0};
	if ((t->GetStatus()&status) == status) {
	  if (t->GetConstrainedChi2()<=20){
	    TParticle *MPart=stack->Particle(label);
	    t->GetConstrainedPxPyPz(p);
	    Double_t pt = TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
	    if (pt >= 6) continue;
	    Float_t P = MPart->Pt();  //transverse momentum from stack
	    Int_t mpart  = TMath::Abs(MPart->GetPdgCode());
	    if (mpart == 211){
	      Int_t ibin = (Int_t)((P-pl)/delta)+1;
	      histospi[ibin]->Fill(pt);
	      piR2->Fill(P);
	      ptrespi->Fill((pt-P)/P/P);
	    }
	    else if (mpart == 321){
	      Int_t ibin = (Int_t)((P-pl)/delta)+1;
	      histoska[ibin]->Fill(pt);
	      kaR2->Fill(P);
	      ptreska->Fill((pt-P)/P/P);
	    }
	    else if (mpart == 2212){
	      Int_t ibin = (Int_t)((P-pl)/delta)+1;
	      histospr[ibin]->Fill(pt);
	      prR2->Fill(P);
	      ptrespr->Fill((pt-P)/P/P);
	    }
	  }
	}
      }
    }
    rl->UnloadHeader();
    rl->UnloadKinematics();
    rl->UnloadgAlice();
    if (gAlice) {
      delete gAlice->GetRunLoader();
      delete gAlice;
      gAlice = 0x0;
   }
   delete rl;
   ef->Close();
  }

  
  TH1F *heffpi = new TH1F("heffpi","reco eff histo for pions", nbins, pl, pu);
  TH1F *heffka = new TH1F("heffka","reco eff histo for kaons", nbins, pl, pu);
  TH1F *heffpr = new TH1F("heffpr","reco eff histo for protons", nbins, pl, pu);  
  heffpi->Sumw2();
  heffka->Sumw2();
  heffpr->Sumw2();

  TH1F *heff2pi = new TH1F("heff2pi","reco eff histo for pions", nbins, pl, pu);
  TH1F *heff2ka = new TH1F("heff2ka","reco eff histo for kaons", nbins, pl, pu);
  TH1F *heff2pr = new TH1F("heff2pr","reco eff histo for protons", nbins, pl, pu);  
  heff2pi->Sumw2();
  heff2ka->Sumw2();
  heff2pr->Sumw2();

  TH1F *hothpi = new TH1F("hothpi","reco oth histo for pions", nbins, pl, pu);
  TH1F *hothka = new TH1F("hothka","reco oth histo for kaons", nbins, pl, pu);
  TH1F *hothpr = new TH1F("hothpr","reco oth histo for protons", nbins, pl, pu);  
  hothpi->Sumw2();
  hothka->Sumw2();
  hothpr->Sumw2();

  TH1F *hpireco = new TH1F("hpireco","reco histo for pions", nbins, pl, pu);
  TH1F *hkareco = new TH1F("hkareco","reco histo for kaons", nbins, pl, pu);
  TH1F *hprreco = new TH1F("hprreco","reco histo for protons", nbins, pl, pu);
  hpireco->Sumw2();
  hkareco->Sumw2();
  hprreco->Sumw2();

  for (Int_t i=1;i<=nbins;i++){
    histospi[i]->SetBinContent(1,0);
    histoska[i]->SetBinContent(1,0);
    histospr[i]->SetBinContent(1,0);
    if (histospi[i]->GetBinContent(i) != 0 && piP->GetBinContent(i) != 0){
      Float_t cont = (piP->GetBinContent(i))/(histospi[i]->Integral());
      Float_t cont1 = (histospi[i]->Integral())/(histospi[i]->GetBinContent(i));
      heffpi->SetBinContent(i,cont);
      hcorrpi->SetBinContent(i,cont1);
    }
    if (histoska[i]->GetBinContent(i) != 0 && kaP->GetBinContent(i) != 0){
      Float_t cont = (kaP->GetBinContent(i))/(histoska[i]->Integral());
      Float_t cont1 = (histoska[i]->Integral())/(histoska[i]->GetBinContent(i));
      heffka->SetBinContent(i,cont);
      hcorrka->SetBinContent(i,cont1);
   }
    if (histospr[i]->GetBinContent(i) != 0 && prP->GetBinContent(i) != 0){
      Float_t cont = (prP->GetBinContent(i))/(histospr[i]->Integral());
      Float_t cont1 = (histospr[i]->Integral())/(histospr[i]->GetBinContent(i));
      heffpr->SetBinContent(i,cont);
      hcorrpr->SetBinContent(i,cont1);
    }
  
    Float_t othpi = 0;
    Float_t othka = 0;
    Float_t othpr = 0;
    for (Int_t j=1;j<=nbins;j++){
	othpi += histospi[j]->GetBinContent(i);
	othka += histoska[j]->GetBinContent(i);
	othpr += histospr[j]->GetBinContent(i);
    }
    Float_t cont =0;
    if ((othpi+histospi[i]->GetBinContent(i))!=0){
      cont = histospi[i]->GetBinContent(i)/othpi;
      hothpi->SetBinContent(i,cont);
    }
    if ((othka+histoska[i]->GetBinContent(i))!=0){
      cont = histoska[i]->GetBinContent(i)/othka;
      hothka->SetBinContent(i,cont);
    }
    if ((othpr+histospr[i]->GetBinContent(i))!=0){
      cont = histospr[i]->GetBinContent(i)/othpr;
      hothpr->SetBinContent(i,cont);
    }
  }
  
  heff2pi->Divide(piP,piR2);   
  heff2ka->Divide(kaP,kaR2);   
  heff2pr->Divide(prP,prR2);   
  
  for (Int_t ifile = 19;ifile<20;ifile++){
    sprintf(namedir,"/exports1/aliprod1/zampolli/ProdNewHEAD/Hij_5fm.Run1/Events/Hij_5fm.Evt%i",ifile);    
    cout << " ------------------------------------ " << endl;
    cout << " ----------FILE = " << namedir << "---------------------- " << endl;
    cout << " ------------------------------------ " << endl;
    gSystem->ChangeDirectory(namedir);
    TArrayF vertex(3);
    
    AliRunLoader *rl = AliRunLoader::Open("galice.root");
    
    if (rl == 0x0) {
      cerr<<"Can not open session"<<endl;
      return;
    }
    rl->LoadgAlice();
    gAlice = rl->GetAliRun();
    rl->LoadKinematics();
    rl->LoadHeader();
    
    sprintf(fname,"AliESDs.root");
    TFile *ef=TFile::Open(fname);
    if (!ef || !ef->IsOpen()) {
      ::Error("chisquare","Can't AliESDs.root !"); 
      delete rl;
      return;
    }
    AliESD* event = new AliESD;
    TTree* tree = (TTree*) ef->Get("esdTree");
    if (!tree) {
      ::Error("AliESDComparison.C", "no ESD tree found");
      delete rl;
      return;
    }
    tree->SetBranchAddress("ESD", &event);
    Int_t e=0;
    while (tree->GetEvent(e)) {
      cout<<endl<<endl<<"********* Processing event number: "<<e<<"*******\n";
      
      rl->GetEvent(e);
      
      e++;
      
      Int_t ntrk=event->GetNumberOfTracks();
      AliStack *stack = rl->Stack();
      
      TArrayF vertex(3);
      rl->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
      while (ntrk--) {
	AliESDtrack *t=event->GetTrack(ntrk);
	Int_t label=TMath::Abs(t->GetLabel());
	UInt_t status=AliESDtrack::kESDpid;
	Double_t p[3];
	if ((t->GetStatus()&status) == status) {
	  if (t->GetConstrainedChi2()<=20){
	    TParticle *MPart=stack->Particle(label);
	    t->GetConstrainedPxPyPz(p);
	    Double_t pt = TMath::Sqrt(p[0]*p[0] + p[1]*p[1]);
	    if (pt >= 6) continue;
	    Float_t P = MPart->Pt();  //transverse momentum from stack
	    Int_t mpart  = TMath::Abs(MPart->GetPdgCode());
	    if (mpart == 211){
	      piPgen->Fill(P);
	      piRpt->Fill(pt);
	    }
	    else if (mpart == 321){
	      kaPgen->Fill(P);
	      kaRpt->Fill(pt);
	    }
	    else if (mpart == 2212){
	      prPgen->Fill(P);
	      prRpt->Fill(pt);
	    }
	  }
	}
      }
    }
  

    hpireco->Add(piRpt);
    hpireco->Multiply(hothpi);
    hpireco->Multiply(hcorrpi);
    //  hpireco->Multiply(hcorr2pi);
    hpireco->Multiply(heff2pi);
    hpireco->SetBinContent(1,0);

    hkareco->Add(kaRpt);
    hkareco->Multiply(hothka);
    hkareco->Multiply(hcorrka);
    hkareco->Multiply(heff2ka);
    // hkareco->Multiply(hcorr2ka);
    hkareco->SetBinContent(1,0);

    hprreco->Add(prRpt);
    hprreco->Multiply(hothpr);
    hprreco->Multiply(hcorrpr);
    hprreco->Multiply(heff2pr);
    //hprreco->Multiply(hcorr2pr);
    hprreco->SetBinContent(1,0);

    rl->UnloadHeader();
    rl->UnloadKinematics();
    rl->UnloadgAlice();
    delete rl;
    ef->Close();
  }

  TCanvas *cpi = new TCanvas("cpi","cpi",-2,30,500,500);
  cpi->cd();
  //  heffpi->Draw();
  hcorr2pi->Draw();
  TCanvas *cka = new TCanvas("cka","cka",-2,30,500,500);
  cka->cd();
  hcorr2ka->Draw();
  //heffka->Draw();
  TCanvas *cpr = new TCanvas("cpr","cpr",-2,30,500,500);
  cpr->cd();
  hcorr2pr->Draw();
  //heffpr->Draw();
  
  TCanvas *cpiR = new TCanvas("cpiR","cpiR",-2,30,500,500);
  cpiR->cd();
  hpireco->Draw("hist");
  TCanvas *ckaR = new TCanvas("ckaR","ckaR",-2,30,500,500);
  ckaR->cd();
  hkareco->Draw("hist");
  TCanvas *cprR = new TCanvas("cprR","cprR",-2,30,500,500);
  cprR->cd();
  hprreco->Draw("hist");

  
  TCanvas *cpiPgen = new TCanvas("cpiPgen","cpiPgen",-2,30,500,500);
  cpiPgen->cd();
  //  piPgen->Draw("hist");
  cout << "piPcheck = " << piPcheck << endl;
  piPcheck->Draw("hist");
  TCanvas *ckaPgen = new TCanvas("ckaPgen","ckaPgen",-2,30,500,500);
  ckaPgen->cd();
  //  kaPgen->Draw("hist");
  kaPcheck->Draw("hist");
  TCanvas *cprPgen = new TCanvas("cprPgen","cprPgen",-2,30,500,500);
  cprPgen->cd();
  //  prPgen->Draw("hist");
  prPcheck->Draw("hist");
  //  prP->Draw();
  

  TCanvas *cptrespi = new TCanvas("cptrespi","cptrespi",-2,30,500,500);
  cptrespi->cd();
  ptrespi->Draw();
  TCanvas *cptreska = new TCanvas("cptreska","cptreska",-2,30,500,500);
  cptreska->cd();
  ptreska->Draw();
  TCanvas *cptrespr = new TCanvas("cptrespr","cptrespr",-2,30,500,500);
  cptrespr->cd();
  ptrespr->Draw();
  
  
  gSystem->ChangeDirectory("/exports1/aliprod1/zampolli/HQ06/Systematics/");

  return;
}