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