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