PWGLF/SPECTRA/PiKaPr/TestAOD/MainAnalysis.C

   1 class AliPID;
   2
   3 TString Charge[]={"Pos","Neg"};
   4 TString Sign[]={"Plus","Minus"};
   5 TString Particle[]={"Pion","Kaon","Proton"};
   6 Int_t Color[3]={1,2,4};
   7 Int_t Marker[6]={20,21,22,24,25,26};
   8 Double_t projl[3]={0,20,80};
   9 Double_t proju[3]={10,40,90};
  10 Double_t Range[3]={0.3,0.3,0.5}; // LowPt range for pi k p
  11 enum ECharge_t {
  12   kPositive,
  13   kNegative,
  14   kNCharges
  15 };
  16
  17 void MainAnalysis()  {
  18
  19   gSystem->Load("libCore.so");
  20   //gSystem->Load("libGeom.so");
  21   gSystem->Load("libPhysics.so");
  22   //gSystem->Load("libVMC");
  23   gSystem->Load("libTree");
  24   //gSystem->Load("libProof");
  25   gSystem->Load("libMatrix");
  26   gSystem->Load("libSTEERBase");
  27   gSystem->Load("libESD");
  28   gSystem->Load("libAOD");
  29   gSystem->Load("libANALYSIS");
  30   gSystem->Load("libOADB");
  31   gSystem->Load("libANALYSISalice");
  32   gSystem->Load("libTENDER");
  33   gSystem->Load("libCORRFW");
  34   //gSystem->Load("libPWG0base");
  35   //gSystem->Load("libMinuit");
  36   gSystem->Load("libPWGTools");
  37   gSystem->Load("libPWGLFSPECTRA");
  38   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  39   gROOT->LoadMacro("QAPlots.C");
  40   gStyle->SetOptStat(000000);
  41   // Set Masses
  42   Double_t mass[3];
  43   mass[0]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  44   mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
  45   mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
  46
  47   //TString fold="test";
  48   TString fold="LEGOTrain_20120625";
  49
  50   //LOOP OVER SELECTION
  51   //                            0    1    2    3    4    5    6    7    8    9
  52   Double_t CentCutMin[4]= {     0,  30,  30,  30};
  53   Double_t CentCutMax[4]= {     5,  40,  40,  40};
  54   Double_t QvecCutMin[4]={      0,   0,   0, 1.5};
  55   Double_t QvecCutMax[4]={    100, 100, 0.4, 100};
  56   Double_t EtaMin[4]={       -0.8,-0.8,-0.8,-0.8};
  57   Double_t EtaMax[4]={        0.8, 0.8, 0.8, 0.8};
  58   Double_t Nsigmapid=3.;
  59   Double_t pt=10.;
  60   Double_t p=10.;
  61   Double_t y=.5;
  62   Double_t ptTofMatch=.6;
  63   UInt_t trkbit=1;
  64   UInt_t trkbitQVector=1;
  65   Bool_t UseCentPatchAOD049=kFALSE;
  66   Double_t DCA=100000;
  67   UInt_t minNclsTPC=70;
  68   Int_t nrebin=0;
  69   TString opt="";
  70
  71   Int_t icut=1;
  72   if(icut==0)Int_t ibinToCompare=0;
  73   else Int_t ibinToCompare=4;
  74
  75   TString sname_data=Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
  76   //MC we do not use the cut on Qvector
  77   TString sname_mc=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
  78
  79   TString dataFile = Form("output/%s/AnalysisResultsDATA.root",fold.Data());
  80   TString mcFile = Form("output/%s/AnalysisResultsMC.root",fold.Data());
  81   //TString mcFile = Form("output/%s/AnalysisResultsAMPT.root",fold.Data());
  82
  83   Printf("\n\n-> Creating Output file Res_%s_%s.root",sname_data.Data(),fold.Data());
  84   TFile * fout=new TFile(Form("results/Res_%s_%s.root",sname_data.Data(),fold.Data()),"RECREATE");
  85
  86   gStyle->SetPalette(1);
  87   // Open root MC file and get classes
  88   cout << "Analysis Macro" << endl;
  89   cout << "  > Reading MC data" << endl;
  90   TFile *f_mc = TFile::Open(mcFile.Data());
  91   f_mc->Print("");
  92   Printf("sname_MC: %s",sname_mc.Data());
  93   TDirectoryFile *Dir_mc=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc.Data()));
  94   AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) Dir_mc->Get("SpectraHistos");
  95   AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) Dir_mc->Get("Event Cuts");
  96   AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) Dir_mc->Get("Track Cuts");
  97   // print info about mc track and Event cuts
  98   cout << " -- Info about MC -- "<< endl;
  99   ecuts_mc->PrintCuts();
 100   tcuts_mc->PrintCuts();
 101   // proceed likewise for data
 102   // Open root DATA file and get classes
 103   cout << "  > Reading DATA data" << endl;
 104   TFile *f_data = TFile::Open(dataFile.Data());
 105   f_data->Print("");
 106   Printf("sname_DATA: %s",sname_data.Data());
 107   TDirectoryFile *Dir_data=(TDirectoryFile*)f_data->Get(Form("%s",sname_data.Data()));
 108   AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) Dir_data->Get("SpectraHistos");
 109   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) Dir_data->Get("Event Cuts");
 110   AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) Dir_data->Get("Track Cuts");
 111   // print info about data track and Event cuts
 112   cout << " -- Info about DATA -- "<< endl;
 113   ecuts_data->PrintCuts();
 114   tcuts_data->PrintCuts();
 115
 116
 117   QAPlots(hman_data,hman_mc,ecuts_data,ecuts_mc,tcuts_data,tcuts_mc,fout);
 118
 119   //efficiencies
 120   Printf("\n\n-> Calculating MC Correction Factors");
 121   TH1F *CorrFact[6];
 122   TCanvas *ceff=new TCanvas("ceff","ceff",700,500);
 123   Bool_t UseMCDCACorrection=kTRUE;
 124   for(Int_t icharge=0;icharge<2;icharge++){
 125     for(Int_t ipart=0;ipart<3;ipart++){
 126       Int_t index=ipart+3*icharge;
 127       //BE CAREFUL! depending on the efficiency you choose, you must change the DCA correction (data or data/mc)
 128
 129       TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()); //MC correction for Prim+Cont+Eff
 130       //TString hname=Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());//MC correction for Prim+Eff if (idRec == idGen)
 131       //TString hname=Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorretion=kFALSE; //MC correction for Cont+Eff. filled with idRec
 132       //TString hname=Form("hHistPtRecPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; //MC correction for Cont+Eff. filled with idGen
 133       //TString hname=Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());UseMCDCACorrection=kFALSE; // Pure MC efficiency for Prim. BE CAREFUL WITH MUONS!!! (idRec == idGen)
 134       Printf("Getting %s",hname.Data());
 135       CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 136       CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 137       CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 138       CorrFact[index]->SetMarkerStyle(Marker[index]);
 139       CorrFact[index]->SetMarkerColor(Color[ipart]);
 140       CorrFact[index]->SetLineColor(Color[ipart]);
 141       hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 142       Printf("... and divide it by %s",hname.Data());
 143       CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
 144       if(index==0)CorrFact[index]->DrawClone();
 145       else CorrFact[index]->DrawClone("same");
 146     }
 147   }
 148   TFile *fESD=new TFile("EffAlex/pionEffPbPb.root");
 149   TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0");
 150   hEffESD->DrawClone("same");
 151   gPad->BuildLegend()->SetFillColor(0);
 152
 153   TFile*f=new TFile("Eff.root","recreate");
 154   for(Int_t icharge=0;icharge<2;icharge++){
 155     for(Int_t ipart=0;ipart<3;ipart++){
 156           Int_t index=ipart+3*icharge;
 157           CorrFact[index]->Write();
 158     }
 159   }
 160
 161   //Normalization
 162   printf("\n\n-> Spectra Normalization");
 163   Double_t events_data =  ecuts_data->NumberOfEvents();
 164   Printf(": accepted events: %.0f",events_data);
 165   Double_t events_mc =  ecuts_mc->NumberOfEvents();
 166   Printf(": accepted events (MC): %.0f",events_mc);
 167
 168   //divide RAW for Correction Factor
 169   Printf("\n\n-> Using MC correction factor to correct RAW spectra");
 170   TH1F *Spectra[6];
 171   for(Int_t icharge=0;icharge<2;icharge++){
 172     for(Int_t ipart=0;ipart<3;ipart++){
 173       Int_t index=ipart+3*icharge;
 174       TString hname=Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 175       printf("Getting %s",hname.Data());
 176       Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 177       Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 178       Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 179       Spectra[index]->SetMarkerStyle(Marker[index]);
 180       Spectra[index]->SetMarkerColor(Color[ipart]);
 181       Spectra[index]->SetLineColor(Color[ipart]);
 182       Printf("... and divide it by %s",hname.Data());
 183       Spectra[index]->Divide(CorrFact[index]);
 184       Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION
 185     }
 186   }
 187
 188
 189
 190   //Put Bin Content = 0 for bin below the the range
 191   for(Int_t icharge=0;icharge<2;icharge++){
 192     for(Int_t ipart=0;ipart<3;ipart++){
 193       Int_t index=ipart+3*icharge;
 194       for(Int_t ibin=0;ibin<Spectra[index]->GetNbinsX();ibin++){
 195         if(Spectra[index]->GetBinCenter(ibin)<Range[ipart]){
 196           Spectra[index]->SetBinContent(ibin,0);
 197           Spectra[index]->SetBinError(ibin,0);
 198         }
 199       }
 200     }
 201   }
 202
 203   //DCA Correction with the "right" DCA sample
 204   DCACorrection(Spectra,hman_data,hman_mc,UseMCDCACorrection,fout);
 205
 206   //DCA Correction forcing loose DCA
 207   // TString fold_LooseDCA="5SigmaPID_AOD046_FilterBit5";
 208   // TString dataFile_LooseDCA = Form("output/%s/Pt.AOD.1._data_ptcut_%s.root",fold_LooseDCA.Data(),sname.Data());
 209   // TString mcFile_LooseDCA =Form("output/%s/Pt.AOD.1._MC_%s.root",fold_LooseDCA.Data(),sname.Data());
 210   // TFile *_mc_LooseDCA = TFile::Open(mcFile_LooseDCA.Data());
 211   // AliSpectraAODHistoManager* hman_mc_LooseDCA = (AliSpectraAODHistoManager*) _mc_LooseDCA->Get("SpectraHistos");
 212   // TFile *_data_LooseDCA = TFile::Open(dataFile_LooseDCA.Data());
 213   // AliSpectraAODHistoManager* hman_data_LooseDCA = (AliSpectraAODHistoManager*) _data_LooseDCA->Get("SpectraHistos");
 214   // DCACorrection(Spectra,hman_data_LooseDCA,hman_mc_LooseDCA,UseMCDCACorrection);
 215
 216   //GFCorrection
 217   GFCorrection(Spectra,tcuts_data,fout);
 218
 219   TH1F *MCTruth[6];
 220   for(Int_t icharge=0;icharge<2;icharge++){
 221     for(Int_t ipart=0;ipart<3;ipart++){
 222       Int_t index=ipart+3*icharge;
 223       hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 224       MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone();
 225       MCTruth[index]->SetName(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 226       MCTruth[index]->SetTitle(Form("MCTruth_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 227       MCTruth[index]->Scale(1./events_mc,"width");//NORMALIZATION
 228     }
 229   }
 230
 231
 232   //Matching efficiency in data and Monte Carlo
 233   TCanvas *cMatchingEff=new TCanvas("MatchingEff","MatchingEff",700,500);
 234   TH1F *hMatcEffPos_data=(TH1F*)tcuts_data->GetHistoNMatchedPos();
 235   hMatcEffPos_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedPos());
 236   hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
 237   TH1F *hMatcEffNeg_data=(TH1F*)tcuts_data->GetHistoNMatchedNeg();
 238   hMatcEffNeg_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedNeg());
 239   hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
 240   hMatcEffNeg_data->SetLineColor(2);
 241   TH1F *hMatcEffPos_mc=(TH1F*)tcuts_mc->GetHistoNMatchedPos();
 242   hMatcEffPos_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedPos());
 243   hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
 244   hMatcEffPos_mc->SetLineStyle(2);
 245   TH1F *hMatcEffNeg_mc=(TH1F*)tcuts_mc->GetHistoNMatchedNeg();
 246   hMatcEffNeg_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedNeg());
 247   hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
 248   hMatcEffNeg_mc->SetLineColor(2);
 249   hMatcEffNeg_mc->SetLineStyle(2);
 250   cMatchingEff->Divide(1,2);
 251   cMatchingEff->cd(1);
 252   gPad->SetGridy();
 253   gPad->SetGridx();
 254   hMatcEffPos_data->DrawClone("lhist");
 255   hMatcEffNeg_data->DrawClone("lhistsame");
 256   hMatcEffPos_mc->DrawClone("lhistsame");
 257   hMatcEffNeg_mc->DrawClone("lhistsame");
 258   gPad->BuildLegend()->SetFillColor(0);
 259   hMatcEffPos_data->Divide(hMatcEffPos_mc);
 260   hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
 261   cMatchingEff->cd(2);
 262   gPad->SetGridy();
 263   gPad->SetGridx();
 264   hMatcEffPos_data->DrawClone("lhist");
 265   hMatcEffNeg_data->DrawClone("lhistsame");
 266   TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",2.5,5);
 267   hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
 268   pol0MatchPos_data->DrawClone("same");
 269   TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",2.5,5);
 270   hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
 271   pol0MatchNeg_data->SetLineColor(2);
 272   pol0MatchNeg_data->DrawClone("same");
 273   Float_t ScalingMatchingPos=pol0MatchPos_data->GetParameter(0);
 274   Float_t ScalingMatchingNeg=pol0MatchNeg_data->GetParameter(0);
 275
 276   //Correction spectra for matching efficiency
 277   //For the moment I'm using the inclusive correction
 278   for(Int_t ipart=0;ipart<3;ipart++){
 279     for(Int_t ibin=1;ibin<Spectra[ipart]->GetNbinsX();ibin++){
 280       Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
 281       if(ptspectra<tcuts_data->GetPtTOFMatching())continue;
 282       //Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/hMatcEffPos_data->GetBinContent(hMatcEffPos_data->FindBin(ptspectra))));
 283       //Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/hMatcEffNeg_data->GetBinContent(hMatcEffNeg_data->FindBin(ptspectra))));
 284       Spectra[ipart]->SetBinContent(ibin,( Spectra[ipart]->GetBinContent(ibin)/ScalingMatchingPos));
 285       Spectra[ipart+3]->SetBinContent(ibin,( Spectra[ipart+3]->GetBinContent(ibin)/ScalingMatchingNeg));
 286     }
 287   }
 288
 289   //Drawing Final Spectra
 290   TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500);
 291   gPad->SetGridy();
 292   gPad->SetGridy();
 293   gPad->SetLogy();
 294   for(Int_t icharge=0;icharge<2;icharge++){
 295     for(Int_t ipart=0;ipart<3;ipart++){
 296       Int_t index=ipart+3*icharge;
 297       if(index==0)Spectra[index]->DrawClone();
 298       else Spectra[index]->DrawClone("same");
 299     }
 300   }
 301   gPad->BuildLegend()->SetFillColor(0);
 302
 303
 304   //if Bin 0-5% with no cut ratio with combined analysis
 305   TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500);
 306   CratioComb->Divide(3,2);
 307   TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
 308                        Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
 309   TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120412.root");
 310   TH1F *Spectra_copy[6]=0x0;
 311   for(Int_t icharge=0;icharge<2;icharge++){
 312     for(Int_t ipart=0;ipart<3;ipart++){
 313       Int_t index=ipart+3*icharge;
 314       TH1F *htmp=(TH1F*)((TH1F*)Spectra[index])->Clone("");
 315       htmp->GetXaxis()->SetRangeUser(0,5);
 316       TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone();
 317       CratioComb->cd(ipart+1);
 318       gPad->SetGridy();
 319       gPad->SetGridx();
 320       //for(Int_t ibin=1;ibin<hcomb->GetNbinsX();ibin++)hcomb->SetBinError(ibin,0);
 321
 322       if(icharge==0)htmp->DrawClone();
 323       else htmp->DrawClone("same");
 324       //MCTruth[index]->DrawClone("same");
 325       hcomb->DrawClone("same");
 326       htmp->Divide(hcomb);
 327       htmp->SetMaximum(1.3);
 328       htmp->SetMinimum(0.7);
 329       CratioComb->cd(ipart+4);
 330       gPad->SetGridy();
 331       gPad->SetGridx();
 332       if(icharge==0)htmp->DrawClone();
 333       else htmp->DrawClone("same");
 334     }
 335   }
 336   //comparison with charged hadron
 337   Printf("\n\n-> ChargedHadron comparison");
 338   TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
 339   TH1F* hMatchCorrectionAllCh=(TH1F*)hMatcEffPos_data->Clone("hMatchCorrectionAllCh");
 340   hMatchCorrectionAllCh->Add(hMatcEffNeg_data); //correction for Matching efficiency!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 341   hMatchCorrectionAllCh->Scale(0.5);
 342   for(Int_t ibin=1;ibin<hChHad_data->GetNbinsX();ibin++){
 343     Float_t ptch=hChHad_data->GetBinCenter(ibin);
 344     if(ptch<tcuts_data->GetPtTOFMatching())continue;
 345     hChHad_data->SetBinContent(ibin,2*(hChHad_data->GetBinContent(ibin)/(ScalingMatchingPos+ScalingMatchingNeg)));
 346   }
 347   //fraction of sec in MC
 348   TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
 349   TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
 350   for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
 351     Double_t en_data=hChHad_data->GetBinContent(ibin);
 352     Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
 353     Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
 354     if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
 355     //Printf("Before: %.0f After: %.0f  fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
 356   }
 357   hPrimRec_mc->Divide(hAllRec_mc);
 358   //efficiency for primaries
 359   TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
 360   hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
 361   TCanvas *cAllChFactors=new TCanvas("cAllChFactors","cAllChFactors",700,500);
 362   cAllChFactors->Divide(1,2);
 363   cAllChFactors->cd(1);
 364   gPad->SetGridy();
 365   gPad->SetGridx();
 366   hPrimRec_mc->SetTitle("Prim/All, charged hadron pure MC");
 367   hPrimRec_mc->DrawClone("lhist");
 368   gPad->BuildLegend()->SetFillColor(0);
 369   cAllChFactors->cd(2);
 370   gPad->SetGridy();
 371   gPad->SetGridx();
 372   hEff_mc->SetTitle("Efficiency for Primaries, charged hadron pure MC");
 373   hEff_mc->DrawClone("lhist");
 374   gPad->BuildLegend()->SetFillColor(0);
 375   //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
 376   hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
 377   hChHad_data->Divide(hEff_mc);//Efficiency
 378   hChHad_data->Scale(1./(TMath::Abs(tcuts_data->GetEtaMin())+TMath::Abs(tcuts_data->GetEtaMax())));
 379   hChHad_data->SetTitle("All Ch from AOD");
 380   hChHad_data->SetName("AllCh");
 381   TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
 382   cAllCh->Divide(1,2);
 383   cAllCh->cd(1);
 384   gPad->SetGridy();
 385   gPad->SetGridx();
 386   hChHad_data->DrawClone();
 387   TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
 388   TH1F *hCh=fCh->Get("hpt_c0_5");
 389   hCh->SetTitle("All Ch from Jacek");
 390   hCh->SetMarkerColor(2);
 391   hCh->SetLineColor(2);
 392   //invariant yield
 393   for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
 394     hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
 395     hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
 396   }
 397   hCh->DrawClone("same");
 398   gPad->BuildLegend()->SetFillColor(0);
 399   TH1F *gRatio=AliPWGHistoTools::MyDivideHistosDifferentBins(hChHad_data,hCh);
 400   gRatio->SetMaximum(1.3);
 401   gRatio->SetMinimum(.7);
 402   cAllCh->cd(2);
 403   gPad->SetGridy();
 404   gPad->SetGridx();
 405   gRatio->Print("all");
 406   gRatio->SetBinContent(48,1.02);
 407   gRatio->SetBinContent(49,1.022);
 408   gRatio->SetBinContent(50,1.021);
 409   gRatio->SetBinContent(51,1.026);
 410   gRatio->SetBinContent(52,1.027);
 411   gRatio->SetBinError(48,0.056);
 412   gRatio->SetBinError(49,0.056);
 413   gRatio->SetBinError(50,0.056);
 414   gRatio->SetBinError(51,0.056);
 415   gRatio->SetBinError(52,0.056);
 416   gRatio->DrawClone("");
 417
 418   //fitting
 419   TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
 420   gPad->SetGridy();
 421   gPad->SetGridx();
 422   hChHad_data->DrawClone();
 423   //Fitting the sum of all particles
 424   AliPWGFunc * fm = new AliPWGFunc;
 425   fm->SetVarType(AliPWGFunc::kdNdpt);
 426   Float_t fitmin = 0.3;
 427   Float_t fitmax = 3;
 428   TF1 * func = 0;
 429   Int_t normPar = 0;
 430   func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
 431   func->SetParLimits(1, 0.1, 0.99);
 432   func->SetParLimits(2, 0.01, 1);
 433   func->SetParLimits(3, 0.01, 2);
 434   TH1F * hToFit = hChHad_data;
 435   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
 436   // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
 437   // //   cout << " FIT ERROR " << endl;
 438   // //   return;
 439   // // }
 440   Double_t yieldTools, yieldETools;
 441   Double_t partialYields[3],partialYieldsErrors[3];
 442   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
 443   func->DrawClone("same");
 444   Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
 445   //Fit All Charged
 446   hToFit = hCh;
 447   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
 448   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
 449   func->SetLineColor(2);
 450   hCh->DrawClone("same");
 451   func->DrawClone("same");
 452   gPad->BuildLegend()->SetFillColor(0);
 453   Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
 454   //sumID vs AllCh
 455   //Convert spectra to dNdeta and sum
 456   TH1F * hsum = (TH1F*) Spectra[0]->Clone();
 457   hsum->Reset("all");
 458   Double_t epsilon = 0.0001;
 459   for(Int_t icharge = 0; icharge < 2; icharge++){
 460     for(Int_t ipart = 0; ipart < 3; ipart++){
 461       Int_t index=ipart+3*icharge;
 462       TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
 463       Int_t nbin = htmp->GetNbinsX();
 464       for(Int_t ibin = 1; ibin <= nbin; ibin++){
 465         Double_t pt = htmp->GetBinCenter(ibin);
 466         Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
 467         Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
 468         //Printf("jacobian: %f pt:%f   BinContent:%f  mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
 469         htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
 470         htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
 471         Int_t ibinSum = hsum->FindBin(pt);
 472         if ( htmp->GetBinContent(ibin) > 0 &&
 473              (TMath::Abs(htmp->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
 474               TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
 475              ) {
 476           cout << "DISCREPANCY IN BIN RANGES" << endl;
 477           cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
 478                << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
 479                << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
 480           cout << "" << endl;
 481         }
 482       }
 483       hsum->Add(htmp);
 484     }
 485   }
 486   hsum->SetTitle("Sum ID from AOD");
 487   TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",700,500);
 488   cChargHadComp->Divide(1,2);
 489   cChargHadComp->cd(1);
 490   gPad->SetGridy();
 491   gPad->SetGridx();
 492   hsum->DrawClone();
 493   hToFit = hsum;
 494   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
 495   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
 496   func->SetLineColor(2);
 497   Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
 498   hChHad_data->SetMarkerColor(2);
 499   hChHad_data->SetLineColor(2);
 500   hChHad_data->DrawClone("same");
 501   gPad->BuildLegend()->SetFillColor(0);
 502   cChargHadComp->cd(2);
 503   gPad->SetGridy();
 504   gPad->SetGridx();
 505   TH1F *hRatio=AliPWGHistoTools::DivideHistosDifferentBins(hsum,hChHad_data);
 506   hRatio->SetMaximum(1.2);
 507   hRatio->SetMinimum(.8);
 508   hRatio->DrawClone("");
 509
 510   //saving spectra
 511   fout->cd();
 512   for(Int_t icharge=0;icharge<2;icharge++){
 513     for(Int_t ipart=0;ipart<3;ipart++){
 514       Int_t index=ipart+3*icharge;
 515       Spectra[index]->Write();
 516       CorrFact[index]->Write();
 517       MCTruth[index]->Write();
 518     }
 519   }
 520   hChHad_data->Write();
 521   cRaw->Write();
 522   ceff->Write();
 523   cMatchingEff->Write();
 524   cspectra->Write();
 525   CratioComb->Write();
 526   cAllChFactors->Write();
 527   cAllCh->Write();
 528   cFitChargHad->Write();
 529   cChargHadComp->Write();
 530
 531   //fout->Close();
 532   //delete fout;
 533
 534 }
 535
 536
 537
 538
 539 void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc,Bool_t UseMCDCACorrection,TFile *fout){
 540   printf("\n\n-> DCA Correction");
 541
 542   Double_t FitRange[2]={-3.,3.};
 543   Int_t nrebin=2;
 544   Printf("\DCACorr");
 545   TH1F *hcorrection[2];
 546   TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
 547   TCanvas *cRatiocorrection=new TCanvas("DCARatiocorrection","DCARatiocorrection",700,500);
 548   cRatiocorrection->Divide(2,1);
 549   ccorrection->Divide(2,1);
 550   TString sample[2]={"data","mc"};
 551   for(Int_t icharge=0;icharge<2;icharge++){
 552     for(Int_t ipart=0;ipart<3;ipart++){
 553       Int_t index=ipart+3*icharge;
 554       for(Int_t isample=0;isample<2;isample++){
 555         TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500);
 556         hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
 557         hcorrection[isample]->Reset("all");
 558         cDCA->Divide(16,4);
 559         for(Int_t ibin_data=6;ibin_data<53;ibin_data++){
 560           Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data);
 561           Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data);
 562           if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 563           if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 564           TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 565           TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 566           TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("hHistPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 567           Double_t minentries=1;
 568           if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue;
 569           // hmc1->Rebin(nrebin);
 570           // hmc2->Rebin(nrebin);
 571           // hmc3->Rebin(nrebin);
 572           // hToFit->Rebin(nrebin);
 573           //Data and MC can have different stat
 574           hToFit->Sumw2();
 575           hmc1->Sumw2();
 576           hmc2->Sumw2();
 577           hmc3->Sumw2();
 578           hToFit->Scale(1./hToFit->GetEntries());
 579           Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries();
 580           hmc1->Scale(1./normMC);
 581           hmc2->Scale(1./normMC);
 582           hmc3->Scale(1./normMC);
 583           cDCA->cd(ibin_data-5);
 584           gPad->SetGridy();
 585           gPad->SetGridx();
 586           gPad->SetLogy();
 587           hToFit->DrawClone("lhist");
 588           hmc3->DrawClone("lhist");
 589           TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
 590           mc->Add(hmc1);
 591           mc->Add(hmc2);
 592           mc->Add(hmc3);
 593           TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
 594           fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 595           fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 596           fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 597           fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
 598           hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
 599           hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data()));
 600           Int_t status = fit->Fit();               // perform the fit
 601           cout << "fit status: " << status << endl;
 602           if (status == 0) {                       // check on fit status
 603             TH1F* result = (TH1F*) fit->GetPlot();
 604             TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
 605             TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
 606             TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2);
 607
 608             Double_t v1=0,v2=0,v3=0;
 609             Double_t ev1=0,ev2=0,ev3=0;
 610             //first method, use directly the fit result
 611             fit->GetResult(0,v1,ev1);
 612             fit->GetResult(1,v2,ev2);
 613             fit->GetResult(2,v3,ev3);
 614             result->Scale(hToFit->GetSumOfWeights()/result->GetSumOfWeights());
 615             PrimMCPred->Scale(hToFit->GetSumOfWeights()*v1/PrimMCPred->GetSumOfWeights());
 616             secStMCPred->Scale(hToFit->GetSumOfWeights()*v2/secStMCPred->GetSumOfWeights());
 617             secMCPred->Scale(hToFit->GetSumOfWeights()*v3/secMCPred->GetSumOfWeights());
 618             //second method, integrated the MC predisction, it should give the same as the first method
 619             // v1=PrimMCPred->Integral();
 620             // v2=secStMCPred->Integral();
 621             // v3=secMCPred->Integral();
 622             //Printf("\n\n\n\n\nv1:%f  v2:%f  v3:%f   ev1:%f  ev2:%f  ev3:%f   ",v1,v2,v3,ev1,ev2,ev3);
 623             hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3));
 624             hcorrection[isample]->SetBinError(ibin_data,0);
 625             //Drawing section
 626             PrimMCPred->SetLineColor(2);
 627             secStMCPred->SetLineColor(6);
 628             secMCPred->SetLineColor(4);
 629             hToFit->SetMinimum(0.0001);
 630             hToFit->DrawClone("lhist");
 631             result->SetTitle("Fit result");
 632             result->DrawClone("lhistsame");
 633             PrimMCPred->DrawClone("lhistsame");
 634             secStMCPred->DrawClone("lhistsame");
 635             secMCPred->DrawClone("lhistsame");
 636           }
 637           else{
 638             hcorrection[isample]->SetBinContent(ibin_data,1);
 639             hcorrection[isample]->SetBinError(ibin_data,0);
 640           }
 641           Printf("deleting");
 642           delete hToFit;
 643         }
 644
 645         ccorrection->cd(icharge+1);
 646         gPad->SetGridy();
 647         gPad->SetGridx();
 648         hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data()));
 649         hcorrection[isample]->SetLineWidth(2);
 650         hcorrection[isample]->SetLineColor(Color[ipart]);
 651         hcorrection[isample]->SetLineStyle(isample+1);
 652         hcorrection[isample]->SetMarkerColor(Color[ipart]);
 653         hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,5);
 654         if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist");
 655         else hcorrection[isample]->DrawClone("lhistsame");
 656         // smooth the DCA correction
 657         // TF1 *fitfun = new TF1("fitfun","[0]+[1]*x^[2]+[3]*x^[4]",0.2,2.5);
 658         // hcorrection[isample]->Fit(fitfun,"WRN","N",0.35,2);
 659         // fitfun->SetLineWidth(1.5);
 660         // fitfun->SetLineColor(Color[ipart]);
 661         // fitfun->SetLineStyle(isample);
 662         // fitfun->DrawClone("same");
 663         // for(Int_t ibin=1;ibin<30;ibin++){
 664         //   hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin)));
 665         // }
 666         fout->cd();
 667         cDCA->Write();
 668       }
 669       Spectra[index]->Multiply(hcorrection[0]);//multiply for data
 670       Printf("DCACorrection for DATA used: Spectra[index]->Multiply(hcorrection[0])")
 671         if(UseMCDCACorrection){
 672           Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo
 673           Printf("DCACorrection for MC used: Spectra[index]->Divide(hcorrection[1]")
 674             }
 675       cRatiocorrection->cd(icharge+1);
 676       gPad->SetGridy();
 677       gPad->SetGridx();
 678       hcorrection[0]->Divide(hcorrection[1]);
 679       if(ipart==0)hcorrection[0]->DrawClone("lhist");
 680       else hcorrection[0]->DrawClone("lhistsame");
 681     }
 682   }
 683   ccorrection->cd(1);
 684   gPad->BuildLegend()->SetFillColor(0);
 685   ccorrection->cd(2);
 686   gPad->BuildLegend()->SetFillColor(0);
 687   fout->cd();
 688   cRatiocorrection->Write();
 689   ccorrection->Write();
 690 }
 691
 692
 693 //////////////
 694
 695 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
 696   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
 697   return TMath::ASinH(pt / mt * TMath::SinH(eta));
 698 }
 699 ///////////////////////
 700 void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data, TFile * fout){
 701   //Geant/Fluka Correction
 702   Printf("\nGF correction for Kaons");
 703   //Getting GF For Kaons in TPC
 704   TGraph *gGFCorrectionKaonPlus=new TGraph();
 705   gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
 706   gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
 707   TGraph *gGFCorrectionKaonMinus=new TGraph();
 708   gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
 709   gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
 710   TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
 711   TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
 712   TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
 713   TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
 714   //getting GF func for Kaons with TOF
 715   TF1 *fGFKPosTracking;
 716   fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
 717   TF1 *fGFKNegTracking;
 718   fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
 719   TF1 *fGFKPosMatching;
 720   fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
 721   TF1 *fGFKNegMatching;
 722   fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
 723   for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
 724     if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
 725       Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
 726       Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
 727       Printf("TPC Geant/Fluka: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
 728       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
 729       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
 730       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
 731       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
 732       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
 733       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
 734     }else{
 735       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
 736       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
 737       Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
 738       Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
 739       Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
 740       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
 741       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
 742       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
 743       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
 744       Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
 745       Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
 746       Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
 747       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
 748       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
 749       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
 750       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
 751     }
 752   }
 753
 754   //Geant Fluka for P in TPC
 755   Printf("\nGF correction for Protons");
 756   const Int_t kNCharge=2;
 757   Int_t kPos=0;
 758   Int_t kNeg=1;
 759   TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
 760   TH2D * hCorrFluka[kNCharge];
 761   TH2D * hCorrFluka[2];
 762   hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
 763   hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
 764   TGraph *gGFCorrectionProtonPlus=new TGraph();
 765   gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
 766   gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
 767   TGraph *gGFCorrectionProtonMinus=new TGraph();
 768   gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
 769   gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
 770   //getting GF func for Protons with TPCTOF
 771   TF1 *fGFpPosTracking;
 772   fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
 773   TF1 *fGFpNegTracking;
 774   fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
 775   TF1 *fGFpPosMatching;
 776   fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
 777   TF1 *fGFpNegMatching;
 778   fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
 779
 780   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 781     Int_t nbins = Spectra[2]->GetNbinsX();
 782     Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
 783     for(Int_t ibin = 0; ibin < nbins; ibin++){
 784       if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
 785         Float_t pt = Spectra[2]->GetBinCenter(ibin);
 786         Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
 787         Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 788         if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 789         if (pt > maxPtCorrection) pt = maxPtCorrection;
 790         Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
 791         //cout<<correction<<"     charge "<<icharge<<endl;
 792         if(icharge==0){
 793           if (correction != 0) {// If the bin is empty this is a  0
 794             Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
 795             Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError  (ibin)*correction);
 796             gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
 797           }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 798             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
 799             cout << " Bin content: " << Spectra[2]->GetBinContent(ibin)  << endl;
 800           }
 801         }
 802         if(icharge==1){
 803           if (correction != 0) {// If the bin is empty this is a  0
 804             Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
 805             Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError  (ibin)*correction);
 806             gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
 807           }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 808             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
 809             cout << " Bin content: " << Spectra[5]->GetBinContent(ibin)  << endl;
 810           }
 811         }
 812       }else{
 813         gGFCorrectionProtonPlus->SetPoint(ibin,Spectra[2]->GetBinCenter(ibin),0);
 814         gGFCorrectionProtonMinus->SetPoint(ibin,Spectra[5]->GetBinCenter(ibin),0);
 815         Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
 816         Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
 817         Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
 818         if(icharge==0){
 819           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
 820           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
 821         }else{
 822           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
 823           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
 824         }
 825         Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
 826         Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
 827         Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
 828         if(icharge==0){
 829           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
 830           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
 831         }else{
 832           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
 833           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
 834         }
 835       }
 836     }//end loop on bins
 837   }
 838   gGFCorrectionKaonPlus->SetLineColor(kRed);
 839   gGFCorrectionKaonMinus->SetLineColor(kRed+2);
 840   gGFCorrectionProtonPlus->SetLineColor(kGreen);
 841   gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
 842   fGFKPosTracking->SetLineColor(kRed);
 843   fGFKNegTracking->SetLineColor(kRed+2);
 844   fGFKPosMatching->SetLineColor(kRed);
 845   fGFKNegMatching->SetLineColor(kRed+2);
 846   fGFpPosTracking->SetLineColor(kGreen);
 847   fGFpNegTracking->SetLineColor(kGreen+2);
 848   fGFpPosMatching->SetLineColor(kGreen);
 849   fGFpNegMatching->SetLineColor(kGreen+2);
 850   fGFKPosTracking->SetLineStyle(2);
 851   fGFKNegTracking->SetLineStyle(2);
 852   fGFKPosMatching->SetLineStyle(3);
 853   fGFKNegMatching->SetLineStyle(3);
 854   fGFpPosTracking->SetLineStyle(2);
 855   fGFpNegTracking->SetLineStyle(2);
 856   fGFpPosMatching->SetLineStyle(3);
 857   fGFpNegMatching->SetLineStyle(3);
 858   fGFKPosTracking->SetRange(.6,5);
 859   fGFKNegTracking->SetRange(.6,5);
 860   fGFKPosMatching->SetRange(.6,5);
 861   fGFKNegMatching->SetRange(.6,5);
 862   fGFpPosTracking->SetRange(.6,5);
 863   fGFpNegTracking->SetRange(.6,5);
 864   fGFpPosMatching->SetRange(.6,5);
 865   fGFpNegMatching->SetRange(.6,5);
 866
 867   TCanvas * cGFCorrection = new TCanvas ("GFCorrection","GFCorrection",700,500);
 868   gPad->SetGridx();
 869   gPad->SetGridy();
 870   gGFCorrectionKaonPlus->DrawClone("al");
 871   gGFCorrectionKaonMinus->DrawClone("lsame");
 872   gGFCorrectionProtonPlus->DrawClone("lsame");
 873   gGFCorrectionProtonMinus->DrawClone("lsame");
 874   fGFKPosTracking->DrawClone("lsame");
 875   fGFKNegTracking->DrawClone("lsame");
 876   fGFKPosMatching->DrawClone("lsame");
 877   fGFKNegMatching->DrawClone("lsame");
 878   fGFpPosTracking->DrawClone("lsame");
 879   fGFpNegTracking->DrawClone("lsame");
 880   fGFpPosMatching->DrawClone("lsame");
 881   fGFpNegMatching->DrawClone("lsame");
 882   gPad->BuildLegend()->SetFillColor(0);
 883   fout->cd();
 884   cGFCorrection->Write();
 885 }
 886
 887
 888
 889 ///////////
 890 TF1 *
 891 TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
 892 {
 893
 894   if (ipart == 3 && icharge == kNegative) {
 895     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
 896     return f;
 897   }
 898   else if (ipart == 4 && icharge == kNegative) {
 899     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
 900   }
 901   else
 902     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
 903
 904   return f;
 905 }
 906
 907 Double_t
 908 TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
 909 {
 910   return 1.;
 911 }
 912
 913 Double_t
 914 TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
 915 {
 916   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
 917 }
 918
 919 Double_t
 920 TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
 921 {
 922   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
 923 }
 924 ///////////////////////////////////////////
 925 TF1 *
 926 TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
 927 {
 928
 929   if (ipart == 3 && icharge == kNegative) {
 930     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
 931     return f;
 932   }
 933   else if (ipart == 4 && icharge == kNegative) {
 934     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
 935   }
 936   else
 937     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
 938
 939   return f;
 940 }
 941
 942
 943 Double_t
 944 MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
 945 {
 946   return 1.;
 947 }
 948
 949 Double_t
 950 MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
 951 {
 952   Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
 953   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
 954 }
 955
 956 Double_t
 957 MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
 958 {
 959   Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
 960   return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
 961 }