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