PWG2/SPECTRA/Fit/CombineSpectra.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <iostream>
   3 #include "TH1F.h"
   4 #include "TGraphErrors.h"
   5 #include "AliBWFunc.h"
   6 #include "AliBWTools.h"
   7 #include "TF1.h"
   8 #include "TFile.h"
   9 #include "TDatabasePDG.h"
  10 #include "TROOT.h"
  11 #include "TCanvas.h"
  12 #include "TFolder.h"
  13 #include "TStyle.h"
  14 #include "AliLatexTable.h"
  15 #include "TLegend.h"
  16 #include "TVirtualFitter.h"
  17 #include "TMath.h"
  18 #include "TH2F.h"
  19 #include "TSystem.h"
  20 #include "TLine.h"
  21 #include "TLatex.h"
  22 #include "TMath.h"
  23
  24 #include "TASImage.h"
  25 #include "TPaveText.h"
  26 #include "StarPPSpectra.C"
  27 #include "GetE735Ratios.C"
  28 #include "TString.h"
  29 #include "TObjString.h"
  30 #include "TGraphAsymmErrors.h"
  31
  32 #endif
  33
  34
  35 using namespace std;
  36
  37
  38 // A bunch of useful enums and constants
  39 enum {kPion=0,kKaon,kProton,kNPart};
  40 enum {kTPC=0,kTOF,kITS,kITSTPC,kK0,kKinks,kCombTOFTPC,kCombAll,kNHist};// "k0" listed here as a kind of PID method...
  41 const Int_t kNDet = kITS+2;
  42 enum {kPos=0,kNeg,kNCharge};
  43 enum {kPhojet=0,kPyTuneAtlasCSC, kPyTuneCMS6D6T, kPyTunePerugia0, kNTunes} ;
  44 enum {kFitLevi=0, kFitUA1, kFitPowerLaw,
  45       kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0,
  46       kNFit};
  47 enum {kDoFits=0, kDoRatios, kDoSuperposition, kDoDrawWithModels, kDoCompareToStar, kDoDrawSyst, kDoCompareToAllCharged, kDoHelp};
  48 enum {kStatErrors = 0, kSystErrors, kStatSystErrors}; // which errors do we put in the histo that we fit? stat,syst or stat+syst?
  49
  50 // flags, labels and names
  51 const char * partFlag[] = {"Pion", "Kaon", "Proton"};
  52 const char * detFlag[]  = {"TPC", "TOF", "ITS standalone", "ITS-TPC", "K0", "Kinks", "Combined TOF + TPC", "Combined TOF + TPC + ITS"};
  53 const char * chargeFlag[]  = {"Pos", "Neg"};
  54 const char * chargeLabel[]  = {"Positive", "Negative"};
  55 const char * partLabel[kNPart][kNCharge] = {{"#pi^{+}", "#pi^{-}"},
  56                                             //                                      {"K^{+} (#times 2)", "K^{-} (#times 2)"},
  57                                             {"K^{+}", "K^{-}"},
  58                                             {"p" ,  "#bar{p}"}};
  59 const char * partLatex[kNPart][kNCharge] = {{"$\\pi^{+}$", "$\\pi^{-}$"},
  60                                             //                                      {"K^{+} (#times 2)", "K^{-} (#times 2)"},
  61                                             {"$K^{+}$", "$K^{-}$"},
  62                                             {"$p$" ,  "$\\bar{p}$"}};
  63 const char * mcTuneName[] = {"Phojet",
  64                              "Pythia - CSC 306",
  65                              "Pythia - D6T 109",
  66                              "Pythia - Perugia0 - 320", };
  67 const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
  68
  69 // Style
  70 //const Int_t marker[] = {25,24,28,20,21}; // set for highlithining marek
  71 //const Int_t marker[] = {20,24,25,28,21}; // standard set
  72 const Int_t marker[] = {24,25,28,27,21}; // No full symbols
  73 const Int_t color [] = {1,2,4};
  74
  75 const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
  76 const Int_t mcLineStyle[] = {1,2,3,4};
  77
  78
  79 // Template needed to combine different detectors
  80 const Float_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,2,2.2,2.4,2.6};
  81 Int_t nbinsTempl=34;
  82
  83 TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
  84
  85 //  Globals
  86 TH1F * hSpectra[kNHist][kNPart][kNCharge];
  87 TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
  88 TH1F * hSystError[kNHist][kNPart][kNCharge];
  89 Double_t mass[kNPart];
  90
  91 // Functions:
  92 // Loading
  93 void LoadSpectra() ;
  94 void LoadMC() ;
  95
  96 // Additional tasks (may be uncommented below)
  97 void LoadLibs();
  98 void DrawStar(Int_t icharge);
  99 void GetITSResiduals();
 100 void DrawWithModels() ;
 101 void DrawAllAndKaons();
 102 void DrawWithJacek();
 103 void DrawRatioToStar();
 104 void DrawRatios();
 105 void DrawWithSyst();
 106 void FitCombined();
 107 void PrintCanvas(TCanvas* c,const TString formats) ;
 108 void RebinHistosForRatios() ;
 109 void Help();
 110
 111 // External stuff
 112 void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance"){
 113
 114   TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.72,0.89,0.89);
 115   //TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.62,0.85,0.75);
 116   myPadLogo->SetFillColor(0);
 117   myPadLogo->SetBorderMode(0);
 118   myPadLogo->SetBorderSize(2);
 119   myPadLogo->SetFrameBorderMode(0);
 120   myPadLogo->SetLeftMargin(0.0);
 121   myPadLogo->SetTopMargin(0.0);
 122   myPadLogo->SetBottomMargin(0.0);
 123   myPadLogo->SetRightMargin(0.0);
 124   myPadLogo->SetFillStyle(0);
 125   myPadLogo->Draw();
 126   myPadLogo->cd();
 127   TASImage *myAliceLogo = new TASImage("alice_logo.png");
 128   myAliceLogo->Draw();
 129   c->cd();
 130   TPaveText* t1=new TPaveText(0.418103, 0.837798, 0.656609, 0.888393,"NDC");
 131   t1->SetFillStyle(0);
 132   t1->SetBorderSize(0);
 133   t1->AddText(0.,0.,label);
 134   t1->SetTextColor(kRed);
 135   t1->SetTextFont(42);
 136   t1->SetTextSize(0.04);
 137   t1->Draw();
 138   TPaveText* t2=new TPaveText(0.418103, 0.80, 0.656609, 0.84,"NDC");
 139   t2->SetFillStyle(0);
 140   t2->SetBorderSize(0);
 141   t2->AddText(0.,0.,"pp at #sqrt{s} = 900 GeV (2009 data)");
 142   t2->SetTextColor(kRed);
 143   t2->SetTextFont(42);
 144   t2->SetTextSize(0.027);
 145   t2->Draw();
 146   TPaveText* t3=new TPaveText(0.418103, 0.76, 0.656609, 0.80,"NDC");
 147   t3->SetFillStyle(0);
 148   t3->SetBorderSize(0);
 149   t3->AddText(0.,0.,"Statistical and systematic errors");
 150   t3->SetTextColor(kRed);
 151   t3->SetTextFont(42);
 152   t3->SetTextSize(0.027);
 153   t3->Draw();
 154   // TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
 155   // t2->SetFillStyle(0);
 156   // t2->SetBorderSize(0);
 157   // t2->SetTextColor(kRed);
 158   // t2->SetTextFont(52);
 159   // t2->AddText(0.,0.,today.Data());
 160   // t2->Draw();
 161 }
 162
 163 // Used to tag plots
 164 TDatime dt;
 165 TString today = "";
 166
 167
 168
 169 // Switches
 170 Bool_t convertToMT = 0;
 171 Bool_t sumCharge = 0;
 172 Int_t whatToFit = kStatErrors;
 173 Bool_t scaleKaons =  kFALSE;
 174 Bool_t drawStar =  kFALSE; // Overlay star when doing fits
 175 Bool_t correctSecondaries  = 1;
 176 Bool_t correctGeantFlukaXS = 1;
 177 Int_t iCombInStudy = kCombAll; //kCombTOFTPC
 178 Int_t fitFuncID = kFitLevi;
 179 Bool_t showMC=kTRUE;
 180 Bool_t showE735=kTRUE;
 181 Bool_t useSecCorrFromDCA=kTRUE;
 182 Bool_t flagPreliminary=kFALSE; // Add "preliminary" flag and logo to the plots
 183 Bool_t doPrint = 0;
 184 TString printFormats = "C,eps,root"; // format in which canvases will be printed, if PrintCanvas is called (not all prints are based on printcanvas at the moment). This is a comma separated list.
 185
 186
 187 void CombineSpectra(Int_t analysisType=kDoHelp, Int_t  locfitFuncID = kFitLevi) {  //kDoSuperposition;//kDoDrawWithModels;// kDoFits; //kDoRatios;
 188
 189   // This macro is used to combine the 900 GeV spectra from 2009
 190
 191   fitFuncID=locfitFuncID;
 192
 193   // Load stuff
 194   LoadLibs();
 195
 196   // Print Help and quit
 197   if (analysisType == kDoHelp) {
 198     Help();
 199     return;
 200   }
 201
 202
 203   // Set date
 204   today = today +  long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
 205
 206
 207   // Set Masses
 208   mass[kPion]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
 209   mass[kKaon]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
 210   mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
 211
 212   // Load histos
 213   LoadSpectra();
 214   LoadMC();
 215
 216   // Additional tasks
 217   //DrawStar(icharge);
 218   //  GetITSResiduals();
 219   if(analysisType==kDoSuperposition) DrawAllAndKaons();
 220   else if(analysisType==kDoDrawWithModels)  DrawWithModels() ;
 221   else if(analysisType==kDoCompareToAllCharged) DrawWithJacek();
 222   else if(analysisType==kDoCompareToStar) DrawRatioToStar();
 223   else if(analysisType==kDoRatios) DrawRatios();
 224   else if(analysisType==kDoDrawSyst) DrawWithSyst();
 225   else if(analysisType==kDoFits) FitCombined();
 226   return;
 227 }
 228
 229
 230 void FitCombined() {
 231   gStyle->SetOptTitle(0);
 232   gStyle->SetOptStat(0);
 233
 234   // Draw combined & Fit
 235   AliBWFunc * fm = new AliBWFunc;
 236   fm->SetVarType(AliBWFunc::kdNdpt);
 237   if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
 238
 239   // table to print results
 240   AliLatexTable table(10,"c|ccccccccc");
 241   if (fitFuncID == kFitLevi) {
 242     table.InsertCustomRow("Part & Yield & Yield (FIT) &  T Slope & n & $\\chi^2$/NDF & Min X & Extrap (%) & $\\langle p_{t} \\rangle (Fit)$  & $\\langle p_{t} \\rangle$ (FIT+DATA)\\\\");
 243   }  else if (fitFuncID == kFitPowerLaw) {
 244     table.InsertCustomRow("Part & Yield & Norm &  n & pt0 & $\\chi^2$/NDF & Min X & Frac Above  & $\\langle p_{t} \\rangle (Fit)$  & $\\langle p_{t} \\rangle$ (FIT+DATA)  \\\\");
 245   } else {
 246     table.InsertCustomRow("Part & Yield & Par0 & Par2  & Par1 & $\\chi^2$/NDF & Min X & Frac Above  & $\\langle p_{t} \\rangle (Fit)$  & $\\langle p_{t} \\rangle$ (FIT+DATA)  \\\\");
 247
 248   }
 249   table.InsertHline();
 250   AliLatexTable tempTable(3,"c|cc");
 251   tempTable.InsertCustomRow("Part & Yield &  $\\langle p_{t} \\rangle$ \\\\");
 252   tempTable.InsertHline();
 253
 254   TH1F* hRatiosToFit[kNPart][kNCharge]; // Ratio data/fit function, stat error by default
 255   TH1F* hRatiosToFitSyst[kNPart][kNCharge]; // Ratio data/fit, stat + syst
 256   //  Fit all
 257   Int_t chargeLoop = sumCharge ? 1 : 2;
 258   for(Int_t icharge = 0; icharge < chargeLoop; icharge++){
 259
 260     TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge],700,700);
 261     c2->SetTickx();
 262     c2->SetTicky();
 263     c2->SetLeftMargin(0.14);
 264     TCanvas * c2r = new TCanvas(TString("cCombinedRatio")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombinedRatio")+chargeFlag[icharge],700,700);
 265     c2->cd();
 266     gPad->SetLogy();
 267     TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
 268     hempty->SetXTitle("p_{t} (GeV/c)");
 269     hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
 270     hempty->GetYaxis()->SetTitleOffset(1.35);
 271     hempty->GetXaxis()->SetTitleOffset(1.1);
 272     hempty->Draw();
 273     c2r->cd();
 274     gPad->SetGridy();
 275     TH2F * hemptyR = new TH2F(TString("hemptyR")+long(icharge),"hemptyR",100,0.,2.9, 100, 0.5,1.5);
 276     hemptyR->SetXTitle("p_{t} (GeV/c)");
 277     hemptyR->SetYTitle("Data/Fit");
 278     hemptyR->Draw();
 279
 280
 281     TLegend * l = new TLegend(0.176724, 0.153274, 0.475575, 0.434524,chargeLabel[icharge]);
 282     l->SetFillColor(kWhite);
 283     l->SetTextSize(0.035);
 284
 285
 286
 287     TPaveText* tf=new TPaveText(0.321839, 0.175595, 0.685345, 0.299107,"NDC");
 288     if(fitFuncID == kFitLevi){
 289       tf->AddText("#frac{dN}{dp_{t}} #propto p_{t} #left(1+#frac{#sqrt{m^{2}+p_{t}^{2}} -m}{nT} #right)^{-n}");
 290       //      tf->SetNDC();
 291       tf->SetTextFont(12);
 292       tf->SetTextSize(0.032);
 293     }
 294     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 295       printf(" ----- Fit %s %s ------\n",partFlag[ipart],chargeFlag[icharge]);
 296       Float_t fitmin = 0;
 297       Float_t fitmax = 3;
 298
 299       // Get functions
 300       TF1 * func = 0;
 301       Int_t normPar = -1;
 302       if(fitFuncID == kFitLevi)          {
 303         normPar = 0; // The levi is normalized by parameter 0
 304         if (ipart == kPion)
 305           func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
 306         if (ipart == kKaon)
 307           func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
 308         if (ipart == kProton)
 309           func = fm->GetLevi(mass[ipart], 0.15, 8.5, 0.09);
 310       }
 311       else if(fitFuncID == kFitUA1)      func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
 312       else if(fitFuncID == kFitPowerLaw) {
 313         // fitmin=0;
 314         // fitmax=1;
 315         if (ipart == kPion)
 316           func = fm->GetPowerLaw(1.0, 8.6, 7);
 317         if (ipart == kKaon)
 318           func = fm->GetPowerLaw(3.0, 12, 2.6);
 319         if (ipart == kProton)
 320           func = fm->GetPowerLaw(24, 72, 0.8);
 321       }
 322       else if(fitFuncID == kFitPhojet)   func = fm->GetHistoFunc(hSpectraMC[kPhojet]        [ipart][icharge]);
 323       else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
 324       else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
 325       else if(fitFuncID == kFitCMS6D6T)  func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
 326       else {
 327         cout << "Unknown fit ID " << fitFuncID << endl;
 328         return;
 329       }
 330
 331       if(fitFuncID >= kFitPhojet){
 332         fitmin = 0.0;
 333         fitmax = 1.0;
 334       }
 335
 336       // Temp: fit histo with sist errors
 337       TH1F * hsyst     = new TH1F(*htemplate);
 338       TH1F * hsyststat = 0;
 339       hsyst->SetFillColor(kGray);
 340       AliBWTools::GetValueAndError(hsyst,hSpectra[iCombInStudy][ipart][icharge],hSystError[iCombInStudy][ipart][icharge],kTRUE);
 341       hsyststat= new TH1F(*hsyst);
 342       AliBWTools::GetHistoCombinedErrors(hsyststat,hSpectra[iCombInStudy][ipart][icharge]); // combine syst and stat
 343       hsyststat->SetFillColor(kGray);
 344
 345       TH1F * hToFit = 0;
 346       if (whatToFit == kStatErrors) hToFit = hSpectra[iCombInStudy][ipart][icharge]; // Shorthand
 347       if (sumCharge) hToFit->Add(hSpectra[iCombInStudy][ipart][1]);
 348       if (whatToFit == kStatSystErrors) {
 349         hToFit = hsyststat;// Shorthand
 350       }
 351       if (whatToFit == kSystErrors) hToFit = hsyst;
 352
 353
 354       if(!AliBWTools::Fit(hToFit,func,fitmin,fitmax)) {
 355         cout << " FIT ERROR " << endl;
 356         return;
 357       }
 358       cout << "DRAWING" << endl;
 359       c2->cd();
 360       hsyststat->Draw("same,e2");
 361       // TH1F * hsyststat2 = new TH1F(*hsyststat);
 362       //      hsyststat2->SetFillColor(kWhite);
 363       //      hsyststat2->SetLineColor(kBlack);
 364       //      hsyststat2->Draw("same,e2");
 365       hToFit->Draw("same");
 366       TF1* fitfunc=(TF1*)hToFit->GetListOfFunctions()->At(0);
 367       fitfunc->Draw("same");
 368       fitfunc->SetRange(0,4);
 369       fitfunc->SetLineColor(hSpectra[iCombInStudy][ipart][icharge]->GetLineColor());
 370       if(drawStar)    DrawStar(icharge);
 371       hRatiosToFit[ipart][icharge]=(TH1F*)hToFit->Clone(Form("hRatio%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
 372       hRatiosToFitSyst[ipart][icharge]=(TH1F*)hsyststat->Clone(Form("hRatioSyst%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
 373       // Compute ratio data/function integrating the function in the bin width
 374       for(Int_t iBin=1; iBin<hToFit->GetNbinsX(); iBin++){
 375         Double_t lowLim=hToFit->GetBinLowEdge(iBin);
 376         Double_t highLim=hToFit->GetBinLowEdge(iBin+1);
 377         Double_t contFunc=fitfunc->Integral(lowLim,highLim)/(highLim-lowLim);
 378         Double_t ratio=hToFit->GetBinContent(iBin)/contFunc;
 379         Double_t eratio=hToFit->GetBinError(iBin)/contFunc;
 380         Double_t eratioSyst=hsyststat->GetBinError(iBin)/contFunc;
 381         hRatiosToFitSyst[ipart][icharge]->SetBinContent(iBin,ratio);
 382         hRatiosToFitSyst[ipart][icharge]->SetBinError(iBin,eratioSyst);
 383         hRatiosToFit[ipart][icharge]->SetBinContent(iBin,ratio);
 384         hRatiosToFit[ipart][icharge]->SetBinError(iBin,eratio);
 385       }
 386       c2->Update();
 387       l->AddEntry(hToFit,
 388                   scaleKaons && ipart == kKaon ?
 389                   (TString(partLabel[ipart][icharge])+" #times 2").Data()
 390                   : partLabel[ipart][icharge]);
 391
 392       // populate table
 393       //     Float_t yield  = func->Integral(0.45,1.05);
 394       //     Float_t yieldE = func->IntegralError(0.45,1.05);
 395
 396       Float_t yield  = func->Integral(0.,100);
 397       //Float_t yieldE = func->IntegralError(0.,100);
 398
 399       Double_t yieldTools, yieldETools;
 400       Double_t partialYields[3],partialYieldsErrors[3];
 401       AliBWTools::GetYield(hToFit, func, yieldTools, yieldETools,
 402                            0, 100, partialYields,partialYieldsErrors);
 403       Double_t tslope   = func->GetParameter(2);
 404       Double_t tslopeE  = func->GetParError(2);
 405
 406       table.SetNextCol(partLatex[ipart][icharge]);
 407       //table.SetNextCol(yield,yieldE,-4);
 408       table.SetNextCol(yieldTools, yieldETools,-4);
 409       table.SetNextCol(func->GetParameter(0));
 410       table.SetNextCol(tslope,tslopeE,-4);
 411       table.SetNextCol(func->GetParameter(1),func->GetParError(1),-2);
 412       table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
 413       Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hToFit);
 414       //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
 415       Float_t yieldAbove  = func->Integral(lowestPoint,100);
 416       table.SetNextCol(lowestPoint,-2);
 417       table.SetNextCol((1 - yieldAbove/yield)*100,-1);
 418       Float_t mean, meane;  // Mean using only fit function
 419       Double_t meanDF, meanDFe; // Mean from data + fit
 420       Float_t mean2, mean2e;
 421       if (fitFuncID < kFitPhojet) {
 422         AliBWTools::GetMean      (func, mean,  meane , 0.,100., normPar);
 423         AliBWTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
 424       }
 425       //      AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
 426       AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
 427       table.SetNextCol(mean,   meane  ,-4);
 428       table.SetNextCol(meanDF, meanDFe,-4);
 429
 430       //                         fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
 431       table.InsertRow();
 432
 433
 434       /// TEMP TABLE
 435       tempTable.SetNextCol(partLatex[ipart][icharge]);
 436       tempTable.SetNextCol(yieldTools, yieldETools, -4);
 437       tempTable.SetNextCol(meanDF,  meanDFe ,-4);
 438       //      tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
 439       //      tempTable.SetNextCol(yieldAbove/yield,-2);
 440       tempTable.InsertRow();
 441       c2r->cd();
 442       hRatiosToFitSyst[ipart][icharge]->Draw("e2same");
 443       hRatiosToFit[ipart][icharge]->Draw("esame");
 444
 445     }
 446     c2->cd();
 447     if (flagPreliminary) ALICEWorkInProgress(c2,"","ALICE Preliminary");
 448     l->Draw();
 449     if (flagPreliminary) tf->Draw();
 450     c2r->cd();
 451     if(flagPreliminary) ALICEWorkInProgress(c2r,"","ALICE Preliminary");
 452     l->Draw();
 453     PrintCanvas(c2,printFormats);
 454     PrintCanvas(c2r,printFormats);
 455
 456
 457   }
 458
 459
 460   table.PrintTable("ASCII");
 461
 462   cout << "" << endl;
 463   tempTable.PrintTable("ASCII");
 464
 465
 466
 467 }
 468
 469 void LoadSpectra() {
 470
 471   // Loads spectra for different detectors and corresponding systematic errors.
 472
 473   TFile * f=0;
 474
 475   // Systematic errors: initialize histos
 476   gROOT->cd();
 477   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 478     for(Int_t idet = 0; idet < kNHist; idet++){
 479       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 480         hSystError[idet][ipart][icharge] = new TH1F (TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
 481                                                      TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
 482                                                      nbinsTempl,templBins);
 483         hSystError[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
 484         hSystError[idet][ipart][icharge]->SetLineColor   (color[ipart] );
 485
 486       }
 487
 488     }
 489
 490   }
 491
 492   // Load
 493
 494
 495   // TOF
 496   // Load Efficiencies
 497   f =  new TFile("./Files/effTOFhistosSmooth.root");
 498   TFile * ftrk =  new TFile("./Files/effhistos.root");
 499   TH1D * hEffTrackTOF[kNPart][kNCharge];
 500   TH1D * hEffMatchTOF[kNPart][kNCharge];
 501   hEffTrackTOF[kPion]  [kPos] = (TH1D*) ftrk->Get("hpitrk_pos");
 502   hEffTrackTOF[kKaon]  [kPos] = (TH1D*) ftrk->Get("hkatrk_pos");
 503   hEffTrackTOF[kProton][kPos] = (TH1D*) ftrk->Get("hprtrk_pos");
 504   hEffMatchTOF[kPion]  [kPos] = (TH1D*) f->Get("hpieff_pos");
 505   hEffMatchTOF[kKaon]  [kPos] = (TH1D*) f->Get("hkaeff_pos");
 506   hEffMatchTOF[kProton][kPos] = (TH1D*) f->Get("hpreff_pos");
 507   hEffTrackTOF[kPion]  [kNeg] = (TH1D*) ftrk->Get("hpitrk_neg");
 508   hEffTrackTOF[kKaon]  [kNeg] = (TH1D*) ftrk->Get("hkatrk_neg");
 509   hEffTrackTOF[kProton][kNeg] = (TH1D*) ftrk->Get("hprtrk_neg");
 510   hEffMatchTOF[kPion]  [kNeg] = (TH1D*) f->Get("hpieff_neg");
 511   hEffMatchTOF[kKaon]  [kNeg] = (TH1D*) f->Get("hkaeff_neg");
 512   hEffMatchTOF[kProton][kNeg] = (TH1D*) f->Get("hpreff_neg");
 513
 514   //  f = new TFile("./Files/spectra-pos-y.root");
 515   f = new TFile("./Files/spectraRaw-pos-y.root");
 516   hSpectra[kTOF][kPion]  [kPos]= (TH1F*) f->Get("hpi");
 517   hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
 518   hSpectra[kTOF][kKaon]  [kPos]= (TH1F*) f->Get("hka");
 519   hSpectra[kTOF][kPion]  [kPos]->SetName("hpiPos");
 520   hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
 521   hSpectra[kTOF][kKaon]  [kPos]->SetName("hkaPos");
 522   //f = new TFile("./Files/spectra-neg-y.root");
 523   f = new TFile("./Files/spectraRaw-neg-y.root");
 524   hSpectra[kTOF][kPion]  [kNeg]= (TH1F*) f->Get("hpi");
 525   hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
 526   hSpectra[kTOF][kKaon]  [kNeg]= (TH1F*) f->Get("hka");
 527   hSpectra[kTOF][kPion]  [kNeg]->SetName("hpiNeg");
 528   hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
 529   hSpectra[kTOF][kKaon]  [kNeg]->SetName("hkaNeg");
 530
 531   // Divide for efficiency
 532   hSpectra[kTOF][kPion]  [kPos]->Divide(hEffTrackTOF[kPion]  [kPos]);
 533   hSpectra[kTOF][kKaon]  [kPos]->Divide(hEffTrackTOF[kKaon]  [kPos]);
 534   hSpectra[kTOF][kProton][kPos]->Divide(hEffTrackTOF[kProton][kPos]);
 535   hSpectra[kTOF][kPion]  [kPos]->Divide(hEffMatchTOF[kPion]  [kPos]);
 536   hSpectra[kTOF][kKaon]  [kPos]->Divide(hEffMatchTOF[kKaon]  [kPos]);
 537   hSpectra[kTOF][kProton][kPos]->Divide(hEffMatchTOF[kProton][kPos]);
 538   hSpectra[kTOF][kPion]  [kNeg]->Divide(hEffTrackTOF[kPion]  [kNeg]);
 539   hSpectra[kTOF][kKaon]  [kNeg]->Divide(hEffTrackTOF[kKaon]  [kNeg]);
 540   hSpectra[kTOF][kProton][kNeg]->Divide(hEffTrackTOF[kProton][kNeg]);
 541   hSpectra[kTOF][kPion]  [kNeg]->Divide(hEffMatchTOF[kPion]  [kNeg]);
 542   hSpectra[kTOF][kKaon]  [kNeg]->Divide(hEffMatchTOF[kKaon]  [kNeg]);
 543   hSpectra[kTOF][kProton][kNeg]->Divide(hEffMatchTOF[kProton][kNeg]);
 544
 545
 546   // Clean UP TOF spectra, removing unwanted points
 547   cout << "Cleaning Up TOF spectra" << endl;
 548   Int_t nbin =  hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
 549   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 550     Float_t pt =  hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
 551     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 552       if (pt > 2.35) {
 553         hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
 554         hSpectra[kTOF][kKaon][icharge]->SetBinError  (ibin,0);
 555         hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
 556         hSpectra[kTOF][kProton][icharge]->SetBinError  (ibin,0);
 557       }
 558       if (pt < 0.7) {
 559         hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
 560         hSpectra[kTOF][kProton][icharge]->SetBinError  (ibin,0);
 561       }
 562     }
 563   }
 564 //   cout << "Scaling TOF to TPC" << endl;
 565 //   // Scale protons to TPC
 566 //   hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
 567 //   // Scale pbar so that pbar/p is compatible with Panos
 568 //   hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
 569
 570   // TOF: systematics
 571   // Load TOF systematic errors:
 572   f = new TFile ("./Files/systMatchingPos.root");
 573   AliBWTools::AddHisto(hSystError[kTOF][kPion][kPos]  ,(TH1*)gDirectory->Get("hErrPionMatch")  );
 574   AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hErrProtonMatch"));
 575   AliBWTools::AddHisto(hSystError[kTOF][kKaon][kPos]  ,(TH1*)gDirectory->Get("hErrKaonMatch")  );
 576   f = new TFile ("./Files/systMatchingNeg.root");
 577   AliBWTools::AddHisto(hSystError[kTOF][kPion]  [kNeg],(TH1*)gDirectory->Get("hErrPionMatch"));
 578   AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hErrProtonMatch"));
 579   AliBWTools::AddHisto(hSystError[kTOF][kKaon]  [kNeg],(TH1*)gDirectory->Get("hErrKaonMatch"));
 580   f = new TFile ("./Files/systPIDall.root");
 581   AliBWTools::AddHisto(hSystError[kTOF][kPion]  [kPos],(TH1*)gDirectory->Get("hpiCorr"));
 582   AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hpCorr"));
 583   AliBWTools::AddHisto(hSystError[kTOF][kKaon]  [kPos],(TH1*)gDirectory->Get("hkCorr"));
 584   AliBWTools::AddHisto(hSystError[kTOF][kPion]  [kNeg],(TH1*)gDirectory->Get("hpiCorr"));
 585   AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hpCorr"));
 586   AliBWTools::AddHisto(hSystError[kTOF][kKaon]  [kNeg],(TH1*)gDirectory->Get("hkCorr"));
 587
 588
 589   // ITS SA
 590   f = new TFile("./Files/ITSsaSpectraCorr-28Sep10.root");
 591   hSpectra[kITS][kPion]  [kPos]= (TH1F*) f->Get("hSpectraPos0");
 592   hSpectra[kITS][kKaon]  [kPos]= (TH1F*) f->Get("hSpectraPos1");
 593   hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
 594   hSpectra[kITS][kPion]  [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
 595   hSpectra[kITS][kKaon]  [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
 596   hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
 597
 598   // Remove unwanted bins
 599   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 600     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 601       Int_t nbinITS = hSpectra[kITS][ipart][icharge]->GetNbinsX();
 602       for(Int_t ibin = 1; ibin <= nbinITS; ibin++){
 603         if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
 604           hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 605           hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 606         }
 607         if(ipart == kProton && ibin==9){
 608           printf("Kill bin %d (%f - %f GeV/c)for ITS protons\n",ibin,hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin),hSpectra[kITS][ipart][icharge]->GetBinLowEdge(ibin)+hSpectra[kITS][ipart][icharge]->GetBinWidth(ibin));
 609           hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 610           hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 611         }
 612 //      if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
 613 //        hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 614 //        hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 615 //      }
 616       }
 617
 618     }
 619   }
 620
 621   if(useSecCorrFromDCA){
 622     TFile* fseccorr = new TFile("./Files/CorrFac-SecProtons3Meth-ITSsa.root");
 623     TH1F* hcorrp=(TH1F*)fseccorr->Get("hSecPCorrFromDCAexpo");
 624     TH1F* hcorrpbar=(TH1F*)fseccorr->Get("hSecPbarCorrFromDCAexpo");
 625     hSpectra[kITS][kProton][kPos]->Multiply(hcorrp);
 626     hSpectra[kITS][kProton][kNeg]->Multiply(hcorrpbar);
 627     fseccorr->Close();
 628   }
 629
 630   // Load ITS sa systematics, only pt dependent ones
 631   //  f = TFile::Open("./Files/ITSsa-systematics_20100930.root");
 632   f = TFile::Open("./Files/ITSsa-systematics-20101014.root");
 633   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 634       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 635         AliBWTools::AddHisto(hSystError[kITS][ipart][icharge], (TH1*) gDirectory->Get(Form("hSystTot%s%s",chargeFlag[icharge],partFlag[ipart]))); // Using total error
 636       }
 637     }
 638
 639   // ITS + TPC (Marek)
 640   //f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons20100720.root");
 641   f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons14092010new.root");
 642   TList * list = (TList*) gDirectory->Get("output");
 643   hSpectra[kITSTPC][kPion]  [kPos]= (TH1F*) list->FindObject("Pions");
 644   hSpectra[kITSTPC][kKaon]  [kPos]= (TH1F*) list->FindObject("Kaons");
 645   hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
 646   hSpectra[kITSTPC][kPion]  [kNeg]= (TH1F*) list->FindObject("AntiPions");
 647   hSpectra[kITSTPC][kKaon]  [kNeg]= (TH1F*) list->FindObject("AntiKaons");
 648   hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
 649
 650   // TPC
 651   f = new TFile("./Files/protonSpectra_20100615.root");
 652   hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
 653   hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
 654   f = new TFile("./Files/pionSpectra_20100615.root");
 655   hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
 656   hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
 657   //  f = new TFile("./Files/kaonSpectra_20100615.root");
 658   f = new TFile("./Files/kaonsTPCcorr_2010_08_31.root");
 659   hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
 660   hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
 661
 662   // Clean UP TPC spectra, removing unwanted points
 663   cout << "Cleaning Up TPC spectra" << endl;
 664   nbin =  hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
 665   for(Int_t ibin = 0; ibin < nbin; ibin++){
 666     Float_t pt =  hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
 667     if (pt > 0.45){  // || pt<0.25) {
 668       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 669         hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
 670         hSpectra[kTPC][kKaon][icharge]->SetBinError  (ibin,0);
 671       }
 672     }
 673     if (pt < 0.25) {
 674       hSpectra[kTPC][kKaon][kNeg]->SetBinContent(ibin,0);
 675       hSpectra[kTPC][kKaon][kNeg]->SetBinError  (ibin,0);
 676     }
 677     if (pt < 0.45) {
 678       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 679         hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
 680         hSpectra[kTPC][kProton][icharge]->SetBinError  (ibin,0);
 681       }
 682     }
 683   }
 684
 685   // Load TPC systematics
 686   cout << "WARNING: TPC SYST FOR NEGATIVES TO BE CORRECTED" << endl;
 687   f = TFile ::Open("./Files/pionsSystSum.root");
 688   AliBWTools::AddHisto(hSystError[kTPC][kPion][kPos],(TH1*) gDirectory->Get("pionsSyst"));
 689   AliBWTools::AddHisto(hSystError[kTPC][kPion][kNeg],(TH1*) gDirectory->Get("pionsSyst"));
 690   f = TFile ::Open("./Files/kaonsSystSum.root");
 691   AliBWTools::AddHisto(hSystError[kTPC][kKaon][kPos],(TH1*) gDirectory->Get("kaonsSyst"));
 692   AliBWTools::AddHisto(hSystError[kTPC][kKaon][kNeg],(TH1*) gDirectory->Get("kaonsSyst"));
 693   f = TFile ::Open("./Files/ProtonsSystSum.root");
 694   AliBWTools::AddHisto(hSystError[kTPC][kProton][kPos],(TH1*) gDirectory->Get("ProtonsSyst"));
 695   AliBWTools::AddHisto(hSystError[kTPC][kProton][kNeg],(TH1*) gDirectory->Get("ProtonsSyst"));
 696
 697   // K0s
 698   f = new TFile ("./Files/PtSpectraCorrectedK0sOff_20100803.root");
 699   //  hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
 700   hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
 701   //  hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
 702   //  hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
 703   hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
 704
 705   // Kinks
 706   //  f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
 707   //  f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
 708   //  f = new TFile ("./Files/KaonKinkJun16PhySel2N.root");
 709   f = new TFile ("./Files/KaonKinkJun16PhySel2NRebPass6.root");
 710   hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
 711   hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
 712   // hSpectra[kKinks][kKaon][kPos]->Scale(0.5/0.7); // different rapidity range for kinks
 713   // hSpectra[kKinks][kKaon][kNeg]->Scale(0.5/0.7); // different rapidity range for kinks
 714   // hSpectra[kKinks][kKaon][kPos]->Scale(276004./263345.); // different N of events
 715   // hSpectra[kKinks][kKaon][kNeg]->Scale(276004./263345.); // different N of events
 716   // This normalization if for KaonKinkJun16PhySel2N
 717   // hSpectra[kKinks][kKaon][kPos]->Scale(1./303214); // different N of events
 718   // hSpectra[kKinks][kKaon][kNeg]->Scale(1./303214); // different N of events
 719   // This normalization is for KaonKinkJun16PhySel2NRebPass6
 720   hSpectra[kKinks][kKaon][kPos]->Scale(1./260432.26); // different N of events
 721   hSpectra[kKinks][kKaon][kNeg]->Scale(1./260432.26); // different N of events
 722
 723   // Apply correction factors
 724   // Secondaries for protons
 725
 726   //f = new TFile ("./Files/corrFactorProtons_20100615.root");
 727   f = new TFile ("./Files/corrFactorProtons_2010_09_24.root");
 728   TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
 729   if(correctSecondaries) {
 730     cout << "CORRECTING SECONDARIES" << endl;
 731
 732     for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
 733       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 734         Int_t ipart = kProton;
 735         TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
 736         if (h){
 737           Int_t nbins = h->GetNbinsX();
 738           for(Int_t ibin = 0; ibin < nbins; ibin++){
 739             //      Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
 740             Float_t pt = h->GetBinCenter(ibin);
 741             if (icharge == kNeg) pt = -pt;
 742             Int_t binCorrection = hCorrSecondaries->FindBin(pt);
 743             Float_t correction    = hCorrSecondaries->GetBinContent(binCorrection);
 744             //      cout << pt << " " << correction << endl;
 745
 746             if (correction) {// If the bin is empty this is a  0
 747               h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
 748               h->SetBinError  (ibin,h->GetBinError  (ibin)/correction);
 749             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 750               cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
 751               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 752
 753             }
 754           }
 755         }
 756       }
 757     }
 758   }
 759   // geant/fluka absorption
 760   if(correctGeantFlukaXS) {
 761     cout << "CORRECTING GEANT3/FLUKA" << endl;
 762     // common file for Kaons
 763     TFile *fFlukakaons = TFile::Open("./Files/correctionForCrossSection.321.root");
 764     TH1D * hCorrFlukakaon[kNCharge];
 765     hCorrFlukakaon[kPos] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionParticles");
 766     hCorrFlukakaon[kNeg] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionAntiParticles");
 767
 768     for(Int_t idet = 0; idet < kNDet; idet++){
 769       if( idet != kITS) continue; // comment to use fluka for kaons for all dets
 770       if (idet == kTOF) continue; // TOF already corrected
 771       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 772         Int_t ipart = kKaon;
 773         TH1 * h = hSpectra[idet][ipart][icharge]; // only ITS sa
 774         if (h){
 775           Int_t nbins = h->GetNbinsX();
 776           Int_t nbinsy=hCorrFlukakaon[icharge]->GetNbinsY();
 777           for(Int_t ibin = 0; ibin < nbins; ibin++){
 778             Float_t pt = h->GetBinCenter(ibin);
 779             Float_t minPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(1);
 780             Float_t maxPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(nbinsy+1);
 781             if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 782             if (pt > maxPtCorrection) pt = maxPtCorrection;
 783             Float_t correction = hCorrFlukakaon[icharge]->GetBinContent(hCorrFlukakaon[icharge]->GetXaxis()->FindBin(pt));
 784             if (correction != 0) {// If the bin is empty this is a  0
 785               h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 786               h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 787             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 788               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
 789               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 790             }
 791           }
 792         }
 793       }
 794     }
 795
 796     // PROTONS
 797
 798     // ITS specific file for protons/antiprotons
 799     TFile* fITS = new TFile ("./Files/correctionForCrossSectionITS_20100719.root");
 800     TH2D * hCorrFlukaITS[kNCharge];
 801     hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
 802     hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
 803
 804     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 805       Int_t ipart = kProton;
 806       TH1 * h = hSpectra[kITS][ipart][icharge]; // only ITS sa
 807       if (h){
 808         Int_t nbins = h->GetNbinsX();
 809         Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
 810         for(Int_t ibin = 0; ibin < nbins; ibin++){
 811           Float_t pt = h->GetBinCenter(ibin);
 812           Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
 813           Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 814           if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 815           if (pt > maxPtCorrection) pt = maxPtCorrection;
 816           Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
 817           if (correction != 0) {// If the bin is empty this is a  0
 818             h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 819             h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 820           } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 821             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
 822             cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 823           }
 824         }
 825       }
 826     }
 827
 828
 829
 830     //f = new TFile ("./Files/correctionForCrossSection_20100615.root");
 831     f = new TFile ("./Files/correctionForCrossSection_20100920.root");
 832     TH2D * hCorrFluka[kNCharge];
 833     hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
 834     hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
 835     for(Int_t idet = 0; idet < kNDet; idet++){
 836       if (idet == kITS) continue; // skip ITS sa
 837       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 838         Int_t ipart = kProton;
 839         TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
 840         if (h){
 841           Int_t nbins = h->GetNbinsX();
 842           for(Int_t ibin = 0; ibin < nbins; ibin++){
 843 //          Float_t pt = convertToMT ?
 844 //            TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
 845 //            h->GetBinCenter(ibin);
 846             Float_t pt = h->GetBinCenter(ibin);
 847             Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
 848             //      hCorrFluka[icharge]->Draw();
 849             TH1D * htmpFluka=hCorrFluka[icharge]->ProjectionY();
 850             Float_t maxPtCorrection = AliBWTools::GetHighestNotEmptyBinEdge(htmpFluka);//->GetYaxis()->GetBinLowEdge(hCorrFluka[icharge]->GetNbinsY()+1);
 851             if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 852             if (pt > maxPtCorrection) pt = maxPtCorrection-0.0001;
 853             Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
 854
 855             // already in the efficiency correction (F. Noferini)
 856             if (idet == kTOF) {
 857               if (icharge == kNeg)  correction = 1; // antiprotons already corrected in efficiency
 858               // Scale correction for the different material budget. Recipe by Francesco Noferini
 859               else {
 860                 correction = TMath::Power(correction,0.07162/0.03471);
 861               }
 862             }
 863             if (correction != 0) {// If the bin is empty this is a  0
 864
 865               h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 866               h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 867             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 868               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
 869               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 870             }
 871
 872           }
 873
 874         }
 875       }
 876     }
 877   }
 878
 879
 880   // Set style and scale
 881   for(Int_t idet = 0; idet < kNDet; idet++){
 882     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 883       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 884         if (hSpectra[idet][ipart][icharge]){
 885           hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
 886           hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
 887           hSpectra[idet][ipart][icharge]->SetLineColor   (color[ipart] );
 888           hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
 889           hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
 890           hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
 891           if (convertToMT) {
 892             TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
 893             hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 894             hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
 895             hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
 896           }
 897 //        if (idet == kTOF || idet == kTPC) {
 898 //            hSpectra[idet][ipart][icharge]->Scale(1./236763);
 899 //        }
 900 //        if (idet == kITS){
 901 //          hSpectra[idet][ipart][icharge]->Scale(202945./236763);
 902 //        }
 903           if (scaleKaons && ipart == kKaon) {
 904             hSpectra[idet][ipart][icharge]->Scale(2.);
 905           }
 906         } else {
 907           printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
 908         }
 909       }
 910     }
 911   }
 912
 913
 914   // Create fake weights for the mean; To be update once we have syst errors
 915   // Using syste from table in paper. It would be better to have this as a function of pt.
 916   TH1F * hWeights[3][kNPart];
 917   const Double_t kWeights[3][kNPart] =
 918     // {{4,  3,  10.2},  // TPC
 919     //  {4.1,8.8,7.0 },  //TOF
 920     //  {4.5,5.6,7.0 }}; // ITS
 921     {{0.1,0.1,0.1},  // TPC
 922      {0.1,0.1,0.1},  //TOF
 923      {0.1,0.1,0.1}}; // ITS
 924   for(Int_t ipart = 0; ipart <= kNPart ; ipart++){
 925     for(Int_t idet = 0; idet <= kITS ; idet++){
 926       hWeights[idet][ipart] = (TH1F*) hSpectra[idet][ipart][kPos]->Clone();
 927       Int_t nbinW = hWeights[idet][ipart]->GetNbinsX();
 928       for(Int_t ibin = 1; ibin <= nbinW; ibin++){
 929         hWeights[idet][ipart]->SetBinError(ibin, kWeights[idet][ipart]);
 930       }
 931     }
 932   }
 933   const Double_t scaleDet[] = {1.00,1.00,1.00}; // During the weekly meeting 19/08/2010 we decided not to apply any scaling.
 934   //  const Double_t scaleDet[] = {0.98,1,1.00}; // Scaling factor for the different detectors. Estimated from ratios, it has an estimated uncertainty of ~2%
 935   //  const Double_t scaleDet[] = {0.88,1,0.88}; // Scaling factor for the different detectors. Estimated from ratios, it has an estimated uncertainty of ~2%
 936
 937
 938   // Combine detectors
 939   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 940     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 941       TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
 942
 943       if (convertToMT) {
 944         TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
 945         htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 946
 947       }
 948       hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
 949                                                                         hSpectra[kTPC][ipart][icharge],
 950                                                                         htemplLocal,1.);;
 951
 952       hSpectra[kCombAll][ipart][icharge]   =
 953         AliBWTools::Combine3HistosWithErrors(hSpectra[kITS][ipart][icharge],  // Histos to combine
 954                                              hSpectra[kTPC][ipart][icharge],
 955                                              hSpectra[kTOF][ipart][icharge],
 956                                              hSystError[kITS][ipart][icharge], // Errors (weights) used for the average
 957                                              hSystError[kTPC][ipart][icharge],
 958                                              hSystError[kTOF][ipart][icharge],
 959                                              // hWeights[kITS][ipart],
 960                                              // hWeights[kTPC][ipart],
 961                                              // hWeights[kTOF][ipart],
 962                                              htemplLocal,1,   // template, take statistical error from TPC in overlap
 963                                              scaleDet[kITS],  // relative scaling
 964                                              scaleDet[kTPC],
 965                                              scaleDet[kTOF],
 966                                              (TH1**)&hSystError[kCombAll][ipart][icharge], // put combined syst error here
 967                                              1 // weights histos contain error in bin content
 968                                              );
 969 //       if (convertToMT) {
 970 //      TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
 971 //      hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 972 //      hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
 973 //      hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
 974 //       }
 975     }
 976   }
 977
 978
 979   // Scale for the number of inelastic collisions and correct for
 980   // efficiency losses due to physics selection:
 981
 982   Double_t effPhysSel[kNPart];
 983   effPhysSel[kPion]   = 1.012;
 984   effPhysSel[kKaon]   = 1.013;
 985   effPhysSel[kProton] = 1.014;
 986
 987
 988   for(Int_t idet = 0; idet < kNHist; idet++){
 989     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 990       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 991         if(hSpectra[idet][ipart][icharge]) {
 992           //      cout << "Scaling!" << endl;
 993           if(idet != kKinks && idet != kK0){ // Kinks use a different run list, k0s normalized by Helene
 994             hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel efficiency (computed by Alexander)
 995           }
 996         }
 997       }
 998     }
 999   }
