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 }