PWGLF/SPECTRA/PiKaPr/COMBINED/CombineSpectra.C

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