Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / calibrate_de_dx.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 <TString.h>
15#include <TSystem.h>
16#include <TLatex.h>
17#include <TF1.h>
18#include <TF2.h>
19#include <TFitResultPtr.h>
20#include <TFitResult.h>
21
22#include "AliHighPtDeDxCalib.h"
23
24#include "my_tools.C"
25#include "my_functions.C"
26
27#include <iostream>
28#include <fstream>
29#include <string>
30
31using namespace std;
32
33/*
34 Ideas to improve:
35 =================
36
37 Use the real mean p -> Might give some effect
38 Effect: Push the <dE/dx> down (let'a see if it is still needed in the fits)
39
40 Use the real <nCl>
41 Effect: Increase sigma for p<15. For p>15 seems to saturate.
42
43 To use:
44
45 =======
46 root is enough
47
48 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC")
49 gSystem->AddIncludePath("-I../lib")
50 gSystem->AddIncludePath("-I../grid")
51 gSystem->AddIncludePath("-I../macros")
52 gROOT->SetMacroPath(".:../macros:../grid:../lib/")
53 .L libMyDeDxAnalysis.so
54 .L my_functions.C+
55 .L my_tools.C+
56 .L DebugClasses.C+
57 .L calibrate_de_dx.C+
58
59 // Step 2: test calibrations done in step 1
60 DrawStep2("calib_eta/7tev_b.root", kFALSE);
61 DrawStep2("calib_eta/7tev_b.root", kTRUE);
62 TestResolutionVsDeDx("calib_eta/7tev_b.root");
63 TestResolutionVsEta("calib_eta/7tev_b.root", kTRUE);
64
65 // Here we want to see that the data is symmetric in eta and roughly symmetric for high vs low eta
66 CompareYields("calib_eta/7tev_b.root", 0, 3.0, 30.0, "eta-80", "eta08")
67 CompareYields("calib_eta/7tev_b.root", 0, 3.0, 30.0, "etaabs04", "etaabs48")
68
69 // Step 3: extract
70 FitDeDxVsP("calib_eta/7tev_b.root", 3.0, 30.0, 4, 2, kTRUE)
71
72
73 // Step 4: inspect the results in results/calibdedx and results/calibplots
74 // using e.g. gwenview. Event the fit did not converge it is still ok if it
75 // describes the data!
76
77
78
79 //
80 // Test
81 //
82 // Step 2: test calibrations done in step 1
83 DrawStep2("calib_eta/7tev_b_test.root", kFALSE);
84 DrawStep2("calib_eta/7tev_b_test.root", kTRUE);
85 TestResolutionVsDeDx("calib_eta/7tev_b_test.root");
86 TestResolutionVsEta("calib_eta/7tev_b_test.root", kTRUE);
87 // Step 3: extract
88 FitDeDxVsP("calib_eta/7tev_b_test.root", 2.0, 10.0, 4, 2, kTRUE)
89
90
91 */
92
93
94//____________________________________________________________________________
95void FitDeDxVsP(const Char_t* calibFileName,
96 Double_t pStart, Double_t pStop,
97 Int_t optionDeDx, Int_t optionSigma,
98 Bool_t performFit = kFALSE,
99 Int_t run = 0,
100 Int_t filter = 1,
101 Bool_t usePhiCut = kTRUE,
102 const Char_t* v0FileName=0,
103 Bool_t fixAllPar=kFALSE)
104{
105 gStyle->SetOptStat(0);
106
107 TString baseName(gSystem->BaseName(calibFileName));
108 baseName.ReplaceAll(".root", "");
109
110 TFile* calibFile = FindFileFresh(calibFileName);
111 if(!calibFile)
112 return;
113 AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run);
114 calib->Print();
115
116 fixMIP = calib->GetHistDeDx(kTRUE, 0)->GetMean();
117 fixPlateau = calib->GetHistDeDx(kFALSE, 0)->GetMean();
118 hDeDxVsP = calib->GetHistDeDxVsP(0);
119 hMeanP = calib->GetHistMeanP();
120
121 AliHighPtDeDxCalib* lambdaCalib = 0;
122 AliHighPtDeDxCalib* kaonCalib = 0;
123 TH2D* hDeDxVsPV0Lambda = 0;
124 TH2D* hDeDxVsPV0Kaon = 0;
125 if(v0FileName) {
126 TFile* v0File = FindFileFresh(v0FileName);
127 if(!v0File)
128 return;
129 lambdaCalib = (AliHighPtDeDxCalib*)GetObject(v0File, 0, usePhiCut, run, "lambda");
130 lambdaCalib->Print();
131 hDeDxVsPV0Lambda = lambdaCalib->GetHistDeDxVsP(0);
132 kaonCalib = (AliHighPtDeDxCalib*)GetObject(v0File, 0, usePhiCut, run, "kaon");
133 kaonCalib->Print();
134 hDeDxVsPV0Kaon = kaonCalib->GetHistDeDxVsP(0);
135 }
136
137 CreateDir("old");
138 gSystem->Exec("mv results/calibdedx/* old/");
139 if(calib->IsMc())
140 gSystem->Exec("mv debugmc/* old/");
141
142 if(hDeDxVsPV0Lambda) {
143 gSystem->Exec("mv debuglambda/* old/");
144 gSystem->Exec("mv debugkaon/* old/");
145 }
146
147 TH2D* hDeDxVsPPi = 0;
148 TH2D* hDeDxVsPK = 0;
149 TH2D* hDeDxVsPP = 0;
150 TH2D* hDeDxVsPMu = 0;
151
152 if(calib->IsMc()) {
153
154 hDeDxVsPPi = calib->GetHistDeDxVsP(1);
155 hDeDxVsPK = calib->GetHistDeDxVsP(2);
156 hDeDxVsPP = calib->GetHistDeDxVsP(3);
157 hDeDxVsPMu = calib->GetHistDeDxVsP(5);
158 }
159
160 TCanvas* cDeDxVsP = new TCanvas("cDeDxVsP", "dE/dx vs p", 400, 300);
161 cDeDxVsP->Clear();
162 cDeDxVsP->cd();
163 cDeDxVsP->SetLogz();
164 hDeDxVsP->SetTitle(0);
165 hDeDxVsP->Draw("COLZ");
166
167 TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "dE/dx vs p", 400, 300);
168 cDeDxVsPLogX->Clear();
169 cDeDxVsPLogX->cd();
170 cDeDxVsPLogX->SetLogz();
171 cDeDxVsPLogX->SetLogx();
172 hDeDxVsP->Draw("COLZ");
173
174 // Root is a bit stupid with finidng bins so we have to add and subtract a
175 // little to be sure we get the right bin as we typically put edges as
176 // limits
177 const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(pStart+0.01);
178 pStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
179 const Int_t binStop = hDeDxVsP->GetXaxis()->FindBin(pStop-0.01);
180 pStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
181 const Int_t nBins = binStop - binStart + 1;
182
183 cout << "Doing 2d fit from pTlow = " << pStart << " (bin: " << binStart
184 << ") to pThigh = " << pStop << " (bin: " << binStop << ")" << endl;
185
186 // Double_t binSize = (histdEdxvsP->GetXaxis()->GetXmax() - histdEdxvsP->GetXaxis()->GetXmin())/ histdEdxvsP->GetXaxis()->GetNbins();
187
188 Double_t parDeDx[3] = {0, 0, 0};
189 Double_t parSigma[3] = {0, 0, 0};
190
191 const Int_t nLocalParams = 3; // pi, K, p yields
192 Int_t nDeDxPar = 0;
193 Int_t nSigmaPar = 0;
194
195 switch(optionDeDx) {
196
197 case 1:
198 nDeDxPar = 2;
199 parDeDx[0] = 39.7;
200 parDeDx[1] = 6.3;
201 break;
202 case 2:
203 nDeDxPar = 1;
204 parDeDx[0] = 7.3;
205 break;
206 case 3:
207 nDeDxPar = 2;
208 parDeDx[0] = 6.85097;
209 parDeDx[1] = 29.5933;
210 break;
211 // case 4:
212 // nDeDxPar = 3;
213 // parDeDx[0] = 35.5471;
214 // parDeDx[1] = 6.85097;
215 // parDeDx[2] = 29.5933;
216 // break;
217 case 4:
218 nDeDxPar = 3;
219 parDeDx[0] = 35.6694;
220 parDeDx[1] = 6.80835;
221 parDeDx[2] = 7.6542;
222 break;
223 default:
224
225 cout << "optionDeDx does not support option: " << optionDeDx << endl;
226 return;
227 break;
228 }
229
230 switch(optionSigma) {
231
232 case 1:
233 nSigmaPar = 1;
234 parSigma[0] = 3.0;
235 break;
236 case 2:
237 nSigmaPar = 1;
238 parSigma[0] = 0.0655688;
239 break;
240 case 3:
241 nSigmaPar = 2;
242 parSigma[0] = 0.06;
243 parSigma[1] = -0.001;
244 case 4:
245 nSigmaPar = 2;
246 parSigma[0] = 0.06;
247 parSigma[1] = 1.0;
248 break;
249 case 5:
250 nSigmaPar = 2;
251 parSigma[0] = 0.06;
252 parSigma[1] = 1.0;
253 break;
254 default:
255
256 cout << "optionSigma does not support option: " << optionSigma << endl;
257 return;
258 break;
259 }
260
261 const Int_t nGlobalParams = 2 // binStart, binStop,
262 + 2 + nDeDxPar // optionDeDx, nDeDxPar, dedxpar0 ....
263 + 2 + nSigmaPar; // nSigmaPar, sigmapar0 .....
264
265 const Int_t nParams = nBins*nLocalParams + nGlobalParams;
266
267
268 TF2* fitAll = new TF2("fitAll", fitf3G, pStart, pStop, 30, 90, nParams);
269 Double_t parametersIn[nParams];
270
271 parametersIn[0] = binStart;
272 fitAll->SetParName(0,"binStart");
273 fitAll->FixParameter(0, parametersIn[0]);
274
275 parametersIn[1] = binStop;
276 fitAll->SetParName(1,"binStop");
277 fitAll->FixParameter(1, parametersIn[1]);
278
279 // dE/dx parameters
280 parametersIn[2] = nDeDxPar;
281 fitAll->SetParName(2,"nDeDxPar");
282 fitAll->FixParameter(2, parametersIn[2]);
283
284 parametersIn[3] = optionDeDx;
285 fitAll->SetParName(3,"optionDeDx");
286 fitAll->FixParameter(3, parametersIn[3]);
287
288 for(Int_t i = 0; i < nDeDxPar; i++) {
289
290 parametersIn[4+i] = parDeDx[i];
291 fitAll->SetParName(4+i,Form("dEdXpar%d", i));
292 // if(optionDeDx == 4 && i <2) {
293
294 // fitAll->FixParameter(4+i, parametersIn[4+i]);
295 // }
296 if(fixAllPar) {
297
298 fitAll->FixParameter(4+i, parametersIn[4+i]);
299 }
300 }
301
302 // sigma parameters
303 parametersIn[4+nDeDxPar] = nSigmaPar;
304 fitAll->SetParName(4+nDeDxPar,"nSigmaPar");
305 fitAll->FixParameter(4+nDeDxPar, parametersIn[4+nDeDxPar]);
306
307 parametersIn[5+nDeDxPar] = optionSigma;
308 fitAll->SetParName(5+nDeDxPar,"optionSigma");
309 fitAll->FixParameter(5+nDeDxPar, parametersIn[5+nDeDxPar]);
310
311 for(Int_t i = 0; i < nSigmaPar; i++) {
312
313 parametersIn[6+nDeDxPar+i] = parSigma[i];
314 fitAll->SetParName(6+nDeDxPar+i,Form("sigmapar%d", i));
315 if(fixAllPar) {
316
317 fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
318 }
319 }
320
321 // Initial yields
322
323 for(Int_t bin = binStart; bin <= binStop; bin++) {
324
325 TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY("hDeDxVsPProj", bin, bin);
326
327 const Int_t offset = nGlobalParams + (bin - binStart)*nLocalParams;
328 const Double_t all = hDeDxVsPProj->Integral();
329 parametersIn[offset + 0] = (all)*0.6;
330 parametersIn[offset + 1] = (all)*0.3;
331 parametersIn[offset + 2] = (all)*0.1;
332
333 fitAll->SetParName(offset + 0, Form("piYield%d", bin));
334 fitAll->SetParLimits(offset + 0, 0, 10*all);
335 fitAll->SetParName(offset + 1, Form("kYield%d", bin));
336 fitAll->SetParLimits(offset + 1, 0, 10*all);
337 fitAll->SetParName(offset + 2, Form("pYield%d", bin));
338 fitAll->SetParLimits(offset + 2, 0, 10*all);
339 // fitAll->SetParLimits(offset + 0, 0., histdEdxvsPpy->GetEntries());
340 // fitAll->SetParLimits(offset + 1, 0., histdEdxvsPpy->GetEntries());
341 // fitAll->SetParLimits(offset + 2, 0., histdEdxvsPpy->GetEntries());
342 }
343
344 fitAll->SetParameters(parametersIn);
345 fitAll->Print();
346
347 Bool_t converged = kFALSE;
348
349 if(performFit) {
350 for(Int_t i = 0; i < 4; i++) {
351 TFitResultPtr fitResultPtr = hDeDxVsP->Fit(fitAll, "0S", "", pStart+0.01, pStop-0.01);
352 if(!fitResultPtr->Status()) {
353 // if(!fitResultPtr->Status() && !fitResultPtr->CovMatrixStatus()) {
354
355 converged = kTRUE;
356 break;
357 }
358 }
359 }
360 // else we just draw how the results looks with the input parameters
361
362 Double_t parametersOut[nParams];
363 fitAll->GetParameters(parametersOut);
364 const Double_t* parameterErrorsOut = fitAll->GetParErrors();
365
366 // overlay the fitfunction
367
368
369 TF1* fDeDxPi = new TF1("fDeDxPi", FitFunc, 0, 50, nDeDxPar+1); // +1 because of option!
370 fDeDxPi->SetParameters(&parametersOut[3]);
371 fDeDxPi->SetLineWidth(2);
372 cDeDxVsP->cd();
373 fDeDxPi->Draw("SAME");
374 cDeDxVsPLogX->cd();
375 fDeDxPi->Draw("SAME");
376
377 TF1* fDeDxK = new TF1("fDeDxK", FitFunc, 0, 50, nDeDxPar+1);
378 fDeDxK->SetParameters(&parametersOut[3]);
379 fDeDxK->SetParameter(0, fDeDxK->GetParameter(0)+10);
380 fDeDxK->SetLineWidth(2);
381 cDeDxVsP->cd();
382 fDeDxK->Draw("SAME");
383 cDeDxVsPLogX->cd();
384 fDeDxK->Draw("SAME");
385
386 TF1* fDeDxP = new TF1("fDeDxP", FitFunc, 0, 50, nDeDxPar+1);
387 fDeDxP->SetParameters(&parametersOut[3]);
388 fDeDxP->SetParameter(0, fDeDxP->GetParameter(0)+20);
389 fDeDxP->SetLineWidth(2);
390 cDeDxVsP->cd();
391 fDeDxP->Draw("SAME");
392 cDeDxVsPLogX->cd();
393 fDeDxP->Draw("SAME");
394
395 TF1* fDeDxE = new TF1("fDeDxE", FitFunc, 0, 50, nDeDxPar+1);
396 fDeDxE->SetParameters(&parametersOut[3]);
397 fDeDxE->SetParameter(0, fDeDxE->GetParameter(0)+30);
398 fDeDxE->SetLineWidth(2);
399 cDeDxVsP->cd();
400 fDeDxE->Draw("SAME");
401 cDeDxVsPLogX->cd();
402 fDeDxE->Draw("SAME");
403
404 TF1* fSigma = new TF1("fSigma", SigmaFunc, 0, 50, nSigmaPar+1);
405 fSigma->SetParameters(&parametersOut[5 + nDeDxPar]);
406
407 // fitAll->Draw("same cont3");
408
409 CreateDir("results/calibdedx");
410
411 cDeDxVsP->cd();
412 cDeDxVsP->Modified();
413 cDeDxVsP->Update();
414 gROOT->ProcessLine(".x drawText.C");
415 cDeDxVsP->SaveAs("results/calibdedx/dedx_vs_p.gif");
416 cDeDxVsP->SaveAs("results/calibdedx/dedx_vs_p.pdf");
417
418 cDeDxVsPLogX->cd();
419 cDeDxVsPLogX->Modified();
420 cDeDxVsPLogX->Update();
421 gROOT->ProcessLine(".x drawText.C");
422 cDeDxVsPLogX->SaveAs("results/calibdedx/dedx_vs_p_logx.gif");
423 cDeDxVsPLogX->SaveAs("results/calibdedx/dedx_vs_p_logx.pdf");
424
425 //cross check
426 TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
427 cFits->Clear();
428 cFits->Divide(7, 5);
429
430 TF1* pion = new TF1("pion", "gausn", 30, 90);
431 pion->SetLineWidth(2);
432 pion->SetLineColor(kRed);
433 TF1* kaon = new TF1("kaon", "gausn", 30, 90);
434 kaon->SetLineWidth(2);
435 kaon->SetLineColor(kGreen);
436 TF1* proton = new TF1("proton", "gausn", 30, 90);
437 proton->SetLineWidth(2);
438 proton->SetLineColor(kBlue);
439 TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", 30, 90);
440 total->SetLineColor(kBlack);
441 total->SetLineWidth(2);
442 total->SetLineStyle(2);
443
444 TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
445 legend->SetBorderSize(0);
446 legend->SetFillColor(0);
447 legend->AddEntry(total, "3-Gauss fit", "L");
448 legend->AddEntry(pion, "#pi", "L");
449 legend->AddEntry(kaon, "K", "L");
450 legend->AddEntry(proton, "p", "L");
451
452
453 TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
454
455 TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
456 hPionRatio->Reset();
457 hPionRatio->GetXaxis()->SetRangeUser(pStart+0.001, pStop-0.001);
458 hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
459 hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
460 TH1D* hKaonRatio = (TH1D*)hPionRatio->Clone("hKaonRatio");
461 TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
462
463 TH1D* hPionRatioMc = (TH1D*)hPionRatio->Clone("hPionRatioMc");
464 TH1D* hKaonRatioMc = (TH1D*)hPionRatio->Clone("hKaonRatioMc");
465 TH1D* hProtonRatioMc = (TH1D*)hPionRatio->Clone("hProtonRatioMc");
466 TH1D* hMuonRatioMc = (TH1D*)hPionRatio->Clone("hMuonRatioMc");
467
468 for(Int_t bin = binStart; bin <= binStop; bin++){
469
470 cout << "Making projection for bin: " << bin << endl;
471
472 const Int_t j = bin-binStart;
473 cFits->cd();
474 cFits->cd(j + 1);
475
476 TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeDxVsPProj%d", bin), bin, bin);
477 hDeDxVsPProj->GetXaxis()->SetRangeUser(40, 90);
478 hDeDxVsPProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
479 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
480 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
481 hDeDxVsPProj->Draw();
482
483 const Int_t offset = nGlobalParams + j*nLocalParams;
484 const Double_t p = hDeDxVsP->GetXaxis()->GetBinCenter(bin);
485 const Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
486 const Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
487 const Double_t meanDeDxPi = fDeDxPi->Eval(p);
488 const Double_t meanDeDxK = fDeDxPi->Eval(pKeff);
489 const Double_t meanDeDxP = fDeDxPi->Eval(pPeff);
490 Double_t gausParams[9] = {
491 parametersOut[offset + 0],
492 meanDeDxPi,
493 fSigma->Eval(meanDeDxPi) ,
494 parametersOut[offset + 1],
495 meanDeDxK,
496 fSigma->Eval(meanDeDxK) ,
497 parametersOut[offset + 2],
498 meanDeDxP,
499 fSigma->Eval(meanDeDxP) ,
500 };
501
502 for(Int_t i = 0; i < 9; i++)
503 cout << gausParams[i] << ", ";
504
505 cout << endl;
506
507 pion->SetParameters(&gausParams[0]);
508 pion->DrawCopy("same");
509 Double_t all = hDeDxVsPProj->Integral();
510 hPionRatio->SetBinContent(bin, parametersOut[offset + 0]/all);
511 hPionRatio->SetBinError(bin, parameterErrorsOut[offset + 0]/all);
512
513 kaon->SetParameters(&gausParams[3]);
514 kaon->DrawCopy("same");
515 hKaonRatio->SetBinContent(bin, parametersOut[offset + 1]/all);
516 hKaonRatio->SetBinError(bin, parameterErrorsOut[offset + 1]/all);
517
518 proton->SetParameters(&gausParams[6]);
519 proton->DrawCopy("same");
520 hProtonRatio->SetBinContent(bin, parametersOut[offset + 2]/all);
521 hProtonRatio->SetBinError(bin, parameterErrorsOut[offset + 2]/all);
522
523 total->SetParameters(gausParams);
524 total->DrawCopy("same");
525
526 cSingleFit->cd();
527 cSingleFit->Clear();
528 // cSingleFit->SetLogy();
529 hDeDxVsPProj->Draw();
530 pion->DrawCopy("same");
531 kaon->DrawCopy("same");
532 proton->DrawCopy("same");
533 total->DrawCopy("same");
534
535 gROOT->ProcessLine(".x drawText.C(2)");
536 cSingleFit->SaveAs(Form("results/calibdedx/ptspectrum_bin%d.gif", bin));
537 cSingleFit->SaveAs(Form("results/calibdedx/ptspectrum_bin%d.pdf", bin));
538 // legend->Draw();
539
540 if(calib->IsMc()) {
541
542 cSingleFit->cd();
543 cSingleFit->Clear();
544 TH1D* hDeDxVsPPiProj =(TH1D*)hDeDxVsPPi->ProjectionY(Form("hDeDxVsPPiProj%d", bin), bin, bin);
545 hDeDxVsPPiProj->SetMarkerStyle(20);
546 hDeDxVsPPiProj->SetMarkerColor(2);
547 hDeDxVsPPiProj->GetXaxis()->SetRangeUser(40, 90);
548 hDeDxVsPPiProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
549 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
550 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
551 hDeDxVsPPiProj->Draw("P");
552 hPionRatioMc->SetBinContent(bin, hDeDxVsPPiProj->Integral()/all);
553 TH1D* hDeDxVsPKProj =(TH1D*)hDeDxVsPK->ProjectionY(Form("hDeDxVsPKProj%d", bin), bin, bin);
554 hDeDxVsPKProj->SetMarkerStyle(21);
555 hDeDxVsPKProj->SetMarkerColor(3);
556 hDeDxVsPKProj->Draw("SAME P");
557 hKaonRatioMc->SetBinContent(bin, hDeDxVsPKProj->Integral()/all);
558 TH1D* hDeDxVsPPProj =(TH1D*)hDeDxVsPP->ProjectionY(Form("hDeDxVsPPProj%d", bin), bin, bin);
559 hDeDxVsPPProj->SetMarkerStyle(22);
560 hDeDxVsPPProj->SetMarkerColor(4);
561 hDeDxVsPPProj->Draw("SAME P");
562 hProtonRatioMc->SetBinContent(bin, hDeDxVsPPProj->Integral()/all);
563 TH1D* hDeDxVsPMuProj =(TH1D*)hDeDxVsPMu->ProjectionY(Form("hDeDxVsPMuProj%d", bin), bin, bin);
564 hDeDxVsPMuProj->SetMarkerStyle(22);
565 hDeDxVsPMuProj->SetMarkerColor(6);
566 // hDeDxVsPMuProj->Draw("SAME P");
567 hMuonRatioMc->SetBinContent(bin, hDeDxVsPMuProj->Integral()/all);
568 // cSingleFit->SetLogy();
569 pion->DrawCopy("same");
570 kaon->DrawCopy("same");
571 proton->DrawCopy("same");
572 CreateDir("results/calibdedx/debugmc");
573 cSingleFit->SaveAs(Form("results/calibdedx/debugmc/ptspectrum_bin%d.gif", bin));
574 }
575
576 if(hDeDxVsPV0Lambda) {
577 cSingleFit->cd();
578 cSingleFit->Clear();
579 TH1D* hDeDxVsPV0LambdaProj =(TH1D*)hDeDxVsPV0Lambda->ProjectionY(Form("hDeDxVsPV0LambdaProj%d", bin), bin, bin);
580 hDeDxVsPV0LambdaProj->SetMarkerStyle(20);
581 hDeDxVsPV0LambdaProj->SetMarkerColor(4);
582 hDeDxVsPV0LambdaProj->GetXaxis()->SetRangeUser(40, 90);
583 hDeDxVsPV0LambdaProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
584 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
585 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
586 hDeDxVsPV0LambdaProj->Draw("P");
587 proton->SetParameter(0, hDeDxVsPV0LambdaProj->Integral());
588 proton->DrawCopy("same");
589
590 CreateDir("results/calibdedx/debuglambda");
591 cSingleFit->SaveAs(Form("results/calibdedx/debuglambda/ptspectrum_bin%d.gif", bin));
592 }
593
594 if(hDeDxVsPV0Kaon) {
595 cSingleFit->cd();
596 cSingleFit->Clear();
597 TH1D* hDeDxVsPV0KaonProj =(TH1D*)hDeDxVsPV0Kaon->ProjectionY(Form("hDeDxVsPV0KaonProj%d", bin), bin, bin);
598 hDeDxVsPV0KaonProj->SetMarkerStyle(20);
599 hDeDxVsPV0KaonProj->SetMarkerColor(2);
600 hDeDxVsPV0KaonProj->GetXaxis()->SetRangeUser(40, 90);
601 hDeDxVsPV0KaonProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
602 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
603 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
604 hDeDxVsPV0KaonProj->Draw("P");
605 pion->SetParameter(0, hDeDxVsPV0KaonProj->Integral());
606 pion->DrawCopy("same");
607 CreateDir("results/calibdedx/debugkaon");
608 cSingleFit->SaveAs(Form("results/calibdedx/debugkaon/ptspectrum_bin%d.gif", bin));
609 }
610 }
611
612 TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p", 600, 400);
613 cRatio->Clear();
614 hPionRatio->SetMaximum(0.8);
615 hPionRatio->SetMarkerStyle(20);
616 hPionRatio->SetMarkerColor(2);
617 hPionRatio->Draw("P E");
618
619 hKaonRatio->SetMarkerStyle(20);
620 hKaonRatio->SetMarkerColor(3);
621 hKaonRatio->Draw("SAME P E");
622
623 hProtonRatio->SetMarkerStyle(20);
624 hProtonRatio->SetMarkerColor(4);
625 hProtonRatio->Draw("SAME P E");
626 gROOT->ProcessLine(".x drawText.C(2)");
627 cRatio->SaveAs("results/calibdedx/particle_ratios.gif");
628 cRatio->SaveAs("results/calibdedx/particle_ratios.pdf");
629
630 if(calib->IsMc()) {
631
632 hPionRatioMc->SetMarkerStyle(24);
633 hPionRatioMc->SetMarkerColor(2);
634 hPionRatioMc->Draw("SAME P");
635
636 hKaonRatioMc->SetMarkerStyle(24);
637 hKaonRatioMc->SetMarkerColor(3);
638 hKaonRatioMc->Draw("SAME P");
639
640 hProtonRatioMc->SetMarkerStyle(24);
641 hProtonRatioMc->SetMarkerColor(4);
642 hProtonRatioMc->Draw("SAME P");
643
644 hMuonRatioMc->SetMarkerStyle(24);
645 hMuonRatioMc->SetMarkerColor(6);
646 // hMuonRatioMc->Draw("SAME P");
647 cRatio->SaveAs("results/calibdedx/debugmc/particle_ratios_mc.gif");
648 }
649
650
651 cout << "MIP was fixed to: " << fixMIP << endl
652 << "Plateau was fixed to: " << fixPlateau << endl;
653
654 //
655 // Store the <dE/dx> parameters in a file that we can get them back to use for the Delta-pi!
656 //
657 DeDxFitInfo* fitInfo = new DeDxFitInfo();
658 fitInfo->MIP = fixMIP;
659 fitInfo->plateau = fixPlateau;
660 fitInfo->optionDeDx = optionDeDx;
661 fitInfo->nDeDxPar = nDeDxPar;
662 for(Int_t i = 0; i < nDeDxPar; i++) {
663 fitInfo->parDeDx[i] = fDeDxPi->GetParameter(i+1); // 1st parameter is option
664 }
665 fitInfo->optionSigma = optionSigma;
666 fitInfo->nSigmaPar = nSigmaPar;
667 for(Int_t i = 0; i < nSigmaPar; i++) {
668 fitInfo->parSigma[i] = fSigma->GetParameter(i+1); // 1st parameter is option
669 }
670 if(!strstr(calibFileName, "no_calib_eta"))
671 fitInfo->calibFileName = calibFileName;
672
673 fitInfo->Print();
674
675 CreateDir("fitparameters");
676 TFile* outFile = new TFile(Form("fitparameters/%s", gSystem->BaseName(calibFileName)), "RECREATE");
677 fitInfo->Write("fitInfo");
678 outFile->Close();
679
680 if(converged) {
681
682 cout << "Fit converged and error matrix was ok" << endl;
683 } else {
684
685 cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was not ok!" << endl;
686 }
687}
688
689
690
691//____________________________________________________________________________
692void DrawStep2(const Char_t* calibFileName,
693 Bool_t forMIP = kTRUE,
694 Int_t run = 0,
695 Int_t filter = 1,
696 Bool_t usePhiCut = kTRUE)
697{
698 // option:
699 // 0 - compare to ALICE INEL
700 // 1 - compare to MC input spectrum
701 // 2 - compare to MC correct spectrum
702 gStyle->SetOptStat(0);
703
704 CreateDir("results/calibplots");
705 TString baseName(gSystem->BaseName(calibFileName));
706 baseName.ReplaceAll(".root", "");
707
708 TFile* calibFile = FindFileFresh(calibFileName);
709 if(!calibFile)
710 return;
711 AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run);
712 calib->Print();
713
714 TCanvas* cPhiCut = calib->DrawPhiCutHistograms();
715 gROOT->ProcessLine(".x drawText.C(2)");
716 cPhiCut->SaveAs(Form("results/calibplots/phicut_%s.gif", baseName.Data()));
717 cPhiCut->SaveAs(Form("results/calibplots/phicut_%s.pdf", baseName.Data()));
718 TCanvas* cEta = calib->DrawEta(forMIP);
719 gROOT->ProcessLine(".x drawText.C(2)");
720 TCanvas* cEtaCal = calib->DrawEtaCalibrated(forMIP);
721 gROOT->ProcessLine(".x drawText.C(2)");
722 if(forMIP) {
723 cEta->cd();
724 cEta->SaveAs(Form("results/calibplots/MIP_dedx_vs_eta_nocal_%s.gif", baseName.Data()));
725 cEta->SaveAs(Form("results/calibplots/MIP_dedx_vs_eta_nocal_%s.pdf", baseName.Data()));
726 cEtaCal->cd();
727 cEtaCal->SaveAs(Form("results/calibplots/MIP_dedx_vs_eta_final_%s.gif", baseName.Data()));
728 cEtaCal->SaveAs(Form("results/calibplots/MIP_dedx_vs_eta_final_%s.pdf", baseName.Data()));
729 } else {
730 cEta->cd();
731 cEta->SaveAs(Form("results/calibplots/Electrons_dedx_vs_eta_nocal_%s.gif", baseName.Data()));
732 cEta->SaveAs(Form("results/calibplots/Electrons_dedx_vs_eta_nocal_%s.pdf", baseName.Data()));
733 cEtaCal->cd();
734 cEtaCal->SaveAs(Form("results/calibplots/Electrons_dedx_vs_eta_final_%s.gif", baseName.Data()));
735 cEtaCal->SaveAs(Form("results/calibplots/Electrons_dedx_vs_eta_final_%s.pdf", baseName.Data()));
736 }
737 TCanvas* cEtaSel = calib->DrawSelectionHistograms();
738 gROOT->ProcessLine(".x drawText.C(2)");
739 cEtaSel->SaveAs(Form("results/calibplots/selection_%s.gif", baseName.Data()));
740 cEtaSel->SaveAs(Form("results/calibplots/selection_%s.pdf", baseName.Data()));
741
742 // TCanvas* cNcl = calib->DrawNclCal();
743 // gROOT->ProcessLine(".x drawText.C");
744 // cNcl->SaveAs(Form("results/calibplots/dedx_vs_ncl_calib__%s.gif", baseName.Data()));
745 // cNcl->SaveAs(Form("results/calibplots/dedx_vs_ncl_calib__%s.pdf", baseName.Data()));
746}
747
748//____________________________________________________________________________
749void TestResolutionVsDeDx(const Char_t* calibFileName,
750 Int_t run = 0,
751 Int_t filter = 1,
752 Bool_t usePhiCut = kTRUE)
753{
754 gStyle->SetOptStat(0);
755
756 TFile* calibFile = FindFileFresh(calibFileName);
757 if(!calibFile)
758 return;
759 AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run);
760 calib->Print();
761
762 TH2D* hDeDxVsNcl = calib->GetHistDeDxVsNcl(kTRUE, 0);
763 TH2D* hDeDxVsNclE = calib->GetHistDeDxVsNcl(kFALSE, 0);
764 TH1D* hDeDx = calib->GetHistDeDx(kTRUE, 0);
765 TH1D* hDeDxE = calib->GetHistDeDx(kFALSE, 0);
766
767 Double_t mdedx = hDeDx->GetMean();
768 Double_t mdedxE = hDeDxE->GetMean();
769
770 TObjArray* helpArray = new TObjArray();
771
772 hDeDxVsNcl->FitSlicesY(0, 0, -1, 0, "QNR", helpArray);
773 TH1D* hMeanVsNcl = (TH1D*)helpArray->At(1);
774 hMeanVsNcl->SetName("hMeanVsNcl");
775 hMeanVsNcl->SetMarkerStyle(29);
776 TH1D* hSigmaVsNcl = (TH1D*)helpArray->At(2);
777 hSigmaVsNcl->SetName("hSigmaVsNcl");
778 hSigmaVsNcl->SetTitle("#sigma vs Ncl; Ncl; #sigma");
779 hSigmaVsNcl->SetMarkerStyle(29);
780 hSigmaVsNcl->SetMarkerColor(2);
781 hSigmaVsNcl->Scale(1.0/mdedx);
782
783 hDeDxVsNclE->FitSlicesY(0, 0, -1, 0, "QNR", helpArray);
784 TH1D* hMeanVsNclE = (TH1D*)helpArray->At(1);
785 hMeanVsNclE->SetName("hSigmaVsNclE");
786 hMeanVsNclE->SetMarkerStyle(29);
787 TH1D* hSigmaVsNclE = (TH1D*)helpArray->At(2);
788 hSigmaVsNclE->SetName("hSigmaVsNclE");
789 hSigmaVsNclE->SetTitle("#sigma vs Ncl; Ncl; #sigma");
790 hSigmaVsNclE->SetMarkerStyle(29);
791 hSigmaVsNclE->SetMarkerColor(kMagenta);
792 hSigmaVsNclE->Scale(1.0/mdedxE);
793
794 TLegend* legend = new TLegend(0.54, 0.67, 0.84, 0.87);
795 legend->SetFillColor(0);
796 legend->SetBorderSize(0);
797 legend->SetTextSize(0.05);
798 legend->AddEntry(hSigmaVsNcl, "MIP pions", "P");
799 legend->AddEntry(hSigmaVsNclE, "electrons", "P");
800
801 TCanvas* cTestDeDx = new TCanvas("cTestDeDx", "Comparing resolution for MIPs and electrons",
802 1200, 400);
803 cTestDeDx->Clear();
804 cTestDeDx->Divide(3, 1);
805
806 cTestDeDx->cd(1);
807 hDeDxVsNcl->Draw("COL");
808 hMeanVsNcl->Draw("SAME");
809
810 cTestDeDx->cd(2);
811 hDeDxVsNclE->Draw("COL");
812 hMeanVsNclE->Draw("SAME");
813
814 cTestDeDx->cd(3);
815 TF1* fSDedxVsNclSqrt = new TF1("fSDedxVsNclSqrt", "[0]*sqrt(159.0/x)", 0, 160);
816 fSDedxVsNclSqrt->SetParameter(0, 0.05);
817 hSigmaVsNcl->SetMinimum(0.0);
818 hSigmaVsNcl->SetMaximum(0.15);
819 hSigmaVsNcl->SetTitle("#sigma/<dE/dx> vs Ncl; Ncl; rel. #sigma");
820 // hSigmaVsNcl->Fit(fSDedxVsNclSqrt);
821 hSigmaVsNcl->Draw();
822 hSigmaVsNclE->Draw("SAME");
823 legend->Draw();
824 gROOT->ProcessLine(".x drawText.C(2)");
825 TString baseName(gSystem->BaseName(calibFileName));
826 baseName.ReplaceAll(".root", "");
827 cTestDeDx->SaveAs(Form("results/calibplots/test_resolution_vs_dedx_%s.gif", baseName.Data()));
828 cTestDeDx->SaveAs(Form("results/calibplots/test_resolution_vs_dedx_%s.pdf", baseName.Data()));
829}
830
831//___________________________________________________________________________________________
832void TestResolutionVsEta(const Char_t* calibFileName,
833 Bool_t forMIP,
834 Int_t run = 0,
835 Int_t filter = 1,
836 Bool_t usePhiCut = kTRUE)
837{
838 gStyle->SetOptStat(0);
839
840 TFile* calibFile = FindFileFresh(calibFileName);
841 if(!calibFile)
842 return;
843 AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run);
844 calib->Print();
845
846 const Int_t nEta = 4;
847 Int_t color[nEta] = {kBlack, kBlue, kGreen+2, kRed};
848 TH1D* hDeDx[nEta] = {0, 0, 0, 0};
849 TH1D* hMeanVsNcl[nEta] = {0, 0, 0, 0};
850 TH1D* hSigmaVsNcl[nEta] = {0, 0, 0, 0};
851
852 const Double_t mdedx = (calib->GetHistDeDx(forMIP, 0))->GetMean();
853
854 TObjArray* helpArray = new TObjArray();
855
856 for(Int_t bin = 1; bin <= nEta; bin++) {
857
858 const Int_t i = bin - 1;
859 hDeDx[i] = calib->GetHistDeDx(forMIP, bin);
860 hDeDx[i]->Scale(1.0/hDeDx[i]->GetEntries());
861 hDeDx[i]->SetMarkerStyle(24+i);
862 hDeDx[i]->SetMarkerColor(color[i]);
863
864 TH2D* hDeDxVsNcl = calib->GetHistDeDxVsNcl(forMIP, bin);
865 hDeDxVsNcl->FitSlicesY(0, 0, -1, 0, "QNR", helpArray);
866
867 hMeanVsNcl[i] = (TH1D*)helpArray->At(1);
868 hMeanVsNcl[i]->SetTitle("<dE/dx> vs Ncl; Ncl; #sigma");
869 hMeanVsNcl[i]->SetName(Form("hMeanVsNcl%d", bin));
870 hMeanVsNcl[i]->SetMarkerStyle(24+i);
871 hMeanVsNcl[i]->SetMarkerColor(color[i]);
872 hMeanVsNcl[i]->SetMinimum(mdedx*0.9);
873 hMeanVsNcl[i]->SetMaximum(mdedx*1.1);
874
875 hSigmaVsNcl[i] = (TH1D*)helpArray->At(2);
876 hSigmaVsNcl[i]->SetName("hSigmaVsNcl");
877 hSigmaVsNcl[i]->SetTitle("#sigma/<dE/dx> vs Ncl; Ncl; rel. #sigma");
878 hSigmaVsNcl[i]->SetMarkerStyle(24+i);
879 hSigmaVsNcl[i]->SetMarkerColor(color[i]);
880 hSigmaVsNcl[i]->Scale(1.0/mdedx);
881 hSigmaVsNcl[i]->SetMinimum(0.0);
882 hSigmaVsNcl[i]->SetMaximum(0.15);
883 }
884
885 TLegend* legend = new TLegend(0.54, 0.67, 0.84, 0.87);
886 legend->SetFillColor(0);
887 legend->SetBorderSize(0);
888 legend->SetTextSize(0.05);
889
890 legend->AddEntry(hDeDx[0], "|#eta|<0.2", "P");
891 legend->AddEntry(hDeDx[1], "0.2<|#eta|<0.4", "P");
892 legend->AddEntry(hDeDx[2], "0.4<|#eta|<0.6", "P");
893 legend->AddEntry(hDeDx[3], "0.6<|#eta|<0.8", "P");
894
895 TCanvas* cTestEta = new TCanvas("cTestEta", "Comparing resolution for MIPs and electrons",
896 1200, 400);
897 cTestEta->Clear();
898 cTestEta->Divide(3, 1);
899
900 for(Int_t i = 0; i < nEta; i++) {
901
902 cTestEta->cd(1);
903 if(i==0)
904 hDeDx[i]->Draw("P");
905 else
906 hDeDx[i]->Draw("SAME P");
907
908 cTestEta->cd(2);
909 if(i==0)
910 hMeanVsNcl[i]->Draw("P");
911 else
912 hMeanVsNcl[i]->Draw("SAME P");
913
914 cTestEta->cd(3);
915 if(i==0) {
916 hSigmaVsNcl[i]->Draw("P");
917 legend->Draw();
918 } else
919 hSigmaVsNcl[i]->Draw("SAME P");
920 }
921
922 gROOT->ProcessLine(".x drawText.C(2)");
923 TString baseName(gSystem->BaseName(calibFileName));
924 baseName.ReplaceAll(".root", "");
925 cTestEta->SaveAs(Form("results/calibplots/test_resolution_vs_eta_%s.gif", baseName.Data()));
926 cTestEta->SaveAs(Form("results/calibplots/test_resolution_vs_eta_%s.pdf", baseName.Data()));
927}
928
929 // hMeanVsNcl = (TH1D*)helpArray->At(1);
930 // hMeanVsNcl->SetName("hMeanVsNcl");
931 // hMeanVsNcl->SetTitle("Mean vs Ncl; Ncl; Mean");
932 // hMeanVsNcl->SetMarkerStyle(29);
933 // hMeanVsNcl->SetMinimum(45);
934 // hMeanVsNcl->SetMaximum(55);
935
936 // if(isMC) {
937 // hDeDxVsNclElectrons->FitSlicesY(0, 0, -1, 0, "QNR", helpArray);
938 // hSigmaVsNclElectrons = (TH1D*)helpArray->At(2);
939 // hSigmaVsNclElectrons->SetName("hSigmaVsNcl");
940 // hSigmaVsNclElectrons->SetTitle("#sigma vs Ncl; Ncl; #sigma");
941 // }
942
943 // TCanvas* cNcl2 = new TCanvas("cNcl2", "sigma dE/dx vs Ncl", 600, 400);
944 // cNcl2->Clear();
945 // cNcl2->cd();
946
947 // hSigmaVsNcl->DrawCopy();
948 // hSigmaVsNcl->Fit(fSDeDxVsNcl, "0");
949 // hSigmaVsNcl->Fit(fSDeDxVsNclSqrt, "0");
950
951 // fSDeDxVsNcl->DrawCopy("SAME");
952 // fSDeDxVsNclSqrt->SetLineStyle(2);
953 // fSDeDxVsNclSqrt->DrawCopy("SAME");
954
955 // TLegend* legFit = new TLegend(0.64, 0.67, 0.84, 0.87);
956 // legFit->SetFillColor(0);
957 // legFit->SetBorderSize(0);
958 // legFit->SetTextSize(0.05);
959 // legFit->AddEntry(fSDeDxVsNcl, "Linear fit", "L");
960 // legFit->AddEntry(fSDeDxVsNclSqrt, "1/Sqrt fit", "L");
961 // legFit->Draw();
962
963 // cNcl2->SaveAs(Form("%s/sigma_vs_ncl.gif", baseName.Data()));
964
965 // if(isMC) {
966 // TCanvas* cNcl2Electrons = new TCanvas("cNcl2Electrons", "sigma dE/dx vs Ncl (comparison)", 600, 400);
967 // cNcl2Electrons->Clear();
968 // cNcl2Electrons->cd();
969
970 // hSigmaVsNcl->SetMarkerStyle(29);
971 // hSigmaVsNcl->DrawCopy();
972 // hSigmaVsNclElectrons->Scale(50.0 / hDeDxVsEtaCalibratedElectrons->GetMean(2));
973 // hSigmaVsNclElectrons->SetMarkerStyle(29);
974 // hSigmaVsNclElectrons->SetMarkerColor(6);
975 // hSigmaVsNclElectrons->DrawCopy("SAME");
976 // fSDeDxVsNcl->DrawCopy("SAME");
977 // fSDeDxVsNclSqrt->DrawCopy("SAME");
978
979 // TLegend* legSigma = new TLegend(0.64, 0.67, 0.84, 0.87);
980 // legSigma->SetFillColor(0);
981 // legSigma->SetBorderSize(0);
982 // legSigma->SetTextSize(0.05);
983 // legSigma->AddEntry(hSigmaVsNcl, "MIP pions", "P");
984 // legSigma->AddEntry(fSDeDxVsNcl, "Linear fit", "L");
985 // legSigma->AddEntry(fSDeDxVsNclSqrt, "1/Sqrt fit", "L");
986 // legSigma->AddEntry(hSigmaVsNclElectrons, "Scaled electrons", "P");
987 // legSigma->Draw();
988
989 // cNcl2Electrons->SaveAs(Form("%s/electrons_sigma_vs_ncl_comparison.gif", baseName.Data()));
990 // }
991
992 // TCanvas* cDeDx = new TCanvas("cDeDx", "dE/dx", 600, 400);
993 // cDeDx->Clear();
994 // hDeDx->Scale(1.0/hDeDx->GetEntries());
995 // hDeDx->SetMarkerStyle(29);
996 // hDeDx->DrawCopy();
997 // TCanvas* cNcl3 = new TCanvas("cNcl3", "#sigma vs Ncl", 600, 400);
998 // cNcl3->Clear();
999 // hSigmaVsNcl->DrawCopy();
1000 // TCanvas* cNcl3Mean = new TCanvas("cNcl3Mean", "Mean vs Ncl", 600, 400);
1001 // cNcl3Mean->Clear();
1002 // hMeanVsNcl->DrawCopy();
1003
1004 // TLegend* legDeDx = new TLegend(0.7, 0.67, 0.9, 0.87);
1005 // legDeDx->SetFillColor(0);
1006 // legDeDx->SetBorderSize(0);
1007 // legDeDx->SetTextSize(0.05);
1008 // legDeDx->AddEntry(hDeDx, "No #eta cut", "P");
1009
1010 // TLegend* legScan = new TLegend(0.44, 0.67, 0.64, 0.87);
1011 // legScan->SetFillColor(0);
1012 // legScan->SetBorderSize(0);
1013 // legScan->SetTextSize(0.05);
1014 // legScan->AddEntry(hSigmaVsNcl, "No #eta cut", "P");
1015
1016 // TH2D* hArray[4] = {hDeDxVsNcl0, hDeDxVsNcl1, hDeDxVsNcl2, hDeDxVsNcl3};
1017 // TH1D* hArrayDeDx[4] = {hDeDx0, hDeDx1, hDeDx3, hDeDx4};
1018
1019 // Int_t color[4] = {2, kGreen+1, 4, 6};
1020 // Double_t x[4];
1021 // Double_t x_err[4];
1022 // Double_t slope[4];
1023 // Double_t slope_err[4];
1024 // Double_t offset[4];
1025 // Double_t offset_err[4];
1026 // for(Int_t i = 0; i < 4; i++) {
1027
1028 // hArrayDeDx[i]->SetMarkerStyle(20);
1029 // hArrayDeDx[i]->SetMarkerColor(color[i]);
1030 // hArrayDeDx[i]->Scale(1.0/hArrayDeDx[i]->GetEntries());
1031 // cDeDx->cd();
1032 // hArrayDeDx[i]->DrawCopy("SAME");
1033
1034 // x[i] = hMeanEta->GetBinContent(i+1);
1035 // x_err[i] = hMeanEta->GetBinError(i+1);
1036
1037 // hArray[i]->FitSlicesY(0, 0, -1, 0, "QNR", helpArray);
1038 // hSigmaVsNclHelp = (TH1D*)helpArray->At(2);
1039 // hSigmaVsNclHelp->SetName(Form("hSigmaVsNcl%d", i));
1040 // hSigmaVsNclHelp->SetTitle("#sigma vs Ncl; Ncl; #sigma");
1041 // hSigmaVsNclHelp->SetMarkerStyle(20);
1042 // hSigmaVsNclHelp->SetMarkerColor(color[i]);
1043 // hSigmaVsNclHelp->Fit(fSDeDxVsNclHelp, "0");
1044 // slope[i] = fSDeDxVsNclHelp->GetParameter(0);
1045 // slope_err[i] = fSDeDxVsNclHelp->GetParError(0);
1046 // offset[i] = fSDeDxVsNclHelp->GetParameter(1);
1047 // offset_err[i] = fSDeDxVsNclHelp->GetParError(1);
1048 // cNcl3->cd();
1049 // hSigmaVsNclHelp->DrawCopy("SAME");
1050 // // fSDeDxVsNclHelp->DrawCopy("SAME");
1051
1052 // if(i==0) {
1053 // legDeDx->AddEntry(hArrayDeDx[i], "|#eta|<0.2", "P");
1054 // legScan->AddEntry(hSigmaVsNclHelp, "|#eta|<0.2", "P");
1055 // } else if(i==1) {
1056 // legDeDx->AddEntry(hArrayDeDx[i], "0.2<|#eta|<0.4", "P");
1057 // legScan->AddEntry(hSigmaVsNclHelp, "0.2<|#eta|<0.4", "P");
1058 // } else if(i==2) {
1059 // legDeDx->AddEntry(hArrayDeDx[i], "0.4<|#eta|<0.6", "P");
1060 // legScan->AddEntry(hSigmaVsNclHelp, "0.4<|#eta|<0.6", "P");
1061 // } else if(i==3) {
1062 // legDeDx->AddEntry(hArrayDeDx[i], "0.6<|#eta|<0.8", "P");
1063 // legScan->AddEntry(hSigmaVsNclHelp, "0.6<|#eta|<0.8", "P");
1064 // }
1065
1066 // cNcl3Mean->cd();
1067 // hMeanVsNclHelp = (TH1D*)helpArray->At(1);
1068 // hMeanVsNclHelp->SetName(Form("hMeanVsNcl%d", i));
1069 // hMeanVsNclHelp->SetTitle("Mean vs Ncl; Ncl; Mean");
1070 // hMeanVsNclHelp->SetMarkerStyle(20);
1071 // hMeanVsNclHelp->SetMarkerColor(color[i]);
1072 // hMeanVsNclHelp->DrawCopy("SAME");
1073 // }
1074
1075 // cDeDx->cd();
1076 // hDeDx->DrawCopy("SAME");
1077 // legDeDx->Draw();
1078 // cDeDx->SaveAs(Form("%s/dedx_resolution_vs_eta.gif", baseName.Data()));
1079
1080 // cNcl3->cd();
1081 // hSigmaVsNcl->DrawCopy("SAME");
1082 // legScan->Draw();
1083 // cNcl3->SaveAs(Form("%s/sigma_vs_ncl_vs_eta.gif", baseName.Data()));
1084
1085 // cNcl3Mean->cd();
1086 // hMeanVsNcl->DrawCopy("SAME");
1087 // legScan->Draw();
1088 // cNcl3Mean->SaveAs(Form("%s/mean_vs_ncl_vs_eta.gif", baseName.Data()));
1089
1090 // TGraphErrors* graphSlope = new TGraphErrors(4, x, slope, x_err, slope_err);
1091 // graphSlope->SetTitle("slope vs <#eta>; <#eta>; slope");
1092 // graphSlope->SetMarkerStyle(29);
1093 // TCanvas* cNcl4 = new TCanvas("cNcl4", "Slope vs <eta>", 600, 400);
1094 // cNcl4->Clear();
1095 // graphSlope->Draw("AP");
1096 // cNcl4->SaveAs(Form("%s/fit_slope_vs_eta.gif", baseName.Data()));
1097
1098 // TGraphErrors* graphOffset = new TGraphErrors(4, x, offset, x_err, offset_err);
1099 // graphOffset->SetMarkerStyle(29);
1100 // graphOffset->SetTitle("offset vs <#eta>; <#eta>; offset");
1101 // TCanvas* cNcl5 = new TCanvas("cNcl5", "Offset vs <eta>", 600, 400);
1102 // cNcl5->Clear();
1103 // graphOffset->Draw("AP");
1104 // cNcl5->SaveAs(Form("%s/fit_offset_vs_eta.gif", baseName.Data()));
1105
1106 // }
1107
1108//____________________________________________________________________________
1109void CompareYields(const Char_t* dataFileName1,
1110 const Char_t* dataFileName2,
1111 Double_t ptStart, Double_t ptStop,
1112 const Char_t* endName1=0,
1113 const Char_t* endName2=0,
1114 Int_t run = 0,
1115 Int_t filter = 1,
1116 Bool_t usePhiCut = kTRUE)
1117{
1118 gStyle->SetOptStat(0);
1119
1120
1121 TFile* dataFile1 = FindFileFresh(dataFileName1);
1122 if(!dataFile1)
1123 return;
1124 AliHighPtDeDxCalib* data1 = (AliHighPtDeDxCalib*)GetObject(dataFile1, filter, usePhiCut, run, "filter", endName1);
1125 data1->Print();
1126
1127 // gSystem->Exec("mv debug/* olddebug/");
1128
1129
1130 TH2D* hDeltaPiVsPt1 = data1->GetHistDeDxVsP(0);
1131
1132 AliHighPtDeDxCalib* data2 = data1;
1133 if(dataFileName2) {
1134 TFile* dataFile2 = FindFileFresh(dataFileName2);
1135 if(!dataFile2)
1136 return;
1137 data2 = (AliHighPtDeDxCalib*)GetObject(dataFile2, filter, usePhiCut, run, "filter", endName2);
1138 data2->Print();
1139 } else if (endName2) {
1140
1141 data2 = (AliHighPtDeDxCalib*)GetObject(dataFile1, filter, usePhiCut, run, "filter", endName2);
1142 }
1143
1144 TH2D* hDeltaPiVsPt2 = data2->GetHistDeDxVsP(0);
1145
1146 // Root is a bit stupid with finidng bins so we have to add and subtract a
1147 // little to be sure we get the right bin as we typically put edges as
1148 // limits
1149 const Int_t binStart = hDeltaPiVsPt1->GetXaxis()->FindBin(ptStart+0.01);
1150 ptStart = hDeltaPiVsPt1->GetXaxis()->GetBinLowEdge(binStart);
1151 const Int_t binStop = hDeltaPiVsPt1->GetXaxis()->FindBin(ptStop-0.01);
1152 ptStop = hDeltaPiVsPt1->GetXaxis()->GetBinUpEdge(binStop);
1153 // const Int_t nBins = binStop - binStart + 1;
1154
1155 cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
1156 << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
1157
1158
1159 //cross check
1160 TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
1161 cFits->Clear();
1162 cFits->Divide(7, 4);
1163
1164 TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
1165 legend->SetBorderSize(0);
1166 legend->SetFillColor(0);
1167
1168 TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
1169
1170 for(Int_t bin = binStart; bin <= binStop; bin++){
1171
1172 cout << "Making projection for bin: " << bin << endl;
1173
1174 const Int_t j = bin-binStart;
1175
1176 TH1D* hDeltaPiVsPtProj1 =(TH1D*)hDeltaPiVsPt1->ProjectionY(Form("hDeltaPiVsPtProj1%d", bin), bin, bin);
1177 // hDeltaPiVsPtProj1->GetXaxis()->SetRangeUser(-20, 20);
1178 hDeltaPiVsPtProj1->SetTitle(Form("%.1f<p<%.1f GeV/c",
1179 hDeltaPiVsPt1->GetXaxis()->GetBinLowEdge(bin),
1180 hDeltaPiVsPt1->GetXaxis()->GetBinUpEdge(bin)));
1181
1182 TH1D* hDeltaPiVsPtProj2 =(TH1D*)hDeltaPiVsPt2->ProjectionY(Form("hDeltaPiVsPtProj2%d", bin), bin, bin);
1183
1184 hDeltaPiVsPtProj1->SetLineColor(4);
1185 hDeltaPiVsPtProj2->SetLineColor(2);
1186 hDeltaPiVsPtProj1->SetNormFactor();
1187 hDeltaPiVsPtProj2->SetNormFactor();
1188
1189 cFits->cd();
1190 cFits->cd(j + 1);
1191 hDeltaPiVsPtProj1->DrawCopy();
1192 hDeltaPiVsPtProj2->DrawCopy("SAME");
1193
1194 cSingleFit->cd();
1195 cSingleFit->Clear();
1196 hDeltaPiVsPtProj1->DrawCopy();
1197 hDeltaPiVsPtProj2->DrawCopy("SAME");
1198
1199 CreateDir("results/debug");
1200 cSingleFit->SaveAs(Form("results/debug/ptspectrum_bin%d.gif", bin));
1201 }
1202}