TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / HighLevelQA / JCGAnalysis.C
1 /////////////////////////////////////////////////////
2 // JCGAnalysis.C                                   //
3 // based on MainAnalysis.C                         //
4 // uses a new method of calculating the efficiency //
5 /////////////////////////////////////////////////////
6
7 class AliPID;
8
9 // global constants
10 const Int_t nCharge = 2;
11 const Int_t nPart = 3;
12 TString Charge[nCharge] = {"Pos","Neg"};
13 TString Sign[nCharge] = {"Plus","Minus"};
14 TString Particle[nPart] = {"Pion","Kaon","Proton"};
15 Int_t Color[nPart] = {1,2,4};
16 Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
17 Double_t Range[nPart] = {0.3,0.3,0.5}; // LowPt range for pi k p
18 enum ECharge_t {kPositive,kNegative,kNCharges};
19 TString Names[nCharge*nPart] = {"#pi^{+}",
20                                 "K^{+}",
21                                 "p",
22                                 "#pi^{-}",
23                                 "K^{-}",
24                                 "#bar{p}"};
25
26 void JCGAnalysis()
27 {
28   // load libraries 'n such
29   gSystem->Load("libCore.so");  
30   gSystem->Load("libPhysics.so");
31   gSystem->Load("libTree");
32   gSystem->Load("libMatrix");
33   gSystem->Load("libSTEERBase");
34   gSystem->Load("libESD");
35   gSystem->Load("libAOD");
36   gSystem->Load("libANALYSIS");
37   gSystem->Load("libOADB");
38   gSystem->Load("libANALYSISalice");
39   gSystem->Load("libTender");
40   gSystem->Load("libCORRFW");
41   gSystem->Load("libPWGTools");
42   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
43   gStyle->SetOptStat(000000);
44   gStyle->SetPalette(1);
45
46   // Set masses
47   Double_t mass[3];
48   mass[0]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
49   mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
50   mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
51   
52   // cuts 'n such
53   //                            0    1    2    3
54   Double_t CentCutMin[4]= {     0,  30,  30,  30};
55   Double_t CentCutMax[4]= {     5,  40,  40,  40};
56   Double_t QvecCutMin[4]={      0,   0,   0, 1.5};
57   Double_t QvecCutMax[4]={    100, 100, 0.4, 100};
58   Double_t EtaMin[4]={       -0.8,-0.8,-0.8,-0.8};
59   Double_t EtaMax[4]={        0.8, 0.8, 0.8, 0.8};
60   Double_t Nsigmapid=3.;
61   Double_t pt=10.;
62   Double_t p=10.;
63   Double_t y=.5;
64   Double_t ptTofMatch=.6;
65   UInt_t trkbit=1024;
66   UInt_t trkbitQVector=1;
67   Bool_t UseCentPatchAOD049=kFALSE;
68   Double_t DCA=100000;
69   UInt_t minNclsTPC=70;
70   Int_t nrebin=0;
71   TString opt="";
72   Int_t icut=1; // chooses which cut to use
73   if(icut==0)Int_t ibinToCompare=0;
74   else Int_t ibinToCompare=4;
75   
76   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);
77   TString sname_data_1sig=Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
78   // for MC we do not use the cut on Qvector
79   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);
80   TString sname_mc_5sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_5.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
81   TString sname_mc_1sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
82
83   // input and output files
84   TString fold="No_Outliers";
85   TString dataFile = Form("output/%s/AnResDATATrain12_NoOutliers.root",fold.Data());
86   TString mcFile = Form("output/%s/AnResMCTrain27_NoOutliers.root",fold.Data());
87   TFile * fout=new TFile(Form("results/%s/Res_%s_%s.root",fold.Data(),sname_data.Data(),fold.Data()),"RECREATE");
88
89   // Open root MC file and get classes
90   TFile *f_mc = TFile::Open(mcFile.Data());
91   TDirectoryFile *Dir_mc=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc.Data()));
92   AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) Dir_mc->Get("SpectraHistos");
93   AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) Dir_mc->Get("Event Cuts");
94   AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) Dir_mc->Get("Track Cuts");
95   TDirectoryFile *Dir_mc_5sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_5sig.Data()));
96   AliSpectraAODHistoManager* hman_mc_5sig = (AliSpectraAODHistoManager*) Dir_mc_5sig->Get("SpectraHistos");
97   AliSpectraAODEventCuts* ecuts_mc_5sig = (AliSpectraAODEventCuts*) Dir_mc_5sig->Get("Event Cuts");
98   AliSpectraAODTrackCuts* tcuts_mc_5sig = (AliSpectraAODTrackCuts*) Dir_mc_5sig->Get("Track Cuts");
99   TDirectoryFile *Dir_mc_1sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_1sig.Data()));
100   AliSpectraAODHistoManager* hman_mc_1sig = (AliSpectraAODHistoManager*) Dir_mc_1sig->Get("SpectraHistos");
101   AliSpectraAODEventCuts* ecuts_mc_1sig = (AliSpectraAODEventCuts*) Dir_mc_1sig->Get("Event Cuts");
102   AliSpectraAODTrackCuts* tcuts_mc_1sig = (AliSpectraAODTrackCuts*) Dir_mc_1sig->Get("Track Cuts");
103   // proceed likewise for data:
104   // Open root DATA file and get classes
105   TFile *f_data = TFile::Open(dataFile.Data());
106   TDirectoryFile *Dir_data=(TDirectoryFile*)f_data->Get(Form("%s",sname_data.Data()));
107   AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) Dir_data->Get("SpectraHistos");
108   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) Dir_data->Get("Event Cuts");
109   AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) Dir_data->Get("Track Cuts");
110   TDirectoryFile *Dir_data_1sig=(TDirectoryFile*)f_data->Get(Form("%s",sname_data_1sig.Data()));
111   AliSpectraAODHistoManager* hman_data_1sig = (AliSpectraAODHistoManager*) Dir_data_1sig->Get("SpectraHistos");
112   AliSpectraAODEventCuts* ecuts_data_1sig = (AliSpectraAODEventCuts*) Dir_data_1sig->Get("Event Cuts");
113   AliSpectraAODTrackCuts* tcuts_data_1sig = (AliSpectraAODTrackCuts*) Dir_data_1sig->Get("Track Cuts");
114
115   //---------------------------------------------
116   // plot nsigma distributions for presentation -
117   //---------------------------------------------
118   // TPC
119   TH2F * hTPCnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPC"))->Clone());
120   hTPCnsigPion->SetTitle("TPC NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
121   TCanvas * cTPCnsigPion = new TCanvas("cTPCnsigPion","cTPCnsigPion");
122   gPad->SetLogz();
123   hTPCnsigPion->GetXaxis()->SetRangeUser(0,5);
124   hTPCnsigPion->GetYaxis()->SetRangeUser(-10,40);
125   hTPCnsigPion->DrawCopy("COLZ");
126   // TOF
127   TH2F * hTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTOF"))->Clone());
128   hTOFnsigPion->SetTitle("TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
129   TCanvas * cTOFnsigPion = new TCanvas("cTOFnsigPion","cTOFnsigPion");
130   gPad->SetLogz();
131   hTOFnsigPion->GetYaxis()->SetRangeUser(-10,40);
132   hTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
133   hTOFnsigPion->DrawCopy("COLZ");
134   // TPCTOF
135   TH2F * hTPCTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPCTOF"))->Clone());
136   hTPCTOFnsigPion->SetTitle("TPC+TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
137   TCanvas * cTPCTOFnsigPion = new TCanvas("cTPCTOFnsigPion","cTPCTOFnsigPion");
138   gPad->SetLogz();
139   hTPCTOFnsigPion->GetYaxis()->SetRangeUser(-5,40);
140   hTPCTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
141   hTPCTOFnsigPion->DrawCopy("COLZ");
142
143
144   //------------------------------------------------------
145   // plot the MC True Generated Primaries for comparison -
146   //------------------------------------------------------
147   TH1F * MCTruthSpectra[nCharge*nPart];
148   Double_t events_mc = ecuts_mc->NumberOfEvents();
149   for (Int_t icharge=0; icharge<nCharge; icharge++)
150     {
151       for (Int_t ipart=0; ipart<nPart; ipart++)
152         {
153           Int_t index = ipart + 3*icharge;
154           TString hname = Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
155           MCTruthSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1)->Clone());
156           MCTruthSpectra[index]->SetTitle("Generated Primaries Spectra (MC Truth);p_{T} (GeV/c);");
157           MCTruthSpectra[index]->SetMarkerStyle(Marker[index]);
158           MCTruthSpectra[index]->SetMarkerColor(Color[ipart]);
159           MCTruthSpectra[index]->SetLineColor(Color[ipart]);
160           MCTruthSpectra[index]->Scale(1./events_mc,"width");
161         }
162     }
163   // neg/pos ratios
164   TH1F * NegPosRatio[nPart];
165   for (Int_t ipart=0; ipart<nPart; ipart++)
166     {
167       NegPosRatio[ipart] = (TH1F*)MCTruthSpectra[ipart+nPart]->Clone();
168       NegPosRatio[ipart]->Divide(MCTruthSpectra[ipart]);
169       NegPosRatio[ipart]->SetTitle("Neg/Pos;p_{T} (GeV/c);");
170     }
171   // drawing
172   TCanvas * cMCTruthSpectra = new TCanvas("cMCTruthSpectra","cMCTruthSpectra");
173   cMCTruthSpectra->Divide(1,2);
174   TLegend * lMCTruthSpectra = new TLegend(.69,.69,.89,.99);
175   lMCTruthSpectra->SetFillColor(0);
176   TLegend * lPosNegRatios = new TLegend (.69,.69,.99,.99);
177   lPosNegRatios->SetFillColor(0);
178   for (Int_t icharge=0; icharge<nCharge; icharge++)
179     {
180       for (Int_t ipart=0; ipart<nPart; ipart++)
181         {
182           Int_t index = ipart + 3*icharge;
183           cMCTruthSpectra->cd(1);
184           gPad->SetLogy();
185           if (index == 0) MCTruthSpectra[index]->DrawCopy("EP");
186           else MCTruthSpectra[index]->DrawCopy("EPsame");
187           lMCTruthSpectra->AddEntry(MCTruthSpectra[index],Names[index].Data(),"p");
188
189           cMCTruthSpectra->cd(2);
190           NegPosRatio[ipart]->GetYaxis()->SetRangeUser(0.4,1.1);
191           if (index == 0) NegPosRatio[ipart]->DrawCopy("hist][");
192           else NegPosRatio[ipart]->DrawCopy("hist][same");
193           if (icharge == 0) lPosNegRatios->AddEntry(NegPosRatio[ipart],Particle[ipart].Data(),"l");
194         }
195     }
196   cMCTruthSpectra->cd(1);
197   lMCTruthSpectra->DrawClone();
198   cMCTruthSpectra->cd(2);
199   lPosNegRatios->DrawClone();
200
201
202   //-----------------------------------
203   // Get the (normalized) raw spectra -
204   //-----------------------------------
205   TH1F * Spectra[nCharge*nPart];
206   Double_t events_data =  ecuts_data->NumberOfEvents();
207   Double_t events_data_1sig =  ecuts_data_1sig->NumberOfEvents();
208   Double_t events_mc =  ecuts_mc->NumberOfEvents();
209   Double_t events_mc_1sig = ecuts_mc_1sig->NumberOfEvents();
210   for (Int_t icharge=0; icharge<nCharge; icharge++)
211     {
212       for (Int_t ipart=0; ipart<nPart; ipart++)
213         {
214           Int_t index = ipart + 3*icharge;
215           TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
216           Spectra[index] = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
217           Spectra[index]->SetTitle(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
218           Spectra[index]->SetMarkerStyle(Marker[index]);
219           Spectra[index]->SetMarkerColor(Color[ipart]); 
220           Spectra[index]->SetLineColor(Color[ipart]);
221           Spectra[index]->Scale(1./events_data,"width"); // normalization
222         }
223     }
224   // Set bins outisde the range to contain zeros
225   for (Int_t icharge=0; icharge<nCharge; icharge++)
226     {
227       for (Int_t ipart=0; ipart<nPart; ipart++)
228         {
229           Int_t index = ipart + 3*icharge;
230           for (Int_t ibin=0; ibin<Spectra[index]->GetNbinsX(); ibin++)
231             {
232               if (Spectra[index]->GetBinCenter(ibin) < Range[ipart])
233                 {
234                   Spectra[index]->SetBinContent(ibin, 0);
235                   Spectra[index]->SetBinError(ibin, 0);
236                 }
237             }
238         }
239     }
240   // save the raw spectra for use later
241   TH1F * RawSpectra[nCharge*nPart];
242   for (Int_t icharge=0; icharge<nCharge; icharge++)
243     {
244       for (Int_t ipart=0; ipart<nPart; ipart++)
245         {
246           Int_t index = ipart + 3*icharge;
247           RawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
248           RawSpectra[index]->SetTitle("Raw Spectra");
249         }
250     }
251
252   // plot the raw spectra from DATA
253   // also include allch
254   TCanvas * cRawSpectra = new TCanvas("cRawSpectra","cRawSpectra");
255   TLegend * lRawSpectra = new TLegend(.69,.69,.99,.99);
256   lRawSpectra->SetFillColor(0);
257   TH1F * hRawAllCh = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
258   hRawAllCh->Scale(1./ecuts_data->NumberOfEvents(),"width");
259   hRawAllCh->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
260   hRawAllCh->SetMarkerColor(kGreen+3);
261   hRawAllCh->SetLineColor(kGreen+3);
262   hRawAllCh->SetMarkerStyle(34);
263   hRawAllCh->DrawCopy();
264   lRawSpectra->AddEntry(hRawAllCh,"All Particles","p");
265   for (Int_t icharge=0; icharge<nCharge; icharge++)
266     {
267       for (Int_t ipart=0; ipart<nPart; ipart++)
268         {
269           Int_t index = ipart + 3*icharge;
270           RawSpectra[index]->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
271           RawSpectra[index]->SetMarkerColor(Color[ipart]);
272           RawSpectra[index]->SetLineColor(Color[ipart]);
273           RawSpectra[index]->SetMarkerStyle(Marker[index]);
274           RawSpectra[index]->DrawCopy("same");
275           lRawSpectra->AddEntry(RawSpectra[index],Names[index].Data(),"p");
276         }
277     }
278   lRawSpectra->DrawClone();
279
280   // also get MC raw for comparison purposes
281   TH1F * MCSpectra[nCharge*nPart];
282   for (Int_t icharge=0; icharge<nCharge; icharge++)
283     {
284       for (Int_t ipart=0; ipart<nPart; ipart++)
285         {
286           Int_t index = ipart + 3*icharge;
287           TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
288           MCSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
289           MCSpectra[index]->SetTitle(Form("MCSpectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
290           MCSpectra[index]->SetMarkerStyle(Marker[index]);
291           MCSpectra[index]->SetMarkerColor(Color[ipart]);
292           MCSpectra[index]->SetLineColor(Color[ipart]);
293           MCSpectra[index]->Scale(1./events_mc,"width"); // normalization
294         }
295     }
296   // Set bins outisde the range to contain zeros
297   for (Int_t icharge=0; icharge<nCharge; icharge++)
298     {
299       for (Int_t ipart=0; ipart<nPart; ipart++)
300         {
301           Int_t index = ipart + 3*icharge;
302           for (Int_t ibin=0; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
303             {
304               if (MCSpectra[index]->GetBinCenter(ibin) < Range[ipart])
305                 {
306                   MCSpectra[index]->SetBinContent(ibin, 0);
307                   MCSpectra[index]->SetBinError(ibin, 0);
308                 }
309             }
310         }
311     }
312   // save the raw spectra for use later
313   TH1F * MCRawSpectra[nCharge*nPart];
314   for (Int_t icharge=0; icharge<nCharge; icharge++)
315     {
316       for (Int_t ipart=0; ipart<nPart; ipart++)
317         {
318           Int_t index = ipart + 3*icharge;
319           MCRawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
320           MCRawSpectra[index]->SetTitle("MC Raw Spectra");
321         }
322     }
323
324
325   //----------------
326   // MC Correction -
327   //----------------
328   // calculate the percentage of primaries reconstructed w/ 5 sigma PID
329   Printf("\n\n-> Calculating MC Correction Factors");
330   TH1F * CorrFact[nPart*nCharge];
331   for (Int_t icharge=0; icharge<nCharge; icharge++)
332     {
333       for (Int_t ipart=0; ipart<nPart; ipart++)
334         {
335           Int_t index = ipart + 3*icharge;
336           TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
337           Printf("Getting %s",hname.Data());
338           CorrFact[index] = new TH1F(*(TH1F*)((TH1F*) hman_mc_5sig->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
339           CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
340           CorrFact[index]->SetTitle("MC Correction Factor");
341           CorrFact[index]->SetMarkerStyle(Marker[index]);
342           CorrFact[index]->SetMarkerColor(Color[ipart]);
343           CorrFact[index]->SetLineColor(Color[ipart]);
344           hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
345           Printf("... and divide it by %s",hname.Data());
346           CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc_5sig->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
347           if (index == 0)
348             {
349               TCanvas * cCorrFactMC = new TCanvas("cCorrFactMC","cCorrFactMC");
350               TLegend * lCorrFactMC = new TLegend(.69,.69,.99,.99);
351               lCorrFactMC->SetFillColor(0);
352               CorrFact[index]->DrawCopy("EP");
353             }
354           else CorrFact[index]->DrawCopy("EPsame");
355           lCorrFactMC->AddEntry(CorrFact[index],Names[index].Data(),"p");
356         } // end loop over ipart
357     } // end loop over icharge
358   lCorrFactMC->DrawClone();
359   // Correct spectra with MC correction factor
360   for (Int_t icharge=0; icharge<nCharge; icharge++)
361     {
362       for (Int_t ipart=0; ipart<nPart; ipart++)
363         {
364           Int_t index = ipart + 3*icharge;
365           Spectra[index]->Divide(CorrFact[index]);
366         }
367     }
368   // save this version of the spectra for use later
369   TH1F * SpectraStep1[nCharge*nPart];
370   for (Int_t icharge=0; icharge<nCharge; icharge++)
371     {
372       for (Int_t ipart=0; ipart<nPart; ipart++)
373         {
374           Int_t index = ipart + 3*icharge;
375           SpectraStep1[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
376           SpectraStep1[index]->SetTitle("Spectra After MC Correction");
377         }
378     }
379   // do the same with MCSpectra
380   for (Int_t icharge=0; icharge<nCharge; icharge++)
381     {
382       for (Int_t ipart=0; ipart<nPart; ipart++)
383         {
384           Int_t index = ipart + 3*icharge;
385           MCSpectra[index]->Divide(CorrFact[index]);
386         }
387     }
388   // save this version of the spectra for use later
389   TH1F * MCSpectraStep1[nCharge*nPart];
390   for (Int_t icharge=0; icharge<nCharge; icharge++)
391     {
392       for (Int_t ipart=0; ipart<nPart; ipart++)
393         {
394           Int_t index = ipart + 3*icharge;
395           MCSpectraStep1[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
396           MCSpectraStep1[index]->SetTitle("MCSpectra After MC Correction");
397         }
398     }
399
400
401   //--------------------------------------
402   // Correct spectra for PID nsigma cuts -
403   //--------------------------------------
404   // for the TPC only, assume the nsigma distribution is purely gaussian, but account for the actual means and sigmas
405   // the correction factor is the integral of a normal distribution from -3sqrt(2) to 3sqrt(2)
406   TF1 * fGaussTPC[nPart];
407   for (Int_t ipart=0; ipart<nPart; ipart++)
408     {
409       fGaussTPC[ipart] = new TF1(Form("fGaussTPC%s",Particle[ipart].Data()),"gaus",-25,25);
410       fGaussTPC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
411     }
412   // pions
413   fGaussTPC[0]->SetParameter(1,0.04);
414   fGaussTPC[0]->SetParameter(2,1.02);
415   // kaons
416   fGaussTPC[1]->SetParameter(1,-0.06);
417   fGaussTPC[1]->SetParameter(2,0.97);
418   // protons
419   fGaussTPC[2]->SetParameter(1,0.25);
420   fGaussTPC[2]->SetParameter(2,0.72);
421   // normalization
422   for (Int_t ipart=0; ipart<nPart; ipart++)
423     {
424       fGaussTPC[ipart]->SetParameter(0, 1./fGaussTPC[ipart]->Integral(-20,20));
425     }
426   // calculate correction factors
427   Float_t PIDCorrFactTPC[nPart];
428   for (Int_t ipart=0; ipart<nPart; ipart++)
429     {
430       PIDCorrFactTPC[ipart] = fGaussTPC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
431     }
432   cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS: ---\n"
433        << "Pions: " << PIDCorrFactTPC[0] << endl
434        << "Kaons: " << PIDCorrFactTPC[1] << endl
435        << "Protons: " << PIDCorrFactTPC[2] << endl << endl << endl;
436   // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
437   TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
438   TH1F * hTPCTOFnsigma[nPart];
439   Float_t PIDCorrFactTPCTOF[nPart];
440   for (Int_t ipart=0; ipart<nPart; ipart++)
441     {
442       hTPCTOFnsigma[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredicted"))->FindObject(Form("hTPCTOFnsigma%s",Particle[ipart].Data())))->Clone();
443       PIDCorrFactTPCTOF[ipart] = hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(3));
444       PIDCorrFactTPCTOF[ipart] /= hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(20));
445     }
446   cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS: ---\n"
447        << "Pions: " << PIDCorrFactTPCTOF[0] << endl
448        << "Kaons: " << PIDCorrFactTPCTOF[1] << endl
449        << "Protons: " << PIDCorrFactTPCTOF[2] << endl << endl << endl;
450
451   // calculate seperate correction factors for MC b/c the parameters will be different...
452   TF1 * fGaussTPCMC[nPart];
453   for (Int_t ipart=0; ipart<nPart; ipart++)
454     {
455       fGaussTPCMC[ipart] = new TF1(Form("fGaussTPCMC%s",Particle[ipart].Data()),"gaus",-25,25);
456       fGaussTPCMC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
457     }
458   // pions
459   fGaussTPCMC[0]->SetParameter(1,-0.03);
460   fGaussTPCMC[0]->SetParameter(2,1.35);
461   // kaons
462   fGaussTPCMC[1]->SetParameter(1,0.14);
463   fGaussTPCMC[1]->SetParameter(2,1.27);
464   // protons
465   fGaussTPCMC[2]->SetParameter(1,-0.17);
466   fGaussTPCMC[2]->SetParameter(2,1.03);
467   // normalization
468   for (Int_t ipart=0; ipart<nPart; ipart++)
469     {
470       fGaussTPCMC[ipart]->SetParameter(0, 1./fGaussTPCMC[ipart]->Integral(-20,20));
471     }
472   // calculate correction factors
473   Float_t PIDCorrFactTPCMC[nPart];
474   for (Int_t ipart=0; ipart<nPart; ipart++)
475     {
476       PIDCorrFactTPCMC[ipart] = fGaussTPCMC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
477     }
478   cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS (MC): ---\n"
479        << "Pions: " << PIDCorrFactTPCMC[0] << endl
480        << "Kaons: " << PIDCorrFactTPCMC[1] << endl
481        << "Protons: " << PIDCorrFactTPCMC[2] << endl << endl << endl;
482   // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
483   TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
484   TH1F * hTPCTOFnsigmaMC[nPart];
485   Float_t PIDCorrFactTPCTOFMC[nPart];
486   for (Int_t ipart=0; ipart<nPart; ipart++)
487     {
488       hTPCTOFnsigmaMC[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredictedMC"))->FindObject(Form("hTPCTOFnsigmaMC%s",Particle[ipart].Data())))->Clone();
489       PIDCorrFactTPCTOFMC[ipart] = hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(3));
490       PIDCorrFactTPCTOFMC[ipart] /= hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(20));
491     }
492   cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS (MC): ---\n"
493        << "Pions: " << PIDCorrFactTPCTOFMC[0] << endl
494        << "Kaons: " << PIDCorrFactTPCTOFMC[1] << endl
495        << "Protons: " << PIDCorrFactTPCTOFMC[2] << endl << endl << endl;
496
497
498   // correct spectra, taking into account the range where the TOF is valid
499   for (Int_t icharge=0; icharge<nCharge; icharge++)
500     {
501       for (Int_t ipart=0; ipart<nPart; ipart++)
502         {
503           Int_t index = ipart + 3*icharge;
504           for (Int_t ibin=1; ibin<Spectra[index]->GetNbinsX(); ibin++)
505             {
506               if (ibin < Spectra[index]->FindBin(0.6))
507                 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPC[ipart]);                  
508               else
509                 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOF[ipart]);
510             }  
511         } // end for loop over ipart
512     } // end for loop over icharge
513   // save this version of the spectra for use later
514   TH1F * SpectraStep2[nCharge*nPart];
515   for (Int_t icharge=0; icharge<nCharge; icharge++)
516     {
517       for (Int_t ipart=0; ipart<nPart; ipart++)
518         {
519           Int_t index = ipart + 3*icharge;
520           SpectraStep2[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
521           SpectraStep2[index]->SetTitle("Spectra After PID Correction");
522         }
523     }
524   // also do this for MCSpectra
525   for (Int_t icharge=0; icharge<nCharge; icharge++)
526     {
527       for (Int_t ipart=0; ipart<nPart; ipart++)
528         {
529           Int_t index = ipart + 3*icharge;
530           for (Int_t ibin=1; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
531             {
532               if (ibin < MCSpectra[index]->FindBin(0.6))
533                 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCMC[ipart]);            
534               else
535                 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOFMC[ipart]);
536             }  
537         } // end for loop over ipart
538     } // end for loop over icharge
539   // save this version of the spectra for use later
540   TH1F * MCSpectraStep2[nCharge*nPart];
541   for (Int_t icharge=0; icharge<nCharge; icharge++)
542     {
543       for (Int_t ipart=0; ipart<nPart; ipart++)
544         {
545           Int_t index = ipart + 3*icharge;
546           MCSpectraStep2[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
547           MCSpectraStep2[index]->SetTitle("MCSpectra After PID Correction");
548         }
549     }
550
551   
552   //-------------------------------------------------------------------------
553   // Correct spectra for contamination from other Pions, Kaons, and Protons -
554   //-------------------------------------------------------------------------
555   // calculate the correction factors first
556   TH1F * ContCorrFact[nPart*nCharge];
557   for (Int_t icharge=0; icharge<nCharge; icharge++)
558     {
559       for (Int_t ipart=0; ipart<nPart; ipart++)
560         {
561           Int_t index = ipart + 3*icharge;
562           TString hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
563           ContCorrFact[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
564           ContCorrFact[index]->SetTitle("Contamination Correction Factor");
565           ContCorrFact[index]->SetMarkerStyle(Marker[index]);
566           ContCorrFact[index]->SetMarkerColor(Color[ipart]);
567           ContCorrFact[index]->SetLineColor(Color[ipart]);
568           hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
569           ContCorrFact[index]->Divide(ContCorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(),-1,-1,"B"); // binomial error
570           if (index == 0)
571             {
572               TCanvas * cContCorrFact = new TCanvas("cContCorrFact","cContCorrFact");
573               TLegend * lContCorrFact = new TLegend(.69,.69,.99,.99);
574               ContCorrFact[index]->DrawCopy("EP");
575             }
576           else ContCorrFact[index]->DrawCopy("EPsame");
577           lContCorrFact->AddEntry(ContCorrFact[index],Names[index].Data(),"p");
578         } // end for loop over ipart
579     } // end for loop over icharge
580   lContCorrFact->DrawClone();
581   // then use them to correct the spectra
582   for (Int_t icharge=0; icharge<nCharge; icharge++)
583     {
584       for(Int_t ipart=0; ipart<nPart; ipart++)
585         {
586           Int_t index = ipart + 3*icharge;
587           Spectra[index]->Multiply(ContCorrFact[index]);
588         }
589     }
590   // save this version of the spectra for use later
591   TH1F * SpectraStep3[nCharge*nPart];
592   for (Int_t icharge=0; icharge<nCharge; icharge++)
593     {
594       for (Int_t ipart=0; ipart<nPart; ipart++)
595         {
596           Int_t index = ipart + 3*icharge;
597           SpectraStep3[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
598           SpectraStep3[index]->SetTitle("Spectra After Contamination Correction");
599         }
600     }
601   // also do this for MCSpectra
602   for (Int_t icharge=0; icharge<nCharge; icharge++)
603     {
604       for(Int_t ipart=0; ipart<nPart; ipart++)
605         {
606           Int_t index = ipart + 3*icharge;
607           MCSpectra[index]->Multiply(ContCorrFact[index]);
608         }
609     }
610   // save this version of the spectra for use later
611   TH1F * MCSpectraStep3[nCharge*nPart];
612   for (Int_t icharge=0; icharge<nCharge; icharge++)
613     {
614       for (Int_t ipart=0; ipart<nPart; ipart++)
615         {
616           Int_t index = ipart + 3*icharge;
617           MCSpectraStep3[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
618           MCSpectraStep3[index]->SetTitle("MCSpectra After Contamination Correction");
619         }
620     }
621
622
623   //--------------------------------------------------
624   // Correct spectra for the Geant-Fluka discrepancy -
625   //--------------------------------------------------
626   GFCorrection(Spectra,tcuts_data,fout);
627   // save this version of the spectra for use later
628   TH1F * SpectraStep4[nCharge*nPart];
629   for (Int_t icharge=0; icharge<nCharge; icharge++)
630     {
631       for (Int_t ipart=0; ipart<nPart; ipart++)
632         {
633           Int_t index = ipart + 3*icharge;
634           SpectraStep4[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
635           SpectraStep4[index]->SetTitle("Spectra After Geant-Fluka Correction");
636         }
637     }
638   // NOTE! DO NOT USE GF FOR MCSPECTRA!
639
640
641   //--------------------------------------------------
642   // Correct spectra for the TOF matching efficiency -
643   //--------------------------------------------------
644   //Matching efficiency in data and Monte Carlo
645   TCanvas *cMatchingEff=new TCanvas("cMatchingEff","cMatchingEff",700,500);
646   TH1F *hMatcEffPos_data=(TH1F*)tcuts_data->GetHistoNMatchedPos();
647   hMatcEffPos_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedPos());
648   hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
649   TH1F *hMatcEffNeg_data=(TH1F*)tcuts_data->GetHistoNMatchedNeg();
650   hMatcEffNeg_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedNeg());
651   hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
652   hMatcEffNeg_data->SetLineColor(2);
653   /*  
654   // for my (John Groh) report, plot the matching efficiency just for data
655   TCanvas * cMatchEffDataOnly = new TCanvas("cMatchEffDataOnly","cMatchEffDataOnly");
656   hMatcEffPos_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
657   hMatcEffPos_data->DrawCopy("lhist");
658   hMatcEffNeg_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
659   hMatcEffNeg_data->DrawCopy("lhistsame");
660   TLegend * lMatcEffDataOnly = new TLegend(.69,.69,.99,.99);
661   lMatcEffDataOnly->SetFillColor(0);
662   lMatcEffDataOnly->AddEntry(hMatcEffPos_data,"Positive Particles","l");
663   lMatcEffDataOnly->AddEntry(hMatcEffNeg_data,"Negative Particles","l");
664   lMatcEffDataOnly->DrawClone();
665   */
666
667   TH1F *hMatcEffPos_mc=(TH1F*)tcuts_mc->GetHistoNMatchedPos();
668   hMatcEffPos_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedPos());
669   hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
670   hMatcEffPos_mc->SetLineStyle(2);
671   TH1F *hMatcEffNeg_mc=(TH1F*)tcuts_mc->GetHistoNMatchedNeg();
672   hMatcEffNeg_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedNeg());
673   hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
674   hMatcEffNeg_mc->SetLineColor(2);
675   hMatcEffNeg_mc->SetLineStyle(2);
676   cMatchingEff->Divide(1,2);
677   cMatchingEff->cd(1);
678   gPad->SetGridy();
679   gPad->SetGridx();
680   hMatcEffPos_data->DrawClone("lhist");
681   hMatcEffNeg_data->DrawClone("lhistsame");
682   hMatcEffPos_mc->DrawClone("lhistsame");
683   hMatcEffNeg_mc->DrawClone("lhistsame");
684   gPad->BuildLegend()->SetFillColor(0);
685   hMatcEffPos_data->Divide(hMatcEffPos_mc);
686   hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
687   cMatchingEff->cd(2);
688   gPad->SetGridy();
689   gPad->SetGridx();
690   hMatcEffPos_data->DrawClone("lhist");
691   hMatcEffNeg_data->DrawClone("lhistsame");
692   TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",2.5,5);
693   hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
694   pol0MatchPos_data->DrawClone("same");
695   TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",2.5,5);
696   hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
697   pol0MatchNeg_data->SetLineColor(2);
698   pol0MatchNeg_data->DrawClone("same");
699   Float_t ScalingMatchingPos=pol0MatchPos_data->GetParameter(0);
700   Float_t ScalingMatchingNeg=pol0MatchNeg_data->GetParameter(0);
701   //Correction spectra for matching efficiency
702   for (Int_t ipart=0; ipart<nPart; ipart++)
703     {
704       for (Int_t ibin=1; ibin<=Spectra[ipart]->GetNbinsX(); ibin++)
705         {
706           Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
707           if (ptspectra < tcuts_data->GetPtTOFMatching()) continue;
708           Spectra[ipart]->SetBinContent(ibin, (Spectra[ipart]->GetBinContent(ibin) / ScalingMatchingPos));
709           Spectra[ipart+3]->SetBinContent(ibin, (Spectra[ipart+3]->GetBinContent(ibin) / ScalingMatchingNeg));
710         }
711     }
712   // save this version of the spectra for use later
713   TH1F * SpectraStep5[nCharge*nPart];
714   for (Int_t icharge=0; icharge<nCharge; icharge++)
715     {
716       for (Int_t ipart=0; ipart<nPart; ipart++)
717         {
718           Int_t index = ipart + 3*icharge;
719           SpectraStep5[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
720           SpectraStep5[index]->SetTitle("Spectra After Matching Efficiency Correction");
721         }
722     }
723   // NOTE! DO NOT USE THE MATCHING EFFICIENCY CORRECTION FOR MCSPECTRA!
724
725
726   //-----------------------------------------------------
727   // Secondaries correction - use MC, not DCA from data -
728   //------------------------------------------------------
729   // (correction factor for secondaries) = (true primaries) / (primaries)
730   TH1F * SecondaryCorrFact[nCharge*nPart];
731   for (Int_t icharge=0; icharge<nCharge; icharge++)
732     {
733       for (Int_t ipart=0; ipart<nPart; ipart++)
734         {
735           Int_t index = ipart + 3*icharge;
736           cout << "index = " << index << endl;
737           TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
738           SecondaryCorrFact[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
739           hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
740           SecondaryCorrFact[index]->Divide(SecondaryCorrFact[index], (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(), 1, 1, "B"); // binomial error!
741           SecondaryCorrFact[index]->SetName(Form("SecondaryCorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
742           SecondaryCorrFact[index]->SetTitle("Correction Factor for Secondaries");
743           SecondaryCorrFact[index]->SetStats(kFALSE);
744           SecondaryCorrFact[index]->SetMarkerStyle(Marker[index]);
745           SecondaryCorrFact[index]->SetMarkerColor(Color[ipart]);
746           SecondaryCorrFact[index]->SetLineColor(Color[ipart]);
747           if (index == 0)
748             {
749               TCanvas * cSecCorrFactMC = new TCanvas("cSecCorrFactMC","cSecCorrFactMC");
750               TLegend * lSecCorrFactMC = new TLegend(.69,.69,.99,.99);
751               lSecCorrFactMC->SetFillColor(0);
752               SecondaryCorrFact[index]->DrawCopy("EP");           
753             }
754           else SecondaryCorrFact[index]->DrawCopy("EPsame");
755           lSecCorrFactMC->AddEntry(SecondaryCorrFact[index],Names[index],"ep");
756         } // end for loop over ipart
757     } // end for loop over icharge
758   lSecCorrFactMC->DrawClone();
759
760   // correct Spectra with this
761   for (Int_t icharge=0; icharge<nCharge; icharge++)
762     {
763       for (Int_t ipart=0; ipart<nPart; ipart++)
764         {
765           Int_t index = ipart + 3*icharge;
766           Spectra[index]->Multiply(SecondaryCorrFact[index]);
767         }
768     }
769
770   // also correct MCSpectra with this
771   for (Int_t icharge=0; icharge<nCharge; icharge++)
772     {
773       for (Int_t ipart=0; ipart<nPart; ipart++)
774         {
775           Int_t index = ipart + 3*icharge;
776           MCSpectra[index]->Multiply(SecondaryCorrFact[index]);
777         }
778     }
779   // save this version of the spectra for use later
780   TH1F * SpectraStep6[nCharge*nPart];
781   for (Int_t icharge=0; icharge<nCharge; icharge++)
782     {
783       for (Int_t ipart=0; ipart<nPart; ipart++)
784         {
785           Int_t index = ipart + 3*icharge;
786           SpectraStep6[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
787           SpectraStep6[index]->SetTitle("Spectra After Secondaries Correction");
788         }
789     }
790   // also save this version of MCSpectra
791   TH1F * MCSpectraStep6[nCharge*nPart];
792   for (Int_t icharge=0; icharge<nCharge; icharge++)
793     {
794       for (Int_t ipart=0; ipart<nPart; ipart++)
795         {
796           Int_t index = ipart + 3*icharge;
797           MCSpectraStep6[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
798           MCSpectraStep6[index]->SetTitle("MCSpectra After Secondaries Correction");
799         }
800     }
801   
802   //-------------------------
803   // draw the final spectra -
804   //-------------------------
805   TCanvas * cSpectra = new TCanvas("cSpectra","cSpectra",100,100,700,500);
806   gPad->SetGridy();
807   gPad->SetLogy();
808   TLegend * lSpectra = new TLegend(.69,.69,.99,.99);
809   for (Int_t icharge=0; icharge<nCharge; icharge++)
810     {
811       for (Int_t ipart=0; ipart<nPart; ipart++)
812         {
813           Int_t index = ipart + 3*icharge;
814           Spectra[index]->SetTitle("Final Spectra");
815           if (index == 0) Spectra[index]->DrawClone("EP");
816           else Spectra[index]->DrawClone("EPsame");
817           lSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
818         }
819     }
820   lSpectra->DrawClone();
821
822   //-----------------------------------------------
823   // compare MCSpectra with hHistPtGenTruePrimary -
824   //-----------------------------------------------
825   TCanvas * cCompareMCRecWithMCTrue = new TCanvas("cCompareMCRecWithMCTrue","cCompareMCRecWithMCTrue",250,175,700,500);
826   cCompareMCRecWithMCTrue->Divide(3,2);
827   TH1F * hMCTruth[nPart*nCharge];
828   for (Int_t icharge=0; icharge<nCharge; icharge++)
829     {
830       for (Int_t ipart=0; ipart<nPart; ipart++)
831         {
832           Int_t index = ipart + 3*icharge;
833           hMCTruth[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),1,1))->Clone());
834           hMCTruth[index]->Scale(1./ecuts_mc->NumberOfEvents(),"width"); // normalize it first!
835           hMCTruth[index]->SetMarkerStyle(Marker[index]);
836           hMCTruth[index]->SetMarkerColor(Color[ipart]);
837           hMCTruth[index]->SetLineColor(Color[ipart]);
838         }
839     }
840   // set bins outside range to zero
841   for (Int_t icharge=0; icharge<nCharge; icharge++)
842     {
843       for (Int_t ipart=0; ipart<nPart; ipart++)
844         {
845           Int_t index = ipart + 3*icharge;
846           for (Int_t ibin=0; ibin<hMCTruth[index]->GetNbinsX(); ibin++)
847             {
848               if (hMCTruth[index]->GetBinCenter(ibin) < Range[ipart])
849                 {
850                   hMCTruth[index]->SetBinContent(ibin, 0);
851                   hMCTruth[index]->SetBinError(ibin, 0);
852                 }
853             }
854         }
855     }
856
857
858   TH1F * hRatioMCRecTrue[nPart*nCharge];  
859   for (Int_t ipart=0; ipart<nPart; ipart++)
860     {
861       for (Int_t icharge=0; icharge<nCharge; icharge++)
862         {
863           Int_t index = ipart + 3*icharge;        
864           hRatioMCRecTrue[index] = (TH1F*)MCSpectra[index]->Clone();
865           hRatioMCRecTrue[index]->Divide(hMCTruth[index]);
866           hRatioMCRecTrue[index]->SetTitle("(Corrected MC Spectra) / (MC Truth Spectra);p_{T} (GeV/c);");
867           if (icharge == 1) hRatioMCRecTrue[index]->SetLineStyle(2);
868         }
869     }
870   for (Int_t ipart=0; ipart<nPart; ipart++)
871     {
872       cCompareMCRecWithMCTrue->cd(ipart+1);
873       gPad->SetLogy();
874       MCSpectra[ipart]->SetMarkerColor(kGreen+3);
875       MCSpectra[ipart+nPart]->SetMarkerColor(kGreen+2);
876       MCSpectra[ipart]->SetMarkerSize(0.75);
877       MCSpectra[ipart+nPart]->SetMarkerSize(0.75);
878       MCSpectra[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
879       MCSpectra[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
880       hMCTruth[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
881       hMCTruth[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
882       hMCTruth[ipart]->SetMarkerSize(0.75);
883       hMCTruth[ipart]->SetMarkerSize(0.75);
884       hMCTruth[ipart]->DrawCopy("EP");
885       hMCTruth[ipart+nPart]->DrawCopy("EPsame");
886       MCSpectra[ipart]->DrawCopy("EPsame");
887       MCSpectra[ipart+nPart]->DrawCopy("EPsame");
888       TLegend * lCompare = new TLegend(.69,.59,.99,.91);
889       lCompare->SetFillColor(0);
890       lCompare->AddEntry(hMCTruth[ipart],Form("%s, Truth",Names[ipart].Data()),"p");
891       lCompare->AddEntry(hMCTruth[ipart+nPart],Form("%s, Truth",Names[ipart+nPart].Data()),"p");
892       lCompare->AddEntry(MCSpectra[ipart],Form("%s, Corrected",Names[ipart].Data()),"p");
893       lCompare->AddEntry(MCSpectra[ipart+nPart],Form("%s, Corrected",Names[ipart+nPart].Data()),"p");
894       lCompare->DrawClone();
895
896       cCompareMCRecWithMCTrue->cd(ipart+1+nPart);
897       hRatioMCRecTrue[ipart]->GetYaxis()->SetRangeUser(.97,1.03);
898       hRatioMCRecTrue[ipart]->DrawCopy("hist][");
899       hRatioMCRecTrue[ipart+nPart]->DrawCopy("hist][same");
900       TLegend * lCompareRatio = new TLegend(.79,.79,.99,.99);
901       lCompareRatio->SetFillColor(0);
902       lCompareRatio->AddEntry(hRatioMCRecTrue[ipart],Names[ipart].Data(),"l");
903       lCompareRatio->AddEntry(hRatioMCRecTrue[ipart+nPart],Names[ipart+nPart].Data(),"l");
904       lCompareRatio->DrawClone();
905     }
906
907
908   //---------------------------------------------------------
909   // plot the evolution of the spectra afte each correction -
910   //---------------------------------------------------------
911   TCanvas * cEvolutionOfSpectra = new TCanvas("cEvolutionOfSpectra","cEvolutionOfSpectra",300,200,700,500);
912   cEvolutionOfSpectra->Divide(4,2);
913   TLegend * lEvolutionOfSpectra = new TLegend(.69,.69,.99,.99);
914   for (Int_t icharge=0; icharge<nCharge; icharge++)
915     {
916       for (Int_t ipart=0; ipart<nPart; ipart++)
917         {
918           Int_t index = ipart + 3*icharge;
919           lEvolutionOfSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
920         }
921     }
922   cEvolutionOfSpectra->cd(1);
923   gPad->SetLogy();
924   for (Int_t icharge=0; icharge<nCharge; icharge++)
925     {
926       for (Int_t ipart=0; ipart<nPart; ipart++)
927         {
928           Int_t index = ipart + 3*icharge;
929           RawSpectra[index]->SetMarkerSize(.5);
930           if (index == 0) RawSpectra[index]->DrawClone("EP");
931           else RawSpectra[index]->DrawClone("EPsame");
932         }
933     }
934   lEvolutionOfSpectra->DrawClone();
935   cEvolutionOfSpectra->cd(2);
936   gPad->SetLogy();
937   for (Int_t icharge=0; icharge<nCharge; icharge++)
938     {
939       for (Int_t ipart=0; ipart<nPart; ipart++)
940         {
941           Int_t index = ipart + 3*icharge;
942           SpectraStep1[index]->SetMarkerSize(.5);
943           if (index == 0) SpectraStep1[index]->DrawClone("EP");
944           else SpectraStep1[index]->DrawClone("EPsame");
945         }
946     }
947   lEvolutionOfSpectra->DrawClone();
948   cEvolutionOfSpectra->cd(3);
949   gPad->SetLogy();
950   for (Int_t icharge=0; icharge<nCharge; icharge++)
951     {
952       for (Int_t ipart=0; ipart<nPart; ipart++)
953         {
954           Int_t index = ipart + 3*icharge;
955           SpectraStep2[index]->SetMarkerSize(.5);
956           if (index == 0) SpectraStep2[index]->DrawClone("EP");
957           else SpectraStep2[index]->DrawClone("EPsame");
958         }
959     }
960   lEvolutionOfSpectra->DrawClone();
961   cEvolutionOfSpectra->cd(4);
962   gPad->SetLogy();
963   for (Int_t icharge=0; icharge<nCharge; icharge++)
964     {
965       for (Int_t ipart=0; ipart<nPart; ipart++)
966         {
967           Int_t index = ipart + 3*icharge;
968           SpectraStep3[index]->SetMarkerSize(.5);
969           if (index == 0) SpectraStep3[index]->DrawClone("EP");
970           else SpectraStep3[index]->DrawClone("EPsame");
971         }
972     }
973   lEvolutionOfSpectra->DrawClone();
974   cEvolutionOfSpectra->cd(5);
975   gPad->SetLogy();
976   for (Int_t icharge=0; icharge<nCharge; icharge++)
977     {
978       for (Int_t ipart=0; ipart<nPart; ipart++)
979         {
980           Int_t index = ipart + 3*icharge;
981           SpectraStep4[index]->SetMarkerSize(.5);
982           if (index == 0) SpectraStep4[index]->DrawClone("EP");
983           else SpectraStep4[index]->DrawClone("EPsame");
984         }
985     }
986   lEvolutionOfSpectra->DrawClone();
987   cEvolutionOfSpectra->cd(6);
988   gPad->SetLogy();
989   for (Int_t icharge=0; icharge<nCharge; icharge++)
990     {
991       for (Int_t ipart=0; ipart<nPart; ipart++)
992         {
993           Int_t index = ipart + 3*icharge;
994           SpectraStep5[index]->SetMarkerSize(.5);
995           if (index == 0) SpectraStep5[index]->DrawClone("EP");
996           else SpectraStep5[index]->DrawClone("EPsame");
997         }
998     }
999   lEvolutionOfSpectra->DrawClone();
1000   cEvolutionOfSpectra->cd(7);
1001   gPad->SetLogy();
1002   for (Int_t icharge=0; icharge<nCharge; icharge++)
1003     {
1004       for (Int_t ipart=0; ipart<nPart; ipart++)
1005         {
1006           Int_t index = ipart + 3*icharge;
1007           SpectraStep6[index]->SetMarkerSize(.5);
1008           if (index == 0) SpectraStep6[index]->DrawClone("EP");
1009           else SpectraStep6[index]->DrawClone("EPsame");
1010         }
1011     }
1012   lEvolutionOfSpectra->DrawClone();
1013   cEvolutionOfSpectra->cd(8);
1014   gPad->SetLogy();
1015   for (Int_t icharge=0; icharge<nCharge; icharge++)
1016     {
1017       for (Int_t ipart=0; ipart<nPart; ipart++)
1018         {
1019           Int_t index = ipart + 3*icharge;
1020           Spectra[index]->SetMarkerSize(.5);
1021           if (index == 0) Spectra[index]->DrawClone("EP");
1022           else Spectra[index]->DrawClone("EPsame");
1023         }
1024     }
1025   lEvolutionOfSpectra->DrawClone();
1026   
1027
1028   //-----------------------------------------
1029   // Compare spectra with published spectra -
1030   //-----------------------------------------
1031   // also add another canvas that just plots the ratios b/c I can't use the spectra in my presentation (hasn't been published yet)
1032   TCanvas * cRatioCombined = new TCanvas("cRatioCombined","cRatioCombined");
1033   cRatioCombined->Divide(3,1);
1034   TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",150,125,700,500);
1035   CratioComb->Divide(3,2);
1036   TString nameComb[6] = {Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
1037                          Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
1038   TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120730.root");
1039   TLegend * lComp[nPart];
1040   TLegend * lRatio[nPart];
1041   for(Int_t icharge=0; icharge<nCharge; icharge++)
1042     {
1043       for(Int_t ipart=0; ipart<nPart; ipart++)
1044         {
1045           Int_t index = ipart + 3*icharge;
1046           TH1F * htmp = (TH1F*)((TH1F*)Spectra[index])->Clone("");
1047           htmp->SetMarkerSize(1);
1048           htmp->GetXaxis()->SetRangeUser(0,5);
1049           TH1F * hcomb = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1050           CratioComb->cd(ipart+1);
1051           gPad->SetGridy();
1052           gPad->SetGridx();
1053           htmp->SetTitle("Comparison with Combined Analysis;p_{T} (GeV/c);");
1054           if(icharge == 0) htmp->DrawCopy();
1055           else htmp->DrawCopy("same");
1056           hcomb->SetMarkerStyle(Marker[index]);
1057           hcomb->DrawCopy("same");
1058           if (icharge == 0) {lComp[ipart] = new TLegend(.49,.49,.89,.89); lComp[ipart]->SetFillColor(0);}
1059           lComp[ipart]->AddEntry(htmp,Form("%s, This Analysis",Names[index].Data()),"p");
1060           lComp[ipart]->AddEntry(hcomb,Form("%s, Combined Analysis",Names[index].Data()),"lep");
1061           if (icharge == 1) lComp[ipart]->DrawClone();
1062  
1063           CratioComb->cd(ipart+nPart+1);
1064           gPad->SetGridy();
1065           gPad->SetGridx();
1066           htmp->Divide(hcomb);
1067           htmp->SetTitle("(This Analysis) / (Combined Analysis);p_{T} (GeV/c);");
1068           htmp->GetYaxis()->SetRangeUser(0.8, 1.2);
1069           if (icharge == 0)
1070             {
1071               CratioComb->cd(ipart+nPart+1);
1072               htmp->DrawClone("hist][");
1073               cRatioCombined->cd(ipart+1);
1074               htmp->DrawClone("hist][");
1075               lRatio[ipart] = new TLegend(.69,.69,.89,.89);
1076               lRatio[ipart]->SetFillColor(0);
1077             }
1078           else 
1079             {
1080               htmp->SetLineStyle(2);
1081               CratioComb->cd(ipart+nPart+1);
1082               htmp->DrawClone("histsame][");
1083               cRatioCombined->cd(ipart+1);
1084               htmp->DrawClone("histsame][");
1085             }
1086           lRatio[ipart]->AddEntry(htmp,Names[index].Data(),"l");
1087           if (icharge == 1)
1088             {
1089               CratioComb->cd(ipart+nPart+1);
1090               lRatio[ipart]->DrawClone();
1091               cRatioCombined->cd(ipart+1);
1092               lRatio[ipart]->DrawClone();    
1093             }
1094         }
1095     }
1096   //comparison with charged hadron
1097   Printf("\n\n-> ChargedHadron comparison");
1098   TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1099   TH1F* hMatchCorrectionAllCh=(TH1F*)hMatcEffPos_data->Clone("hMatchCorrectionAllCh");
1100   hMatchCorrectionAllCh->Add(hMatcEffNeg_data); //correction for Matching efficiency!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1101   hMatchCorrectionAllCh->Scale(0.5);
1102   for(Int_t ibin=1;ibin<hChHad_data->GetNbinsX();ibin++){
1103     Float_t ptch=hChHad_data->GetBinCenter(ibin);
1104     if(ptch<tcuts_data->GetPtTOFMatching())continue;
1105     hChHad_data->SetBinContent(ibin,2*(hChHad_data->GetBinContent(ibin)/(ScalingMatchingPos+ScalingMatchingNeg)));
1106   }
1107   //fraction of sec in MC
1108   TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1109   TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1110   for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
1111     Double_t en_data=hChHad_data->GetBinContent(ibin);
1112     Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
1113     Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
1114     if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
1115     //Printf("Before: %.0f After: %.0f  fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
1116   }
1117   hPrimRec_mc->Divide(hAllRec_mc);
1118   //efficiency for primaries
1119   TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1120   hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
1121   TCanvas *cAllChFactors=new TCanvas("cAllChFactors","cAllChFactors",700,500);
1122   cAllChFactors->Divide(1,2);
1123   cAllChFactors->cd(1);
1124   gPad->SetGridy();
1125   gPad->SetGridx();
1126   hPrimRec_mc->SetTitle("Prim/All, charged hadron pure MC");
1127   hPrimRec_mc->DrawClone("lhist");
1128   gPad->BuildLegend()->SetFillColor(0);
1129   cAllChFactors->cd(2);
1130   gPad->SetGridy();
1131   gPad->SetGridx();
1132   hEff_mc->SetTitle("Efficiency for Primaries, charged hadron pure MC");
1133   hEff_mc->DrawClone("lhist");
1134   gPad->BuildLegend()->SetFillColor(0);
1135   //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
1136   hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
1137   hChHad_data->Divide(hEff_mc);//Efficiency
1138   hChHad_data->Scale(1./(TMath::Abs(tcuts_data->GetEtaMin())+TMath::Abs(tcuts_data->GetEtaMax())));
1139   hChHad_data->SetTitle("All Ch from AOD");
1140   hChHad_data->SetName("AllCh");
1141   TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
1142   cAllCh->Divide(1,2);
1143   cAllCh->cd(1);
1144   gPad->SetGridy();
1145   gPad->SetGridx();
1146   hChHad_data->GetYaxis()->SetRangeUser(0,2500);
1147   hChHad_data->DrawClone();
1148   TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
1149   TH1F *hCh=fCh->Get("hpt_c0_5");
1150   hCh->SetTitle("All Ch from Jacek");
1151   hCh->SetMarkerColor(2);
1152   hCh->SetLineColor(2);
1153   //invariant yield
1154   for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
1155     hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1156     hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1157   }
1158   hCh->DrawClone("same");
1159   gPad->BuildLegend()->SetFillColor(0);
1160   TH1F *gRatio=AliPWGHistoTools::DivideHistosDifferentBins(hChHad_data,hCh);
1161   gRatio->SetMaximum(1.3);
1162   gRatio->SetMinimum(.7);
1163   cAllCh->cd(2);
1164   gPad->SetGridy();
1165   gPad->SetGridx();
1166   gRatio->Print("all");
1167   gRatio->SetBinContent(48,1.02);
1168   gRatio->SetBinContent(49,1.022);
1169   gRatio->SetBinContent(50,1.021);
1170   gRatio->SetBinContent(51,1.026);
1171   gRatio->SetBinContent(52,1.027);
1172   gRatio->SetBinError(48,0.056);
1173   gRatio->SetBinError(49,0.056);
1174   gRatio->SetBinError(50,0.056);
1175   gRatio->SetBinError(51,0.056);
1176   gRatio->SetBinError(52,0.056);
1177   gRatio->GetYaxis()->SetRangeUser(-1.5,1.5);
1178   gRatio->DrawClone("");
1179   
1180   //fitting
1181   TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
1182   gPad->SetGridy();
1183   gPad->SetGridx();
1184   hChHad_data->DrawClone();
1185   //Fitting the sum of all particles
1186   AliPWGFunc * fm = new AliPWGFunc;
1187   fm->SetVarType(AliPWGFunc::kdNdpt);
1188   Float_t fitmin = 0.3;
1189   Float_t fitmax = 3;
1190   TF1 * func = 0;
1191   Int_t normPar = 0;
1192   func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
1193   func->SetParLimits(1, 0.1, 0.99);
1194   func->SetParLimits(2, 0.01, 1);
1195   func->SetParLimits(3, 0.01, 2);
1196   TH1F * hToFit = hChHad_data;
1197   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1198   // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
1199   // //   cout << " FIT ERROR " << endl;
1200   // //   return;      
1201   // // }
1202   Double_t yieldTools, yieldETools;
1203   Double_t partialYields[3],partialYieldsErrors[3]; 
1204   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1205   func->DrawClone("same");   
1206   Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
1207   //Fit All Charged
1208   hToFit = hCh;
1209   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1210   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1211   func->SetLineColor(2);
1212   hCh->DrawClone("same");
1213   func->DrawClone("same");   
1214   gPad->BuildLegend()->SetFillColor(0);
1215   Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
1216   //sumID vs AllCh
1217   //Convert spectra to dNdeta and sum
1218   TH1F * hsum = (TH1F*) Spectra[0]->Clone();
1219   hsum->Reset("all");
1220   hsum->SetMarkerSize(1);
1221   Double_t epsilon = 0.0001;
1222   for(Int_t icharge = 0; icharge < 2; icharge++){
1223     for(Int_t ipart = 0; ipart < 3; ipart++){
1224       Int_t index=ipart+3*icharge;
1225       TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
1226       Int_t nbin = htmp->GetNbinsX();
1227       for(Int_t ibin = 1; ibin <= nbin; ibin++){
1228         Double_t pt = htmp->GetBinCenter(ibin);
1229         Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
1230         Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
1231         //Printf("jacobian: %f pt:%f   BinContent:%f  mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
1232         htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
1233         htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
1234         Int_t ibinSum = hsum->FindBin(pt);
1235         if ( htmp->GetBinContent(ibin) > 0 && 
1236              (TMath::Abs(htmp->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon || 
1237               TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1238              ) {
1239           cout << "DISCREPANCY IN BIN RANGES" << endl;
1240           cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
1241                << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
1242                << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
1243           cout << "" << endl;       
1244         }
1245       }
1246       hsum->Add(htmp);
1247     }
1248   }
1249   hsum->SetTitle("Sum ID from AOD");
1250   TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",200,150,700,500);
1251   cChargHadComp->Divide(1,2);
1252   cChargHadComp->cd(1);
1253   gPad->SetGridy();
1254   gPad->SetGridx();
1255   hsum->SetTitle("Charged Hadron Comparison;p_{T} (GeV/c);");
1256   hsum->DrawClone();
1257   hToFit = hsum;
1258   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1259   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1260   func->SetLineColor(2);
1261   Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
1262   hChHad_data->SetMarkerColor(2);
1263   hChHad_data->SetLineColor(2);
1264   hChHad_data->DrawClone("same");
1265   TLegend * lChargHadComp = new TLegend(.69,.69,.89,.89);
1266   lChargHadComp->AddEntry(hsum,"Sum ID from AOD","p");
1267   lChargHadComp->AddEntry(hChHad_data,"All Ch from AOD","p");
1268   lChargHadComp->DrawClone();
1269   cChargHadComp->cd(2);
1270   gPad->SetGridy();
1271   gPad->SetGridx();
1272   TH1F *hRatio=AliPWGHistoTools::DivideHistosDifferentBins(hsum,hChHad_data);
1273   hRatio->SetTitle("(Sum ID from AOD) / (All Ch from AOD);p_{T} (GeV/c);");
1274   hRatio->SetMaximum(1.2);
1275   hRatio->SetMinimum(.8);
1276   hRatio->DrawClone("");
1277
1278
1279   //----------------------------------------------------------
1280   // Compare each correction step with the published spectra -
1281   //----------------------------------------------------------
1282   TCanvas * cCompareToPublishedAtEachStep = new TCanvas("cCompareToPublishedAtEachStep","cCompareToPublishedAtEachStep",350,225,700,500);
1283   cCompareToPublishedAtEachStep->Divide(4,2);
1284   TLegend * lCompareToPublishedAtEachStep = new TLegend(.69,.69,.99,.99);
1285   for (Int_t icharge=0; icharge<nCharge; icharge++)
1286     {
1287       for (Int_t ipart=0; ipart<nPart; ipart++)
1288         {
1289           Int_t index = ipart + 3*icharge;
1290           lCompareToPublishedAtEachStep->AddEntry(Spectra[index],Names[index].Data(),"p");
1291         }
1292     }
1293   TH1F * hComb[nPart*nCharge];
1294   for (Int_t icharge=0; icharge<nCharge; icharge++)
1295     {
1296       for (Int_t ipart=0; ipart<nPart; ipart++)
1297         {
1298           Int_t index = ipart + 3*icharge;
1299           fComb->cd();
1300           hComb[index] = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1301         }
1302     }
1303   cCompareToPublishedAtEachStep->cd(1);
1304   for (Int_t icharge=0; icharge<nCharge; icharge++)
1305     {
1306       for (Int_t ipart=0; ipart<nPart; ipart++)
1307         {
1308           Int_t index = ipart + 3*icharge;
1309           RawSpectra[index]->Divide(hComb[index]);
1310           if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1311           if (index == 0)
1312             {
1313               RawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1314               RawSpectra[index]->DrawCopy("hist][");
1315             }
1316           else RawSpectra[index]->DrawCopy("hist][same");
1317         }
1318     }
1319   lCompareToPublishedAtEachStep->DrawClone();
1320   cCompareToPublishedAtEachStep->cd(2);
1321   for (Int_t icharge=0; icharge<nCharge; icharge++)
1322     {
1323       for (Int_t ipart=0; ipart<nPart; ipart++)
1324         {
1325           Int_t index = ipart + 3*icharge;
1326           SpectraStep1[index]->Divide(hComb[index]);
1327           if (icharge == 1) SpectraStep1[index]->SetLineStyle(2);
1328           if (index == 0)
1329             {
1330               SpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1331               SpectraStep1[index]->DrawCopy("hist][");
1332             }
1333           else SpectraStep1[index]->DrawCopy("hist][same");
1334         }
1335     }
1336   lCompareToPublishedAtEachStep->DrawClone();
1337   cCompareToPublishedAtEachStep->cd(3);
1338   for (Int_t icharge=0; icharge<nCharge; icharge++)
1339     {
1340       for (Int_t ipart=0; ipart<nPart; ipart++)
1341         {
1342           Int_t index = ipart + 3*icharge;
1343           SpectraStep2[index]->Divide(hComb[index]);
1344           if (icharge == 1) SpectraStep2[index]->SetLineStyle(2);
1345           if (index == 0)
1346             {
1347               SpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1348               SpectraStep2[index]->DrawCopy("hist][");
1349             }
1350           else SpectraStep2[index]->DrawCopy("hist][same");
1351         }
1352     }
1353   lCompareToPublishedAtEachStep->DrawClone();
1354   cCompareToPublishedAtEachStep->cd(4);
1355   for (Int_t icharge=0; icharge<nCharge; icharge++)
1356     {
1357       for (Int_t ipart=0; ipart<nPart; ipart++)
1358         {
1359           Int_t index = ipart + 3*icharge;
1360           SpectraStep3[index]->Divide(hComb[index]);
1361           if (icharge == 1) SpectraStep3[index]->SetLineStyle(2);
1362           if (index == 0)
1363             {
1364               SpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1365               SpectraStep3[index]->DrawCopy("hist][");
1366             }
1367           else SpectraStep3[index]->DrawCopy("hist][same");
1368         }
1369     }
1370   lCompareToPublishedAtEachStep->DrawClone();
1371   cCompareToPublishedAtEachStep->cd(5);
1372   for (Int_t icharge=0; icharge<nCharge; icharge++)
1373     {
1374       for (Int_t ipart=0; ipart<nPart; ipart++)
1375         {
1376           Int_t index = ipart + 3*icharge;
1377           SpectraStep4[index]->Divide(hComb[index]);
1378           if (icharge == 1) SpectraStep4[index]->SetLineStyle(2);
1379           if (index == 0)
1380             {
1381               SpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1382               SpectraStep4[index]->DrawCopy("hist][");
1383             }
1384           else SpectraStep4[index]->DrawCopy("hist][same");
1385         }
1386     }
1387   lCompareToPublishedAtEachStep->DrawClone();
1388   cCompareToPublishedAtEachStep->cd(6);
1389   for (Int_t icharge=0; icharge<nCharge; icharge++)
1390     {
1391       for (Int_t ipart=0; ipart<nPart; ipart++)
1392         {
1393           Int_t index = ipart + 3*icharge;
1394           SpectraStep5[index]->Divide(hComb[index]);
1395           if (icharge == 1) SpectraStep5[index]->SetLineStyle(2);
1396           if (index == 0)
1397             {
1398               SpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1399               SpectraStep5[index]->DrawCopy("hist][");
1400             }
1401           else SpectraStep5[index]->DrawCopy("hist][same");
1402         }
1403     }
1404   lCompareToPublishedAtEachStep->DrawClone();
1405   cCompareToPublishedAtEachStep->cd(7);
1406   for (Int_t icharge=0; icharge<nCharge; icharge++)
1407     {
1408       for (Int_t ipart=0; ipart<nPart; ipart++)
1409         {
1410           Int_t index = ipart + 3*icharge;
1411           SpectraStep6[index]->Divide(hComb[index]);
1412           if (icharge == 1) SpectraStep6[index]->SetLineStyle(2);
1413           if (index == 0)
1414             {
1415               SpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1416               SpectraStep6[index]->DrawCopy("hist][");
1417             }
1418           else SpectraStep6[index]->DrawCopy("hist][same");
1419         }
1420     }
1421   lCompareToPublishedAtEachStep->DrawClone();
1422   cCompareToPublishedAtEachStep->cd(8);
1423   for (Int_t icharge=0; icharge<nCharge; icharge++)
1424     {
1425       for (Int_t ipart=0; ipart<nPart; ipart++)
1426         {
1427           Int_t index = ipart + 3*icharge;
1428           Spectra[index]->Divide(hComb[index]);
1429           if (icharge == 1) Spectra[index]->SetLineStyle(2);
1430           if (index == 0)
1431             {
1432               Spectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1433               Spectra[index]->DrawCopy("hist][");
1434             }
1435           else Spectra[index]->DrawCopy("hist][same");
1436         }
1437     }
1438   lCompareToPublishedAtEachStep->DrawClone();
1439
1440
1441
1442   //---------------------------------------------------------
1443   // For MC, compare each correction step with the MC truth -
1444   //---------------------------------------------------------
1445   TCanvas * cCompareMCtoMCTruthAtEachStep = new TCanvas("cCompareMCtoMCTruthAtEachStep","cCompareMCtoMCTruthAtEachStep",350,225,700,500);
1446   cCompareMCtoMCTruthAtEachStep->Divide(4,2);
1447   TLegend * lCompareMCtoMCTruthAtEachStep = new TLegend(.69,.69,.99,.99);
1448   for (Int_t icharge=0; icharge<nCharge; icharge++)
1449     {
1450       for (Int_t ipart=0; ipart<nPart; ipart++)
1451         {
1452           Int_t index = ipart + 3*icharge;
1453           lCompareMCtoMCTruthAtEachStep->AddEntry(MCSpectra[index],Names[index].Data(),"p");
1454         }
1455     }
1456   cCompareMCtoMCTruthAtEachStep->cd(1);
1457   for (Int_t icharge=0; icharge<nCharge; icharge++)
1458     {
1459       for (Int_t ipart=0; ipart<nPart; ipart++)
1460         {
1461           Int_t index = ipart + 3*icharge;
1462           MCRawSpectra[index]->Divide(hMCTruth[index]);
1463           if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1464           if (index == 0)
1465             {
1466               MCRawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1467               MCRawSpectra[index]->DrawCopy("hist][");
1468             }
1469           else MCRawSpectra[index]->DrawCopy("hist][same");
1470         }
1471     }
1472   lCompareMCtoMCTruthAtEachStep->DrawClone();
1473   cCompareMCtoMCTruthAtEachStep->cd(2);
1474   for (Int_t icharge=0; icharge<nCharge; icharge++)
1475     {
1476       for (Int_t ipart=0; ipart<nPart; ipart++)
1477         {
1478           Int_t index = ipart + 3*icharge;
1479           MCSpectraStep1[index]->Divide(hMCTruth[index]);
1480           if (icharge == 1) MCSpectraStep1[index]->SetLineStyle(2);
1481           if (index == 0)
1482             {
1483               MCSpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1484               MCSpectraStep1[index]->DrawCopy("hist][");
1485             }
1486           else MCSpectraStep1[index]->DrawCopy("hist][same");
1487         }
1488     }
1489   lCompareMCtoMCTruthAtEachStep->DrawClone();
1490   cCompareMCtoMCTruthAtEachStep->cd(3);
1491   for (Int_t icharge=0; icharge<nCharge; icharge++)
1492     {
1493       for (Int_t ipart=0; ipart<nPart; ipart++)
1494         {
1495           Int_t index = ipart + 3*icharge;
1496           MCSpectraStep2[index]->Divide(hMCTruth[index]);
1497           if (icharge == 1) MCSpectraStep2[index]->SetLineStyle(2);
1498           if (index == 0)
1499             {
1500               MCSpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1501               MCSpectraStep2[index]->DrawCopy("hist][");
1502             }
1503           else MCSpectraStep2[index]->DrawCopy("hist][same");
1504         }
1505     }
1506   lCompareMCtoMCTruthAtEachStep->DrawClone();
1507   cCompareMCtoMCTruthAtEachStep->cd(4);
1508   for (Int_t icharge=0; icharge<nCharge; icharge++)
1509     {
1510       for (Int_t ipart=0; ipart<nPart; ipart++)
1511         {
1512           Int_t index = ipart + 3*icharge;
1513           MCSpectraStep3[index]->Divide(hMCTruth[index]);
1514           if (icharge == 1) MCSpectraStep3[index]->SetLineStyle(2);
1515           if (index == 0)
1516             {
1517               MCSpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1518               MCSpectraStep3[index]->DrawCopy("hist][");
1519             }
1520           else MCSpectraStep3[index]->DrawCopy("hist][same");
1521         }
1522     }
1523   lCompareMCtoMCTruthAtEachStep->DrawClone();
1524   /*  cCompareMCtoMCTruthAtEachStep->cd(5);
1525   for (Int_t icharge=0; icharge<nCharge; icharge++)
1526     {
1527       for (Int_t ipart=0; ipart<nPart; ipart++)
1528         {
1529           Int_t index = ipart + 3*icharge;
1530           MCSpectraStep4[index]->Divide(hMCTruth[index]);
1531           if (icharge == 1) MCSpectraStep4[index]->SetLineStyle(2);
1532           if (index == 0)
1533             {
1534               MCSpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1535               MCSpectraStep4[index]->DrawCopy("hist][");
1536             }
1537           else MCSpectraStep4[index]->DrawCopy("hist][same");
1538         }
1539     }
1540     lCompareMCtoMCTruthAtEachStep->DrawClone();
1541   cCompareMCtoMCTruthAtEachStep->cd(6);
1542   for (Int_t icharge=0; icharge<nCharge; icharge++)
1543     {
1544       for (Int_t ipart=0; ipart<nPart; ipart++)
1545         {
1546           Int_t index = ipart + 3*icharge;
1547           MCSpectraStep5[index]->Divide(hMCTruth[index]);
1548           if (icharge == 1) MCSpectraStep5[index]->SetLineStyle(2);
1549           if (index == 0)
1550             {
1551               MCSpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1552               MCSpectraStep5[index]->DrawCopy("hist][");
1553             }
1554           else MCSpectraStep5[index]->DrawCopy("hist][same");
1555         }
1556     }
1557     lCompareMCtoMCTruthAtEachStep->DrawClone();*/
1558   cCompareMCtoMCTruthAtEachStep->cd(7);
1559   for (Int_t icharge=0; icharge<nCharge; icharge++)
1560     {
1561       for (Int_t ipart=0; ipart<nPart; ipart++)
1562         {
1563           Int_t index = ipart + 3*icharge;
1564           MCSpectraStep6[index]->Divide(hMCTruth[index]);
1565           if (icharge == 1) MCSpectraStep6[index]->SetLineStyle(2);
1566           if (index == 0)
1567             {
1568               MCSpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1569               MCSpectraStep6[index]->DrawCopy("hist][");
1570             }
1571           else MCSpectraStep6[index]->DrawCopy("hist][same");
1572         }
1573     }
1574   lCompareMCtoMCTruthAtEachStep->DrawClone();
1575   cCompareMCtoMCTruthAtEachStep->cd(8);
1576   for (Int_t icharge=0; icharge<nCharge; icharge++)
1577     {
1578       for (Int_t ipart=0; ipart<nPart; ipart++)
1579         {
1580           Int_t index = ipart + 3*icharge;
1581           MCSpectra[index]->Divide(hMCTruth[index]);
1582           if (icharge == 1) MCSpectra[index]->SetLineStyle(2);
1583           if (index == 0)
1584             {
1585               MCSpectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1586               MCSpectra[index]->DrawCopy("hist][");
1587             }
1588           else MCSpectra[index]->DrawCopy("hist][same");
1589         }
1590     }
1591   lCompareMCtoMCTruthAtEachStep->DrawClone();
1592
1593
1594
1595   // close unwanted canvases
1596   //  cGFCorrection->Close();
1597   //cMatchingEff->Close();
1598   //cContCorrFact->Close();
1599   //cCorrFactMC->Close();
1600   //cSecCorrFactMC->Close();
1601   cAllCh->Close();
1602   cFitChargHad->Close();
1603   cAllChFactors->Close();
1604
1605   // save wanted canvases
1606   fout->cd();
1607   cSpectra->Write();
1608   CratioComb->Write();
1609   cMatchingEff->Write();
1610   cEvolutionOfSpectra->Write();  
1611   cCompareToPublishedAtEachStep->Write();
1612 }
1613
1614
1615
1616 ///////////////////////
1617 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
1618   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1619   return TMath::ASinH(pt / mt * TMath::SinH(eta));
1620 }
1621
1622 ///////////////////////
1623 void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data, TFile * fout){
1624   //Geant/Fluka Correction
1625   Printf("\nGF correction for Kaons");
1626   //Getting GF For Kaons in TPC
1627   TGraph *gGFCorrectionKaonPlus=new TGraph();
1628   gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
1629   gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
1630   TGraph *gGFCorrectionKaonMinus=new TGraph();
1631   gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
1632   gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
1633   TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
1634   TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
1635   TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
1636   TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
1637   //getting GF func for Kaons with TOF
1638   TF1 *fGFKPosTracking;
1639   fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
1640   TF1 *fGFKNegTracking;
1641   fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
1642   TF1 *fGFKPosMatching;
1643   fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
1644   TF1 *fGFKNegMatching;
1645   fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
1646   for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
1647     if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1648       Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
1649       Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
1650       Printf("TPC Geant/Fluka: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
1651       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
1652       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
1653       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
1654       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
1655       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
1656       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
1657     }else{
1658       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1659       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1660       Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1661       Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1662       Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1663       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1664       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1665       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1666       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1667       Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1668       Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1669       Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1670       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1671       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1672       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1673       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1674     }
1675   }
1676   
1677   //Geant Fluka for P in TPC
1678   Printf("\nGF correction for Protons");
1679   const Int_t kNCharge=2;
1680   Int_t kPos=0;
1681   Int_t kNeg=1;
1682   TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
1683   TH2D * hCorrFluka[kNCharge];
1684   TH2D * hCorrFluka[2];
1685   hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1686   hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1687   TGraph *gGFCorrectionProtonPlus=new TGraph();
1688   gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
1689   gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
1690   TGraph *gGFCorrectionProtonMinus=new TGraph();
1691   gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
1692   gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
1693   //getting GF func for Protons with TPCTOF
1694   TF1 *fGFpPosTracking;
1695   fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1696   TF1 *fGFpNegTracking;
1697   fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1698   TF1 *fGFpPosMatching;
1699   fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1700   TF1 *fGFpNegMatching;
1701   fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1702   
1703   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1704     Int_t nbins = Spectra[2]->GetNbinsX();
1705     Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1706     for(Int_t ibin = 0; ibin < nbins; ibin++){
1707       if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1708         Float_t pt = Spectra[2]->GetBinCenter(ibin);
1709         Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1710         Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1711         if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
1712         if (pt > maxPtCorrection) pt = maxPtCorrection;
1713         Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1714         //cout<<correction<<"     charge "<<icharge<<endl;
1715         if(icharge==0){
1716           if (correction != 0) {// If the bin is empty this is a  0
1717             Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
1718             Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError  (ibin)*correction);
1719             gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
1720           }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1721             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
1722             cout << " Bin content: " << Spectra[2]->GetBinContent(ibin)  << endl;
1723           }
1724         }
1725         if(icharge==1){
1726           if (correction != 0) {// If the bin is empty this is a  0
1727             Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
1728             Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError  (ibin)*correction);
1729             gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
1730           }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1731             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
1732             cout << " Bin content: " << Spectra[5]->GetBinContent(ibin)  << endl;
1733           }
1734         }
1735       }else{
1736         gGFCorrectionProtonPlus->SetPoint(ibin,Spectra[2]->GetBinCenter(ibin),0);
1737         gGFCorrectionProtonMinus->SetPoint(ibin,Spectra[5]->GetBinCenter(ibin),0);
1738         Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1739         Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1740         Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1741         if(icharge==0){
1742           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1743           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1744         }else{
1745           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1746           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1747         }
1748         Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1749         Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1750         Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1751         if(icharge==0){
1752           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1753           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1754         }else{
1755           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1756           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1757         }
1758       }
1759     }//end loop on bins 
1760   }
1761   gGFCorrectionKaonPlus->SetLineColor(kRed);
1762   gGFCorrectionKaonMinus->SetLineColor(kRed+2);
1763   gGFCorrectionProtonPlus->SetLineColor(kGreen);
1764   gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
1765   fGFKPosTracking->SetLineColor(kRed);
1766   fGFKNegTracking->SetLineColor(kRed+2);
1767   fGFKPosMatching->SetLineColor(kRed);
1768   fGFKNegMatching->SetLineColor(kRed+2);
1769   fGFpPosTracking->SetLineColor(kGreen);
1770   fGFpNegTracking->SetLineColor(kGreen+2);
1771   fGFpPosMatching->SetLineColor(kGreen);
1772   fGFpNegMatching->SetLineColor(kGreen+2);
1773   fGFKPosTracking->SetLineStyle(2);
1774   fGFKNegTracking->SetLineStyle(2);
1775   fGFKPosMatching->SetLineStyle(3);
1776   fGFKNegMatching->SetLineStyle(3);
1777   fGFpPosTracking->SetLineStyle(2);
1778   fGFpNegTracking->SetLineStyle(2);
1779   fGFpPosMatching->SetLineStyle(3);
1780   fGFpNegMatching->SetLineStyle(3);
1781   fGFKPosTracking->SetRange(.6,5);
1782   fGFKNegTracking->SetRange(.6,5);
1783   fGFKPosMatching->SetRange(.6,5);
1784   fGFKNegMatching->SetRange(.6,5);
1785   fGFpPosTracking->SetRange(.6,5);
1786   fGFpNegTracking->SetRange(.6,5);
1787   fGFpPosMatching->SetRange(.6,5);
1788   fGFpNegMatching->SetRange(.6,5);
1789  
1790   TCanvas * cGFCorrection = new TCanvas ("cGFCorrection","cGFCorrection",700,500);
1791   gPad->SetGridx();
1792   gPad->SetGridy();
1793   gGFCorrectionKaonPlus->DrawClone("al");
1794   gGFCorrectionKaonMinus->DrawClone("lsame");
1795   gGFCorrectionProtonPlus->DrawClone("lsame");
1796   gGFCorrectionProtonMinus->DrawClone("lsame");
1797   fGFKPosTracking->DrawClone("lsame");
1798   fGFKNegTracking->DrawClone("lsame");
1799   fGFKPosMatching->DrawClone("lsame");
1800   fGFKNegMatching->DrawClone("lsame");
1801   fGFpPosTracking->DrawClone("lsame");
1802   fGFpNegTracking->DrawClone("lsame");
1803   fGFpPosMatching->DrawClone("lsame");
1804   fGFpNegMatching->DrawClone("lsame");
1805   gPad->BuildLegend()->SetFillColor(0);
1806   fout->cd();
1807   cGFCorrection->Write();
1808 }
1809
1810 ///////////
1811 TF1 * TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1812 {
1813
1814   if (ipart == 3 && icharge == kNegative) {
1815     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1816     return f;
1817   }
1818   else if (ipart == 4 && icharge == kNegative) {
1819     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1820   }
1821   else
1822     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1823
1824   return f;
1825 }
1826
1827 Double_t TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1828 {
1829   return 1.;
1830 }
1831
1832 Double_t TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1833 {
1834   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1835 }
1836
1837 Double_t TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1838 {
1839   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1840 }
1841
1842 ///////////////////////////////////////////
1843 TF1 * TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1844 {
1845   if (ipart == 3 && icharge == kNegative) {
1846     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1847     return f;
1848   }
1849   else if (ipart == 4 && icharge == kNegative) {
1850     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1851   }
1852   else
1853     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1854
1855   return f;
1856 }
1857
1858 Double_t MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1859 {
1860   return 1.;
1861 }
1862
1863 Double_t MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1864 {
1865   Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1866   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1867 }
1868
1869 Double_t MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1870 {
1871   Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1872   return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1873 }