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