]> 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
From C.Loizides: Handle manual branch loading if requested by the user (via SetAutoBr...
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / Fit / CombineSpectra.C
CommitLineData
4c0d7bc7 1#include <iostream>
2#include "TH1F.h"
3#include "TGraphErrors.h"
4#include "AliBWFunc.h"
5#include "AliBWTools.h"
6#include "TF1.h"
7#include "TFile.h"
8#include "TDatabasePDG.h"
9#include "TROOT.h"
10#include "TCanvas.h"
11#include "TFolder.h"
12#include "TStyle.h"
13#include "AliLatexTable.h"
14#include "TLegend.h"
15#include "TVirtualFitter.h"
16#include "TMath.h"
17#include "TH2F.h"
18#include "TSystem.h"
19#include "TLine.h"
20#include "TLatex.h"
21#include "TMath.h"
22
23#include "TASImage.h"
24#include "TPaveText.h"
25#include "ALICEWorkInProgress.C"
26#include "StarPPSpectra.C"
27#include "GetE735Ratios.C"
28
29using namespace std;
30
31
32// A bunch of useful enums and constants
33enum {kPion=0,kKaon,kProton,kNPart};
ee08b77d 34enum {kTPC=0,kTOF,kITS,kITSTPC,kK0,kKinks,kCombTOFTPC,kCombAll,kNHist};// "k0" listed here as a kind of PID method...
4c0d7bc7 35const Int_t kNDet = kITS+2;
36enum {kPos=0,kNeg,kNCharge};
37enum {kPhojet=0,kPyTuneAtlasCSC, kPyTuneCMS6D6T, kPyTunePerugia0, kNTunes} ;
38enum {kFitLevi=0, kFitUA1, kFitPowerLaw,
39 kFitPhojet, kFitAtlasCSC, kFitCMS6D6T, kFitPerugia0,
40 kNFit};
41
42
43// flags, labels and names
44const char * partFlag[] = {"Pion", "Kaon", "Proton"};
b7018bc7 45const char * detFlag[] = {"TPC", "TOF", "ITS", "ITS Global", "K0", "Kinks", "Combined TOF + TPC", "Combined TOF + TPC + ITS"};
4c0d7bc7 46const char * chargeFlag[] = {"Pos", "Neg"};
47const char * chargeLabel[] = {"Positive", "Negative"};
48const char * partLabel[kNPart][kNCharge] = {{"#pi^{+}", "#pi^{-}"},
49 // {"K^{+} (#times 2)", "K^{-} (#times 2)"},
50 {"K^{+}", "K^{-}"},
51 {"p" , "#bar{p}"}};
52const char * mcTuneName[] = {"Phojet",
53 "Pythia - CSC 306",
54 "Pythia - D6T 109",
55 "Pythia - Perugia0 - 320", };
56const char * funcName[] = { "Levi", "UA1", "PowerLaw", "Phojet", "AtlasCSC", "CMS6D6T", "Perugia0"};
57
58// Style
59//const Int_t marker[] = {25,24,28,20,21}; // set for highlithining marek
60const Int_t marker[] = {20,24,25,28,21}; // standard set
61const Int_t color [] = {1,2,4};
62
63const Int_t mcLineColor[] = {kGreen,kRed,kBlue,kBlack};
64const Int_t mcLineStyle[] = {1,2,3,4};
65
66
67// Template needed to combine different detectors
ee08b77d 68const 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};
69Int_t nbinsTempl=34;
4c0d7bc7 70
71TH1F * htemplate = new TH1F("htemplate", "htemplate",nbinsTempl, templBins );
72
73// Globals
74TH1F * hSpectra[kNHist][kNPart][kNCharge];
75TH1F * hSpectraMC[kNTunes][kNPart][kNCharge];
76Double_t mass[kNPart];
77
78// Functions:
79// Loading
80void LoadSpectra() ;
81void LoadMC() ;
82
83// Additional tasks (may be uncommented below)
84void DrawStar(Int_t icharge);
85void GetITSResiduals();
86void DrawWithModels() ;
87void DrawAllAndKaons();
88void DrawWithJacek();
89void DrawRatioToStar();
90void DrawRatios();
91
92// External stuff
93//void ALICEWorkInProgress(TCanvas *c,TString today, TString label);
94
95// Used to tag plots
96TDatime dt;
97TString today = "";
98
99
100
101// Switches
102Bool_t convertToMT = 0;
b7018bc7 103Bool_t doPrint = 1;
4c0d7bc7 104Int_t fitFuncID = kFitLevi;
105Bool_t scaleKaons = kFALSE;
106Bool_t correctSecondaries = 1;
107Bool_t correctGeantFlukaXS = 1;
ee08b77d 108Int_t iCombInStudy = kCombAll;
4c0d7bc7 109
110void CombineSpectra() {
111
112 // This macro is used to combine the 900 GeV spectra from 2009
113
114 // Load stuff
115 gSystem->Load("libTree.so");
116 gSystem->Load("libVMC.so");
117 gSystem->Load("libMinuit.so");
118 gSystem->Load("libSTEERBase.so");
119 gSystem->Load("libESD.so");
120 gSystem->Load("libAOD.so");
121 gSystem->Load("libANALYSIS.so");
122 gSystem->Load("libANALYSISalice.so");
123 gSystem->Load("libCORRFW.so");
124 gSystem->Load("libPWG2spectra.so");
125
126 // Set date
127 today = today + long(dt.GetDay()) +"/" + long(dt.GetMonth()) +"/"+ long(dt.GetYear());
128
129
130 // Set Masses
131 mass[kPion] = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
132 mass[kKaon] = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
133 mass[kProton] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
134
135 // Load histos
136 LoadSpectra();
137 LoadMC();
138
139 // Additional tasks
140 //DrawStar(icharge);
141 // GetITSResiduals();
142 //DrawAllAndKaons();
143 // DrawWithModels() ;
144 //DrawWithJacek();
145 //DrawRatioToStar();
b7018bc7 146 DrawRatios();
147 return;
4c0d7bc7 148
149
150 // Draw combined & Fit
151 AliBWFunc * fm = new AliBWFunc;
152 fm->SetVarType(AliBWFunc::kdNdpt);
153 if (convertToMT) fm->SetVarType(AliBWFunc::kOneOverMtdNdmt);
154
155 // table to print results
156 AliLatexTable table(10,"c|cccccccc");
157 if (fitFuncID == kFitLevi) {
158 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");
159 } else if (fitFuncID == kFitPowerLaw) {
160 table.InsertCustomRow("Part & Yield & Norm & n & pt0 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle & \\langle p_{t}^{2} \\rangle");
161 } else {
162 table.InsertCustomRow("Part & Yield & Par0 & Par2 & Par1 & $\\Chi^2$/NDF & Min X & Frac Above & \\langle p_{t} \\rangle & \\langle p_{t}^{2} \\rangle");
163
164 }
165
166 AliLatexTable tempTable(4,"c|ccc");
167 tempTable.InsertCustomRow("Part & Yield & Yield Below & Frac Above\\\\");
168
169
170 // Fit all
171 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
172
173 TCanvas * c2 = new TCanvas(TString("cCombined")+chargeFlag[icharge]+"_"+funcName[fitFuncID], TString("cCombined")+chargeFlag[icharge]);
174 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,2.9, 100, 0.001,5);
175 hempty->SetXTitle("p_{t} (GeV/c)");
176 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
177 hempty->Draw();
178 // DrawStar(icharge);
179 c2->SetLogy();
180 TLegend * l = new TLegend(0.516779, 0.729021, 0.89094 ,0.916084, chargeLabel[icharge]);
181 l->SetFillColor(kWhite);
182
183 for(Int_t ipart = 0; ipart < kNPart; ipart++){
184 Float_t fitmin = 0;
185 Float_t fitmax = 3;
186
187 // Get functions
188 TF1 * func = 0;
189 if(fitFuncID == kFitLevi) {
190 if (ipart == kPion)
191 func = fm->GetLevi(mass[ipart], 0.12, 7, 1.5);
192 if (ipart == kKaon)
193 func = fm->GetLevi(mass[ipart], 0.17, 7, 0.17);
194 if (ipart == kProton)
195 func = fm->GetLevi(mass[ipart], 0.15, 8, 0.08);
196 }
197 else if(fitFuncID == kFitUA1) func = fm->GetUA1(mass[ipart],0.2,1.25,1.25,0.2,1.5);
198 else if(fitFuncID == kFitPowerLaw) {
199 if (ipart == kPion)
200 func = fm->GetPowerLaw(1.0, 8.6, 7);
201 if (ipart == kKaon)
202 func = fm->GetPowerLaw(3.0, 12, 2.6);
203 if (ipart == kProton)
204 func = fm->GetPowerLaw(24, 72, 0.8);
205 }
206 else if(fitFuncID == kFitPhojet) func = fm->GetHistoFunc(hSpectraMC[kPhojet] [ipart][icharge]);
207 else if(fitFuncID == kFitAtlasCSC) func = fm->GetHistoFunc(hSpectraMC[kPyTuneAtlasCSC][ipart][icharge]);
208 else if(fitFuncID == kFitPerugia0) func = fm->GetHistoFunc(hSpectraMC[kPyTunePerugia0][ipart][icharge]);
209 else if(fitFuncID == kFitCMS6D6T) func = fm->GetHistoFunc(hSpectraMC[kPyTuneCMS6D6T] [ipart][icharge]);
210 else {
211 cout << "Unknown fit ID " << fitFuncID << endl;
212 return;
213 }
214
215 if(fitFuncID >= kFitPhojet){
216 fitmin = 0.0;
217 fitmax = 1.0;
218 }
219
ee08b77d 220 if(!AliBWTools::Fit(hSpectra[iCombInStudy][ipart][icharge],func,fitmin,fitmax)) {
4c0d7bc7 221 cout << " FIT ERROR " << endl;
222 return;
223 }
ee08b77d 224 hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
225 hSpectra[iCombInStudy][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
226 ((TF1*)hSpectra[iCombInStudy][ipart][icharge]->GetListOfFunctions()->At(0))->SetRange(0,4);
227 l->AddEntry(hSpectra[iCombInStudy][ipart][icharge],
4c0d7bc7 228 scaleKaons && ipart == kKaon ?
229 (TString(partLabel[ipart][icharge])+" #times 2").Data()
230 : partLabel[ipart][icharge]);
ee08b77d 231// TF1 * fClone = (TF1*) hSpectra[iCombInStudy][ipart][icharge]->GetListOfFunctions()->At(0)->Clone();
232// hSpectra[iCombInStudy][ipart][icharge]->GetListOfFunctions()->Add(fClone);
4c0d7bc7 233// fClone->SetLineStyle(kDashed);
234// fClone->SetRange(0,100);
235// fClone->Draw("same");
236
237 // populate table
238 // Float_t yield = func->Integral(0.45,1.05);
239 // Float_t yieldE = func->IntegralError(0.45,1.05);
240
241 Float_t yield = func->Integral(0.,100);
242 //Float_t yieldE = func->IntegralError(0.,100);
243
244 Double_t yieldTools, yieldETools;
245 Double_t partialYields[3],partialYieldsErrors[3];
ee08b77d 246 AliBWTools::GetYield(hSpectra[iCombInStudy][ipart][icharge], func, yieldTools, yieldETools,
4c0d7bc7 247 0, 100, partialYields,partialYieldsErrors);
248 Double_t tslope = func->GetParameter(2);
249 Double_t tslopeE = func->GetParError(2);
250
251 table.SetNextCol(partLabel[ipart][icharge]);
252 //table.SetNextCol(yield,yieldE,-4);
253 table.SetNextCol(yieldTools, yieldETools,-4);
254 table.SetNextCol(func->GetParameter(0));
255 table.SetNextCol(tslope,tslopeE,-4);
256 table.SetNextCol(func->GetParameter(1));
257 table.SetNextCol(Form("%2.2f/%d",func->GetChisquare(),func->GetNDF()));
ee08b77d 258 Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[iCombInStudy][ipart][icharge]);
4c0d7bc7 259 //Float_t lowestPoint = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kITS][ipart][icharge]);
260 Float_t yieldAbove = func->Integral(lowestPoint,100);
261 table.SetNextCol(lowestPoint,-2);
262 table.SetNextCol(yieldAbove/yield,-2);
263 Float_t mean, meane;
264 Float_t mean2, mean2e;
265 AliBWTools::GetMean (func, mean, meane );
266 AliBWTools::GetMeanSquare(func, mean2, mean2e);
267 table.SetNextCol(mean, meane ,-4);
268 table.SetNextCol(mean2, mean2e,-4);
269
270 // fMean2->IntegralError(0,100)/func->Integral(0,100),-7);
271 table.InsertRow();
272
273
274 /// TEMP TABLE
275 tempTable.SetNextCol(partLabel[ipart][icharge]);
276 tempTable.SetNextCol(yieldTools, yieldETools, -4);
277 tempTable.SetNextCol(partialYields[1], partialYieldsErrors[1], -4);
278 tempTable.SetNextCol(yieldAbove/yield,-2);
279 tempTable.InsertRow();
280 }
281 l->Draw();
282 if (doPrint) {
283 c2->Update();
284 gSystem->ProcessEvents();
285 c2->Print(TString(c2->GetName()) + ".eps");
286 ALICEWorkInProgress(c2,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
287 c2->Print(TString(c2->GetName()) + ".png");
288 }
289
290 }
291
292
293 table.PrintTable("ASCII");
294
295 cout << "" << endl;
296 tempTable.PrintTable("ASCII");
297
298}
299
300void LoadSpectra() {
301
302 TFile * f=0;
303
304 // Load
305
306
307 // TOF
308 // f = new TFile("./Files/spectra-pos-y.root");
309 f = new TFile("./Files/spectra-pos-y_20100615.root");
310 hSpectra[kTOF][kPion] [kPos]= (TH1F*) f->Get("hpi");
311 hSpectra[kTOF][kProton][kPos]= (TH1F*) f->Get("hpr");
312 hSpectra[kTOF][kKaon] [kPos]= (TH1F*) f->Get("hka");
313 hSpectra[kTOF][kPion] [kPos]->SetName("hpiPos");
314 hSpectra[kTOF][kProton][kPos]->SetName("hprPos");
315 hSpectra[kTOF][kKaon] [kPos]->SetName("hkaPos");
316 //f = new TFile("./Files/spectra-neg-y.root");
317 f = new TFile("./Files/spectra-neg-y_20100615.root");
318 hSpectra[kTOF][kPion] [kNeg]= (TH1F*) f->Get("hpi");
319 hSpectra[kTOF][kProton][kNeg]= (TH1F*) f->Get("hpr");
320 hSpectra[kTOF][kKaon] [kNeg]= (TH1F*) f->Get("hka");
321 hSpectra[kTOF][kPion] [kNeg]->SetName("hpiNeg");
322 hSpectra[kTOF][kProton][kNeg]->SetName("hprNeg");
323 hSpectra[kTOF][kKaon] [kNeg]->SetName("hkaNeg");
324
325
326 // Clean UP TPC spectra, removing unwanted points
327 cout << "Cleaning Up TOF spectra" << endl;
328 Int_t nbin = hSpectra[kTOF][kKaon][kPos]->GetNbinsX();
329 for(Int_t ibin = 1; ibin <= nbin; ibin++){
330 Float_t pt = hSpectra[kTOF][kKaon][kPos]->GetBinCenter(ibin);
331 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
332 if (pt > 2.35) {
333 hSpectra[kTOF][kKaon][icharge]->SetBinContent(ibin,0);
334 hSpectra[kTOF][kKaon][icharge]->SetBinError (ibin,0);
335 hSpectra[kTOF][kProton][icharge]->SetBinContent(ibin,0);
336 hSpectra[kTOF][kProton][icharge]->SetBinError (ibin,0);
337 }
338 }
339 }
340// cout << "Scaling TOF to TPC" << endl;
341// // Scale protons to TPC
342// hSpectra[kTOF][kProton][kPos]->Scale(1./1.05);
343// // Scale pbar so that pbar/p is compatible with Panos
344// hSpectra[kTOF][kProton][kNeg]->Scale(1./1.1);
345
346
347
348
349
350 // ITS SA (Emanuele)
351 // f = new TFile("./Files/ITSsaSPECTRA_3clusters20100619.root");
352 f = new TFile("./Files/ITSsaSPECTRAperMICHELE_20100703.root");
353 hSpectra[kITS][kPion] [kPos]= (TH1F*) f->Get("hSpectraPos0");
354 hSpectra[kITS][kKaon] [kPos]= (TH1F*) f->Get("hSpectraPos1");
355 hSpectra[kITS][kProton][kPos]= (TH1F*) f->Get("hSpectraPos2");
356 hSpectra[kITS][kPion] [kNeg]= (TH1F*) f->Get("hSpectraNeg0");
357 hSpectra[kITS][kKaon] [kNeg]= (TH1F*) f->Get("hSpectraNeg1");
358 hSpectra[kITS][kProton][kNeg]= (TH1F*) f->Get("hSpectraNeg2");
359
360 for(Int_t ipart = 0; ipart < kNPart; ipart++){
361 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
362 Int_t nbin = hSpectra[kITS][ipart][icharge]->GetNbinsX();
363 // hSpectra[kITS][ipart][icharge]->Scale(276004.);// Emanule divided his spectra...
364 for(Int_t ibin = 1; ibin <= nbin; ibin++){
365 if(hSpectra[kITS][ipart][icharge]->GetBinContent(ibin) < 0 ) {
366 hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
367 hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
368 }
369// if ((ipart == kKaon && ibin >= 12) || (ipart == kProton && ibin >= 20)) {
370// hSpectra[kITS][ipart][icharge]->SetBinContent(ibin,0);
371// hSpectra[kITS][ipart][icharge]->SetBinError (ibin,0);
372// }
373 }
374
375 }
376 }
377
378
379 // ITS + TPC (Marek)
380 f = TFile::Open("./Files/SpectraCorrectedITSBeforeProtons_20100618.root");
381 TList * list = (TList*) gDirectory->Get("output");
382 hSpectra[kITSTPC][kPion] [kPos]= (TH1F*) list->FindObject("Pions");
383 hSpectra[kITSTPC][kKaon] [kPos]= (TH1F*) list->FindObject("Kaons");
384 hSpectra[kITSTPC][kProton][kPos]= (TH1F*) list->FindObject("Protons");
385 hSpectra[kITSTPC][kPion] [kNeg]= (TH1F*) list->FindObject("AntiPions");
386 hSpectra[kITSTPC][kKaon] [kNeg]= (TH1F*) list->FindObject("AntiKaons");
387 hSpectra[kITSTPC][kProton][kNeg]= (TH1F*) list->FindObject("AntiProtons");
388
389 // TPC
390 // htemplate = hSpectra[kITS][kProton][kNeg]; //FIXME
391 f = new TFile("./Files/protonSpectra_20100615.root");
392 hSpectra[kTPC][kProton][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonPosClassic"),htemplate);
393 hSpectra[kTPC][kProton][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("protonNegClassic"),htemplate);
394 f = new TFile("./Files/pionSpectra_20100615.root");
395 hSpectra[kTPC][kPion][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionPosClassic"),htemplate);
396 hSpectra[kTPC][kPion][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("pionNegClassic"),htemplate);
397 f = new TFile("./Files/kaonSpectra_20100615.root");
398 hSpectra[kTPC][kKaon][kPos]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonPosClassic"),htemplate);
399 hSpectra[kTPC][kKaon][kNeg]= AliBWTools::GetHistoFromGraph((TGraphErrors*)f->Get("kaonNegClassic"),htemplate);
400
401 // Clean UP TPC spectra, removing unwanted points
402 cout << "Cleaning Up TPC spectra" << endl;
403 nbin = hSpectra[kTPC][kKaon][kPos]->GetNbinsX();
404 for(Int_t ibin = 0; ibin < nbin; ibin++){
405 Float_t pt = hSpectra[kTPC][kKaon][kPos]->GetBinCenter(ibin);
406 if (pt > 0.45) {
407 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
408 hSpectra[kTPC][kKaon][icharge]->SetBinContent(ibin,0);
409 hSpectra[kTPC][kKaon][icharge]->SetBinError (ibin,0);
410 }
411 }
412 if (pt < 0.45) {
413 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
414 hSpectra[kTPC][kProton][icharge]->SetBinContent(ibin,0);
415 hSpectra[kTPC][kProton][icharge]->SetBinError (ibin,0);
416 }
417 }
418 }
419
420
421
422
423
424 // K0s
425 f = new TFile ("./Files/PtSpectraCorrectedK0sOff-dNdy-Jun02.root");
426 // hSpectra[kK0][kKaon][kPos] = (TH1F*) AliBWTools::GetdNdPtFromOneOverPt((TH1*) gDirectory->Get("hSpectraOff"));
427 hSpectra[kK0][kKaon][kPos] = (TH1F*) gDirectory->Get("hSpectraOff");
428 // hSpectra[kK0][kKaon][kPos]->Scale(2*TMath::Pi());
429 // hSpectra[kK0][kKaon][kPos]->Scale(1./272463);
430 hSpectra[kK0][kKaon][kNeg] = hSpectra[kK0][kKaon][kPos];
431
432 // Kinks: TO BE FIXED WITH POSITIVES AND NEGATIVES
433 // f = new TFile ("./Files/PtAllKaonKinkRap6Apr24.root");
434 f = new TFile ("./Files/PtKaonKinkJune13AllPN_20100615.root");
435 hSpectra[kKinks][kKaon][kPos] = (TH1F*)gDirectory->Get("fptallKPA");
436 hSpectra[kKinks][kKaon][kNeg] = (TH1F*)gDirectory->Get("fptallKNA");
437
438
439 // Apply correction factrs
440 // Secondaries for protons
441 f = new TFile ("./Files/corrFactorProtons_20100615.root");
442 TH1F * hCorrSecondaries = (TH1F*) gDirectory->Get("corrFactorProtons");
443 if(correctSecondaries) {
444 cout << "CORRECTING SECONDARIES" << endl;
445
446 for(Int_t idet = 0; idet <= kTOF; idet++){ // TPC and TOF only
447 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
448 Int_t ipart = kProton;
449 TH1* h = hSpectra[idet][ipart][icharge]; // lighter notation below
450 if (h){
451 Int_t nbins = h->GetNbinsX();
452 for(Int_t ibin = 0; ibin < nbins; ibin++){
453 // Float_t pt = convertToMT ? TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) : h->GetBinCenter(ibin);
454 Float_t pt = h->GetBinCenter(ibin);
455 if (icharge == kNeg) pt = -pt;
456 Int_t binCorrection = hCorrSecondaries->FindBin(pt);
457 Float_t correction = hCorrSecondaries->GetBinContent(binCorrection);
458 // cout << pt << " " << correction << endl;
459
460 if (correction != 0) {// If the bin is empty this is a 0
461 h->SetBinContent(ibin,h->GetBinContent(ibin)/correction);
462 h->SetBinError (ibin,h->GetBinError (ibin)/correction);
463 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
464 cout << "Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
465 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
466
467 }
468 }
469 }
470 }
471 }
472 }
473 // geant/fluka absorption
474 if(correctGeantFlukaXS) {
475 cout << "CORRECTING GEANT3/FLUKA" << endl;
476
477 f = new TFile ("./Files/correctionForCrossSection_20100615.root");
478 TH2D * hCorrFluka[kNCharge];
479 hCorrFluka[kPos] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionProtons");
480 hCorrFluka[kNeg] = (TH2D*)gDirectory->Get("gHistCorrectionForCrossSectionAntiProtons");
481 for(Int_t idet = 0; idet < kNDet; idet++){
482 if (idet == kITS) continue; // skip ITS sa
483 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
484 Int_t ipart = kProton;
485 TH1 * h = hSpectra[idet][ipart][icharge]; // lighter notation below
486 if (h){
487 Int_t nbins = h->GetNbinsX();
488 for(Int_t ibin = 0; ibin < nbins; ibin++){
489// Float_t pt = convertToMT ?
490// TMath::Sqrt(h->GetBinCenter(ibin)*h->GetBinCenter(ibin)-mass[kProton]*mass[kProton]) :
491// h->GetBinCenter(ibin);
492 Float_t pt = h->GetBinCenter(ibin);
493 Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
494 Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(h->GetNbinsY()+1);
495 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
496 if (pt > maxPtCorrection) pt = maxPtCorrection;
497 Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
498 if (idet == kTOF && icharge == kNeg) {
499
500 // Apply parametrized correction computed by francesco
501 // Fitted panos correction and using momentum at the outer radius of the TPC
502
503 Float_t ptav = pt; // Just to use the same name francesco uses...
504
505 // from pT constrained at P.V. (ptav) to pT TPC outer (ptTPCout)
506 Float_t ptTPCout=ptav*(1-6.81059e-01*TMath::Exp(-ptav*4.20094));
507
508 // traking correction (fit to Panos)
509 Float_t antiprotonEC = 1 - 0.129758 *TMath::Exp(-ptav*0.679612);
510
511 // TOF matching efficiency correction (derived from Panos one scaled for M.B.(TOF)/M.B.(TPC)).
512 Float_t antiprotonEC2 = TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCout*0.679612),0.07162/0.03471);
513 correction = antiprotonEC * antiprotonEC2;
514 }
515 // cout << icharge<< " " << h->GetBinCenter(ibin) << " " << pt << " " << correction << endl;
516 if (correction != 0) {// If the bin is empty this is a 0
517 h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
518 h->SetBinError (ibin,h->GetBinError (ibin)*correction);
519// if (idet == kTOF) {
520// cout << "CORRECTING TOF TWICE" << endl;
521// h->SetBinContent(ibin,h->GetBinContent(ibin)*correction);
522// h->SetBinError (ibin,h->GetBinError (ibin)*correction);
523// }
524 } else if (h->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
525 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, " << detFlag[idet] << " " << chargeFlag[icharge] << endl;
526 cout << " Bin content: " << h->GetBinContent(ibin) << endl;
527 }
528
529 }
530
531 }
532 }
533 }
534 }
535
536
537 // Set style and scale
538 for(Int_t idet = 0; idet < kNDet; idet++){
539 for(Int_t ipart = 0; ipart < kNPart; ipart++){
540 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
541 if (hSpectra[idet][ipart][icharge]){
542 hSpectra[idet][ipart][icharge]->SetStats(0); // disable stats
543 hSpectra[idet][ipart][icharge]->SetMarkerColor (color[ipart] );
544 hSpectra[idet][ipart][icharge]->SetLineColor (color[ipart] );
545 hSpectra[idet][ipart][icharge]->SetMarkerStyle (marker[ipart]);
546 hSpectra[idet][ipart][icharge]->SetXTitle("p_{t} (GeV/c)");
547 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} dN/dp_{t} (GeV/c)^{-1}");
548 if (convertToMT) {
549 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[idet][ipart][icharge]);
550 hSpectra[idet][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
551 hSpectra[idet][ipart][icharge]->SetXTitle("m_{t} (GeV)");
552 hSpectra[idet][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
553 }
554// if (idet == kTOF || idet == kTPC) {
555// hSpectra[idet][ipart][icharge]->Scale(1./236763);
556// }
557// if (idet == kITS){
558// hSpectra[idet][ipart][icharge]->Scale(202945./236763);
559// }
560 if (scaleKaons && ipart == kKaon) {
561 hSpectra[idet][ipart][icharge]->Scale(2.);
562 }
563 } else {
564 printf("Cannot load %s,%s,%s\n",detFlag[idet],partFlag[ipart],chargeFlag[icharge]);
565 }
566 }
567 }
568 }
569
570
571
572
ee08b77d 573
574 // Create fake weights for the mean; To be update once we have syst errors
575 TH1F * hWeights[3];
576 const Double_t kWeights[3] = {0.1,0.1,0.2}; // TPC, TOF, ITS
577 for(Int_t idet = 0; idet <= kITS ; idet++){
578 hWeights[idet] = (TH1F*) hSpectra[idet][kPion][kPos]->Clone();
579 Int_t nbin = hWeights[idet]->GetNbinsX();
580 for(Int_t ibin = 0; ibin < nbin; ibin++){
581 hWeights[idet]->SetBinError(ibin, kWeights[idet]);
582 }
583 }
584
585
586
4c0d7bc7 587 // Combine detectors
588 for(Int_t ipart = 0; ipart < kNPart; ipart++){
589 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
590 TH1F * htemplLocal = htemplate; // If we are converting to 1/mt dNdmt we need to convert the template as well...
591
592 if (convertToMT) {
593 TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(htemplate);
594 htemplLocal = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
595
596 }
597 hSpectra[kCombTOFTPC][ipart][icharge] = AliBWTools::CombineHistos(hSpectra[kTOF][ipart][icharge],
598 hSpectra[kTPC][ipart][icharge],
599 htemplLocal,1.);;
ee08b77d 600 hSpectra[kCombAll][ipart][icharge] =
601 AliBWTools::Combine3HistosWithErrors(hSpectra[kITS][ipart][icharge],
602 hSpectra[kTPC][ipart][icharge],
603 hSpectra[kTOF][ipart][icharge],
604 hWeights[kITS],
605 hWeights[kTPC],
606 hWeights[kTOF],
607 htemplLocal,1);;
4c0d7bc7 608// if (convertToMT) {
609// TH1F * htmp = (TH1F*) AliBWTools::GetOneOverPtdNdPt(hSpectra[kCombTOFTPC][ipart][icharge]);
610// hSpectra[kCombTOFTPC][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
611// hSpectra[kCombTOFTPC][ipart][icharge]->SetXTitle("m_{t} (GeV)");
612// hSpectra[kCombTOFTPC][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
613// }
614 }
615 }
616
617
618 // Scale for the number of inelastic collisions and correct for
619 // efficiency losses due to physics selection:
620
621 Double_t effPhysSel[kNPart];
622 effPhysSel[kPion] = 1.012;
623 effPhysSel[kKaon] = 1.013;
624 effPhysSel[kProton] = 1.014;
625
626
627 for(Int_t idet = 0; idet < kNHist; idet++){
628 for(Int_t ipart = 0; ipart < kNPart; ipart++){
629 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
630 if(hSpectra[idet][ipart][icharge]) {
631 // cout << "Scaling!" << endl;
632 hSpectra[idet][ipart][icharge]->Scale(1.*effPhysSel[ipart]/278366.15); // Scale PhysSel tutti? // FIXME
633 }
634 }
635 }
636 }
637
638
639}
640
641void LoadMC() {
642
643 TFile * f = 0;
644 const char * evClass= "INEL";
645 const char * files[] = {"./Files/dndeta_Phojet_10M_900GeV.root",
646 "./Files/dndeta_AtlasCSC306_10M_900GeV.root",
647 "./Files/dndeta_CMSD6T109_10M_900GeV.root",
648 "./Files/dndeta_Perugia0320_10M_900GeV.root", };
649
650 // Phojet
651 for(Int_t itune = 0; itune < kNTunes; itune++){
652 f = new TFile(files[itune]);
653 for(Int_t ipart = 0; ipart < kNPart; ipart++){
654 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
655 hSpectraMC[itune][ipart][icharge] = (TH1F*) f->Get(Form("fHistPtID_%s_%s%s",evClass,partFlag[ipart],icharge==0 ? "Pos" : "Neg"));
656 }
657 }
658 }
659
660
661 // Set style
662 for(Int_t itune = 0; itune < kNTunes; itune++){
663 for(Int_t ipart = 0; ipart < kNPart; ipart++){
664 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
665 if (hSpectraMC[itune][ipart][icharge]){
666 hSpectraMC[itune][ipart][icharge]->SetName(TString(hSpectraMC[itune][ipart][icharge]->GetName())+mcTuneName[itune]);
667 hSpectraMC[itune][ipart][icharge]->SetMarkerColor (mcLineColor[itune] );
668 hSpectraMC[itune][ipart][icharge]->SetLineColor (mcLineColor[itune] );
669 hSpectraMC[itune][ipart][icharge]->SetLineStyle (mcLineStyle[itune] );
670 hSpectraMC[itune][ipart][icharge]->SetMarkerStyle (1);
671 if (convertToMT) {
672 TH1F * htmp = (TH1F*)AliBWTools::GetOneOverPtdNdPt(hSpectraMC[itune][ipart][icharge]);
673 hSpectraMC[itune][ipart][icharge] = (TH1F*)AliBWTools::GetdNdmtFromdNdpt(htmp,mass[ipart]);
674 hSpectraMC[itune][ipart][icharge]->SetXTitle("m_{t} (GeV)");
675 hSpectraMC[itune][ipart][icharge]->SetYTitle("1/N_{ev} 1/m_{t} dN/dp_{t} (GeV)^{-1}");
676 }
677
678 } else {
679 printf("Cannot load MC # %d,%s,%s\n",itune,partFlag[ipart],chargeFlag[icharge]);
680 }
681 }
682 }
683 }
684
685}
686
687void DrawStar(Int_t icharge) {
688 // cout << icharge << endl;
689
690 // gROOT->LoadMacro("StarPPSpectra.C");
691 TGraphErrors ** gStar = StarPPSpectra();
692
693 for(Int_t istar = 0; istar < 6; istar++){
694 gStar[istar]->SetMarkerStyle(kOpenStar);
695 if (icharge==kPos && (istar%2) ) gStar[istar]->Draw("P");
696 else if (icharge==kNeg && !(istar%2) ) gStar[istar]->Draw("P");
697 else cout << "Skipping Star " << istar << endl;
698 }
699}
700
701void GetITSResiduals() {
702
703
704 for(Int_t ipart = 0; ipart < kNPart; ipart++){
705 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
706 // cout << "1 " << ipart << " " << icharge << endl;
707
708// gSystem->ProcessEvents();
709// gSystem->Sleep(1000);
710 TF1* f = (TF1*) hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0);
711 TH1F * hres1, *hres2;
712 AliBWTools::GetResiduals(hSpectra[kITS][ipart][icharge], f, &hres1, &hres2);
713
714 TCanvas * c1 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge],TString(partFlag[ipart])+"_"+chargeFlag[icharge]);
715 c1->SetLogy();
716 hSpectra[kCombTOFTPC][ipart][icharge]->Draw();
717 hSpectra[kITS][ipart][icharge]->SetMarkerStyle(24);
718 hSpectra[kCombTOFTPC][ipart][icharge]->SetMarkerStyle(20);
719 hSpectra[kITS][ipart][icharge]->Draw("same");
720 hSpectra[kCombTOFTPC][ipart][icharge]->GetListOfFunctions()->At(0)->Draw("same");
721 TLegend* l = new TLegend( 0.182886, 0.192308, 0.505034,0.384615, TString(partLabel[ipart][icharge])+" "+chargeFlag[icharge]);
722 l->AddEntry(hSpectra[kCombTOFTPC][ipart][icharge], "TOF+TPC");
723 l->AddEntry(hSpectra[kITS][ipart][icharge], "ITS");
724 l->AddEntry(f, "Fit to TOF+TPC");
725 l->Draw();
726
727 TCanvas * c2 = new TCanvas(TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res",
728 TString(partFlag[ipart])+"_"+chargeFlag[icharge]+"_res");
729 c2->SetGridy();
730 hres2->SetMinimum(-1);
731 hres2->SetMaximum(1);
732 hres2->Draw();
733 hres2->GetYaxis()->SetTitleOffset(1.2);
734 Float_t x = AliBWTools::GetLowestNotEmptyBinEdge(hSpectra[kCombTOFTPC][ipart][icharge]);
735 TLine * line = new TLine(x,-1,x,1);
736 line->SetLineStyle(kDashed);
737 line->Draw("same");
738
739 if (doPrint) {
740 c1->Update();
741 c2->Update();
742 gSystem->ProcessEvents();
743 c1->Print(TString(c1->GetName()) + ".png");
744 c2->Print(TString(c2->GetName()) + ".png");
745 }
746 }
747 }
748}
749
750void DrawWithModels() {
751
752 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
753
754 // Draw with models
755 for(Int_t ipart = 0; ipart < kNPart; ipart++){
756 // Multipad canvas
757 TCanvas * c1 = new TCanvas(TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge],TString("cSpectra")+partFlag[ipart]+chargeFlag[icharge] );
758 TPad *p1 = new TPad(TString("p1")+partFlag[ipart]+chargeFlag[icharge], "p1", 0.0, 0.3, 1.0, 0.95, 0, 0, 0);
759 p1->SetBottomMargin(0);
760 p1->Draw();
761
762 TPad *p2 = new TPad(TString("p2")+partFlag[ipart]+chargeFlag[icharge], "p2", 0.0, 0.05, 1.0, 0.3, 0, 0, 0);
763 p2->SetTopMargin(0);
764 p2->SetBottomMargin(0.4);
765 p2->Draw();
766
767 Float_t scaleFonts = (0.95-0.3)/(0.3-0.05);
768
769 // Draw spectra
770 p1->cd();
771 p1->SetLogy();
772 TH2F * hempty = new TH2F(TString("hempty")+long(ipart),"hempty",100,0.,4, 100, 0.0015,5);
773 hempty->SetXTitle("p_{t} (GeV/c)");
774 hempty->SetYTitle("1/N_{ev} d^{2}N / dydp_{t} (GeV/c)^{-1}");
775 hempty->Draw();
776 c1->SetLogy();
777
778
779 TLegend * l =new TLegend( 0.543624, 0.431818, 0.892617,0.629371);
780 l->SetFillColor(kWhite);
ee08b77d 781 hSpectra[iCombInStudy][ipart][icharge]->Draw("same");
4c0d7bc7 782 l->AddEntry(hSpectra[kTOF][ipart][icharge],TString ("Data: ")+partLabel[ipart][icharge]);
783 for(Int_t itune = 0; itune < kNTunes; itune++){
784 l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
785 hSpectraMC[itune][ipart][icharge]->SetLineWidth(2);
786 hSpectraMC[itune][ipart][icharge]->Draw("same,chist");
787 }
788 l->Draw("same");
789
790 // Draw ratio
791 p2->cd();
792 TH2F * hemptyr = new TH2F(TString("hemptyratio")+long(ipart),"hempty",100,0.,4, 100, 0.01,2.99);
793 hemptyr->SetXTitle("p_{t} (GeV/c)");
794 hemptyr->SetYTitle("Data/MC");
795 hemptyr->GetXaxis()->SetLabelSize(0.04*scaleFonts);
796 hemptyr->GetYaxis()->SetLabelSize(0.04*scaleFonts);
797 hemptyr->GetYaxis()->SetTitleSize(0.05*scaleFonts);
798 hemptyr->GetYaxis()->SetTitleOffset(1.4/scaleFonts);
799 hemptyr->GetXaxis()->SetTitleSize(0.05*scaleFonts);
800 hemptyr->SetTickLength(0.03*scaleFonts, "X");
801 hemptyr->SetTickLength(0.02*scaleFonts, "Y");
802 // hemptyr->GetXaxis()->SetTitleOffset(1.4/scaleFonts);
803 hemptyr->GetYaxis()->SetNdivisions(5);
804 hemptyr->Draw("");
805
806 AliBWFunc fm;
807 for(Int_t itune = 0; itune < kNTunes; itune++){
808 TF1* f = fm.GetHistoFunc(hSpectraMC[itune][ipart][icharge], TString("f")+mcTuneName[itune]);
809
810 // l->AddEntry(hSpectraMC[itune][ipart][icharge],mcTuneName[itune]);
ee08b77d 811 TH1F* hRatio = AliBWTools::DivideHistoByFunc(hSpectra[iCombInStudy][ipart][icharge],f);
4c0d7bc7 812 hRatio->SetLineStyle(hSpectraMC[itune][ipart][icharge]->GetLineStyle());
813 hRatio->SetLineColor(hSpectraMC[itune][ipart][icharge]->GetLineColor());
814 hRatio->SetLineWidth(hSpectraMC[itune][ipart][icharge]->GetLineWidth());
815 hRatio->Draw("lhist,same");
816 }
817
818
819 // print
820 if(doPrint) {
821 c1->Update();
822 gSystem->ProcessEvents();
823 c1->Print(TString(c1->GetName())+".eps");
824 ALICEWorkInProgress(c1,"","#splitline{ALICE Preliminary}{Statistical Error Only}");
825 c1->Print(TString(c1->GetName())+".png");
826
827 }
828 }
829 }
830
831
832
833}
834
835void DrawAllAndKaons() {
836
837
838 // gROOT->LoadMacro("ALICEWorkInProgress.C");
839
840 gStyle->SetOptFit(0);
841
ee08b77d 842 TH1F * hKaonsAllTPCTOF = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
843 hKaonsAllTPCTOF->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
4c0d7bc7 844
845 TH1F * hK0Scaled = (TH1F*) hSpectra[kK0][kKaon][kPos]->Clone();
846 hK0Scaled->Add(hSpectra[kK0][kKaon][kPos]);
847
848 hSpectra[kKinks][kKaon][kPos]->SetMarkerStyle(25);
849 hSpectra[kKinks][kKaon][kPos]->SetStats(0);
850 TH1F * hKinksAll = (TH1F*) hSpectra[kKinks][kKaon][kPos]->Clone();
851 hKinksAll->Add(hSpectra[kKinks][kKaon][kNeg]);
852
853 TCanvas * c1 = new TCanvas("cKaons","cKaons");
854 c1->SetLogy();
855 TH2F * hempty = new TH2F("hempty_allkaons","hempty",100,0.,3, 100, 1e-4,10);
856 hempty->SetXTitle("p_{t} (GeV/c)");
857 hempty->SetYTitle("dN / dp_{t} (A.U.)");
858 hempty->Draw();
859 hKaonsAllTPCTOF->Draw("same");
860 hK0Scaled->Draw("same");
861 hKinksAll->Draw("same");
862
863 TLegend * leg = new TLegend(0.2013423,0.2255245,0.5503356,0.4335664,NULL,"brNDC");
864// leg->SetBorderSize(0);
865// leg->SetLineColor(1);
866// leg->SetLineStyle(1);
867// leg->SetLineWidth(1);
868// leg->SetFillColor(19);
869// leg->SetFillStyle(1001);
870 TLegendEntry *entry=leg->AddEntry("hkaPos_h_tpcKaonPos","K^{+} + K^{-}, TPC+TOF ","lpf");
871 entry=leg->AddEntry("h1PtSpectraOff_inv","K^{0} #times 2","lpf");
872 entry=leg->AddEntry("fptallK","K^{+} + K ^{-}, Kinks","lpf");
873 leg->Draw();
874
875 ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not fully corrected}");
876 TLatex * tex = new TLatex(0.2120805,0.01288336,"Statistical error only");
877 tex->SetTextColor(2);
878 tex->SetTextFont(42);
879 tex->SetTextSize(0.03496503);
880 tex->Draw();
881
882 c1->Update();
883 if(doPrint) c1->Print(TString(c1->GetName())+".png");
884
885 // Draw all "stable" hadrons
886 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
887 TCanvas * c1 = new TCanvas(TString("cAll_")+chargeFlag[icharge],TString("cAll_")+chargeFlag[icharge]);
888 c1->SetLogy();
889 TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.,4, 100, 1e-4,10);
890 hempty->SetXTitle("p_{t} (GeV/c)");
891 hempty->SetYTitle("dN / dp_{t} (A.U.)");
892 hempty->Draw();
893 leg = new TLegend( 0.645973, 0.325175, 0.892617,0.636364, NULL,"brNDC");
894 for(Int_t ipart = 0; ipart < kNPart; ipart++) {
895 for(Int_t idet = 0; idet <= kITSTPC; idet++){
896 // if (idet == kITS) continue;
897 // if (idet == kITSTPC) hSpectra[idet][ipart][icharge]->SetMarkerColor(kGreen);
898 hSpectra[idet][ipart][icharge]->SetMarkerStyle(marker[idet]);
899 hSpectra[idet][ipart][icharge]->Draw("same");
900 leg->AddEntry(hSpectra[idet][ipart][icharge],TString(partLabel[ipart][icharge])+" (" + detFlag[idet] + ")","lpf");
901 }
902 // leg->AddLine();
903 }
904 leg->Draw();
905 ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Performance}{Not Fully Corrected}");
906 c1->Update();
907 if(doPrint) c1->Print(TString(c1->GetName())+".png");
908 }
909
910
911 // Draw ratios (tmp)
912
913// new TCanvas;
914// hSpectra[kTPC][kKaon][kNeg]->Divide(hSpectra[kTPC][kKaon][kPos]);
915// hSpectra[kTPC][kKaon][kNeg]->Draw();
916// new TCanvas;
917// hSpectra[kITS][kKaon][kNeg]->Divide(hSpectra[kITS][kKaon][kPos]);
918// hSpectra[kITS][kKaon][kNeg]->Draw("");
919 // hSpectra[kTOF][kProton][kPos]->Draw("same");
920
921// for(Int_t icharge = 0; icharge < kNCharge; icharge++){
922// TCanvas * c1 = new TCanvas(TString("cAllRatio_")+chargeFlag[icharge],TString("cAllRatio_")+chargeFlag[icharge]);
923// c1->SetGridy();
924// TH2F * hempty = new TH2F(TString("hempty")+long(icharge),"hempty",100,0.2,1, 100, 0.0,2.0);
925// hempty->SetXTitle("p_{t} (GeV/c)");
926// hempty->SetYTitle("ITSsa / TPC");
927// hempty->Draw();
928
929
930// for(Int_t ipart = 0; ipart < kNPart; ipart++) {
931// hSpectra[kITS][ipart][icharge]->Divide(hSpectra[kTPC][ipart][icharge]);
932// hSpectra[kITS][ipart][icharge]->Draw("same");
933// }
934// if(doPrint) c1->Print(TString(c1->GetName())+".png");
935// }
936
937// end of tmp
938
939}
940
941void DrawWithJacek() {
942
943 //1. Convert spectra to dNdeta and sum
944 TH1F * hsum = (TH1F*) htemplate->Clone();
945 hsum->Reset();
ee08b77d 946 Int_t idet= iCombInStudy;
4c0d7bc7 947 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
948 for(Int_t ipart = 0; ipart < kNPart; ipart++){
949 TH1 * h = hSpectra[idet][ipart][icharge];
950 Int_t nbin = h->GetNbinsX();
951 for(Int_t ibin = 1; ibin <= nbin; ibin++){
952 Double_t pt = h->GetBinCenter(ibin);
953 Double_t mt = TMath::Sqrt(pt*pt + mass[ipart]*mass[ipart]);
954 Double_t jacobian = pt/mt;
955 h->SetBinContent(ibin,h->GetBinContent(ibin)*jacobian);
956 h->SetBinError (ibin,h->GetBinError (ibin)*jacobian);
957 Int_t ibinSum = hsum->FindBin(pt);
958 Double_t epsilon = 0.0001;
959 if ( h->GetBinContent(ibin) > 0 &&
960 (TMath::Abs(h->GetBinLowEdge(ibin) - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
961 TMath::Abs(h->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
962 ) {
963 cout << "DISCREPANCY IN BIN RANGES" << endl;
964 cout << pt << " " << ibinSum << " " << ibin << "; " << h->GetBinContent(ibin) << endl
965 << h->GetBinLowEdge(ibin) << "-" << h->GetBinLowEdge(ibin+1) << endl
966 << hsum->GetBinLowEdge(ibinSum) << "-" << hsum->GetBinLowEdge(ibinSum+1) << endl;
967 cout << "" << endl;
968 }
969
970
971 hsum->SetBinContent(ibinSum,hsum->GetBinContent(ibinSum)+h->GetBinContent(ibin)); // EROOR FIXME
972 hsum->SetBinError (ibinSum,0);
973 }
974 // hsum->Add(h);
975 }
976 }
977
978
979 // Load Jacek and Draw both:
980 new TFile ("./Files/dNdPt_Data_Points_ALICE_900GeV.root");
981 TGraphErrors * gJacek = (TGraphErrors*) gDirectory->Get("inel");
982 gJacek->Draw("AP");
983 hsum->Draw("same");
984
985 TGraphErrors * gRatio = AliBWTools::DivideGraphByHisto(gJacek,hsum);
986
987
988 new TCanvas();
989 gRatio->Draw("AP");
990
991
992
993}
994
995
996void DrawRatioToStar() {
997
998 // Star data
999 // gROOT->LoadMacro("StarPPSpectra.C");
1000 TGraphErrors ** gStar = StarPPSpectra();
1001
1002 // ALICE, INEL -> NSD
1003 Double_t scaleYield = 3.58/3.02; // from paper 2
1004 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1005 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
ee08b77d 1006 hSpectra[iCombInStudy][ipart][icharge]->Scale(scaleYield);
4c0d7bc7 1007 }
1008 }
1009
1010
1011 TCanvas * c1 = new TCanvas("cRatioToStarNeg","cRatioToStarNeg");
1012 TH2F * hempty = new TH2F(TString("hemptyNeg"),"hemptyNeg",100,0.,1.5, 100, 0.001,1.8);
1013 hempty->SetXTitle("p_{t} (GeV/c)");
1014 hempty->SetYTitle("ALICE/STAR (NSD)");
1015
1016 TLegend *leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Negative","brNDC");
1017 leg->SetBorderSize(0);
1018 leg->SetLineColor(1);
1019 leg->SetLineStyle(1);
1020 leg->SetLineWidth(1);
1021 leg->SetFillColor(0);
1022 leg->SetFillStyle(1001);
1023
1024
1025 hempty->Draw();
1026 TGraphErrors * g ;
ee08b77d 1027 g = AliBWTools::DivideGraphByHisto(gStar[0],hSpectra[iCombInStudy][kPion][kNeg],1);
4c0d7bc7 1028 g->SetMarkerStyle(kFullCircle);
1029 g->SetMarkerColor(kBlack);
1030 g->Draw("p");
1031 leg->AddEntry(g,"#pi^{-}","lp");
ee08b77d 1032 g = AliBWTools::DivideGraphByHisto(gStar[2],hSpectra[iCombInStudy][kKaon][kNeg],1);
4c0d7bc7 1033 g->SetMarkerStyle(kOpenCircle);
1034 g->SetMarkerColor(kRed);
1035 g->Draw("p");
1036 leg->AddEntry(g,"K^{-}","lp");
ee08b77d 1037 g = AliBWTools::DivideGraphByHisto(gStar[4],hSpectra[iCombInStudy][kProton][kNeg],1);
4c0d7bc7 1038 g->SetMarkerStyle(kOpenSquare);
1039 g->SetMarkerColor(kBlue);
1040 g->Draw("p");
1041 leg->AddEntry(g,"#bar{p}","lp");
1042 leg->Draw();
1043
1044
1045
1046
1047 TCanvas * c2 = new TCanvas("cRatioToStarPos","cRatioToStarPos");
1048 hempty->Draw();
1049 leg = new TLegend(0.6510067,0.1853147,0.8892617,0.4178322,"Positive","brNDC");
1050 leg->SetBorderSize(0);
1051 leg->SetLineColor(1);
1052 leg->SetLineStyle(1);
1053 leg->SetLineWidth(1);
1054 leg->SetFillColor(0);
1055 leg->SetFillStyle(1001);
1056 // TGraphErrors * g ;
ee08b77d 1057 g = AliBWTools::DivideGraphByHisto(gStar[1],hSpectra[iCombInStudy][kPion][kPos],1);
4c0d7bc7 1058 g->SetMarkerStyle(kFullCircle);
1059 g->SetMarkerColor(kBlack);
1060 g->Draw("p");
1061 leg->AddEntry(g,"#pi^{+}","lp");
ee08b77d 1062 g = AliBWTools::DivideGraphByHisto(gStar[3],hSpectra[iCombInStudy][kKaon][kPos],1);
4c0d7bc7 1063 g->SetMarkerStyle(kOpenCircle);
1064 g->SetMarkerColor(kRed);
1065 g->Draw("p");
1066 leg->AddEntry(g,"K^{+}","lp");
ee08b77d 1067 g = AliBWTools::DivideGraphByHisto(gStar[5],hSpectra[iCombInStudy][kProton][kPos],1);
4c0d7bc7 1068 g->SetMarkerStyle(kOpenSquare);
1069 g->SetMarkerColor(kBlue);
1070 g->Draw("p");
1071 leg->AddEntry(g,"p","lp");
1072
1073
1074 leg->Draw();
1075
1076 c1->Update();
1077 c2->Update();
1078 gSystem->ProcessEvents();
1079 c1->Print(TString(c1->GetName()) + ".eps");
1080 c2->Print(TString(c2->GetName()) + ".eps");
1081 ALICEWorkInProgress(c1,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1082 ALICEWorkInProgress(c2,today.Data(),"#splitline{ALICE Preliminary}{Statistical Error Only}");
1083 c1->Update();
1084 c2->Update();
1085 c1->Print(TString(c1->GetName()) + ".png");
1086 c2->Print(TString(c2->GetName()) + ".png");
1087
1088
1089
1090
1091}
1092
1093
1094
1095void DrawRatios() {
1096
1097 // Draws ratios of combined spectra
1098
1099 // Compute ratios
1100 TH1F * hPosNegRatio[kNPart];
1101
1102 for(Int_t ipart = 0; ipart < kNPart; ipart++){
ee08b77d 1103 hPosNegRatio[ipart] = (TH1F*) hSpectra[iCombInStudy][ipart][kPos]->Clone();
1104 hPosNegRatio[ipart]->Divide(hSpectra[iCombInStudy][ipart][kNeg]);
4c0d7bc7 1105 hPosNegRatio[ipart]->SetYTitle(TString(partLabel[ipart][kPos])+"/"+partLabel[ipart][kNeg]);
1106 hPosNegRatio[ipart]->SetMinimum(0.5);
1107 hPosNegRatio[ipart]->SetMaximum(1.5);
1108 }
1109
ee08b77d 1110 TH1F * hKPiRatio = (TH1F*) hSpectra[iCombInStudy][kKaon][kPos]->Clone();
1111 hKPiRatio->Add(hSpectra[iCombInStudy][kKaon][kNeg]);
1112 TH1F * htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1113 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
4c0d7bc7 1114 hKPiRatio->Divide(htmp);
1115 hKPiRatio->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1116
1117 TH1F * hKPiRatioMC[kNTunes];
1118 for(Int_t itune = 0; itune < kNTunes; itune++){
1119 hKPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kKaon][kPos]->Clone();
1120 hKPiRatioMC[itune]->Add(hSpectraMC[itune][kKaon][kNeg]);
1121 TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1122 htmp->Add(hSpectraMC[itune][kPion][kNeg]);
1123 hKPiRatioMC[itune]->Divide(htmp);
1124 hKPiRatioMC[itune]->SetYTitle("(K^{+}+K^{-})/(#pi^{+}+#pi^{-})");
1125 }
1126
1127
1128
ee08b77d 1129 TH1F * hPPiRatio = (TH1F*) hSpectra[iCombInStudy][kProton][kPos]->Clone();
1130 hPPiRatio->Add(hSpectra[iCombInStudy][kProton][kNeg]);
1131 htmp = (TH1F*) hSpectra[iCombInStudy][kPion][kPos]->Clone();
1132 htmp->Add(hSpectra[iCombInStudy][kPion][kNeg]);
4c0d7bc7 1133 hPPiRatio->Divide(htmp);
1134 hPPiRatio->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1135
1136 TH1F * hPPiRatioMC[kNTunes];
1137 for(Int_t itune = 0; itune < kNTunes; itune++){
1138 hPPiRatioMC[itune] = (TH1F*) hSpectraMC[itune][kProton][kPos]->Clone();
1139 hPPiRatioMC[itune]->Add(hSpectraMC[itune][kProton][kNeg]);
1140 TH1F * htmp = (TH1F*) hSpectraMC[itune][kPion][kPos]->Clone();
1141 htmp->Add(hSpectraMC[itune][kPion][kNeg]);
1142 hPPiRatioMC[itune]->Divide(htmp);
1143 hPPiRatioMC[itune]->SetYTitle("(p+#bar{p})/(#pi^{+}+#pi^{-})");
1144 }
1145
1146 // Draw
1147// TH2F * hempty = new TH2F(TString("hempty"),"hempty",100,0.,1.5, 100, 0.001,1.8);
1148// hempty->SetXTitle("p_{t} (GeV/c)");
1149 // hempty->SetYTitle("");
1150
1151 // tmp: overlay levi fits
1152 AliBWFunc * fm2 = new AliBWFunc;
1153 fm2->SetVarType(AliBWFunc::kdNdpt);
1154 TF1 * fLevi[kNPart] [kNCharge];
1155 fLevi[kPion] [kPos] = fm2->GetLevi (mass[0], 0.1243, 7.614785, 1.524167, "fLeviPiPlus");
1156 fLevi[kKaon] [kPos] = fm2->GetLevi (mass[1], 0.1625, 5.869318, 0.186361, "fLeviKPlus");
1157 fLevi[kProton][kPos] = fm2->GetLevi (mass[2], 0.1773, 6.918065, 0.086389, "fLeviPPlus");
1158 fLevi[kPion] [kNeg] = fm2->GetLevi (mass[0], 0.1267, 7.979582, 1.515908, "fLeviPiNeg");
1159 fLevi[kKaon] [kNeg] = fm2->GetLevi (mass[1], 0.1721, 6.927956, 0.191140, "fLeviKNeg");
1160 fLevi[kProton][kNeg] = fm2->GetLevi (mass[2], 0.1782, 8.160362, 0.082091, "fLeviPNeg");
1161 for(Int_t ipart = 0; ipart < kNPart; ipart++){
1162 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1163 fLevi[ipart][icharge]->SetRange(0,4);
1164 }
1165 }
1166
1167
1168 for(Int_t ipart = 0; ipart < kNPart; ipart++){
b7018bc7 1169 TString detName = detFlag[iCombInStudy];
1170 detName.ReplaceAll(" ", "_");
1171 detName.ReplaceAll("+", "");
1172
1173 TCanvas * c1 = new TCanvas(TString("cRatio_")+detName+partFlag[ipart], TString("cRatio_")+detName+partFlag[ipart]);
4c0d7bc7 1174 hPosNegRatio[ipart]->Draw();
1175 TF1 * fRatio = new TF1 (TString("fRatio")+partFlag[ipart], TString(fLevi[ipart][kPos]->GetName())+"/"+fLevi[ipart][kNeg]->GetName());
1176 // fRatio->Draw("same");
1177 fRatio->SetRange(0,5);
1178 if (doPrint) {
1179 c1->Update();
1180 gSystem->ProcessEvents();
b7018bc7 1181 c1->Print(TString(c1->GetName()) + ".png");
4c0d7bc7 1182 }
1183
1184 }
1185
1186
1187 TCanvas * c2 = new TCanvas(TString("cRatio_KPi"),TString("cRatio_KPi"));
1188 hKPiRatio->Draw();
1189 TLegend * lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1190 lMC->SetFillColor(kWhite);
1191
1192 // gROOT->LoadMacro("GetE735Ratios.C");
1193 GetE735Ratios(0,0)->Draw("EX0,same");
1194 GetE735Ratios(0,1)->Draw("EX0,same");
1195 GetE735Ratios(0,2)->Draw("EX0,same");
1196 GetE735Ratios(0,3)->Draw("EX0,same");
1197 hKPiRatio->SetMarkerStyle(20);
1198 hKPiRatio->Draw("same");
1199
1200 for(Int_t itune = 0; itune < kNTunes; itune++){
1201 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1202 hKPiRatioMC[itune]->SetLineWidth(2);
1203 hKPiRatioMC[itune]->Draw("same,chist");
1204 }
1205
1206 lMC->Draw();
1207
1208
1209 TLegend * l = new TLegend( 0.1879, 0.68, 0.54,0.92);
1210 l->SetFillColor(kWhite);
1211 l->AddEntry(hKPiRatio, "ALICE, #sqrt{s} = 900 GeV","lpf");
1212 l->AddEntry(GetE735Ratios(0,0), "E735, #sqrt{s} = 300 GeV","lpf");
1213 l->AddEntry(GetE735Ratios(0,1), "E735, #sqrt{s} = 540 GeV","lpf");
1214 l->AddEntry(GetE735Ratios(0,2), "E735, #sqrt{s} = 1000 GeV","lpf");
1215 l->AddEntry(GetE735Ratios(0,3), "E735, #sqrt{s} = 1800 GeV","lpf");
1216 l->Draw();
1217
1218
1219 TCanvas * c3 = new TCanvas(TString("cRatio_PPi"),TString("cRatio_PPi"));
1220 hPPiRatio->Draw();
1221 hPPiRatio->SetMaximum(0.39);
1222 lMC = new TLegend(0.526846, 0.18007, 0.887584,0.407343);
1223 lMC->SetFillColor(kWhite);
1224 for(Int_t itune = 0; itune < kNTunes; itune++){
1225 lMC->AddEntry(hKPiRatioMC[itune],mcTuneName[itune]);
1226 hKPiRatioMC[itune]->SetLineWidth(2);
1227 hKPiRatioMC[itune]->Draw("same,chist");
1228 }
1229
1230 lMC->Draw();
1231
1232 if (doPrint) {
1233 c2->Update();
1234 gSystem->ProcessEvents();
1235 c2->Print(TString(c2->GetName()) + ".eps");
1236 c3->Update();
1237 gSystem->ProcessEvents();
1238 c3->Print(TString(c3->GetName()) + ".eps");
1239 }
1240
1241
1242}
1243