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