TENDER becomes Tender
[u/mrichter/AliRoot.git] / 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 }