Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / fit_yields_old.C
1
2 #include <TFile.h>
3 #include <TH1.h>
4 #include <TProfile.h>
5 #include <TF1.h>
6 #include <TGraphErrors.h>
7 #include <TCanvas.h>
8 #include <TAxis.h>
9 #include <TLegend.h>
10 #include <TStyle.h>
11 #include <TGraphAsymmErrors.h>
12 #include <TList.h>
13 #include <TMath.h>
14 #include <TSystem.h>
15 #include <TLatex.h>
16 #include <TF1.h>
17 #include <TF2.h>
18
19 #include "AliHighPtDeDxData.h"
20
21 #include "my_tools.C"
22 #include "my_functions.C"
23
24 #include <iostream>
25 #include <fstream>
26 #include <string>
27
28 using namespace std;
29
30 /*
31   Ideas to improve:
32   =================
33
34   Use the real mean p -> Might give some effect
35   Effect: Push the <dE/dx> down (let'a see if it is still needed in the fits)
36
37   Use the real <nCl> 
38   Effect: Increase sigma for p<15. For p>15 seems to saturate.
39   
40   To use:
41   =======
42   root is enough
43   .L libMyDeDxAnalysis.so 
44   .L my_tools.C+
45   .L my_functions.C+
46   .L fit_yields.C+
47
48
49   mkdir results
50   mkdir results/fits
51   mkdir results/eta04
52   mkdir results/eta48
53   mkdir results/etaneg
54   mkdir results/etapos
55   mkdir results/neg
56   mkdir results/pos
57   mkdir results/comparison
58   mkdir results/kaon
59   mkdir results/lambda
60
61
62   *********************
63        ALL ETAS
64   *********************
65   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all.root", 0, "results/fits") 
66
67   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, -1, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_neg.root", 0, "results/neg") 
68
69   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, +1, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_pos.root", 0, "results/pos") 
70
71   Compare("fit_yields_results/lhc10h_aod_all.root", "fit_yields_results/lhc10h_aod_all_neg.root", "fit_yields_results/lhc10h_aod_all_pos.root", "q<0", "q>0", "results/comparison/neg_vs_pos")
72
73   CompareYields("data/lhc10h_aod_all.root", 0, 3.0, 20.0, -1, +1, 0, 0, kFALSE, 0, 1, kFALSE)
74
75
76
77   *********************
78      ETA LOW & HIGH
79   *********************
80
81   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_low.root", "etaabs04", "results/eta04") 
82
83   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_high.root", "etaabs48", "results/eta48") 
84
85   Compare("fit_yields_results/lhc10h_aod_all.root", "fit_yields_results/lhc10h_aod_all_eta_low.root", "fit_yields_results/lhc10h_aod_all_eta_high.root", "|#eta|<0.4", "0.4<|#eta|<0.8", "results/comparison/low_vs_high")
86
87   CompareYields("data/lhc10h_aod_all.root", 0, 3.0, 20.0, 0, 0, "etaabs04", "etaabs48")
88
89
90   *********************
91      ETA NEG & POS
92   *********************
93   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_neg.root", "eta-80", "results/etaneg") 
94
95   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_pos.root", "eta08", "results/etapos") 
96
97   Compare("fit_yields_results/lhc10h_aod_all.root", "fit_yields_results/lhc10h_aod_all_eta_neg.root", "fit_yields_results/lhc10h_aod_all_eta_pos.root", "#eta < 0", "#eta > 0", "results/comparison/etaneg_vs_etapos")
98
99   CompareYields("data/lhc10h_aod_all.root", 0, 3.0, 20.0, 0, 0, "eta-80", "eta08")
100
101
102   *********************
103      V0 ETA LOW & HIGH
104   *********************
105   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "kaon", 3.0, 20.0, 0, kTRUE, 1, 0, kTRUE, "etaabs04") 
106   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "kaon", 3.0, 20.0, 0, kTRUE, 1, 0, kTRUE, "etaabs48") 
107   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "lambda", 3.0, 20.0, 0, kTRUE, 1, 0, kTRUE, "etaabs04") 
108   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "lambda", 3.0, 20.0, 0, kTRUE, 1, 0, kTRUE, "etaabs48") 
109
110   *********************
111
112   *********************
113      ETA FINE
114   *********************
115   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_86.root", "eta-8-6", "results/eta86") 
116   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_64.root", "eta-6-4", "results/eta64") 
117   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_42.root", "eta-4-2", "results/eta42") 
118   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_20.root", "eta-20", "results/eta20") 
119   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_02.root", "eta02", "results/eta02") 
120   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_24.root", "eta24", "results/eta24") 
121   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_46.root", "eta46", "results/eta46") 
122   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE, 0, 1, kTRUE, "lhc10h_aod_all_eta_68.root", "eta68", "results/eta68") 
123
124   Compare(3.1)
125   Compare(5.1)
126   Compare(7.1)
127   Compare(9.1)
128
129  FitYields("data/lhc10h_mc_aod_all.root", 3.0, 20.0, 0, kTRUE)   
130
131   FitYields("data/lhc10h_aod_test.root", 3.0, 20.0, 0, kTRUE)
132
133   FitYields("data/lhc10h_aod_all.root", 3.0, 20.0, 0, kTRUE)
134
135   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "kaon", 3.0, 20.0, 0, kTRUE)
136   FitYieldsV0("data/lhc10h_aod_all.root", "data/lhc10h_v0_all_loose.root", "lambda", 3.0, 20.0, 0, kTRUE)
137
138
139
140   MakeRatios("../ptspectraPP/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPP/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPP/results/comparison/")
141
142   MakeRatios("../ptspectraPbPb_60_80/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPbPb_60_80/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPbPb_60_80/results/comparison/")
143
144   MakeRatios("../ptspectraPbPb_40_60/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPbPb_40_60/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPbPb_40_60/results/comparison/")
145
146   MakeRatios("../ptspectraPbPb_20_40/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPbPb_20_40/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPbPb_20_40/results/comparison/")
147
148   MakeRatios("../ptspectraPbPb_10_20/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPbPb_10_20/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPbPb_10_20/results/comparison/")
149
150   MakeRatios("../ptspectraPbPb_0_5/fit_yields_results/lhc10h_aod_all_neg.root", "../ptspectraPbPb_0_5/fit_yields_results/lhc10h_aod_all_pos.root", "../ptspectraPbPb_0_5/results/comparison/")
151
152   MakeRatios("fit_yields_results/lhc10h_aod_all_neg.root", "fit_yields_results/lhc10h_aod_all_pos.root", "results/comparison/")
153
154
155  */
156
157 void MakeRatios(const Char_t* file1Name, const Char_t* file2Name, 
158                 Bool_t drawFractionRatios,
159                 Bool_t drawYieldRatios);
160
161
162 //___________________________________________________________________________________________
163 void MakeRatios(const Char_t* fileNeg, const Char_t* filePos, 
164                 const Char_t* outdirname)
165 {
166   /*
167     For yields we assume that file 1 is negative charge and file
168     2 is positive charge.
169    */
170   
171   TFile* file1 = FindFileFresh(fileNeg);
172   if(!file1)
173     return;
174   
175   TH1D* hPionRatio1   = (TH1D*)file1->Get("hPionYield");
176   SetHistError(hPionRatio1, 0.0);
177   // TH1D* hKaonRatio1   = (TH1D*)file1->Get("hKaonRatio");
178   // TH1D* hProtonRatio1 = (TH1D*)file1->Get("hProtonRatio");
179
180
181   TFile* file2 = FindFileFresh(filePos);
182   if(!file2)
183     return;
184   
185   TH1D* hPionRatio2   = (TH1D*)file2->Get("hPionYield");
186   // TH1D* hKaonRatio2   = (TH1D*)file2->Get("hKaonRatio");
187   // TH1D* hProtonRatio2 = (TH1D*)file2->Get("hProtonRatio");
188
189   
190   TH1D* hPionAssymetry = (TH1D*)hPionRatio1->Clone("hPionAssymetry");
191   hPionAssymetry->Divide(hPionRatio2);
192   TF1 f1("f1", "-1.0", 0.0, 50.0);
193   hPionAssymetry->Add(&f1); // correct for muons and electrons
194   CutHistogram(hPionAssymetry, 3.0, 20.0);
195   hPionAssymetry->GetYaxis()->SetRangeUser(-0.25, 0.25);
196   hPionAssymetry->SetTitle("Syst. error evaluation: Assymetry; p_{T} [GeV/c]; pion assymmetry: (neg-pos)/pos");
197   
198   
199   hPionRatio1->Divide(hPionRatio2);
200   hPionRatio1->GetYaxis()->SetRangeUser(0.75, 1.25);
201   hPionRatio1->SetTitle("Syst. error evaluation: Ratio; p_{T} [GeV/c]; pion fraction ratio: neg/pos");
202
203
204   TCanvas* cPionFractionRatio = new TCanvas("cPionFractionRatio", "pion fraction ratio", 400, 300);
205   cPionFractionRatio->Clear();
206   cPionFractionRatio->SetGridy();
207   cPionFractionRatio->cd();
208   hPionRatio1->Draw();
209   if(strstr(fileNeg, "PP"))
210     gROOT->ProcessLine(".x ../ptspectraPP/drawText.C");
211   else if(strstr(fileNeg, "0_5"))
212     gROOT->ProcessLine(".x ../ptspectraPbPb_0_5/drawText.C");
213   else if(strstr(fileNeg, "5_10"))
214     gROOT->ProcessLine(".x ../ptspectraPbPb_5_10/drawText.C");
215   else if(strstr(fileNeg, "10_20"))
216     gROOT->ProcessLine(".x ../ptspectraPbPb_10_20/drawText.C");
217   else if(strstr(fileNeg, "20_40"))
218     gROOT->ProcessLine(".x ../ptspectraPbPb_20_40/drawText.C");
219   else if(strstr(fileNeg, "40_60"))
220     gROOT->ProcessLine(".x ../ptspectraPbPb_40_60/drawText.C");
221   else if(strstr(fileNeg, "60_80"))
222     gROOT->ProcessLine(".x ../ptspectraPbPb_60_80/drawText.C");
223
224   cPionFractionRatio->SaveAs(Form("%s/pion_frac_ratio_neg_over_pos.gif", outdirname));
225   cPionFractionRatio->SaveAs(Form("%s/pion_frac_ratio_neg_over_pos.pdf", outdirname));
226
227
228   // TCanvas* cPionAssymetry = new TCanvas("cPionAssymetry", "pion fraction assymetry", 400, 300);
229   // cPionAssymetry->Clear();
230   // cPionAssymetry->SetGridy();
231   // cPionAssymetry->cd();
232   // hPionAssymetry->Draw();
233   // if(strstr(fileNeg, "PP"))
234   //   gROOT->ProcessLine(".x ../ptspectraPP/drawText.C");
235   // else if(strstr(fileNeg, "0_5"))
236   //   gROOT->ProcessLine(".x ../ptspectraPbPb_0_5/drawText.C");
237   // else if(strstr(fileNeg, "5_10"))
238   //   gROOT->ProcessLine(".x ../ptspectraPbPb_5_10/drawText.C");
239   // else if(strstr(fileNeg, "10_20"))
240   //   gROOT->ProcessLine(".x ../ptspectraPbPb_10_20/drawText.C");
241   // else if(strstr(fileNeg, "20_40"))
242   //   gROOT->ProcessLine(".x ../ptspectraPbPb_20_40/drawText.C");
243   // else if(strstr(fileNeg, "40_60"))
244   //   gROOT->ProcessLine(".x ../ptspectraPbPb_40_60/drawText.C");
245   // else if(strstr(fileNeg, "60_80"))
246   //   gROOT->ProcessLine(".x ../ptspectraPbPb_60_80/drawText.C");
247   // cPionAssymetry->SaveAs(Form("%s/pion_assymetry_between_neg_and_pos.gif", outdirname));
248   // cPionAssymetry->SaveAs(Form("%s/pion_assymetry_between_neg_and_pos.pdf", outdirname));  
249
250
251 }
252
253 //____________________________________________________________________________
254 void FitYields(const Char_t* dataFileName,
255                Double_t ptStart, Double_t ptStop,
256                Int_t charge,
257                Bool_t performFit = kFALSE,
258                Int_t run    = 0,
259                Int_t filter = 1,
260                Bool_t usePhiCut = kTRUE,
261                const Char_t* outFileName=0,
262                const Char_t* endName=0,
263                const Char_t* dirName="debugfits")
264 {
265   gStyle->SetOptStat(0);
266
267   
268   TFile* dataFile = FindFileFresh(dataFileName);
269   if(!dataFile)
270     return;
271   AliHighPtDeDxData* data = (AliHighPtDeDxData*)GetObject(dataFile, filter, usePhiCut, run, "filter", endName);
272   data->Print();
273
274   gSystem->Exec(Form("mv %s/* old/", dirName));
275   if(data->IsMc())
276     gSystem->Exec("mv debugfitsmc/* olddebugfitsmc/");
277
278
279   TH2D* hDeltaPiVsPt = data->GetHistDeltaPiVsPt(0, charge);
280   hDeDxVsP = hDeltaPiVsPt; // for the 2d fit to pick up the right bin
281
282   TH2D* hDeltaPiVsPtPiGen = data->GetHistDeltaPiVsPt(1, charge);
283   TH2D* hDeltaPiVsPtKGen  = data->GetHistDeltaPiVsPt(2, charge);
284   TH2D* hDeltaPiVsPtPGen  = data->GetHistDeltaPiVsPt(3, charge);
285
286   TH2D* hDeltaPiVsPtPi = 0;
287   TH2D* hDeltaPiVsPtK  = 0;
288   TH2D* hDeltaPiVsPtP  = 0;
289
290   if(data->IsMc()) {
291
292     hDeltaPiVsPtPi = data->GetHistDeltaPiVsPtMc(1, charge);
293     hDeltaPiVsPtK  = data->GetHistDeltaPiVsPtMc(2, charge);
294     hDeltaPiVsPtP  = data->GetHistDeltaPiVsPtMc(3, charge);
295   }
296
297
298
299   TProfile* hPiGenProfile = hDeltaPiVsPtPiGen->ProfileX();
300   hPiGenProfile->SetMarkerStyle(29);
301   TProfile* hKGenProfile = hDeltaPiVsPtKGen->ProfileX();
302   hKGenProfile->SetMarkerStyle(29);
303   TProfile* hPGenProfile = hDeltaPiVsPtPGen->ProfileX();
304   hPGenProfile->SetMarkerStyle(29);
305
306   TCanvas* cDeltaPiVsPt = new TCanvas("cDeltaPiVsPt", "dE/dx vs p", 400, 300);
307   cDeltaPiVsPt->Clear();
308   cDeltaPiVsPt->cd();
309   cDeltaPiVsPt->SetLogz();
310   hDeltaPiVsPt->Draw("COLZ");
311   hPiGenProfile->Draw("SAME P");
312   hKGenProfile->Draw("SAME P");
313   hPGenProfile->Draw("SAME P");
314   gROOT->ProcessLine(".x drawText.C");
315   cDeltaPiVsPt->SaveAs(Form("%s/deltapi_vs_pt.gif", dirName));
316   cDeltaPiVsPt->SaveAs(Form("%s/deltapi_vs_pt.pdf", dirName));
317
318   TCanvas* cDeltaPiVsPtLogX = new TCanvas("cDeltaPiVsPtLogX", "dE/dx vs p", 400, 300);
319   cDeltaPiVsPtLogX->Clear();
320   cDeltaPiVsPtLogX->cd();
321   cDeltaPiVsPtLogX->SetLogz();
322   cDeltaPiVsPtLogX->SetLogx();
323   hDeltaPiVsPt->Draw("COLZ");
324   hPiGenProfile->Draw("SAME P");
325   hKGenProfile->Draw("SAME P");
326   hPGenProfile->Draw("SAME P");
327   gROOT->ProcessLine(".x drawText.C");
328   cDeltaPiVsPtLogX->SaveAs(Form("%s/deltapi_vs_pt_logx.gif", dirName));
329   cDeltaPiVsPtLogX->SaveAs(Form("%s/deltapi_vs_pt_logx.pdf", dirName));
330
331   // Root is a bit stupid with finidng bins so we have to add and subtract a
332   // little to be sure we get the right bin as we typically put edges as
333   // limits
334   const Int_t binStart = hDeltaPiVsPt->GetXaxis()->FindBin(ptStart+0.01);
335   ptStart = hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(binStart);
336   const Int_t binStop  = hDeltaPiVsPt->GetXaxis()->FindBin(ptStop-0.01);
337   ptStop = hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(binStop);
338   //  const Int_t nBins    = binStop - binStart + 1;
339
340   cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
341        << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
342   
343
344   //cross check
345   TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
346   cFits->Clear();
347   cFits->Divide(7, 4);
348
349   // TF1* pionGen = new TF1("pionGen", "[0]/sqrt(6.2832*[2]*[2])*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", -30, 20);
350   // pionGen->SetLineWidth(2);
351   // pionGen->SetLineColor(kRed);
352
353   TF1* pion = new TF1("pion", "gausn", -30, 20);
354   pion->SetLineWidth(2);
355   pion->SetLineColor(kRed);
356   TF1* kaon = new TF1("kaon", "gausn", -30, 20);
357   kaon->SetLineWidth(2);
358   kaon->SetLineColor(kGreen);
359   TF1* proton = new TF1("proton", "gausn", -30, 20);
360   proton->SetLineWidth(2);
361   proton->SetLineColor(kBlue);
362   TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", -30, 20);
363   total->SetLineColor(kBlack);
364   total->SetLineWidth(2);
365   total->SetLineStyle(2);
366
367   TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);    
368   legend->SetBorderSize(0);
369   legend->SetFillColor(0);
370   legend->AddEntry(total, "3-Gauss fit", "L");
371   legend->AddEntry(pion, "#pi", "L");
372   legend->AddEntry(kaon, "K", "L");
373   legend->AddEntry(proton, "p", "L");
374
375   TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
376   //  cSingleFit->SetLogy();
377
378   TH1D* hPionRatio =(TH1D*)hDeltaPiVsPt->ProjectionX("hPionRatio", 1, 1);
379   hPionRatio->SetTitle("particle fractions; p_{T} [GeV/c]; particle fraction");
380   hPionRatio->Reset();
381   TH1D* hKaonRatio   = (TH1D*)hPionRatio->Clone("hKaonRatio");
382   TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
383   TH1D* hPionRatioMc = (TH1D*)hPionRatio->Clone("hPionRatioMc");
384   TH1D* hKaonRatioMc = (TH1D*)hPionRatio->Clone("hKaonRatioMc");
385   TH1D* hProtonRatioMc = (TH1D*)hPionRatio->Clone("hProtonRatioMc");
386
387   TH1D* hPionYield =(TH1D*)hDeltaPiVsPt->ProjectionX("hPionYield", 1, 1);
388   hPionYield->SetTitle("particle fractions; p_{T} [GeV/c]; particle fraction");
389   hPionYield->Reset();
390   TH1D* hKaonYield   = (TH1D*)hPionYield->Clone("hKaonYield");
391   TH1D* hProtonYield = (TH1D*)hPionYield->Clone("hProtonYield");
392   TH1D* hPionYieldMc = (TH1D*)hPionYield->Clone("hPionYieldMc");
393   TH1D* hKaonYieldMc = (TH1D*)hPionYield->Clone("hKaonYieldMc");
394   TH1D* hProtonYieldMc = (TH1D*)hPionYield->Clone("hProtonYieldMc");
395   
396   for(Int_t bin = binStart; bin <= binStop; bin++){
397     
398     cout << "Making projection for bin: " << bin << endl;
399     
400     const Int_t j = bin-binStart;
401     
402     TH1D* hDeltaPiVsPtProj =(TH1D*)hDeltaPiVsPt->ProjectionY(Form("hDeltaPiVsPtProj%d", bin), bin, bin);
403     hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(-25, 20);
404     hDeltaPiVsPtProj->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c", 
405                                 hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(bin),
406                                 hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(bin)));
407
408     Double_t all =  hDeltaPiVsPtProj->GetEntries();
409
410     TH1D* hDeltaPiVsPtPiGenProj =(TH1D*)hDeltaPiVsPtPiGen->ProjectionY(Form("hDeltaPiVsPtPiGenProj%d", bin), bin, bin);
411     TH1D* hDeltaPiVsPtKGenProj =(TH1D*)hDeltaPiVsPtKGen->ProjectionY(Form("hDeltaPiVsPtKGenProj%d", bin), bin, bin);
412     TH1D* hDeltaPiVsPtPGenProj =(TH1D*)hDeltaPiVsPtPGen->ProjectionY(Form("hDeltaPiVsPtPGenProj%d", bin), bin, bin);
413     
414     Double_t gausParams[9] = { 
415       0.6*all,
416       hDeltaPiVsPtPiGenProj->GetMean(), 
417       hDeltaPiVsPtPiGenProj->GetRMS(), 
418       0.2*all,
419       hDeltaPiVsPtKGenProj->GetMean(), 
420       hDeltaPiVsPtKGenProj->GetRMS(), 
421       0.2*all,
422       hDeltaPiVsPtPGenProj->GetMean(), 
423       hDeltaPiVsPtPGenProj->GetRMS(), 
424     };
425
426     cFits->cd();
427     cFits->cd(j + 1);
428
429     total->SetParameters(gausParams);
430     for(Int_t i = 0; i < 9; i++) {
431
432       if((i%3) > 0)
433         total->FixParameter(i, gausParams[i]);
434       // else {
435       //        if(bin>56)
436       //          total->SetParLimits(i, 0, 10*all);
437       // }
438     }
439     
440     if(performFit) {
441
442       hDeltaPiVsPtProj->Fit(total, "0L");
443
444     } 
445
446     hDeltaPiVsPtProj->DrawCopy();
447     total->DrawCopy("same");    
448     
449     Double_t parametersOut[9];
450     total->GetParameters(parametersOut);
451     const Double_t* parameterErrorsOut = total->GetParErrors();
452
453     for(Int_t i = 0; i < 9; i++) 
454       cout << parametersOut[i] << ", ";
455     cout << endl;
456
457
458     pion->SetParameters(&parametersOut[0]);
459     pion->DrawCopy("same");
460     hPionRatio->SetBinContent(bin, parametersOut[0]/all);
461     hPionRatio->SetBinError(bin, parameterErrorsOut[0]/all);
462     hPionYield->SetBinContent(bin, parametersOut[0]);
463     hPionYield->SetBinError(bin, parameterErrorsOut[0]);
464
465     kaon->SetParameters(&parametersOut[3]);
466     kaon->DrawCopy("same");
467     hKaonRatio->SetBinContent(bin, parametersOut[3]/all);
468     hKaonRatio->SetBinError(bin, parameterErrorsOut[3]/all);
469     hKaonYield->SetBinContent(bin, parametersOut[3]);
470     hKaonYield->SetBinError(bin, parameterErrorsOut[3]);
471     
472     proton->SetParameters(&parametersOut[6]);
473     proton->DrawCopy("same");
474     hProtonRatio->SetBinContent(bin, parametersOut[6]/all);
475     hProtonRatio->SetBinError(bin, parameterErrorsOut[6]/all);
476     hProtonYield->SetBinContent(bin, parametersOut[6]);
477     hProtonYield->SetBinError(bin, parameterErrorsOut[6]);
478
479
480     cSingleFit->cd();
481     cSingleFit->Clear();
482     //    cSingleFit->SetLogy();
483     hDeltaPiVsPtProj->Draw();
484     pion->DrawCopy("same");
485     kaon->DrawCopy("same");
486     proton->DrawCopy("same");
487     total->DrawCopy("same");
488     
489     if(bin==42 || bin==49) {
490
491       TFile* fileOut = new TFile(Form("%s/ptspectrum_bin%d.root", dirName, bin), "RECREATE");
492
493       TH1D* hist = (TH1D*)hDeltaPiVsPtProj->Clone("hist");
494
495       TF1* sumFit    = (TF1*)total->Clone("sumFit");
496       TF1* pionFit   = (TF1*)pion->Clone("pionFit");
497       TF1* kaonFit   = (TF1*)kaon->Clone("kaonFit");
498       TF1* protonFit = (TF1*)proton->Clone("protonFit");
499
500       hist->Write();
501
502       sumFit->Write();
503       pionFit->Write();
504       kaonFit->Write();
505       protonFit->Write();
506
507       fileOut->Close();
508     }
509
510     gROOT->ProcessLine(".x drawText.C");
511     cSingleFit->SaveAs(Form("%s/ptspectrum_bin%d.gif", dirName, bin));
512     cSingleFit->SaveAs(Form("%s/ptspectrum_bin%d.pdf", dirName, bin));
513     //    legend->Draw();
514
515     if(data->IsMc()) {
516
517       cSingleFit->cd();
518       cSingleFit->Clear();
519       TH1D* hDeltaPiVsPtPiProj =(TH1D*)hDeltaPiVsPtPi->ProjectionY(Form("hDeltaPiVsPtPiProj%d", bin), bin, bin);
520       hDeltaPiVsPtPiProj->SetMarkerStyle(20);
521       hDeltaPiVsPtPiProj->SetMarkerColor(2);
522       hDeltaPiVsPtPiProj->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c", 
523                                         hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(bin),
524                                         hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(bin)));
525       hDeltaPiVsPtPiProj->Draw("P");
526       Double_t integral = hDeltaPiVsPtPiProj->Integral();
527       hPionRatioMc->SetBinContent(bin, integral/all);
528       hPionRatioMc->SetBinError(bin, TMath::Sqrt(integral)/all);
529       hPionYieldMc->SetBinContent(bin, integral);
530       hPionYieldMc->SetBinError(bin, TMath::Sqrt(integral));
531       TH1D* hDeltaPiVsPtKProj =(TH1D*)hDeltaPiVsPtK->ProjectionY(Form("hDeltaPiVsPtKProj%d", bin), bin, bin);
532       hDeltaPiVsPtKProj->SetMarkerStyle(21);
533       hDeltaPiVsPtKProj->SetMarkerColor(3);
534       hDeltaPiVsPtKProj->Draw("SAME P");
535       integral = hDeltaPiVsPtKProj->Integral();
536       hKaonRatioMc->SetBinContent(bin, integral/all);
537       hKaonRatioMc->SetBinError(bin, TMath::Sqrt(integral)/all);
538       hKaonYieldMc->SetBinContent(bin, integral);
539       hKaonYieldMc->SetBinError(bin, TMath::Sqrt(integral));
540       TH1D* hDeltaPiVsPtPProj =(TH1D*)hDeltaPiVsPtP->ProjectionY(Form("hDeltaPiVsPtPProj%d", bin), bin, bin);
541       hDeltaPiVsPtPProj->SetMarkerStyle(22);
542       hDeltaPiVsPtPProj->SetMarkerColor(4);
543       hDeltaPiVsPtPProj->Draw("SAME P");
544       integral = hDeltaPiVsPtPProj->Integral();
545       hProtonRatioMc->SetBinContent(bin, integral/all);
546       hProtonRatioMc->SetBinError(bin, TMath::Sqrt(integral)/all);
547       hProtonYieldMc->SetBinContent(bin, integral);
548       hProtonYieldMc->SetBinError(bin, TMath::Sqrt(integral));
549
550       pion->DrawCopy("same");
551       kaon->DrawCopy("same");
552       proton->DrawCopy("same");
553       cSingleFit->SaveAs(Form("debugfitsmc/ptspectrum_bin%d.gif", bin));
554     }
555   }
556
557   TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p", 600, 400);
558   cRatio->Clear();
559   hPionRatio->SetMarkerStyle(20);
560   hPionRatio->SetMarkerColor(2);
561   hPionRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
562   hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
563   hPionRatio->DrawCopy("P E");
564   hPionYield->SetMarkerStyle(20);
565   hPionYield->SetMarkerColor(2);
566   hPionYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
567
568   hKaonRatio->SetMarkerStyle(20);
569   hKaonRatio->SetMarkerColor(3);
570   hKaonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
571   hKaonRatio->DrawCopy("SAME P E");
572   hKaonYield->SetMarkerStyle(20);
573   hKaonYield->SetMarkerColor(3);
574   hKaonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
575
576   hProtonRatio->SetMarkerStyle(20);
577   hProtonRatio->SetMarkerColor(4);
578   hProtonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
579   hProtonRatio->DrawCopy("SAME P E");
580   hProtonYield->SetMarkerStyle(20);
581   hProtonYield->SetMarkerColor(4);
582   hProtonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
583
584   gROOT->ProcessLine(".x drawText.C");
585   cRatio->SaveAs(Form("%s/particle_ratios.gif", dirName));
586   cRatio->SaveAs(Form("%s/particle_ratios.pdf", dirName));
587
588   if(data->IsMc()) {
589     
590     hPionRatioMc->SetMarkerStyle(24);
591     hPionRatioMc->SetMarkerColor(2);
592     hPionRatioMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
593     hPionYieldMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
594     hPionYieldMc->SetMarkerStyle(24);
595     hPionYieldMc->SetMarkerColor(2);
596     hPionRatioMc->DrawCopy("SAME P");
597     
598     hKaonRatioMc->SetMarkerStyle(24);
599     hKaonRatioMc->SetMarkerColor(3);
600     hKaonRatioMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
601     hKaonYieldMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
602     hKaonYieldMc->SetMarkerStyle(24);
603     hKaonYieldMc->SetMarkerColor(3);
604     hKaonRatioMc->DrawCopy("SAME P");
605     
606     hProtonRatioMc->SetMarkerStyle(24);
607     hProtonRatioMc->SetMarkerColor(4);
608     hProtonRatioMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
609     hProtonYieldMc->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
610     hProtonYieldMc->SetMarkerStyle(24);
611     hProtonYieldMc->SetMarkerColor(4);
612     hProtonRatioMc->DrawCopy("SAME P");
613
614     cRatio->SaveAs("debugfitsmc/particle_ratios.gif");
615   }
616
617   if(outFileName) {
618
619     TFile* fileOut = new TFile(Form("fit_yields_results/%s", outFileName), "RECREATE");
620
621     // Histograms for normalization
622     data->GetHistVtxStatus()->Write();
623     TH1D* hPhiCutEff = data->GetHistNclVsPhiVsPtAfter()->ProjectionX("hPhiCutEff");
624     TH1D* hHelp = data->GetHistNclVsPhiVsPtBefore()->ProjectionX("hHelp");
625     hPhiCutEff->Divide(hPhiCutEff, hHelp, 1.0, 1.0, "B");
626     hPhiCutEff->Write();
627     delete hHelp;
628     delete hPhiCutEff;
629
630     hPionRatio->Write();
631     hKaonRatio->Write();
632     hProtonRatio->Write();
633
634     hPionYield->Write();
635     hKaonYield->Write();
636     hProtonYield->Write();
637
638     if(data->IsMc()) {
639       hPionRatioMc->Write();
640       hKaonRatioMc->Write();
641       hProtonRatioMc->Write();
642
643       hPionYieldMc->Write();
644       hKaonYieldMc->Write();
645       hProtonYieldMc->Write();                         
646     }
647
648     fileOut->Close();
649   }
650 }
651
652 //____________________________________________________________________________
653 void FitYieldsV0(const Char_t* dataFileName,
654                  const Char_t* v0FileName,
655                  const Char_t* v0Name,
656                  Double_t ptStart, Double_t ptStop,
657                  Int_t charge,
658                  Bool_t performFit = kFALSE,
659                  Int_t filter = 1,
660                  Int_t run    = 0,
661                  Bool_t usePhiCut = kTRUE,
662                  const Char_t* endName=0)
663 {
664   gStyle->SetOptStat(0);
665   
666   TFile* dataFile = FindFileFresh(dataFileName);
667   if(!dataFile)
668     return;
669   AliHighPtDeDxData* data = (AliHighPtDeDxData*)GetObject(dataFile, filter, usePhiCut, run, "filter", endName);
670   data->Print();
671
672   TFile* v0File = FindFileFresh(v0FileName);
673   if(!v0File)
674     return;
675   AliHighPtDeDxData* v0 = (AliHighPtDeDxData*)GetObject(v0File, 0, usePhiCut, run, v0Name, endName);
676   v0->Print();
677
678   TH2D* hDeltaPiVsPt = v0->GetHistDeltaPiVsPt(0, charge);
679   hDeDxVsP = hDeltaPiVsPt; // for the 2d fit to pick up the right bin
680
681   TH2D* hDeltaPiVsPtPiGen = data->GetHistDeltaPiVsPt(1, charge);
682   TH2D* hDeltaPiVsPtKGen  = data->GetHistDeltaPiVsPt(2, charge);
683   TH2D* hDeltaPiVsPtPGen  = data->GetHistDeltaPiVsPt(3, charge);
684
685   TProfile* hPiGenProfile = hDeltaPiVsPtPiGen->ProfileX();
686   hPiGenProfile->SetMarkerStyle(29);
687   TProfile* hKGenProfile = hDeltaPiVsPtKGen->ProfileX();
688   hKGenProfile->SetMarkerStyle(29);
689   TProfile* hPGenProfile = hDeltaPiVsPtPGen->ProfileX();
690   hPGenProfile->SetMarkerStyle(29);
691
692   TCanvas* cDeltaPiVsPt = new TCanvas("cDeltaPiVsPt", "dE/dx vs p", 400, 300);
693   cDeltaPiVsPt->Clear();
694   cDeltaPiVsPt->cd();
695   cDeltaPiVsPt->SetLogz();
696   hDeltaPiVsPt->Draw("COLZ");
697   hPiGenProfile->Draw("SAME P");
698   hKGenProfile->Draw("SAME P");
699   hPGenProfile->Draw("SAME P");
700
701   TCanvas* cDeltaPiVsPtLogX = new TCanvas("cDeltaPiVsPtLogX", "dE/dx vs p", 400, 300);
702   cDeltaPiVsPtLogX->Clear();
703   cDeltaPiVsPtLogX->cd();
704   cDeltaPiVsPtLogX->SetLogz();
705   cDeltaPiVsPtLogX->SetLogx();
706   hDeltaPiVsPt->Draw("COLZ");
707   hPiGenProfile->Draw("SAME P");
708   hKGenProfile->Draw("SAME P");
709   hPGenProfile->Draw("SAME P");
710
711   // Root is a bit stupid with finidng bins so we have to add and subtract a
712   // little to be sure we get the right bin as we typically put edges as
713   // limits
714   const Int_t binStart = hDeltaPiVsPt->GetXaxis()->FindBin(ptStart+0.01);
715   ptStart = hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(binStart);
716   const Int_t binStop  = hDeltaPiVsPt->GetXaxis()->FindBin(ptStop-0.01);
717   ptStop = hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(binStop);
718   //  const Int_t nBins    = binStop - binStart + 1;
719
720   cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
721        << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
722   
723
724   //cross check
725   TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
726   cFits->Clear();
727   cFits->Divide(7, 4);
728
729   TF1* pion = new TF1("pion", "[0]/sqrt(6.2832*[2]*[2])*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", -30, 20);
730     //  TF1* pion = new TF1("pion", "gausn", -30, 20);
731   pion->SetLineWidth(2);
732   pion->SetLineColor(kRed);
733   // TF1* kaon = new TF1("kaon", "gausn", -30, 20);
734   // kaon->SetLineWidth(2);
735   // kaon->SetLineColor(kGreen);
736   TF1* proton = new TF1("proton", "[0]/sqrt(6.2832*[2]*[2])*exp(-(x-[1])*(x-[1])/2.0/[2]/[2])", -30, 20);
737   //  TF1* proton = new TF1("proton", "gausn", -30, 20);
738   proton->SetLineWidth(2);
739   proton->SetLineColor(kBlue);
740
741   TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);    
742   legend->SetBorderSize(0);
743   legend->SetFillColor(0);
744   legend->AddEntry(pion, "#pi", "L");
745   //  legend->AddEntry(kaon, "K", "L");
746   legend->AddEntry(proton, "p", "L");
747
748   if(strstr(v0Name, "lambda")) {
749     gSystem->Exec("mv debugfitslambda/* olddebugfitslambda/");
750   } else {
751     gSystem->Exec("mv debugfitskaon/* olddebugfitskaon/");
752   }
753   TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
754
755   for(Int_t bin = binStart; bin <= binStop; bin++){
756     
757     cout << "Making projection for bin: " << bin << endl;
758     
759     const Int_t j = bin-binStart;
760     
761     TH1D* hDeltaPiVsPtProj =(TH1D*)hDeltaPiVsPt->ProjectionY(Form("hDeltaPiVsPtProj%d", bin), bin, bin);
762     //    hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40, 85);
763     hDeltaPiVsPtProj->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c", 
764                                 hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(bin),
765                                 hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(bin)));
766
767     Double_t all =  hDeltaPiVsPtProj->GetEntries();
768
769     TH1D* hDeltaPiVsPtPiGenProj =(TH1D*)hDeltaPiVsPtPiGen->ProjectionY(Form("hDeltaPiVsPtPiGenProj%d", bin), bin, bin);
770     //    TH1D* hDeltaPiVsPtKGenProj =(TH1D*)hDeltaPiVsPtKGen->ProjectionY(Form("hDeltaPiVsPtKGenProj%d", bin), bin, bin);
771     TH1D* hDeltaPiVsPtPGenProj =(TH1D*)hDeltaPiVsPtPGen->ProjectionY(Form("hDeltaPiVsPtPGenProj%d", bin), bin, bin);
772     
773     cFits->cd();
774     cFits->cd(j + 1);
775
776     hDeltaPiVsPtProj->Draw();
777
778     if(strstr(v0Name, "lambda")) {
779     
780       Double_t mean  = hDeltaPiVsPtPGenProj->GetMean();
781       Double_t sigma = hDeltaPiVsPtPGenProj->GetRMS();
782
783       cout << "Bin: " << bin 
784            <<   " (" << hDeltaPiVsPtProj->GetTitle()                            
785            << "), sigma: " << sigma << endl;
786       proton->SetParameters(all, mean, sigma);
787       proton->FixParameter(1, mean);
788       proton->FixParameter(2, sigma);
789
790       if(performFit) {
791         
792         hDeltaPiVsPtProj->Fit(proton, "L0", "", mean-sigma, mean+sigma);
793         
794         cout << "Bin: " << bin 
795              <<   " (" << hDeltaPiVsPtProj->GetTitle()                          
796              << "), sigma: " << proton->GetParameter(2) << endl;
797       } 
798
799       proton->DrawCopy("same");
800
801       cSingleFit->cd();
802       cSingleFit->Clear();
803       hDeltaPiVsPtProj->Draw();
804       proton->DrawCopy("same");
805     } else {
806
807       Double_t mean  = hDeltaPiVsPtPiGenProj->GetMean();
808       Double_t sigma = hDeltaPiVsPtPiGenProj->GetRMS();
809       cout << "Bin: " << bin 
810            <<   " (" << hDeltaPiVsPtProj->GetTitle()                            
811            << "), sigma: " << sigma << endl;
812       pion->SetParameters(all, mean, sigma);
813       pion->FixParameter(1, mean);
814       pion->FixParameter(2, sigma);
815       
816       if(performFit) {
817         
818         hDeltaPiVsPtProj->Fit(pion, "L0", "", mean-sigma, mean+sigma);
819
820         cout << "Bin: " << bin 
821              <<   " (" << hDeltaPiVsPtProj->GetTitle()                          
822              << "), sigma: " << pion->GetParameter(2) << endl;  
823       } 
824       pion->DrawCopy("same");
825
826       cSingleFit->cd();
827       cSingleFit->Clear();
828       hDeltaPiVsPtProj->Draw();
829       pion->DrawCopy("same");
830     }
831     
832     //    cSingleFit->SetLogy();
833     
834     gROOT->ProcessLine(".x drawText.C");
835
836     if(strstr(v0Name, "lambda")) {
837
838       if(!endName) {
839         cSingleFit->SaveAs(Form("results/lambda/ptspectrum_bin%d.gif", bin));
840         cSingleFit->SaveAs(Form("results/lambda/ptspectrum_bin%d.pdf", bin));
841       } else {
842         cSingleFit->SaveAs(Form("results/lambda%s/ptspectrum_bin%d.gif", endName, bin));
843         cSingleFit->SaveAs(Form("results/lambda%s/ptspectrum_bin%d.pdf", endName, bin));
844       }
845     } else {
846       if(!endName) {
847         cSingleFit->SaveAs(Form("results/kaon/ptspectrum_bin%d.gif", bin));
848         cSingleFit->SaveAs(Form("results/kaon/ptspectrum_bin%d.pdf", bin));
849       }else {
850         cSingleFit->SaveAs(Form("results/kaon%s/ptspectrum_bin%d.gif", endName, bin));
851         cSingleFit->SaveAs(Form("results/kaon%s/ptspectrum_bin%d.pdf", endName, bin));
852       }
853     }
854     //    legend->Draw();
855   }
856
857 }
858
859 //____________________________________________________________________________
860 void CompareYields(const Char_t* dataFileName1,
861                    const Char_t* dataFileName2,
862                    Double_t ptStart, Double_t ptStop,
863                    Int_t charge1,
864                    Int_t charge2,
865                    const Char_t* endName1=0,
866                    const Char_t* endName2=0,
867                    Bool_t performFit = kFALSE,
868                    Int_t run    = 0,
869                    Int_t filter = 1,
870                    Bool_t usePhiCut = kTRUE)
871 {
872   gStyle->SetOptStat(0);
873
874   
875   TFile* dataFile1 = FindFileFresh(dataFileName1);
876   if(!dataFile1)
877     return;
878   AliHighPtDeDxData* data1 = (AliHighPtDeDxData*)GetObject(dataFile1, filter, usePhiCut, run, "filter", endName1);
879   data1->Print();
880
881   gSystem->Exec("mv debugfits/* olddebugfits/");
882   gSystem->Exec("mv debugfitsratio/* olddebugfitsratio/");
883   // if(data1->IsMc())
884   //   gSystem->Exec("mv debugfitsmc/* olddebugfitsmc/");
885
886
887   TH2D* hDeltaPiVsPt1 = data1->GetHistDeltaPiVsPt(0, charge1);
888   hDeDxVsP = hDeltaPiVsPt1; // for the 2d fit to pick up the right bin
889
890   TH2D* hDeltaPiVsPtPiGen1 = data1->GetHistDeltaPiVsPt(1, charge1);
891   TH2D* hDeltaPiVsPtKGen1  = data1->GetHistDeltaPiVsPt(2, charge1);
892   TH2D* hDeltaPiVsPtPGen1  = data1->GetHistDeltaPiVsPt(3, charge1);
893
894   // TH2D* hDeltaPiVsPtPi1 = 0;
895   // TH2D* hDeltaPiVsPtK1  = 0;
896   // TH2D* hDeltaPiVsPtP1  = 0;
897
898   // if(data1->IsMc()) {
899
900   //   hDeltaPiVsPtPi1 = data1->GetHistDeltaPiVsPtPiMc();
901   //   hDeltaPiVsPtK1  = data1->GetHistDeltaPiVsPtKMc();
902   //   hDeltaPiVsPtP1  = data1->GetHistDeltaPiVsPtPMc();
903   // }
904
905   AliHighPtDeDxData* data2 = data1;
906   if(dataFileName2) {
907     TFile* dataFile2 = FindFileFresh(dataFileName2);
908     if(!dataFile2)
909       return;
910     data2 = (AliHighPtDeDxData*)GetObject(dataFile2, filter, usePhiCut, run, "filter", endName2);
911     data2->Print();
912   } else if (endName2) {
913
914     data2 = (AliHighPtDeDxData*)GetObject(dataFile1, filter, usePhiCut, run, "filter", endName2);
915   }
916
917   TH2D* hDeltaPiVsPt2 = data2->GetHistDeltaPiVsPt(0, charge2);
918   hDeDxVsP = hDeltaPiVsPt2; // for the 2d fit to pick up the right bin
919
920   // TH2D* hDeltaPiVsPtPiGen2 = data2->GetHistDeltaPiVsPt(1, charge2);
921   // TH2D* hDeltaPiVsPtKGen2  = data2->GetHistDeltaPiVsPt(2, charge2);
922   // TH2D* hDeltaPiVsPtPGen2  = data2->GetHistDeltaPiVsPt(3, charge2);
923
924   // TH2D* hDeltaPiVsPtPi2 = 0;
925   // TH2D* hDeltaPiVsPtK2  = 0;
926   // TH2D* hDeltaPiVsPtP2  = 0;
927
928   // if(data2->IsMc()) {
929
930   //   hDeltaPiVsPtPi2 = data2->GetHistDeltaPiVsPtPiMc();
931   //   hDeltaPiVsPtK2  = data2->GetHistDeltaPiVsPtKMc();
932   //   hDeltaPiVsPtP2  = data2->GetHistDeltaPiVsPtPMc();
933   // }
934
935
936   // Root is a bit stupid with finidng bins so we have to add and subtract a
937   // little to be sure we get the right bin as we typically put edges as
938   // limits
939   const Int_t binStart = hDeltaPiVsPt1->GetXaxis()->FindBin(ptStart+0.01);
940   ptStart = hDeltaPiVsPt1->GetXaxis()->GetBinLowEdge(binStart);
941   const Int_t binStop  = hDeltaPiVsPt1->GetXaxis()->FindBin(ptStop-0.01);
942   ptStop = hDeltaPiVsPt1->GetXaxis()->GetBinUpEdge(binStop);
943   //  const Int_t nBins    = binStop - binStart + 1;
944
945   cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
946        << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
947   
948
949   //cross check
950   TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
951   cFits->Clear();
952   cFits->Divide(7, 4);
953
954   TF1* pion = new TF1("pion", "gausn", -30, 20);
955   pion->SetLineWidth(2);
956   pion->SetLineColor(kRed);
957   TF1* kaon = new TF1("kaon", "gausn", -30, 20);
958   kaon->SetLineWidth(2);
959   kaon->SetLineColor(kGreen);
960   TF1* proton = new TF1("proton", "gausn", -30, 20);
961   proton->SetLineWidth(2);
962   proton->SetLineColor(kBlue);
963   TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", -30, 20);
964   total->SetLineColor(kBlack);
965   total->SetLineWidth(2);
966   total->SetLineStyle(2);
967
968   TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);    
969   legend->SetBorderSize(0);
970   legend->SetFillColor(0);
971   legend->AddEntry(total, "3-Gauss fit", "L");
972   legend->AddEntry(pion, "#pi", "L");
973   legend->AddEntry(kaon, "K", "L");
974   legend->AddEntry(proton, "p", "L");
975
976   TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
977
978   for(Int_t bin = binStart; bin <= binStop; bin++){
979     
980     cout << "Making projection for bin: " << bin << endl;
981     
982     const Int_t j = bin-binStart;
983     
984     TH1D* hDeltaPiVsPtProj1 =(TH1D*)hDeltaPiVsPt1->ProjectionY(Form("hDeltaPiVsPtProj1%d", bin), bin, bin);
985     hDeltaPiVsPtProj1->GetXaxis()->SetRangeUser(-25, 20);
986     hDeltaPiVsPtProj1->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c", 
987                                 hDeltaPiVsPt1->GetXaxis()->GetBinLowEdge(bin),
988                                 hDeltaPiVsPt1->GetXaxis()->GetBinUpEdge(bin)));
989
990     TH1D* hDeltaPiVsPtProj2 =(TH1D*)hDeltaPiVsPt2->ProjectionY(Form("hDeltaPiVsPtProj2%d", bin), bin, bin);
991
992     Double_t all1 =  hDeltaPiVsPtProj1->GetEntries();
993
994     TH1D* hDeltaPiVsPtPiGenProj1 =(TH1D*)hDeltaPiVsPtPiGen1->ProjectionY(Form("hDeltaPiVsPtPiGenProj1%d", bin), bin, bin);
995     TH1D* hDeltaPiVsPtKGenProj1 =(TH1D*)hDeltaPiVsPtKGen1->ProjectionY(Form("hDeltaPiVsPtKGenProj1%d", bin), bin, bin);
996     TH1D* hDeltaPiVsPtPGenProj1 =(TH1D*)hDeltaPiVsPtPGen1->ProjectionY(Form("hDeltaPiVsPtPGenProj1%d", bin), bin, bin);
997     
998     Double_t gausParams[9] = { 
999       0.6*all1,
1000       hDeltaPiVsPtPiGenProj1->GetMean(), 
1001       hDeltaPiVsPtPiGenProj1->GetRMS(), 
1002       0.2*all1,
1003       hDeltaPiVsPtKGenProj1->GetMean(), 
1004       hDeltaPiVsPtKGenProj1->GetRMS(), 
1005       0.2*all1,
1006       hDeltaPiVsPtPGenProj1->GetMean(), 
1007       hDeltaPiVsPtPGenProj1->GetRMS(), 
1008     };
1009
1010     cFits->cd();
1011     cFits->cd(j + 1);
1012
1013     total->SetParameters(gausParams);
1014     for(Int_t i = 0; i < 9; i++) {
1015
1016       if((i%3) > 0)
1017         total->FixParameter(i, gausParams[i]);
1018     }
1019     
1020     if(performFit) {
1021
1022       hDeltaPiVsPtProj1->Fit(total, "0L");
1023
1024     } 
1025
1026     hDeltaPiVsPtProj1->SetLineColor(4);
1027     hDeltaPiVsPtProj2->SetLineColor(2);
1028     hDeltaPiVsPtProj1->DrawCopy();
1029     hDeltaPiVsPtProj2->DrawCopy("SAME");
1030     if(performFit)
1031       total->DrawCopy("same");    
1032     
1033     Double_t parametersOut[9];
1034     total->GetParameters(parametersOut);
1035     //    const Double_t* parameterErrorsOut = total->GetParErrors();
1036
1037     for(Int_t i = 0; i < 9; i++) 
1038       cout << parametersOut[i] << ", ";
1039     cout << endl;
1040
1041
1042     if(performFit) {
1043       pion->SetParameters(&parametersOut[0]);
1044       pion->DrawCopy("same");
1045       
1046       kaon->SetParameters(&parametersOut[3]);
1047       kaon->DrawCopy("same");
1048       
1049       proton->SetParameters(&parametersOut[6]);
1050       proton->DrawCopy("same");
1051     }
1052
1053     cSingleFit->cd();
1054     cSingleFit->Clear();
1055     //    cSingleFit->SetLogy();
1056     hDeltaPiVsPtProj1->DrawCopy();
1057     hDeltaPiVsPtProj2->DrawCopy("SAME");
1058     if(performFit) {
1059       pion->DrawCopy("same");
1060       kaon->DrawCopy("same");
1061       proton->DrawCopy("same");
1062       total->DrawCopy("same");
1063     }
1064
1065     cSingleFit->SaveAs(Form("debugfits/ptspectrum_bin%d.gif", bin));
1066
1067     cSingleFit->cd();
1068     cSingleFit->Clear();
1069     //    cSingleFit->SetLogy();
1070     hDeltaPiVsPtProj1->Divide(hDeltaPiVsPtProj2);
1071     hDeltaPiVsPtProj1->GetYaxis()->SetRangeUser(0.8, 1.2);
1072     hDeltaPiVsPtProj1->DrawCopy();
1073
1074     cSingleFit->SaveAs(Form("debugfitsratio/ptspectrum_bin%d.gif", bin));
1075     //    legend->Draw();
1076
1077     // if(data->IsMc()) {
1078
1079     //   cSingleFit->cd();
1080     //   cSingleFit->Clear();
1081     //   TH1D* hDeltaPiVsPtPiProj =(TH1D*)hDeltaPiVsPtPi->ProjectionY(Form("hDeltaPiVsPtPiProj%d", bin), bin, bin);
1082     //   hDeltaPiVsPtPiProj->SetMarkerStyle(20);
1083     //   hDeltaPiVsPtPiProj->SetMarkerColor(2);
1084     //   hDeltaPiVsPtPiProj->SetTitle(Form("%.1f<p_{T}<%.1f GeV/c", 
1085     //                                  hDeltaPiVsPt->GetXaxis()->GetBinLowEdge(bin),
1086     //                                  hDeltaPiVsPt->GetXaxis()->GetBinUpEdge(bin)));
1087     //   hDeltaPiVsPtPiProj->Draw("P");
1088     //   hPionRatioMc->SetBinContent(bin, hDeltaPiVsPtPiProj->Integral()/all);
1089     //   hPionYieldMc->SetBinContent(bin, hDeltaPiVsPtPiProj->Integral());
1090     //   TH1D* hDeltaPiVsPtKProj =(TH1D*)hDeltaPiVsPtK->ProjectionY(Form("hDeltaPiVsPtKProj%d", bin), bin, bin);
1091     //   hDeltaPiVsPtKProj->SetMarkerStyle(21);
1092     //   hDeltaPiVsPtKProj->SetMarkerColor(3);
1093     //   hDeltaPiVsPtKProj->Draw("SAME P");
1094     //   hKaonRatioMc->SetBinContent(bin, hDeltaPiVsPtKProj->Integral()/all);
1095     //   hKaonYieldMc->SetBinContent(bin, hDeltaPiVsPtKProj->Integral());
1096     //   TH1D* hDeltaPiVsPtPProj =(TH1D*)hDeltaPiVsPtP->ProjectionY(Form("hDeltaPiVsPtPProj%d", bin), bin, bin);
1097     //   hDeltaPiVsPtPProj->SetMarkerStyle(22);
1098     //   hDeltaPiVsPtPProj->SetMarkerColor(4);
1099     //   hDeltaPiVsPtPProj->Draw("SAME P");
1100     //   hProtonRatioMc->SetBinContent(bin, hDeltaPiVsPtPProj->Integral()/all);
1101     //   hProtonYieldMc->SetBinContent(bin, hDeltaPiVsPtPProj->Integral());
1102
1103     //   pion->DrawCopy("same");
1104     //   kaon->DrawCopy("same");
1105     //   proton->DrawCopy("same");
1106     //   cSingleFit->SaveAs(Form("debugfitsmc/ptspectrum_bin%d.gif", bin));
1107     // }
1108   }
1109 }
1110
1111 //___________________________________________________________________________________________
1112 void MakeRatios(const Char_t* file1Name, const Char_t* file2Name, 
1113                 Bool_t drawFractionRatios,
1114                 Bool_t drawYieldRatios)
1115 {
1116   /*
1117     For yields we assume that file 1 is negative charge and file
1118     2 is positive charge.
1119    */
1120   
1121   TFile* file1 = FindFileFresh(file1Name);
1122   if(!file1)
1123     return;
1124   
1125   TH1D* hPionRatio1   = (TH1D*)file1->Get("hPionRatio");
1126   TH1D* hKaonRatio1   = (TH1D*)file1->Get("hKaonRatio");
1127   TH1D* hProtonRatio1 = (TH1D*)file1->Get("hProtonRatio");
1128
1129   TH1D* hPionYield1   = (TH1D*)file1->Get("hPionYield");
1130   TH1D* hKaonYield1   = (TH1D*)file1->Get("hKaonYield");
1131   TH1D* hProtonYield1 = (TH1D*)file1->Get("hProtonYield");
1132
1133   TH1D* hPionRatioMc1   = (TH1D*)file1->Get("hPionRatioMc");
1134   TH1D* hKaonRatioMc1   = (TH1D*)file1->Get("hKaonRatioMc");
1135   TH1D* hProtonRatioMc1 = (TH1D*)file1->Get("hProtonRatioMc");
1136
1137   TH1D* hPionYieldMc1   = (TH1D*)file1->Get("hPionYieldMc");
1138   TH1D* hKaonYieldMc1   = (TH1D*)file1->Get("hKaonYieldMc");
1139   TH1D* hProtonYieldMc1 = (TH1D*)file1->Get("hProtonYieldMc");
1140
1141   TFile* file2 = FindFileFresh(file2Name);
1142   if(!file2)
1143     return;
1144   
1145   TH1D* hPionRatio2   = (TH1D*)file2->Get("hPionRatio");
1146   TH1D* hKaonRatio2   = (TH1D*)file2->Get("hKaonRatio");
1147   TH1D* hProtonRatio2 = (TH1D*)file2->Get("hProtonRatio");
1148
1149   TH1D* hPionYield2   = (TH1D*)file2->Get("hPionYield");
1150   TH1D* hKaonYield2   = (TH1D*)file2->Get("hKaonYield");
1151   TH1D* hProtonYield2 = (TH1D*)file2->Get("hProtonYield");
1152
1153   TH1D* hPionRatioMc2   = (TH1D*)file2->Get("hPionRatioMc");
1154   TH1D* hKaonRatioMc2   = (TH1D*)file2->Get("hKaonRatioMc");
1155   TH1D* hProtonRatioMc2 = (TH1D*)file2->Get("hProtonRatioMc");
1156
1157   TH1D* hPionYieldMc2   = (TH1D*)file2->Get("hPionYieldMc");
1158   TH1D* hKaonYieldMc2   = (TH1D*)file2->Get("hKaonYieldMc");
1159   TH1D* hProtonYieldMc2 = (TH1D*)file2->Get("hProtonYieldMc");
1160   
1161   hPionRatio1->Divide(hPionRatio2);
1162   hPionRatio1->GetYaxis()->SetRangeUser(0.8, 1.2);
1163   hPionRatio1->SetTitle("Xcheck: pion fraction ratios; p_{T} [GeV/c]; pion fraction ratio");
1164
1165   hKaonRatio1->Divide(hKaonRatio2);
1166   hKaonRatio1->GetYaxis()->SetRangeUser(0.8, 1.2);
1167   hKaonRatio1->SetTitle("Xcheck: kaon fraction ratios; p_{T} [GeV/c]; kaon fraction ratio");
1168
1169   hProtonRatio1->Divide(hProtonRatio2);
1170   hProtonRatio1->GetYaxis()->SetRangeUser(0.8, 1.2);
1171   hProtonRatio1->SetTitle("Xcheck: proton fraction ratios; p_{T} [GeV/c]; proton fraction ratio");
1172
1173   hPionYield1->Divide(hPionYield2);
1174   hPionYield1->GetYaxis()->SetRangeUser(0.8, 1.2);
1175   hPionYield1->SetTitle("#pi^{-}/#pi^{+} vs p_{T}; p_{T} [GeV/c]; #pi^{-}/#pi^{+}");
1176
1177   hKaonYield1->Divide(hKaonYield2);
1178   hKaonYield1->GetYaxis()->SetRangeUser(0.8, 1.2);
1179   hKaonYield1->SetTitle("K^{-}/K^{+} vs p_{T}; p_{T} [GeV/c]; K^{-}/K^{+}");
1180
1181   hProtonYield1->Divide(hProtonYield2);
1182   hProtonYield1->GetYaxis()->SetRangeUser(0.8, 1.2);
1183   hProtonYield1->SetTitle("#bar{p}/p vs p_{T}; p_{T} [GeV/c]; #bar{p}/p");
1184
1185   if(hPionRatioMc1) {
1186     hPionRatioMc1->Divide(hPionRatioMc2);
1187     hPionRatioMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1188     hPionRatioMc1->SetTitle("Xcheck: pion fraction ratios (MC TRUTH); p_{T} [GeV/c]; pion fraction ratio");
1189     
1190     hKaonRatioMc1->Divide(hKaonRatioMc2);
1191     hKaonRatioMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1192     hKaonRatioMc1->SetTitle("Xcheck: kaon fraction ratios (MC TRUTH); p_{T} [GeV/c]; kaon fraction ratio");
1193     
1194     hProtonRatioMc1->Divide(hProtonRatioMc2);
1195     hProtonRatioMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1196     hProtonRatioMc1->SetTitle("Xcheck: proton fraction ratios (MC TRUTH); p_{T} [GeV/c]; proton fraction ratio");
1197     
1198     hPionYieldMc1->Divide(hPionYieldMc2);
1199     hPionYieldMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1200     hPionYieldMc1->SetTitle("#pi^{-}/#pi^{+} vs p_{T} (MC TRUTH); p_{T} [GeV/c]; #pi^{-}/#pi^{+}");
1201     
1202     hKaonYieldMc1->Divide(hKaonYieldMc2);
1203     hKaonYieldMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1204     hKaonYieldMc1->SetTitle("K^{-}/K^{+} vs p_{T} (MC TRUTH); p_{T} [GeV/c]; K^{-}/K^{+}");
1205     
1206     hProtonYieldMc1->Divide(hProtonYieldMc2);
1207     hProtonYieldMc1->GetYaxis()->SetRangeUser(0.8, 1.2);
1208     hProtonYieldMc1->SetTitle("#bar{p}/p vs p_{T} (MC TRUTH); p_{T} [GeV/c]; #bar{p}/p");
1209   }
1210
1211   if(drawFractionRatios) {
1212     TCanvas* cPionFractionRatio = new TCanvas("cPionFractionRatio", "pion fraction ratio", 400, 300);
1213     cPionFractionRatio->Clear();
1214     cPionFractionRatio->SetGridy();
1215     cPionFractionRatio->cd();
1216     hPionRatio1->Draw();
1217     cPionFractionRatio->SaveAs("pion_frac_ratio.gif");
1218
1219     TCanvas* cKaonFractionRatio = new TCanvas("cKaonFractionRatio", "kaon fraction ratio", 400, 300);
1220     cKaonFractionRatio->Clear();
1221     cKaonFractionRatio->SetGridy();
1222     cKaonFractionRatio->cd();
1223     hKaonRatio1->Draw();
1224     cKaonFractionRatio->SaveAs("kaon_frac_ratio.gif");
1225
1226     TCanvas* cProtonFractionRatio = new TCanvas("cProtonFractionRatio", "proton fraction ratio", 400, 300);
1227     cProtonFractionRatio->Clear();
1228     cProtonFractionRatio->SetGridy();
1229     cProtonFractionRatio->cd();
1230     hProtonRatio1->Draw();
1231     cProtonFractionRatio->SaveAs("proton_frac_ratio.gif");
1232
1233     if(hPionRatioMc1) {
1234       TCanvas* cPionFractionRatioMc = new TCanvas("cPionFractionRatioMc", "pion fraction ratio", 400, 300);
1235       cPionFractionRatioMc->Clear();
1236       cPionFractionRatioMc->SetGridy();
1237       cPionFractionRatioMc->cd();
1238       hPionRatioMc1->Draw();
1239       cPionFractionRatioMc->SaveAs("pion_frac_ratio_mc.gif");
1240       
1241       TCanvas* cKaonFractionRatioMc = new TCanvas("cKaonFractionRatioMc", "kaon fraction ratio", 400, 300);
1242       cKaonFractionRatioMc->Clear();
1243       cKaonFractionRatioMc->SetGridy();
1244       cKaonFractionRatioMc->cd();
1245       hKaonRatioMc1->Draw();
1246       cKaonFractionRatioMc->SaveAs("kaon_frac_ratio_mc.gif");
1247       
1248       TCanvas* cProtonFractionRatioMc = new TCanvas("cProtonFractionRatioMc", "proton fraction ratio", 400, 300);
1249       cProtonFractionRatioMc->Clear();
1250       cProtonFractionRatioMc->SetGridy();
1251       cProtonFractionRatioMc->cd();
1252       hProtonRatioMc1->Draw();
1253       cProtonFractionRatioMc->SaveAs("proton_frac_ratio_mc.gif");
1254     }
1255   }
1256   
1257   if(drawYieldRatios) {
1258     TCanvas* cPionRatio = new TCanvas("cPionRatio", "pion yields ratio", 400, 300);
1259     cPionRatio->Clear();
1260     cPionRatio->cd();
1261     cPionRatio->SetGridy();
1262     hPionYield1->Draw();
1263     cPionRatio->SaveAs("pion_ratio.gif");
1264
1265     TCanvas* cKaonRatio = new TCanvas("cKaonRatio", "kaon yields ratio", 400, 300);
1266     cKaonRatio->Clear();
1267     cKaonRatio->cd();
1268     cKaonRatio->SetGridy();
1269     hKaonYield1->Draw();
1270     cKaonRatio->SaveAs("kaon_ratio.gif");
1271
1272     TCanvas* cProtonRatio = new TCanvas("cProtonRatio", "proton yields ratio", 400, 300);
1273     cProtonRatio->Clear();
1274     cProtonRatio->cd();
1275     cProtonRatio->SetGridy();
1276     hProtonYield1->Draw();
1277     cProtonRatio->SaveAs("proton_ratio.gif");
1278
1279     if(hPionRatioMc1) {
1280       
1281       TCanvas* cPionRatioMc = new TCanvas("cPionRatioMc", "pion yields ratio", 400, 300);
1282       cPionRatioMc->Clear();
1283       cPionRatioMc->cd();
1284       cPionRatioMc->SetGridy();
1285       hPionYieldMc1->Draw();
1286       cPionRatioMc->SaveAs("pion_ratio_mc.gif");
1287       
1288       TCanvas* cKaonRatioMc = new TCanvas("cKaonRatioMc", "kaon yields ratio", 400, 300);
1289       cKaonRatioMc->Clear();
1290       cKaonRatioMc->cd();
1291       cKaonRatioMc->SetGridy();
1292       hKaonYieldMc1->Draw();
1293       cKaonRatioMc->SaveAs("kaon_ratio_mc.gif");
1294       
1295       TCanvas* cProtonRatioMc = new TCanvas("cProtonRatioMc", "proton yields ratio", 400, 300);
1296       cProtonRatioMc->Clear();
1297       cProtonRatioMc->cd();
1298       cProtonRatioMc->SetGridy();
1299       hProtonYieldMc1->Draw();
1300       cProtonRatioMc->SaveAs("proton_ratio_mc.gif");
1301     }
1302   }
1303 }
1304
1305 //___________________________________________________________________________________________
1306 void Compare(const Char_t* file1Name, const Char_t* file2Name, const Char_t* file3Name, 
1307              const Char_t* legend2, const Char_t* legend3, const Char_t* outfilename)
1308 {
1309   /*
1310     filename1 is the default
1311    */
1312   
1313   TLegend* legend = new TLegend(0.11, 0.68, 0.35, 0.88);    
1314   legend->SetBorderSize(0);
1315   legend->SetFillColor(0);
1316
1317   TFile* file1 = FindFileFresh(file1Name);
1318   if(!file1)
1319     return;
1320   
1321   TH1D* hPionRatio1   = (TH1D*)file1->Get("hPionRatio");
1322   hPionRatio1->SetMarkerStyle(29);
1323   TH1D* hPionRatio1Clone = (TH1D*)hPionRatio1->Clone("hPionsRatio1Clone");
1324   hPionRatio1Clone->SetMarkerColor(1);
1325   legend->AddEntry(hPionRatio1Clone, "default", "P");
1326   TH1D* hKaonRatio1   = (TH1D*)file1->Get("hKaonRatio");
1327   hKaonRatio1->SetMarkerStyle(29);
1328   TH1D* hProtonRatio1 = (TH1D*)file1->Get("hProtonRatio");
1329   hProtonRatio1->SetMarkerStyle(29);
1330
1331   TFile* file2 = FindFileFresh(file2Name);
1332   if(!file2)
1333     return;
1334   
1335   TH1D* hPionRatio2   = (TH1D*)file2->Get("hPionRatio");
1336   hPionRatio2->SetMarkerStyle(20);
1337   TH1D* hPionRatio2Clone = (TH1D*)hPionRatio2->Clone("hPionsRatio2Clone");
1338   hPionRatio2Clone->SetMarkerColor(1);
1339   legend->AddEntry(hPionRatio2Clone, legend2, "P");
1340   TH1D* hKaonRatio2   = (TH1D*)file2->Get("hKaonRatio");
1341   hKaonRatio2->SetMarkerStyle(20);
1342   TH1D* hProtonRatio2 = (TH1D*)file2->Get("hProtonRatio");
1343   hProtonRatio2->SetMarkerStyle(20);
1344
1345   TFile* file3 = FindFileFresh(file3Name);
1346   if(!file3)
1347     return;
1348   
1349   TH1D* hPionRatio3   = (TH1D*)file3->Get("hPionRatio");
1350   hPionRatio3->SetMarkerStyle(24);
1351   TH1D* hPionRatio3Clone = (TH1D*)hPionRatio3->Clone("hPionsRatio3Clone");
1352   hPionRatio3Clone->SetMarkerColor(1);
1353   legend->AddEntry(hPionRatio3Clone, legend3, "P");
1354   TH1D* hKaonRatio3   = (TH1D*)file3->Get("hKaonRatio");
1355   hKaonRatio3->SetMarkerStyle(24);
1356   TH1D* hProtonRatio3 = (TH1D*)file3->Get("hProtonRatio");
1357   hProtonRatio3->SetMarkerStyle(24);
1358
1359   TCanvas* cRatios = new TCanvas("cRatios", "pion fraction ratio", 400, 300);
1360   cRatios->Clear();
1361   cRatios->SetGridy();
1362   cRatios->cd();
1363   hPionRatio1->Draw("EP");
1364   hKaonRatio1->Draw("SAME EP");
1365   hProtonRatio1->Draw("SAME EP");
1366   hPionRatio2->Draw("HIST SAME P");
1367   hKaonRatio2->Draw("HIST SAME P");
1368   hProtonRatio2->Draw("HIST SAME P");
1369   hPionRatio3->Draw("HIST SAME P");
1370   hKaonRatio3->Draw("HIST SAME P");
1371   hProtonRatio3->Draw("HIST SAME P");
1372   legend->Draw();
1373   gROOT->ProcessLine(".x drawText.C");
1374   cRatios->SaveAs(Form("%s.gif", outfilename));
1375   cRatios->SaveAs(Form("%s.pdf", outfilename));
1376   
1377 }
1378
1379 //___________________________________________________________________________________________
1380 void Compare(Double_t x)
1381 {
1382   /*
1383     filename1 is the default
1384    */
1385   gStyle->SetOptStat(0);
1386
1387   TCanvas* cRatios = new TCanvas("cRatios", "pion fraction ratio", 400, 300);
1388
1389   const Int_t nEtaBins = 8;
1390   Double_t etaLimits[nEtaBins+1] = {-0.8, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8};
1391
1392   TGraphErrors* graphPionFractions = new TGraphErrors(nEtaBins);
1393   graphPionFractions->SetMarkerStyle(29);
1394   graphPionFractions->SetMarkerColor(2);
1395
1396   TGraphErrors* graphKaonFractions = new TGraphErrors(nEtaBins);
1397   graphKaonFractions->SetMarkerStyle(29);
1398   graphKaonFractions->SetMarkerColor(3);
1399
1400   TGraphErrors* graphProtonFractions = new TGraphErrors(nEtaBins);
1401   graphProtonFractions->SetMarkerStyle(29);
1402   graphProtonFractions->SetMarkerColor(4);
1403
1404   TH1F* hHelp = new TH1F("hHelp", "particle fractions vs #eta; #eta; particle fraction",
1405                          nEtaBins, etaLimits[0], etaLimits[nEtaBins]);
1406   hHelp->SetMinimum(0.0);
1407   hHelp->SetMaximum(1.0);
1408   hHelp->SetDirectory(0);
1409
1410   for(Int_t i = 0; i < nEtaBins; i++) {
1411
1412     Double_t etaCenter = (etaLimits[i] + etaLimits[i+1])/2.0;
1413     Double_t etaWidth  = TMath::Abs(etaCenter - etaLimits[i]);
1414
1415     Int_t etaLow  = Int_t(TMath::Abs(etaLimits[i]*10.0));
1416     Int_t etaHigh = Int_t(TMath::Abs(etaLimits[i+1]*10.0));
1417
1418     cout << etaCenter << ", " << etaWidth << ", " << etaLow << ", " << etaHigh << endl;
1419
1420     TFile* file = FindFileFresh(Form("fit_yields_results/lhc10h_aod_all_eta_%d%d.root", 
1421                                      etaLow, etaHigh));
1422     if(!file)
1423       return;
1424     
1425     TH1D* hPionRatio   = (TH1D*)file->Get("hPionRatio");
1426
1427     Int_t bin = hPionRatio->FindBin(x);
1428     if(i==0)
1429       hHelp->SetTitle(Form("%s (%.1f<p_{T}<%.1f GeV/c)", 
1430                            hHelp->GetTitle(),
1431                            hPionRatio->GetXaxis()->GetBinLowEdge(bin),
1432                            hPionRatio->GetXaxis()->GetBinUpEdge(bin)));
1433
1434     
1435     graphPionFractions->SetPoint(i, etaCenter, hPionRatio->GetBinContent(bin)); 
1436     graphPionFractions->SetPointError(i, etaWidth, hPionRatio->GetBinError(bin)); 
1437     
1438     TH1D* hKaonRatio   = (TH1D*)file->Get("hKaonRatio");
1439     graphKaonFractions->SetPoint(i, etaCenter, hKaonRatio->GetBinContent(bin)); 
1440     graphKaonFractions->SetPointError(i, etaWidth, hKaonRatio->GetBinError(bin)); 
1441
1442     TH1D* hProtonRatio = (TH1D*)file->Get("hProtonRatio");
1443     graphProtonFractions->SetPoint(i, etaCenter, hProtonRatio->GetBinContent(bin)); 
1444     graphProtonFractions->SetPointError(i, etaWidth, hProtonRatio->GetBinError(bin));     
1445   }
1446   
1447
1448   cRatios->Clear();
1449   cRatios->cd();
1450   hHelp->DrawCopy();
1451   graphPionFractions->Draw("P");
1452   graphKaonFractions->Draw("P");
1453   graphProtonFractions->Draw("P");
1454   gROOT->ProcessLine(".x drawText.C");
1455   cRatios->SaveAs(Form("results/comparison/ratios_vs_eta_%.1f.gif", x));
1456   cRatios->SaveAs(Form("results/comparison/ratios_vs_eta_%.1f.pdf", x));
1457   
1458   delete hHelp;
1459 }