PWGLF/QATasks/post/PostProcessQAHighPtDeDx.C

   1 /*
   2
   3   Use AliRoot because of AliXRDPROOFtoolkit:
   4
   5
   6
   7   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/Base")
   8   gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros")
   9   gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid")
  10   gSystem->AddIncludePath("-I$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib")
  11   gROOT->SetMacroPath(".:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/macros:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/grid:$ALICE_ROOT/PWGLF/SPECTRA/IdentifiedHighPt/lib/")
  12   .L my_functions.C+
  13   .L my_tools.C+
  14   .L PostProcessQAHighPtDeDx.C+
  15
  16    PlotQA("FileRoot/AnalysisResults.root")
  17    MakeFitsExternalData("FileRoot/AnalysisResults.root", "HistosForBB")
  18
  19   MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",0)
  20   MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",1)
  21   MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",2)
  22   MakeFitsV0s("FileRoot/AnalysisResults.root", "HistosForBB/PrimaryElectrons.root", "HistosForBB",3)
  23
  24
  25   PlotParametrizations("HistosForBB")
  26
  27
  28
  29
  30   FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  0,  2, 0,1, 0, "HistosForBB/hres_0_5_02.root",27);
  31   FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  2,  4, 0,1, 0, "HistosForBB/hres_0_5_24.root",27);
  32   FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  4,  6, 0,1, 0, "HistosForBB/hres_0_5_46.root",27);
  33   FitDeDxVsP("FileRoot/AnalysisResults.root", 3.0, 10.0, 0, 6, 13, kTRUE,  6,  8, 0,1, 0, "HistosForBB/hres_0_5_68.root",27);
  34
  35
  36   MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/02_dataPbPb.root",2,50, 0,  "02_dataPbPb.root");
  37   MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/24_dataPbPb.root",2,50, 1,  "24_dataPbPb.root");
  38   MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/46_dataPbPb.root",2,50, 2,  "46_dataPbPb.root");
  39   MakeNSigmaPlot("FileRoot/AnalysisResults.root","fitparameters/MB/68_dataPbPb.root",2,50, 3,  "68_dataPbPb.root");
  40
  41
  42   PlotNSigma("nsigma_results")
  43
  44
  45
  46   //add particle fractions vs p
  47
  48
  49   ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,0, "results/eta02","fractions");
  50   ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,1, "results/eta24","fractions");
  51   ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,2, "results/eta46","fractions");
  52   ExtractUncPartFractvsP("FileRoot/AnalysisResults.root", 2, 50,0,3, "results/eta68","fractions");
  53
  54
  55 */
  56 #include <TLegend.h>
  57 #include <TFile.h>
  58 #include <TList.h>
  59 #include <TTree.h>
  60 #include <TH2.h>
  61 #include <TF1.h>
  62 #include <TF2.h>
  63 #include <TH1.h>
  64 #include <TCanvas.h>
  65 #include <TProfile.h>
  66 #include <TClonesArray.h>
  67 #include <TObjString.h>
  68 #include <TString.h>
  69 #include <TROOT.h>
  70 #include <TStyle.h>
  71 #include <TLatex.h>
  72 #include <TGraphErrors.h>
  73 #include <TSpectrum.h>
  74 #include <TFitResultPtr.h>
  75 #include <TFitResult.h>
  76
  77 #include "my_tools.C"
  78 #include "my_functions.C"
  79
  80 #include <iostream>
  81 #include <fstream>
  82 #include <string>
  83
  84 using namespace std;
  85
  86
  87
  88
  89
  90 //const Char_t* ending[4] = {"02", "24", "46", "68"};
  91
  92 const Char_t * legEtaCut[4]=
  93   {"|#eta_{lab}|<0.2",  "0.2<|#eta_{lab}|<0.4",   "0.4<|#eta_{lab}|<0.6",   "0.6<|#eta_{lab}|<0.8"};
  94
  95
  96 const Char_t* endingCent[1] =
  97   {"MB"};
  98 const Char_t* CentLatex[1] =
  99   {"MB"};
 100
 101   const Char_t *Pid[4]={"Pion", "Kaon", "Proton", "Electron"};
 102 //const Char_t *Pid[3]={"Pion","Kaon","Proton"};
 103 const Char_t *PidLatex[3]={"#pi^{+} + #pi^{-}","K^{+} + K^{-}","p + #bar{p}"};
 104
 105
 106 const Color_t PidColor[3]={2,kGreen,4};
 107
 108 TF1* piFunc = 0;
 109 TF1* kFunc  = 0;
 110 TF1* pFunc = 0;
 111 TF1* eFunc = 0;
 112 TF1* sigmaFunc = 0;
 113
 114 const Int_t nPtBins = 68;
 115
 116 Double_t xBins[nPtBins+1] = {0. ,  0.05, 0.1,  0.15, 0.2,  0.25, 0.3,  0.35, 0.4,  0.45,
 117                              0.5,  0.55, 0.6,  0.65, 0.7,  0.75, 0.8,  0.85, 0.9,  0.95,
 118                              1.0,  1.1 , 1.2,  1.3 , 1.4,  1.5 , 1.6,  1.7 , 1.8,  1.9 ,
 119                              2.0,  2.2 , 2.4,  2.6 , 2.8,  3.0 , 3.2,  3.4 , 3.6,  3.8 ,
 120                              4.0,  4.5 , 5.0,  5.5 , 6.0,  6.5 , 7.0,  8.0 , 9.0,  10.0,
 121                              11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
 122                              26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
 123
 124
 125
 126
 127 const Int_t nPtBinsV0s = 25;
 128 Double_t ptBinsV0s[nPtBinsV0s+1] = { 0.0 , 0.1 , 0.2 , 0.3 , 0.4 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 , 1.0 ,
 129                                      1.2 , 1.4 , 1.6 , 1.8 , 2.0 , 2.5 , 3.0 , 3.5 , 4.0 , 5.0 , 7.0 ,
 130                                      9.0 , 12.0, 15.0, 20.0 };
 131
 132
 133 TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist);
 134
 135
 136 void PlotNSigma(const Char_t *path){
 137
 138   gStyle->SetCanvasColor(10);
 139   gStyle->SetFrameFillColor(10);
 140   gStyle->SetOptStat(0);
 141   gStyle->SetOptTitle(0);
 142   gStyle->SetOptFit(0);
 143
 144
 145
 146   TFile *fin[2]={0,0};
 147
 148   fin[0] = FindFileFresh(Form("%s/02_dataPbPb.root",path));
 149   if(!fin[0])
 150     return;
 151
 152   fin[1] = FindFileFresh(Form("%s/68_dataPbPb.root",path));
 153   if(!fin[1])
 154     return;
 155
 156   TH1D *hNPiP[2]={0,0};
 157   TH1D *hNPiK[2]={0,0};
 158   TH1D *hNPK[2]={0,0};
 159   Int_t colorestmp[2]={2,4};
 160   for(Int_t i=0; i<2; ++i){
 161
 162     hNPiP[i]=(TH1D *)fin[i]->Get("hPionProtonSeparation");
 163     hNPiP[i]->SetLineColor(colorestmp[i]);
 164     hNPiP[i]->SetLineWidth(2);
 165     hNPiK[i]=(TH1D *)fin[i]->Get("hPionKaonSeparation");
 166     hNPiK[i]->SetLineColor(colorestmp[i]);
 167     hNPiK[i]->SetLineWidth(2);
 168     hNPK[i]=(TH1D *)fin[i]->Get("hKaonProtonSeparation");
 169     hNPK[i]->SetLineColor(colorestmp[i]);
 170     hNPK[i]->SetLineWidth(2);
 171
 172   }
 173
 174
 175   TH1D *hframe[3]={0,0,0};
 176   for(Int_t i=0; i<3; ++i){
 177
 178     hframe[i] = new TH1D(Form("hframe%d",i),";#it{p} (GeV/#it{c});",10,3,15);
 179     hframe[i]->GetYaxis()->SetRangeUser(-0.01,7.01);
 180     hframe[i]->GetXaxis()->SetTitleSize(0.05);
 181     hframe[i]->GetYaxis()->SetTitleSize(0.05);
 182   }
 183
 184
 185   TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,1050,400);
 186   TPad *pad[3]={0,0,0};
 187   for(Int_t i=0; i<3; ++i){
 188
 189     pad[i]=new TPad(Form("pad%d",i),"pad1",0.01+0.33*i,0.01,0.33+0.33*i,0.99,0);
 190     pad[i]->SetLeftMargin(0.12);
 191     pad[i]->SetRightMargin(0.01);
 192     pad[i]->SetTopMargin(0.01);
 193     pad[i]->SetBottomMargin(0.1);
 194     pad[i]->SetBorderSize(0);
 195     pad[i]->SetBorderMode(0);
 196
 197   }
 198
 199
 200
 201
 202   for(Int_t i =0 ; i<3; ++i){
 203
 204
 205
 206
 207     vC1->cd();
 208     pad[i]->Draw();
 209
 210
 211     pad[i]->cd();
 212
 213
 214     pad[i]->SetTickx(1);
 215     pad[i]->SetTicky(1);
 216
 217     if(i==0)
 218       hframe[i]->GetYaxis()->SetTitle("Separation #pi/p, N_{#sigma}");
 219     if(i==1)
 220       hframe[i]->GetYaxis()->SetTitle("Separation #pi/K, N_{#sigma}");
 221     if(i==2)
 222       hframe[i]->GetYaxis()->SetTitle("Separation p/K, N_{#sigma}");
 223
 224
 225     hframe[i]->Draw();
 226
 227     if(i==0){
 228       for(Int_t j =0 ; j<2; ++j){
 229
 230         hNPiP[j]->Draw("samel");
 231
 232       }
 233       TF1 *fpip1=new TF1("fpip1","5.0+pol0",5,15);
 234       fpip1->SetLineColor(1);
 235       fpip1->SetLineWidth(2);
 236       fpip1->SetLineStyle(2);
 237       fpip1->Draw("same");
 238
 239       TF1 *fpip2=new TF1("fpip2","3.5+pol0",5,15);
 240       fpip2->SetLineColor(1);
 241       fpip2->SetLineWidth(2);
 242       fpip2->SetLineStyle(4);
 243       fpip2->Draw("same");
 244
 245       TLatex* latex0 = new TLatex();
 246       latex0->SetNDC();
 247       latex0->SetTextAlign(22);
 248       latex0->SetTextFont(42);
 249       latex0->SetTextSize(0.05);
 250       latex0->DrawLatex(0.5,0.77,"pp @ 2.76 TeV, 0.6<|#eta|<0.8 (final)");
 251
 252       latex0->DrawLatex(0.5,0.57,"Pb-Pb, 0-5%, 0.6<|#eta|<0.8 (final)");
 253
 254     }
 255
 256     if(i==1){
 257       for(Int_t j =0 ; j<2; ++j){
 258         hNPiK[j]->Draw("samel");
 259
 260       }
 261
 262
 263       TF1 *fpik1=new TF1("fpik1","3.5+pol0",5,15);
 264       fpik1->SetLineColor(1);
 265       fpik1->SetLineWidth(2);
 266       fpik1->SetLineStyle(2);
 267       fpik1->Draw("same");
 268
 269       TF1 *fpik2=new TF1("fpik2","2.2+pol0",5,15);
 270       fpik2->SetLineColor(1);
 271       fpik2->SetLineWidth(2);
 272       fpik2->SetLineStyle(4);
 273       fpik2->Draw("same");
 274
 275       TLegend* legend = new TLegend(0.51, 0.65, 0.85, 0.85);
 276       legend->SetBorderSize(0);
 277       legend->SetFillColor(0);
 278       legend->SetTextFont(42);
 279       legend->SetTextSize(0.05);
 280       legend->SetLineColor(kWhite);
 281       legend->SetLineStyle(3);
 282       legend->SetShadowColor(kWhite);
 283       legend->SetFillColor(kWhite);
 284       legend->SetFillStyle(0);
 285
 286       legend->AddEntry(hNPiK[0], "|#eta|<0.2", "L");
 287       legend->AddEntry(hNPiK[1], "0.6<|#eta|<0.8", "L");
 288       legend->Draw();
 289
 290
 291     }
 292
 293     if(i==2){
 294       for(Int_t j =0 ; j<2; ++j){
 295         hNPK[j]->Draw("samel");
 296
 297       }
 298
 299       TF1 *fpk1=new TF1("fpk1","1.8+pol0",5,15);
 300       fpk1->SetLineColor(1);
 301       fpk1->SetLineWidth(2);
 302       fpk1->SetLineStyle(2);
 303       fpk1->Draw("same");
 304
 305       TF1 *fpk2=new TF1("fpk2","1.2+pol0",5,15);
 306       fpk2->SetLineColor(1);
 307       fpk2->SetLineWidth(2);
 308       fpk2->SetLineStyle(4);
 309       fpk2->Draw("same");
 310
 311     }
 312
 313   }
 314
 315
 316
 317   vC1->SaveAs("NSigma.png");
 318   vC1->SaveAs("NSigma.eps");
 319
 320
 321
 322
 323
 324
 325
 326 }
 327
 328
 329
 330 //____________________________________________________________________________
 331 void MakeNSigmaPlot(const Char_t* inFile, const Char_t* fitFileName,
 332                     Double_t ptStart, Double_t ptStop, Int_t i_eta,
 333                     const Char_t* outFileName=0)
 334 {
 335
 336   const Char_t* ending[4] = {"02", "24", "46", "68"};
 337   gStyle->SetOptStat(0);
 338
 339   if(fitFileName) {
 340
 341     TFile* fitFile = FindFileFresh(fitFileName);
 342     if(!fitFile)
 343       return;
 344     DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
 345     fitPar->Print();
 346
 347     fixMIP      = fitPar->MIP;
 348     fixPlateau  = fitPar->plateau;
 349
 350     Double_t dedxPar[6]  = {0, 0, 0, 0, 0, 0};
 351     Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
 352
 353     dedxPar[0] = fitPar->optionDeDx;
 354     for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
 355       dedxPar[i+1] = fitPar->parDeDx[i];
 356     }
 357
 358     sigmaPar[0] = fitPar->optionSigma;
 359     for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
 360       sigmaPar[i+1] = fitPar->parSigma[i];
 361     }
 362
 363     piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 364     piFunc->SetParameters(dedxPar);
 365
 366
 367     kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 368     kFunc->SetParameters(dedxPar);
 369     kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
 370
 371     pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 372     pFunc->SetParameters(dedxPar);
 373     pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
 374
 375     eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 376     eFunc->SetParameters(dedxPar);
 377     eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
 378
 379     sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1);
 380     sigmaFunc->SetParameters(sigmaPar);
 381   }
 382
 383   kFunc->SetLineColor(kGreen+2);
 384   pFunc->SetLineColor(4);
 385   eFunc->SetLineColor(kMagenta);
 386
 387
 388   /*
 389   TString baseName(gSystem->BaseName(calibFileName));
 390   baseName.ReplaceAll(".root", "");
 391
 392   TFile* calibFile = FindFileFresh(calibFileName);
 393   if(!calibFile)
 394     return;
 395   AliHighPtDeDxCalib* calib = (AliHighPtDeDxCalib*)GetObject(calibFile, filter, usePhiCut, run, etaAbs, etaLow, etaHigh);
 396   if(calib) {
 397     fixMIP      = calib->GetHistDeDx(kTRUE, 0)->GetMean();
 398     fixPlateau  = calib->GetHistDeDx(kFALSE, 0)->GetMean();
 399     cout << "Plateau: " << fixPlateau << endl;
 400   } else {
 401     calib = (AliHighPtDeDxCalib*)calibFile->Get("v0phicut");
 402     fixMIP = 50.0;
 403     fixPlateau  = 83.461;
 404   }
 405   calib->Print();
 406
 407   hDeDxVsP = calib->GetHistDeDxVsP(0);
 408   */
 409
 410   TFile* dedxFile = FindFileFresh(inFile);
 411   if(!dedxFile)
 412     return;
 413
 414   TList * list = (TList *)dedxFile->Get("outputdedx");
 415   TH2D *hDeDxVsPPlus = 0;
 416   hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
 417   hDeDxVsPPlus->Sumw2();
 418
 419   TH2D *hDeDxVsPMinus = 0;
 420   hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
 421   hDeDxVsPMinus->Sumw2();
 422
 423   gSystem->Exec("rm *.root");
 424   TFile *fplus = new TFile("hist2Dplus.root","recreate");
 425   fplus->cd();
 426   hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
 427   hDeDxVsPPlus->Write();
 428   fplus->Close();
 429
 430   TFile *fminus = new TFile("hist2Dminus.root","recreate");
 431   fminus->cd();
 432   hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
 433   hDeDxVsPMinus->Write();
 434   fminus->Close();
 435
 436   gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
 437
 438   TFile *ffinal = TFile::Open("hist2Ddedx.root");
 439   hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
 440
 441
 442
 443
 444    //in this case, sigma is not the relative resolution
 445
 446   // Root is a bit stupid with finidng bins so we have to add and subtract a
 447   // little to be sure we get the right bin as we typically put edges as
 448   // limits
 449   const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
 450   ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
 451   const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
 452   ptStop  = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
 453   //  const Int_t nBins    = binStop - binStart + 1;
 454
 455   cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
 456        << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
 457
 458
 459   TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
 460   hPionRatio->Reset();
 461
 462   TH1D* hPionKaonSeparation = (TH1D*)hPionRatio->Clone("hPionKaonSeparation");
 463   hPionKaonSeparation->SetTitle(Form("#pi/K separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
 464   TH1D* hPionProtonSeparation = (TH1D*)hPionRatio->Clone("hPionProtonSeparation");
 465   hPionProtonSeparation->SetTitle(Form("#pi/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
 466   TH1D* hKaonProtonSeparation = (TH1D*)hPionRatio->Clone("hKaonProtonSeparation");
 467   hKaonProtonSeparation->SetTitle(Form("K/p separation, %s; #it{p} (GeV/#it{c}); N_{#sigma}",legEtaCut[i_eta]));
 468
 469
 470
 471
 472   for(Int_t bin = binStart; bin <= binStop; bin++){
 473
 474     Double_t pt=hPionRatio->GetBinCenter(bin);
 475
 476     Double_t dEdx_pi = piFunc->Eval(pt);
 477     Double_t Sigma_pi = sigmaFunc->Eval(piFunc->Eval(pt));
 478
 479     Double_t dEdx_k = kFunc->Eval(pt);
 480     Double_t Sigma_k = sigmaFunc->Eval(kFunc->Eval(pt));
 481
 482     Double_t dEdx_p = pFunc->Eval(pt);
 483     Double_t Sigma_p = sigmaFunc->Eval(pFunc->Eval(pt));
 484
 485     Double_t N1 = (dEdx_pi - dEdx_k)/( (Sigma_pi + Sigma_k)/2.0 );
 486     Double_t N2 = (dEdx_pi - dEdx_p)/( (Sigma_pi + Sigma_p)/2.0 );
 487     Double_t N3 = (dEdx_k - dEdx_p)/( (Sigma_k + Sigma_p)/2.0 );
 488
 489     hPionKaonSeparation->SetBinContent(bin,N1);
 490     hPionProtonSeparation->SetBinContent(bin,N2);
 491     hKaonProtonSeparation->SetBinContent(bin,N3);
 492
 493   }
 494
 495
 496
 497   if(outFileName) {
 498
 499     CreateDir("nsigma_results");
 500     TFile* fileOut = new TFile(Form("nsigma_results/%s", outFileName), "RECREATE");
 501
 502     hPionKaonSeparation->Write();
 503     hPionProtonSeparation->Write();
 504     hKaonProtonSeparation->Write();
 505
 506
 507
 508     fileOut->Close();
 509   }
 510
 511
 512
 513 }
 514
 515
 516
 517 //________________________________________________________________________
 518 void PlotParametrizations(const Char_t *path){
 519
 520
 521   gStyle->SetCanvasColor(10);
 522   gStyle->SetFrameFillColor(10);
 523   gStyle->SetOptStat(0);
 524   gStyle->SetOptTitle(0);
 525   gStyle->SetOptFit(0);
 526
 527   TFile *fin[2]={0,0};
 528
 529
 530   fin[0] = FindFileFresh(Form("%s/hres_0_5_02.root",path));
 531   if(!fin[0])
 532     return;
 533
 534   fin[1] = FindFileFresh(Form("%s/hres_0_5_68.root",path));
 535   if(!fin[1])
 536     return;
 537
 538
 539   TH1D *hframeRes = new TH1D("hframeRes","; #LT dE/dx #GT; #sigma/#LT dE/dx #GT",50,40,79);
 540   hframeRes->GetYaxis()->SetRangeUser(0,0.105);
 541   hframeRes->GetYaxis()->SetTitleSize(0.05);
 542   hframeRes->GetXaxis()->SetTitleSize(0.05);
 543
 544   TH1D *hframeBB = new TH1D("hframeBB","; #beta#gamma;#LT dE/dx #GT",50,2,59);
 545   hframeBB->GetYaxis()->SetRangeUser(40,81);
 546   hframeBB->GetYaxis()->SetTitleSize(0.05);
 547   hframeBB->GetXaxis()->SetTitleSize(0.05);
 548
 549   TGraphErrors *gRes[2] = {0,0};
 550
 551   gRes[0] = (TGraphErrors *)fin[0]->Get("res_allpions");
 552   gRes[1] = (TGraphErrors *)fin[1]->Get("res_allpions");
 553
 554
 555   TGraphErrors *gBB[2] = {0,0};
 556
 557   gBB[0] = (TGraphErrors *)fin[0]->Get("gBB");
 558   gBB[1] = (TGraphErrors *)fin[1]->Get("gBB");
 559
 560   Int_t colorestmp[2]={2,4};
 561
 562
 563
 564   TCanvas* vC1 = new TCanvas("vC1","vC1",200,10,700,400);
 565   TPad * pad1=new TPad("pad1","pad1",0.01,0.01,0.5,0.99,0);
 566   TPad * pad2=new TPad("pad2","pad2",0.5,0.01,0.99,0.99,0);
 567
 568   pad1->SetLeftMargin(0.12);
 569   pad1->SetRightMargin(0.01);
 570   pad1->SetTopMargin(0.01);
 571   pad1->SetBottomMargin(0.1);
 572   pad1->SetBorderSize(0);
 573   pad1->SetBorderMode(0);
 574
 575   pad2->SetLeftMargin(0.12);
 576   pad2->SetRightMargin(0.01);
 577   pad2->SetBottomMargin(0.1);
 578   pad2->SetTopMargin(0);
 579   pad2->SetBorderSize(0);
 580   vC1->cd();
 581   pad1->Draw();
 582
 583
 584   pad1->cd();
 585
 586
 587   pad1->SetTickx(1);
 588   pad1->SetTicky(1);
 589
 590
 591
 592
 593
 594
 595   hframeRes->Draw();
 596
 597   for(Int_t i =0 ; i<2; ++i){
 598     gRes[i]->SetMarkerStyle(20);
 599     gRes[i]->SetLineWidth(2);
 600     gRes[i]->SetLineColor(colorestmp[i]);
 601     gRes[i]->SetMarkerColor(colorestmp[i]);
 602     gRes[i]->Draw("samep");
 603   }
 604
 605
 606   TF1 *f0005 = new TF1("f0005","pol2",50,80);
 607   //  f0005->SetParameter(0,1.20000000000000000e+01);
 608   f0005->SetParameter(0,6.08803411659195118e-02);
 609   f0005->SetParameter(1,2.34760597152924448e-04);
 610   f0005->SetParameter(2,-2.81841683363273653e-06);
 611
 612   f0005->SetLineColor(1);
 613   f0005->SetLineWidth(2);
 614   f0005->SetLineStyle(2);
 615   f0005->Draw("samel");
 616
 617
 618   TF1 *fpp = new TF1("fpp","pol2",50,80);
 619   fpp->SetParameter(0,1.06773399595441187e-01);
 620   fpp->SetParameter(1,-1.55087701882609280e-03);
 621   fpp->SetParameter(2,1.01579672586848698e-05);
 622
 623   fpp->SetLineColor(1);
 624   fpp->SetLineWidth(2);
 625   fpp->SetLineStyle(3);
 626   fpp->Draw("samel");
 627
 628
 629   TLegend* legendRes = new TLegend(0.21, 0.15, 0.65, 0.35);
 630   legendRes->SetBorderSize(0);
 631   legendRes->SetFillColor(0);
 632   legendRes->SetTextFont(42);
 633   legendRes->SetTextSize(0.05);
 634   legendRes->SetLineColor(kWhite);
 635   legendRes->SetLineStyle(3);
 636   legendRes->SetShadowColor(kWhite);
 637   legendRes->SetFillColor(kWhite);
 638   legendRes->SetFillStyle(0);
 639   legendRes->SetHeader("Final, 0.6<|#eta|<0.8");
 640   legendRes->AddEntry(f0005, "Pb-Pb 0-5%", "L");
 641   legendRes->AddEntry(fpp, "pp @ 2.76 TeV", "L");
 642   legendRes->Draw();
 643
 644   vC1->cd();
 645   pad2->Draw();
 646
 647
 648   pad2->cd();
 649
 650
 651   pad2->SetTickx(1);
 652   pad2->SetTicky(1);
 653
 654
 655
 656
 657
 658
 659   hframeBB->Draw();
 660
 661   for(Int_t i =0 ; i<2; ++i){
 662     gBB[i]->SetMarkerStyle(20);
 663     gBB[i]->SetLineWidth(2);
 664     gBB[i]->SetLineColor(colorestmp[i]);
 665     gBB[i]->SetMarkerColor(colorestmp[i]);
 666     gBB[i]->Draw("samep");
 667   }
 668
 669   TLegend* legend = new TLegend(0.51, 0.8, 0.85, 0.95);
 670   legend->SetBorderSize(0);
 671   legend->SetFillColor(0);
 672   legend->SetTextFont(42);
 673   legend->SetTextSize(0.05);
 674   legend->SetLineColor(kWhite);
 675   legend->SetLineStyle(3);
 676   legend->SetShadowColor(kWhite);
 677   legend->SetFillColor(kWhite);
 678   legend->SetFillStyle(0);
 679
 680   legend->AddEntry(gBB[0], "|#eta|<0.2", "P");
 681   legend->AddEntry(gBB[1], "0.6<|#eta|<0.8", "P");
 682   legend->Draw();
 683
 684
 685
 686
 687   vC1->SaveAs("Parametrizations.png");
 688   vC1->SaveAs("Parametrizations.eps");
 689
 690
 691
 692 }
 693 //____________________________________________________________________________
 694 void ExtractUncPartFractvsP(const Char_t * inFile,  Double_t ptStart, Double_t ptStop,
 695                             Int_t i_cent=1,
 696                             Int_t i_eta=1,
 697                             const Char_t* dirName="debugfits", const Char_t* outFileName=0)
 698
 699
 700 {
 701
 702
 703   const Char_t* ending[4] = {"02", "24", "46", "68"};
 704   const Double_t etaHigh[4]={2,4,6,8};
 705   gStyle->SetOptStat(0);
 706
 707
 708
 709   TFile* fitFile = FindFileFresh(Form("fitparameters/MB/%s_dataPbPb.root",ending[i_eta]));
 710   if(!fitFile)
 711     return;
 712
 713
 714   DeDxFitInfo* fitPar = (DeDxFitInfo*)fitFile->Get("fitInfo");
 715   fitPar->Print();
 716
 717   fixMIP      = fitPar->MIP;
 718
 719
 720
 721   //return;
 722
 723   Double_t dedxPar[6]  = {0, 0, 0, 0, 0, 0};
 724   Double_t sigmaPar[8] = {0, 0, 0, 0, 0, 0, 0, 0};
 725
 726   dedxPar[0] = fitPar->optionDeDx;
 727   for(Int_t i = 0; i < fitPar->nDeDxPar; i++) {
 728     dedxPar[i+1] = fitPar->parDeDx[i];
 729     cout<<"idedx_par"<<i+1<<"="<<dedxPar[i+1]<<endl;
 730
 731   }
 732   fixPlateau  = dedxPar[5];
 733   cout<<"fixPlateau="<<fixPlateau<<"  fixMIP="<<fixMIP<<endl;
 734   //return;
 735
 736   sigmaPar[0] = fitPar->optionSigma;
 737   for(Int_t i = 0; i < fitPar->nSigmaPar; i++) {
 738     sigmaPar[i+1] = fitPar->parSigma[i];
 739   }
 740
 741   piFunc = new TF1("piFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 742   piFunc->SetParameters(dedxPar);
 743
 744
 745   kFunc = new TF1("kFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 746   kFunc->SetParameters(dedxPar);
 747   kFunc->SetParameter(0, kFunc->GetParameter(0)+10);
 748
 749   pFunc = new TF1("pFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 750   pFunc->SetParameters(dedxPar);
 751   pFunc->SetParameter(0, pFunc->GetParameter(0)+20);
 752
 753   eFunc = new TF1("eFunc", FitFunc, 0, 100, fitPar->nDeDxPar+1);
 754   eFunc->SetParameters(dedxPar);
 755   eFunc->SetParameter(0, eFunc->GetParameter(0)+30);
 756
 757   sigmaFunc = new TF1("sigmaFunc", SigmaFunc, 0, 100, fitPar->nSigmaPar+1);
 758   sigmaFunc->SetParameters(sigmaPar);
 759
 760
 761   kFunc->SetLineColor(kGreen+2);
 762   pFunc->SetLineColor(4);
 763   eFunc->SetLineColor(kMagenta);
 764
 765   /*
 766     TFile* calibFile = FindFileFresh(inFile);
 767     if(!calibFile)
 768     return;
 769     hDeDxVsP = (TH2D *)calibFile->Get(Form("hDeDxVsP%s_%s", ending[i_eta], endingCent[i_cent]));
 770     hDeDxVsP->Sumw2();
 771   */
 772
 773   TFile* dedxFile = FindFileFresh(inFile);
 774   if(!dedxFile)
 775     return;
 776
 777   TList * list = (TList *)dedxFile->Get("outputdedx");
 778   TH2D *hDeDxVsPPlus = 0;
 779   hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
 780   hDeDxVsPPlus->Sumw2();
 781
 782   TH2D *hDeDxVsPMinus = 0;
 783   hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%s",ending[i_eta]));
 784   hDeDxVsPMinus->Sumw2();
 785
 786   gSystem->Exec("rm *.root");
 787   TFile *fplus = new TFile("hist2Dplus.root","recreate");
 788   fplus->cd();
 789   hDeDxVsPPlus->SetName(Form("histAllCh%s", ending[i_eta]));
 790   hDeDxVsPPlus->Write();
 791   fplus->Close();
 792
 793   TFile *fminus = new TFile("hist2Dminus.root","recreate");
 794   fminus->cd();
 795   hDeDxVsPMinus->SetName(Form("histAllCh%s", ending[i_eta]));
 796   hDeDxVsPMinus->Write();
 797   fminus->Close();
 798
 799   gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
 800
 801   TFile *ffinal = TFile::Open("hist2Ddedx.root");
 802   hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%s", ending[i_eta]));
 803
 804
 805
 806
 807
 808
 809
 810   CreateDir(Form("fit_results/%s_%s",dirName,endingCent[i_cent]));
 811
 812   TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}", 500, 400);
 813   cDeDxVsPLogX->Clear();
 814   cDeDxVsPLogX->cd();
 815   cDeDxVsPLogX->SetLogz();
 816
 817   TH2D *hDeDxVsP2=(TH2D *)hDeDxVsP->Clone("h2DP");
 818   hDeDxVsP2->SetTitle("; #it{p} (GeV/#it{c}) ;d#it{E}/d#it{x}");
 819   hDeDxVsP2->GetXaxis()->SetRangeUser(2,19.5);
 820   hDeDxVsP2->Draw("COLZ");
 821
 822   piFunc->SetLineColor(1);
 823   piFunc->Draw("samel");
 824
 825   eFunc->SetLineColor(1);
 826   eFunc->Draw("samel");
 827
 828   kFunc->SetLineColor(1);
 829   kFunc->Draw("samel");
 830
 831   pFunc->SetLineColor(1);
 832   pFunc->Draw("samel");
 833
 834   TLatex* latex0 = new TLatex();
 835   latex0->SetNDC();
 836   latex0->SetTextAlign(22);
 837   latex0->SetTextSize(0.05);
 838   latex0->DrawLatex(0.7,0.94+0.035,Form("p-Pb, %s",legEtaCut[i_eta]));
 839   latex0->DrawLatex(0.71,0.89+0.035,Form("#sqrt{s_{NN}}=5.02 TeV, %s",CentLatex[i_cent]));
 840
 841
 842
 843   cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.png",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
 844   cDeDxVsPLogX->SaveAs(Form("fit_results/%s_%s/%s.eps",dirName,endingCent[i_cent],cDeDxVsPLogX->GetName()));
 845   //in this case, sigma is not the relative resolution
 846
 847   // Root is a bit stupid with finidng bins so we have to add and subtract a
 848   // little to be sure we get the right bin as we typically put edges as
 849   // limits
 850   const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(ptStart+0.01);
 851   ptStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
 852   const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(ptStop-0.01);
 853   ptStop  = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
 854   //  const Int_t nBins    = binStop - binStart + 1;
 855
 856   cout << "Doing fits from pTlow = " << ptStart << " (bin: " << binStart
 857        << ") to pThigh = " << ptStop << " (bin: " << binStop << ")" << endl;
 858
 859
 860   //cross check
 861   TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
 862   cFits->Clear();
 863   cFits->Divide(7, 4);
 864
 865
 866   TF1* pion = new TF1("pion", "gausn", 40, 100);
 867
 868   TF1* kaon = new TF1("kaon", "gausn", 40, 100);
 869
 870   TF1* proton = new TF1("proton", "gausn", 40, 100);
 871   //proton->SetLineWidth(2);
 872   //proton->SetLineColor(kBlue);
 873   TF1* electron = new TF1("electron", "gausn", 40, 100);
 874   //electron->SetLineWidth(2);
 875   //electron->SetLineColor(kMagenta);
 876   //  TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)+gausn(9)", -30, 30);
 877   TF1* total = 0;
 878   total = new TF1("total", "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
 879
 880   total->SetLineColor(kBlack);
 881   total->SetLineWidth(2);
 882   total->SetLineStyle(2);
 883
 884   TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
 885   legend->SetBorderSize(0);
 886   legend->SetFillColor(0);
 887   legend->AddEntry(total, "4-Gauss fit", "L");
 888   legend->AddEntry(pion, "#pi", "L");
 889   legend->AddEntry(kaon, "K", "L");
 890   legend->AddEntry(proton, "p", "L");
 891   legend->AddEntry(electron, "e", "L");
 892
 893   TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
 894   cSingleFit->SetLeftMargin(0.11);
 895   cSingleFit->SetRightMargin(0.03);
 896   cSingleFit->SetTopMargin(0.01);
 897   cSingleFit->SetBottomMargin(0.1);
 898   cSingleFit->SetBorderSize(0);
 899   cSingleFit->SetBorderMode(0);
 900
 901
 902   TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
 903   hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
 904   hPionRatio->Reset();
 905   TH1D* hKaonRatio   = (TH1D*)hPionRatio->Clone("hKaonRatio");
 906   TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
 907   TH1D* hElectronRatio = (TH1D*)hPionRatio->Clone("hElectronRatio");
 908
 909
 910
 911   TH1D* hPionYield =(TH1D*)hDeDxVsP->ProjectionX("hPionYield", 1, 1);
 912   hPionYield->SetTitle("particle fractions; p [GeV/c]; particle fraction");
 913   hPionYield->Reset();
 914   TH1D* hKaonYield   = (TH1D*)hPionYield->Clone("hKaonYield");
 915   TH1D* hProtonYield = (TH1D*)hPionYield->Clone("hProtonYield");
 916   TH1D* hElectronYield = (TH1D*)hPionYield->Clone("hElectronYield");
 917
 918
 919   TF1 *fmip=new TF1("fmip","pol0",0,1);
 920   fmip->SetParameter(0,fixMIP);
 921
 922
 923
 924   TF1* fElectronFraction = 0;
 925   if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
 926     fElectronFraction = new TF1("fElectronFraction", "[0]+(x<9.0)*[1]*(x-9.0)", 0, ptStop);
 927   if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
 928     fElectronFraction = new TF1("fElectronFraction", "[0]+(x<8.0)*[1]*(x-8.0)", 0, ptStop);
 929   if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
 930     fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.5)*[1]*(x-7.5)", 0, ptStop);
 931   if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
 932     fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
 933
 934   //TF1* fElectronFraction = new TF1("fElectronFraction", "[0]+(x<7.0)*[1]*(x-7.0)", 0, ptStop);
 935   fElectronFraction->SetParameters(0.01, 0.0);
 936
 937
 938   TH1D *dEdxpoints[binStop-binStart];
 939   TF1 *fPion[binStop-binStart];
 940   TF1 *fKaon[binStop-binStart];
 941   TF1 *fProton[binStop-binStart];
 942   TF1 *fElectron[binStop-binStart];
 943   TF1 *fTotal[binStop-binStart];
 944
 945
 946   for(Int_t bin = binStart; bin <= binStop; bin++){
 947
 948     Double_t pt=hPionRatio->GetBinCenter(bin);
 949
 950     cout << "Making projection for bin: " << bin << endl;
 951
 952     const Int_t j = bin-binStart;
 953
 954     TH1D* hDeltaPiVsPtProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeltaPiVsPtProj%d", bin), bin, bin);
 955     //    hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(-25, 20);
 956     hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40, 90);
 957     hDeltaPiVsPtProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
 958                                     hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
 959                                     hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
 960
 961     dEdxpoints[j]=0;
 962     dEdxpoints[j]=hDeltaPiVsPtProj;
 963     fPion[j]=0;
 964     fPion[j]=new TF1(Form("fPiongauss_%d",bin), "gausn", 40, 100);
 965
 966     fKaon[j]=0;
 967     fKaon[j]=new TF1(Form("fKaongauss_%d",bin), "gausn", 40, 100);
 968
 969     fProton[j]=0;
 970     fProton[j]=new TF1(Form("fProtongauss_%d",bin), "gausn", 40, 100);
 971
 972     fElectron[j]=0;
 973     fElectron[j]=new TF1(Form("fElectrongauss_%d",bin), "gausn", 40, 100);
 974
 975     fTotal[j]=0;
 976     fTotal[j]=new TF1(Form("fTotalgauss_%d",bin), "([0]+[12])*gausn(1)+[4]*gausn(5)+[8]*gausn(9)+[12]*gausn(13)", 40, 100);
 977
 978
 979     Double_t all =  hDeltaPiVsPtProj->GetEntries();
 980
 981
 982     const Int_t nPar = 16;
 983     Double_t gausParams[nPar] = {
 984       0.59,
 985       all,
 986       piFunc->Eval(pt),
 987       //fpi->Eval(pt),
 988       sigmaFunc->Eval(piFunc->Eval(pt)),
 989
 990       0.3,
 991       all,
 992       kFunc->Eval(pt),
 993       //fp->Eval(pt),
 994       sigmaFunc->Eval(kFunc->Eval(pt)),
 995
 996       0.1,
 997       all,
 998       pFunc->Eval(pt),
 999       //fp->Eval(pt),
1000       sigmaFunc->Eval(pFunc->Eval(pt)),
1001
1002       0.01,
1003       all,
1004       eFunc->Eval(pt),
1005       //fpi->Eval(pt),
1006       sigmaFunc->Eval(eFunc->Eval(pt)),
1007     };
1008
1009
1010     for(Int_t ipar=0;ipar<nPar;++ipar)
1011       cout<<gausParams[ipar]<<endl;
1012
1013
1014
1015     cFits->cd();
1016     cFits->cd(j + 1);
1017
1018     total->SetParameters(gausParams);
1019
1020
1021     for(Int_t i = 0; i < nPar; i++) {
1022       if(i!=0 && i!=4 && i!=8 && i!=12)
1023         total->FixParameter(i, gausParams[i]);
1024       else
1025         continue;
1026     }
1027
1028     total->SetParLimits(7, gausParams[7]-0.05*gausParams[7],gausParams[7]+0.05*gausParams[7]);
1029     total->SetParLimits(8, 0.005,0.4);
1030     total->SetParLimits(4, 0.005,0.6);
1031
1032     if(bin==48) {
1033       if(etaHigh[i_eta]==8||etaHigh[i_eta]==-6)
1034         hElectronRatio->Fit(fElectronFraction, "N", "", 4.0, 10.0);
1035       if(etaHigh[i_eta]==6||etaHigh[i_eta]==-4)
1036         hElectronRatio->Fit(fElectronFraction, "N", "", 3.66, 10.0);
1037       if(etaHigh[i_eta]==4||etaHigh[i_eta]==-2)
1038         hElectronRatio->Fit(fElectronFraction, "N", "", 3.33, 10.0);
1039       if(etaHigh[i_eta]==2||etaHigh[i_eta]==0)
1040         hElectronRatio->Fit(fElectronFraction, "N", "", 3.0, 10.0);
1041       fElectronFraction->SetRange(0, ptStop);
1042     }
1043
1044     if(bin>=48) {
1045       total->FixParameter(12, fElectronFraction->Eval(hDeDxVsP->GetXaxis()->GetBinCenter(bin)));
1046     }
1047
1048
1049     hDeltaPiVsPtProj->Fit(total, "0L");
1050
1051
1052     hDeltaPiVsPtProj->Draw();
1053     total->DrawCopy("same");
1054
1055     Double_t parametersOut[nPar];
1056     total->GetParameters(parametersOut);
1057     const Double_t* parameterErrorsOut = total->GetParErrors();
1058
1059     for(Int_t i = 0; i < nPar; i++)
1060       cout << parametersOut[i] << ", ";
1061     cout << endl;
1062     fTotal[j]->SetParameters(parametersOut);
1063
1064
1065
1066     pion->SetParameters(&parametersOut[1]);
1067     pion->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
1068
1069     fPion[j]->SetParameters(&parametersOut[1]);;
1070     fPion[j]->SetParameter(0,all*(parametersOut[0]+parametersOut[12]));
1071
1072
1073     hPionRatio->SetBinContent(bin, parametersOut[0]);
1074     hPionRatio->SetBinError(bin, parameterErrorsOut[0]);
1075     hPionYield->SetBinContent(bin, parametersOut[0]*all);
1076     hPionYield->SetBinError(bin, parameterErrorsOut[0]*all);
1077
1078     kaon->SetParameters(&parametersOut[5]);
1079     kaon->SetParameter(0,all*parametersOut[4]);
1080
1081     fKaon[j]->SetParameters(&parametersOut[5]);
1082     fKaon[j]->SetParameter(0,all*parametersOut[4]);
1083
1084
1085     hKaonRatio->SetBinContent(bin, parametersOut[4]);
1086     hKaonRatio->SetBinError(bin, parameterErrorsOut[4]);
1087     hKaonYield->SetBinContent(bin, parametersOut[4]*all);
1088     hKaonYield->SetBinError(bin, parameterErrorsOut[4]*all);
1089
1090     proton->SetParameters(&parametersOut[9]);
1091     proton->SetParameter(0,all*parametersOut[8]);
1092     //proton->DrawCopy("same");
1093
1094     fProton[j]->SetParameters(&parametersOut[9]);
1095     fProton[j]->SetParameter(0,all*parametersOut[8]);
1096
1097
1098     hProtonRatio->SetBinContent(bin, parametersOut[8]);
1099     hProtonRatio->SetBinError(bin, parameterErrorsOut[8]);
1100     hProtonYield->SetBinContent(bin, parametersOut[8]*all);
1101     hProtonYield->SetBinError(bin, parameterErrorsOut[8]*all);
1102
1103     electron->SetParameters(&parametersOut[13]);
1104     electron->SetParameter(0,all*parametersOut[12]);
1105
1106     fElectron[j]->SetParameters(&parametersOut[13]);
1107     fElectron[j]->SetParameter(0,all*parametersOut[12]);
1108
1109
1110     //electron->DrawCopy("same");
1111     hElectronRatio->SetBinContent(bin, parametersOut[12]);
1112     hElectronRatio->SetBinError(bin, parameterErrorsOut[12]);
1113     hElectronYield->SetBinContent(bin, parametersOut[12]*all);
1114     hElectronYield->SetBinError(bin, parameterErrorsOut[12]*all);
1115
1116     //DrawALICELogo(kFALSE, 0.71, 0.76, 0.82, 0.98);
1117     cSingleFit->cd();
1118     cSingleFit->Clear();
1119     cSingleFit->SetTickx(1);
1120     cSingleFit->SetTicky(1);
1121     cSingleFit->SetLogy(0);
1122
1123     //    cSingleFit->SetLogy();
1124     hDeltaPiVsPtProj->GetXaxis()->SetRangeUser(40,90);
1125     hDeltaPiVsPtProj->SetMarkerStyle(25);
1126     hDeltaPiVsPtProj->SetMarkerColor(1);
1127     hDeltaPiVsPtProj->SetLineColor(1);
1128     hDeltaPiVsPtProj->SetLineWidth(2);
1129     hDeltaPiVsPtProj->GetYaxis()->SetLabelFont(42);
1130     hDeltaPiVsPtProj->GetXaxis()->SetLabelFont(42);
1131     hDeltaPiVsPtProj->GetZaxis()->SetLabelFont(42);
1132     hDeltaPiVsPtProj->GetYaxis()->SetTitleFont(42);
1133     hDeltaPiVsPtProj->GetXaxis()->SetTitleFont(42);
1134     hDeltaPiVsPtProj->GetZaxis()->SetTitleFont(42);
1135     hDeltaPiVsPtProj->GetYaxis()->SetTitle("Entries");
1136     hDeltaPiVsPtProj->GetXaxis()->SetTitle("d#it{E}/d#it{x} in TPC (arb. units)");
1137     hDeltaPiVsPtProj->GetXaxis()->SetTitleSize(0.05);
1138     hDeltaPiVsPtProj->GetXaxis()->SetTitleOffset(0.9);
1139     hDeltaPiVsPtProj->GetYaxis()->SetTitleSize(0.05);
1140     hDeltaPiVsPtProj->GetYaxis()->SetTitleOffset(0.9);
1141     hDeltaPiVsPtProj->Draw();
1142
1143     total->SetLineColor(kGray+1);
1144     total->SetLineWidth(3);
1145     total->SetLineStyle(1);
1146     total->DrawCopy("same");
1147
1148     pion->GetXaxis()->SetRangeUser(40,90);
1149     pion->SetLineColor(2);
1150     pion->SetLineWidth(2);
1151     //pion->SetLineStyle(2);
1152     pion->DrawCopy("same");
1153
1154     kaon->SetLineColor(kGreen+2);
1155     kaon->SetLineWidth(2);
1156     //kaon->SetLineStyle(2);
1157     kaon->DrawCopy("same");
1158
1159
1160     proton->SetLineColor(4);
1161     proton->SetLineWidth(2);
1162     //proton->SetLineStyle(2);
1163     proton->DrawCopy("same");
1164
1165
1166     TLatex* latex = new TLatex();
1167     //  latex->SetTextFont(132);
1168     latex->SetNDC();
1169     latex->SetTextAlign(22);
1170
1171     latex->SetTextSize(0.05);
1172     latex->SetTextFont(42);
1173     latex->SetTextColor(kRed+2);
1174     //latex->DrawLatex(0.8,0.3,"pp");
1175     //latex->DrawLatex(0.81,0.25,"#sqrt{s}=2.76 TeV");
1176     latex->DrawLatex(0.8,0.25+0.035,Form("p-Pb, %s",CentLatex[i_cent]));
1177     latex->DrawLatex(0.81,0.2+0.035,"#sqrt{s_{NN}}=5.02 TeV");
1178
1179     latex->SetTextColor(1);
1180     latex->SetTextSize(0.04);
1181     //latex->DrawLatex(0.81,0.2,"25/07/2012");
1182     //latex->DrawLatex(0.80,0.3+0.036,"31/10/2012");
1183
1184     latex->SetTextSize(0.05);
1185     latex->DrawLatex(0.8,0.71,legEtaCut[i_eta]);
1186     latex->DrawLatex(0.8,0.64,Form("%.1f<#it{p}<%.1f GeV/#it{c}", hDeDxVsP->GetXaxis()->GetBinLowEdge(bin), hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
1187
1188
1189
1190     //cSingleFit->SaveAs(Form("%s/ptspectrum_bin%d.gif", dirName, bin));
1191     cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.png", dirName, endingCent[i_cent],bin));
1192     cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.pdf", dirName, endingCent[i_cent],bin));
1193     cSingleFit->SaveAs(Form("fit_results/%s_%s/ptspectrum_bin%d.eps", dirName, endingCent[i_cent],bin));
1194     //    legend->Draw();
1195
1196
1197
1198
1199   }
1200
1201
1202   TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p; pt", 600, 400);
1203
1204   cRatio->Clear();
1205   cRatio->SetGridy();
1206   hElectronRatio->GetYaxis()->SetRangeUser(-0.05, 0.1);
1207   hElectronRatio->DrawCopy("P E");
1208   fElectronFraction->DrawCopy("SAME");
1209
1210   gROOT->ProcessLine(".x drawText.C");
1211   cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.gif", dirName,endingCent[i_cent]));
1212   cRatio->SaveAs(Form("fit_results/%s_%s/electron_ratio.pdf", dirName,endingCent[i_cent]));
1213
1214   cRatio->Clear();
1215   cRatio->SetGridy(kFALSE);
1216   hPionRatio->SetMarkerStyle(20);
1217   hPionRatio->SetMarkerColor(2);
1218   hPionRatio->SetLineColor(2);
1219   hPionRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1220   hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
1221   hPionRatio->DrawCopy("P E");
1222   hPionYield->SetMarkerStyle(20);
1223   hPionYield->SetMarkerColor(2);
1224   hPionYield->SetLineColor(2);
1225   hPionYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1226
1227   hKaonRatio->SetMarkerStyle(20);
1228   hKaonRatio->SetMarkerColor(kGreen);
1229   hKaonRatio->SetLineColor(kGreen);
1230   hKaonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1231   hKaonYield->SetMarkerStyle(20);
1232   hKaonYield->SetMarkerColor(kGreen);
1233   hKaonYield->SetLineColor(kGreen);
1234   hKaonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1235   hKaonRatio->DrawCopy("SAME P E");
1236
1237   hProtonRatio->SetMarkerStyle(20);
1238   hProtonRatio->SetMarkerColor(4);
1239   hProtonRatio->SetLineColor(4);
1240   hProtonRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1241   hProtonYield->SetMarkerStyle(20);
1242   hProtonYield->SetMarkerColor(4);
1243   hProtonYield->SetLineColor(4);
1244   hProtonYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1245   hProtonRatio->DrawCopy("SAME P E");
1246
1247   hElectronRatio->SetMarkerStyle(20);
1248   hElectronRatio->SetMarkerColor(6);
1249   hElectronRatio->SetLineColor(6);
1250   hElectronRatio->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1251   hElectronRatio->DrawCopy("SAME P E");
1252   hElectronYield->SetMarkerStyle(20);
1253   hElectronYield->SetMarkerColor(6);
1254   hElectronYield->SetLineColor(6);
1255   hElectronYield->GetXaxis()->SetRangeUser(0.0, ptStop-0.001);
1256
1257   TH1D *hallRec=0;
1258
1259   hallRec=(TH1D *)hPionRatio->Clone("hRecAll");
1260   hallRec->Add(hKaonRatio);
1261   hallRec->Add(hProtonRatio);
1262   hallRec->SetLineColor(kYellow-3);
1263   hallRec->SetMarkerColor(kYellow-3);
1264   hallRec->DrawCopy("SAME P E");
1265
1266
1267   gROOT->ProcessLine(".x drawText.C");
1268   cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.gif", dirName,endingCent[i_cent]));
1269   cRatio->SaveAs(Form("fit_results/%s_%s/particle_ratios.pdf", dirName,endingCent[i_cent]));
1270
1271
1272
1273   if(outFileName) {
1274
1275     CreateDir(Form("fit_results/fit_yields_results_%s",endingCent[i_cent]));
1276     TFile* fileOut = new TFile(Form("fit_results/fit_yields_results_%s/%s%s.root", endingCent[i_cent], outFileName, ending[i_eta]), "RECREATE");
1277
1278
1279     hPionRatio->Write();
1280     hKaonRatio->Write();
1281     hProtonRatio->Write();
1282     hElectronRatio->Write();
1283
1284     hPionYield->Write();
1285     hKaonYield->Write();
1286     hProtonYield->Write();
1287     hElectronYield->Write();
1288
1289
1290
1291     fElectronFraction->Write();
1292
1293     for(Int_t bin = binStart; bin <= binStop; bin++){
1294
1295       Int_t j = bin-binStart;
1296       dEdxpoints[j]->Write();
1297       fPion[j]->Write();
1298       fKaon[j]->Write();
1299       fProton[j]->Write();
1300       fElectron[j]->Write();
1301       fTotal[j]->Write();
1302
1303
1304     }
1305     fmip->Write();
1306     fileOut->Close();
1307   }
1308
1309
1310
1311
1312 }
1313
1314
1315
1316
1317
1318
1319 //____________________________________________________________________________
1320 void FitDeDxVsP(const Char_t * dedxFileName,
1321                 Double_t pStart, Double_t pStop,  Int_t i_cent,
1322                 Int_t optionDeDx, Int_t optionSigma,
1323                 Bool_t performFit = kFALSE,
1324                 Int_t etaLow=0, Int_t etaHigh=8,
1325                 Bool_t fixAllParBB=kFALSE,
1326                 Bool_t fixAllParSigma=kFALSE,
1327                 Bool_t fixKtoZero=kFALSE,
1328                 const Char_t* initialFitFileName=0,
1329                 Int_t fixParametersDedx=-1,
1330                 Int_t fixParametersSigma=-1)
1331 {
1332   //gStyle->SetOptStat(0);
1333
1334
1335
1336   TFile* dedxFile = FindFileFresh(dedxFileName);
1337   if(!dedxFile)
1338     return;
1339
1340   TList * list = (TList *)dedxFile->Get("outputdedx");
1341   TH2D *hDeDxVsPPlus = 0;
1342   hDeDxVsPPlus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaLow, etaHigh));
1343   hDeDxVsPPlus->Sumw2();
1344
1345   TH2D *hDeDxVsPMinus = 0;
1346   hDeDxVsPMinus = (TH2D *)list->FindObject(Form("histAllCh%d%d", etaHigh, etaLow));
1347   hDeDxVsPMinus->Sumw2();
1348
1349   gSystem->Exec("rm *.root");
1350   TFile *fplus = new TFile("hist2Dplus.root","recreate");
1351   fplus->cd();
1352   hDeDxVsPPlus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
1353   hDeDxVsPPlus->Write();
1354   fplus->Close();
1355
1356   TFile *fminus = new TFile("hist2Dminus.root","recreate");
1357   fminus->cd();
1358   hDeDxVsPMinus->SetName(Form("histAllCh%d%d", etaLow, etaHigh));
1359   hDeDxVsPMinus->Write();
1360   fminus->Close();
1361
1362   gSystem->Exec("hadd hist2Ddedx.root hist2Dplus.root hist2Dminus.root");
1363
1364   TFile *ffinal = TFile::Open("hist2Ddedx.root");
1365   hDeDxVsP = (TH2D *)ffinal->Get(Form("histAllCh%d%d", etaLow, etaHigh));
1366
1367
1368
1369   CreateDir("old");
1370   gSystem->Exec(Form("mv results/calibdedx_%d%d/* old/",etaLow, etaHigh));
1371
1372
1373   TCanvas* cDeDxVsP = new TCanvas("cDeDxVsP", "dE/dx vs p", 400, 300);
1374   cDeDxVsP->Clear();
1375   cDeDxVsP->cd();
1376   cDeDxVsP->SetLogz();
1377   hDeDxVsP->SetTitle(0);
1378   hDeDxVsP->Draw("COLZ");
1379
1380   TCanvas* cDeDxVsPLogX = new TCanvas("cDeDxVsPLogX", "dE/dx vs p", 400, 300);
1381   cDeDxVsPLogX->Clear();
1382   cDeDxVsPLogX->cd();
1383   cDeDxVsPLogX->SetLogz();
1384   cDeDxVsPLogX->SetLogx();
1385   hDeDxVsP->Draw("COLZ");
1386
1387   // Root is a bit stupid with finidng bins so we have to add and subtract a
1388   // little to be sure we get the right bin as we typically put edges as
1389   // limits
1390   const Int_t binStart = hDeDxVsP->GetXaxis()->FindBin(pStart+0.01);
1391   pStart = hDeDxVsP->GetXaxis()->GetBinLowEdge(binStart);
1392   const Int_t binStop  = hDeDxVsP->GetXaxis()->FindBin(pStop-0.01);
1393   pStop = hDeDxVsP->GetXaxis()->GetBinUpEdge(binStop);
1394   const Int_t nBins    = binStop - binStart + 1;
1395
1396   cout << "Doing 2d fit from pTlow = " << pStart << " (bin: " << binStart
1397        << ") to pThigh = " << pStop << " (bin: " << binStop << ")" << endl;
1398
1399   // Double_t binSize = (histdEdxvsP->GetXaxis()->GetXmax() - histdEdxvsP->GetXaxis()->GetXmin())/ histdEdxvsP->GetXaxis()->GetNbins();
1400
1401   Double_t parDeDx[6]  = {0, 0, 0, 0, 0, 0};
1402   Double_t parSigma[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1403
1404   const Int_t nLocalParams  = 3; // pi, K, p yields
1405   Int_t nDeDxPar      = 0;
1406   Int_t nSigmaPar     = 0;
1407
1408   switch(optionDeDx) {
1409
1410   case 1:
1411     nDeDxPar = 2;
1412     parDeDx[0] = 39.7;
1413     parDeDx[1] =  6.3;
1414     break;
1415   case 2:
1416     nDeDxPar = 1;
1417     parDeDx[0] =  7.3;
1418     break;
1419   case 3:
1420     nDeDxPar = 2;
1421     parDeDx[0] =  6.85097;
1422     parDeDx[1] =  29.5933;
1423     break;
1424   case 4:
1425     nDeDxPar = 3;
1426     parDeDx[0] =  31.2435;
1427     parDeDx[1] =  9.73037;
1428     parDeDx[2] =  1.95832;
1429     break;
1430   case 5:
1431     nDeDxPar = 4;
1432     parDeDx[0] =  31.384;
1433     parDeDx[1] =  9.5253;
1434     parDeDx[2] =  7.3542;
1435     parDeDx[3] =  1.4;
1436     break;
1437   case 6:
1438     nDeDxPar = 5;
1439     parDeDx[0] =  31.384;
1440     parDeDx[1] =  9.5253;
1441     parDeDx[2] =  7.3542;
1442     parDeDx[3] =  1.5;
1443     parDeDx[4] =  81;
1444     break;
1445   case 7:
1446     nDeDxPar = 3;
1447     parDeDx[0] =  31.3;
1448     parDeDx[1] =  9.5;
1449     parDeDx[2] =  1.4;
1450     break;
1451   default:
1452
1453     cout << "optionDeDx does not support option: " << optionDeDx << endl;
1454     return;
1455     break;
1456   }
1457
1458   switch(optionSigma) {
1459
1460   case 1:
1461     nSigmaPar = 1;
1462     parSigma[0] = 3.0;
1463     break;
1464   case 2:
1465     nSigmaPar = 1;
1466     parSigma[0] = 0.0655688;
1467     break;
1468   case 3:
1469     nSigmaPar = 2;
1470     parSigma[0] = 0.06;
1471     parSigma[1] = -0.001;
1472   case 4:
1473     nSigmaPar = 2;
1474     parSigma[0] = 0.06;
1475     parSigma[1] = 1.0;
1476     break;
1477   case 5:
1478     nSigmaPar = 2;
1479     parSigma[0] = 0.06;
1480     parSigma[1] = 1.0;
1481     break;
1482   case 6:
1483     nSigmaPar = 2;
1484     parSigma[0] = 0.06;
1485     parSigma[1] = 0.0;
1486     break;
1487   case 7:
1488     nSigmaPar = 3;
1489     parSigma[0] = 0.06;
1490     parSigma[1] = 0.0;
1491     parSigma[2] = 2.0;
1492     break;
1493   case 8:
1494     nSigmaPar = 3;
1495     parSigma[0] = 0.06;
1496     parSigma[1] = 0.0;
1497     parSigma[2] = 2.0;
1498     break;
1499   case 9://works for 0-5 %
1500     nSigmaPar = 3;
1501     parSigma[0] = 1.88409e-01;
1502     parSigma[1] = -3.84073e-03;
1503     parSigma[2] = 3.03110e-05;
1504     break;
1505   case 10://works for 0-5 %
1506     nSigmaPar = 6;
1507     parSigma[0]=7.96380e-03;
1508     parSigma[1]=8.09869e-05;
1509     parSigma[2]=6.29707e-06;
1510     parSigma[3]=-3.27295e-08;
1511     parSigma[4]=-1.20200e+03;
1512     parSigma[5]=-3.97089e+01;
1513     break;
1514   case 11://works for 0-5 %
1515     nSigmaPar = 6;
1516     parSigma[0]=7.96380e-03;
1517     parSigma[1]=8.09869e-05;
1518     parSigma[2]=6.29707e-06;
1519     parSigma[3]=-3.27295e-08;
1520     parSigma[4]=-1.20200e+03;
1521     parSigma[5]=-3.97089e+01;
1522     break;
1523   case 12://works for 0-5 %
1524     nSigmaPar = 3;
1525     parSigma[0]=7.96380e-03;
1526     parSigma[1]=8.09869e-05;
1527     parSigma[2]=6.29707e-06;
1528     break;
1529   case 13://works for 0-5 %
1530     nSigmaPar = 3;
1531     parSigma[0]=7.96380e-03;
1532     parSigma[1]=8.09869e-05;
1533     parSigma[2]=6.29707e-06;
1534     break;
1535   case 14:
1536     nSigmaPar = 2;
1537     parSigma[0]=7.96380e-03;
1538     parSigma[1]=8.09869e-05;
1539     break;
1540   case 15:
1541     nSigmaPar = 3;
1542     parSigma[0]=7.96380e-03;
1543     parSigma[1]=8.09869e-05;
1544     break;
1545
1546
1547
1548
1549   default:
1550
1551     cout << "optionSigma does not support option: " << optionSigma << endl;
1552     return;
1553     break;
1554   }
1555
1556   if(initialFitFileName) {
1557
1558
1559
1560     TFile* fresBB = FindFileFresh(initialFitFileName);
1561     if(!fresBB)
1562       return;
1563     TF1 *fres=(TF1 *)fresBB->Get("sigmaRes");
1564     TF1 *fBB=(TF1 *)fresBB->Get("dedxFunc");
1565
1566     Int_t nparres=fres->GetNpar();
1567     Int_t nparBB=fBB->GetNpar();
1568     for(Int_t ires=0;ires<nparres;++ires){
1569       parSigma[ires] = fres->GetParameter(ires+1);
1570       //cout<<fres->GetParameter(ires+1)<<endl;
1571     }
1572     for(Int_t iBB=0;iBB<nparBB;++iBB){
1573       parDeDx[iBB] = fBB->GetParameter(iBB+1);
1574       cout<<parDeDx[iBB]<<endl;
1575
1576     }
1577     //Esto lo puse apenas
1578     fixPlateau  = parDeDx[4];
1579     fixMIP = fBB->Eval(3.5);
1580     cout<<"Hello!!"<<"   fixPlateau="<<fixPlateau<<"  fixMIP="<<fixMIP<<endl;
1581
1582     //fixMIP      = fitPar->MIP;
1583     //return;
1584   }
1585
1586
1587
1588   const Int_t nGlobalParams = 2  // binStart, binStop,
1589     + 2 + nDeDxPar               // optionDeDx, nDeDxPar, dedxpar0 ....
1590     + 2 + nSigmaPar;             // nSigmaPar, sigmapar0 .....
1591
1592   const Int_t nParams = nBins*nLocalParams + nGlobalParams;
1593
1594
1595   TF2* fitAll = new TF2("fitAll", fitf3G, pStart, pStop, 30, 90, nParams);
1596   Double_t parametersIn[nParams];
1597
1598   parametersIn[0] = binStart;
1599   fitAll->SetParName(0,"binStart");
1600   fitAll->FixParameter(0, parametersIn[0]);
1601
1602   parametersIn[1] = binStop;
1603   fitAll->SetParName(1,"binStop");
1604   fitAll->FixParameter(1, parametersIn[1]);
1605
1606   // dE/dx parameters
1607   parametersIn[2] = nDeDxPar;
1608   fitAll->SetParName(2,"nDeDxPar");
1609   fitAll->FixParameter(2, parametersIn[2]);
1610
1611   parametersIn[3] = optionDeDx;
1612   fitAll->SetParName(3,"optionDeDx");
1613   fitAll->FixParameter(3, parametersIn[3]);
1614
1615   for(Int_t i = 0; i < nDeDxPar; i++) {
1616
1617     parametersIn[4+i] = parDeDx[i];
1618     fitAll->SetParName(4+i,Form("dEdXpar%d", i));
1619     Int_t bit = TMath::Nint(TMath::Power(2, i));
1620     if(fixParametersDedx>=0)
1621       if (fixParametersDedx==0 || fixParametersDedx&bit) {
1622
1623         cout << "Fixing dE/dx parameter " << i << endl;
1624         fitAll->FixParameter(4+i, parametersIn[4+i]);
1625       }
1626
1627     if(fixAllParBB) {
1628
1629       fitAll->FixParameter(4+i, parametersIn[4+i]);
1630       }
1631   }
1632
1633   // sigma parameters
1634   parametersIn[4+nDeDxPar] = nSigmaPar;
1635   fitAll->SetParName(4+nDeDxPar,"nSigmaPar");
1636   fitAll->FixParameter(4+nDeDxPar, parametersIn[4+nDeDxPar]);
1637
1638   parametersIn[5+nDeDxPar] = optionSigma;
1639   fitAll->SetParName(5+nDeDxPar,"optionSigma");
1640   fitAll->FixParameter(5+nDeDxPar, parametersIn[5+nDeDxPar]);
1641
1642   for(Int_t i = 0; i < nSigmaPar; i++) {
1643
1644     parametersIn[6+nDeDxPar+i] = parSigma[i];
1645     fitAll->SetParName(6+nDeDxPar+i,Form("sigmapar%d", i));
1646     Int_t bit = TMath::Nint(TMath::Power(2, i));
1647     if(fixParametersSigma>=0)
1648       if (fixParametersSigma==0 || fixParametersSigma&bit) {
1649
1650         cout << "Fixing sigma parameter " << i << endl;
1651         fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
1652       }
1653
1654     if(fixAllParSigma) {
1655
1656       fitAll->FixParameter(6+nDeDxPar+i, parametersIn[6+nDeDxPar+i]);
1657       }
1658   }
1659
1660   // Initial yields
1661
1662   for(Int_t bin = binStart; bin <= binStop; bin++) {
1663
1664     TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY("hDeDxVsPProj", bin, bin);
1665
1666     const Int_t offset = nGlobalParams + (bin - binStart)*nLocalParams;
1667     const Double_t all = hDeDxVsPProj->Integral();
1668
1669     fitAll->SetParName(offset + 0, Form("piYield%d", bin));
1670     fitAll->SetParName(offset + 1, Form("kYield%d", bin));
1671     fitAll->SetParName(offset + 2, Form("pYield%d", bin));
1672     fitAll->SetParLimits(offset + 0, 0, 10*all);
1673     fitAll->SetParLimits(offset + 2, 0, 10*all);
1674
1675     if(fixKtoZero) {
1676       parametersIn[offset + 0] = (all)*0.5;
1677       parametersIn[offset + 1] = (all)*0.0;
1678       fitAll->FixParameter(offset + 1, 0.0);
1679       parametersIn[offset + 2] = (all)*0.5;
1680     } else {
1681       parametersIn[offset + 0] = (all)*0.6;
1682       parametersIn[offset + 1] = (all)*0.3;
1683       fitAll->SetParLimits(offset + 1, 0, 10*all);
1684       parametersIn[offset + 2] = (all)*0.1;
1685     }
1686     // fitAll->SetParLimits(offset + 0, 0., histdEdxvsPpy->GetEntries());
1687     // fitAll->SetParLimits(offset + 1, 0., histdEdxvsPpy->GetEntries());
1688     // fitAll->SetParLimits(offset + 2, 0., histdEdxvsPpy->GetEntries());
1689   }
1690
1691   fitAll->SetParameters(parametersIn);
1692   fitAll->Print();
1693
1694   Bool_t converged = kFALSE;
1695
1696   if(performFit) {
1697     for(Int_t i = 0; i < 4; i++) {
1698       TFitResultPtr fitResultPtr =  hDeDxVsP->Fit(fitAll, "0S", "", pStart+0.01, pStop-0.01);
1699       if(!fitResultPtr->Status()) {
1700         //      if(!fitResultPtr->Status() && !fitResultPtr->CovMatrixStatus()) {
1701
1702         converged = kTRUE;
1703         break;
1704       }
1705     }
1706   }
1707   // else we just draw how the results looks with the input parameters
1708
1709   Double_t parametersOut[nParams];
1710   fitAll->GetParameters(parametersOut);
1711   const Double_t* parameterErrorsOut = fitAll->GetParErrors();
1712
1713   // overlay the fitfunction
1714
1715
1716   TF1* fDeDxPi = new TF1("fDeDxPi", FitFunc, 0, 50, nDeDxPar+1); // +1 because of option!
1717   fDeDxPi->SetParameters(&parametersOut[3]);
1718   fDeDxPi->SetLineWidth(2);
1719   cDeDxVsP->cd();
1720   fDeDxPi->Draw("SAME");
1721   cDeDxVsPLogX->cd();
1722   fDeDxPi->Draw("SAME");
1723
1724   TF1* fDeDxK = new TF1("fDeDxK", FitFunc, 0, 50, nDeDxPar+1);
1725   fDeDxK->SetParameters(&parametersOut[3]);
1726   fDeDxK->SetParameter(0, fDeDxK->GetParameter(0)+10);
1727   fDeDxK->SetLineWidth(2);
1728   cDeDxVsP->cd();
1729   fDeDxK->Draw("SAME");
1730   cDeDxVsPLogX->cd();
1731   fDeDxK->Draw("SAME");
1732
1733   TF1* fDeDxP = new TF1("fDeDxP", FitFunc, 0, 50, nDeDxPar+1);
1734   fDeDxP->SetParameters(&parametersOut[3]);
1735   fDeDxP->SetParameter(0, fDeDxP->GetParameter(0)+20);
1736   fDeDxP->SetLineWidth(2);
1737   cDeDxVsP->cd();
1738   fDeDxP->Draw("SAME");
1739   cDeDxVsPLogX->cd();
1740   fDeDxP->Draw("SAME");
1741
1742   TF1* fDeDxE = new TF1("fDeDxE", FitFunc, 0, 50, nDeDxPar+1);
1743   fDeDxE->SetParameters(&parametersOut[3]);
1744   fDeDxE->SetParameter(0, fDeDxE->GetParameter(0)+30);
1745   fDeDxE->SetLineWidth(2);
1746   cDeDxVsP->cd();
1747   fDeDxE->Draw("SAME");
1748   cDeDxVsPLogX->cd();
1749   fDeDxE->Draw("SAME");
1750
1751   TF1* fSigma = new TF1("fSigma", SigmaFunc, 0, 50, nSigmaPar+1);
1752   fSigma->SetParameters(&parametersOut[5 + nDeDxPar]);
1753
1754   //  fitAll->Draw("same cont3");
1755
1756   CreateDir(Form("results/%s/calibdedx_%d%d",endingCent[i_cent],etaLow,etaHigh));
1757
1758   //CreateDir(Form("results/calibdedx_%d%d",etaLow,etaHigh));
1759
1760   cDeDxVsP->cd();
1761   cDeDxVsP->Modified();
1762   cDeDxVsP->Update();
1763   gROOT->ProcessLine(".x drawText.C");
1764   cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.gif",endingCent[i_cent],etaLow,etaHigh));
1765   cDeDxVsP->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p.pdf",endingCent[i_cent],etaLow,etaHigh));
1766
1767   cDeDxVsPLogX->cd();
1768   cDeDxVsPLogX->Modified();
1769   cDeDxVsPLogX->Update();
1770   gROOT->ProcessLine(".x drawText.C");
1771   cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.gif",endingCent[i_cent],etaLow,etaHigh));
1772   cDeDxVsPLogX->SaveAs(Form("results/%s/calibdedx_%d%d/dedx_vs_p_logx.pdf",endingCent[i_cent],etaLow,etaHigh));
1773
1774   //cross check
1775   TCanvas* cFits = new TCanvas("cFits", "Fit comparison to data", 1200, 800);
1776
1777   cFits->Clear();
1778   cFits->Divide(7, 5);
1779
1780   TF1* pion = new TF1("pion", "gausn", 30, 90);
1781   pion->SetLineWidth(2);
1782   pion->SetLineColor(kRed);
1783   TF1* kaon = new TF1("kaon", "gausn", 30, 90);
1784   kaon->SetLineWidth(2);
1785   kaon->SetLineColor(kGreen);
1786   TF1* proton = new TF1("proton", "gausn", 30, 90);
1787   proton->SetLineWidth(2);
1788   proton->SetLineColor(kBlue);
1789   TF1* total = new TF1("total", "gausn(0)+gausn(3)+gausn(6)", 30, 90);
1790   total->SetLineColor(kBlack);
1791   total->SetLineWidth(2);
1792   total->SetLineStyle(2);
1793
1794   TLegend* legend = new TLegend(0.11, 0.6, 0.35, 0.85);
1795   legend->SetBorderSize(0);
1796   legend->SetFillColor(0);
1797   legend->AddEntry(total, "3-Gauss fit", "L");
1798   legend->AddEntry(pion, "#pi", "L");
1799   legend->AddEntry(kaon, "K", "L");
1800   legend->AddEntry(proton, "p", "L");
1801
1802
1803   TCanvas* cSingleFit = new TCanvas("cSingleFit", "single fit");
1804   cSingleFit->SetLogy(1);
1805
1806   TH1D* hPionRatio =(TH1D*)hDeDxVsP->ProjectionX("hPionRatio", 1, 1);
1807   hPionRatio->Reset();
1808   hPionRatio->GetXaxis()->SetRangeUser(pStart+0.001, pStop-0.001);
1809   hPionRatio->GetYaxis()->SetRangeUser(0.0, 1.0);
1810   hPionRatio->SetTitle("particle fractions; p [GeV/c]; particle fraction");
1811   TH1D* hKaonRatio   = (TH1D*)hPionRatio->Clone("hKaonRatio");
1812   TH1D* hProtonRatio = (TH1D*)hPionRatio->Clone("hProtonRatio");
1813
1814   for(Int_t bin = binStart; bin <= binStop; bin++){
1815
1816     cout << "Making projection for bin: " << bin << endl;
1817
1818     const Int_t j = bin-binStart;
1819     cFits->cd();
1820     cFits->cd(j + 1);
1821
1822     TH1D* hDeDxVsPProj =(TH1D*)hDeDxVsP->ProjectionY(Form("hDeDxVsPProj%d", bin), bin, bin);
1823     //    hDeDxVsPProj->GetXaxis()->SetRangeUser(40, 90);
1824     hDeDxVsPProj->SetTitle(Form("%.1f<p<%.1f GeV/c",
1825                                 hDeDxVsP->GetXaxis()->GetBinLowEdge(bin),
1826                                 hDeDxVsP->GetXaxis()->GetBinUpEdge(bin)));
1827     hDeDxVsPProj->Draw();
1828
1829     const Int_t offset = nGlobalParams + j*nLocalParams;
1830     const Double_t p = hDeDxVsP->GetXaxis()->GetBinCenter(bin);
1831     const Double_t pKeff = p*piMass/kMass; // corresponding p of a pion with same dE/dx
1832     const Double_t pPeff = p*piMass/pMass; // corresponding p of a pion with same dE/dx
1833     const Double_t meanDeDxPi = fDeDxPi->Eval(p);
1834     const Double_t meanDeDxK  = fDeDxPi->Eval(pKeff);
1835     const Double_t meanDeDxP  = fDeDxPi->Eval(pPeff);
1836     Double_t gausParams[9] = {
1837       parametersOut[offset + 0],
1838       meanDeDxPi,
1839       fSigma->Eval(meanDeDxPi) ,
1840       parametersOut[offset + 1],
1841       meanDeDxK,
1842       fSigma->Eval(meanDeDxK) ,
1843       parametersOut[offset + 2],
1844       meanDeDxP,
1845       fSigma->Eval(meanDeDxP) ,
1846     };
1847
1848     for(Int_t i = 0; i < 9; i++)
1849       cout << gausParams[i] << ", ";
1850
1851     cout << endl;
1852
1853     pion->SetParameters(&gausParams[0]);
1854     pion->DrawCopy("same");
1855     Double_t all =  hDeDxVsPProj->Integral();
1856     hPionRatio->SetBinContent(bin, parametersOut[offset + 0]/all);
1857     hPionRatio->SetBinError(bin, parameterErrorsOut[offset + 0]/all);
1858
1859     kaon->SetParameters(&gausParams[3]);
1860     kaon->DrawCopy("same");
1861     hKaonRatio->SetBinContent(bin, parametersOut[offset + 1]/all);
1862     hKaonRatio->SetBinError(bin, parameterErrorsOut[offset + 1]/all);
1863
1864     proton->SetParameters(&gausParams[6]);
1865     proton->DrawCopy("same");
1866     hProtonRatio->SetBinContent(bin, parametersOut[offset + 2]/all);
1867     hProtonRatio->SetBinError(bin, parameterErrorsOut[offset + 2]/all);
1868
1869     total->SetParameters(gausParams);
1870     total->DrawCopy("same");
1871
1872     cSingleFit->cd();
1873     cSingleFit->Clear();
1874     //    cSingleFit->SetLogy();
1875     hDeDxVsPProj->Draw();
1876     pion->DrawCopy("same");
1877     kaon->DrawCopy("same");
1878     proton->DrawCopy("same");
1879     total->DrawCopy("same");
1880
1881     gROOT->ProcessLine(".x drawText.C(2)");
1882
1883
1884     cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.gif",endingCent[i_cent], etaLow,etaHigh, bin));
1885     cSingleFit->SaveAs(Form("results/%s/calibdedx_%d%d/ptspectrum_bin%d.pdf",endingCent[i_cent],etaLow, etaHigh, bin));
1886     //    legend->Draw();
1887
1888
1889   }
1890
1891   TCanvas* cRatio = new TCanvas("cRatio", "ratios/all vs p", 600, 400);
1892   cRatio->Clear();
1893   hPionRatio->SetMaximum(0.8);
1894   hPionRatio->SetMarkerStyle(20);
1895   hPionRatio->SetMarkerColor(2);
1896   hPionRatio->Draw("P E");
1897
1898   hKaonRatio->SetMarkerStyle(20);
1899   hKaonRatio->SetMarkerColor(3);
1900   hKaonRatio->Draw("SAME P E");
1901
1902   hProtonRatio->SetMarkerStyle(20);
1903   hProtonRatio->SetMarkerColor(4);
1904   hProtonRatio->Draw("SAME P E");
1905   gROOT->ProcessLine(".x drawText.C(2)");
1906   cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.gif",endingCent[i_cent],etaLow,etaHigh));
1907   cRatio->SaveAs(Form("results/%s/calibdedx_%d%d/particle_ratios.pdf",endingCent[i_cent],etaLow,etaHigh));
1908
1909
1910
1911   //
1912   // Store the <dE/dx> parameters in a file that we can get them back to use for the Delta-pi!
1913   //
1914   DeDxFitInfo* fitInfo = new DeDxFitInfo();
1915   fitInfo->MIP     = fixMIP;
1916   //fitInfo->plateau = 80;
1917   fitInfo->optionDeDx = optionDeDx;
1918   fitInfo->nDeDxPar = nDeDxPar;
1919   for(Int_t i = 0; i < nDeDxPar; i++) {
1920     fitInfo->parDeDx[i] = fDeDxPi->GetParameter(i+1); // 1st parameter is option
1921   }
1922   //fitInfo->plateau = fDeDxPi->GetParameter(5);//lo puse 16/09/13
1923   fitInfo->optionSigma = optionSigma;
1924   fitInfo->nSigmaPar = nSigmaPar;
1925   for(Int_t i = 0; i < nSigmaPar; i++) {
1926     cout<<"my sugma param="<<fSigma->GetParameter(i+1)<<endl;
1927     fitInfo->parSigma[i] = fSigma->GetParameter(i+1); // 1st parameter is option
1928   }
1929
1930
1931   fitInfo->Print();
1932
1933
1934   if(converged) {
1935
1936     cout << "Fit converged and error matrix was ok" << endl;
1937   } else {
1938
1939     cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was  not ok!" << endl;
1940   }
1941
1942
1943   CreateDir(Form("fitparameters/%s",endingCent[i_cent]));
1944   //cout<<"flag 0"<<"    gSystem->BaseName(calibFileName))="<< gSystem->BaseName(calibFileName)<<endl;
1945   TFile* outFile = new TFile(Form("fitparameters/%s/%d%d_dataPbPb.root", endingCent[i_cent], etaLow, etaHigh), "RECREATE");
1946   outFile->cd();
1947   fitInfo->Write("fitInfo");
1948   cout<<"flag 1"<<endl;
1949
1950   outFile->Close();
1951
1952   if(converged) {
1953
1954     cout << "Fit converged and error matrix was ok" << endl;
1955   } else {
1956
1957     cout << "WARNING!!!!!!!!!!!!!!! Fit did not converge, or error matrix was  not ok!" << endl;
1958   }
1959 }
1960
1961
1962
1963
1964
1965 void MakeFitsV0s(const Char_t * fileInputName, const Char_t * fileE, const Char_t *outDir, Int_t index_eta){
1966   /*
1967   gStyle->SetPalette(1);
1968   gStyle->SetCanvasColor(10);
1969   gStyle->SetFrameFillColor(10);
1970   gStyle->SetOptStat(0);
1971   gStyle->SetOptTitle(0);
1972   gStyle->SetOptFit(0);
1973   */
1974   TFile *f12=TFile::Open(fileE);
1975   TF1 *fPlateau = (TF1 *)f12->Get(Form("fGauss%d",index_eta));
1976   /*
1977   TH1D * hEV0s=(TH1D *)f12->Get("hdedx_vs_p_e");
1978   TGraphErrors  *gResE=(TGraphErrors    *)f12->Get("gRes_Electrons");
1979   */
1980
1981   const Char_t * endNamePos[4]={"02","24","46","68"};
1982   const Char_t * endNameNeg[4]={"20","42","64","86"};
1983
1984
1985
1986   TFile *f1=TFile::Open(fileInputName);
1987   TList * list = (TList *)f1->Get("outputdedx");
1988
1989
1990   //positive eta
1991   TH2D *hPi2DPos=(TH2D *)list->FindObject(Form("histPiTof%s",endNamePos[index_eta]));
1992   TH2D *hPi2DV0sPos=(TH2D *)list->FindObject(Form("histPiV0%s",endNamePos[index_eta]));
1993   TH2D *hP2DPos=(TH2D *)list->FindObject(Form("histPV0%s",endNamePos[index_eta]));
1994   TH1D *hPi2DpPos=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNamePos[index_eta]));
1995   TH1D *hPi2DV0spPos=(TH1D *)list->FindObject(Form("histPiV01D%s",endNamePos[index_eta]));
1996   TH1D *hP2DpPos=(TH1D *)list->FindObject(Form("histPV01D%s",endNamePos[index_eta]));
1997
1998
1999   gSystem->Exec("rm *.root");
2000
2001   TFile *fpos = new TFile("histtmppos.root","recreate");
2002   fpos->cd();
2003   hPi2DPos->SetName(Form("histPiTof%s",endNamePos[index_eta]));
2004   hPi2DV0sPos->SetName(Form("histPiV0%s",endNamePos[index_eta]));
2005   hP2DPos->SetName(Form("histPV0%s",endNamePos[index_eta]));
2006   hPi2DpPos->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
2007   hPi2DV0spPos->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
2008   hP2DpPos->SetName(Form("histPV01D%s",endNamePos[index_eta]));
2009   hPi2DPos->Write();
2010   hPi2DV0sPos->Write();
2011   hP2DPos->Write();
2012   hPi2DpPos->Write();
2013   hPi2DV0spPos->Write();
2014   hP2DpPos->Write();
2015   fpos->Close();
2016
2017
2018   //negative eta
2019   TH2D *hPi2DNeg=(TH2D *)list->FindObject(Form("histPiTof%s",endNameNeg[index_eta]));
2020   TH2D *hPi2DV0sNeg=(TH2D *)list->FindObject(Form("histPiV0%s",endNameNeg[index_eta]));
2021   TH2D *hP2DNeg=(TH2D *)list->FindObject(Form("histPV0%s",endNameNeg[index_eta]));
2022   TH1D *hPi2DpNeg=(TH1D *)list->FindObject(Form("histPiTof1D%s",endNameNeg[index_eta]));
2023   TH1D *hPi2DV0spNeg=(TH1D *)list->FindObject(Form("histPiV01D%s",endNameNeg[index_eta]));
2024   TH1D *hP2DpNeg=(TH1D *)list->FindObject(Form("histPV01D%s",endNameNeg[index_eta]));
2025
2026
2027
2028   TFile *fneg = new TFile("histtmpneg.root","recreate");
2029   fneg->cd();
2030   hPi2DNeg->SetName(Form("histPiTof%s",endNamePos[index_eta]));
2031   hPi2DV0sNeg->SetName(Form("histPiV0%s",endNamePos[index_eta]));
2032   hP2DNeg->SetName(Form("histPV0%s",endNamePos[index_eta]));
2033   hPi2DpNeg->SetName(Form("histPiTof1D%s",endNamePos[index_eta]));
2034   hPi2DV0spNeg->SetName(Form("histPiV01D%s",endNamePos[index_eta]));
2035   hP2DpNeg->SetName(Form("histPV01D%s",endNamePos[index_eta]));
2036   hPi2DNeg->Write();
2037   hPi2DV0sNeg->Write();
2038   hP2DNeg->Write();
2039   hPi2DpNeg->Write();
2040   hPi2DV0spNeg->Write();
2041   hP2DpNeg->Write();
2042   fneg->Close();
2043
2044   //pos+neg eta
2045   TH2D *hPi2D=0;
2046   TH2D *hPi2DV0s=0;
2047   TH2D *hP2D=0;
2048   TH1D *hPi2Dp=0;
2049   TH1D *hPi2DV0sp=0;
2050   TH1D *hP2Dp=0;
2051
2052
2053   gSystem->Exec(Form("hadd tmp%d.root histtmppos.root histtmpneg.root",index_eta));
2054
2055   TFile *ffinalv0 = TFile::Open(Form("tmp%d.root",index_eta));
2056   hPi2D=(TH2D *)ffinalv0->Get(Form("histPiTof%s",endNamePos[index_eta]));
2057   hPi2DV0s=(TH2D *)ffinalv0->Get(Form("histPiV0%s",endNamePos[index_eta]));
2058   hP2D=(TH2D *)ffinalv0->Get(Form("histPV0%s",endNamePos[index_eta]));
2059   hPi2Dp=(TH1D *)ffinalv0->Get(Form("histPiTof1D%s",endNamePos[index_eta]));
2060   hPi2DV0sp=(TH1D *)ffinalv0->Get(Form("histPiV01D%s",endNamePos[index_eta]));
2061   hP2Dp=(TH1D *)ffinalv0->Get(Form("histPV01D%s",endNamePos[index_eta]));
2062
2063
2064
2065   TList *l1=new TList();
2066   l1->SetOwner();
2067
2068
2069   TGraphErrors* graphRes = new TGraphErrors();
2070   Int_t nERes = 0;
2071
2072   TGraphErrors* graphBB = new TGraphErrors();
2073   Int_t nBB = 0;
2074
2075   TGraphErrors* graphBBpitof = new TGraphErrors();
2076   Int_t nBBpitof = 0;
2077   TGraphErrors* graphBBpiv0s = new TGraphErrors();
2078   Int_t nBBpiv0s = 0;
2079   TGraphErrors* graphBBpv0s = new TGraphErrors();
2080   Int_t nBBpv0s = 0;
2081
2082
2083
2084
2085   TF1* fFit1D = new TF1("fFit1D", "gausn(0)", 40, 90);
2086
2087
2088
2089   //TF1* pionsv0s = new TF1("pionsv0s", "gausn", 56, 80);
2090   TF1* pionsv0s = new TF1("pionsv0s", "gausn", 40, 80);
2091   pionsv0s->SetLineColor(kRed);
2092   pionsv0s->SetLineWidth(2);
2093   pionsv0s->SetLineStyle(2);
2094
2095   //TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 55);
2096   TF1* protonsv0s = new TF1("protonsv0s", "gausn", 40, 80);
2097   protonsv0s->SetLineColor(kBlue);
2098   protonsv0s->SetLineWidth(2);
2099   protonsv0s->SetLineStyle(2);
2100
2101
2102   TF1* total = new TF1("total", "gausn(0)+gausn(3)", 40, 80);
2103   total->SetLineColor(1);
2104   total->SetLineWidth(2);
2105   total->SetLineStyle(2);
2106
2107
2108
2109   Int_t nbins=hPi2D->GetXaxis()->GetNbins();
2110
2111
2112
2113   for(Int_t bin=1;bin<=nbins;++bin){
2114
2115
2116
2117     TH1D *hproy_pion=0;
2118     TH1D *hproy_proton=0;
2119
2120     Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
2121     Double_t lowp = hPi2D->GetXaxis()->GetBinLowEdge(bin);
2122     Double_t widthp = hPi2D->GetXaxis()->GetBinWidth(bin);
2123     hPi2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2124     Double_t meanp= hPi2Dp->GetMean();
2125     Double_t meanpE= hPi2Dp->GetMeanError();
2126
2127
2128     cout<<"p="<<p<<endl;
2129
2130     if(p>2.0&&p<9.0){
2131
2132       cout<<bin<<endl;
2133
2134       hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pionTof_%d",bin),bin,bin);
2135
2136       cout<<bin<<endl;
2137
2138       hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2139                                 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2140                                 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2141       cout<<bin<<endl;
2142       hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2143
2144       cout<<"nbinsproj="<<hproy_proton->GetNbinsX()<<endl;
2145       //hproy_pion->Draw();
2146
2147
2148       TSpectrum *s = new TSpectrum(2);
2149       Int_t nfound = s->Search(hproy_pion,2,"",0.15);
2150       cout<<"!!!!!!!!!!!!!!         nfound="<<nfound<<endl;
2151       Int_t npeaks = 0;
2152       Float_t *xpeaks = s->GetPositionX();
2153
2154
2155       for (Int_t p1=0;p1<nfound;p1++) {
2156         //Float_t xp = xpeaks[p1];
2157         npeaks++;
2158       }
2159
2160       if(npeaks==0)
2161         continue;
2162       if(npeaks==1){
2163
2164         fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
2165         hproy_pion->Fit(fFit1D, "Q");
2166         hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
2167                         xpeaks[0]+hproy_pion->GetRMS());
2168
2169         hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2170                         fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2171
2172
2173
2174
2175           graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2176           graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2177           nBB++;
2178
2179           graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2180           graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2181           nBBpitof++;
2182
2183           if(p>3&&p<9){
2184             graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2185             graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2186             nERes++;
2187           }
2188
2189
2190
2191
2192       }
2193       else if(npeaks==2){
2194         Int_t bin1 = hproy_pion->GetXaxis()->FindBin(xpeaks[0]);
2195         Int_t bin2 = hproy_pion->GetXaxis()->FindBin(xpeaks[1]);
2196
2197         Int_t bin1C = hproy_pion->GetBinContent(bin1);
2198         Int_t bin2C = hproy_pion->GetBinContent(bin2);
2199
2200         cout<<"bin1="<<bin1<<"  content="<<bin1C<<endl;
2201         cout<<"bin2="<<bin2<<"  content="<<bin2C<<endl;
2202
2203
2204         if(bin==24){
2205
2206           hproy_pion->Fit(fFit1D, "R", "", 65, 85);
2207           graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2208           graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2209           nBB++;
2210
2211           graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2212           graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2213           nBBpitof++;
2214
2215
2216
2217         }
2218
2219
2220         else if(bin1C>bin2C){//pion signal larger
2221
2222
2223           fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[0], hproy_pion->GetRMS());
2224           hproy_pion->Fit(fFit1D, "Q");
2225           hproy_pion->Fit(fFit1D, "Q", "", xpeaks[0]-hproy_pion->GetRMS(),
2226                           xpeaks[0]+hproy_pion->GetRMS());
2227
2228
2229
2230
2231           hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2232                           fFit1D->GetParameter(1)+2.5*fFit1D->GetParameter(2));
2233
2234
2235
2236
2237           cout<<"hello"<<endl;
2238           //protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], hproy_proton->GetRMS());
2239
2240
2241           hproy_proton->Fit(protonsv0s, "Q");
2242
2243
2244           hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[1]-hproy_proton->GetRMS(),
2245                             xpeaks[1]+hproy_proton->GetRMS());
2246
2247
2248           hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2*protonsv0s->GetParameter(2),
2249                           xpeaks[1]+2*protonsv0s->GetParameter(2));
2250
2251           protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[1], protonsv0s->GetParameter(2));
2252           hproy_proton->Fit(protonsv0s, "R", "", xpeaks[1]-2.0*protonsv0s->GetParameter(2),
2253                             xpeaks[1]+2.0*protonsv0s->GetParameter(2));
2254
2255           total->SetParameter(1,protonsv0s->GetParameter(1));
2256           total->SetParameter(2,protonsv0s->GetParameter(2));
2257           total->SetParameter(4,fFit1D->GetParameter(1));
2258           total->SetParameter(5,fFit1D->GetParameter(2));
2259           hproy_pion->Fit(total,"R");
2260
2261           Float_t reducedCh2=total->GetChisquare()/total->GetNDF();
2262           if(reducedCh2>1.0){
2263             hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
2264                             fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2265             //hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-0.5*fFit1D->GetParameter(2),
2266             hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2267                             fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2268
2269
2270             graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2271             graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2272             nBB++;
2273
2274
2275             graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2276             graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2277             nBBpitof++;
2278
2279             if(p>3.0&&p<9){
2280               graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2281               graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2282               nERes++;
2283             }
2284
2285
2286
2287           }else{
2288
2289             graphBB->SetPoint(nBB,meanp/0.14, total->GetParameter(4));
2290             graphBB->SetPointError(nBB,meanpE, total->GetParError(4));
2291             nBB++;
2292
2293             graphBBpitof->SetPoint(nBBpitof,meanp/0.14, total->GetParameter(4));
2294             graphBBpitof->SetPointError(nBBpitof,meanpE, total->GetParError(4));
2295             nBBpitof++;
2296
2297             if(p>3.0&&p<9){
2298               graphRes->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
2299
2300               graphRes->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
2301               nERes++;
2302             }
2303
2304
2305
2306
2307           }
2308
2309
2310
2311
2312         }else{
2313
2314           fFit1D->SetParameters(hproy_pion->GetEntries(), xpeaks[1], hproy_pion->GetRMS());
2315           hproy_pion->Fit(fFit1D, "Q");
2316           hproy_pion->Fit(fFit1D, "Q", "", xpeaks[1]-hproy_pion->GetRMS(),
2317                           xpeaks[0]+hproy_pion->GetRMS());
2318
2319           hproy_pion->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-2*fFit1D->GetParameter(2),
2320                           fFit1D->GetParameter(1)+2*fFit1D->GetParameter(2));
2321
2322
2323           protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], hproy_proton->GetRMS());
2324           hproy_proton->Fit(protonsv0s, "Q");
2325           hproy_proton->Fit(protonsv0s, "Q", "", xpeaks[0]-hproy_proton->GetRMS(),
2326                           xpeaks[0]+hproy_proton->GetRMS());
2327
2328           hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
2329                           xpeaks[0]+2*protonsv0s->GetParameter(2));
2330
2331           protonsv0s->SetParameters(hproy_proton->GetEntries(), xpeaks[0], protonsv0s->GetParameter(2));
2332           hproy_proton->Fit(protonsv0s, "R", "", xpeaks[0]-2*protonsv0s->GetParameter(2),
2333                           xpeaks[0]+2*protonsv0s->GetParameter(2));
2334
2335
2336           //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
2337           //graphBB->SetPointError(nBB,0, fFit1D->GetParError(1));
2338           //nBB++;
2339           total->SetParameter(1,protonsv0s->GetParameter(1));
2340           total->SetParameter(2,protonsv0s->GetParameter(2));
2341           total->SetParameter(4,fFit1D->GetParameter(1));
2342           total->SetParameter(5,fFit1D->GetParameter(2));
2343           hproy_pion->Fit(total,"R");
2344
2345           //graphBB->SetPoint(nBB,p/0.14, xpeaks[1]);
2346           //graphBB->SetPointError(nBB,0, total->GetParError(4));
2347           graphBB->SetPoint(nBB,meanp/0.14, fFit1D->GetParameter(1));
2348           graphBB->SetPointError(nBB,meanpE, fFit1D->GetParError(1));
2349           nBB++;
2350
2351           graphBBpitof->SetPoint(nBBpitof,meanp/0.14, fFit1D->GetParameter(1));
2352           graphBBpitof->SetPointError(nBBpitof,meanpE, fFit1D->GetParError(1));
2353           nBBpitof++;
2354
2355
2356           if(p>3.0&&p<9){
2357             graphRes->SetPoint(nERes,fFit1D->GetParameter(1),fFit1D->GetParameter(2)/fFit1D->GetParameter(1));
2358             graphRes->SetPointError(nERes,fFit1D->GetParError(1),fFit1D->GetParError(2)/fFit1D->GetParameter(1));
2359             nERes++;
2360           }
2361
2362
2363
2364         }
2365
2366
2367
2368
2369
2370
2371
2372
2373
2374
2375
2376       }
2377
2378       l1->Add(hproy_pion);
2379
2380
2381     }
2382     else
2383       continue;
2384
2385   }
2386
2387
2388
2389   graphRes->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
2390   graphRes->SetName("gRes_PionsTof");
2391
2392   graphBBpitof->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
2393   graphBBpitof->SetName("BB_piTOF");
2394   graphBBpitof->SetMarkerColor(kGreen+2);
2395   graphBBpitof->SetLineColor(kGreen+2);
2396   graphBBpitof->SetMarkerStyle(21);
2397
2398   l1->Add(graphRes);
2399   l1->Add(graphBBpitof);
2400
2401
2402
2403   //pi by V0s
2404
2405   TGraphErrors* graphResV0s = new TGraphErrors();
2406   nERes = 0;
2407
2408   nbins=hPi2DV0s->GetXaxis()->GetNbins();
2409
2410
2411
2412   for(Int_t bin=1;bin<=nbins;++bin){
2413
2414
2415
2416     TH1D *hproy_pionV0s=0;
2417     TH1D *hproy_proton=0;
2418
2419     Double_t p=hPi2DV0s->GetXaxis()->GetBinCenter(bin);
2420
2421     Double_t lowp = hPi2DV0s->GetXaxis()->GetBinLowEdge(bin);
2422     Double_t widthp = hPi2DV0s->GetXaxis()->GetBinWidth(bin);
2423     hPi2DV0sp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2424     Double_t meanp= hPi2DV0sp->GetMean();
2425     Double_t meanpE= hPi2DV0sp->GetMeanError();
2426
2427
2428
2429     if(p>2.0&&p<9){
2430
2431       hproy_pionV0s=(TH1D *)hPi2DV0s->ProjectionY(Form("hproy_pionV0s_%d",bin),bin,bin);
2432       hproy_pionV0s->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2433                                 hPi2DV0s->GetXaxis()->GetBinLowEdge(bin),
2434                                 hPi2DV0s->GetXaxis()->GetBinUpEdge(bin)));
2435
2436      fFit1D->SetParameters(hproy_pionV0s->GetEntries(), hproy_pionV0s->GetMean(), hproy_pionV0s->GetRMS());
2437       hproy_pionV0s->Fit(fFit1D, "Q");
2438       hproy_pionV0s->Fit(fFit1D, "Q", "", hproy_pionV0s->GetMean()-hproy_pionV0s->GetRMS(),
2439                          hproy_pionV0s->GetMean()+hproy_pionV0s->GetRMS());
2440
2441
2442       //      hproy_pionV0s->Fit(fFit1D, "Q", "", fFit1D->GetParameter(1)-1.5*fFit1D->GetParameter(2),
2443       //              fFit1D->GetParameter(1)+3.0*fFit1D->GetParameter(2));
2444
2445       hproy_pionV0s->Fit(fFit1D, "R", "", fFit1D->GetParameter(1)-1.0*fFit1D->GetParameter(2),
2446                          fFit1D->GetParameter(1)+2.0*fFit1D->GetParameter(2));
2447
2448
2449       hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2450       hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx; Entries",
2451                                   hPi2D->GetXaxis()->GetBinLowEdge(bin),
2452                                   hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2453
2454       //protonsv0s->SetParameters(hproy_proton->GetEntries(), hproy_proton->GetMean(), hproy_proton->GetRMS());
2455       hproy_proton->Fit(protonsv0s,"R");
2456       //hproy_proton->Fit(protonsv0s,"R", "",protonsv0s->GetParameter(1)-3.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+1.0*protonsv0s->GetParameter(2));
2457       hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
2458
2459       total->SetParameter(0,0.1*hproy_pionV0s->GetEntries());
2460       total->SetParameter(1,protonsv0s->GetParameter(1));
2461       total->SetParameter(2,protonsv0s->GetParameter(2));
2462       total->SetParameter(3,0.9*hproy_pionV0s->GetEntries());
2463       total->SetParameter(4,fFit1D->GetParameter(1));
2464       total->SetParameter(5,fFit1D->GetParameter(2));
2465
2466       hproy_pionV0s->Fit(total,"0R");
2467
2468       cout<<"p="<<p<<"   yield0="<<total->GetParameter(0)<<"  yield1="<<total->GetParameter(3)<<endl;
2469       cout<<"p="<<p<<"   mean0="<<total->GetParameter(1)<<"  mean1="<<total->GetParameter(4)<<endl;
2470       cout<<"p="<<p<<"   sigma0="<<total->GetParameter(1)<<"  sigma1="<<total->GetParameter(5)<<endl;
2471
2472       total->SetParameter(1,total->GetParameter(1));
2473       total->SetParameter(2,total->GetParameter(2));
2474       total->SetParameter(4,total->GetParameter(4));
2475       total->SetParameter(5,total->GetParameter(5));
2476
2477       hproy_pionV0s->Fit(total,"R");
2478
2479
2480
2481       if(p>4.0){
2482
2483
2484         //return;
2485
2486         graphResV0s->SetPoint(nERes,total->GetParameter(4),total->GetParameter(5)/total->GetParameter(4));
2487         graphResV0s->SetPointError(nERes,total->GetParError(4),total->GetParError(5)/total->GetParameter(4));
2488         nERes++;
2489
2490         //graphBB->SetPoint(nBB,p/0.14, total->GetParameter(4));
2491         //graphBB->SetPointError(nBB,0, total->GetParError(4));
2492         //nBB++;
2493
2494
2495         graphBBpiv0s->SetPoint(nBBpiv0s,meanp/0.14, total->GetParameter(4));
2496         graphBBpiv0s->SetPointError(nBBpiv0s,meanpE,total->GetParError(4));
2497         nBBpiv0s++;
2498
2499         l1->Add(hproy_pionV0s);
2500
2501       }
2502
2503
2504
2505     }
2506
2507   }
2508
2509
2510   graphResV0s->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
2511   graphResV0s->SetName("gRes_PionsV0sArmenteros");
2512   l1->Add(graphResV0s);
2513
2514   graphBBpiv0s->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
2515   graphBBpiv0s->SetName("BB_piV0s");
2516   graphBBpiv0s->SetMarkerColor(2);
2517   graphBBpiv0s->SetLineColor(2);
2518   graphBBpiv0s->SetMarkerStyle(24);
2519   l1->Add(graphBBpiv0s);
2520
2521   //protons
2522
2523
2524
2525
2526   TGraphErrors* graphPRes = new TGraphErrors();
2527   nERes = 0;
2528   TGraphErrors* graphdEdxvsP = new TGraphErrors();
2529   TGraphErrors* graphSigmavsP = new TGraphErrors();
2530   graphdEdxvsP->SetName("graphdEdxvsP_p");
2531   graphSigmavsP->SetName("graphSigmavsP_p");
2532
2533   nbins=hPi2D->GetXaxis()->GetNbins();
2534   for(Int_t bin=1;bin<=nbins;++bin){
2535
2536     TH1D *hproy_pion=0;
2537     TH1D *hproy_proton=0;
2538
2539
2540     Double_t p=hPi2D->GetXaxis()->GetBinCenter(bin);
2541     Double_t lowp = hP2D->GetXaxis()->GetBinLowEdge(bin);
2542     Double_t widthp = hP2D->GetXaxis()->GetBinWidth(bin);
2543     hP2Dp->GetXaxis()->SetRangeUser(lowp,lowp+widthp);
2544     Double_t meanp= hP2Dp->GetMean();
2545     Double_t meanpE= hP2Dp->GetMeanError();
2546
2547
2548     if(p<2.0||p>5.0)
2549       continue;
2550
2551
2552     if(p>2){
2553
2554       cout<<"p="<<hPi2D->GetXaxis()->GetBinCenter(bin)<<endl;
2555
2556       hproy_pion=(TH1D *)hPi2D->ProjectionY(Form("hproy_pion_%d",bin),bin,bin);
2557       hproy_pion->SetTitle(Form("pionsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2558                                 hPi2D->GetXaxis()->GetBinLowEdge(bin),
2559                                 hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2560
2561       hproy_pion->Fit(pionsv0s,"R");
2562
2563       hproy_proton=(TH1D *)hP2D->ProjectionY(Form("hproy_proton_%d",bin),bin,bin);
2564       hproy_proton->SetTitle(Form("protonsv0s: %.1f < p < %.1f GeV/c; dE/dx (GeV/c); Entries",
2565                                   hPi2D->GetXaxis()->GetBinLowEdge(bin),
2566                                   hPi2D->GetXaxis()->GetBinUpEdge(bin)));
2567
2568
2569
2570       hproy_proton->Fit(protonsv0s,"R");
2571
2572       hproy_proton->Fit(protonsv0s,"R","",protonsv0s->GetParameter(1)-2.0*protonsv0s->GetParameter(2),protonsv0s->GetParameter(1)+0.8*protonsv0s->GetParameter(2));
2573
2574
2575       Float_t reducedCh2=protonsv0s->GetChisquare()/protonsv0s->GetNDF();
2576       //if(reducedCh2<3)
2577       if(reducedCh2>3){
2578         total->SetParameter(1,protonsv0s->GetParameter(1));
2579         total->SetParameter(2,protonsv0s->GetParameter(2));
2580         total->SetParameter(4,pionsv0s->GetParameter(1));
2581         total->SetParameter(5,pionsv0s->GetParameter(2));
2582
2583         hproy_proton->Fit(total,"R");
2584       }
2585       //if(total->GetParameter(1)<=0)continue;
2586
2587       if(p>2&&p<7){
2588         if(reducedCh2<5){
2589
2590
2591
2592           if(p>3.0&&p<5){
2593             graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
2594             graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
2595             nERes++;
2596           }
2597           graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
2598           graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
2599
2600           graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
2601           graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
2602           nBBpv0s++;
2603
2604           graphdEdxvsP->SetPoint(nBB,p, protonsv0s->GetParameter(1));
2605           graphdEdxvsP->SetPointError(nBB,0.001, protonsv0s->GetParError(1));
2606           graphSigmavsP->SetPoint(nBB,p, protonsv0s->GetParameter(2));
2607           graphSigmavsP->SetPointError(nBB,0.0001, protonsv0s->GetParError(2));
2608           nBB++;
2609
2610         }
2611         else{
2612
2613           if(p>3.0&&p<5){
2614             graphPRes->SetPoint(nERes,total->GetParameter(1),total->GetParameter(2)/total->GetParameter(1));
2615             graphPRes->SetPointError(nERes,total->GetParError(1),total->GetParError(2)/total->GetParameter(1));
2616             nERes++;
2617           }
2618
2619           //graphBB->SetPoint(nBB,p/0.943, total->GetParameter(1));
2620           //graphBB->SetPointError(nBB,0,total->GetParError(1));
2621
2622           graphBBpv0s->SetPoint(nBBpv0s,meanp/0.943, protonsv0s->GetParameter(1));
2623           graphBBpv0s->SetPointError(nBBpv0s,meanpE,protonsv0s->GetParError(1));
2624           nBBpv0s++;
2625
2626
2627           graphdEdxvsP->SetPoint(nBB,p, total->GetParameter(1));
2628           graphdEdxvsP->SetPointError(nBB,0.001, total->GetParError(1));
2629           graphSigmavsP->SetPoint(nBB,p, total->GetParameter(2));
2630           graphSigmavsP->SetPointError(nBB,0.0001, total->GetParError(2));
2631
2632           graphBB->SetPoint(nBB,meanp/0.943, protonsv0s->GetParameter(1));
2633           graphBB->SetPointError(nBB,meanpE,protonsv0s->GetParError(1));
2634
2635           nBB++;
2636         }
2637
2638
2639       }
2640       /*else{
2641
2642         graphPRes->SetPoint(nERes,protonsv0s->GetParameter(1),protonsv0s->GetParameter(2)/protonsv0s->GetParameter(1));
2643         graphPRes->SetPointError(nERes,protonsv0s->GetParError(1),protonsv0s->GetParError(2)/protonsv0s->GetParameter(1));
2644
2645
2646         nERes++;
2647
2648         }*/
2649
2650       l1->Add(hproy_proton);
2651
2652
2653     }
2654
2655
2656
2657   }
2658
2659   graphPRes->SetTitle(";#LT dE/dx  #GT;#sigma/#LT dE/dx  #GT");
2660   graphPRes->SetName("gRes_Protons");
2661
2662   l1->Add(graphPRes);
2663
2664   graphBBpv0s->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
2665   graphBBpv0s->SetName("BB_pv0s");
2666   graphBBpv0s->SetMarkerColor(4);
2667   graphBBpv0s->SetLineColor(4);
2668   graphBBpv0s->SetMarkerStyle(24);
2669   l1->Add(graphBBpv0s);
2670
2671
2672
2673   graphBB->SetTitle(";#beta*#gamma; #LT dE/dx  #GT");
2674   graphBB->SetName("gBB");
2675
2676
2677   Double_t dedxPar[10]  = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
2678
2679   dedxPar[0] =  0;
2680   dedxPar[1] =  31.384;
2681   dedxPar[2] =  9.5253;
2682   dedxPar[3] =  7.3542;
2683   dedxPar[4] =  1.5;
2684
2685   TF1 *dedxFunc = new TF1("dedxFunc", FitFunc, 2.5, 60, 6);
2686   dedxFunc->SetParameters(dedxPar);
2687   // Fit betagamma (option+40)
2688   dedxFunc->SetParameter(0, 6+40);
2689   dedxFunc->FixParameter(0, 6+40);
2690   //parameters for MB
2691   dedxFunc->SetParameter(1, 33.7538);
2692   dedxFunc->SetParameter(2, 8.05926);
2693   dedxFunc->SetParameter(3, 1.69954);
2694   dedxFunc->SetParameter(4, 1.37793);
2695   dedxFunc->SetParameter(5, 78.0465);
2696
2697
2698   cout << "!!!!!!!!!!!!!!!!!                    Setting new plateau=" << fPlateau->GetParameter(1) <<  endl;
2699   dedxFunc->SetParameter(5, fPlateau->GetParameter(1));
2700   dedxFunc->FixParameter(5, fPlateau->GetParameter(1));
2701
2702
2703   graphBB->Fit(dedxFunc, "0R", "SAME");
2704   graphBB->Fit(dedxFunc, "0R", "SAME");
2705   graphBB->Fit(dedxFunc, "0R", "SAME");
2706
2707   graphBB->Draw();
2708
2709   for(Int_t i=0;i<dedxFunc->GetNpar();++i)
2710     dedxFunc->SetParameter(i,dedxFunc->GetParameter(i));
2711
2712   l1->Add(dedxFunc);
2713   l1->Add(graphBB);
2714
2715
2716
2717   //Fit to get resolution
2718   TGraphErrors *g1=0;
2719   g1=new TGraphErrors();
2720   g1->SetName("res_allpions");
2721
2722   Int_t npionts=0;
2723
2724
2725   for(Int_t ipoint = 0; ipoint < graphPRes->GetN(); ipoint++){
2726     g1->SetPoint(npionts,graphPRes->GetX()[ipoint],graphPRes->GetY()[ipoint]);
2727     g1->SetPointError(npionts,graphPRes->GetEX()[ipoint],graphPRes->GetEY()[ipoint]);
2728     npionts++;
2729   }
2730
2731   //  pions TOF
2732
2733   for(Int_t ipoint = 0; ipoint < graphRes->GetN(); ipoint++){
2734
2735     g1->SetPoint(npionts,graphRes->GetX()[ipoint],graphRes->GetY()[ipoint]);
2736     g1->SetPointError(npionts,graphRes->GetEX()[ipoint],graphRes->GetEY()[ipoint]);
2737     npionts++;
2738     }
2739   //electrons
2740   g1->SetPoint(npionts,fPlateau->GetParameter(1),fPlateau->GetParameter(2)/fPlateau->GetParameter(1));
2741   g1->SetPointError(npionts,fPlateau->GetParError(1),fPlateau->GetParError(2)/fPlateau->GetParameter(1));
2742
2743
2744   Float_t sigmaPar[7];
2745   sigmaPar[0]=12;
2746   sigmaPar[1]=7.96380e-03;
2747   sigmaPar[2]=8.09869e-05;
2748   sigmaPar[3]=6.29707e-06;
2749   sigmaPar[4]=-3.27295e-08;
2750   sigmaPar[5]=-1.20200e+03;
2751   sigmaPar[6]=-3.97089e+01;
2752
2753
2754
2755   TF1 *sigmarrrrr=new TF1("sigmaRes",SigmaFunc,47.5, 85,4);
2756   sigmarrrrr->FixParameter(0,sigmaPar[0]);
2757   g1->Fit(sigmarrrrr, "0R", "SAME");
2758   g1->Fit(sigmarrrrr, "0R", "SAME");
2759   g1->Fit(sigmarrrrr, "0R", "SAME");
2760
2761   TF1 *fdedxvspP=new TF1("fdedxvspP","pol0",3.0,5.5);
2762   graphdEdxvsP->Fit(fdedxvspP, "R", "SAME");
2763   TF1 *fsigmavspP=new TF1("fsigmavspP","pol0",3.0,5.5);
2764   graphSigmavsP->Fit(fsigmavspP, "R", "SAME");
2765
2766
2767
2768   for(Int_t i=0;i<6;++i)
2769     sigmarrrrr->SetParameter(i,sigmarrrrr->GetParameter(i));
2770
2771   l1->Add(graphdEdxvsP);
2772   l1->Add(graphSigmavsP);
2773
2774   l1->Add(fdedxvspP);
2775   l1->Add(fsigmavspP);
2776
2777   l1->Add(sigmarrrrr);
2778   l1->Add(g1);
2779
2780   TFile *fout=new TFile(Form("%s/hres_0_5_%s.root",outDir,endNamePos[index_eta]),"recreate");
2781   fout->cd();
2782   l1->Write();
2783   fout->Close();
2784
2785
2786
2787
2788 }
2789
2790
2791
2792 void MakeFitsExternalData(const Char_t * inFile, const Char_t * outDir){
2793
2794
2795   TFile *fIn = TFile::Open(inFile);
2796   TList * list = (TList *)fIn->Get("outputdedx");
2797   //plateau, electrons 0.4<p<0.6 GeV/c
2798   TH2D *hEVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
2799   hEVsEta->Sumw2();
2800   //four th1, projections
2801   TH1D *hdEdxE[4]={0,0,0,0};// |eta|
2802   TH1D *hdEdxEpos[4]={0,0,0,0};// eta>0
2803   TF1 *fGaussE[4]={0,0,0,0};
2804   Int_t index_max = 1;
2805   for(Int_t i=0; i<4; ++i){
2806
2807     Int_t index_min_negeta = i+index_max;
2808     Int_t index_max_negeta = i+index_max+1;
2809     Int_t index_min_poseta = 17-(i+index_max+1);
2810     Int_t index_max_poseta = 17-(i+index_max);
2811
2812     hdEdxE[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxNegEta%d",i),index_min_negeta,index_max_negeta);
2813     hdEdxEpos[i] = (TH1D *)hEVsEta->ProjectionY(Form("hdEdxPosEta%d",i),index_min_poseta,index_max_poseta);
2814     hdEdxE[i]->Add(hdEdxEpos[i]);
2815
2816     //Fitting the signal
2817     fGaussE[i]=new TF1(Form("fGauss%d",i),"gaus");
2818     hdEdxE[i]->Fit(fGaussE[i],"0","");
2819     hdEdxE[i]->Fit(fGaussE[i],"","", fGaussE[i]->GetParameter(1)-1.5*fGaussE[i]->GetParameter(2), fGaussE[i]->GetParameter(1)+1.5*fGaussE[i]->GetParameter(2));
2820
2821     index_max++;
2822
2823   }
2824
2825   CreateDir(outDir);
2826
2827   TFile *fout = 0;
2828
2829   if(outDir){
2830
2831     fout = new TFile(Form("%s/PrimaryElectrons.root",outDir),"recreate");
2832     fout->cd();
2833     for(Int_t i=0; i<4; ++i){
2834       hdEdxE[i]->Write();
2835       fGaussE[i]->Write();
2836     }
2837     fout->Close();
2838
2839
2840   }
2841
2842
2843
2844 }
2845 //__________________________________________
2846 void PlotQA(const Char_t * inFile){
2847
2848   //gStyle->SetPalette(1);
2849   gStyle->SetCanvasColor(10);
2850   gStyle->SetFrameFillColor(10);
2851   gStyle->SetOptStat(0);
2852   gStyle->SetOptTitle(0);
2853   gStyle->SetOptFit(0);
2854
2855   const Char_t* ending[8] = { "86", "64", "42", "20", "02", "24", "46", "68" };
2856   const Char_t* LatexEta[8] = { "-0.8<#eta<-0.6", "-0.6<#eta<-0.4", "-0.4<#eta<-0.2", "-0.2<#eta<0",
2857                                 "0<#eta<0.2", "0.2<#eta<0.4", "0.4<#eta<0.6", "0.6<#eta<0.8" };
2858   TFile *fIn = TFile::Open(inFile);
2859   if(!fIn)
2860     return;
2861
2862   TList * list = (TList *)fIn->Get("outputdedx");
2863
2864
2865   TH2D *hMIP[8];
2866   TH2D *hPlateau[8];
2867   TProfile *pMIP[8];
2868   TProfile *pPlateau[8];
2869
2870   TH2D *hMIPVsNch[8];
2871   TProfile *pMIPVsNch[8];
2872
2873
2874   //Secondary pions
2875   TH2D *hPiV0s2D[8];
2876   TH1D *hPion3[8];
2877   TH1D *hPion2GeV = new TH1D("hPion2GeV","hPion2GeV",8,-0.8,0.8);
2878
2879
2880   for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
2881     hMIP[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsPhi%s",ending[index_eta]));
2882     hPlateau[index_eta] = (TH2D *)list->FindObject(Form("hPlateauVsPhi%s",ending[index_eta]));
2883     pMIP[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsPhi%s",ending[index_eta]));
2884     pPlateau[index_eta] = (TProfile *)list->FindObject(Form("pPlateauVsPhi%s",ending[index_eta]));
2885     hPiV0s2D[index_eta] = (TH2D *)list->FindObject(Form("histPiV0%s",ending[index_eta]));
2886     hPion3[index_eta] = (TH1D *)hPiV0s2D[index_eta]->ProjectionY(Form("Pion3%s",ending[index_eta]),16,16);
2887
2888     hMIPVsNch[index_eta] = (TH2D *)list->FindObject(Form("hMIPVsNch%s",ending[index_eta]));
2889     pMIPVsNch[index_eta] = (TProfile *)list->FindObject(Form("pMIPVsNch%s",ending[index_eta]));
2890
2891
2892
2893     TF1 *fgaus = 0;
2894     fgaus = new TF1("fgaus","gausn",40,80);
2895     hPion3[index_eta]->Fit(fgaus,"0","",40,80);
2896
2897     hPion2GeV->SetBinContent(index_eta+1,fgaus->GetParameter(1));
2898     hPion2GeV->SetBinError(index_eta+1,fgaus->GetParError(1));
2899   }
2900
2901
2902
2903   TH2D *hMIPVsEta = (TH2D *)list->FindObject("hMIPVsEta");
2904   TProfile *pMIPVsEta = (TProfile *)list->FindObject("pMIPVsEta");
2905   TProfile *pMIPVsEtaV0s = (TProfile *)list->FindObject("pMIPVsEtaV0s");
2906   TProfile *pPlateauVsEta = (TProfile *)list->FindObject("pPlateauVsEta");
2907   TH2D *hPlateauVsEta = (TH2D *)list->FindObject("hPlateauVsEta");
2908
2909
2910
2911   const Int_t nPadX = 4;
2912   const Int_t nPadY = 2;
2913   Float_t noMargin    = 0.000;
2914   Float_t topMargin   = 0.01;
2915   Float_t botMargin   = 0.1;
2916   Float_t leftMargin  = 0.05;
2917   Float_t rightMargin = 0.01;
2918   Float_t width       = (1-rightMargin-leftMargin)/nPadX;
2919   Float_t height      = (1-botMargin-topMargin)/nPadY;
2920   Float_t shift = 0.05;
2921
2922
2923   TCanvas* cMIP = new TCanvas("cMIP2", "Raa Pions", 900, 500);
2924   TPad* padMIP[nPadX*nPadY];
2925   //Float_t shift = 0.1;
2926   //Float_t shift = 0.05;
2927   const char* yTitleMIP = "d#it{E}/d#it{x}_{MIP}";
2928   const char* xTitleMIP = "#phi (rad)";
2929   cMIP->cd();
2930
2931
2932   TLatex* latex = new TLatex();
2933   latex->SetNDC();
2934   latex->SetTextAlign(22);
2935   latex->SetTextAngle(0);
2936   latex->SetTextFont(62);
2937
2938   latex->DrawLatex(0.53,0.04,xTitleMIP);
2939   latex->SetTextAngle(90);
2940   latex->SetTextSize(0.055);
2941   latex->DrawLatex(0.025,0.5,yTitleMIP);
2942   latex->SetTextSize(0.05);
2943
2944
2945   for (Int_t i = 1; i < 9; i++) {
2946
2947     Int_t iX = (i-1)%nPadX;
2948     Int_t iY = (i-1)/nPadX;
2949     Float_t x1 = leftMargin + iX*width;
2950     if(iX==2)
2951       x1-=0.01;
2952     else if(iX==1)
2953       x1+=0.01;
2954     Float_t x2 = leftMargin + (iX + 1)*width;
2955      if(iX==0)
2956        x2+=0.01;
2957      else if(iX==1)
2958        x2-=0.01;
2959      Float_t y1 = 1 - topMargin - (iY +1)*height;
2960      Float_t y2 = 1 - topMargin - iY*height;
2961      padMIP[i-1] = new TPad(Form("padRaa%d",i),"", x1, y1, x2, y2, 0, 0, 0);
2962      Float_t mBot = noMargin;
2963      Float_t mTop = noMargin;
2964      Float_t mLeft = noMargin;
2965      Float_t mRight = noMargin;
2966      if(iY==0)       mTop   = shift;
2967      if(iY==nPadY-1) mBot   = shift;
2968      if(iX==0)       mLeft  = 0.08;
2969      if(iX==nPadX-1) mRight = 0.08;
2970      padMIP[i-1]->SetBottomMargin(mBot);
2971      padMIP[i-1]->SetTopMargin(mTop);
2972      padMIP[i-1]->SetLeftMargin(mLeft);
2973      padMIP[i-1]->SetRightMargin(mRight);
2974
2975      cMIP->cd();
2976
2977      padMIP[i-1]->Draw();
2978
2979   }
2980
2981   latex->SetTextAngle(0);
2982   latex->SetTextSize(0.085);
2983
2984   Int_t index_dec = 7;
2985
2986   for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
2987     padMIP[index_eta]->cd();
2988     padMIP[index_eta]->SetTickx(1);
2989     padMIP[index_eta]->SetTicky(1);
2990     padMIP[index_eta]->SetLogy(0);
2991     padMIP[index_eta]->SetLogx(0);
2992     hMIP[index_eta]->GetXaxis()->SetLabelSize(0.06);
2993     hMIP[index_eta]->GetYaxis()->SetLabelSize(0.06);
2994     hMIP[index_eta]->GetXaxis()->SetTitleSize(0);
2995     hMIP[index_eta]->GetYaxis()->SetTitleSize(0);
2996
2997     if(index_eta<4){
2998       hMIP[index_eta]->Draw("colz");
2999       pMIP[index_eta]->Draw("samep");
3000       latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3001     }
3002     else{
3003       hMIP[index_dec]->Draw("colz");
3004       pMIP[index_dec]->Draw("samep");
3005       latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
3006       index_dec--;
3007     }
3008
3009   }
3010
3011
3012
3013   cMIP->SaveAs("MIPvsPhi.png");
3014
3015   //Plateau vs phi, different eta intervals
3016   TCanvas* cPlateau = new TCanvas("cPlateau2", "Raa Pions", 900, 500);
3017   TPad* padPlateau[nPadX*nPadY];
3018
3019   const char* yTitlePlateau = "d#it{E}/d#it{x}_{Plateau}";
3020   const char* xTitlePlateau = "#phi (rad)";
3021   cPlateau->cd();
3022
3023
3024   //TLatex* latex = new TLatex();
3025   latex->SetNDC();
3026   latex->SetTextAlign(22);
3027   latex->SetTextAngle(0);
3028   latex->SetTextFont(62);
3029
3030   latex->DrawLatex(0.53,0.04,xTitlePlateau);
3031   latex->SetTextAngle(90);
3032   latex->SetTextSize(0.055);
3033   latex->DrawLatex(0.025,0.5,yTitlePlateau);
3034   latex->SetTextSize(0.05);
3035
3036
3037   for (Int_t i = 1; i < 9; i++) {
3038
3039     Int_t iX = (i-1)%nPadX;
3040     Int_t iY = (i-1)/nPadX;
3041     Float_t x1 = leftMargin + iX*width;
3042     if(iX==2)
3043       x1-=0.01;
3044     else if(iX==1)
3045       x1+=0.01;
3046     Float_t x2 = leftMargin + (iX + 1)*width;
3047      if(iX==0)
3048        x2+=0.01;
3049      else if(iX==1)
3050        x2-=0.01;
3051      Float_t y1 = 1 - topMargin - (iY +1)*height;
3052      Float_t y2 = 1 - topMargin - iY*height;
3053      padPlateau[i-1] = new TPad(Form("padPlateau%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3054      Float_t mBot = noMargin;
3055      Float_t mTop = noMargin;
3056      Float_t mLeft = noMargin;
3057      Float_t mRight = noMargin;
3058      if(iY==0)       mTop   = shift;
3059      if(iY==nPadY-1) mBot   = shift;
3060      if(iX==0)       mLeft  = 0.08;
3061      if(iX==nPadX-1) mRight = 0.08;
3062      padPlateau[i-1]->SetBottomMargin(mBot);
3063      padPlateau[i-1]->SetTopMargin(mTop);
3064      padPlateau[i-1]->SetLeftMargin(mLeft);
3065      padPlateau[i-1]->SetRightMargin(mRight);
3066
3067      cPlateau->cd();
3068
3069      padPlateau[i-1]->Draw();
3070
3071   }
3072
3073   latex->SetTextAngle(0);
3074   latex->SetTextSize(0.085);
3075
3076   index_dec = 7;
3077
3078   for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
3079     padPlateau[index_eta]->cd();
3080     padPlateau[index_eta]->SetTickx(1);
3081     padPlateau[index_eta]->SetTicky(1);
3082     padPlateau[index_eta]->SetLogy(0);
3083     padPlateau[index_eta]->SetLogx(0);
3084     hPlateau[index_eta]->GetXaxis()->SetLabelSize(0.06);
3085
3086     hPlateau[index_eta]->GetYaxis()->SetLabelSize(0.06);
3087     hPlateau[index_eta]->GetXaxis()->SetTitleSize(0);
3088     hPlateau[index_eta]->GetYaxis()->SetTitleSize(0);
3089
3090     if(index_eta<4){
3091       hPlateau[index_eta]->Draw("colz");
3092       pPlateau[index_eta]->Draw("samep");
3093       latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3094     }
3095     else{
3096       hPlateau[index_dec]->Draw("colz");
3097       pPlateau[index_dec]->Draw("samep");
3098       latex->DrawLatex(0.5,0.1,LatexEta[index_dec]);
3099       index_dec--;
3100     }
3101
3102   }
3103
3104
3105
3106   cPlateau->SaveAs("PlateauVsPhi.png");
3107
3108
3109   //MIP vs Nch
3110   TCanvas* cMIPVsNch = new TCanvas("cMIPVsNch2", "Raa Pions", 900, 500);
3111   TPad* padMIPVsNch[nPadX*nPadY];
3112
3113   const char* yTitleMIPVsNch = "d#it{E}/d#it{x}_{MIP}";
3114   const char* xTitleMIPVsNch = "TPC-track multiplicity";
3115   cMIPVsNch->cd();
3116
3117
3118   latex->SetNDC();
3119   latex->SetTextAlign(22);
3120   latex->SetTextAngle(0);
3121   latex->SetTextFont(62);
3122
3123   latex->DrawLatex(0.53,0.04,xTitleMIPVsNch);
3124   latex->SetTextAngle(90);
3125   latex->SetTextSize(0.055);
3126   latex->DrawLatex(0.025,0.5,yTitleMIPVsNch);
3127   latex->SetTextSize(0.05);
3128
3129
3130   for (Int_t i = 1; i < 9; i++) {
3131
3132     Int_t iX = (i-1)%nPadX;
3133     Int_t iY = (i-1)/nPadX;
3134     Float_t x1 = leftMargin + iX*width;
3135     if(iX==2)
3136       x1-=0.01;
3137     else if(iX==1)
3138       x1+=0.01;
3139     Float_t x2 = leftMargin + (iX + 1)*width;
3140      if(iX==0)
3141        x2+=0.01;
3142      else if(iX==1)
3143        x2-=0.01;
3144      Float_t y1 = 1 - topMargin - (iY +1)*height;
3145      Float_t y2 = 1 - topMargin - iY*height;
3146      padMIPVsNch[i-1] = new TPad(Form("padMIPVsNch%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3147      Float_t mBot = noMargin;
3148      Float_t mTop = noMargin;
3149      Float_t mLeft = noMargin;
3150      Float_t mRight = noMargin;
3151      if(iY==0)       mTop   = shift;
3152      if(iY==nPadY-1) mBot   = shift;
3153      if(iX==0)       mLeft  = 0.08;
3154      if(iX==nPadX-1) mRight = 0.08;
3155      padMIPVsNch[i-1]->SetBottomMargin(mBot);
3156      padMIPVsNch[i-1]->SetTopMargin(mTop);
3157      padMIPVsNch[i-1]->SetLeftMargin(mLeft);
3158      padMIPVsNch[i-1]->SetRightMargin(mRight);
3159
3160      if(i>=5&&i<9)
3161        padMIPVsNch[i-1]->SetBottomMargin(mBot+0.05);
3162
3163      cMIPVsNch->cd();
3164
3165      padMIPVsNch[i-1]->Draw();
3166
3167   }
3168
3169   latex->SetTextAngle(0);
3170   latex->SetTextSize(0.085);
3171
3172   index_dec = 7;
3173
3174   for(Int_t index_eta = 0; index_eta < 8; ++index_eta){
3175     padMIPVsNch[index_eta]->cd();
3176
3177     padMIPVsNch[index_eta]->SetTickx(1);
3178     padMIPVsNch[index_eta]->SetTicky(1);
3179     padMIPVsNch[index_eta]->SetLogy(0);
3180     padMIPVsNch[index_eta]->SetLogx(1);
3181     hMIPVsNch[index_eta]->GetXaxis()->SetLabelSize(0.08);
3182
3183     hMIPVsNch[index_eta]->GetXaxis()->SetRangeUser(10,2000);
3184     hMIPVsNch[index_eta]->GetYaxis()->SetRangeUser(43.5,54.5);
3185     hMIPVsNch[index_eta]->GetYaxis()->SetLabelSize(0.06);
3186     hMIPVsNch[index_eta]->GetXaxis()->SetTitleSize(0);
3187     hMIPVsNch[index_eta]->GetYaxis()->SetTitleSize(0);
3188
3189     if(index_eta<4){
3190       hMIPVsNch[index_eta]->Draw("colz");
3191       pMIPVsNch[index_eta]->Draw("samep");
3192       latex->DrawLatex(0.5,0.1,LatexEta[index_eta]);
3193     }
3194     else{
3195       hMIPVsNch[index_dec]->Draw("colz");
3196       pMIPVsNch[index_dec]->Draw("samep");
3197       latex->DrawLatex(0.5,0.2,LatexEta[index_dec]);
3198       index_dec--;
3199     }
3200
3201   }
3202
3203
3204
3205
3206
3207
3208
3209
3210   cMIPVsNch->SaveAs("MIPvsNch.png");
3211
3212
3213
3214
3215
3216
3217
3218
3219
3220   //eta dependence, scaling
3221   const Int_t nPadXS = 3;
3222   const Int_t nPadYS = 1;
3223   Float_t noMarginS    = 0.0;
3224   Float_t topMarginS   = 0.01;
3225   Float_t botMarginS   = 0.1;
3226   Float_t leftMarginS  = 0.05;
3227   Float_t rightMarginS = 0.01;
3228   Float_t widthS       = (1-rightMarginS-leftMarginS)/nPadXS;
3229   Float_t heightS      = (1-botMarginS-topMarginS)/nPadYS;
3230   Float_t shiftS = 0.05;
3231
3232
3233   TCanvas* cMIPEta = new TCanvas("cMIPEta", "MIP eta", 900, 500);
3234   TPad* padMIPEta[nPadXS*nPadYS];
3235   //Float_t shift = 0.1;
3236   //Float_t shift = 0.05;
3237   const char* yTitleMIPS = "d#it{E}/d#it{x}";
3238   const char* xTitleMIPS = "#eta";
3239   cMIPEta->cd();
3240
3241
3242   //TLatexlatex = new TLatex();
3243   latex->SetNDC();
3244   latex->SetTextAlign(22);
3245   latex->SetTextAngle(0);
3246   latex->SetTextFont(42);
3247
3248   latex->DrawLatex(0.53,0.04,xTitleMIPS);
3249   latex->SetTextAngle(90);
3250   latex->SetTextSize(0.055);
3251   latex->DrawLatex(0.025,0.5,yTitleMIPS);
3252   latex->SetTextSize(0.05);
3253
3254
3255   for (Int_t i = 1; i < 4; i++) {
3256
3257     Int_t iX = (i-1)%nPadXS;
3258     Int_t iY = (i-1)/nPadXS;
3259     Float_t x1 = leftMarginS + iX*widthS;
3260     if(iX==2)
3261       x1-=0.01;
3262     else if(iX==1)
3263       x1+=0.01;
3264     Float_t x2 = leftMarginS + (iX + 1)*widthS;
3265      if(iX==0)
3266        x2+=0.01;
3267      else if(iX==1)
3268        x2-=0.01;
3269      Float_t y1 = 1 - topMarginS - (iY +1)*heightS;
3270      Float_t y2 = 1 - topMarginS - iY*heightS;
3271      padMIPEta[i-1] = new TPad(Form("padMIPEta%d",i),"", x1, y1, x2, y2, 0, 0, 0);
3272      Float_t mBot = noMarginS;
3273      Float_t mTop = noMarginS;
3274      Float_t mLeft = noMarginS;
3275      Float_t mRight = noMarginS;
3276      if(iY==0)       mTop   = shiftS;
3277      if(iY==nPadYS-1) mBot   = shiftS;
3278      if(iX==0)       mLeft  = 0.08;
3279      if(iX==nPadXS-1) mRight = 0.08;
3280      padMIPEta[i-1]->SetBottomMargin(mBot);
3281      padMIPEta[i-1]->SetTopMargin(mTop);
3282      padMIPEta[i-1]->SetLeftMargin(mLeft);
3283      padMIPEta[i-1]->SetRightMargin(mRight);
3284
3285      cMIPEta->cd();
3286
3287      padMIPEta[i-1]->Draw();
3288
3289   }
3290
3291   latex->SetTextAngle(0);
3292   latex->SetTextSize(0.085);
3293
3294
3295   TF1 *fetapos = new TF1("fetapos","pol3",0.0,0.8);
3296   TF1 *fetaneg = new TF1("fetaneg","pol3",-0.8,0.0);
3297
3298   pMIPVsEta->Fit(fetaneg,"0","",-0.8,0.0);
3299   pMIPVsEta->Fit(fetapos,"0","",0.0,0.8);
3300   TLegend *leg1=new TLegend(0.15,0.5,0.4,0.7);
3301
3302   fetapos->SetLineColor(1);
3303   fetapos->SetLineStyle(1);
3304   fetapos->SetLineWidth(2);
3305
3306   fetaneg->SetLineColor(1);
3307   fetaneg->SetLineStyle(1);
3308   fetaneg->SetLineWidth(2);
3309
3310   TH1D *hframe20 = new TH1D("hframe20","hframe20",8,-0.8,0.8);
3311
3312   TH1D *hMIPPiPrim = new TH1D("hMIPPiPrim","",16,-0.8,0.8);
3313   for(Int_t bin =1; bin <= hMIPPiPrim->GetNbinsX(); ++bin){
3314     Double_t eta = pMIPVsEta->GetBinCenter(bin);
3315     Double_t dedx = pMIPVsEta->GetBinContent(bin);
3316     Double_t e_dedx = pMIPVsEta->GetBinError(bin);
3317     if(eta<0){
3318       dedx *= 50/fetaneg->Eval(eta);
3319       e_dedx *= 50/fetaneg->Eval(eta);
3320     }else{
3321       dedx *= 50/fetapos->Eval(eta);
3322       e_dedx *= 50/fetapos->Eval(eta);
3323     }
3324     hMIPPiPrim->SetBinContent(bin,dedx);
3325     hMIPPiPrim->SetBinError(bin,e_dedx);
3326
3327   }
3328
3329   TH1D *hPlateauPiPrim = new TH1D("hPlateauPiPrim","",16,-0.8,0.8);
3330   for(Int_t bin =1; bin <= hPlateauPiPrim->GetNbinsX(); ++bin){
3331     Double_t eta = pPlateauVsEta->GetBinCenter(bin);
3332     Double_t dedx = pPlateauVsEta->GetBinContent(bin);
3333     Double_t e_dedx = pPlateauVsEta->GetBinError(bin);
3334     if(eta<0){
3335       dedx *= 50/fetaneg->Eval(eta);
3336       e_dedx *= 50/fetaneg->Eval(eta);
3337     }else{
3338       dedx *= 50/fetapos->Eval(eta);
3339       e_dedx *= 50/fetapos->Eval(eta);
3340     }
3341     hPlateauPiPrim->SetBinContent(bin,dedx);
3342     hPlateauPiPrim->SetBinError(bin,e_dedx);
3343
3344   }
3345
3346
3347   TH1D *h2GeVPiSec = new TH1D("h2GeVPiSec","",8,-0.8,0.8);
3348   for(Int_t bin =1; bin <= h2GeVPiSec->GetNbinsX(); ++bin){
3349     Double_t eta = hPion2GeV->GetBinCenter(bin);
3350     Double_t dedx = hPion2GeV->GetBinContent(bin);
3351     Double_t e_dedx = hPion2GeV->GetBinError(bin);
3352     if(eta<0){
3353       dedx *= 50/fetaneg->Eval(eta);
3354       e_dedx *= 50/fetaneg->Eval(eta);
3355     }else{
3356       dedx *= 50/fetapos->Eval(eta);
3357       e_dedx *= 50/fetapos->Eval(eta);
3358     }
3359     h2GeVPiSec->SetBinContent(bin,dedx);
3360     h2GeVPiSec->SetBinError(bin,e_dedx);
3361
3362   }
3363
3364
3365
3366
3367
3368
3369
3370
3371
3372
3373   for(Int_t index_eta = 0; index_eta < 3; ++index_eta){
3374     padMIPEta[index_eta]->cd();
3375     padMIPEta[index_eta]->SetTickx(1);
3376     padMIPEta[index_eta]->SetTicky(1);
3377     padMIPEta[index_eta]->SetLogy(0);
3378     padMIPEta[index_eta]->SetLogx(0);
3379
3380     if(index_eta==0){
3381       hframe20->GetYaxis()->SetRangeUser(30,95);
3382       hframe20->GetXaxis()->SetLabelSize(0.06);
3383       hframe20->GetXaxis()->SetLabelOffset(0.004);
3384       hframe20->GetYaxis()->SetLabelSize(0.06);
3385       hframe20->GetXaxis()->SetTitleSize(0);
3386       hframe20->GetYaxis()->SetTitleSize(0);
3387       hframe20->Draw();
3388       hMIPVsEta->Draw("colz same");
3389       hPlateauVsEta->Draw("colz same");
3390
3391       latex->SetTextSize(0.05);
3392       latex->DrawLatex(0.55,0.45,"MIP Pions, 0.4<#it{p}<0.6 GeV/c");
3393       latex->DrawLatex(0.55,0.85,"Electrons, 0.4<#it{p}<0.6 GeV/c");
3394
3395       }
3396
3397
3398     if(index_eta==1){
3399       pMIPVsEtaV0s->GetXaxis()->SetLabelSize(0.06);
3400       pMIPVsEtaV0s->GetYaxis()->SetLabelSize(0.06);
3401       pMIPVsEtaV0s->GetXaxis()->SetLabelOffset(0.004);
3402       pMIPVsEtaV0s->GetXaxis()->SetTitleSize(0);
3403       pMIPVsEtaV0s->GetYaxis()->SetTitleSize(0);
3404       pMIPVsEtaV0s->GetYaxis()->SetRangeUser(30,95);
3405       pMIPVsEtaV0s->SetMarkerStyle(29);
3406       pMIPVsEtaV0s->SetMarkerColor(4);
3407       pMIPVsEtaV0s->SetLineColor(4);
3408       pMIPVsEtaV0s->Draw();
3409
3410       pMIPVsEta->SetMarkerStyle(25);
3411       pMIPVsEta->SetMarkerColor(1);
3412       pMIPVsEta->SetLineColor(1);
3413       pMIPVsEta->Draw("samep");
3414       pPlateauVsEta->SetMarkerStyle(21);
3415       pPlateauVsEta->SetMarkerColor(kBlue);
3416       pPlateauVsEta->SetLineColor(4);
3417       pPlateauVsEta->Draw("samep");
3418
3419       //fetapos->Draw("same");
3420       //fetaneg->Draw("same");
3421
3422
3423       pMIPVsEta->Draw("samep");
3424       hPion2GeV->SetMarkerStyle(20);
3425       hPion2GeV->SetMarkerColor(4);
3426       hPion2GeV->SetLineColor(4);
3427       hPion2GeV->Draw("samep");
3428       pPlateauVsEta->Draw("samep");
3429
3430
3431       leg1->SetTextFont(42);
3432       leg1->SetTextSize(0.045);
3433       leg1->SetLineColor(kWhite);
3434       leg1->SetLineStyle(3);
3435       leg1->SetShadowColor(kWhite);
3436       leg1->SetFillColor(kWhite);
3437       leg1->SetFillStyle(0);
3438       leg1->AddEntry(pMIPVsEta,"Primary pions at MIP","P");
3439       leg1->AddEntry(pMIPVsEtaV0s,"Secondary pions at MIP","P");
3440       leg1->AddEntry(hPion2GeV,"Secondary pions at #beta#gamma ~ 14","P");
3441       leg1->AddEntry(pPlateauVsEta,"Electrons,  at #beta#gamma ~ 800","P");
3442       leg1->Draw();
3443
3444
3445
3446     }
3447
3448
3449
3450     if(index_eta==2){
3451       hMIPPiPrim->GetXaxis()->SetLabelSize(0.06);
3452       hMIPPiPrim->GetYaxis()->SetLabelSize(0.06);
3453       hMIPPiPrim->GetXaxis()->SetLabelOffset(0.004);
3454       hMIPPiPrim->GetXaxis()->SetTitleSize(0);
3455       hMIPPiPrim->GetYaxis()->SetTitleSize(0);
3456       hMIPPiPrim->GetYaxis()->SetRangeUser(30,95);
3457       hMIPPiPrim->SetMarkerStyle(29);
3458       hMIPPiPrim->SetMarkerColor(4);
3459       hMIPPiPrim->SetLineColor(4);
3460       hMIPPiPrim->Draw();
3461
3462
3463       hPlateauPiPrim->SetMarkerStyle(21);
3464       hPlateauPiPrim->SetLineColor(4);
3465       hPlateauPiPrim->SetMarkerColor(4);
3466       hPlateauPiPrim->Draw("samep");
3467       h2GeVPiSec->SetMarkerStyle(20);
3468       h2GeVPiSec->SetLineColor(4);
3469       h2GeVPiSec->SetMarkerColor(4);
3470       h2GeVPiSec->Draw("samep");
3471
3472       latex->DrawLatex(0.45,0.85,"#LT d#it{E}/d#it{x} #GT #times 50 / #LT d#it{E}/d#it{x} #GT_{primary #pi at MIP}");
3473
3474     }
3475
3476
3477
3478
3479
3480
3481
3482
3483
3484
3485
3486   }
3487
3488
3489   cMIPEta->SaveAs("MIPvsEta.png");
3490
3491
3492
3493
3494 }
3495 TH2D * AddTwoSameBinningTH2D(TH2D *hPos, TH2D *hNeg, const Char_t *nameHist){
3496
3497   TH2D *hout = (TH2D *)hPos->Clone(nameHist);
3498   hout->Reset();
3499
3500
3501   for(Int_t binx = 1; binx <= hPos->GetNbinsX(); ++binx){
3502
3503     for(Int_t biny = 1; biny <= hPos->GetNbinsY(); ++biny){
3504
3505       Double_t yield1 =  hPos->GetBinContent(binx,biny);
3506       Double_t yield2 =  hNeg->GetBinContent(binx,biny);
3507       Double_t e_yield1 =  hPos->GetBinError(binx,biny);
3508       Double_t e_yield2 =  hNeg->GetBinError(binx,biny);
3509
3510       Double_t yield = yield1 + yield2;
3511       Double_t e_yield = TMath::Sqrt( e_yield1*e_yield1 + e_yield2*e_yield2 );
3512
3513       hout->SetBinContent(binx,biny,yield);
3514       hout->SetBinError(binx,biny,e_yield);
3515
3516     }
3517
3518   }
3519
3520   hout->SetEntries( hPos->GetEntries()+hNeg->GetEntries() );
3521
3522   return hout;
3523
3524 }