Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / fit_yields_old.C
CommitLineData
4ebdd20e 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
28using 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
157void MakeRatios(const Char_t* file1Name, const Char_t* file2Name,
158 Bool_t drawFractionRatios,
159 Bool_t drawYieldRatios);
160
161
162//___________________________________________________________________________________________
163void 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//____________________________________________________________________________
254void 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//____________________________________________________________________________
653void 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//____________________________________________________________________________
860void 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//___________________________________________________________________________________________
1112void 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//___________________________________________________________________________________________
1306void 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//___________________________________________________________________________________________
1380void 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}