1000
1001
1002 }
1003
1004 void LoadMC() {
1005
1006   TFile * f = 0;
1007   const char * evClass= "INEL";
1008   const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
1009                           "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
1010                           "./Files/dndeta_CMSD6T109_10M_900GeV.root",
1011                           "./Files/dndeta_Perugia0320_10M_900GeV.root", };
1012
1013   // Phojet
1014   for(Int_t itune = 0; itune < kNTunes; itune++){
1015     f = new TFile(files[itune]);
1016     for(Int_t ipart = 0; ipart < kNPart; ipart++){
1017       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1018         hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
1019       }
1020     }
1021   }
1022
1023
1024   // Set style
1025   for(Int_t itune = 0; itune < kNTunes; itune++){
1026     for(Int_t ipart = 0; ipart < kNPart; ipart++){
1027       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1028         if (hSpectraMC[itune][ipart][icharge]){
1029           hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
1030           hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
1031           hSpectraMC[itune][ipart][icharge]->SetLineColor   (mcLineColor[itune] );
1032           hSpectraMC[itune][ipart][icharge]->SetLineStyle   (mcLineStyle[itune] );
1033           hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
1034           if (convertToMT) {
1035             TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
1036             hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
1037             hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
1038             hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
1039           }
1040
1041         } else {
1042           printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
1043         }
1044       }
1045     }
1046   }
1047
1048 }
1049
1050 void DrawStar(Int_t icharge) {
1051   //  cout << icharge << endl;
1052
1053   //  gROOT->LoadMacro("StarPPSpectra.C");
1054   TGraphErrors ** gStar = StarPPSpectra();
1055
1056   for(Int_t istar = 0; istar < 6; istar++){
1057     gStar[istar]->SetMarkerStyle(kOpenStar);
1058     if      (icharge==kPos &&  (istar%2) ) gStar[istar]->Draw("P");
1059     else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
1060     else cout << "Skipping Star " << istar << endl;
1061   }
1062 }
1063
1064 void GetITSResiduals() {
1065
1066
1067   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1068     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1069       //      cout << "1 " << ipart << " " << icharge << endl;
1070
1071 //       gSystem->ProcessEvents();
1072 //       gSystem->Sleep(1000);
1073       TF1* f = (TF1*)   hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
1074       TH1F * hres1, *hres2;
1075       AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
1076
1077       TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
1078       c1->SetLogy();
1079       hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
1080       hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
1081       hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
1082       hSpectra[kITS][ipart][icharge]->Draw("same");
1083       hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
1084       TLegend* l = new TLegend(    0.182886,    0.192308,    0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
1085       l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
1086       l->AddEntry(hSpectra[kITS][ipart][icharge],        "ITS");
1087       l->AddEntry(f,        "Fit to TOF+TPC");
1088       l->Draw();
1089
1090       TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
1091                                  TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
1092       c2->SetGridy();
1093       hres2->SetMinimum(-1);
1094       hres2->SetMaximum(1);
1095       hres2->Draw();
1096       hres2->GetYaxis()->SetTitleOffset(1.2);
1097       Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
1098       TLine * line = new TLine(x,-1,x,1);
1099       line->SetLineStyle(kDashed);
1100       line->Draw("same");
1101
1102       if (doPrint) {
1103         c1->Update();
1104         c2->Update();
1105         gSystem->ProcessEvents();
1106         c1->Print(TString(c1->GetName()) + ".png");
1107         c2->Print(TString(c2->GetName()) + ".png");
1108       }
1109     }
1110   }
1111 }
1112
1113 void DrawWithModels() {
1114
1115   Int_t chargeLoop = sumCharge ? 1 : 2;
1116   for(Int_t icharge = 0; icharge < chargeLoop; icharge++){
1117
1118     // Draw with models
1119     for(Int_t ipart = 0; ipart < kNPart; ipart++){
1120       // Multipad canvas
1121       TString chargeFlagLocal = chargeFlag[icharge];
1122       if(sumCharge) chargeFlagLocal += chargeFlag[1];
1123       TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlagLocal,TString("cSpectra")+partFlag[ipart]+chargeFlagLocal,700,700);
1124       TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlagLocal, "p1", 0.0, 0.35, 1.0,  0.95, 0, 0, 0);
1125       p1->SetBottomMargin(0);
1126       p1->Draw();
1127
1128       TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlagLocal, "p2", 0.0, 0.05, 1.0,  0.35, 0, 0, 0);
1129       p2->SetTopMargin(0);
1130       p2->SetBottomMargin(0.3);
1131       p2->Draw();
1132
1133       Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
1134
1135       // Draw spectra
1136       p1->cd();
1137       p1->SetLogy();
1138       Float_t maxy = sumCharge ? 10 : 5;
1139       TH2F * hempty = new TH2F(TString("hempty")+long(icharge)+long(ipart),"hempty",100,0.,3, 100, 0.0015,maxy);
1140       hempty->SetXTitle("p_{t} (GeV/c)");
1141       hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1142       hempty->Draw();
1143       c1->SetLogy();
1144
1145
1146
1147       TLegend * l =new TLegend(0.502874, 0.493056, 0.892241, 0.904762);
1148       l->SetFillColor(kWhite);
1149       if(sumCharge)     {
1150         hSpectra[iCombInStudy][ipart][kPos]->Add(hSpectra[iCombInStudy][ipart][kNeg]); // Draw pos+neg
1151       }
1152       hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
1153       l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data"));
1154       for(Int_t itune = 0; itune < kNTunes; itune++){
1155         if(sumCharge)   hSpectraMC[itune][ipart][icharge]->Add(hSpectraMC[itune][ipart][kNeg]); // Draw pos+neg;
1156         l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1157         hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
1158         AliBWTools::GetGraphFromHisto(hSpectraMC[itune][ipart][icharge])->Draw("CX");
1159       }
1160       l->Draw("same");
1161
1162       TLatex * tex = new TLatex(0.6712643,2.353486,sumCharge ? Form("%s+%s",partLabel[ipart][icharge],partLabel[ipart][1]) : partLabel[ipart][icharge]);
1163       tex->SetTextFont(42);
1164       tex->SetTextSize(0.07936508);
1165       tex->SetLineWidth(2);
1166       tex->Draw();
1167
1168       // Draw ratio
1169       p2->cd();
1170       Float_t maxyr = sumCharge ? 3.99 : 2.99;
1171       TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(icharge)+long(ipart),"hempty",100,0.,3, 100, 0.01,maxyr);
1172       hemptyr->SetXTitle("p_{t} (GeV/c)");
1173       hemptyr->SetYTitle("Data/MC");
1174       hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
1175       hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
1176       hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
1177       hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
1178       hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
1179       hemptyr->GetXaxis()->SetTitleOffset(1.05);
1180       hemptyr->SetTickLength(0.03*scaleFonts, "X");
1181       hemptyr->SetTickLength(0.02*scaleFonts, "Y");
1182       //      hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
1183       hemptyr->GetYaxis()->SetNdivisions(5);
1184       hemptyr->Draw("");
1185
1186       AliBWFunc fm;
1187       for(Int_t itune = 0; itune < kNTunes; itune++){
1188         TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
1189
1190         //      l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1191         TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[iCombInStudy][ipart][icharge],f);
1192         hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
1193         hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
1194         hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
1195         hRatio->SetMarkerStyle(0);
1196         hRatio->Draw("same");
1197         //      AliBWTools::GetGraphFromHisto(hRatio)->Draw("CX");
1198       }
1199
1200
1201       // print
1202       PrintCanvas(c1,printFormats);
1203     }
1204   }
1205
1206
1207
1208 }
1209
1210 void DrawAllAndKaons() {
1211
1212
1213
1214   //  gROOT->LoadMacro("ALICEWorkInProgress.C");
1215
1216   //  gStyle->SetOptFit(0);
1217   gStyle->SetStatX(0.9);
1218   gStyle->SetStatY(0.88);
1219   gStyle->SetStatW(0.4);
1220   gStyle->SetStatH(0.1);
1221
1222   TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1223   hKaonsAllTPCTOF->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1224
1225   TH1F * hK0Scaled    = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
1226   hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
1227
1228   hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
1229   hSpectra[kKinks][kKaon][kPos]->SetLineColor(4);
1230   hSpectra[kKinks][kKaon][kPos]->SetStats(0);
1231   TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
1232   hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
1233
1234   TCanvas * c1 = new TCanvas("cKaons","cKaons",700,700);
1235   c1->SetLogy();
1236   TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-3,6);
1237   hempty->SetXTitle("p_{t} (GeV/c)");
1238   //  hempty->SetYTitle("dN / dp_{t} (A.U.)");
1239   hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1240   hempty->Draw();
1241   hK0Scaled->Draw("same");
1242   hKaonsAllTPCTOF->Draw("same");
1243   hKinksAll->Draw("same");
1244
1245
1246   TLegend * leg = new TLegend(0.456897, 0.71875, 0.892241, 0.936012,NULL,"brNDC");
1247   //    leg->SetBorderSize(0);
1248 //    leg->SetLineColor(1);
1249 //    leg->SetLineStyle(1);
1250 //    leg->SetLineWidth(1);
1251 //    leg->SetFillColor(19);
1252     leg->SetFillColor(0);
1253    TLegendEntry *entry=leg->AddEntry(hKaonsAllTPCTOF,"K^{+} + K^{-}, ITS+TPC+TOF ","lpf");
1254    entry=leg->AddEntry(hK0Scaled,"K^{0} #times 2","lpf");
1255    entry=leg->AddEntry(hKinksAll,"K^{+} + K ^{-}, Kinks","lpf");
1256    leg->Draw();
1257
1258    //   ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Prelimiary}{Statistical Error Only}");
1259    // TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
1260    // tex->SetTextColor(2);
1261    // tex->SetTextFont(42);
1262    // tex->SetTextSize(0.03496503);
1263    // tex->Draw();
1264
1265    c1->Update();
1266    PrintCanvas(c1, printFormats);
1267
1268   // Draw all "stable" hadrons
1269    Bool_t divideCanvas=kTRUE;
1270   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1271     TCanvas * c1h = 0;
1272     if(divideCanvas) c1h = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],1200,500);
1273     else c1h = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],700,700);
1274     //    c1h->SetLogy();
1275     c1h->SetLeftMargin(0.14);
1276     if(divideCanvas) {
1277       c1h->Divide(3,1);
1278       c1h->cd(1);
1279     }
1280     TH2F * hemptyLoc = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-3,10);
1281     hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1282     hemptyLoc->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1283     hemptyLoc->GetYaxis()->SetTitleOffset(1.35);
1284     hemptyLoc->GetXaxis()->SetTitleOffset(1.1);
1285     hemptyLoc->Draw();
1286     leg = new TLegend(0.482759, 0.489583, 0.896552, 0.925595, NULL,"brNDC");
1287     leg->SetFillColor(0);
1288     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1289       if(divideCanvas) {
1290         cout << "CD!!" << ipart+1 << endl;
1291         c1h->cd(ipart+1);
1292         TH2F * hClone = new TH2F(*hemptyLoc);
1293         hClone->GetXaxis()->SetRangeUser(0,1);
1294         if(ipart == kPion)   hClone->GetYaxis()->SetRangeUser(0,3.5);
1295         if(ipart == kKaon)   hClone->GetYaxis()->SetRangeUser(0,0.28);
1296         if(ipart == kProton) hClone->GetYaxis()->SetRangeUser(0,0.1);
1297         hClone->Draw();
1298       }
1299       for(Int_t idet = 0; idet <= kITSTPC; idet++){
1300         //      if (idet == kITS) continue;
1301         //      if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
1302         hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
1303         hSpectra[idet][ipart][icharge]->Draw("same");
1304         leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet]  + ")","lpf");
1305       }
1306       //      leg->AddLine();
1307     }
1308     if(divideCanvas) c1h->cd(1);
1309     leg->Draw();
1310     //    ALICEWorkInProgress(c1h,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1311     PrintCanvas(c1h,printFormats);
1312   }
1313
1314
1315   //  Draw ratios
1316
1317   // K-/K+ in the different detectors
1318   TCanvas * cpm=new TCanvas("cpm","Kminus/Kplus",700,700);
1319   cpm->Divide(2,2);
1320   cpm->cd(1);
1321   TH1F* hRatioKPKM_TPC=new TH1F(*(hSpectra[kTPC][kKaon][kNeg]));
1322   hRatioKPKM_TPC->SetMinimum(0.5);
1323   hRatioKPKM_TPC->SetMaximum(1.5);
1324   hRatioKPKM_TPC->Divide(hSpectra[kTPC][kKaon][kPos]);
1325   hRatioKPKM_TPC->GetYaxis()->SetTitle("K-/K+ (TPC)");
1326   hRatioKPKM_TPC->Draw();
1327   cpm->cd(2);
1328   TH1F* hRatioKPKM_ITS=new TH1F(*(hSpectra[kITS][kKaon][kNeg]));
1329   hRatioKPKM_ITS->Divide(hSpectra[kITS][kKaon][kPos]);
1330   hRatioKPKM_ITS->SetMinimum(0.5);
1331   hRatioKPKM_ITS->SetMaximum(1.5);
1332   hRatioKPKM_ITS->GetYaxis()->SetTitle("K-/K+ (ITSsa)");
1333   hRatioKPKM_ITS->Draw("");
1334   cpm->cd(3);
1335   TH1F* hRatioKPKM_TOF=new TH1F(*(hSpectra[kTOF][kKaon][kNeg]));
1336   hRatioKPKM_TOF->Divide(hSpectra[kTOF][kKaon][kPos]);
1337   hRatioKPKM_TOF->SetMinimum(0.5);
1338   hRatioKPKM_TOF->SetMaximum(1.5);
1339   hRatioKPKM_TOF->GetYaxis()->SetTitle("K-/K+ (TOF)");
1340   hRatioKPKM_TOF->Draw("");
1341   cpm->cd(4);
1342   TH1F* hRatioKPKM_ITSTPC=new TH1F(*(hSpectra[kITSTPC][kKaon][kNeg]));
1343   hRatioKPKM_ITSTPC->Divide(hSpectra[kITSTPC][kKaon][kPos]);
1344   hRatioKPKM_ITSTPC->SetMinimum(0.5);
1345   hRatioKPKM_ITSTPC->SetMaximum(1.5);
1346   hRatioKPKM_ITSTPC->GetYaxis()->SetTitle("K-/K+ (ITS Global)");
1347   hRatioKPKM_ITSTPC->Draw("");
1348
1349
1350   // ITS/TPC
1351   TH1F * hRatioITSTPC[kNPart][kNCharge];
1352   for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1353     // Create canvas
1354     TCanvas * c1h = new TCanvas(TString("cITSTPCRatio_")+chargeFlag[icharge],TString("cITSTPCRatio_")+chargeFlag[icharge],1200,500);
1355     c1h->Divide(3,1);
1356     c1h->SetGridy();
1357     TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITSsa/TPC ",100,0.,1., 100, 0.5,1.5);
1358     hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1359     hemptyLoc->SetYTitle("ITSsa / TPC");
1360     // Loop over particles
1361     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1362       // Clone histo
1363       hRatioITSTPC[ipart][icharge]=new TH1F(*hSpectra[kITS][ipart][icharge]);
1364       Int_t nBinsITS=hSpectra[kITS][ipart][icharge]->GetNbinsX();
1365       Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1366       // Loop over ITS bins,
1367       for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1368         hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1369         hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1370         Float_t lowPtITS=hSpectra[kITS][ipart][icharge]->GetBinLowEdge(iBin);
1371         Float_t binWidITS=hSpectra[kITS][ipart][icharge]->GetBinWidth(iBin);
1372         // Loop over TPC bins and find overlapping bins to ITS
1373         for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1374           Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1375           Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1376           if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1377             Float_t numer=hSpectra[kITS][ipart][icharge]->GetBinContent(iBin);
1378             Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1379             Float_t enumer=hSpectra[kITS][ipart][icharge]->GetBinError(iBin);
1380             Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1381             Double_t ratio=0.;
1382             Double_t eratio=0.;
1383             if(numer>0. && denom>0.){
1384               ratio=numer/denom;
1385               eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1386             }
1387             hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1388             hRatioITSTPC[ipart][icharge]->SetBinError(iBin,eratio);
1389             break;
1390           }
1391         }
1392       }
1393       c1h->cd(ipart+1);
1394       // hemptyLoc->SetStats(1);
1395       // hemptyLoc->Draw();
1396       hRatioITSTPC[ipart][icharge]->SetYTitle("ITSsa/TPC");
1397       hRatioITSTPC[ipart][icharge]->SetStats(1);
1398       hRatioITSTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1399       hRatioITSTPC[ipart][icharge]->Draw("");
1400       hRatioITSTPC[ipart][icharge]->Fit("pol0","","same");
1401
1402     }
1403     PrintCanvas(c1h,printFormats);
1404   }
1405
1406   // ITS/TPC
1407   TH1F * hRatioITSGlobalTPC[kNPart][kNCharge];
1408   for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1409     // Create canvas
1410     TCanvas * c1h = new TCanvas(TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],1200,500);
1411     c1h->Divide(3,1);
1412     c1h->SetGridy();
1413     TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITS Global/TPC ",100,0.,1., 100, 0.5,1.5);
1414     hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1415     hemptyLoc->SetYTitle("ITS Global / TPC");
1416     // Loop over particles
1417     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1418       // Clone histo
1419       hRatioITSGlobalTPC[ipart][icharge]=new TH1F(*hSpectra[kITSTPC][ipart][icharge]);
1420       Int_t nBinsITS=hSpectra[kITSTPC][ipart][icharge]->GetNbinsX();
1421       Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1422       // Loop over ITS bins,
1423       for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1424         hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1425         hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1426         Float_t lowPtITS=hSpectra[kITSTPC][ipart][icharge]->GetBinLowEdge(iBin);
1427         Float_t binWidITS=hSpectra[kITSTPC][ipart][icharge]->GetBinWidth(iBin);
1428         // Loop over TPC bins and find overlapping bins to ITS
1429         for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1430           Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1431           Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1432           if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1433             Float_t numer=hSpectra[kITSTPC][ipart][icharge]->GetBinContent(iBin);
1434             Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1435             Float_t enumer=hSpectra[kITSTPC][ipart][icharge]->GetBinError(iBin);
1436             Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1437             Double_t ratio=0.;
1438             Double_t eratio=0.;
1439             if(numer>0. && denom>0.){
1440               ratio=numer/denom;
1441               eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1442             }
1443             hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1444             hRatioITSGlobalTPC[ipart][icharge]->SetBinError(iBin,eratio);
1445             break;
1446           }
1447         }
1448       }
1449       c1h->cd(ipart+1);
1450       // hemptyLoc->SetStats(1);
1451       // hemptyLoc->Draw();
1452       hRatioITSGlobalTPC[ipart][icharge]->SetYTitle("ITS Global/TPC");
1453       hRatioITSGlobalTPC[ipart][icharge]->SetStats(1);
1454       hRatioITSGlobalTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1455       hRatioITSGlobalTPC[ipart][icharge]->Draw("");
1456       hRatioITSGlobalTPC[ipart][icharge]->Fit("pol0","","same");
1457
1458     }
1459     PrintCanvas(c1h,printFormats);
1460   }
1461
1462   // TOF/TPC
1463   TH1F * hRatioTOFTPC[kNPart][kNCharge];
1464   for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1465     // create canvas
1466     TCanvas * c1t = new TCanvas(TString("cTOFTPCRatio_")+chargeFlag[icharge],TString("cTOFTPCRatio_")+chargeFlag[icharge],1200,500);
1467     c1t->SetGridy();
1468     c1t->Divide(3,1);
1469     TH2F * hemptyt = new TH2F(TString("hemptyRt")+long(icharge),"TOF/TPC ",100,0.,1., 100, 0.5,1.5);
1470     hemptyt->SetXTitle("p_{t} (GeV/c)");
1471     hemptyt->SetYTitle("TOF / TPC");
1472     //    hemptyt->Draw();
1473     for(Int_t ipart = 0; ipart < kNPart; ipart++) { // loop over particles
1474       // Clone histo
1475       hRatioTOFTPC[ipart][icharge]=new TH1F(*hSpectra[kTOF][ipart][icharge]);
1476       Int_t nBinsTOF=hSpectra[kTOF][ipart][icharge]->GetNbinsX();
1477       Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1478       // Loop over TOF bins
1479       for(Int_t iBin=1; iBin<=nBinsTOF; iBin++){
1480         hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1481         hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1482         Float_t lowPtTOF=hSpectra[kTOF][ipart][icharge]->GetBinLowEdge(iBin);
1483         Float_t binWidTOF=hSpectra[kTOF][ipart][icharge]->GetBinWidth(iBin);
1484         // Loop over TPC bins and find overlapping bins to ITS
1485         for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1486           Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1487           Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1488           if(TMath::Abs(lowPtTOF-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidTOF)<0.001){
1489             Float_t numer=hSpectra[kTOF][ipart][icharge]->GetBinContent(iBin);
1490             Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1491             Float_t enumer=hSpectra[kTOF][ipart][icharge]->GetBinError(iBin);
1492             Float_t edenom=0;//hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1493             Double_t ratio=0.;
1494             Double_t eratio=0.;
1495             if(numer>0. && denom>0.){
1496               ratio=numer/denom;
1497               eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1498             }
1499             hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1500             hRatioTOFTPC[ipart][icharge]->SetBinError(iBin,eratio);
1501             break;
1502           }
1503         }
1504       }
1505       c1t->cd(ipart+1);
1506       hRatioTOFTPC[ipart][icharge]->SetYTitle("TOF/TPC");
1507       hRatioTOFTPC[ipart][icharge]->SetStats(1);
1508       hRatioTOFTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1509       hRatioTOFTPC[ipart][icharge]->Draw("");
1510       hRatioTOFTPC[ipart][icharge]->Fit("pol0","","same");
1511     }
1512     PrintCanvas(c1t,printFormats);
1513   }
1514
1515 }
1516
1517 void DrawWithJacek() {
1518
1519   //1. Convert spectra to dNdeta and sum
1520   TH1F * hsum = (TH1F*) htemplate->Clone();
1521   hsum->Reset();
1522   Int_t idet= iCombInStudy;
1523   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1524     for(Int_t ipart = 0; ipart < kNPart; ipart++){
1525       TH1 * h = hSpectra[idet][ipart][icharge];
1526       Int_t nbin = h->GetNbinsX();
1527       for(Int_t ibin = 1; ibin <= nbin; ibin++){
1528         Double_t pt = h->GetBinCenter(ibin);
1529         Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
1530         Double_t jacobian = pt/mt;
1531         h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
1532         h->SetBinError  (ibin,h->GetBinError  (ibin)*jacobian);
1533         Int_t ibinSum = hsum->FindBin(pt);
1534         Double_t epsilon = 0.0001;
1535         if ( h->GetBinContent(ibin) > 0 &&
1536              (TMath::Abs(h->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1537               TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1538              ) {
1539           cout << "DISCREPANCY IN BIN RANGES" << endl;
1540           cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
1541                << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
1542                << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
1543           cout << "" << endl;
1544         }
1545
1546
1547         hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
1548         hsum->SetBinError  (ibinSum,0);
1549       }
1550       //        hsum->Add(h);
1551     }
1552   }
1553
1554
1555   // Load Jacek and Draw both:
1556   new TFile ("./Files/pt_allcharged_published.root");
1557   TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("dndpt");
1558   gJacek->SetMarkerStyle(20);
1559   gJacek->SetMarkerColor(kGray);
1560   gJacek->Draw("AP");
1561   gJacek->GetHistogram()->SetXTitle("p_{t} (GeV)");
1562   gJacek->GetHistogram()->SetYTitle("d^{2}N/dp_{t}d#eta (GeV^{-1})");
1563   hsum->Draw("same");
1564
1565   TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
1566   new TCanvas();
1567   gRatio->Draw("AP");
1568
1569
1570
1571 }
1572
1573
1574 void DrawRatioToStar() {
1575
1576   // Star data
1577   //  gROOT->LoadMacro("StarPPSpectra.C");
1578   TGraphErrors ** gStar = StarPPSpectra();
1579   gStar[0]->SetMarkerStyle(kOpenStar);
1580   gStar[1]->SetMarkerStyle(kFullStar);
1581   gStar[2]->SetMarkerStyle(kOpenStar);
1582   gStar[3]->SetMarkerStyle(kFullStar);
1583   gStar[4]->SetMarkerStyle(kOpenStar);
1584   gStar[5]->SetMarkerStyle(kFullStar);
1585
1586   // ALICE, INEL -> NSD
1587   Double_t scaleYield = 3.58/3.02; // from paper 2
1588   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1589     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1590       hSpectra[iCombInStudy][ipart][icharge]->Scale(scaleYield);
1591     }
1592   }
1593
1594
1595   TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
1596   TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
1597   hempty->SetXTitle("p_{t} (GeV/c)");
1598   hempty->SetYTitle("ALICE/STAR (NSD)");
1599   hempty->Draw();
1600
1601   TCanvas * c1Comp = new TCanvas("cCompToStarNeg","cCompToStarNeg");
1602   c1Comp->SetLogy();
1603   TH2F * hempty2 = new TH2F(TString("hemptyCompNeg"),"hemptyCompNeg",100,0.,1.5, 100, 0.001,10);
1604   hempty2->SetXTitle("p_{t} (GeV/c)");
1605   hempty2->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1} [NSD]");
1606   hempty2->Draw();
1607
1608   TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
1609   leg->SetBorderSize(0);
1610   leg->SetLineColor(1);
1611   leg->SetLineStyle(1);
1612   leg->SetLineWidth(1);
1613   leg->SetFillColor(0);
1614   leg->SetFillStyle(1001);
1615
1616
1617
1618   c1->cd();
1619   TGraphErrors * g ;
1620   g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[iCombInStudy][kPion][kNeg],1);
1621   g->SetMarkerStyle(kFullCircle);
1622   g->SetMarkerColor(kBlack);
1623   g->Draw("p");
1624   leg->AddEntry(g,"#pi^{-}","lp");
1625   g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[iCombInStudy][kKaon][kNeg],1);
1626   g->SetMarkerStyle(kOpenCircle);
1627   g->SetMarkerColor(kRed);
1628   g->Draw("p");
1629   leg->AddEntry(g,"K^{-}","lp");
1630   g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[iCombInStudy][kProton][kNeg],1);
1631   g->SetMarkerStyle(kOpenSquare);
1632   g->SetMarkerColor(kBlue);
1633   g->Draw("p");
1634   leg->AddEntry(g,"#bar{p}","lp");
1635   leg->Draw();
1636
1637   c1Comp->cd();
1638   gStar[0]->Draw("p");
1639   hSpectra[iCombInStudy][kPion][kNeg]->Draw("same");
1640   gStar[2]->Draw("p");
1641   hSpectra[iCombInStudy][kKaon][kNeg]->Draw("same");
1642   gStar[4]->Draw("p");
1643   hSpectra[iCombInStudy][kProton][kNeg]->Draw("same");
1644
1645
1646
1647   TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
1648   hempty->Draw();
1649   leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
1650   leg->SetBorderSize(0);
1651   leg->SetLineColor(1);
1652   leg->SetLineStyle(1);
1653   leg->SetLineWidth(1);
1654   leg->SetFillColor(0);
1655   leg->SetFillStyle(1001);
1656
1657   TCanvas * c2Comp = new TCanvas("cCompToStarPos","cCompToStarPos");
1658   c2Comp->SetLogy();
1659   hempty2->Draw();
1660
1661   c2->cd();
1662  //  TGraphErrors * g ;
1663   g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[iCombInStudy][kPion][kPos],1);
1664   g->SetMarkerStyle(kFullCircle);
1665   g->SetMarkerColor(kBlack);
1666   g->Draw("p");
1667   leg->AddEntry(g,"#pi^{+}","lp");
1668   g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[iCombInStudy][kKaon][kPos],1);
1669   g->SetMarkerStyle(kOpenCircle);
1670   g->SetMarkerColor(kRed);
1671   g->Draw("p");
1672   leg->AddEntry(g,"K^{+}","lp");
1673   g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[iCombInStudy][kProton][kPos],1);
1674   g->SetMarkerStyle(kOpenSquare);
1675   g->SetMarkerColor(kBlue);
1676   g->Draw("p");
1677   leg->AddEntry(g,"p","lp");
1678   leg->Draw();
1679
1680
1681   c2Comp->cd();
1682   gStar[1]->Draw("p");
1683   hSpectra[iCombInStudy][kPion][kPos]->Draw("same");
1684   gStar[3]->Draw("p");
1685   hSpectra[iCombInStudy][kKaon][kPos]->Draw("same");
1686   gStar[5]->Draw("p");
1687   hSpectra[iCombInStudy][kProton][kPos]->Draw("same");
1688
1689
1690   c1->Update();
1691   c2->Update();
1692   gSystem->ProcessEvents();
1693   c1->Print(TString(c1->GetName()) + ".eps");
1694   c2->Print(TString(c2->GetName()) + ".eps");
1695   ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1696   ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1697   c1->Update();
1698   c2->Update();
1699   c1->Print(TString(c1->GetName()) + ".png");
1700   c2->Print(TString(c2->GetName()) + ".png");
1701
1702
1703
1704
1705 }
1706
1707
1708
1709 void DrawRatios() {
1710
1711   // Draws ratios of combined spectra
1712
1713   // Rebin histograms (hi pt was fluctuating too much in the TOF)
1714   RebinHistosForRatios();
1715
1716   // Compute ratios
1717   TH1F * hPosNegRatio[kNPart];
1718
1719   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1720     hPosNegRatio[ipart] = (TH1F*) hSpectra[iCombInStudy][ipart][kPos]->Clone();
1721     hPosNegRatio[ipart]->Divide(hSpectra[iCombInStudy][ipart][kNeg]);
1722     hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
1723     hPosNegRatio[ipart]->SetMinimum(0.5);
1724     hPosNegRatio[ipart]->SetMaximum(1.5);
1725   }
1726
1727
1728   TH1F * hKPiRatio = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1729   hKPiRatio->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1730   TH1F * htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1731   htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1732   hKPiRatio->Divide(htmp);
1733   hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1734   // Get ratio of syst errors
1735   // FIXME: commented combined error
1736   TH1F * hKPiRatioSyst = new TH1F(*htemplate);;
1737   AliBWTools::GetValueAndError(hKPiRatioSyst,hSpectra[iCombInStudy][kKaon][kPos],hSystError[iCombInStudy][kKaon][kPos],kTRUE);
1738   //  AliBWTools::GetHistoCombinedErrors(hKPiRatioSyst,hSpectra[iCombInStudy][kKaon][kPos]);
1739   AliBWTools::GetValueAndError(htmp,hSpectra[iCombInStudy][kKaon][kNeg],hSystError[iCombInStudy][kKaon][kNeg],kTRUE);
1740   // //  AliBWTools::GetHistoCombinedErrors(htmp,hSpectra[iCombInStudy][kKaon][kNeg]);
1741   hKPiRatioSyst->Add(htmp);
1742   AliBWTools::GetValueAndError(htmp,hSpectra[iCombInStudy][kPion][kNeg],hSystError[iCombInStudy][kPion][kNeg],kTRUE);
1743   //  AliBWTools::GetHistoCombinedErrors(htmp,hSpectra[iCombInStudy][kPion][kNeg]);
1744   TH1F * htmp2 = new TH1F(*htemplate);
1745   AliBWTools::GetValueAndError(htmp2,hSpectra[iCombInStudy][kPion][kPos],hSystError[iCombInStudy][kPion][kPos],kTRUE);
1746   //  AliBWTools::GetHistoCombinedErrors(htmp2,hSpectra[iCombInStudy][kPion][kPos]);
1747   htmp->Add(htmp2);
1748   hKPiRatioSyst->Divide(htmp);
1749   hKPiRatioSyst->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1750   // Multiply error for sqrt(2): the syst is not reduced by combining positive and negative charges
1751   Int_t nbin = hKPiRatioSyst->GetNbinsX();
1752   for(Int_t ibin = 1; ibin <= nbin; ibin++){
1753     hKPiRatioSyst->SetBinError(ibin, hKPiRatioSyst->GetBinError(ibin)*TMath::Sqrt(2));
1754   }
1755
1756
1757
1758
1759   TH1F * hKPiRatioMC[kNTunes];
1760   if(showMC){
1761     for(Int_t itune = 0; itune < kNTunes; itune++){
1762       hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
1763       hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
1764       TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1765       htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1766       hKPiRatioMC[itune]->Divide(htmpLoc);
1767       hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1768     }
1769   }
1770
1771
1772
1773   TH1F * hPPiRatio = (TH1F*) hSpectra[iCombInStudy][kProton][kPos]->Clone();
1774   hPPiRatio->Add(hSpectra[iCombInStudy][kProton][kNeg]);
1775   htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1776   htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1777   hPPiRatio->Divide(htmp);
1778   hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1779
1780   TH1F * hPPiRatioMC[kNTunes];
1781   if(showMC){
1782     for(Int_t itune = 0; itune < kNTunes; itune++){
1783       hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
1784       hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
1785       TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1786       htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1787       hPPiRatioMC[itune]->Divide(htmpLoc);
1788       hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1789     }
1790   }
1791
1792   // Phenix ratio
1793   TFile * fPhenix = new TFile("./Files/phenixdata.root");
1794   TH1F * hPPiRatioPhenix = (TH1F*) fPhenix->Get("ProtonPHNX");
1795   hPPiRatioPhenix->Add((TH1F*) fPhenix->Get("PbarPHNX"));
1796   htmp = (TH1F*) fPhenix->Get("PiPlusPHNX");
1797   htmp->Add((TH1F*) fPhenix->Get("PiMinusPHNX"));
1798   hPPiRatioPhenix = AliBWTools::DivideHistosDifferentBins(hPPiRatioPhenix,htmp);
1799   hPPiRatioPhenix->SetMarkerStyle(24);
1800
1801   // Draw
1802 //   TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
1803 //   hempty->SetXTitle("p_{t} (GeV/c)");
1804   //  hempty->SetYTitle("");
1805
1806   // tmp: overlay levi fits
1807   AliBWFunc * fm2 = new AliBWFunc;
1808   fm2->SetVarType(AliBWFunc::kdNdpt);
1809   TF1 * fLevi[kNPart]  [kNCharge];
1810   fLevi[kPion]  [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
1811   fLevi[kKaon]  [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
1812   fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
1813   fLevi[kPion]  [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
1814   fLevi[kKaon]  [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
1815   fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
1816   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1817     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1818       fLevi[ipart][icharge]->SetRange(0,4);
1819     }
1820   }
1821
1822
1823   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1824     TString detName = detFlag[iCombInStudy];
1825     detName.ReplaceAll(" ", "_");
1826     detName.ReplaceAll("+", "");
1827
1828     TCanvas * c1 = new TCanvas(TString("cRatio_")+detName+TString("_")+partFlag[ipart], TString("cRatio_")+detName+partFlag[ipart]);
1829     c1->SetGridy();
1830     hPosNegRatio[ipart]->Draw();
1831     TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
1832     //    fRatio->Draw("same");
1833     fRatio->SetRange(0,5);
1834     PrintCanvas(c1,printFormats);
1835
1836   }
1837
1838
1839   TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
1840   //  c2->SetGridy();
1841
1842   hKPiRatioSyst->SetMaximum(0.62);
1843   hKPiRatioSyst->SetFillColor(kGray);
1844   hKPiRatioSyst->SetMarkerStyle(24);
1845   //  hKPiRatioSyst->Draw("e2");
1846   // TMP FIXME: ERROR FROM FP:
1847   TFile * f = TFile::Open("./Files/ITSsa-systematics-20101014.root");
1848   TH1F * hsystRatio = (TH1F*) f->Get("hSystRatioKpi");
1849   TH1F * hsyst = new TH1F(*htemplate);
1850   AliBWTools::GetValueAndError(hsyst,hKPiRatio,hsystRatio,1);
1851   //  AliBWTools::GetHistoCombinedErrors(hsyst,hKPiRatio);
1852   //  hsyst->Draw("samee1");
1853   // END TMP
1854   hKPiRatio->SetMaximum(0.62);
1855   hKPiRatio->Draw("");
1856   TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1857   lMC->SetFillColor(kWhite);
1858   TGraphErrors ** gStar = StarPPSpectra();
1859   cout << gStar[6] << endl;
1860
1861   gStar[6]->SetMarkerStyle(kOpenStar);
1862   gStar[6]->Draw("P");
1863
1864   if(showE735){
1865     gROOT->LoadMacro("GetE735Ratios.C");
1866     // GetE735Ratios(0,0)->Draw("EX0,same");
1867     // GetE735Ratios(0,1)->Draw("EX0,same");
1868     // GetE735Ratios(0,2)->Draw("EX0,same");
1869     GetE735Ratios(0,3)->Draw("EX0,same");
1870   }
1871   hKPiRatio->SetMarkerStyle(20);
1872   hKPiRatio->SetMarkerColor(kRed);
1873   hKPiRatio->SetLineColor(kRed);
1874   //hKPiRatio->Draw("same");
1875
1876   if(showMC){
1877     for(Int_t itune = 0; itune < kNTunes; itune++){
1878       lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1879       hKPiRatioMC[itune]->SetLineWidth(2);
1880       hKPiRatioMC[itune]->Draw("same,chist");
1881     }
1882
1883     lMC->Draw();
1884   }
1885   cout << "12" << endl;
1886
1887   if(showE735){
1888     TLegend * l = new TLegend(  0.1879,  0.68,  0.54,0.92);
1889     l->SetFillColor(kWhite);
1890     l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1891     l->AddEntry(gStar[6],  "STAR, #sqrt{s} = 200 GeV", "lp");
1892     // l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
1893     // l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
1894     // l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
1895     l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
1896     l->Draw();
1897   }
1898
1899
1900   TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
1901   cout << "DRAWING" << endl;
1902   //  c3->SetGridy();
1903   hPPiRatio->SetMarkerStyle(20);
1904   hPPiRatio->SetMarkerColor(kRed);
1905   hPPiRatio->SetLineColor(kRed);
1906   hPPiRatio->Draw();
1907   hPPiRatio->SetMaximum(0.39);
1908   if(showMC){
1909     lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1910     lMC->SetFillColor(kWhite);
1911
1912     for(Int_t itune = 0; itune < kNTunes; itune++){
1913       lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1914       hPPiRatioMC[itune]->SetLineWidth(2);
1915       hPPiRatioMC[itune]->Draw("same,chist");
1916     }
1917
1918     lMC->Draw();
1919   }
1920   hPPiRatioPhenix->Draw("same");
1921
1922
1923   TLegend * l = new TLegend(0.186242, 0.781469, 0.538591, 0.921329);
1924   l->SetFillColor(kWhite);
1925   l->AddEntry(hPPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1926   l->AddEntry(hPPiRatioPhenix, "PHENIX, #sqrt{s} = 200 GeV","lpf");
1927   //  TLegend * l = new TLegend(0.206376, 0.77972, 0.600671, 0.909091);
1928   l->Draw();
1929
1930   PrintCanvas(c2,printFormats);
1931   PrintCanvas(c3,printFormats);
1932
1933 }
1934
1935 void DrawWithSyst() {
1936
1937   // Draws detector and combined with syst errors.
1938
1939   for(Int_t idet = 0; idet < kNHist; idet++){
1940     if(idet > kITS && idet != iCombInStudy) continue;
1941     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1942       TCanvas * c = new TCanvas(TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge],TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge]);
1943       TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
1944       hempty->SetXTitle("p_{t} (GeV/c)");
1945       hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1946       hempty->GetYaxis()->SetTitleOffset(1.35);
1947       hempty->GetXaxis()->SetTitleOffset(1.1);
1948       hempty->Draw();
1949       c->SetLogy();
1950       for(Int_t ipart = 0; ipart < kNPart; ipart++){
1951         //      cout << detFlag[idet] << " " << chargeFlag[icharge] << " " << partFlag[ipart] << endl;
1952
1953         TString opt = ipart ? "" : "same";
1954         TH1F * hsyst = new TH1F(*htemplate);
1955         AliBWTools::GetValueAndError(hsyst,hSpectra[idet][ipart][icharge],hSystError[idet][ipart][icharge],kTRUE);
1956         hsyst->SetFillColor(kGray);
1957         hsyst->Draw("e5,same");
1958         hSpectra[idet][ipart][icharge]->Draw("same");
1959         // Do not draw syst error outside of spectra range
1960         Float_t lowEdge = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1961         Float_t hiEdge  = AliBWTools::GetHighestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1962         Int_t lowBin = hSystError[idet][ipart][icharge]->FindBin(lowEdge);
1963         Int_t hiBin = hSystError[idet][ipart][icharge]->FindBin(hiEdge);
1964         Int_t nbin =    hSystError[idet][ipart][icharge]->GetNbinsX();
1965         for(Int_t ibin = 1; ibin < lowBin; ibin++){
1966           hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1967           hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1968         }
1969         for(Int_t ibin = hiBin; ibin <= nbin; ibin++){
1970           hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1971           hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1972         }
1973         // Draw relative syst error, convert to graph to avoid vertical bars at the edges (empty bins in histo)
1974         AliBWTools::GetGraphFromHisto(hSystError[idet][ipart][icharge])->Draw("LX");
1975
1976
1977
1978       }
1979       PrintCanvas(c,"png,eps");
1980     }
1981   }
1982 }
1983
1984 void Help() {
1985
1986   cout << "Macro: void CombineSpectra(Int_t analysisType=kDoFits, Int_t  fitFuncID = kFitLevi) " << endl;
1987   cout << "" << endl;
1988   cout << "WARNING: some function modify the histos upon execution. You are advised to quit root" << endl;
1989   cout << "         before running a different function" << endl;
1990   cout << "" << endl;
1991   cout << "Possible Arguments" << endl;
1992   cout << "- analysisType:" << endl;
1993   cout << "    kDoFits:           Fit Combined Spectra " << endl;
1994   cout << "    kDoRatios:         Particle ratios K/pi and p/pi" << endl;
1995   cout << "    kDoSuperposition:  Compare different detectors (superimpose and ratios)" << endl;
1996   cout << "    kDoCompareStar:    Compare combined spectra to star results" << endl;
1997   cout << "    kDoDrawWithModels: Compare combined spectra and models" << endl;
1998   cout << "    kDoDrawSyst:       Draws spectra from individual detectors with their systematic error" << endl;
1999   cout << "    kDoCompareToAllCharged:       convert to dn/deta/dpt, sum up and compare to all charged" << endl;
2000   cout << "    kDoHelp:           This help" << endl;
2001   cout << "- fitFuncID, function used to extrapolate and compute yields" << endl;
2002   cout << "    An analitic fit function [kFitLevi, kFitUA1, kFitPowerLaw]" << endl;
2003   cout << "    Or a shape from a MC moder [kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0]" << endl;
2004   cout << "    Which is fitted to the data at low pt and used to extrapolate at low pt" << endl;
2005
2006
2007 }
2008
2009 void PrintCanvas(TCanvas* c,const TString formats) {
2010   // print a canvas in every of the given comma-separated formats
2011   // ensure the canvas is updated
2012   if(!doPrint) return;
2013   c->Update();
2014   gSystem->ProcessEvents();
2015   TObjArray * arr = formats.Tokenize(",");
2016   TIterator * iter = arr->MakeIterator();
2017   TObjString * element = 0;
2018   TString name  =c ->GetName();
2019   name.ReplaceAll(" ","_");
2020   name.ReplaceAll("+","Plus");
2021   name.ReplaceAll("-","");
2022   while ((element = (TObjString*) iter->Next())) {
2023     c->Print(name+ "."+element->GetString());
2024   }
2025 }
2026
2027
2028 void LoadLibs(){
2029
2030   gSystem->Load("libTree.so");
2031   gSystem->Load("libVMC.so");
2032   gSystem->Load("libMinuit.so");
2033   gSystem->Load("libSTEERBase.so");
2034   gSystem->Load("libESD.so");
2035   gSystem->Load("libAOD.so");
2036   gSystem->Load("libANALYSIS.so");
2037   gSystem->Load("libANALYSISalice.so");
2038   gSystem->Load("libCORRFW.so");
2039   gSystem->Load("libPWG2spectra.so");
2040
2041 }
2042
2043 void RebinHistosForRatios() {
2044   // rebins all spectra histograms
2045
2046   cout << "WARNING: REBINNING HISTOS BEFORE RATIOS" << endl;
2047
2048
2049   const Float_t rebinBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.6,1.8,2.4};
2050   Int_t nrebinBins=29;
2051
2052
2053
2054   Int_t idet = iCombInStudy;
2055   for(Int_t ipart = 0; ipart < kNPart; ipart++){
2056     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
2057
2058       TH1F * h =  hSpectra[idet][ipart][icharge];
2059       TString hname = h->GetName();
2060       hname += "_rebin";
2061
2062       TH1F *hrebin   = new TH1F(hname, hname,nrebinBins,rebinBins);
2063       hrebin->Sumw2();
2064
2065       for(Int_t i=1;i <= h->GetNbinsX();i++){
2066         Float_t x = h->GetBinCenter(i);
2067         Int_t j=0;
2068         while(x > rebinBins[j+1]) j++;
2069         j++;
2070         hrebin->AddBinContent(j,h->GetBinContent(i));
2071
2072         Float_t err = hrebin->GetBinError(j) + h->GetBinError(i)*h->GetBinError(i);
2073         hrebin->SetBinError(j,err);
2074       }
2075
2076       for(Int_t i=1;i <= hrebin->GetNbinsX();i++){
2077         hrebin->SetBinError(i,TMath::Sqrt(hrebin->GetBinError(i)));
2078       }
2079       hSpectra[idet][ipart][icharge] = hrebin;
2080     }
2081   }
2082
2083 }