Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxBase.cxx
CommitLineData
4ebdd20e 1#include "AliHighPtDeDxBase.h"
2
3#include <TMath.h>
4#include <TStyle.h>
5#include <TROOT.h>
6
7#ifndef __IOSTREAM__
8#include <iostream>
9#endif
10
11using namespace std;
12
13ClassImp(AliHighPtDeDxBase);
14
15//
16// AliHighPtDeDxBase class
17//
18// This class contains the AliHighPtDeDxBase information
19//
20
21//_________________________________________________________
22AliHighPtDeDxBase::AliHighPtDeDxBase():
23 TNamed(),
24 fIsMc(kFALSE),
25 fUseRunCut(kFALSE),
26 fRun(0),
27 fUseEtaCut(kFALSE),
28 fUseEtaCutAbs(kFALSE),
29 fEtaLow(-0.8),
30 fEtaHigh(0.8),
31 fUseFilterCut(kFALSE),
32 fFilter(1),
33 fUsePhiCut(kFALSE),
34 fPhiCutLow(0x0),
35 fPhiCutHigh(0x0),
36 fEventVtxStatus(-999),
37 fEventVtxStatusMc(-999),
38 fEventRun(0),
39 fEventMag(0),
40 fEventTrigger(-1),
41 fTrackCharge(0),
42 fTrackEta(0),
43 fTrackP(0),
44 fTrackPt(0),
45 fTrackFilter(-999),
46 fTrackPhi(-999),
47 fTrackPhiCutVariable(-999),
48 fTrackDeDx(-999),
49 fTrackNcl(-999),
50 fTrackPidMc(0),
51 fTrackPrimaryMc(-1),
52 hVtxStatus(0x0),
53 hNevents(0x0),
54 hPhi(0x0),
55 hEta(0x0),
56 hPt(0x0),
57 hMeanPt(0x0),
58 hNclVsPhiVsPtBefore(0x0),
59 hNclVsPhiVsPtAfter(0x0)
60{
61 // default constructor - do not use
62}
63
64//_________________________________________________________
65AliHighPtDeDxBase::AliHighPtDeDxBase(const char* name, const char* title):
66 TNamed(name, title),
67 fIsMc(kFALSE),
68 fUseRunCut(kFALSE),
69 fRun(0),
70 fUseEtaCut(kFALSE),
71 fUseEtaCutAbs(kFALSE),
72 fEtaLow(-0.8),
73 fEtaHigh(0.8),
74 fUseFilterCut(kFALSE),
75 fFilter(1),
76 fUsePhiCut(kFALSE),
77 fPhiCutLow(0x0),
78 fPhiCutHigh(0x0),
79 fEventVtxStatus(-999),
80 fEventVtxStatusMc(-999),
81 fEventRun(0),
82 fEventMag(0),
83 fEventTrigger(-1),
84 fTrackCharge(0),
85 fTrackEta(0),
86 fTrackP(0),
87 fTrackPt(0),
88 fTrackFilter(-999),
89 fTrackPhi(-999),
90 fTrackPhiCutVariable(-999),
91 fTrackDeDx(-999),
92 fTrackNcl(-999),
93 fTrackPidMc(0),
94 fTrackPrimaryMc(-1),
95 hVtxStatus(0x0),
96 hNevents(0x0),
97 hPhi(0x0),
98 hEta(0x0),
99 hPt(0x0),
100 hMeanPt(0x0),
101 hNclVsPhiVsPtBefore(0x0),
102 hNclVsPhiVsPtAfter(0x0)
103{
104 // named constructor
105}
106
107//_________________________________________________________
108AliHighPtDeDxBase::~AliHighPtDeDxBase()
109{
110 delete hVtxStatus;
111 delete hNevents;
112 delete hEta;
113 delete hPhi;
114 delete hPt;
115 delete hMeanPt;
116 delete hNclVsPhiVsPtBefore;
117 delete hNclVsPhiVsPtAfter;
118}
119
120//_________________________________________________________
121void AliHighPtDeDxBase::Init(Int_t nPtBins, Double_t* ptBins)
122{
123 //
124 // Create histograms
125 //
126 hVtxStatus = new TH1D("hVtxStatus", "Vtx status - No Vtx = -1, Vtx outside cut = 0, Vtx inside = 1",
127 3, -1.5, 1.5);
128 hVtxStatus->Sumw2();
129 hVtxStatus->SetDirectory(0);
130
131 hNevents = new TH1D("hNevents", "N events - No Vtx = 0, Vtx OK = 1",
132 2, 0, 2);
133 hNevents->Sumw2();
134 hNevents->SetDirectory(0);
135
136 hEta = new TH1D("hEta", "#eta distribution; #eta; Counts",
137 100, -1.0, 1.0);
138 hEta->Sumw2();
139 hEta->SetDirectory(0);
140
141 hPhi = new TH1D("hPhi", "#varphi distribution; #varphi; Counts",
142 360, 0, TMath::TwoPi());
143 hPhi->Sumw2();
144 hPhi->SetDirectory(0);
145
146 hPt = new TH1D("hPt", "p_{T} spectrum; p_{T} [GeV/c]; Counts",
147 nPtBins, ptBins);
148 hPt->Sumw2();
149 hPt->SetDirectory(0);
150
151 hMeanPt = new TProfile("hMeanPt", "mean p_{T}; p_{T} [GeV/c]; mean p_{T}",
152 nPtBins, ptBins);
153 hMeanPt->SetDirectory(0);
154
155 const Int_t nPhiBins = 50;
156 const Double_t phiBinSize = TMath::Pi()/9.0/nPhiBins;
157 Double_t phiBins[nPhiBins+1];
158
159 for(Int_t i = 0; i <= nPhiBins; i++) {
160
161 phiBins[i] = phiBinSize*i;
162 }
163
164 const Int_t nNclBins = 45;
165 const Double_t nclBinSize = 90.0/nNclBins;
166 Double_t nclBins[nNclBins+1];
167
168 for(Int_t i = 0; i <= nNclBins; i++) {
169
170 nclBins[i] = 69.5 + nclBinSize*i;
171 }
172
173 hNclVsPhiVsPtBefore = new TH3F("hNclVsPhiVsPtBefore", "<Ncl> vs Pt and #phi (before #phi cut); p_{T} [GeV/c]; fmod(#phi, #pi/9.0)",
174 nPtBins, ptBins, nPhiBins, phiBins, nNclBins, nclBins);
175 hNclVsPhiVsPtBefore->SetDirectory(0);
176 hNclVsPhiVsPtAfter = new TH3F("hNclVsPhiVsPtAfter", "<Ncl> vs Pt and #phi (after #phi cut); p_{T} [GeV/c]; fmod(#phi, #pi/9.0)",
177 nPtBins, ptBins, nPhiBins, phiBins, nNclBins, nclBins);
178 hNclVsPhiVsPtAfter->SetDirectory(0);
179}
180
181//_________________________________________________________
182Bool_t AliHighPtDeDxBase::EventAccepted()
183{
184 if(fUseRunCut && fRun!=fEventRun)
185 return kFALSE;
186
187 return kTRUE;
188}
189
190//_________________________________________________________
191Bool_t AliHighPtDeDxBase::TrackAccepted()
192{
193 if(fUseFilterCut && !(fTrackFilter&fFilter))
194 return kFALSE;
195
196 if(fUseEtaCut && (fTrackEta<fEtaLow || fTrackEta>fEtaHigh))
197 return kFALSE;
198
199 if(fUseEtaCutAbs && (TMath::Abs(fTrackEta)<fEtaLow || TMath::Abs(fTrackEta)>fEtaHigh))
200 return kFALSE;
201
202 if(fUsePhiCut) {
203
204 fTrackPhiCutVariable = fTrackPhi;
205 if(fEventMag<0) // for negatve polarity field
206 fTrackPhiCutVariable = TMath::TwoPi() - fTrackPhiCutVariable;
207 if(fTrackCharge<0) // for negatve charge
208 fTrackPhiCutVariable = TMath::TwoPi()-fTrackPhiCutVariable;
209 if(fTrackPhiCutVariable < 0)
210 cout << "Warning!!!!! phi < 0: " << fTrackPhiCutVariable << endl;
211
212 fTrackPhiCutVariable += TMath::Pi()/18.0; // to center gap in the middle
213 fTrackPhiCutVariable = fmod(fTrackPhiCutVariable, TMath::Pi()/9.0);
214
215 hNclVsPhiVsPtBefore->Fill(fTrackPt, fTrackPhiCutVariable, fTrackNcl);
216
217 if(fTrackPt>2.0 && fTrackPhiCutVariable<fPhiCutHigh->Eval(fTrackPt)
218 && fTrackPhiCutVariable>fPhiCutLow->Eval(fTrackPt))
219 return kFALSE; // reject track
220
221 hNclVsPhiVsPtAfter->Fill(fTrackPt, fTrackPhiCutVariable, fTrackNcl);
222 }
223
224 return kTRUE;
225}
226
227//_________________________________________________________
228void AliHighPtDeDxBase::FillEventInfo()
229{
230 if(fEventTrigger==1) {
231 hVtxStatus->Fill(fEventVtxStatus);
232 if(fEventVtxStatus != 0) {
233 hNevents->Fill(1.0+0.5*fEventVtxStatus);
234 }
235 }
236}
237
238//_________________________________________________________
239void AliHighPtDeDxBase::FillTrackInfo(Float_t weight)
240{
241 hEta->Fill(fTrackEta, weight);
242 hPhi->Fill(fTrackPhi, weight);
243 hPt->Fill(fTrackPt, weight);
244 hMeanPt->Fill(fTrackPt, fTrackPt);
245}
246
247//_________________________________________________________
248void AliHighPtDeDxBase::Print(Option_t* option) const
249{
250 cout << ClassName() << " : " << GetName() << endl
251 << "Event cuts: " << endl;
252 if(fUseRunCut)
253 cout << " Run cut: " << fRun << endl;
254 else
255 cout << " Run cut is diabled " << endl;
256
257 cout << "Track cuts: " << endl;
258 if(fUseFilterCut)
259 cout << " Filter cut: " << fFilter << endl;
260 else
261 cout << " Filter cut is diabled " << endl;
262 if(fUseEtaCut)
263 cout << " Eta range: " << fEtaLow << " - " << fEtaHigh << endl;
264 if(fUseEtaCutAbs)
265 cout << " |Eta| range: " << fEtaLow << " - " << fEtaHigh << endl;
266 else
267 cout << " Eta cut is diabled " << endl;
268 if(fUsePhiCut)
269 cout << " Phi cut is ENABLED" << endl;
270 else
271 cout << " Phi cut is diabled " << endl;
272}
273
274//_________________________________________________________
275void AliHighPtDeDxBase::SetUseEtaCut(Bool_t value)
276{
277 fUseEtaCut = value;
278
279 if(value == kTRUE)
280 fUseEtaCutAbs = kFALSE;
281}
282
283//_________________________________________________________
284void AliHighPtDeDxBase::SetUseEtaCutAbs(Bool_t value)
285{
286 fUseEtaCutAbs = value;
287
288 if(value == kTRUE)
289 fUseEtaCut = kFALSE;
290}
291
292//_________________________________________________________
293TF1* AliHighPtDeDxBase::GetStandardPhiCutLow()
294{
295 // TF1* cutLow = new TF1("StandardPhiCutLow", "-0.01/x+pi/18.0-0.015", 0, 50);
296 TF1* cutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
297 return cutLow;
298}
299
300//_________________________________________________________
301TF1* AliHighPtDeDxBase::GetStandardPhiCutHigh()
302{
303 // TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.55/x/x+pi/18.0+0.03", 0, 50);
304 TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
305 return cutHigh;
306}
307
308
309//___________________________________________________________________________
310TCanvas* AliHighPtDeDxBase::DrawPhiCutHistograms()
311{
312 gStyle->SetOptStat(0);
313
314 TCanvas* cPhiCut = FindCanvas("cPhiCut", 1200, 800);
315 cPhiCut->SetTitle("phi cut histograms");
316 cPhiCut->Clear();
317 cPhiCut->Divide(3,2);
318
319 // cPhiCut->cd(1);
320 // TH2D* hPhiVsPtBefore = (TH2D*)hNclVsPhiVsPtBefore->Project3D("yx");
321 // hPhiVsPtBefore->SetName("hPhiVsPtBefore");
322 // hPhiVsPtBefore->SetTitle("Phi vs p_{T} (before cuts)");
323 // MakeNice2dHisto(hPhiVsPtBefore, gPad, kTRUE);
324 // hPhiVsPtBefore->Draw("COL");
325 // fPhiCutHigh->SetRange(2.0, 50.0);
326 // fPhiCutHigh->Draw("SAME");
327 // fPhiCutLow->SetRange(2.0, 50.0);
328 // fPhiCutLow->Draw("SAME");
329
330 cPhiCut->cd(1);
331 TProfile2D* hNclBefore = hNclVsPhiVsPtBefore->Project3DProfile("yx");
332 hNclBefore->SetName("hNclBefore");
333 hNclBefore->SetTitle("<Ncl> (before cuts); p_{T} [GeV/c]; #varphi'");
334 MakeNice2dHisto(hNclBefore, gPad, kTRUE);
335 DrawNice(hNclBefore, 0, 0, 0, "COLZ");
336 fPhiCutHigh->Draw("SAME");
337 fPhiCutLow->Draw("SAME");
338
339 cPhiCut->cd(2);
340 gPad->SetGridy();
341 TH2D* hNclVsPtBefore = (TH2D*)hNclVsPhiVsPtBefore->Project3D("zx");
342 hNclVsPtBefore->SetName("hNclVsPtBefore");
343 hNclVsPtBefore->SetTitle("; p_{T} [GeV/c]; Ncl;");
344 MakeNice2dHisto(hNclVsPtBefore, gPad, kTRUE);
345 hNclVsPtBefore->Draw("COL");
346 TProfile* hNclVsPtBeforeProf = hNclVsPtBefore->ProfileX();
347 hNclVsPtBeforeProf->SetMarkerStyle(29);
348 hNclVsPtBeforeProf->SetMarkerColor(2);
349 hNclVsPtBeforeProf->Draw("SAME");
350
351 // cPhiCut->cd(2);
352 // TH2D* hPhiVsPtAfter = (TH2D*)hNclVsPhiVsPtAfter->Project3D("yx");
353 // MakeNice2dHisto(hNclVsPtAfter, gPad, kTRUE);
354 // hPhiVsPtAfter->SetName("hPhiVsPtAfter");
355 // hPhiVsPtAfter->SetTitle("Phi vs p_{T} (after cuts)");
356 // DrawNice(hPhiVsPtAfter, 0, 0, 0, "COL");
357
358 cPhiCut->cd(4);
359 TProfile2D* hNclAfter = (TProfile2D*)hNclVsPhiVsPtAfter->Project3DProfile("yx");
360 hNclAfter->SetName("hNclAfter");
361 hNclAfter->SetTitle("<Ncl> (after cuts); p_{T} [GeV/c]; #varphi'");
362 MakeNice2dHisto(hNclAfter, gPad, kTRUE);
363 DrawNice(hNclAfter, 0, 0, 0, "COLZ");
364
365 cPhiCut->cd(5);
366 gPad->SetGridy();
367 TH2D* hNclVsPtAfter = (TH2D*)hNclVsPhiVsPtAfter->Project3D("zx");
368 hNclVsPtAfter->SetName("hNclVsPtAfter");
369 hNclVsPtAfter->SetTitle("; p_{T} [GeV/c]; Ncl;");
370 MakeNice2dHisto(hNclVsPtAfter, gPad, kTRUE);
371 hNclVsPtAfter->Draw("COL");
372 TProfile* hNclVsPtAfterProf = hNclVsPtAfter->ProfileX();
373 hNclVsPtAfterProf->SetMarkerStyle(29);
374 hNclVsPtAfterProf->SetMarkerColor(2);
375 hNclVsPtAfterProf->Draw("SAME");
376
377 cPhiCut->cd(3);
378 gPad->SetGridy();
379 TH1D* hEfficiency = (TH1D*)hNclVsPhiVsPtAfter->Project3D("x");
380 hEfficiency->SetName("hEfficiency");
381 hEfficiency->SetTitle("; p_{T} [GeV/c]; Cut efficiency;");
382 TH1D* hHelp = (TH1D*)hNclVsPhiVsPtBefore->Project3D("x");
383 hEfficiency->Divide(hEfficiency, hHelp, 1, 1, "B");
384 delete hHelp;
385 MakeNice1dHisto(hEfficiency, gPad);
386 hEfficiency->SetMarkerStyle(20);
387 hEfficiency->Draw("P");
388
389 return cPhiCut;
390}
391
392/*
393 These methods does not fit in very well here, but for now I put them
394 here. Maybe they should go in some tool class.
395 */
396
397//_______________________________________________________________________
398void AliHighPtDeDxBase::MakeNice1dHisto(TH1* hist, TVirtualPad* c1)
399{
400 if(c1) {
401
402 c1->SetLeftMargin(0.12);
403 c1->SetRightMargin(0.02);
404 c1->SetBottomMargin(0.12);
405 c1->SetTopMargin(0.06);
406 }
407
408 gStyle->SetTitleH(0.08);
409 gStyle->SetTitleW(0.46);
410 gStyle->SetTitleX(0.29);
411 gStyle->SetTitleY(0.99);
412 hist->GetXaxis()->SetTitleSize(0.07);
413 hist->GetYaxis()->SetTitleSize(0.07);
414 hist->GetXaxis()->SetTitleOffset(0.7);
415 hist->GetYaxis()->SetTitleOffset(0.8);
416}
417
418//_______________________________________________________________________
419void AliHighPtDeDxBase::MakeNice2dHisto(TH2* hist, TVirtualPad* c1, Bool_t colz)
420{
421 gStyle->SetTitleH(0.08);
422 gStyle->SetTitleW(0.76);
423 gStyle->SetTitleX(0.14);
424 gStyle->SetTitleY(0.98);
425 if(colz)
426 gStyle->SetTitleW(0.68);
427
428 hist->GetXaxis()->SetTitleSize(0.07);
429 hist->GetYaxis()->SetTitleSize(0.07);
430 hist->GetXaxis()->SetTitleOffset(0.7);
431 hist->GetYaxis()->SetTitleOffset(0.8);
432 if(c1) {
433
434 c1->SetLeftMargin(0.12);
435 if(!colz)
436 c1->SetRightMargin(0.02);
437 else
438 c1->SetRightMargin(0.12);
439 c1->SetBottomMargin(0.12);
440 c1->SetTopMargin(0.06);
441 }
442}
443
444//______________________________________________________________________
445TCanvas* AliHighPtDeDxBase::FindCanvas(const Char_t* canvasName,
446 Int_t xwidth, Int_t ywidth)
447{
448
449 TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(canvasName);
450 if(!c1)
451 c1 = new TCanvas(canvasName, canvasName, xwidth, ywidth);
452 else
453 c1->SetWindowSize(xwidth, ywidth);
454
455 return c1;
456}
457
458//_______________________________________________________________________
459TCanvas* AliHighPtDeDxBase::DrawNice(TH1* hist, const Char_t* canvasName,
460 Int_t xwidth, Int_t ywidth, const Char_t* option)
461{
462 TCanvas* canv = 0;
463
464 if(canvasName) {
465 canv = FindCanvas(canvasName, xwidth, ywidth);
466 canv->Clear();
467 }
468
469 Bool_t is2dHist = hist->InheritsFrom("TH2");
470
471 if(is2dHist) {
472
473 TString opt(option);
474 opt.ToUpper();
475
476 if (opt.Contains("COL"))
477 MakeNice2dHisto((TH2*)hist, canv, kTRUE);
478 else
479 MakeNice2dHisto((TH2*)hist, canv);
480
481 } else
482 MakeNice1dHisto(hist, canv);
483
484 hist->Draw(option);
485
486 return canv;
487}