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