PWG2/SPECTRA/Fit/CombineSpectra.C

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