Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / fit_yields_final.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
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
88void MakeRatios(const Char_t* file1Name, const Char_t* file2Name,
89 Bool_t drawFractionRatios,
90 Bool_t drawYieldRatios);
91
92
93//____________________________________________________________________________
94void 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//____________________________________________________________________________
561void 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//____________________________________________________________________________
772void 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//___________________________________________________________________________________________
1024void 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//___________________________________________________________________________________________
1218void 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//___________________________________________________________________________________________
1292void 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}