1 #if !defined(__CINT__) || defined(__MAKECINT__)
4 #include "TGraphErrors.h"
6 #include "AliBWTools.h"
9 #include "TDatabasePDG.h"
14 #include "AliLatexTable.h"
16 #include "TVirtualFitter.h"
25 #include "TPaveText.h"
26 #include "StarPPSpectra.C"
27 #include "GetE735Ratios.C"
29 #include "TObjString.h"
30 #include "TGraphAsymmErrors.h"
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,
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?
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)"},
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",
66 "Pythia - Perugia0 - 320", };
67 const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
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};
75 const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
76 const Int_t mcLineStyle[] = {1,2,3,4};
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};
83 TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
86 TH1F * hSpectra[kNHist][kNPart][kNCharge];
87 TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
88 TH1F * hSystError[kNHist][kNPart][kNCharge];
89 Double_t mass[kNPart];
96 // Additional tasks (may be uncommented below)
98 void DrawStar(Int_t icharge);
99 void GetITSResiduals();
100 void DrawWithModels() ;
101 void DrawAllAndKaons();
102 void DrawWithJacek();
103 void DrawRatioToStar();
107 void PrintCanvas(TCanvas* c,const TString formats) ;
111 void ALICEWorkInProgress(TCanvas *c,TString today="11/05/2010", TString label = "ALICE performance"){
113 TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.72,0.89,0.89);
114 //TPad *myPadLogo = new TPad("myPadLogo", "Pad for ALICE Logo",0.72,0.62,0.85,0.75);
115 myPadLogo->SetFillColor(0);
116 myPadLogo->SetBorderMode(0);
117 myPadLogo->SetBorderSize(2);
118 myPadLogo->SetFrameBorderMode(0);
119 myPadLogo->SetLeftMargin(0.0);
120 myPadLogo->SetTopMargin(0.0);
121 myPadLogo->SetBottomMargin(0.0);
122 myPadLogo->SetRightMargin(0.0);
123 myPadLogo->SetFillStyle(0);
126 TASImage *myAliceLogo = new TASImage("alice_logo.png");
129 TPaveText* t1=new TPaveText(0.418103, 0.837798, 0.656609, 0.888393,"NDC");
131 t1->SetBorderSize(0);
132 t1->AddText(0.,0.,label);
133 t1->SetTextColor(kRed);
135 t1->SetTextSize(0.04);
137 TPaveText* t2=new TPaveText(0.418103, 0.80, 0.656609, 0.84,"NDC");
139 t2->SetBorderSize(0);
140 t2->AddText(0.,0.,"pp at #sqrt{s} = 900 GeV (2009 data)");
141 t2->SetTextColor(kRed);
143 t2->SetTextSize(0.027);
145 TPaveText* t3=new TPaveText(0.418103, 0.76, 0.656609, 0.80,"NDC");
147 t3->SetBorderSize(0);
148 t3->AddText(0.,0.,"Statistical and systematic errors");
149 t3->SetTextColor(kRed);
151 t3->SetTextSize(0.027);
153 // TPaveText* t2=new TPaveText(0.65,0.65,0.89,0.7,"NDC");
154 // t2->SetFillStyle(0);
155 // t2->SetBorderSize(0);
156 // t2->SetTextColor(kRed);
157 // t2->SetTextFont(52);
158 // t2->AddText(0.,0.,today.Data());
169 Bool_t convertToMT = 0;
170 Bool_t sumCharge = 0;
171 Int_t whatToFit = kStatErrors;
172 Bool_t scaleKaons = kFALSE;
173 Bool_t drawStar = kFALSE; // Overlay star when doing fits
174 Bool_t correctSecondaries = 1;
175 Bool_t correctGeantFlukaXS = 1;
176 Int_t iCombInStudy = kCombAll; //kCombTOFTPC
177 Int_t fitFuncID = kFitLevi;
179 Bool_t showE735=kTRUE;
180 Bool_t useSecCorrFromDCA=kTRUE;
181 Bool_t flagPreliminary=kFALSE; // Add "preliminary" flag and logo to the plots
183 TString printFormats = "C,eps,root"; // format in which canvases will be printed, if PrintCanvas is called (not all prints are based on printcanvas at the moment). This is a comma separated list.
186 void CombineSpectra(Int_t analysisType=kDoHelp, Int_t locfitFuncID = kFitLevi) { //kDoSuperposition;//kDoDrawWithModels;// kDoFits; //kDoRatios;
188 // This macro is used to combine the 900 GeV spectra from 2009
190 fitFuncID=locfitFuncID;
195 // Print Help and quit
196 if (analysisType == kDoHelp) {
203 today = today + long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
207 mass[kPion] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
208 mass[kKaon] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
209 mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
217 // GetITSResiduals();
218 if(analysisType==kDoSuperposition) DrawAllAndKaons();
219 else if(analysisType==kDoDrawWithModels) DrawWithModels() ;
220 else if(analysisType==kDoCompareToAllCharged) DrawWithJacek();
221 else if(analysisType==kDoCompareToStar) DrawRatioToStar();
222 else if(analysisType==kDoRatios) DrawRatios();
223 else if(analysisType==kDoDrawSyst) DrawWithSyst();
224 else if(analysisType==kDoFits) FitCombined();
230 gStyle->SetOptTitle(0);
231 gStyle->SetOptStat(0);
233 // Draw combined & Fit
234 AliBWFunc * fm = new AliBWFunc;
235 fm->SetVarType(AliBWFunc::kdNdpt);
236 if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
238 // table to print results
239 AliLatexTable table(10,"c|ccccccccc");
240 if (fitFuncID == kFitLevi) {
241 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)\\\\");
242 } else if (fitFuncID == kFitPowerLaw) {
243 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 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) \\\\");
249 AliLatexTable tempTable(3,"c|cc");
250 tempTable.InsertCustomRow("Part & Yield & $\\langle p_{t} \\rangle$ \\\\");
251 tempTable.InsertHline();
253 TH1F* hRatiosToFit[kNPart][kNCharge]; // Ratio data/fit function, stat error by default
254 TH1F* hRatiosToFitSyst[kNPart][kNCharge]; // Ratio data/fit, stat + syst
256 Int_t chargeLoop = sumCharge ? 1 : 2;
257 for(Int_t icharge = 0; icharge < chargeLoop; icharge++){
259 TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge],700,700);
262 c2->SetLeftMargin(0.14);
263 TCanvas * c2r = new TCanvas(TString("cCombinedRatio")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombinedRatio")+chargeFlag[icharge],700,700);
266 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
267 hempty->SetXTitle("p_{t} (GeV/c)");
268 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
269 hempty->GetYaxis()->SetTitleOffset(1.35);
270 hempty->GetXaxis()->SetTitleOffset(1.1);
274 TH2F * hemptyR = new TH2F(TString("hemptyR")+long(icharge),"hemptyR",100,0.,2.9, 100, 0.5,1.5);
275 hemptyR->SetXTitle("p_{t} (GeV/c)");
276 hemptyR->SetYTitle("Data/Fit");
280 TLegend * l = new TLegend(0.176724, 0.153274, 0.475575, 0.434524,chargeLabel[icharge]);
281 l->SetFillColor(kWhite);
282 l->SetTextSize(0.035);
286 TPaveText* tf=new TPaveText(0.321839, 0.175595, 0.685345, 0.299107,"NDC");
287 if(fitFuncID == kFitLevi){
288 tf->AddText("#frac{dN}{dp_{t}} #propto p_{t} #left(1+#frac{#sqrt{m^{2}+p_{t}^{2}} -m}{nT} #right)^{-n}");
291 tf->SetTextSize(0.032);
293 for(Int_t ipart = 0; ipart < kNPart; ipart++){
294 printf(" ----- Fit %s %s ------\n",partFlag[ipart],chargeFlag[icharge]);
301 if(fitFuncID == kFitLevi) {
302 normPar = 0; // The levi is normalized by parameter 0
304 func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
306 func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
307 if (ipart == kProton)
308 func = fm->GetLevi(mass[ipart], 0.15, 8.5, 0.09);
310 else if(fitFuncID == kFitUA1) func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
311 else if(fitFuncID == kFitPowerLaw) {
315 func = fm->GetPowerLaw(1.0, 8.6, 7);
317 func = fm->GetPowerLaw(3.0, 12, 2.6);
318 if (ipart == kProton)
319 func = fm->GetPowerLaw(24, 72, 0.8);
321 else if(fitFuncID == kFitPhojet) func = fm->GetHistoFunc(hSpectraMC[kPhojet] [ipart][icharge]);
322 else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
323 else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
324 else if(fitFuncID == kFitCMS6D6T) func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
326 cout << "Unknown fit ID " << fitFuncID << endl;
330 if(fitFuncID >= kFitPhojet){
335 // Temp: fit histo with sist errors
336 TH1F * hsyst = new TH1F(*htemplate);
337 TH1F * hsyststat = 0;
338 hsyst->SetFillColor(kGray);
339 AliBWTools::GetValueAndError(hsyst,hSpectra[iCombInStudy][ipart][icharge],hSystError[iCombInStudy][ipart][icharge],kTRUE);
340 hsyststat= new TH1F(*hsyst);
341 AliBWTools::GetHistoCombinedErrors(hsyststat,hSpectra[iCombInStudy][ipart][icharge]); // combine syst and stat
342 hsyststat->SetFillColor(kGray);
345 if (whatToFit == kStatErrors) hToFit = hSpectra[iCombInStudy][ipart][icharge]; // Shorthand
346 if (sumCharge) hToFit->Add(hSpectra[iCombInStudy][ipart][1]);
347 if (whatToFit == kStatSystErrors) {
348 hToFit = hsyststat;// Shorthand
350 if (whatToFit == kSystErrors) hToFit = hsyst;
353 if(!AliBWTools::Fit(hToFit,func,fitmin,fitmax)) {
354 cout << " FIT ERROR " << endl;
357 cout << "DRAWING" << endl;
359 hsyststat->Draw("same,e2");
360 // TH1F * hsyststat2 = new TH1F(*hsyststat);
361 // hsyststat2->SetFillColor(kWhite);
362 // hsyststat2->SetLineColor(kBlack);
363 // hsyststat2->Draw("same,e2");
364 hToFit->Draw("same");
365 TF1* fitfunc=(TF1*)hToFit->GetListOfFunctions()->At(0);
366 fitfunc->Draw("same");
367 fitfunc->SetRange(0,4);
368 fitfunc->SetLineColor(hSpectra[iCombInStudy][ipart][icharge]->GetLineColor());
369 if(drawStar) DrawStar(icharge);
370 hRatiosToFit[ipart][icharge]=(TH1F*)hToFit->Clone(Form("hRatio%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
371 hRatiosToFitSyst[ipart][icharge]=(TH1F*)hsyststat->Clone(Form("hRatioSyst%s%s",chargeFlag[icharge],partFlag[icharge])); // Ratio data/fit
372 for(Int_t iBin=1; iBin<hToFit->GetNbinsX(); iBin++){
373 Double_t lowLim=hToFit->GetBinLowEdge(iBin);
374 Double_t highLim=hToFit->GetBinLowEdge(iBin+1);
375 Double_t contFunc=fitfunc->Integral(lowLim,highLim)/(highLim-lowLim);
376 Double_t ratio=hToFit->GetBinContent(iBin)/contFunc;
377 Double_t eratio=hToFit->GetBinError(iBin)/contFunc;
378 Double_t eratioSyst=hsyststat->GetBinError(iBin)/contFunc;
379 hRatiosToFitSyst[ipart][icharge]->SetBinContent(iBin,ratio);
380 hRatiosToFitSyst[ipart][icharge]->SetBinError(iBin,eratioSyst);
381 hRatiosToFit[ipart][icharge]->SetBinContent(iBin,ratio);
382 hRatiosToFit[ipart][icharge]->SetBinError(iBin,eratio);
384 // hToFit->GetListOfFunctions()->At(0)->Draw("same");
385 // ((TF1*)hToFit->GetListOfFunctions()->At(0))->SetRange(0,4);
386 // ((TF1*)hToFit->GetListOfFunctions()->At(0))->SetLineColor(hToFit->GetLineColor());
389 scaleKaons && ipart == kKaon ?
390 (TString(partLabel[ipart][icharge])+" #times 2").Data()
391 : partLabel[ipart][icharge]);
392 // TF1 * fClone = (TF1*) hToFit->GetListOfFunctions()->At(0)->Clone();
393 // hToFit->GetListOfFunctions()->Add(fClone);
394 // fClone->SetLineStyle(kDashed);
395 // fClone->SetRange(0,100);
396 // fClone->Draw("same");
399 // Float_t yield = func->Integral(0.45,1.05);
400 // Float_t yieldE = func->IntegralError(0.45,1.05);
402 Float_t yield = func->Integral(0.,100);
403 //Float_t yieldE = func->IntegralError(0.,100);
405 Double_t yieldTools, yieldETools;
406 Double_t partialYields[3],partialYieldsErrors[3];
407 AliBWTools::GetYield(hToFit, func, yieldTools, yieldETools,
408 0, 100, partialYields,partialYieldsErrors);
409 Double_t tslope = func->GetParameter(2);
410 Double_t tslopeE = func->GetParError(2);
412 table.SetNextCol(partLatex[ipart][icharge]);
413 //table.SetNextCol(yield,yieldE,-4);
414 table.SetNextCol(yieldTools, yieldETools,-4);
415 table.SetNextCol(func->GetParameter(0));
416 table.SetNextCol(tslope,tslopeE,-4);
417 table.SetNextCol(func->GetParameter(1),func->GetParError(1),-2);
418 table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
419 Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hToFit);
420 //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
421 Float_t yieldAbove = func->Integral(lowestPoint,100);
422 table.SetNextCol(lowestPoint,-2);
423 table.SetNextCol((1 - yieldAbove/yield)*100,-1);
424 Float_t mean, meane; // Mean using only fit function
425 Double_t meanDF, meanDFe; // Mean from data + fit
426 Float_t mean2, mean2e;
427 if (fitFuncID < kFitPhojet) {
428 AliBWTools::GetMean (func, mean, meane , 0.,100., normPar);
429 AliBWTools::GetMeanSquare(func, mean2, mean2e, 0.,100., normPar);
431 // AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
432 AliBWTools::GetMeanDataAndExtrapolation(hToFit, func, meanDF, meanDFe, 0.,100.);
433 table.SetNextCol(mean, meane ,-4);
434 table.SetNextCol(meanDF, meanDFe,-4);
436 // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
441 tempTable.SetNextCol(partLatex[ipart][icharge]);
442 tempTable.SetNextCol(yieldTools, yieldETools, -4);
443 tempTable.SetNextCol(meanDF, meanDFe ,-4);
444 // tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
445 // tempTable.SetNextCol(yieldAbove/yield,-2);
446 tempTable.InsertRow();
448 hRatiosToFitSyst[ipart][icharge]->Draw("e2same");
449 hRatiosToFit[ipart][icharge]->Draw("esame");
453 if (flagPreliminary) ALICEWorkInProgress(c2,"","ALICE Preliminary");
455 if (flagPreliminary) tf->Draw();
457 if(flagPreliminary) ALICEWorkInProgress(c2r,"","ALICE Preliminary");
459 PrintCanvas(c2,printFormats);
460 PrintCanvas(c2r,printFormats);
466 table.PrintTable("ASCII");
469 tempTable.PrintTable("ASCII");
477 // Loads spectra for different detectors and corresponding systematic errors.
481 // Systematic errors: initialize histos
483 for(Int_t ipart = 0; ipart < kNPart; ipart++){
484 for(Int_t idet = 0; idet < kNHist; idet++){
485 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
486 hSystError[idet][ipart][icharge] = new TH1F (TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
487 TString("hSyst_")+detFlag[idet]+partFlag[ipart]+chargeFlag[icharge],
488 nbinsTempl,templBins);
489 hSystError[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
490 hSystError[idet][ipart][icharge]->SetLineColor (color[ipart] );
503 f = new TFile("./Files/effhistos.root");
504 TH1D * hEffTrackTOF[kNPart][kNCharge];
505 TH1D * hEffMatchTOF[kNPart][kNCharge];
506 hEffTrackTOF[kPion] [kPos] = (TH1D*) f->Get("hpitrk_pos");
507 hEffTrackTOF[kKaon] [kPos] = (TH1D*) f->Get("hkatrk_pos");
508 hEffTrackTOF[kProton][kPos] = (TH1D*) f->Get("hprtrk_pos");
509 hEffMatchTOF[kPion] [kPos] = (TH1D*) f->Get("hpieff_pos");
510 hEffMatchTOF[kKaon] [kPos] = (TH1D*) f->Get("hkaeff_pos");
511 hEffMatchTOF[kProton][kPos] = (TH1D*) f->Get("hpreff_pos");
512 hEffTrackTOF[kPion] [kNeg] = (TH1D*) f->Get("hpitrk_neg");
513 hEffTrackTOF[kKaon] [kNeg] = (TH1D*) f->Get("hkatrk_neg");
514 hEffTrackTOF[kProton][kNeg] = (TH1D*) f->Get("hprtrk_neg");
515 hEffMatchTOF[kPion] [kNeg] = (TH1D*) f->Get("hpieff_neg");
516 hEffMatchTOF[kKaon] [kNeg] = (TH1D*) f->Get("hkaeff_neg");
517 hEffMatchTOF[kProton][kNeg] = (TH1D*) f->Get("hpreff_neg");
519 // f = new TFile("./Files/spectra-pos-y.root");
520 f = new TFile("./Files/spectraRaw-pos-y.root");
521 hSpectra[kTOF][kPion] [kPos]= (TH1F*) f->Get("hpi");
522 hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
523 hSpectra[kTOF][kKaon] [kPos]= (TH1F*) f->Get("hka");
524 hSpectra[kTOF][kPion] [kPos]->SetName("hpiPos");
525 hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
526 hSpectra[kTOF][kKaon] [kPos]->SetName("hkaPos");
527 //f = new TFile("./Files/spectra-neg-y.root");
528 f = new TFile("./Files/spectraRaw-neg-y.root");
529 hSpectra[kTOF][kPion] [kNeg]= (TH1F*) f->Get("hpi");
530 hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
531 hSpectra[kTOF][kKaon] [kNeg]= (TH1F*) f->Get("hka");
532 hSpectra[kTOF][kPion] [kNeg]->SetName("hpiNeg");
533 hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
534 hSpectra[kTOF][kKaon] [kNeg]->SetName("hkaNeg");
536 // Divide for efficiency
537 hSpectra[kTOF][kPion] [kPos]->Divide(hEffTrackTOF[kPion] [kPos]);
538 hSpectra[kTOF][kKaon] [kPos]->Divide(hEffTrackTOF[kKaon] [kPos]);
539 hSpectra[kTOF][kProton][kPos]->Divide(hEffTrackTOF[kProton][kPos]);
540 hSpectra[kTOF][kPion] [kPos]->Divide(hEffMatchTOF[kPion] [kPos]);
541 hSpectra[kTOF][kKaon] [kPos]->Divide(hEffMatchTOF[kKaon] [kPos]);
542 hSpectra[kTOF][kProton][kPos]->Divide(hEffMatchTOF[kProton][kPos]);
543 hSpectra[kTOF][kPion] [kNeg]->Divide(hEffTrackTOF[kPion] [kNeg]);
544 hSpectra[kTOF][kKaon] [kNeg]->Divide(hEffTrackTOF[kKaon] [kNeg]);
545 hSpectra[kTOF][kProton][kNeg]->Divide(hEffTrackTOF[kProton][kNeg]);
546 hSpectra[kTOF][kPion] [kNeg]->Divide(hEffMatchTOF[kPion] [kNeg]);
547 hSpectra[kTOF][kKaon] [kNeg]->Divide(hEffMatchTOF[kKaon] [kNeg]);
548 hSpectra[kTOF][kProton][kNeg]->Divide(hEffMatchTOF[kProton][kNeg]);
551 // Clean UP TOF spectra, removing unwanted points
552 cout << "Cleaning Up TOF spectra" << endl;
553 Int_t nbin = hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
554 for(Int_t ibin = 1; ibin <= nbin; ibin++){
555 Float_t pt = hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
556 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
558 hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
559 hSpectra[kTOF][kKaon][icharge]->SetBinError (ibin,0);
560 hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
561 hSpectra[kTOF][kProton][icharge]->SetBinError (ibin,0);
565 // cout << "Scaling TOF to TPC" << endl;
566 // // Scale protons to TPC
567 // hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
568 // // Scale pbar so that pbar/p is compatible with Panos
569 // hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
572 // Load TOF systematic errors:
573 f = new TFile ("./Files/systMatchingPos.root");
574 AliBWTools::AddHisto(hSystError[kTOF][kPion][kPos] ,(TH1*)gDirectory->Get("hErrPionMatch") );
575 AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hErrProtonMatch"));
576 AliBWTools::AddHisto(hSystError[kTOF][kKaon][kPos] ,(TH1*)gDirectory->Get("hErrKaonMatch") );
577 f = new TFile ("./Files/systMatchingNeg.root");
578 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kNeg],(TH1*)gDirectory->Get("hErrPionMatch"));
579 AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hErrProtonMatch"));
580 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kNeg],(TH1*)gDirectory->Get("hErrKaonMatch"));
581 f = new TFile ("./Files/systPIDall.root");
582 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kPos],(TH1*)gDirectory->Get("hpiCorr"));
583 AliBWTools::AddHisto(hSystError[kTOF][kProton][kPos],(TH1*)gDirectory->Get("hpCorr"));
584 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kPos],(TH1*)gDirectory->Get("hkCorr"));
585 AliBWTools::AddHisto(hSystError[kTOF][kPion] [kNeg],(TH1*)gDirectory->Get("hpiCorr"));
586 AliBWTools::AddHisto(hSystError[kTOF][kProton][kNeg],(TH1*)gDirectory->Get("hpCorr"));
587 AliBWTools::AddHisto(hSystError[kTOF][kKaon] [kNeg],(TH1*)gDirectory->Get("hkCorr"));
591 f = new TFile("./Files/ITSsaSpectraCorr-28Sep10.root");
592 hSpectra[kITS][kPion] [kPos]= (TH1F*) f->Get("hSpectraPos0");
593 hSpectra[kITS][kKaon] [kPos]= (TH1F*) f->Get("hSpectraPos1");
594 hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
595 hSpectra[kITS][kPion] [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
596 hSpectra[kITS][kKaon] [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
597 hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
599 for(Int_t ipart = 0; ipart < kNPart; ipart++){
600 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
601 Int_t nbinITS = hSpectra[kITS][ipart][icharge]->GetNbinsX();
602 for(Int_t ibin = 1; ibin <= nbinITS; ibin++){
603 if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
604 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
605 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
607 if(ipart == kProton && ibin==9){
608 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));
609 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
610 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
612 // if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
613 // hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
614 // hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
621 if(useSecCorrFromDCA){
622 TFile* fseccorr = new TFile("./Files/CorrFac-SecProtons3Meth-ITSsa.root");
623 TH1F* hcorrp=(TH1F*)fseccorr->Get("hSecPCorrFromDCAexpo");
624 TH1F* hcorrpbar=(TH1F*)fseccorr->Get("hSecPbarCorrFromDCAexpo");
625 hSpectra[kITS][kProton][kPos]->Multiply(hcorrp);
626 hSpectra[kITS][kProton][kNeg]->Multiply(hcorrpbar);
630 // Load ITS sa systematics, only pt dependent ones
631 // f = TFile::Open("./Files/ITSsa-systematics.root");
632 f = TFile::Open("./Files/ITSsa-systematics_20100928.root");
633 for(Int_t ipart = 0; ipart < kNPart; ipart++){
634 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
635 AliBWTools::AddHisto(hSystError[kITS][ipart][icharge], (TH1*) gDirectory->Get(Form("hSystTot%s%s",chargeFlag[icharge],partFlag[ipart]))); // Using total error
640 //f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons20100720.root");
641 f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons14092010new.root");
642 TList * list = (TList*) gDirectory->Get("output");
643 hSpectra[kITSTPC][kPion] [kPos]= (TH1F*) list->FindObject("Pions");
644 hSpectra[kITSTPC][kKaon] [kPos]= (TH1F*) list->FindObject("Kaons");
645 hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
646 hSpectra[kITSTPC][kPion] [kNeg]= (TH1F*) list->FindObject("AntiPions");
647 hSpectra[kITSTPC][kKaon] [kNeg]= (TH1F*) list->FindObject("AntiKaons");
648 hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
651 f = new TFile("./Files/protonSpectra_20100615.root");
652 hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
653 hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
654 f = new TFile("./Files/pionSpectra_20100615.root");
655 hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
656 hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
657 // f = new TFile("./Files/kaonSpectra_20100615.root");
658 f = new TFile("./Files/kaonsTPCcorr_2010_08_31.root");
659 hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
660 hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
662 // Clean UP TPC spectra, removing unwanted points
663 cout << "Cleaning Up TPC spectra" << endl;
664 nbin = hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
665 for(Int_t ibin = 0; ibin < nbin; ibin++){
666 Float_t pt = hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
667 if (pt > 0.45){ // || pt<0.25) {
668 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
669 hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
670 hSpectra[kTPC][kKaon][icharge]->SetBinError (ibin,0);
674 hSpectra[kTPC][kKaon][kNeg]->SetBinContent(ibin,0);
675 hSpectra[kTPC][kKaon][kNeg]->SetBinError (ibin,0);
678 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
679 hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
680 hSpectra[kTPC][kProton][icharge]->SetBinError (ibin,0);
685 // Load TPC systematics
686 cout << "WARNING: TPC SYST FOR NEGATIVES TO BE CORRECTED" << endl;
687 f = TFile ::Open("./Files/pionsSystSum.root");
688 AliBWTools::AddHisto(hSystError[kTPC][kPion][kPos],(TH1*) gDirectory->Get("pionsSyst"));
689 AliBWTools::AddHisto(hSystError[kTPC][kPion][kNeg],(TH1*) gDirectory->Get("pionsSyst"));
690 f = TFile ::Open("./Files/kaonsSystSum.root");
691 AliBWTools::AddHisto(hSystError[kTPC][kKaon][kPos],(TH1*) gDirectory->Get("kaonsSyst"));
692 AliBWTools::AddHisto(hSystError[kTPC][kKaon][kNeg],(TH1*) gDirectory->Get("kaonsSyst"));
693 f = TFile ::Open("./Files/ProtonsSystSum.root");
694 AliBWTools::AddHisto(hSystError[kTPC][kProton][kPos],(TH1*) gDirectory->Get("ProtonsSyst"));
695 AliBWTools::AddHisto(hSystError[kTPC][kProton][kNeg],(TH1*) gDirectory->Get("ProtonsSyst"));
698 f = new TFile ("./Files/PtSpectraCorrectedK0sOff_20100803.root");
699 // hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
700 hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
701 // hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
702 // hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
703 hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
706 // f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
707 // f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
708 // f = new TFile ("./Files/KaonKinkJun16PhySel2N.root");
709 f = new TFile ("./Files/KaonKinkJun16PhySel2NRebPass6.root");
710 hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
711 hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
712 // hSpectra[kKinks][kKaon][kPos]->Scale(0.5/0.7); // different rapidity range for kinks
713 // hSpectra[kKinks][kKaon][kNeg]->Scale(0.5/0.7); // different rapidity range for kinks
714 // hSpectra[kKinks][kKaon][kPos]->Scale(276004./263345.); // different N of events
715 // hSpectra[kKinks][kKaon][kNeg]->Scale(276004./263345.); // different N of events
716 // This normalization if for KaonKinkJun16PhySel2N
717 // hSpectra[kKinks][kKaon][kPos]->Scale(1./303214); // different N of events
718 // hSpectra[kKinks][kKaon][kNeg]->Scale(1./303214); // different N of events
719 // This normalization is for KaonKinkJun16PhySel2NRebPass6
720 hSpectra[kKinks][kKaon][kPos]->Scale(1./260432.26); // different N of events
721 hSpectra[kKinks][kKaon][kNeg]->Scale(1./260432.26); // different N of events
723 // Apply correction factors
724 // Secondaries for protons
726 //f = new TFile ("./Files/corrFactorProtons_20100615.root");
727 f = new TFile ("./Files/corrFactorProtons_2010_09_24.root");
728 TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
729 if(correctSecondaries) {
730 cout << "CORRECTING SECONDARIES" << endl;
732 for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
733 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
734 Int_t ipart = kProton;
735 TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
737 Int_t nbins = h->GetNbinsX();
738 for(Int_t ibin = 0; ibin < nbins; ibin++){
739 // Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
740 Float_t pt = h->GetBinCenter(ibin);
741 if (icharge == kNeg) pt = -pt;
742 Int_t binCorrection = hCorrSecondaries->FindBin(pt);
743 Float_t correction = hCorrSecondaries->GetBinContent(binCorrection);
744 // cout << pt << " " << correction << endl;
746 if (correction) {// If the bin is empty this is a 0
747 h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
748 h->SetBinError (ibin,h->GetBinError (ibin)/correction);
749 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
750 cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
751 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
759 // geant/fluka absorption
760 if(correctGeantFlukaXS) {
761 cout << "CORRECTING GEANT3/FLUKA" << endl;
762 // common file for Kaons
763 TFile *fFlukakaons = TFile::Open("./Files/correctionForCrossSection.321.root");
764 TH1D * hCorrFlukakaon[kNCharge];
765 hCorrFlukakaon[kPos] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionParticles");
766 hCorrFlukakaon[kNeg] = (TH1D*)fFlukakaons->Get("gHistCorrectionForCrossSectionAntiParticles");
768 for(Int_t idet = 0; idet < kNDet; idet++){
769 if( idet != kITS) continue; // comment to use fluka for kaons for all dets
770 if (idet == kTOF) continue; // TOF already corrected
771 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
773 TH1 * h = hSpectra[idet][ipart][icharge]; // only ITS sa
775 Int_t nbins = h->GetNbinsX();
776 Int_t nbinsy=hCorrFlukakaon[icharge]->GetNbinsY();
777 for(Int_t ibin = 0; ibin < nbins; ibin++){
778 Float_t pt = h->GetBinCenter(ibin);
779 Float_t minPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(1);
780 Float_t maxPtCorrection = hCorrFlukakaon[icharge]->GetXaxis()->GetBinLowEdge(nbinsy+1);
781 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
782 if (pt > maxPtCorrection) pt = maxPtCorrection;
783 Float_t correction = hCorrFlukakaon[icharge]->GetBinContent(hCorrFlukakaon[icharge]->GetXaxis()->FindBin(pt));
784 if (correction != 0) {// If the bin is empty this is a 0
785 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
786 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
787 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
788 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
789 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
798 // ITS specific file for protons/antiprotons
799 TFile* fITS = new TFile ("./Files/correctionForCrossSectionITS_20100719.root");
800 TH2D * hCorrFlukaITS[kNCharge];
801 hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
802 hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
804 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
805 Int_t ipart = kProton;
806 TH1 * h = hSpectra[kITS][ipart][icharge]; // only ITS sa
808 Int_t nbins = h->GetNbinsX();
809 Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
810 for(Int_t ibin = 0; ibin < nbins; ibin++){
811 Float_t pt = h->GetBinCenter(ibin);
812 Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
813 Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
814 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
815 if (pt > maxPtCorrection) pt = maxPtCorrection;
816 Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
817 if (correction != 0) {// If the bin is empty this is a 0
818 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
819 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
820 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
821 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << chargeFlag[icharge] << endl;
822 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
830 //f = new TFile ("./Files/correctionForCrossSection_20100615.root");
831 f = new TFile ("./Files/correctionForCrossSection_20100920.root");
832 TH2D * hCorrFluka[kNCharge];
833 hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
834 hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
835 for(Int_t idet = 0; idet < kNDet; idet++){
836 if (idet == kITS) continue; // skip ITS sa
837 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
838 Int_t ipart = kProton;
839 TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
841 Int_t nbins = h->GetNbinsX();
842 for(Int_t ibin = 0; ibin < nbins; ibin++){
843 // Float_t pt = convertToMT ?
844 // TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
845 // h->GetBinCenter(ibin);
846 Float_t pt = h->GetBinCenter(ibin);
847 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
848 // hCorrFluka[icharge]->Draw();
849 TH1D * htmpFluka=hCorrFluka[icharge]->ProjectionY();
850 Float_t maxPtCorrection = AliBWTools::GetHighestNotEmptyBinEdge(htmpFluka);//->GetYaxis()->GetBinLowEdge(hCorrFluka[icharge]->GetNbinsY()+1);
851 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
852 if (pt > maxPtCorrection) pt = maxPtCorrection-0.0001;
853 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
855 // already in the efficiency correction (F. Noferini)
857 if (icharge == kNeg) correction = 1; // antiprotons already corrected in efficiency
858 // Scale correction for the different material budget. Recipe by Francesco Noferini
860 correction = TMath::Power(correction,0.07162/0.03471);
863 if (correction != 0) {// If the bin is empty this is a 0
865 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
866 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
867 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
868 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
869 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
880 // Set style and scale
881 for(Int_t idet = 0; idet < kNDet; idet++){
882 for(Int_t ipart = 0; ipart < kNPart; ipart++){
883 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
884 if (hSpectra[idet][ipart][icharge]){
885 hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
886 hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
887 hSpectra[idet][ipart][icharge]->SetLineColor (color[ipart] );
888 hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
889 hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
890 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
892 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
893 hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
894 hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
895 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
897 // if (idet == kTOF || idet == kTPC) {
898 // hSpectra[idet][ipart][icharge]->Scale(1./236763);
900 // if (idet == kITS){
901 // hSpectra[idet][ipart][icharge]->Scale(202945./236763);
903 if (scaleKaons && ipart == kKaon) {
904 hSpectra[idet][ipart][icharge]->Scale(2.);
907 printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
914 // Create fake weights for the mean; To be update once we have syst errors
915 // Using syste from table in paper. It would be better to have this as a function of pt.
916 TH1F * hWeights[3][kNPart];
917 const Double_t kWeights[3][kNPart] =
918 // {{4, 3, 10.2}, // TPC
919 // {4.1,8.8,7.0 }, //TOF
920 // {4.5,5.6,7.0 }}; // ITS
921 {{0.1,0.1,0.1}, // TPC
923 {0.1,0.1,0.1}}; // ITS
924 for(Int_t ipart = 0; ipart <= kNPart ; ipart++){
925 for(Int_t idet = 0; idet <= kITS ; idet++){
926 hWeights[idet][ipart] = (TH1F*) hSpectra[idet][ipart][kPos]->Clone();
927 Int_t nbinW = hWeights[idet][ipart]->GetNbinsX();
928 for(Int_t ibin = 1; ibin <= nbinW; ibin++){
929 hWeights[idet][ipart]->SetBinError(ibin, kWeights[idet][ipart]);
933 const Double_t scaleDet[] = {1.00,1.00,1.00}; // During the weekly meeting 19/08/2010 we decided not to apply any scaling.
934 // 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%
935 // 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%
939 for(Int_t ipart = 0; ipart < kNPart; ipart++){
940 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
941 TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
944 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
945 htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
948 hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
949 hSpectra[kTPC][ipart][icharge],
952 hSpectra[kCombAll][ipart][icharge] =
953 AliBWTools::Combine3HistosWithErrors(hSpectra[kITS][ipart][icharge], // Histos to combine
954 hSpectra[kTPC][ipart][icharge],
955 hSpectra[kTOF][ipart][icharge],
956 hSystError[kITS][ipart][icharge], // Errors (weights) used for the average
957 hSystError[kTPC][ipart][icharge],
958 hSystError[kTOF][ipart][icharge],
959 // hWeights[kITS][ipart],
960 // hWeights[kTPC][ipart],
961 // hWeights[kTOF][ipart],
962 htemplLocal,1, // template, take statistical error from TPC in overlap
963 scaleDet[kITS], // relative scaling
966 (TH1**)&hSystError[kCombAll][ipart][icharge], // put combined syst error here
967 1 // weights histos contain error in bin content
969 // if (convertToMT) {
970 // TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
971 // hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
972 // hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
973 // hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
979 // Scale for the number of inelastic collisions and correct for
980 // efficiency losses due to physics selection:
982 Double_t effPhysSel[kNPart];
983 effPhysSel[kPion] = 1.012;
984 effPhysSel[kKaon] = 1.013;
985 effPhysSel[kProton] = 1.014;
988 for(Int_t idet = 0; idet < kNHist; idet++){
989 for(Int_t ipart = 0; ipart < kNPart; ipart++){
990 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
991 if(hSpectra[idet][ipart][icharge]) {
992 // cout << "Scaling!" << endl;
993 if(idet != kKinks && idet != kK0){ // Kinks use a different run list, k0s normalized by Helene
994 hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel tutti? // FIXME
1007 const char * evClass= "INEL";
1008 const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
1009 "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
1010 "./Files/dndeta_CMSD6T109_10M_900GeV.root",
1011 "./Files/dndeta_Perugia0320_10M_900GeV.root", };
1014 for(Int_t itune = 0; itune < kNTunes; itune++){
1015 f = new TFile(files[itune]);
1016 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1017 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1018 hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
1025 for(Int_t itune = 0; itune < kNTunes; itune++){
1026 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1027 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1028 if (hSpectraMC[itune][ipart][icharge]){
1029 hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
1030 hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
1031 hSpectraMC[itune][ipart][icharge]->SetLineColor (mcLineColor[itune] );
1032 hSpectraMC[itune][ipart][icharge]->SetLineStyle (mcLineStyle[itune] );
1033 hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
1035 TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
1036 hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
1037 hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
1038 hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
1042 printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
1050 void DrawStar(Int_t icharge) {
1051 // cout << icharge << endl;
1053 // gROOT->LoadMacro("StarPPSpectra.C");
1054 TGraphErrors ** gStar = StarPPSpectra();
1056 for(Int_t istar = 0; istar < 6; istar++){
1057 gStar[istar]->SetMarkerStyle(kOpenStar);
1058 if (icharge==kPos && (istar%2) ) gStar[istar]->Draw("P");
1059 else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
1060 else cout << "Skipping Star " << istar << endl;
1064 void GetITSResiduals() {
1067 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1068 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1069 // cout << "1 " << ipart << " " << icharge << endl;
1071 // gSystem->ProcessEvents();
1072 // gSystem->Sleep(1000);
1073 TF1* f = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
1074 TH1F * hres1, *hres2;
1075 AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
1077 TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
1079 hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
1080 hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
1081 hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
1082 hSpectra[kITS][ipart][icharge]->Draw("same");
1083 hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
1084 TLegend* l = new TLegend( 0.182886, 0.192308, 0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
1085 l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
1086 l->AddEntry(hSpectra[kITS][ipart][icharge], "ITS");
1087 l->AddEntry(f, "Fit to TOF+TPC");
1090 TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
1091 TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
1093 hres2->SetMinimum(-1);
1094 hres2->SetMaximum(1);
1096 hres2->GetYaxis()->SetTitleOffset(1.2);
1097 Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
1098 TLine * line = new TLine(x,-1,x,1);
1099 line->SetLineStyle(kDashed);
1105 gSystem->ProcessEvents();
1106 c1->Print(TString(c1->GetName()) + ".png");
1107 c2->Print(TString(c2->GetName()) + ".png");
1113 void DrawWithModels() {
1115 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1118 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1120 TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],700,700);
1121 TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlag[icharge], "p1", 0.0, 0.35, 1.0, 0.95, 0, 0, 0);
1122 p1->SetBottomMargin(0);
1125 TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlag[icharge], "p2", 0.0, 0.05, 1.0, 0.35, 0, 0, 0);
1126 p2->SetTopMargin(0);
1127 p2->SetBottomMargin(0.3);
1130 Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
1135 TH2F * hempty = new TH2F(TString("hempty")+long(icharge)+long(ipart),"hempty",100,0.,4, 100, 0.0015,5);
1136 hempty->SetXTitle("p_{t} (GeV/c)");
1137 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1143 TLegend * l =new TLegend(0.502874, 0.493056, 0.892241, 0.904762);
1144 l->SetFillColor(kWhite);
1145 hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
1146 l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data"));
1147 for(Int_t itune = 0; itune < kNTunes; itune++){
1148 l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1149 hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
1150 AliBWTools::GetGraphFromHisto(hSpectraMC[itune][ipart][icharge])->Draw("CX");
1154 TLatex * tex = new TLatex(0.6712643,2.353486,partLabel[ipart][icharge]);
1155 tex->SetTextFont(42);
1156 tex->SetTextSize(0.07936508);
1157 tex->SetLineWidth(2);
1162 TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(icharge)+long(ipart),"hempty",100,0.,4, 100, 0.01,2.99);
1163 hemptyr->SetXTitle("p_{t} (GeV/c)");
1164 hemptyr->SetYTitle("Data/MC");
1165 hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
1166 hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
1167 hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
1168 hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
1169 hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
1170 hemptyr->GetXaxis()->SetTitleOffset(1.05);
1171 hemptyr->SetTickLength(0.03*scaleFonts, "X");
1172 hemptyr->SetTickLength(0.02*scaleFonts, "Y");
1173 // hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
1174 hemptyr->GetYaxis()->SetNdivisions(5);
1178 for(Int_t itune = 0; itune < kNTunes; itune++){
1179 TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
1181 // l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
1182 TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[iCombInStudy][ipart][icharge],f);
1183 hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
1184 hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
1185 hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
1186 AliBWTools::GetGraphFromHisto(hRatio)->Draw("CX");
1193 gSystem->ProcessEvents();
1194 c1->Print(TString(c1->GetName())+".eps");
1195 // ALICEWorkInProgress(c1,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
1196 c1->Print(TString(c1->GetName())+".png");
1206 void DrawAllAndKaons() {
1210 // gROOT->LoadMacro("ALICEWorkInProgress.C");
1212 // gStyle->SetOptFit(0);
1213 gStyle->SetStatX(0.9);
1214 gStyle->SetStatY(0.88);
1215 gStyle->SetStatW(0.4);
1216 gStyle->SetStatH(0.1);
1218 TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1219 hKaonsAllTPCTOF->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1221 TH1F * hK0Scaled = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
1222 hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
1224 hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
1225 hSpectra[kKinks][kKaon][kPos]->SetLineColor(4);
1226 hSpectra[kKinks][kKaon][kPos]->SetStats(0);
1227 TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
1228 hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
1230 TCanvas * c1 = new TCanvas("cKaons","cKaons",700,700);
1232 TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-3,6);
1233 hempty->SetXTitle("p_{t} (GeV/c)");
1234 // hempty->SetYTitle("dN / dp_{t} (A.U.)");
1235 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1237 hK0Scaled->Draw("same");
1238 hKaonsAllTPCTOF->Draw("same");
1239 hKinksAll->Draw("same");
1242 TLegend * leg = new TLegend(0.456897, 0.71875, 0.892241, 0.936012,NULL,"brNDC");
1243 // leg->SetBorderSize(0);
1244 // leg->SetLineColor(1);
1245 // leg->SetLineStyle(1);
1246 // leg->SetLineWidth(1);
1247 // leg->SetFillColor(19);
1248 leg->SetFillColor(0);
1249 TLegendEntry *entry=leg->AddEntry(hKaonsAllTPCTOF,"K^{+} + K^{-}, ITS+TPC+TOF ","lpf");
1250 entry=leg->AddEntry(hK0Scaled,"K^{0} #times 2","lpf");
1251 entry=leg->AddEntry(hKinksAll,"K^{+} + K ^{-}, Kinks","lpf");
1254 // ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Prelimiary}{Statistical Error Only}");
1255 // TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
1256 // tex->SetTextColor(2);
1257 // tex->SetTextFont(42);
1258 // tex->SetTextSize(0.03496503);
1262 PrintCanvas(c1, printFormats);
1264 // Draw all "stable" hadrons
1265 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1266 TCanvas * c1h = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge],700,700);
1268 c1h->SetLeftMargin(0.14);
1269 TH2F * hemptyLoc = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-3,10);
1270 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1271 hemptyLoc->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1272 hemptyLoc->GetYaxis()->SetTitleOffset(1.35);
1273 hemptyLoc->GetXaxis()->SetTitleOffset(1.1);
1275 leg = new TLegend(0.482759, 0.489583, 0.896552, 0.925595, NULL,"brNDC");
1276 leg->SetFillColor(0);
1277 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1278 for(Int_t idet = 0; idet <= kITSTPC; idet++){
1279 // if (idet == kITS) continue;
1280 // if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
1281 hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
1282 hSpectra[idet][ipart][icharge]->Draw("same");
1283 leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet] + ")","lpf");
1288 // ALICEWorkInProgress(c1h,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1289 PrintCanvas(c1h,printFormats);
1295 // K-/K+ in the different detectors
1296 TCanvas * cpm=new TCanvas("cpm","Kminus/Kplus",700,700);
1299 TH1F* hRatioKPKM_TPC=new TH1F(*(hSpectra[kTPC][kKaon][kNeg]));
1300 hRatioKPKM_TPC->SetMinimum(0.5);
1301 hRatioKPKM_TPC->SetMaximum(1.5);
1302 hRatioKPKM_TPC->Divide(hSpectra[kTPC][kKaon][kPos]);
1303 hRatioKPKM_TPC->GetYaxis()->SetTitle("K-/K+ (TPC)");
1304 hRatioKPKM_TPC->Draw();
1306 TH1F* hRatioKPKM_ITS=new TH1F(*(hSpectra[kITS][kKaon][kNeg]));
1307 hRatioKPKM_ITS->Divide(hSpectra[kITS][kKaon][kPos]);
1308 hRatioKPKM_ITS->SetMinimum(0.5);
1309 hRatioKPKM_ITS->SetMaximum(1.5);
1310 hRatioKPKM_ITS->GetYaxis()->SetTitle("K-/K+ (ITSsa)");
1311 hRatioKPKM_ITS->Draw("");
1313 TH1F* hRatioKPKM_TOF=new TH1F(*(hSpectra[kTOF][kKaon][kNeg]));
1314 hRatioKPKM_TOF->Divide(hSpectra[kTOF][kKaon][kPos]);
1315 hRatioKPKM_TOF->SetMinimum(0.5);
1316 hRatioKPKM_TOF->SetMaximum(1.5);
1317 hRatioKPKM_TOF->GetYaxis()->SetTitle("K-/K+ (TOF)");
1318 hRatioKPKM_TOF->Draw("");
1320 TH1F* hRatioKPKM_ITSTPC=new TH1F(*(hSpectra[kITSTPC][kKaon][kNeg]));
1321 hRatioKPKM_ITSTPC->Divide(hSpectra[kITSTPC][kKaon][kPos]);
1322 hRatioKPKM_ITSTPC->SetMinimum(0.5);
1323 hRatioKPKM_ITSTPC->SetMaximum(1.5);
1324 hRatioKPKM_ITSTPC->GetYaxis()->SetTitle("K-/K+ (ITS Global)");
1325 hRatioKPKM_ITSTPC->Draw("");
1329 TH1F * hRatioITSTPC[kNPart][kNCharge];
1330 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1332 TCanvas * c1h = new TCanvas(TString("cITSTPCRatio_")+chargeFlag[icharge],TString("cITSTPCRatio_")+chargeFlag[icharge],1200,500);
1335 TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITSsa/TPC ",100,0.,1., 100, 0.5,1.5);
1336 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1337 hemptyLoc->SetYTitle("ITSsa / TPC");
1338 // Loop over particles
1339 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1341 hRatioITSTPC[ipart][icharge]=new TH1F(*hSpectra[kITS][ipart][icharge]);
1342 Int_t nBinsITS=hSpectra[kITS][ipart][icharge]->GetNbinsX();
1343 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1344 // Loop over ITS bins,
1345 for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1346 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1347 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,0.);
1348 Float_t lowPtITS=hSpectra[kITS][ipart][icharge]->GetBinLowEdge(iBin);
1349 Float_t binWidITS=hSpectra[kITS][ipart][icharge]->GetBinWidth(iBin);
1350 // Loop over TPC bins and find overlapping bins to ITS
1351 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1352 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1353 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1354 if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1355 Float_t numer=hSpectra[kITS][ipart][icharge]->GetBinContent(iBin);
1356 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1357 Float_t enumer=hSpectra[kITS][ipart][icharge]->GetBinError(iBin);
1358 Float_t edenom=hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1361 if(numer>0. && denom>0.){
1363 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1365 hRatioITSTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1366 hRatioITSTPC[ipart][icharge]->SetBinError(iBin,eratio);
1372 // hemptyLoc->SetStats(1);
1373 // hemptyLoc->Draw();
1374 hRatioITSTPC[ipart][icharge]->SetYTitle("ITSsa/TPC");
1375 hRatioITSTPC[ipart][icharge]->SetStats(1);
1376 hRatioITSTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1377 hRatioITSTPC[ipart][icharge]->Draw("");
1378 hRatioITSTPC[ipart][icharge]->Fit("pol0","","same");
1381 PrintCanvas(c1h,printFormats);
1385 TH1F * hRatioITSGlobalTPC[kNPart][kNCharge];
1386 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1388 TCanvas * c1h = new TCanvas(TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],TString("cITSGlobalTPCRatio_")+chargeFlag[icharge],1200,500);
1391 TH2F * hemptyLoc = new TH2F(TString("hemptyR")+long(icharge),"ITS Global/TPC ",100,0.,1., 100, 0.5,1.5);
1392 hemptyLoc->SetXTitle("p_{t} (GeV/c)");
1393 hemptyLoc->SetYTitle("ITS Global / TPC");
1394 // Loop over particles
1395 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
1397 hRatioITSGlobalTPC[ipart][icharge]=new TH1F(*hSpectra[kITSTPC][ipart][icharge]);
1398 Int_t nBinsITS=hSpectra[kITSTPC][ipart][icharge]->GetNbinsX();
1399 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1400 // Loop over ITS bins,
1401 for(Int_t iBin=1; iBin<=nBinsITS; iBin++){
1402 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1403 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,0.);
1404 Float_t lowPtITS=hSpectra[kITSTPC][ipart][icharge]->GetBinLowEdge(iBin);
1405 Float_t binWidITS=hSpectra[kITSTPC][ipart][icharge]->GetBinWidth(iBin);
1406 // Loop over TPC bins and find overlapping bins to ITS
1407 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1408 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1409 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1410 if(TMath::Abs(lowPtITS-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidITS)<0.001){
1411 Float_t numer=hSpectra[kITSTPC][ipart][icharge]->GetBinContent(iBin);
1412 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1413 Float_t enumer=hSpectra[kITSTPC][ipart][icharge]->GetBinError(iBin);
1414 Float_t edenom=hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1417 if(numer>0. && denom>0.){
1419 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1421 hRatioITSGlobalTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1422 hRatioITSGlobalTPC[ipart][icharge]->SetBinError(iBin,eratio);
1428 // hemptyLoc->SetStats(1);
1429 // hemptyLoc->Draw();
1430 hRatioITSGlobalTPC[ipart][icharge]->SetYTitle("ITS Global/TPC");
1431 hRatioITSGlobalTPC[ipart][icharge]->SetStats(1);
1432 hRatioITSGlobalTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1433 hRatioITSGlobalTPC[ipart][icharge]->Draw("");
1434 hRatioITSGlobalTPC[ipart][icharge]->Fit("pol0","","same");
1437 PrintCanvas(c1h,printFormats);
1441 TH1F * hRatioTOFTPC[kNPart][kNCharge];
1442 for(Int_t icharge = 0; icharge < kNCharge; icharge++){ // loop over charges
1444 TCanvas * c1t = new TCanvas(TString("cTOFTPCRatio_")+chargeFlag[icharge],TString("cTOFTPCRatio_")+chargeFlag[icharge],1200,500);
1447 TH2F * hemptyt = new TH2F(TString("hemptyRt")+long(icharge),"TOF/TPC ",100,0.,1., 100, 0.5,1.5);
1448 hemptyt->SetXTitle("p_{t} (GeV/c)");
1449 hemptyt->SetYTitle("TOF / TPC");
1451 for(Int_t ipart = 0; ipart < kNPart; ipart++) { // loop over particles
1453 hRatioTOFTPC[ipart][icharge]=new TH1F(*hSpectra[kTOF][ipart][icharge]);
1454 Int_t nBinsTOF=hSpectra[kTOF][ipart][icharge]->GetNbinsX();
1455 Int_t nBinsTPC=hSpectra[kTPC][ipart][icharge]->GetNbinsX();
1456 // Loop over TOF bins
1457 for(Int_t iBin=1; iBin<=nBinsTOF; iBin++){
1458 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1459 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,0.);
1460 Float_t lowPtTOF=hSpectra[kTOF][ipart][icharge]->GetBinLowEdge(iBin);
1461 Float_t binWidTOF=hSpectra[kTOF][ipart][icharge]->GetBinWidth(iBin);
1462 // Loop over TPC bins and find overlapping bins to ITS
1463 for(Int_t jBin=1; jBin<=nBinsTPC; jBin++){
1464 Float_t lowPtTPC=hSpectra[kTPC][ipart][icharge]->GetBinLowEdge(jBin);
1465 Float_t binWidTPC=hSpectra[kTPC][ipart][icharge]->GetBinWidth(jBin);
1466 if(TMath::Abs(lowPtTOF-lowPtTPC)<0.001 && TMath::Abs(binWidTPC-binWidTOF)<0.001){
1467 Float_t numer=hSpectra[kTOF][ipart][icharge]->GetBinContent(iBin);
1468 Float_t denom=hSpectra[kTPC][ipart][icharge]->GetBinContent(jBin);
1469 Float_t enumer=hSpectra[kTOF][ipart][icharge]->GetBinError(iBin);
1470 Float_t edenom=hSpectra[kTPC][ipart][icharge]->GetBinError(jBin);
1473 if(numer>0. && denom>0.){
1475 eratio=ratio*TMath::Sqrt((enumer/numer)*(enumer/numer)+(edenom/denom)*(edenom/denom));
1477 hRatioTOFTPC[ipart][icharge]->SetBinContent(iBin,ratio);
1478 hRatioTOFTPC[ipart][icharge]->SetBinError(iBin,eratio);
1484 hRatioTOFTPC[ipart][icharge]->SetYTitle("TOF/TPC");
1485 hRatioTOFTPC[ipart][icharge]->SetStats(1);
1486 hRatioTOFTPC[ipart][icharge]->GetYaxis()->SetRangeUser(0.5,1.5);
1487 hRatioTOFTPC[ipart][icharge]->Draw("");
1488 hRatioTOFTPC[ipart][icharge]->Fit("pol0","","same");
1490 PrintCanvas(c1t,printFormats);
1495 void DrawWithJacek() {
1497 //1. Convert spectra to dNdeta and sum
1498 TH1F * hsum = (TH1F*) htemplate->Clone();
1500 Int_t idet= iCombInStudy;
1501 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1502 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1503 TH1 * h = hSpectra[idet][ipart][icharge];
1504 Int_t nbin = h->GetNbinsX();
1505 for(Int_t ibin = 1; ibin <= nbin; ibin++){
1506 Double_t pt = h->GetBinCenter(ibin);
1507 Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
1508 Double_t jacobian = pt/mt;
1509 h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
1510 h->SetBinError (ibin,h->GetBinError (ibin)*jacobian);
1511 Int_t ibinSum = hsum->FindBin(pt);
1512 Double_t epsilon = 0.0001;
1513 if ( h->GetBinContent(ibin) > 0 &&
1514 (TMath::Abs(h->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1515 TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1517 cout << "DISCREPANCY IN BIN RANGES" << endl;
1518 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
1519 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
1520 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
1525 hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
1526 hsum->SetBinError (ibinSum,0);
1533 // Load Jacek and Draw both:
1534 new TFile ("./Files/pt_allcharged_published.root");
1535 TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("dndpt");
1536 gJacek->SetMarkerStyle(20);
1537 gJacek->SetMarkerColor(kGray);
1539 gJacek->GetHistogram()->SetXTitle("p_{t} (GeV)");
1540 gJacek->GetHistogram()->SetYTitle("d^{2}N/dp_{t}d#eta (GeV^{-1})");
1543 TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
1552 void DrawRatioToStar() {
1555 // gROOT->LoadMacro("StarPPSpectra.C");
1556 TGraphErrors ** gStar = StarPPSpectra();
1557 gStar[0]->SetMarkerStyle(kOpenStar);
1558 gStar[1]->SetMarkerStyle(kFullStar);
1559 gStar[2]->SetMarkerStyle(kOpenStar);
1560 gStar[3]->SetMarkerStyle(kFullStar);
1561 gStar[4]->SetMarkerStyle(kOpenStar);
1562 gStar[5]->SetMarkerStyle(kFullStar);
1564 // ALICE, INEL -> NSD
1565 Double_t scaleYield = 3.58/3.02; // from paper 2
1566 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1567 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1568 hSpectra[iCombInStudy][ipart][icharge]->Scale(scaleYield);
1573 TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
1574 TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
1575 hempty->SetXTitle("p_{t} (GeV/c)");
1576 hempty->SetYTitle("ALICE/STAR (NSD)");
1579 TCanvas * c1Comp = new TCanvas("cCompToStarNeg","cCompToStarNeg");
1581 TH2F * hempty2 = new TH2F(TString("hemptyCompNeg"),"hemptyCompNeg",100,0.,1.5, 100, 0.001,10);
1582 hempty2->SetXTitle("p_{t} (GeV/c)");
1583 hempty2->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1} [NSD]");
1586 TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
1587 leg->SetBorderSize(0);
1588 leg->SetLineColor(1);
1589 leg->SetLineStyle(1);
1590 leg->SetLineWidth(1);
1591 leg->SetFillColor(0);
1592 leg->SetFillStyle(1001);
1598 g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[iCombInStudy][kPion][kNeg],1);
1599 g->SetMarkerStyle(kFullCircle);
1600 g->SetMarkerColor(kBlack);
1602 leg->AddEntry(g,"#pi^{-}","lp");
1603 g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[iCombInStudy][kKaon][kNeg],1);
1604 g->SetMarkerStyle(kOpenCircle);
1605 g->SetMarkerColor(kRed);
1607 leg->AddEntry(g,"K^{-}","lp");
1608 g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[iCombInStudy][kProton][kNeg],1);
1609 g->SetMarkerStyle(kOpenSquare);
1610 g->SetMarkerColor(kBlue);
1612 leg->AddEntry(g,"#bar{p}","lp");
1616 gStar[0]->Draw("p");
1617 hSpectra[iCombInStudy][kPion][kNeg]->Draw("same");
1618 gStar[2]->Draw("p");
1619 hSpectra[iCombInStudy][kKaon][kNeg]->Draw("same");
1620 gStar[4]->Draw("p");
1621 hSpectra[iCombInStudy][kProton][kNeg]->Draw("same");
1625 TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
1627 leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
1628 leg->SetBorderSize(0);
1629 leg->SetLineColor(1);
1630 leg->SetLineStyle(1);
1631 leg->SetLineWidth(1);
1632 leg->SetFillColor(0);
1633 leg->SetFillStyle(1001);
1635 TCanvas * c2Comp = new TCanvas("cCompToStarPos","cCompToStarPos");
1640 // TGraphErrors * g ;
1641 g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[iCombInStudy][kPion][kPos],1);
1642 g->SetMarkerStyle(kFullCircle);
1643 g->SetMarkerColor(kBlack);
1645 leg->AddEntry(g,"#pi^{+}","lp");
1646 g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[iCombInStudy][kKaon][kPos],1);
1647 g->SetMarkerStyle(kOpenCircle);
1648 g->SetMarkerColor(kRed);
1650 leg->AddEntry(g,"K^{+}","lp");
1651 g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[iCombInStudy][kProton][kPos],1);
1652 g->SetMarkerStyle(kOpenSquare);
1653 g->SetMarkerColor(kBlue);
1655 leg->AddEntry(g,"p","lp");
1660 gStar[1]->Draw("p");
1661 hSpectra[iCombInStudy][kPion][kPos]->Draw("same");
1662 gStar[3]->Draw("p");
1663 hSpectra[iCombInStudy][kKaon][kPos]->Draw("same");
1664 gStar[5]->Draw("p");
1665 hSpectra[iCombInStudy][kProton][kPos]->Draw("same");
1670 gSystem->ProcessEvents();
1671 c1->Print(TString(c1->GetName()) + ".eps");
1672 c2->Print(TString(c2->GetName()) + ".eps");
1673 ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1674 ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1677 c1->Print(TString(c1->GetName()) + ".png");
1678 c2->Print(TString(c2->GetName()) + ".png");
1689 // Draws ratios of combined spectra
1692 TH1F * hPosNegRatio[kNPart];
1694 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1695 hPosNegRatio[ipart] = (TH1F*) hSpectra[iCombInStudy][ipart][kPos]->Clone();
1696 hPosNegRatio[ipart]->Divide(hSpectra[iCombInStudy][ipart][kNeg]);
1697 hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
1698 hPosNegRatio[ipart]->SetMinimum(0.5);
1699 hPosNegRatio[ipart]->SetMaximum(1.5);
1702 TH1F * hKPiRatio = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1703 hKPiRatio->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1704 TH1F * htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1705 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1706 hKPiRatio->Divide(htmp);
1707 hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1709 TH1F * hKPiRatioMC[kNTunes];
1711 for(Int_t itune = 0; itune < kNTunes; itune++){
1712 hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
1713 hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
1714 TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1715 htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1716 hKPiRatioMC[itune]->Divide(htmpLoc);
1717 hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1723 TH1F * hPPiRatio = (TH1F*) hSpectra[iCombInStudy][kProton][kPos]->Clone();
1724 hPPiRatio->Add(hSpectra[iCombInStudy][kProton][kNeg]);
1725 htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1726 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
1727 hPPiRatio->Divide(htmp);
1728 hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1730 TH1F * hPPiRatioMC[kNTunes];
1732 for(Int_t itune = 0; itune < kNTunes; itune++){
1733 hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
1734 hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
1735 TH1F * htmpLoc = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1736 htmpLoc->Add(hSpectraMC[itune][kPion][kNeg]);
1737 hPPiRatioMC[itune]->Divide(htmpLoc);
1738 hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1743 TFile * fPhenix = new TFile("./Files/phenixdata.root");
1744 TH1F * hPPiRatioPhenix = (TH1F*) fPhenix->Get("ProtonPHNX");
1745 hPPiRatioPhenix->Add((TH1F*) fPhenix->Get("PbarPHNX"));
1746 htmp = (TH1F*) fPhenix->Get("PiPlusPHNX");
1747 htmp->Add((TH1F*) fPhenix->Get("PiMinusPHNX"));
1748 hPPiRatioPhenix = AliBWTools::DivideHistosDifferentBins(hPPiRatioPhenix,htmp);
1749 hPPiRatioPhenix->SetMarkerStyle(24);
1752 // TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
1753 // hempty->SetXTitle("p_{t} (GeV/c)");
1754 // hempty->SetYTitle("");
1756 // tmp: overlay levi fits
1757 AliBWFunc * fm2 = new AliBWFunc;
1758 fm2->SetVarType(AliBWFunc::kdNdpt);
1759 TF1 * fLevi[kNPart] [kNCharge];
1760 fLevi[kPion] [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
1761 fLevi[kKaon] [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
1762 fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
1763 fLevi[kPion] [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
1764 fLevi[kKaon] [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
1765 fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
1766 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1767 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1768 fLevi[ipart][icharge]->SetRange(0,4);
1773 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1774 TString detName = detFlag[iCombInStudy];
1775 detName.ReplaceAll(" ", "_");
1776 detName.ReplaceAll("+", "");
1778 TCanvas * c1 = new TCanvas(TString("cRatio_")+detName+TString("_")+partFlag[ipart], TString("cRatio_")+detName+partFlag[ipart]);
1780 hPosNegRatio[ipart]->Draw();
1781 TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
1782 // fRatio->Draw("same");
1783 fRatio->SetRange(0,5);
1786 gSystem->ProcessEvents();
1787 c1->Print(TString(c1->GetName()) + ".png");
1793 TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
1795 hKPiRatio->SetMaximum(0.62);
1797 TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1798 lMC->SetFillColor(kWhite);
1801 gROOT->LoadMacro("GetE735Ratios.C");
1802 GetE735Ratios(0,0)->Draw("EX0,same");
1803 GetE735Ratios(0,1)->Draw("EX0,same");
1804 GetE735Ratios(0,2)->Draw("EX0,same");
1805 GetE735Ratios(0,3)->Draw("EX0,same");
1807 hKPiRatio->SetMarkerStyle(20);
1808 hKPiRatio->SetMarkerColor(kRed);
1809 hKPiRatio->SetLineColor(kRed);
1810 hKPiRatio->Draw("same");
1813 for(Int_t itune = 0; itune < kNTunes; itune++){
1814 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1815 hKPiRatioMC[itune]->SetLineWidth(2);
1816 hKPiRatioMC[itune]->Draw("same,chist");
1823 TLegend * l = new TLegend( 0.1879, 0.68, 0.54,0.92);
1824 l->SetFillColor(kWhite);
1825 l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1826 l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
1827 l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
1828 l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
1829 l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
1834 TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
1836 hPPiRatio->SetMarkerStyle(20);
1837 hPPiRatio->SetMarkerColor(kRed);
1838 hPPiRatio->SetLineColor(kRed);
1840 hPPiRatio->SetMaximum(0.39);
1842 lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1843 lMC->SetFillColor(kWhite);
1845 for(Int_t itune = 0; itune < kNTunes; itune++){
1846 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1847 hPPiRatioMC[itune]->SetLineWidth(2);
1848 hPPiRatioMC[itune]->Draw("same,chist");
1853 hPPiRatioPhenix->Draw("same");
1856 TLegend * l = new TLegend(0.186242, 0.781469, 0.538591, 0.921329);
1857 l->SetFillColor(kWhite);
1858 l->AddEntry(hPPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1859 l->AddEntry(hPPiRatioPhenix, "PHENIX, #sqrt{s} = 200 GeV","lpf");
1860 // TLegend * l = new TLegend(0.206376, 0.77972, 0.600671, 0.909091);
1866 gSystem->ProcessEvents();
1867 c2->Print(TString(c2->GetName()) + ".png");
1868 c2->Print(TString(c2->GetName()) + ".eps");
1870 gSystem->ProcessEvents();
1871 c3->Print(TString(c3->GetName()) + ".png");
1872 c3->Print(TString(c3->GetName()) + ".eps");
1878 void DrawWithSyst() {
1880 // Draws detector and combined with syst errors.
1882 for(Int_t idet = 0; idet < kNHist; idet++){
1883 if(idet > kITS && idet != iCombInStudy) continue;
1884 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1885 TCanvas * c = new TCanvas(TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge],TString("cWithSyst")+detFlag[idet]+chargeFlag[icharge]);
1886 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.0005,5);
1887 hempty->SetXTitle("p_{t} (GeV/c)");
1888 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
1889 hempty->GetYaxis()->SetTitleOffset(1.35);
1890 hempty->GetXaxis()->SetTitleOffset(1.1);
1893 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1894 // cout << detFlag[idet] << " " << chargeFlag[icharge] << " " << partFlag[ipart] << endl;
1896 TString opt = ipart ? "" : "same";
1897 TH1F * hsyst = new TH1F(*htemplate);
1898 AliBWTools::GetValueAndError(hsyst,hSpectra[idet][ipart][icharge],hSystError[idet][ipart][icharge],kTRUE);
1899 hsyst->SetFillColor(kGray);
1900 hsyst->Draw("e5,same");
1901 hSpectra[idet][ipart][icharge]->Draw("same");
1902 // Do not draw syst error outside of spectra range
1903 Float_t lowEdge = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1904 Float_t hiEdge = AliBWTools::GetHighestNotEmptyBinEdge(hSpectra[idet][ipart][icharge]);
1905 Int_t lowBin = hSystError[idet][ipart][icharge]->FindBin(lowEdge);
1906 Int_t hiBin = hSystError[idet][ipart][icharge]->FindBin(hiEdge);
1907 Int_t nbin = hSystError[idet][ipart][icharge]->GetNbinsX();
1908 for(Int_t ibin = 1; ibin < lowBin; ibin++){
1909 hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1910 hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1912 for(Int_t ibin = hiBin; ibin <= nbin; ibin++){
1913 hSystError[idet][ipart][icharge]->SetBinContent(ibin,0);
1914 hSystError[idet][ipart][icharge]->SetBinError(ibin,0);
1916 // Draw relative syst error, convert to graph to avoid vertical bars at the edges (empty bins in histo)
1917 AliBWTools::GetGraphFromHisto(hSystError[idet][ipart][icharge])->Draw("LX");
1922 PrintCanvas(c,"png,eps");
1929 cout << "Macro: void CombineSpectra(Int_t analysisType=kDoFits, Int_t fitFuncID = kFitLevi) " << endl;
1932 cout << "Possible Arguments" << endl;
1933 cout << "- analysisType:" << endl;
1934 cout << " kDoFits: Fit Combined Spectra " << endl;
1935 cout << " kDoRatios: Particle ratios K/pi and p/pi" << endl;
1936 cout << " kDoSuperposition: Compare different detectors (superimpose and ratios)" << endl;
1937 cout << " kDoCompareStar: Compare combined spectra to star results" << endl;
1938 cout << " kDoDrawWithModels: Compare combined spectra and models" << endl;
1939 cout << " kDoDrawSyst: Draws spectra from individual detectors with their systematic error" << endl;
1940 cout << " kDoCompareToAllCharged: convert to dn/deta/dpt, sum up and compare to all charged" << endl;
1941 cout << " kDoHelp: This help" << endl;
1942 cout << "- fitFuncID, function used to extrapolate and compute yields" << endl;
1943 cout << " An analitic fit function [kFitLevi, kFitUA1, kFitPowerLaw]" << endl;
1944 cout << " Or a shape from a MC moder [kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0]" << endl;
1945 cout << " Which is fitted to the data at low pt and used to extrapolate at low pt" << endl;
1950 void PrintCanvas(TCanvas* c,const TString formats) {
1951 // print a canvas in every of the given comma-separated formats
1952 // ensure the canvas is updated
1953 if(!doPrint) return;
1955 gSystem->ProcessEvents();
1956 TObjArray * arr = formats.Tokenize(",");
1957 TIterator * iter = arr->MakeIterator();
1958 TObjString * element = 0;
1959 TString name =c ->GetName();
1960 name.ReplaceAll(" ","_");
1961 name.ReplaceAll("+","Plus");
1962 name.ReplaceAll("-","");
1963 while ((element = (TObjString*) iter->Next())) {
1964 c->Print(name+ "."+element->GetString());
1971 gSystem->Load("libTree.so");
1972 gSystem->Load("libVMC.so");
1973 gSystem->Load("libMinuit.so");
1974 gSystem->Load("libSTEERBase.so");
1975 gSystem->Load("libESD.so");
1976 gSystem->Load("libAOD.so");
1977 gSystem->Load("libANALYSIS.so");
1978 gSystem->Load("libANALYSISalice.so");
1979 gSystem->Load("libCORRFW.so");
1980 gSystem->Load("libPWG2spectra.so");