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 =  kTRUE; // 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(4,"c|ccc");
 225   tempTable.InsertCustomRow("Part & Yield & Yield Below & Frac Above\\\\");
 226
 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(partLabel[ipart][icharge]);
 376       tempTable.SetNextCol(yieldTools, yieldETools, -4);
 377       tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
 378       tempTable.SetNextCol(yieldAbove/yield,-2);
 379       tempTable.InsertRow();
 380       c2r->cd();
 381       hRatiosToFit[ipart][icharge]->Draw("esame");
 382
 383     }
 384     c2->cd();
 385     l->Draw();
 386     c2r->cd();
 387     l->Draw();
 388     if (doPrint) {
 389       c2->cd();
 390       c2->Update();
 391       gSystem->ProcessEvents();
 392       tf->Draw();
 393       c2->Print(TString(c2->GetName()) + ".eps");
 394       ALICEWorkInProgress(c2,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
 395       c2->Update();
 396       c2->Print(TString(c2->GetName()) + ".png");
 397       c2r->Update();
 398       gSystem->ProcessEvents();
 399       c2r->Print(TString(c2r->GetName()) + ".eps");
 400       c2r->Print(TString(c2r->GetName()) + ".png");
 401     }
 402
 403
 404   }
 405
 406
 407   table.PrintTable("");
 408
 409   cout << "" << endl;
 410   tempTable.PrintTable("ASCII");
 411
 412
 413
 414 }
 415
 416 void LoadSpectra() {
 417
 418   TFile * f=0;
 419
 420   // Load
 421
 422
 423   // TOF
 424   // Load Efficiencies
 425   f =  new TFile("./Files/effhistos.root");
 426   TH1D * hEffTrackTOF[kNPart][kNCharge];
 427   TH1D * hEffMatchTOF[kNPart][kNCharge];
 428   hEffTrackTOF[kPion]  [kPos] = (TH1D*) f->Get("hpitrk_pos");
 429   hEffTrackTOF[kKaon]  [kPos] = (TH1D*) f->Get("hkatrk_pos");
 430   hEffTrackTOF[kProton][kPos] = (TH1D*) f->Get("hprtrk_pos");
 431   hEffMatchTOF[kPion]  [kPos] = (TH1D*) f->Get("hpieff_pos");
 432   hEffMatchTOF[kKaon]  [kPos] = (TH1D*) f->Get("hkaeff_pos");
 433   hEffMatchTOF[kProton][kPos] = (TH1D*) f->Get("hpreff_pos");
 434   hEffTrackTOF[kPion]  [kNeg] = (TH1D*) f->Get("hpitrk_neg");
 435   hEffTrackTOF[kKaon]  [kNeg] = (TH1D*) f->Get("hkatrk_neg");
 436   hEffTrackTOF[kProton][kNeg] = (TH1D*) f->Get("hprtrk_neg");
 437   hEffMatchTOF[kPion]  [kNeg] = (TH1D*) f->Get("hpieff_neg");
 438   hEffMatchTOF[kKaon]  [kNeg] = (TH1D*) f->Get("hkaeff_neg");
 439   hEffMatchTOF[kProton][kNeg] = (TH1D*) f->Get("hpreff_neg");
 440
 441   //  f = new TFile("./Files/spectra-pos-y.root");
 442   f = new TFile("./Files/spectraRaw-pos-y.root");
 443   hSpectra[kTOF][kPion]  [kPos]= (TH1F*) f->Get("hpi");
 444   hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
 445   hSpectra[kTOF][kKaon]  [kPos]= (TH1F*) f->Get("hka");
 446   hSpectra[kTOF][kPion]  [kPos]->SetName("hpiPos");
 447   hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
 448   hSpectra[kTOF][kKaon]  [kPos]->SetName("hkaPos");
 449   //f = new TFile("./Files/spectra-neg-y.root");
 450   f = new TFile("./Files/spectraRaw-neg-y.root");
 451   hSpectra[kTOF][kPion]  [kNeg]= (TH1F*) f->Get("hpi");
 452   hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
 453   hSpectra[kTOF][kKaon]  [kNeg]= (TH1F*) f->Get("hka");
 454   hSpectra[kTOF][kPion]  [kNeg]->SetName("hpiNeg");
 455   hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
 456   hSpectra[kTOF][kKaon]  [kNeg]->SetName("hkaNeg");
 457
 458   // Divide for efficiency
 459   hSpectra[kTOF][kPion]  [kPos]->Divide(hEffTrackTOF[kPion]  [kPos]);
 460   hSpectra[kTOF][kKaon]  [kPos]->Divide(hEffTrackTOF[kKaon]  [kPos]);
 461   hSpectra[kTOF][kProton][kPos]->Divide(hEffTrackTOF[kProton][kPos]);
 462   hSpectra[kTOF][kPion]  [kPos]->Divide(hEffMatchTOF[kPion]  [kPos]);
 463   hSpectra[kTOF][kKaon]  [kPos]->Divide(hEffMatchTOF[kKaon]  [kPos]);
 464   hSpectra[kTOF][kProton][kPos]->Divide(hEffMatchTOF[kProton][kPos]);
 465   hSpectra[kTOF][kPion]  [kNeg]->Divide(hEffTrackTOF[kPion]  [kNeg]);
 466   hSpectra[kTOF][kKaon]  [kNeg]->Divide(hEffTrackTOF[kKaon]  [kNeg]);
 467   hSpectra[kTOF][kProton][kNeg]->Divide(hEffTrackTOF[kProton][kNeg]);
 468   hSpectra[kTOF][kPion]  [kNeg]->Divide(hEffMatchTOF[kPion]  [kNeg]);
 469   hSpectra[kTOF][kKaon]  [kNeg]->Divide(hEffMatchTOF[kKaon]  [kNeg]);
 470   hSpectra[kTOF][kProton][kNeg]->Divide(hEffMatchTOF[kProton][kNeg]);
 471
 472
 473   // Clean UP TOF spectra, removing unwanted points
 474   cout << "Cleaning Up TOF spectra" << endl;
 475   Int_t nbin =  hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
 476   for(Int_t ibin = 1; ibin <= nbin; ibin++){
 477     Float_t pt =  hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
 478     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 479       if (pt > 2.35) {
 480         hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
 481         hSpectra[kTOF][kKaon][icharge]->SetBinError  (ibin,0);
 482         hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
 483         hSpectra[kTOF][kProton][icharge]->SetBinError  (ibin,0);
 484       }
 485     }
 486   }
 487 //   cout << "Scaling TOF to TPC" << endl;
 488 //   // Scale protons to TPC
 489 //   hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
 490 //   // Scale pbar so that pbar/p is compatible with Panos
 491 //   hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
 492
 493
 494
 495   // ITS SA
 496   f = new TFile("./Files/ITSsaSpectraCorr_20100727.root");
 497   hSpectra[kITS][kPion]  [kPos]= (TH1F*) f->Get("hSpectraPos0");
 498   hSpectra[kITS][kKaon]  [kPos]= (TH1F*) f->Get("hSpectraPos1");
 499   hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
 500   hSpectra[kITS][kPion]  [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
 501   hSpectra[kITS][kKaon]  [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
 502   hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
 503
 504   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 505     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 506       Int_t nbin = hSpectra[kITS][ipart][icharge]->GetNbinsX();
 507       for(Int_t ibin = 1; ibin <= nbin; ibin++){
 508         if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
 509           hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 510           hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 511         }
 512         if(ipart == kProton && ibin==9){
 513           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));
 514           hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 515           hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 516         }
 517 //      if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
 518 //        hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
 519 //        hSpectra[kITS][ipart][icharge]->SetBinError  (ibin,0);
 520 //      }
 521       }
 522
 523     }
 524   }
 525
 526   if(useSecCorrFromDCA){
 527     TFile* fseccorr = new TFile("./Files/CorrFac-SecProtonsFromDCA-ITSsa.root");
 528     TH1F* hcorrp=(TH1F*)fseccorr->Get("hSecPCorrFromDCA");
 529     TH1F* hcorrpbar=(TH1F*)fseccorr->Get("hSecPbarCorrFromDCA");
 530     hSpectra[kITS][kProton][kPos]->Multiply(hcorrp);
 531     hSpectra[kITS][kProton][kNeg]->Multiply(hcorrpbar);
 532     fseccorr->Close();
 533   }
 534
 535   // ITS + TPC (Marek)
 536   f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons20100720.root");
 537   TList * list = (TList*) gDirectory->Get("output");
 538   hSpectra[kITSTPC][kPion]  [kPos]= (TH1F*) list->FindObject("Pions");
 539   hSpectra[kITSTPC][kKaon]  [kPos]= (TH1F*) list->FindObject("Kaons");
 540   hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
 541   hSpectra[kITSTPC][kPion]  [kNeg]= (TH1F*) list->FindObject("AntiPions");
 542   hSpectra[kITSTPC][kKaon]  [kNeg]= (TH1F*) list->FindObject("AntiKaons");
 543   hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
 544
 545   // TPC
 546   f = new TFile("./Files/protonSpectra_20100615.root");
 547   hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
 548   hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
 549   f = new TFile("./Files/pionSpectra_20100615.root");
 550   hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
 551   hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
 552   f = new TFile("./Files/kaonSpectra_20100615.root");
 553   hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
 554   hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
 555
 556   // Clean UP TPC spectra, removing unwanted points
 557   cout << "Cleaning Up TPC spectra" << endl;
 558   nbin =  hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
 559   for(Int_t ibin = 0; ibin < nbin; ibin++){
 560     Float_t pt =  hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
 561     if (pt > 0.45){  // || pt<0.25) {
 562       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 563         hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
 564         hSpectra[kTPC][kKaon][icharge]->SetBinError  (ibin,0);
 565       }
 566     }
 567     if (pt < 0.45) {
 568       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 569         hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
 570         hSpectra[kTPC][kProton][icharge]->SetBinError  (ibin,0);
 571       }
 572     }
 573   }
 574
 575
 576
 577
 578
 579   // K0s
 580   f = new TFile ("./Files/PtSpectraCorrectedK0sOff.root");
 581   //  hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
 582   hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
 583   //  hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
 584   //  hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
 585   hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
 586
 587   // Kinks: TO BE FIXED WITH POSITIVES AND NEGATIVES
 588   //  f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
 589   f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
 590   hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
 591   hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
 592   hSpectra[kKinks][kKaon][kPos]->Scale(0.5/0.7); // different rapidity range for kinks
 593   hSpectra[kKinks][kKaon][kNeg]->Scale(0.5/0.7); // different rapidity range for kinks
 594   hSpectra[kKinks][kKaon][kPos]->Scale(276004./263345.); // different N of events
 595   hSpectra[kKinks][kKaon][kNeg]->Scale(276004./263345.); // different N of events
 596
 597   // Apply correction factors
 598   // Secondaries for protons
 599
 600   f = new TFile ("./Files/corrFactorProtons_20100615.root");
 601   TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
 602   if(correctSecondaries) {
 603     cout << "CORRECTING SECONDARIES" << endl;
 604
 605     for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
 606       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 607         Int_t ipart = kProton;
 608         TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
 609         if (h){
 610           Int_t nbins = h->GetNbinsX();
 611           for(Int_t ibin = 0; ibin < nbins; ibin++){
 612             //      Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
 613             Float_t pt = h->GetBinCenter(ibin);
 614             if (icharge == kNeg) pt = -pt;
 615             Int_t binCorrection = hCorrSecondaries->FindBin(pt);
 616             Float_t correction    = hCorrSecondaries->GetBinContent(binCorrection);
 617             //      cout << pt << " " << correction << endl;
 618
 619             if (correction) {// If the bin is empty this is a  0
 620               h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
 621               h->SetBinError  (ibin,h->GetBinError  (ibin)/correction);
 622             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 623               cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
 624               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 625
 626             }
 627           }
 628         }
 629       }
 630     }
 631   }
 632   // geant/fluka absorption
 633   if(correctGeantFlukaXS) {
 634     cout << "CORRECTING GEANT3/FLUKA" << endl;
 635     // common file for Kaons
 636     TFile *fFlukakaons = TFile::Open("./Files/correctionForCrossSection.321.root");
 637     TH1D * hCorrFlukakaon[kNCharge];
 638     hCorrFlukakaon[kPos] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionParticles");
 639     hCorrFlukakaon[kNeg] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionAntiParticles");
 640
 641     for(Int_t idet = 0; idet < kNDet; idet++){
 642       if( idet != kITS) continue; // comment to use fluka for kaons for all dets
 643       if (idet == kTOF) continue; // TOF already corrected
 644       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 645         Int_t ipart = kKaon;
 646         TH1 * h = hSpectra[idet][ipart][icharge]; // only ITS sa
 647         if (h){
 648           Int_t nbins = h->GetNbinsX();
 649           Int_t nbinsy=hCorrFlukakaon[icharge]->GetNbinsY();
 650           for(Int_t ibin = 0; ibin < nbins; ibin++){
 651             Float_t pt = h->GetBinCenter(ibin);
 652             Float_t minPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(1);
 653             Float_t maxPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(nbinsy+1);
 654             if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 655             if (pt > maxPtCorrection) pt = maxPtCorrection;
 656             Float_t correction = hCorrFlukakaon[icharge]->GetBinContent(hCorrFlukakaon[icharge]->GetXaxis()->FindBin(pt));
 657             if (correction != 0) {// If the bin is empty this is a  0
 658               h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 659               h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 660             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 661               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
 662               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 663             }
 664           }
 665         }
 666       }
 667     }
 668
 669     // PROTONS
 670
 671     // ITS specific file for protons/antiprotons
 672     TFile* fITS = new TFile ("./Files/correctionForCrossSectionITS_20100719.root");
 673     TH2D * hCorrFlukaITS[kNCharge];
 674     hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
 675     hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
 676
 677     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 678       Int_t ipart = kProton;
 679       TH1 * h = hSpectra[kITS][ipart][icharge]; // only ITS sa
 680       if (h){
 681         Int_t nbins = h->GetNbinsX();
 682         Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
 683         for(Int_t ibin = 0; ibin < nbins; ibin++){
 684           Float_t pt = h->GetBinCenter(ibin);
 685           Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
 686           Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 687           if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 688           if (pt > maxPtCorrection) pt = maxPtCorrection;
 689           Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
 690           if (correction != 0) {// If the bin is empty this is a  0
 691             h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 692             h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 693           } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 694             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
 695             cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 696           }
 697         }
 698       }
 699     }
 700
 701
 702     f = new TFile ("./Files/correctionForCrossSection_20100615.root");
 703     TH2D * hCorrFluka[kNCharge];
 704     hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
 705     hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
 706     for(Int_t idet = 0; idet < kNDet; idet++){
 707       if (idet == kITS) continue; // skip ITS sa
 708       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 709         Int_t ipart = kProton;
 710         TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
 711         if (h){
 712           Int_t nbins = h->GetNbinsX();
 713           for(Int_t ibin = 0; ibin < nbins; ibin++){
 714 //          Float_t pt = convertToMT ?
 715 //            TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
 716 //            h->GetBinCenter(ibin);
 717             Float_t pt = h->GetBinCenter(ibin);
 718             Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
 719             Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1);
 720             if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 721             if (pt > maxPtCorrection) pt = maxPtCorrection;
 722             Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
 723
 724             // already in the efficiency correction (F. Noferini)
 725             //      if (idet == kTOF && icharge == kNeg) {
 726             if (idet == kTOF) {
 727               correction = 1;
 728 //            // Apply parametrized correction computed by francesco
 729 //            // Fitted panos correction and using momentum at the outer radius of the TPC
 730
 731 //            Float_t ptav = pt; // Just to use the same name francesco uses...
 732
 733 //            // from pT constrained at P.V. (ptav) to pT TPC outer (ptTPCout)
 734 //            Float_t ptTPCout=ptav*(1-6.81059e-01*TMath::Exp(-ptav*4.20094));
 735
 736 //            // traking correction (fit to Panos)
 737 //            Float_t antiprotonEC = 1 - 0.129758 *TMath::Exp(-ptav*0.679612);
 738
 739 //            // TOF matching efficiency correction (derived from Panos one scaled for M.B.(TOF)/M.B.(TPC)).
 740 //            Float_t antiprotonEC2 = TMath::Power(1 -  0.129758*TMath::Exp(-ptTPCout*0.679612),0.07162/0.03471);
 741 //            correction = antiprotonEC * antiprotonEC2;
 742             }
 743
 744             //      cout << icharge<< " " << h->GetBinCenter(ibin) << " " << pt << " " << correction << endl;
 745             if (correction != 0) {// If the bin is empty this is a  0
 746               h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 747               h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 748 //            if (idet == kTOF) {
 749 //              cout << "CORRECTING TOF TWICE" << endl;
 750 //              h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
 751 //              h->SetBinError  (ibin,h->GetBinError  (ibin)*correction);
 752 //            }
 753             } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 754               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
 755               cout << " Bin content: " << h->GetBinContent(ibin)  << endl;
 756             }
 757
 758           }
 759
 760         }
 761       }
 762     }
 763   }
 764
 765
 766   // Set style and scale
 767   for(Int_t idet = 0; idet < kNDet; idet++){
 768     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 769       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 770         if (hSpectra[idet][ipart][icharge]){
 771           hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
 772           hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
 773           hSpectra[idet][ipart][icharge]->SetLineColor   (color[ipart] );
 774           hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
 775           hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
 776           hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
 777           if (convertToMT) {
 778             TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
 779             hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 780             hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
 781             hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
 782           }
 783 //        if (idet == kTOF || idet == kTPC) {
 784 //            hSpectra[idet][ipart][icharge]->Scale(1./236763);
 785 //        }
 786 //        if (idet == kITS){
 787 //          hSpectra[idet][ipart][icharge]->Scale(202945./236763);
 788 //        }
 789           if (scaleKaons && ipart == kKaon) {
 790             hSpectra[idet][ipart][icharge]->Scale(2.);
 791           }
 792         } else {
 793           printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
 794         }
 795       }
 796     }
 797   }
 798
 799
 800   // Create fake weights for the mean; To be update once we have syst errors
 801   // Using syste from table in paper. It would be better to have this as a function of pt.
 802   TH1F * hWeights[3][kNPart];
 803   const Double_t kWeights[3][kNPart] =
 804     {{4,  3,  10.2},  // TPC
 805      {4.1,8.8,7.0 },  //TOF
 806      {4.5,5.6,7.0 }}; // ITS
 807     // {{0.1,0.1,0.1},  // TPC
 808     //  {0.1,0.1,0.1},  //TOF
 809     //  {0.1,0.1,0.1}}; // ITS
 810   for(Int_t ipart = 0; ipart <= kNPart ; ipart++){
 811     for(Int_t idet = 0; idet <= kITS ; idet++){
 812       hWeights[idet][ipart] = (TH1F*) hSpectra[idet][ipart][kPos]->Clone();
 813       Int_t nbin = hWeights[idet][ipart]->GetNbinsX();
 814       for(Int_t ibin = 1; ibin <= nbin; ibin++){
 815         hWeights[idet][ipart]->SetBinError(ibin, kWeights[idet][ipart]);
 816       }
 817     }
 818   }
 819   const Double_t scaleDet[] = {0.98,1,0.98}; // Scaling factor for the different detectors. Estimated from ratios, it has an estimated uncertainty of ~2%
 820   //  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%
 821
 822
 823   // Combine detectors
 824   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 825     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 826       TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
 827
 828       if (convertToMT) {
 829         TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
 830         htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 831
 832       }
 833       hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
 834                                                                         hSpectra[kTPC][ipart][icharge],
 835                                                                         htemplLocal,1.);;
 836       hSpectra[kCombAll][ipart][icharge]   =
 837         AliBWTools::Combine3HistosWithErrors(hSpectra[kITS][ipart][icharge],
 838                                              hSpectra[kTPC][ipart][icharge],
 839                                              hSpectra[kTOF][ipart][icharge],
 840                                              hWeights[kITS][ipart],
 841                                              hWeights[kTPC][ipart],
 842                                              hWeights[kTOF][ipart],
 843                                              htemplLocal,1,
 844                                              scaleDet[kITS],
 845                                              scaleDet[kTPC],
 846                                              scaleDet[kTOF]
 847                                              );
 848 //       if (convertToMT) {
 849 //      TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
 850 //      hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 851 //      hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
 852 //      hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
 853 //       }
 854     }
 855   }
 856
 857
 858   // Scale for the number of inelastic collisions and correct for
 859   // efficiency losses due to physics selection:
 860
 861   Double_t effPhysSel[kNPart];
 862   effPhysSel[kPion]   = 1.012;
 863   effPhysSel[kKaon]   = 1.013;
 864   effPhysSel[kProton] = 1.014;
 865
 866
 867   for(Int_t idet = 0; idet < kNHist; idet++){
 868     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 869       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 870         if(hSpectra[idet][ipart][icharge]) {
 871           //      cout << "Scaling!" << endl;
 872           if(idet!=kK0){
 873             hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel tutti? // FIXME
 874           }
 875         }
 876       }
 877     }
 878   }
 879
 880
 881 }
 882
 883 void LoadMC() {
 884
 885   TFile * f = 0;
 886   const char * evClass= "INEL";
 887   const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
 888                           "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
 889                           "./Files/dndeta_CMSD6T109_10M_900GeV.root",
 890                           "./Files/dndeta_Perugia0320_10M_900GeV.root", };
 891
 892   // Phojet
 893   for(Int_t itune = 0; itune < kNTunes; itune++){
 894     f = new TFile(files[itune]);
 895     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 896       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 897         hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
 898       }
 899     }
 900   }
 901
 902
 903   // Set style
 904   for(Int_t itune = 0; itune < kNTunes; itune++){
 905     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 906       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 907         if (hSpectraMC[itune][ipart][icharge]){
 908           hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
 909           hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
 910           hSpectraMC[itune][ipart][icharge]->SetLineColor   (mcLineColor[itune] );
 911           hSpectraMC[itune][ipart][icharge]->SetLineStyle   (mcLineStyle[itune] );
 912           hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
 913           if (convertToMT) {
 914             TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
 915             hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
 916             hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
 917             hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
 918           }
 919
 920         } else {
 921           printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
 922         }
 923       }
 924     }
 925   }
 926
 927 }
 928
 929 void DrawStar(Int_t icharge) {
 930   //  cout << icharge << endl;
 931
 932   //  gROOT->LoadMacro("StarPPSpectra.C");
 933   TGraphErrors ** gStar = StarPPSpectra();
 934
 935   for(Int_t istar = 0; istar < 6; istar++){
 936     gStar[istar]->SetMarkerStyle(kOpenStar);
 937     if      (icharge==kPos &&  (istar%2) ) gStar[istar]->Draw("P");
 938     else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
 939     else cout << "Skipping Star " << istar << endl;
 940   }
 941 }
 942
 943 void GetITSResiduals() {
 944
 945
 946   for(Int_t ipart = 0; ipart < kNPart; ipart++){
 947     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 948       //      cout << "1 " << ipart << " " << icharge << endl;
 949
 950 //       gSystem->ProcessEvents();
 951 //       gSystem->Sleep(1000);
 952       TF1* f = (TF1*)   hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
 953       TH1F * hres1, *hres2;
 954       AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
 955
 956       TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
 957       c1->SetLogy();
 958       hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
 959       hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
 960       hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
 961       hSpectra[kITS][ipart][icharge]->Draw("same");
 962       hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
 963       TLegend* l = new TLegend(    0.182886,    0.192308,    0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
 964       l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
 965       l->AddEntry(hSpectra[kITS][ipart][icharge],        "ITS");
 966       l->AddEntry(f,        "Fit to TOF+TPC");
 967       l->Draw();
 968
 969       TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
 970                                  TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
 971       c2->SetGridy();
 972       hres2->SetMinimum(-1);
 973       hres2->SetMaximum(1);
 974       hres2->Draw();
 975       hres2->GetYaxis()->SetTitleOffset(1.2);
 976       Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
 977       TLine * line = new TLine(x,-1,x,1);
 978       line->SetLineStyle(kDashed);
 979       line->Draw("same");
 980
 981       if (doPrint) {
 982         c1->Update();
 983         c2->Update();
 984         gSystem->ProcessEvents();
 985         c1->Print(TString(c1->GetName()) + ".png");
 986         c2->Print(TString(c2->GetName()) + ".png");
 987       }
 988     }
 989   }
 990 }
 991
 992 void DrawWithModels() {
 993
 994   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 995
 996     // Draw with models
 997     for(Int_t ipart = 0; ipart < kNPart; ipart++){
 998       // Multipad canvas
 999       TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],700,700);
