PWGLF/SPECTRA/PiKaPr/TestAOD/MainAnalysis.C

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