Updates to Trains. create a job-script to help
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / IdentifiedHighPt / lib / AliHighPtDeDxBase.cxx
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
11 using namespace std;
12
13 ClassImp(AliHighPtDeDxBase);
14
15 //
16 // AliHighPtDeDxBase class
17 //
18 // This class contains the AliHighPtDeDxBase information 
19 //
20
21 //_________________________________________________________
22 AliHighPtDeDxBase::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 //_________________________________________________________
65 AliHighPtDeDxBase::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 //_________________________________________________________
108 AliHighPtDeDxBase::~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 //_________________________________________________________
121 void 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 //_________________________________________________________
182 Bool_t AliHighPtDeDxBase::EventAccepted()
183 {
184   if(fUseRunCut && fRun!=fEventRun)
185     return kFALSE;
186
187   return kTRUE;
188 }
189
190 //_________________________________________________________
191 Bool_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 //_________________________________________________________
228 void 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 //_________________________________________________________
239 void 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 //_________________________________________________________
248 void 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 //_________________________________________________________
275 void AliHighPtDeDxBase::SetUseEtaCut(Bool_t  value)
276 {
277   fUseEtaCut = value;
278
279   if(value == kTRUE)
280     fUseEtaCutAbs = kFALSE;
281 }
282
283 //_________________________________________________________
284 void AliHighPtDeDxBase::SetUseEtaCutAbs(Bool_t  value)
285 {
286   fUseEtaCutAbs = value;
287
288   if(value == kTRUE)
289     fUseEtaCut = kFALSE;
290 }
291
292 //_________________________________________________________
293 TF1* 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 //_________________________________________________________
301 TF1* 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 //___________________________________________________________________________
310 TCanvas* 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 //_______________________________________________________________________
398 void 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 //_______________________________________________________________________
419 void 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 //______________________________________________________________________
445 TCanvas* 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 //_______________________________________________________________________
459 TCanvas* 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 }