1000       TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlag[icharge], "p1", 0.0, 0.35, 1.0,  0.95, 0, 0, 0);
1001       p1->SetBottomMargin(0);
1002       p1->Draw();
1003
1004       TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlag[icharge], "p2", 0.0, 0.05, 1.0,  0.35, 0, 0, 0);
1005       p2->SetTopMargin(0);
1006       p2->SetBottomMargin(0.3);
1007       p2->Draw();
1008
1009       Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
1010
1011       // Draw spectra
1012       p1->cd();
1013       p1->SetLogy();
1014       TH2F * hempty = new TH2F(TString("hempty")+long(icharge)+long(ipart),"hempty",100,0.,4, 100, 0.0015,5);
1015       hempty->SetXTitle("p_{t} (GeV/c)");
1016       hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1017       hempty->Draw();
1018       c1->SetLogy();
1019
1020
1021       TLegend * l =new TLegend(       0.543624,  0.431818,  0.892617,0.629371);
1022       l->SetFillColor(kWhite);
1023       hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
1024       l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data: ")+partLabel[ipart][icharge]);
1025       for(Int_t itune = 0; itune < kNTunes; itune++){
1026         l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1027         hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
1028         hSpectraMC[itune][ipart][icharge]->Draw("same,chist");
1029       }
1030       l->Draw("same");
1031
1032       // Draw ratio
1033       p2->cd();
1034       TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(icharge)+long(ipart),"hempty",100,0.,4, 100, 0.01,2.99);
1035       hemptyr->SetXTitle("p_{t} (GeV/c)");
1036       hemptyr->SetYTitle("Data/MC");
1037       hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
1038       hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
1039       hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
1040       hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
1041       hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
1042       hemptyr->SetTickLength(0.03*scaleFonts, "X");
1043       hemptyr->SetTickLength(0.02*scaleFonts, "Y");
1044       //      hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
1045       hemptyr->GetYaxis()->SetNdivisions(5);
1046       hemptyr->Draw("");
1047
1048       AliBWFunc fm;
1049       for(Int_t itune = 0; itune < kNTunes; itune++){
1050         TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
1051
1052         //      l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1053         TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[iCombInStudy][ipart][icharge],f);
1054         hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
1055         hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
1056         hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
1057         hRatio->Draw("lhist,same");
1058       }
1059
1060
1061       // print
1062       if(doPrint) {
1063         c1->Update();
1064         gSystem->ProcessEvents();
1065         c1->Print(TString(c1->GetName())+".eps");
1066         ALICEWorkInProgress(c1,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
1067         c1->Print(TString(c1->GetName())+".png");
1068
1069       }
1070     }
1071   }
1072
1073
1074
1075 }
1076
1077 void DrawAllAndKaons() {
1078
1079
1080   //  gROOT->LoadMacro("ALICEWorkInProgress.C");
1081
1082   //  gStyle->SetOptFit(0);
1083   gStyle->SetStatX(0.9);
1084   gStyle->SetStatY(0.88);
1085   gStyle->SetStatW(0.4);
1086   gStyle->SetStatH(0.1);
1087
1088   TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1089   hKaonsAllTPCTOF->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1090
1091   TH1F * hK0Scaled    = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
1092   hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
1093
1094   hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
1095   hSpectra[kKinks][kKaon][kPos]->SetLineColor(4);
1096   hSpectra[kKinks][kKaon][kPos]->SetStats(0);
1097   TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
1098   hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
1099
1100   TCanvas * c1 = new TCanvas("cKaons","cKaons",700,700);
1101   c1->SetLogy();
1102   TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-3,6);
1103   hempty->SetXTitle("p_{t} (GeV/c)");
1104   hempty->SetYTitle("dN / dp_{t} (A.U.)");
1105   hempty->Draw();
1106   hK0Scaled->Draw("same");
1107   hKaonsAllTPCTOF->Draw("same");
1108   hKinksAll->Draw("same");
1109
1110   TLegend * leg = new TLegend(0.2013423,0.2255245,0.5503356,0.4335664,NULL,"brNDC");
1111   //    leg->SetBorderSize(0);
1112 //    leg->SetLineColor(1);
1113 //    leg->SetLineStyle(1);
1114 //    leg->SetLineWidth(1);
1115 //    leg->SetFillColor(19);
1116     leg->SetFillColor(0);
1117    TLegendEntry *entry=leg->AddEntry(hKaonsAllTPCTOF,"K^{+} + K^{-}, ITS+TPC+TOF ","lpf");
1118    entry=leg->AddEntry(hK0Scaled,"K^{0} #times 2","lpf");
1119    entry=leg->AddEntry(hKinksAll,"K^{+} + K ^{-}, Kinks","lpf");
1120    leg->Draw();
1121
1122    ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Prelimiary}{Statistical Error Only}");
1123    TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
1124    tex->SetTextColor(2);
1125    tex->SetTextFont(42);
1126    tex->SetTextSize(0.03496503);
1127    tex->Draw();
1128
1129    c1->Update();
1130    if(doPrint) c1->Print(TString(c1->GetName())+".png");
1131
1132   // Draw all "stable" hadrons
1133   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1134     TCanvas * c1 = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],700,700);
1135     c1->SetLogy();
1136     c1->SetLeftMargin(0.14);
1137     TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-4,10);
1138     hempty->SetXTitle("p_{t} (GeV/c)");
1139     hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1140     hempty->GetYaxis()->SetTitleOffset(1.35);
1141     hempty->GetXaxis()->SetTitleOffset(1.1);
1142     hempty->Draw();
1143     leg = new TLegend(  0.645973,  0.2,  0.892617,0.636364, NULL,"brNDC");
1144     leg->SetFillColor(0);
1145     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1146       for(Int_t idet = 0; idet <= kITSTPC; idet++){
1147         //      if (idet == kITS) continue;
1148         //      if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
1149         hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
1150         hSpectra[idet][ipart][icharge]->Draw("same");
1151         leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet]  + ")","lpf");
1152       }
1153       //      leg->AddLine();
1154     }
1155     leg->Draw();
1156     ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1157     c1->Update();
1158     if(doPrint) c1->Print(TString(c1->GetName())+".png");
1159   }
1160
1161
1162   //  Draw ratios
1163
1164   // K-/K+ in the different detectors
1165   TCanvas * cpm=new TCanvas("cpm","Kminus/Kplus",700,700);
1166   cpm->Divide(2,2);
1167   cpm->cd(1);
1168   TH1F* hRatioKPKM_TPC=new TH1F(*(hSpectra[kTPC][kKaon][kNeg]));
1169   hRatioKPKM_TPC->SetMinimum(0.5);
1170   hRatioKPKM_TPC->SetMaximum(1.5);
1171   hRatioKPKM_TPC->Divide(hSpectra[kTPC][kKaon][kPos]);
1172   hRatioKPKM_TPC->GetYaxis()->SetTitle("K-/K+ (TPC)");
1173   hRatioKPKM_TPC->Draw();
1174   cpm->cd(2);
1175   TH1F* hRatioKPKM_ITS=new TH1F(*(hSpectra[kITS][kKaon][kNeg]));
1176   hRatioKPKM_ITS->Divide(hSpectra[kITS][kKaon][kPos]);
1177   hRatioKPKM_ITS->SetMinimum(0.5);
1178   hRatioKPKM_ITS->SetMaximum(1.5);
1179   hRatioKPKM_ITS->GetYaxis()->SetTitle("K-/K+ (ITSsa)");
1180   hRatioKPKM_ITS->Draw("");
1181   cpm->cd(3);
1182   TH1F* hRatioKPKM_TOF=new TH1F(*(hSpectra[kTOF][kKaon][kNeg]));
1183   hRatioKPKM_TOF->Divide(hSpectra[kTOF][kKaon][kPos]);
1184   hRatioKPKM_TOF->SetMinimum(0.5);
1185   hRatioKPKM_TOF->SetMaximum(1.5);
1186   hRatioKPKM_TOF->GetYaxis()->SetTitle("K-/K+ (TOF)");
1187   hRatioKPKM_TOF->Draw("");
1188   cpm->cd(4);
1189   TH1F* hRatioKPKM_ITSTPC=new TH1F(*(hSpectra[kITSTPC][kKaon][kNeg]));
1190   hRatioKPKM_ITSTPC->Divide(hSpectra[kITSTPC][kKaon][kPos]);
1191   hRatioKPKM_ITSTPC->SetMinimum(0.5);
1192   hRatioKPKM_ITSTPC->SetMaximum(1.5);
1193   hRatioKPKM_ITSTPC->GetYaxis()->SetTitle("K-/K+ (ITS Global)");
1194   hRatioKPKM_ITSTPC->Draw("");
1195
1196
1197   // ITS/TPC
1198   TH1F * hRatioITSTPC[kNPart][kNCharge];
1199   for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1200     // Create canvas
1201     TCanvas * c1 = new TCanvas(TString("cITSTPCRatio_")+chargeFlag[icharge],TString("cITSTPCRatio_")+chargeFlag[icharge],1200,500);
1202     c1->Divide(3,1);
1203     c1->SetGridy();
1204     TH2F * hempty = new TH2F(TString("hemptyR")+long(icharge),"ITSsa/TPC ",100,0.,1., 100, 0.5,1.5);
1205     hempty->SetXTitle("p_{t} (GeV/c)");
1206     hempty->SetYTitle("ITSsa / TPC");
1207     // Loop over particles
1208     for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1209       // Clone histo
1210       hRatioITSTPC[ipart][icharge]=new TH1F(*hSpectra[kITS][ipart][icharge]);
1211       Int_t nBinsITS=hSpectra[kITS][ipart][icharge]->GetNbinsX();
1212       Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1213       // Loop over ITS bins,
1214       for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1215         hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1216         hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1217         Float_t lowPtITS=hSpectra[kITS][ipart][icharge]->GetBinLowEdge(iBin);
1218         Float_t binWidITS=hSpectra[kITS][ipart][icharge]->GetBinWidth(iBin);
1219         // Loop over TPC bins and find overlapping bins to ITS
1220         for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1221           Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1222           Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1223           if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1224             Float_t numer=hSpectra[kITS][ipart][icharge]->GetBinContent(iBin);
1225             Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1226             Float_t enumer=hSpectra[kITS][ipart][icharge]->GetBinError(iBin);
1227             Float_t edenom=hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1228             Double_t ratio=0.;
1229             Double_t eratio=0.;
1230             if(numer>0. && denom>0.){
1231               ratio=numer/denom;
1232               eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1233             }
1234             hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1235             hRatioITSTPC[ipart][icharge]->SetBinError(iBin,eratio);
1236             break;
1237           }
1238         }
1239       }
1240       c1->cd(ipart+1);
1241       // hempty->SetStats(1);
1242       // hempty->Draw();
1243       hRatioITSTPC[ipart][icharge]->SetStats(1);
1244       hRatioITSTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1245       hRatioITSTPC[ipart][icharge]->Draw("");
1246       hRatioITSTPC[ipart][icharge]->Fit("pol0","","same");
1247
1248     }
1249     if(doPrint) c1->Print(TString(c1->GetName())+".png");
1250   }
1251
1252   // TOF/TPC
1253   TH1F * hRatioTOFTPC[kNPart][kNCharge];
1254   for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1255     // create canvas
1256     TCanvas * c1t = new TCanvas(TString("cTOFTPCRatio_")+chargeFlag[icharge],TString("cTOFTPCRatio_")+chargeFlag[icharge],1200,500);
1257     c1t->SetGridy();
1258     c1t->Divide(3,1);
1259     TH2F * hemptyt = new TH2F(TString("hemptyRt")+long(icharge),"TOF/TPC ",100,0.,1., 100, 0.5,1.5);
1260     hemptyt->SetXTitle("p_{t} (GeV/c)");
1261     hemptyt->SetYTitle("TOF / TPC");
1262     //    hemptyt->Draw();
1263     for(Int_t ipart = 0; ipart < kNPart; ipart++) { // loop over particles
1264       // Clone histo
1265       hRatioTOFTPC[ipart][icharge]=new TH1F(*hSpectra[kTOF][ipart][icharge]);
1266       Int_t nBinsTOF=hSpectra[kTOF][ipart][icharge]->GetNbinsX();
1267       Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1268       // Loop over TOF bins
1269       for(Int_t iBin=1; iBin<=nBinsTOF; iBin++){
1270         hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1271         hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1272         Float_t lowPtTOF=hSpectra[kTOF][ipart][icharge]->GetBinLowEdge(iBin);
1273         Float_t binWidTOF=hSpectra[kTOF][ipart][icharge]->GetBinWidth(iBin);
1274         // Loop over TPC bins and find overlapping bins to ITS
1275         for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1276           Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1277           Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1278           if(TMath::Abs(lowPtTOF-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidTOF)<0.001){
1279             Float_t numer=hSpectra[kTOF][ipart][icharge]->GetBinContent(iBin);
1280             Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1281             Float_t enumer=hSpectra[kTOF][ipart][icharge]->GetBinError(iBin);
1282             Float_t edenom=hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1283             Double_t ratio=0.;
1284             Double_t eratio=0.;
1285             if(numer>0. && denom>0.){
1286               ratio=numer/denom;
1287               eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1288             }
1289             hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1290             hRatioTOFTPC[ipart][icharge]->SetBinError(iBin,eratio);
1291             break;
1292           }
1293         }
1294       }
1295       c1t->cd(ipart+1);
1296       hRatioITSTPC[ipart][icharge]->SetStats(1);
1297       hRatioITSTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1298       hRatioITSTPC[ipart][icharge]->Draw("");
1299       hRatioITSTPC[ipart][icharge]->Fit("pol0","","same");
1300     }
1301     if(doPrint) c1t->Print(TString(c1t->GetName())+".png");
1302   }
1303
1304 }
1305
1306 void DrawWithJacek() {
1307
1308   //1. Convert spectra to dNdeta and sum
1309   TH1F * hsum = (TH1F*) htemplate->Clone();
1310   hsum->Reset();
1311   Int_t idet= iCombInStudy;
1312   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1313     for(Int_t ipart = 0; ipart < kNPart; ipart++){
1314       TH1 * h = hSpectra[idet][ipart][icharge];
1315       Int_t nbin = h->GetNbinsX();
1316       for(Int_t ibin = 1; ibin <= nbin; ibin++){
1317         Double_t pt = h->GetBinCenter(ibin);
1318         Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
1319         Double_t jacobian = pt/mt;
1320         h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
1321         h->SetBinError  (ibin,h->GetBinError  (ibin)*jacobian);
1322         Int_t ibinSum = hsum->FindBin(pt);
1323         Double_t epsilon = 0.0001;
1324         if ( h->GetBinContent(ibin) > 0 &&
1325              (TMath::Abs(h->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1326               TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1327              ) {
1328           cout << "DISCREPANCY IN BIN RANGES" << endl;
1329           cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
1330                << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
1331                << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
1332           cout << "" << endl;
1333         }
1334
1335
1336         hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
1337         hsum->SetBinError  (ibinSum,0);
1338       }
1339       //        hsum->Add(h);
1340     }
1341   }
1342
1343
1344   // Load Jacek and Draw both:
1345 //   new TFile ("./Files/dNdPt_Data_Points_ALICE_900GeV.root");
1346 //   TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("inel");
1347 //   gJacek->Draw("AP");
1348 //   hsum->Draw("same");
1349
1350 //   TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
1351
1352
1353 //   new TCanvas();
1354 //   gRatio->Draw("AP");
1355
1356
1357
1358 }
1359
1360
1361 void DrawRatioToStar() {
1362
1363   // Star data
1364   //  gROOT->LoadMacro("StarPPSpectra.C");
1365   TGraphErrors ** gStar = StarPPSpectra();
1366   gStar[0]->SetMarkerStyle(kOpenStar);
1367   gStar[1]->SetMarkerStyle(kFullStar);
1368   gStar[2]->SetMarkerStyle(kOpenStar);
1369   gStar[3]->SetMarkerStyle(kFullStar);
1370   gStar[4]->SetMarkerStyle(kOpenStar);
1371   gStar[5]->SetMarkerStyle(kFullStar);
1372
1373   // ALICE, INEL -> NSD
1374   Double_t scaleYield = 3.58/3.02; // from paper 2
1375   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1376     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1377       hSpectra[iCombInStudy][ipart][icharge]->Scale(scaleYield);
1378     }
1379   }
1380
1381
1382   TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
1383   TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
1384   hempty->SetXTitle("p_{t} (GeV/c)");
1385   hempty->SetYTitle("ALICE/STAR (NSD)");
1386   hempty->Draw();
1387
1388   TCanvas * c1Comp = new TCanvas("cCompToStarNeg","cCompToStarNeg");
1389   c1Comp->SetLogy();
1390   TH2F * hempty2 = new TH2F(TString("hemptyCompNeg"),"hemptyCompNeg",100,0.,1.5, 100, 0.001,10);
1391   hempty2->SetXTitle("p_{t} (GeV/c)");
1392   hempty2->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1} [NSD]");
1393   hempty2->Draw();
1394
1395   TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
1396   leg->SetBorderSize(0);
1397   leg->SetLineColor(1);
1398   leg->SetLineStyle(1);
1399   leg->SetLineWidth(1);
1400   leg->SetFillColor(0);
1401   leg->SetFillStyle(1001);
1402
1403
1404
1405   c1->cd();
1406   TGraphErrors * g ;
1407   g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[iCombInStudy][kPion][kNeg],1);
1408   g->SetMarkerStyle(kFullCircle);
1409   g->SetMarkerColor(kBlack);
1410   g->Draw("p");
1411   leg->AddEntry(g,"#pi^{-}","lp");
1412   g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[iCombInStudy][kKaon][kNeg],1);
1413   g->SetMarkerStyle(kOpenCircle);
1414   g->SetMarkerColor(kRed);
1415   g->Draw("p");
1416   leg->AddEntry(g,"K^{-}","lp");
1417   g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[iCombInStudy][kProton][kNeg],1);
1418   g->SetMarkerStyle(kOpenSquare);
1419   g->SetMarkerColor(kBlue);
1420   g->Draw("p");
1421   leg->AddEntry(g,"#bar{p}","lp");
1422   leg->Draw();
1423
1424   c1Comp->cd();
1425   gStar[0]->Draw("p");
1426   hSpectra[iCombInStudy][kPion][kNeg]->Draw("same");
1427   gStar[2]->Draw("p");
1428   hSpectra[iCombInStudy][kKaon][kNeg]->Draw("same");
1429   gStar[4]->Draw("p");
1430   hSpectra[iCombInStudy][kProton][kNeg]->Draw("same");
1431
1432
1433
1434   TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
1435   hempty->Draw();
1436   leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
1437   leg->SetBorderSize(0);
1438   leg->SetLineColor(1);
1439   leg->SetLineStyle(1);
1440   leg->SetLineWidth(1);
1441   leg->SetFillColor(0);
1442   leg->SetFillStyle(1001);
1443
1444   TCanvas * c2Comp = new TCanvas("cCompToStarPos","cCompToStarPos");
1445   c2Comp->SetLogy();
1446   hempty2->Draw();
1447
1448   c2->cd();
1449  //  TGraphErrors * g ;
1450   g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[iCombInStudy][kPion][kPos],1);
1451   g->SetMarkerStyle(kFullCircle);
1452   g->SetMarkerColor(kBlack);
1453   g->Draw("p");
1454   leg->AddEntry(g,"#pi^{+}","lp");
1455   g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[iCombInStudy][kKaon][kPos],1);
1456   g->SetMarkerStyle(kOpenCircle);
1457   g->SetMarkerColor(kRed);
1458   g->Draw("p");
1459   leg->AddEntry(g,"K^{+}","lp");
1460   g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[iCombInStudy][kProton][kPos],1);
1461   g->SetMarkerStyle(kOpenSquare);
1462   g->SetMarkerColor(kBlue);
1463   g->Draw("p");
1464   leg->AddEntry(g,"p","lp");
1465   leg->Draw();
1466
1467
1468   c2Comp->cd();
1469   gStar[1]->Draw("p");
1470   hSpectra[iCombInStudy][kPion][kPos]->Draw("same");
1471   gStar[3]->Draw("p");
1472   hSpectra[iCombInStudy][kKaon][kPos]->Draw("same");
1473   gStar[5]->Draw("p");
1474   hSpectra[iCombInStudy][kProton][kPos]->Draw("same");
1475
1476
1477   c1->Update();
1478   c2->Update();
1479   gSystem->ProcessEvents();
1480   c1->Print(TString(c1->GetName()) + ".eps");
1481   c2->Print(TString(c2->GetName()) + ".eps");
1482   ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1483   ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1484   c1->Update();
1485   c2->Update();
1486   c1->Print(TString(c1->GetName()) + ".png");
1487   c2->Print(TString(c2->GetName()) + ".png");
1488
1489
1490
1491
1492 }
1493
1494
1495
1496 void DrawRatios() {
1497
1498   // Draws ratios of combined spectra
1499
1500   // Compute ratios
1501   TH1F * hPosNegRatio[kNPart];
1502
1503   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1504     hPosNegRatio[ipart] = (TH1F*) hSpectra[iCombInStudy][ipart][kPos]->Clone();
1505     hPosNegRatio[ipart]->Divide(hSpectra[iCombInStudy][ipart][kNeg]);
1506     hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
1507     hPosNegRatio[ipart]->SetMinimum(0.5);
1508     hPosNegRatio[ipart]->SetMaximum(1.5);
1509   }
1510
1511   TH1F * hKPiRatio = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1512   hKPiRatio->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1513   TH1F * htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1514   htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1515   hKPiRatio->Divide(htmp);
1516   hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1517
1518   TH1F * hKPiRatioMC[kNTunes];
1519   if(showMC){
1520     for(Int_t itune = 0; itune < kNTunes; itune++){
1521       hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
1522       hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
1523       TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1524       htmp->Add(hSpectraMC[itune][kPion][kNeg]);
1525       hKPiRatioMC[itune]->Divide(htmp);
1526       hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1527     }
1528   }
1529
1530
1531
1532   TH1F * hPPiRatio = (TH1F*) hSpectra[iCombInStudy][kProton][kPos]->Clone();
1533   hPPiRatio->Add(hSpectra[iCombInStudy][kProton][kNeg]);
1534   htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1535   htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1536   hPPiRatio->Divide(htmp);
1537   hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1538
1539   if(showMC){
1540     TH1F * hPPiRatioMC[kNTunes];
1541     for(Int_t itune = 0; itune < kNTunes; itune++){
1542       hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
1543       hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
1544       TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1545       htmp->Add(hSpectraMC[itune][kPion][kNeg]);
1546       hPPiRatioMC[itune]->Divide(htmp);
1547       hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1548     }
1549   }
1550
1551   // Draw
1552 //   TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
1553 //   hempty->SetXTitle("p_{t} (GeV/c)");
1554   //  hempty->SetYTitle("");
1555
1556   // tmp: overlay levi fits
1557   AliBWFunc * fm2 = new AliBWFunc;
1558   fm2->SetVarType(AliBWFunc::kdNdpt);
1559   TF1 * fLevi[kNPart]  [kNCharge];
1560   fLevi[kPion]  [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
1561   fLevi[kKaon]  [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
1562   fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
1563   fLevi[kPion]  [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
1564   fLevi[kKaon]  [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
1565   fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
1566   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1567     for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1568       fLevi[ipart][icharge]->SetRange(0,4);
1569     }
1570   }
1571
1572
1573   for(Int_t ipart = 0; ipart < kNPart; ipart++){
1574     TString detName = detFlag[iCombInStudy];
1575     detName.ReplaceAll(" ", "_");
1576     detName.ReplaceAll("+", "");
1577
1578     TCanvas * c1 = new TCanvas(TString("cRatio_")+detName+TString("_")+partFlag[ipart], TString("cRatio_")+detName+partFlag[ipart]);
1579     c1->SetGridy();
1580     hPosNegRatio[ipart]->Draw();
1581     TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
1582     //    fRatio->Draw("same");
1583     fRatio->SetRange(0,5);
1584     if (doPrint) {
1585       c1->Update();
1586       gSystem->ProcessEvents();
1587       c1->Print(TString(c1->GetName()) + ".png");
1588     }
1589
1590   }
1591
1592
1593   TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
1594   c2->SetGridy();
1595   hKPiRatio->Draw();
1596   TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1597   lMC->SetFillColor(kWhite);
1598
1599   if(showE735){
1600     gROOT->LoadMacro("GetE735Ratios.C");
1601     GetE735Ratios(0,0)->Draw("EX0,same");
1602     GetE735Ratios(0,1)->Draw("EX0,same");
1603     GetE735Ratios(0,2)->Draw("EX0,same");
1604     GetE735Ratios(0,3)->Draw("EX0,same");
1605   }
1606   hKPiRatio->SetMarkerStyle(20);
1607   hKPiRatio->Draw("same");
1608
1609   if(showMC){
1610     for(Int_t itune = 0; itune < kNTunes; itune++){
1611       lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1612       hKPiRatioMC[itune]->SetLineWidth(2);
1613       hKPiRatioMC[itune]->Draw("same,chist");
1614     }
1615
1616     lMC->Draw();
1617   }
1618
1619   if(showE735){
1620     TLegend * l = new TLegend(  0.1879,  0.68,  0.54,0.92);
1621     l->SetFillColor(kWhite);
1622     l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1623     l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
1624     l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
1625     l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
1626     l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
1627     l->Draw();
1628   }
1629
1630
1631   TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
1632   c3->SetGridy();
1633   hPPiRatio->Draw();
1634   hPPiRatio->SetMaximum(0.39);
1635   if(showMC){
1636     lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1637     lMC->SetFillColor(kWhite);
1638
1639     for(Int_t itune = 0; itune < kNTunes; itune++){
1640       lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1641       hKPiRatioMC[itune]->SetLineWidth(2);
1642       hKPiRatioMC[itune]->Draw("same,chist");
1643     }
1644
1645     lMC->Draw();
1646   }
1647
1648   if (doPrint) {
1649     c2->Update();
1650     gSystem->ProcessEvents();
1651     c2->Print(TString(c2->GetName()) + ".png");
1652     c2->Print(TString(c2->GetName()) + ".eps");
1653     c3->Update();
1654     gSystem->ProcessEvents();
1655     c3->Print(TString(c3->GetName()) + ".png");
1656     c3->Print(TString(c3->GetName()) + ".eps");
1657   }
1658
1659
1660 }
1661
1662
1663 void Help() {
1664
1665   cout << "Macro: void CombineSpectra(Int_t analysisType=kDoFits, Int_t  fitFuncID = kFitLevi) " << endl;
1666   cout << "" << endl;
1667
1668   cout << "Possible Arguments" << endl;
1669   cout << "- analysisType:" << endl;
1670   cout << "    kDoFits:           Fit Combined Spectra " << endl;
1671   cout << "    kDoRatios:         Particle ratios K/pi and p/pi" << endl;
1672   cout << "    kDoSuperposition:  Compare different detectors (superimpose and ratios)" << endl;
1673   cout << "    kDoCompareStar:    Compare combined spectra to star results" << endl;
1674   cout << "    kDoDrawWithModels: Compare combined spectra and models" << endl;
1675   cout << "    kDoHelp:           This help" << endl;
1676   cout << "- fitFuncID, function used to extrapolate and compute yields" << endl;
1677   cout << "    An analitic fit function [kFitLevi, kFitUA1, kFitPowerLaw]" << endl;
1678   cout << "    Or a shape from a MC moder [kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0]" << endl;
1679   cout << "    Which is fitted to the data at low pt and used to extrapolate at low pt" << endl;
1680
1681
1682 }