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