Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / macros / calibrate_de_dx.C
1
2 #include <TFile.h>
3 #include <TH1.h>
4 #include <TProfile.h>
5 #include <TF1.h>
6 #include <TGraphErrors.h>
7 #include <TCanvas.h>
8 #include <TAxis.h>
9 #include <TLegend.h>
10 #include <TStyle.h>
11 #include <TGraphAsymmErrors.h>
12 #include <TList.h>
13 #include <TMath.h>
14 #include <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
31 using 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 //____________________________________________________________________________
95 void 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 //____________________________________________________________________________
692 void 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 //____________________________________________________________________________
749 void 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 //___________________________________________________________________________________________
832 void 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 //____________________________________________________________________________
1109 void 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 }