1 /*************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
16 // -----------------------------------------------------------------------
17 // Track Cut class for the DiHadronPID analysis.
18 // -----------------------------------------------------------------------
19 // Author: Misha Veldhoen (misha.veldhoen@cern.ch)
21 #include "AliAODTrackCutsDiHadronPID.h"
32 #include "TIterator.h"
35 #include "AliAODTrack.h"
36 #include "AliAODEvent.h"
37 #include "AliAODVertex.h"
38 #include "AliAODMCParticle.h"
42 #include "AliPIDResponse.h"
43 #include "AliTPCPIDResponse.h"
45 #include "AliTrackDiHadronPID.h"
46 #include "AliHistToolsDiHadronPID.h"
47 #include "AliFunctionsDiHadronPID.h"
49 ClassImp(AliAODTrackCutsDiHadronPID);
51 // -----------------------------------------------------------------------
52 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID():
59 fMinimumNumberOfTPCClusters(-999),
61 fMinSPDHitsForPtDeptDCAcut(0),
62 fPtDeptDCAxyCutFormula(0x0),
65 fLowPtNSigmaTOFOnly(kFALSE),
66 fUseNSigmaOnPIDAxes(kFALSE),
68 fTestFilterMask(kFALSE),
70 fTestMaxRapidity(kFALSE),
72 fTestNumberOfTPCClusters(kFALSE),
74 fTestTOFmismatch(kFALSE),
75 fTestPtDeptDCAcut(kFALSE),
76 fDataTrackQAHistos(0x0),
77 fHistAcceptedFilterBits(0x0),
78 fRelevantBitsArray(0x0),
79 fTOFMatchingStat(0x0),
80 fPrimRecMCTrackQAHistos(0x0),
81 fPrimGenMCTrackQAHistos(0x0),
82 fSecRecMCTrackQAHistos(0x0),
83 fSecGenMCTrackQAHistos(0x0),
92 // Default constructor
95 cout<<"AliAODTrackCutsDiHadronPID Default Constructor Called."<<endl;
96 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
98 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
100 // Initialize the data histograms.
101 if (iHistoName < 3) {
102 fHistDataPt[iHistoName] = 0x0;
103 fHistDataPhiEtaPt[iHistoName] = 0x0;
104 fHistDataNTracks[iHistoName] = 0x0;
107 fHistDataDCAxy[iHistoName] = 0x0;
108 fHistDataDCAz[iHistoName] = 0x0;
110 // Initialize PID histograms.
111 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
112 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
113 fHistPrimRecPID[iHistoName][iSpecies][iPtClass] = 0x0;
114 fHistPrimRecMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
115 fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
116 fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
117 fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
122 fNTracks[iHistoName] = 0;
124 // Initialize data DCA for one sigma.
125 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
127 // Initialize MC histograms.
128 fHistPrimGenMCPt[iHistoName] = 0x0;
129 fHistPrimRecPtGenPt[iHistoName] = 0x0;
130 fHistPrimGenMCPhiEtaPt[iHistoName] = 0x0;
131 fHistPrimRecMCPt[iHistoName] = 0x0;
132 fHistPrimRecMCPhiEtaPt[iHistoName] = 0x0;
133 fHistPrimRecNTracks[iHistoName] = 0x0;
135 fHistSecGenMCPt[iHistoName] = 0x0;
136 fHistSecGenMCPhiEtaPt[iHistoName] = 0x0;
137 fHistSecRecMCPt[iHistoName] = 0x0;
138 fHistSecRecMCPhiEtaPt[iHistoName] = 0x0;
140 // Initialze MC DCA histograms
141 fHistPrimRecMCDCA[iHistoName] = 0x0;
142 fHistSecRecMCDCAMat[iHistoName] = 0x0;
143 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
145 // Initialize other stuff.
146 fHistRequested[iHistoName] = kFALSE;
147 f3DSpectraEnabeled[iHistoName] = kFALSE;
148 fPIDHistosEnabeled[iHistoName] = kFALSE;
151 // TODO: It would be a bit neater to initialize all values to zero
152 // instead of the default values...
153 InitializeDefaultHistoNamesAndAxes();
157 // -----------------------------------------------------------------------
158 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID(const char* name):
159 TNamed(name,"AOD Track Cuts"),
165 fMinimumNumberOfTPCClusters(-999),
167 fMinSPDHitsForPtDeptDCAcut(0),
168 fPtDeptDCAxyCutFormula(0x0),
171 fLowPtNSigmaTOFOnly(kFALSE),
172 fUseNSigmaOnPIDAxes(kFALSE),
174 fTestFilterMask(kFALSE),
176 fTestMaxRapidity(kFALSE),
178 fTestNumberOfTPCClusters(kFALSE),
180 fTestTOFmismatch(kFALSE),
181 fTestPtDeptDCAcut(kFALSE),
182 fDataTrackQAHistos(0x0),
183 fHistAcceptedFilterBits(0x0),
184 fRelevantBitsArray(0x0),
185 fTOFMatchingStat(0x0),
186 fPrimRecMCTrackQAHistos(0x0),
187 fPrimGenMCTrackQAHistos(0x0),
188 fSecRecMCTrackQAHistos(0x0),
189 fSecGenMCTrackQAHistos(0x0),
201 cout<<"AliAODTrackCutsDiHadronPID Named Constructor Called."<<endl;
202 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
204 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
206 // Initialize the data histograms.
207 if (iHistoName < 3) {
208 fHistDataPt[iHistoName] = 0x0;
209 fHistDataPhiEtaPt[iHistoName] = 0x0;
210 fHistDataNTracks[iHistoName] = 0x0;
213 fHistDataDCAxy[iHistoName] = 0x0;
214 fHistDataDCAz[iHistoName] = 0x0;
216 // Initialize PID histograms.
217 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
218 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
219 fHistPrimRecPID[iHistoName][iSpecies][iPtClass] = 0x0;
220 fHistPrimRecMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
221 fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
222 fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
223 fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
229 fNTracks[iHistoName] = 0;
231 // Initialize data DCA for one sigma.
232 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
234 // Initialize MC histograms.
235 fHistPrimGenMCPt[iHistoName] = 0x0;
236 fHistPrimRecPtGenPt[iHistoName] = 0x0;
237 fHistPrimGenMCPhiEtaPt[iHistoName] = 0x0;
238 fHistPrimRecMCPt[iHistoName] = 0x0;
239 fHistPrimRecMCPhiEtaPt[iHistoName] = 0x0;
240 fHistPrimRecNTracks[iHistoName] = 0x0;
242 fHistSecGenMCPt[iHistoName] = 0x0;
243 fHistSecGenMCPhiEtaPt[iHistoName] = 0x0;
244 fHistSecRecMCPt[iHistoName] = 0x0;
245 fHistSecRecMCPhiEtaPt[iHistoName] = 0x0;
247 // Initialze MC DCA histograms
248 fHistPrimRecMCDCA[iHistoName] = 0x0;
249 fHistSecRecMCDCAMat[iHistoName] = 0x0;
250 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
252 // Initialize other stuff.
253 fHistRequested[iHistoName] = kFALSE;
254 f3DSpectraEnabeled[iHistoName] = kFALSE;
255 fPIDHistosEnabeled[iHistoName] = kFALSE;
258 InitializeDefaultHistoNamesAndAxes();
262 // -----------------------------------------------------------------------
263 void AliAODTrackCutsDiHadronPID::InitializeDefaultHistoNamesAndAxes() {
265 // Initializes the histogram name conventions and the ranges of all the histograms
266 // by their default values. This method should only be called by the
267 // (default) constructor.
269 // TODO: User should be able to change these standard values at initialization.
270 // TODO: User should be able to retrieve all these variables with appropriate Getters.
272 cout<<"AliAODTrackCutsDiHadronPID - Initializing Default Histogram Names and axes..."<<endl;
273 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
275 // Setting the Pt axis for all histograms except the PID and Mismatch histograms -> Now the same for
277 Double_t ptaxis[52] = {0.500,0.525,0.550,0.575,0.600,0.625,0.650,0.675,0.700,
278 0.75,0.80,0.85,0.90,0.95,1.00,
279 1.05,1.10,1.15,1.20,1.25,1.30,1.35,1.40,1.45,1.50,1.55,1.60,1.65,1.70,
280 1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
281 3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
282 for (Int_t iPtBins = 0; iPtBins < 52; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
285 Double_t ptaxis[57] = {0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.00,
286 1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,
287 3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5,9.0,9.5,10.0};
288 for (Int_t iPtBins = 0; iPtBins < 57; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
292 // Setting the pt range of the five PID histogram pt classes.
293 Double_t ptboundarypid[6] = {0.5,0.7,1.0,1.7,3.0,5.0};
294 for (Int_t iPtBoundary = 0; iPtBoundary < 6; iPtBoundary++) {fPtBoundaryPID[iPtBoundary] = ptboundarypid[iPtBoundary];}
296 // Setting the number of bins in pt for the five PID histogram pt classes.
297 Int_t nptbinspid[5] = {8,6,14,13,10};
298 for (Int_t iPtBins = 0; iPtBins < 5; iPtBins++) {fNPtBinsPID[iPtBins] = nptbinspid[iPtBins];}
300 // Setting the TOF axes for the PID histograms.
301 //Double_t tofsigmaapprox = 80.;
302 Double_t toflowerbound[5][3] = {{-2000.,-6000.,-10000.},{-2000.,-4000.,-10000.},{-1000.,-2000.,-5000.},{-1000.,-1000.,-2500.},{-500.,-500.,-1000.}};
303 Double_t tofupperbound[5][3] = {{10000.,10000.,6000.},{10000.,8000.,6000.},{6000.,6000.,6000.},{6000.,6000.,6000.},{4000.,4000.,6000.}};
304 Int_t tofbins[5][3] = {{120,160,160},{120,120,160},{140,140,165},{140,140,170},{90,90,140}};
305 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
306 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
307 fTOFLowerBound[iPtClass][iSpecies] = toflowerbound[iPtClass][iSpecies];
308 //if (fUseNSigmaOnPIDAxes) {toflowerbound[iPtClass][iSpecies] /= tofsigmaapprox;}
309 fTOFUpperBound[iPtClass][iSpecies] = tofupperbound[iPtClass][iSpecies];
310 //if (fUseNSigmaOnPIDAxes) {tofupperbound[iPtClass][iSpecies] /= tofsigmaapprox;}
311 fTOFbins[iPtClass][iSpecies] = tofbins[iPtClass][iSpecies];
315 // Setting the TPC axes for the PID histograms.
316 //Double_t tpcsigmaaxpprox = 3.5;
317 Double_t tpclowerbound[5][3] = {{-20.,-50.,-100.},{-20.,-30.,-80.},{-25.,-25.,-45.},{-25.,-25.,-45.},{-25.,-20.,-20.}};
318 Double_t tpcupperbound[5][3] = {{60.,30.,20.},{60.,40.,20.},{50.,50.,25.},{45.,45.,25.},{25.,30.,30.}}; // Check highest pT bin boundaries for K,p
319 Int_t tpcbins[5][3] = {{80,80,120},{80,70,100},{75,75,70},{70,70,70},{50,50,50}};
320 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
321 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
322 fTPCLowerBound[iPtClass][iSpecies] = tpclowerbound[iPtClass][iSpecies];
323 //if (fUseNSigmaOnPIDAxes) {tpclowerbound[iPtClass][iSpecies] /= tpcsigmaaxpprox;}
324 fTPCUpperBound[iPtClass][iSpecies] = tpcupperbound[iPtClass][iSpecies];
325 //if (fUseNSigmaOnPIDAxes) {tpcupperbound[iPtClass][iSpecies] /= tpcsigmaaxpprox;}
326 fTPCbins[iPtClass][iSpecies] = tpcbins[iPtClass][iSpecies];
330 // Names for the 12 (species,charge) combinations considered.
331 const char* histoname[12] = {"AllCharged","Pos","Neg","AllPion","PosPion","NegPion","AllKaon","PosKaon","NegKaon","AllProton","PosProton","NegProton"};
332 const char* histolatex[12] = {"ch","ch^{+}","ch^{-}","#pi^{+/-}","#pi^{+}","#pi^{-}","K^{+/-}","K^{+}","K^{-}","p/#bar{p}","p","#bar{p}"};
333 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
334 fHistoName[iHistoName] = histoname[iHistoName];
335 fHistoLatex[iHistoName] = histolatex[iHistoName];
338 // Names of the 3 species considered.
339 const char* particlename[3] = {"Pion","Kaon","Proton"};
340 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {fParticleName[iSpecies] = particlename[iSpecies];}
342 // Names of the 5 pt classes considered for the PID histograms.
343 const char* ptclassname[5] = {"VeryLowPt","LowPt","MedPt","MedHighPt","HighPt"};
344 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {fPtClassName[iPtClass] = ptclassname[iPtClass];}
348 // -----------------------------------------------------------------------
349 AliAODTrackCutsDiHadronPID::~AliAODTrackCutsDiHadronPID() {
355 cout<<"AliAODTrackCutsDiHadronPID Destructor Called."<<endl;
356 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
358 if (fDataTrackQAHistos) delete fDataTrackQAHistos;
359 fDataTrackQAHistos = 0x0;
360 if (fPrimRecMCTrackQAHistos) delete fPrimRecMCTrackQAHistos;
361 fPrimRecMCTrackQAHistos = 0x0;
362 if (fPrimGenMCTrackQAHistos) delete fPrimGenMCTrackQAHistos;
363 fPrimGenMCTrackQAHistos = 0x0;
364 if (fSecRecMCTrackQAHistos) delete fSecRecMCTrackQAHistos;
365 fSecRecMCTrackQAHistos = 0x0;
366 if (fSecGenMCTrackQAHistos) delete fSecGenMCTrackQAHistos;
367 fSecGenMCTrackQAHistos = 0x0;
371 // -----------------------------------------------------------------------
372 Long64_t AliAODTrackCutsDiHadronPID::Merge(TCollection* list) {
378 cout<<"AliAODTrackCutsDiHadronPID Merger Called."<<endl;
379 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
381 Bool_t HistosOK = kTRUE;
384 if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos||!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {HistosOK = kFALSE;}
386 if (!fDataTrackQAHistos) {HistosOK = kFALSE;}
390 cout<<"AliAODTrackCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
395 //cout<<"No list available..."<<endl;
398 if (list->IsEmpty()) {
399 //cout<<"List is empty..."<<endl;
403 //Int_t NEntries = list->GetEntries();
404 //cout<<"Supplied list has "<<NEntries<<" entries."<<endl;
406 TIterator* iter = list->MakeIterator();
409 // List of collections
410 TList collection_fDataTrackQAHistos;
411 TList collection_fPrimRecMCTrackQAHistos;
412 TList collection_fPrimGenMCTrackQAHistos;
413 TList collection_fSecRecMCTrackQAHistos;
414 TList collection_fSecGenMCTrackQAHistos;
418 while ((obj = iter->Next())) {
419 AliAODTrackCutsDiHadronPID* entry = dynamic_cast<AliAODTrackCutsDiHadronPID*> (obj);
420 if (entry == 0) continue;
422 // Check if the object to be merged really has the same name! (FIXME!)
424 // Getting the lists from obj.
425 TList* list_fDataTrackQAHistos = entry->GetListOfDataQAHistos();
426 TList* list_fPrimRecMCTrackQAHistos = entry->GetListOfPrimRecMCTrackQAHistos();
427 TList* list_fPrimGenMCTrackQAHistos = entry->GetListOfPrimGenMCTrackQAHistos();
428 TList* list_fSecRecMCTrackQAHistos = entry->GetListOfSecRecMCTrackQAHistos();
429 TList* list_fSecGenMCTrackQAHistos = entry->GetListOfSecGenMCTrackQAHistos();
431 // Adding the retrieved lists to the collection.
432 if (list_fDataTrackQAHistos) collection_fDataTrackQAHistos.Add(list_fDataTrackQAHistos);
433 if (list_fPrimRecMCTrackQAHistos) collection_fPrimRecMCTrackQAHistos.Add(list_fPrimRecMCTrackQAHistos);
434 if (list_fPrimGenMCTrackQAHistos) collection_fPrimGenMCTrackQAHistos.Add(list_fPrimGenMCTrackQAHistos);
435 if (list_fSecRecMCTrackQAHistos) collection_fSecRecMCTrackQAHistos.Add(list_fSecRecMCTrackQAHistos);
436 if (list_fSecGenMCTrackQAHistos) collection_fSecGenMCTrackQAHistos.Add(list_fSecGenMCTrackQAHistos);
438 //cout<<"Entries intermediate list Data: "<<collection_fDataTrackQAHistos.GetEntries()<<endl;
439 //cout<<"Entries intermediate list RecMC: "<<collection_fPrimRecMCTrackQAHistos.GetEntries()<<endl;
444 // Merging. Note that we require the original list to exist.
445 // * Assume that if the collection happens to be empty, then nothing will happen.
446 // |-> This of course leads to some problems, since if the first file does not
447 // | have such a list, then the merged file will have no results...
448 // | IDEA: if fDataTrackQAHistos does not exist, then create a new one. (TO DO)
449 // * All other variables are taken from the original object.
450 if (fDataTrackQAHistos) fDataTrackQAHistos->Merge(&collection_fDataTrackQAHistos);
451 if (fPrimRecMCTrackQAHistos) fPrimRecMCTrackQAHistos->Merge(&collection_fPrimRecMCTrackQAHistos);
452 if (fPrimGenMCTrackQAHistos) fPrimGenMCTrackQAHistos->Merge(&collection_fPrimGenMCTrackQAHistos);
453 if (fSecRecMCTrackQAHistos) fSecRecMCTrackQAHistos->Merge(&collection_fSecRecMCTrackQAHistos);
454 if (fSecGenMCTrackQAHistos) fSecGenMCTrackQAHistos->Merge(&collection_fSecGenMCTrackQAHistos);
462 // -----------------------------------------------------------------------
463 void AliAODTrackCutsDiHadronPID::PrintCuts() const { /* NOT IMPLEMENTED */ }
465 // -----------------------------------------------------------------------
466 TList* AliAODTrackCutsDiHadronPID::GetListOfDataQAHistos() const {
468 // Returns the list of data histograms.
469 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
471 if (fDataTrackQAHistos) {
472 return fDataTrackQAHistos;
479 // -----------------------------------------------------------------------
480 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimRecMCTrackQAHistos() const {
482 // Returns the list of histograms of reconstructed primary tracks.
483 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
485 if (fPrimRecMCTrackQAHistos) {
486 return fPrimRecMCTrackQAHistos;
493 // -----------------------------------------------------------------------
494 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimGenMCTrackQAHistos() const {
496 // Returns the list of histograms of generator level primary particles.
497 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
499 if (fPrimGenMCTrackQAHistos) {
500 return fPrimGenMCTrackQAHistos;
507 // -----------------------------------------------------------------------
508 TList* AliAODTrackCutsDiHadronPID::GetListOfSecRecMCTrackQAHistos() const {
510 // Returns the list of histograms of reconstructed secondary tracks.
511 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
513 if (fSecRecMCTrackQAHistos) {
514 return fSecRecMCTrackQAHistos;
521 // -----------------------------------------------------------------------
522 TList* AliAODTrackCutsDiHadronPID::GetListOfSecGenMCTrackQAHistos() const {
524 // Returns the list of histograms of generator level secondary particles.
525 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
527 if (fSecGenMCTrackQAHistos) {
528 return fSecGenMCTrackQAHistos;
535 // -----------------------------------------------------------------------
536 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
538 // Returns a projection in TOF of the data histogram.
539 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
541 Int_t ptclass = GetPtClass(ptbin);
542 if (ptclass == -1) return 0x0;
543 Int_t bininptclass = GetBinInPtClass(ptbin);
545 // Retrieve original the 3D histogram (we don't own it).
546 TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
548 // Make the projection.
549 TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojection_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
551 // Some settings of the output histogram.
552 htmp_proj->SetDirectory(0);
553 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
560 // -----------------------------------------------------------------------
561 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFProjection(Int_t charge, Int_t species) {
563 // Returns a TObjArray with all TOF histograms (as needed by AliSpectrumDiHadronPID)
564 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
566 TObjArray* aout = new TObjArray(GetNPtBinsPID());
567 aout->SetOwner(kTRUE);
568 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
569 aout->AddLast((TH1F*)GetHistDataTOFProjection(charge, species, iPtBin));
576 // -----------------------------------------------------------------------
577 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
579 // Returns a projection in TOF of the mismatch histogram.
580 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
582 Int_t ptclass = GetPtClass(ptbin);
583 if (ptclass == -1) return 0x0;
584 Int_t bininptclass = GetBinInPtClass(ptbin);
586 // Retrieve original the 2D histogram (we don't own it).
587 TH2F* htmp = (TH2F*)GetHistData(Form("fHistTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
589 // Make the projection.
590 TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojectionMismatch_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
592 // Some settings of the output histogram.
593 htmp_proj->SetDirectory(0);
594 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
601 // -----------------------------------------------------------------------
602 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFMismatch(Int_t charge, Int_t species) {
604 // Returns a TObjArray with all TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
605 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
607 TObjArray* aout = new TObjArray(GetNPtBinsPID());
608 aout->SetOwner(kTRUE);
609 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
610 aout->AddLast((TH1F*)GetHistDataTOFMismatch(charge, species, iPtBin));
617 // -----------------------------------------------------------------------
618 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
620 // Returns a projection in TOF of the data histogram.
621 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
623 Int_t ptclass = GetPtClass(ptbin);
624 if (ptclass == -1) return 0x0;
625 Int_t bininptclass = GetBinInPtClass(ptbin);
627 // Retrieve original the 3D histogram (we don't own it).
628 TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
630 // Make the projection.
631 TAxis* ptaxis = htmp->GetXaxis();
632 Int_t NbinsPt = ptaxis->GetNbins();
633 //cout<<"bin in pt class: "<<bininptclass<<endl;
634 ptaxis->SetRange(bininptclass, bininptclass);
635 TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
636 htmp_proj->SetName(Form("TPCTOFprojection_%i_%i_%i",charge,species,ptbin));
638 //cout<<"projection: "<<htmp_proj<<endl;
640 // Some settings of the output histogram.
641 htmp_proj->SetDirectory(0);
642 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
645 // Putting back the range on the p_T axis.
646 ptaxis->SetRange(1, NbinsPt);
652 // -----------------------------------------------------------------------
653 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFProjection(Int_t charge, Int_t species) {
655 // TODO: This is basically a copy of GetDataTOFProjection -> Make into one function.
657 // Returns a TObjArray with all TPC/TOF histograms (as needed by AliSpectrumDiHadronPID)
658 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
660 TObjArray* aout = new TObjArray(GetNPtBinsPID());
661 aout->SetOwner(kTRUE);
662 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
664 aout->AddLast((TH2F*)GetHistDataTPCTOFProjection(charge, species, iPtBin));
665 //cout<<"aout entries: "<<aout->GetEntriesFast()<<endl;
672 // -----------------------------------------------------------------------
673 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
675 // Returns a projection in TPC/TOF of the mismatch histogram.
676 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
678 Int_t ptclass = GetPtClass(ptbin);
679 if (ptclass == -1) return 0x0;
680 Int_t bininptclass = GetBinInPtClass(ptbin);
682 // Retrieve original the 3D histogram (we don't own it).
683 TH3F* htmp = (TH3F*)GetHistData(Form("fHistTPCTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
685 // Make the projection.
686 TAxis* ptaxis = htmp->GetXaxis();
687 Int_t NbinsPt = ptaxis->GetNbins();
688 ptaxis->SetRange(bininptclass, bininptclass);
689 TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
690 htmp_proj->SetName(Form("TPCTOFprojectionMismatch_%i_%i_%i",charge,species,ptbin));
692 // Some settings of the output histogram.
693 htmp_proj->SetDirectory(0);
694 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
697 // Putting back the range on the p_T axis.
698 ptaxis->SetRange(1, NbinsPt);
704 // -----------------------------------------------------------------------
705 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFMismatch(Int_t charge, Int_t species) {
707 // Returns a TObjArray with all TPC/TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
708 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
710 TObjArray* aout = new TObjArray(GetNPtBinsPID());
711 aout->SetOwner(kTRUE);
712 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
713 aout->AddLast((TH2F*)GetHistDataTPCTOFMismatch(charge, species, iPtBin));
720 // -----------------------------------------------------------------------
721 Double_t AliAODTrackCutsDiHadronPID::GetPtMin(Int_t bin) const {
723 // Same as: TAxis::GetBinLowEdge()
725 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
727 if ((bin < 1) || (bin > fNPtBins + 1)) {
728 cout<<"Bin is out of range..."<<endl;
731 return fPtAxis[bin - 1];
736 // -----------------------------------------------------------------------
737 Double_t AliAODTrackCutsDiHadronPID::GetPtMax(Int_t bin) const {
739 // Same as: TAxis::GetBinUpEdge()
741 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
743 if ((bin < 1) || (bin > fNPtBins + 1)) {
744 cout<<"Bin is out of range..."<<endl;
752 // -----------------------------------------------------------------------
753 Int_t AliAODTrackCutsDiHadronPID::GetNPtBinsPID(Int_t ptclass) const {
755 // Returns the number of pt bins that are used in every "pt class",
756 // where a "pt class" is a range in pT which have the same range in
757 // TOF/ TPC. If class == -1, then the retrun value is the total number
758 // of bins in all classes together.
760 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
763 Int_t nptbinspid = 0;
764 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
765 nptbinspid += fNPtBinsPID[iPtClass];
768 } else if (ptclass >= 0 && ptclass < 5) {
769 return fNPtBinsPID[ptclass];
776 // -----------------------------------------------------------------------
777 Double_t* AliAODTrackCutsDiHadronPID::GetPtAxisPID() const {
779 // Returns an array representing the pT axis for the PID
782 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
784 const Int_t nptbinspid = GetNPtBinsPID();
785 Double_t* ptaxis = new Double_t[nptbinspid + 1];
787 for (Int_t iPtBin = 0; iPtBin < nptbinspid; iPtBin++) {
788 ptaxis[iPtBin] = GetPtMinPID(iPtBin + 1);
791 ptaxis[nptbinspid] = GetPtMaxPID(nptbinspid);
796 // -----------------------------------------------------------------------
797 Double_t AliAODTrackCutsDiHadronPID::GetPtMinPID(Int_t bin) const {
799 // Same as: TAxis::GetBinLowEdge()
801 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
803 Int_t ptclass = GetPtClass(bin);
805 if (ptclass == -1) {return -999.;}
807 Int_t bininptclass = GetBinInPtClass(bin);
809 Double_t minpt = fPtBoundaryPID[ptclass];
810 Double_t maxpt = fPtBoundaryPID[ptclass+1];
811 Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
813 return (minpt + ptres * ((Double_t)(bininptclass - 1)) );
817 // -----------------------------------------------------------------------
818 Double_t AliAODTrackCutsDiHadronPID::GetPtMaxPID(Int_t bin) const {
820 // Same as: TAxis::GetBinUpEdge()
822 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
824 Int_t ptclass = GetPtClass(bin);
826 if (ptclass == -1) {return -999.;}
828 Int_t bininptclass = GetBinInPtClass(bin);
830 Double_t minpt = fPtBoundaryPID[ptclass];
831 Double_t maxpt = fPtBoundaryPID[ptclass+1];
832 Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
834 return (minpt + ptres * ((Double_t)(bininptclass)) );
838 // -----------------------------------------------------------------------
839 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMin(Int_t ptclass) const {
841 // Returns the minimum p_T of a p_T class.
843 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
845 if (ptclass >= 0 && ptclass < 5) {
846 return fPtBoundaryPID[ptclass];
853 // -----------------------------------------------------------------------
854 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMax(Int_t ptclass) const {
856 // Returns the maximum p_T of a p_T class.
858 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
860 if (ptclass >= 0 && ptclass < 5) {
861 return fPtBoundaryPID[ptclass+1];
868 // -----------------------------------------------------------------------
869 Bool_t AliAODTrackCutsDiHadronPID::RequestQAHistos(Int_t histoclass, Bool_t Enable3DSpectra, Bool_t EnablePIDHistos) {
871 // Request certain histograms to be filled.
873 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
875 if ((histoclass > -1) && (histoclass < 12)) {
876 fHistRequested[histoclass] = kTRUE;
877 f3DSpectraEnabeled[histoclass] = Enable3DSpectra;
878 fPIDHistosEnabeled[histoclass] = EnablePIDHistos;
879 //cout<<"histoclass: "<<histoclass<<" requested: "<<fHistRequested[histoclass]<<endl;
886 // -----------------------------------------------------------------------
887 void AliAODTrackCutsDiHadronPID::SetPtRange(Double_t minpt, Double_t maxpt) {
889 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
896 // -----------------------------------------------------------------------
897 void AliAODTrackCutsDiHadronPID::SetFilterMask(UInt_t filtermask) {
899 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
901 fFilterMask = filtermask;
902 fTestFilterMask = kTRUE;
905 // -----------------------------------------------------------------------
906 void AliAODTrackCutsDiHadronPID::SetMaxEta(Double_t maxeta) {
908 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
914 // -----------------------------------------------------------------------
915 void AliAODTrackCutsDiHadronPID::SetMaxRapidity(Double_t maxrapidity) {
917 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
919 fMaxRapidity = maxrapidity;
920 fTestMaxRapidity = kTRUE;
923 // -----------------------------------------------------------------------
924 void AliAODTrackCutsDiHadronPID::SetDemandNoMismatch() {
926 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
928 fTestTOFmismatch = kTRUE;
931 // -----------------------------------------------------------------------
932 void AliAODTrackCutsDiHadronPID::SetDemandFlags(ULong_t demandedflags) {
934 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
936 fDemandedFlags = demandedflags;
940 // -----------------------------------------------------------------------
941 void AliAODTrackCutsDiHadronPID::SetMinimumNumberOfTPCClusters(Int_t minimumnumberoftpcclusters) {
943 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
945 fMinimumNumberOfTPCClusters = minimumnumberoftpcclusters;
946 fTestNumberOfTPCClusters = kTRUE;
949 // -----------------------------------------------------------------------
950 void AliAODTrackCutsDiHadronPID::SetDemandSPDCluster() {
952 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
957 // -----------------------------------------------------------------------
958 void AliAODTrackCutsDiHadronPID::SetPtDeptDCACut(TFormula* DCAxyCutFormula, Double_t DCAzCut, UInt_t MinSPDHits) {
960 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
962 fPtDeptDCAxyCutFormula = DCAxyCutFormula;
964 fMinSPDHitsForPtDeptDCAcut = MinSPDHits;
965 fTestPtDeptDCAcut = kTRUE;
968 // -----------------------------------------------------------------------
969 void AliAODTrackCutsDiHadronPID::StartNewEvent() {
971 // FIXME: This method is now only suited for Data (3 histo classes only.)
972 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
974 // Resetting the counters.
975 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
976 fNTracks[iHistoClass] = 0;
980 // -----------------------------------------------------------------------
981 void AliAODTrackCutsDiHadronPID::EventIsDone(Bool_t IsMC) {
983 // FIXME: This method is now only suited for Data (3 histo classes only.)
984 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
986 // Fill NTracks histos.
988 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
990 if (fHistRequested[iHistoClass]) {
993 // THE FOLLOWING SHOULD NEVER HAPPEN.
994 //if (!fHistPrimRecMCPt[iHistoClass]) InitializeRecMCHistos(iHistoClass);
995 fHistPrimRecNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
997 // THE FOLLOWING SHOULD NEVER HAPPEN.
998 //if (!fHistDataPt[iHistoClass]) InitializeDataHistos(iHistoClass);
999 fHistDataNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
1005 // -----------------------------------------------------------------------
1006 void AliAODTrackCutsDiHadronPID::CreateHistos() {
1008 // Function should be called by the UserCreateOutput() function of the
1009 // analysis task. This function then generates all the histograms that
1010 // were requested locally on the workernode. Even if case that the
1011 // histograms are not filled, it is imperative that they are still
1012 // created, and in the same order, otherwise problems with merging will
1015 // TODO: In principle it should never happen that if this method is called,
1016 // that the lists of QA histograms already exist, so this may become a fatal error
1017 // instead of a warning.
1019 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1023 fTOFMatchingStat = new TH1F("fTOFMatchingStat","fTOFMatchingStat",3,0,3);
1024 (fTOFMatchingStat->GetXaxis())->SetBinLabel(1,"Match");
1025 (fTOFMatchingStat->GetXaxis())->SetBinLabel(2,"Mismatch");
1026 (fTOFMatchingStat->GetXaxis())->SetBinLabel(3,"No TOF hit");
1028 if (!fPrimGenMCTrackQAHistos) {
1029 cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Gen. MC Track QA TList..."<<endl;
1030 fPrimGenMCTrackQAHistos = new TList();
1031 fPrimGenMCTrackQAHistos->SetName("PrimGenMCTrackQAHistos");
1032 fPrimGenMCTrackQAHistos->SetOwner(kTRUE);
1033 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Gen. MC Track QA TList was already created..."<<endl;}
1035 if (!fSecGenMCTrackQAHistos) {
1036 cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Gen. MC Track QA TList..."<<endl;
1037 fSecGenMCTrackQAHistos = new TList();
1038 fSecGenMCTrackQAHistos->SetName("SecGenMCTrackQAHistos");
1039 fSecGenMCTrackQAHistos->SetOwner(kTRUE);
1040 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {if (fHistRequested[iHistoClass]) {InitializeGenMCHistos(iHistoClass);}}
1041 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Gen. MC Track QA TList was already created..."<<endl;}
1043 if (!fPrimRecMCTrackQAHistos) {
1044 cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Rec. MC Track QA TList..."<<endl;
1045 fPrimRecMCTrackQAHistos = new TList();
1046 fPrimRecMCTrackQAHistos->SetName("PrimRecMCTrackQAHistos");
1047 fPrimRecMCTrackQAHistos->SetOwner(kTRUE);
1048 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Rec. MC Track QA TList was already created..."<<endl;}
1050 if (!fSecRecMCTrackQAHistos) {
1051 cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Rec. MC Track QA TList..."<<endl;
1052 fSecRecMCTrackQAHistos = new TList();
1053 fSecRecMCTrackQAHistos->SetName("SecRecMCTrackQAHistos");
1054 fSecRecMCTrackQAHistos->SetOwner(kTRUE);
1055 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1056 if (fHistRequested[iHistoClass]) {
1057 InitializeRecMCHistos(iHistoClass);
1060 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Rec. MC Track QA TList was already created..."<<endl;}
1064 if (!fDataTrackQAHistos) {
1065 cout<<"AliAODTrackCutsDiHadronPID - Creating Data Track QA TList..."<<endl;
1066 fDataTrackQAHistos = new TList();
1067 fDataTrackQAHistos->SetName("DataTrackQAHistos");
1068 fDataTrackQAHistos->SetOwner(kTRUE);
1070 // Add general histograms.
1071 fHistAcceptedFilterBits = InitializeAcceptedFilterBits("fHistAcceptedFilterBits");
1072 fDataTrackQAHistos->Add(fHistAcceptedFilterBits);
1074 // Add histograms per class.
1075 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1076 fHistDataDCAxyOneSigma[iHistoClass] = InitializeDCASpectrum("fHistDataDCAxyOneSigma",iHistoClass);
1077 fDataTrackQAHistos->Add(fHistDataDCAxyOneSigma[iHistoClass]);
1078 if (fHistRequested[iHistoClass]) {
1079 if (iHistoClass < 3) {InitializeDataHistos(iHistoClass);}
1082 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Data QA TList was already created..."<<endl;}
1087 // -----------------------------------------------------------------------
1088 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime) {
1091 // Checks performed on a data track.
1094 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1096 if (!track) return kFALSE;
1098 if (!fDataTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1100 // Check the track cuts (NB. check function return kTRUE if track is accepted.)
1101 if (!CheckPt(track->Pt())) return kFALSE;
1102 if (!CheckMaxEta(track->Eta())) return kFALSE;
1103 if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
1104 if (!CheckFlags(track->GetFlags())) return kFALSE;
1105 if (!CheckNclsTPC(track->GetNclsTPC())) return kFALSE;
1106 if (!CheckTOFmismatch(track->IsTOFMismatch())) return kFALSE;
1109 if (track->HasPointOnITSLayer(0)) NSPDhits++;
1110 if (track->HasPointOnITSLayer(1)) NSPDhits++;
1111 if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
1112 if (fTestSPDAny) {if (NSPDhits < 1) return kFALSE;}
1114 // Fill the filterbit histogram.
1115 for (Int_t iBin = (fRelevantBitsArray->GetSize() - 1); iBin >= 0; --iBin) {
1116 if ( (((Int_t)track->GetFilterMap())&(1<<fRelevantBitsArray->At(iBin))) == (1<<fRelevantBitsArray->At(iBin))) {
1117 fHistAcceptedFilterBits->Fill(iBin);
1122 // Track has passed the cuts, fill QA histograms.
1123 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
1125 if (fHistRequested[iHistoClass]) {
1127 // Check the charge (could be neater).
1128 if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
1129 if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
1130 if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
1132 FillDataHistos(iHistoClass, track);
1134 // Ignore if random hit is < -1.e20.
1135 if (randomhittime > -1.e20) FillTOFMismatchHistos(iHistoClass,track,randomhittime);
1137 fNTracks[iHistoClass]++;
1142 // Track Has Passed.
1147 // -----------------------------------------------------------------------
1148 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedGeneratedMC(AliAODMCParticle* particle) {
1151 // Checks performed on a generated MC particle.
1154 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1156 // PDG codes for particles:
1157 // pi+ 211, K+ 321, p 2212
1159 if (!particle) return kFALSE;
1161 if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1163 // Check the track cuts.
1164 if (!CheckPt(particle->Pt())) return kFALSE;
1165 if (!CheckMaxEta(particle->Eta())) return kFALSE;
1166 if (!CheckRapidity(particle->Y())) return kFALSE; // NEW: rapidity cut. (not done on data)
1168 // NOT YET IMPLEMENTED! - Histoclasses for specific particles.
1169 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1171 if (fHistRequested[iHistoClass]) {
1173 // Check the charge (could be neater).
1174 if ((iHistoClass == 0) && (particle->Charge() == 0)) continue;
1175 if ((iHistoClass == 1) && (particle->Charge() <= 0)) continue;
1176 if ((iHistoClass == 2) && (particle->Charge() >= 0)) continue;
1177 if ((iHistoClass == 3) && (TMath::Abs(particle->GetPdgCode()) != 211)) continue;
1178 if ((iHistoClass == 4) && (particle->GetPdgCode()) != 211) continue;
1179 if ((iHistoClass == 5) && (particle->GetPdgCode()) != -211) continue;
1180 if ((iHistoClass == 6) && (TMath::Abs(particle->GetPdgCode()) != 321)) continue;
1181 if ((iHistoClass == 7) && (particle->GetPdgCode()) != 321) continue;
1182 if ((iHistoClass == 8) && (particle->GetPdgCode()) != -321) continue;
1183 if ((iHistoClass == 9) && (TMath::Abs(particle->GetPdgCode()) != 2212)) continue;
1184 if ((iHistoClass == 10) && (particle->GetPdgCode()) != 2212) continue;
1185 if ((iHistoClass == 11) && (particle->GetPdgCode()) != -2212) continue;
1187 // Secondary specification not set.
1188 //cout<<"Charge: "<<particle->Charge()<<" PDG code: "<<particle->GetPdgCode()<<endl;
1189 //cout<<particle->IsPhysicalPrimary()<<" "<<particle->IsSecondaryFromWeakDecay()<<" "<<particle->IsSecondaryFromMaterial()<<endl;
1191 // These two functions are not implemented...
1192 //if (particle->IsSecondaryFromWeakDecay()) cout<<"Secondary From Weak Decay!"<<endl;
1193 //if (particle->IsSecondaryFromMaterial()) cout<<"Secondary From Material!"<<endl;
1195 FillGenMCHistos(iHistoClass, particle);
1200 // Particle has Passed.
1205 // -----------------------------------------------------------------------
1206 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedReconstructedMC(AliTrackDiHadronPID* track) {
1208 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1210 if (!track) return kFALSE;
1212 if (!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1214 // Check the track cuts.
1215 if (!CheckPt(track->MCPt())) return kFALSE; // Kinematic cuts are done on the MC particles
1216 if (!CheckMaxEta(track->MCEta())) return kFALSE; // and not on the data.
1217 if (!CheckRapidity(track->MCY())) return kFALSE; // NEW: rapidity cut.
1218 if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
1219 if (!CheckFlags(track->GetFlags())) return kFALSE;
1220 if (!CheckNclsTPC(track->GetNclsTPC())) return kFALSE;
1223 if (track->HasPointOnITSLayer(0)) NSPDhits++;
1224 if (track->HasPointOnITSLayer(1)) NSPDhits++;
1225 if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
1226 if (fTestSPDAny) {if (NSPDhits < 1) return kFALSE;}
1228 // Track has passed the cuts, fill QA histograms.
1229 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1231 if (fHistRequested[iHistoClass]) {
1233 // Check the charge (could be neater).
1234 if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
1235 if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
1236 if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
1237 if ((iHistoClass == 3) && (TMath::Abs(track->GetPdgCode()) != 211)) continue;
1238 if ((iHistoClass == 4) && (track->GetPdgCode()) != 211) continue;
1239 if ((iHistoClass == 5) && (track->GetPdgCode()) != -211) continue;
1240 if ((iHistoClass == 6) && (TMath::Abs(track->GetPdgCode()) != 321)) continue;
1241 if ((iHistoClass == 7) && (track->GetPdgCode()) != 321) continue;
1242 if ((iHistoClass == 8) && (track->GetPdgCode()) != -321) continue;
1243 if ((iHistoClass == 9) && (TMath::Abs(track->GetPdgCode()) != 2212)) continue;
1244 if ((iHistoClass == 10) && (track->GetPdgCode()) != 2212) continue;
1245 if ((iHistoClass == 11) && (track->GetPdgCode()) != -2212) continue;
1247 FillRecMCHistos(iHistoClass, track);
1249 fNTracks[iHistoClass]++;
1253 // Track Has Passed.
1258 // -----------------------------------------------------------------------
1259 Int_t AliAODTrackCutsDiHadronPID::GetPtClass(Int_t ptbin) const {
1261 // Returns a p_T class as a function of bin (PID histos)
1262 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1264 Int_t currentptclass = 0;
1265 Int_t currentptbin = fNPtBinsPID[0];
1267 while (currentptbin < ptbin) {
1269 if (currentptclass == 5) {break;}
1270 currentptbin += fNPtBinsPID[currentptclass];
1273 if (currentptclass == 5) {
1274 cout<<"GetPtClass -> ptbin out of range!"<<endl;
1278 return currentptclass;
1281 // -----------------------------------------------------------------------
1282 Int_t AliAODTrackCutsDiHadronPID::GetBinInPtClass(Int_t ptbin) const {
1284 // Returns the bin withing the p_T class (i.e., 1..Nbins)
1285 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1287 Int_t ptclass = GetPtClass(ptbin);
1288 if (ptclass == -1) {return -1;}
1290 Int_t ptbinout = ptbin;
1291 for (Int_t iPtClass = 0; iPtClass < ptclass; iPtClass++) {ptbinout -= fNPtBinsPID[iPtClass];}
1297 // -----------------------------------------------------------------------
1298 Bool_t AliAODTrackCutsDiHadronPID::CheckPt(Double_t pt) const {
1300 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1302 if (!fTestPt) return kTRUE;
1303 if ((pt > fMinPt) && (pt < fMaxPt)) return kTRUE;
1307 // -----------------------------------------------------------------------
1308 Bool_t AliAODTrackCutsDiHadronPID::CheckMaxEta(Double_t eta) const {
1310 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1312 if (!fTestMaxEta) return kTRUE; // Accepted if there is no check on this parameter.
1313 if (TMath::Abs(eta) < fMaxEta) return kTRUE;
1317 // -----------------------------------------------------------------------
1318 Bool_t AliAODTrackCutsDiHadronPID::CheckRapidity(Double_t rap) const {
1320 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1322 if (!fTestMaxRapidity) return kTRUE;
1323 if (TMath::Abs(rap) < fMaxRapidity) return kTRUE;
1327 // -----------------------------------------------------------------------
1328 Bool_t AliAODTrackCutsDiHadronPID::CheckFilterMask(UInt_t filtermap) const {
1330 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1332 if (!fTestFilterMask) return kTRUE;
1333 if (fFilterMask & filtermap) return kTRUE;
1337 // -----------------------------------------------------------------------
1338 Bool_t AliAODTrackCutsDiHadronPID::CheckFlags(ULong_t flags) const {
1340 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1342 if (!fTestFlags) return kTRUE;
1343 if ((flags & fDemandedFlags) == fDemandedFlags) return kTRUE;
1347 // -----------------------------------------------------------------------
1348 Bool_t AliAODTrackCutsDiHadronPID::CheckNclsTPC(Int_t ncls) const {
1350 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1352 if (!fTestNumberOfTPCClusters) return kTRUE;
1353 if (ncls > fMinimumNumberOfTPCClusters) return kTRUE;
1357 // -----------------------------------------------------------------------
1358 Bool_t AliAODTrackCutsDiHadronPID::CheckTOFmismatch(Bool_t ismismatch) const {
1360 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1362 if (!fTestTOFmismatch) return kTRUE; // if we're not cutting on mismatch, then it's accepted.
1363 if (!ismismatch) return kTRUE; // so if the track is not a mismatch, then it is accepted.
1364 return kFALSE; // if it is a mismatch, then it's not accepted.
1367 // -----------------------------------------------------------------------
1368 Bool_t AliAODTrackCutsDiHadronPID::CheckPtDeptDCACut(Double_t dcaz, Double_t dcaxy, Double_t pt, UInt_t SPDhits) const {
1370 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1372 if (!fTestPtDeptDCAcut) return kTRUE;
1373 if (SPDhits < fMinSPDHitsForPtDeptDCAcut) return kTRUE; // If there are not enough SPD hits to do the cut.
1374 if ((dcaz < fDCAzCut) && (dcaxy < fPtDeptDCAxyCutFormula->Eval(pt))) return kTRUE;
1378 // -----------------------------------------------------------------------
1379 Bool_t AliAODTrackCutsDiHadronPID::FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1381 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1383 // Fill the histograms.
1384 fHistDataPt[histoclass]->Fill(track->Pt());
1385 fHistDataPhiEtaPt[histoclass]->Fill(track->Phi(),track->Eta(),track->Pt());
1387 // Fill DCA histograms.
1388 fHistDataDCAxy[histoclass]->Fill(track->GetXYAtDCA());
1389 fHistDataDCAz[histoclass]->Fill(track->GetZAtDCA());
1394 Introduced a separate selection mechanism here as the TPC is off for low pt for
1395 in particular protons and kaons. Therefor use only the TOF for identifying under 1.8 GeV.
1396 Setting this with fLowPtNSigmaTOFOnly flag.
1398 // histoclass = 0: all charges; histoclass = 1: positive; histoclass = 2: negative;
1399 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(0) * track->GetNumberOfSigmasTOF(0) +
1400 track->GetNumberOfSigmasTPC(0) * track->GetNumberOfSigmasTPC(0)) < 1.) {
1401 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1402 fHistDataDCAxyOneSigma[3 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Pions.
1403 // checkSum++; cout<<"Pion found: nSigTOF: "<<track->GetNumberOfSigmasTOF(0)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(0)<<endl;
1406 // for protons and low pt, only when fLowPtNSigmaTOFOnly set to kTRUE:
1407 if ((TMath::Abs(track->Pt()) < 1.8) && fLowPtNSigmaTOFOnly) {
1409 if (TMath::Abs(track->GetNumberOfSigmasTOF(1)) < 1.) {
1410 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1411 fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Kaons.
1413 if (TMath::Abs(track->GetNumberOfSigmasTOF(2)) < 1.) {
1414 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1415 fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Protons.
1416 // checkSum++; cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
1419 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(1) * track->GetNumberOfSigmasTOF(1) +
1420 track->GetNumberOfSigmasTPC(1) * track->GetNumberOfSigmasTPC(1)) < 1.) {
1421 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1422 fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Kaons.
1423 // checkSum++; cout<<"Kaon found: nSigTOF: "<<track->GetNumberOfSigmasTOF(1)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(1)<<endl;
1425 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(2) * track->GetNumberOfSigmasTOF(2) +
1426 track->GetNumberOfSigmasTPC(2) * track->GetNumberOfSigmasTPC(2)) < 1.) {
1427 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1428 fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Protons.
1429 // checkSum++; cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
1433 // check for double identification (or triple..)
1434 if(checkSum > 1) {AliError("More than one particle identified for the same track!"); }
1437 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1439 // Note that a possible Y cut is only done on the PID cuts!
1440 if (!CheckRapidity(track->Y(iSpecies))) continue;
1442 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1443 fHistDataPID[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),
1444 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)),
1445 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)));
1453 // -----------------------------------------------------------------------
1454 Bool_t AliAODTrackCutsDiHadronPID::FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime) {
1456 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1458 // Fill TOF mismatch.
1459 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1461 // Note that a possible Y cut is only done on the PID cuts!
1462 if (!CheckRapidity(track->Y(iSpecies))) continue;
1464 if (!track->GetNumberOfSigmasTOF(iSpecies)) {
1465 cout<<"NSigma TOF: "<<track->GetNumberOfSigmasTOF(iSpecies)<<endl;
1466 cout<<"TOFsignal exp: "<<track->GetTOFsignalExpected(iSpecies)<<endl;
1467 cout<<"TOFsigma exp: "<<track->GetTOFsigmaExpected(iSpecies)<<endl;
1472 Double_t TOFmismatchSigma = randomhittime - track->GetTOFsignalExpected(iSpecies);
1473 if (fUseNSigmaOnPIDAxes) {
1474 if (track->GetTOFsigmaExpected(iSpecies) < 10e-30) {
1476 cout << "ERROR: division through (almost) zero on the mismatch signal..."<<endl;
1477 cout << "Random time: "<<randomhittime<<endl;
1478 cout << "Exp TOF signal track for species "<<iSpecies<<" "<<track->GetTOFsignalExpected(iSpecies)<<endl;
1479 cout << "Exp Mismatch signal assuming species "<< iSpecies<<": "<<TOFmismatchSigma<<endl;
1480 cout << "Expected TOF sigma: "<<track->GetTOFsigmaExpected(iSpecies)<<endl;
1481 ULong_t flags_kTOFout_kTIME = (UInt_t)(AliAODTrack::kTOFout)|(UInt_t)(AliAODTrack::kTIME);
1482 cout << "Random track TOF status (kTOFout, kTIME): "<<(((track->GetFlags())&flags_kTOFout_kTIME)==flags_kTOFout_kTIME)<<endl;
1485 TOFmismatchSigma /= track->GetTOFsigmaExpected(iSpecies);
1489 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1490 fHistTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(), TOFmismatchSigma);
1491 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(), TOFmismatchSigma,
1492 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)));
1500 // -----------------------------------------------------------------------
1501 Bool_t AliAODTrackCutsDiHadronPID::FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle) {
1503 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1505 //cout << "histoclass: "<<histoclass<<" particle: "<<particle<<endl;
1507 // Fill the histograms.
1508 if (particle->IsPhysicalPrimary()) {
1509 //cout<<"Trying to fill: "<<fHistPrimGenMCPt[histoclass]->GetName()<<" of type: " <<fHistPrimGenMCPt[histoclass]->ClassName()<<" at: "<<fHistPrimGenMCPt[histoclass]<<endl;
1510 fHistPrimGenMCPt[histoclass]->Fill(particle->Pt());
1511 //cout<< "phi: "<<particle->Phi()<< " eta: "<<particle->Eta() << " pt: "<<particle->Pt()<<" ";
1512 fHistPrimGenMCPhiEtaPt[histoclass]->Fill(particle->Phi(), particle->Eta(), particle->Pt());
1513 //cout<<"OK!"<<endl;
1515 //cout<<"Trying to fill: "<<fHistPrimGenMCPt[histoclass]->GetName()<<" of type: " <<fHistPrimGenMCPt[histoclass]->ClassName()<<" at: "<<fHistPrimGenMCPt[histoclass]<<endl;
1516 fHistSecGenMCPt[histoclass]->Fill(particle->Pt());
1517 //cout<< "phi: "<<particle->Phi()<< " eta: "<<particle->Eta() << " pt: "<<particle->Pt()<<" ";
1518 fHistSecGenMCPhiEtaPt[histoclass]->Fill(particle->Phi(), particle->Eta(), particle->Pt());
1519 //cout<<"OK!"<<endl;
1526 // -----------------------------------------------------------------------
1527 Bool_t AliAODTrackCutsDiHadronPID::FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1529 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1531 // Fill the Pt and acceptance histograms.
1532 if (track->IsPhysicalPrimary()) {
1533 fHistPrimRecMCPt[histoclass]->Fill(track->MCPt());
1534 fHistPrimRecPtGenPt[histoclass]->Fill(track->MCPt(), track->Pt());
1535 fHistPrimRecMCPhiEtaPt[histoclass]->Fill(track->MCPhi(), track->MCEta(), track->MCPt());
1537 fHistSecRecMCPt[histoclass]->Fill(track->MCPt());
1538 fHistSecRecMCPhiEtaPt[histoclass]->Fill(track->MCPhi(), track->MCEta(), track->MCPt());
1541 // Fill the DCA histograms.
1542 if (track->IsPhysicalPrimary()) {
1543 fHistPrimRecMCDCA[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1545 if (track->IsSecondaryFromMaterial()) {
1546 fHistSecRecMCDCAMat[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1548 if (track->IsSecondaryFromWeakDecay()) {
1549 fHistSecRecMCDCAWeak[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1553 if (histoclass < 3) {
1555 // Note that the following histogram is only sensible if the fTOFlabel is set properly.
1556 // If not, the histogram will be filled with "no match".
1557 fTOFMatchingStat->Fill(((Double_t)track->GetTOFMatchingStatus())+0.5);
1559 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1561 // Note that a possible Y cut is only done on the PID cuts!
1562 if (!CheckRapidity(track->MCY())) continue;
1564 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1565 //cout << "recpt: " << track->Pt() << " mcpt: "<<track->MCPt() << " DTOF: " <<
1566 //(fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)) << " DTPC: " <<
1567 //(fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTPC(iSpecies) : track->GetTPCsignalMinusExpected(iSpecies)) << endl;
1569 fHistPrimRecPID[histoclass][iSpecies][iPtClass]->Fill(track->MCPt(),
1570 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)));
1572 if (track->IsTOFMismatch()) {
1573 fHistPrimRecMismatch[histoclass][iSpecies][iPtClass]->Fill(track->MCPt(),
1574 (fUseNSigmaOnPIDAxes ? track->GetNumberOfSigmasTOF(iSpecies) : track->GetTOFsignalMinusExpected(iSpecies)));
1584 // -----------------------------------------------------------------------
1585 Bool_t AliAODTrackCutsDiHadronPID::InitializeDataHistos(Int_t histoclass) {
1587 cout<<"AliAODTrackCutsDiHadronPID - Creating Data Histograms of Class "<<histoclass<<endl;
1588 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1590 // Add the Pt spectra.
1591 fHistDataPt[histoclass] = InitializePtSpectrum("fHistDataPt",histoclass);
1592 fDataTrackQAHistos->Add(fHistDataPt[histoclass]);
1593 fHistDataPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistDataPhiEtaPt",histoclass);
1594 fDataTrackQAHistos->Add(fHistDataPhiEtaPt[histoclass]);
1596 // Add the NTrack histograms.
1597 fHistDataNTracks[histoclass] = InitializeNTracksHisto("fHistDataNTracks",histoclass);
1598 fDataTrackQAHistos->Add(fHistDataNTracks[histoclass]);
1600 // Add the DCA histograms.
1601 fHistDataDCAxy[histoclass] = InitializeDCAxyHisto("fHistDCAxy",histoclass);
1602 fDataTrackQAHistos->Add(fHistDataDCAxy[histoclass]);
1603 fHistDataDCAz[histoclass] = InitializeDCAzHisto("fHistDCAz",histoclass);
1604 fDataTrackQAHistos->Add(fHistDataDCAz[histoclass]);
1606 // Add the PID histograms. (FIXME shoudl be able to turn the creation of these histos off.)
1607 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1608 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1610 fHistDataPID[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistDataPID",histoclass,iSpecies,iPtClass);
1611 fDataTrackQAHistos->Add(fHistDataPID[histoclass][iSpecies][iPtClass]);
1613 fHistTOFMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistTOFMismatch",histoclass,iSpecies,iPtClass);
1614 fDataTrackQAHistos->Add(fHistTOFMismatch[histoclass][iSpecies][iPtClass]);
1616 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistTPCTOFMismatch",histoclass,iSpecies,iPtClass);
1617 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->SetTitle(Form("PID %s (Exp: %s)",fHistoName[histoclass].Data(),fParticleName[iSpecies].Data()));
1618 fDataTrackQAHistos->Add(fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]);
1627 // -----------------------------------------------------------------------
1628 Bool_t AliAODTrackCutsDiHadronPID::InitializeGenMCHistos(Int_t histoclass) {
1630 cout<<"AliAODTrackCutsDiHadronPID - Creating Generated MC Histograms of Class "<<histoclass<<endl;
1631 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1633 // Primary Particles.
1634 fHistPrimGenMCPt[histoclass] = InitializePtSpectrum("fHistPrimGenMCPt",histoclass);
1635 fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPt[histoclass]);
1637 fHistPrimGenMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistPrimGenMCPhiEtaPt",histoclass);
1638 fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPhiEtaPt[histoclass]);
1640 // Secondary Particles.
1641 fHistSecGenMCPt[histoclass] = InitializePtSpectrum("fHistSecGenMCPt",histoclass);
1642 fSecGenMCTrackQAHistos->Add(fHistSecGenMCPt[histoclass]);
1644 fHistSecGenMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistSecGenMCPhiEtaPt",histoclass);
1645 fSecGenMCTrackQAHistos->Add(fHistSecGenMCPhiEtaPt[histoclass]);
1651 // -----------------------------------------------------------------------
1652 Bool_t AliAODTrackCutsDiHadronPID::InitializeRecMCHistos(Int_t histoclass) {
1654 cout<<"AliAODTrackCutsDiHadronPID - Creating Reconstructed MC Histograms of Class "<<histoclass<<endl;
1655 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1657 // Primary Particles.
1658 fHistPrimRecMCPt[histoclass] = InitializePtSpectrum("fHistPrimRecMCPt",histoclass);
1659 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPt[histoclass]);
1661 fHistPrimRecPtGenPt[histoclass] = InitializeRecPtGenPt("fHistPrimRecPtGenPt",histoclass);
1662 fPrimRecMCTrackQAHistos->Add(fHistPrimRecPtGenPt[histoclass]);
1664 fHistPrimRecMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistPrimRecMCPhiEtaPt",histoclass);
1665 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPhiEtaPt[histoclass]);
1667 fHistPrimRecNTracks[histoclass] = InitializeNTracksHisto("fHistPrimRecNTracks",histoclass);
1668 fPrimRecMCTrackQAHistos->Add(fHistPrimRecNTracks[histoclass]);
1670 fHistPrimRecMCDCA[histoclass] = InitializeDCASpectrum("fHistPrimRecDCA",histoclass);
1671 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCDCA[histoclass]);
1673 if (histoclass < 3) {
1674 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1675 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1677 fHistPrimRecPID[histoclass][iSpecies][iPtClass] = InitializeTOFHisto("fHistPrimRecPID",histoclass,iSpecies,iPtClass);
1678 fPrimRecMCTrackQAHistos->Add(fHistPrimRecPID[histoclass][iSpecies][iPtClass]);
1680 fHistPrimRecMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistPrimRecMismatch",histoclass,iSpecies,iPtClass);
1681 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMismatch[histoclass][iSpecies][iPtClass]);
1688 // Secondary Particles.
1689 fHistSecRecMCPt[histoclass] = InitializePtSpectrum("fHistSecRecMCPt",histoclass);
1690 fSecRecMCTrackQAHistos->Add(fHistSecRecMCPt[histoclass]);
1692 fHistSecRecMCPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistSecRecMCPhiEtaPt",histoclass);
1693 fSecRecMCTrackQAHistos->Add(fHistSecRecMCPhiEtaPt[histoclass]);
1695 fHistSecRecMCDCAMat[histoclass] = InitializeDCASpectrum("fHistSecRecDCAMat",histoclass);
1696 fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAMat[histoclass]);
1698 fHistSecRecMCDCAWeak[histoclass] = InitializeDCASpectrum("fHistSecRecDCAWeak",histoclass);
1699 fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAWeak[histoclass]);
1705 // -----------------------------------------------------------------------
1706 TH1F* AliAODTrackCutsDiHadronPID::InitializeAcceptedFilterBits(const char* name) {
1708 // This histogram keeps track of the filtermask of all accepted tracks, projected
1709 // onto the requested filtermask. For example, we requested mask 2 or 4, then this
1710 // histogram will have three bins, 2, 4 and 6. Suppose now a track is found which has
1711 // mask 12, then since (12 & 2) = 0, (12 & 4) = 4, (12 & 6) = 4, the track will end up
1714 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1715 if (!fTestFilterMask || fFilterMask == 0) {cout << Form("%s -> ERROR: No filtermask requested.",__func__) << endl; return 0x0;}
1717 // Determine the labels of the X axis.
1718 SetXaxisAcceptedFilterBits();
1720 // Create histogram.
1721 TH1F* hout = new TH1F(name,"Filtermask of accepted track;Mask;N",fRelevantBitsArray->GetSize(),-0.5,fRelevantBitsArray->GetSize()-0.5);
1722 hout->SetDirectory(0);
1725 TAxis* axistmp = hout->GetXaxis();
1726 for (Int_t iBin = 1; iBin <= axistmp->GetNbins(); ++iBin) {
1727 axistmp->SetBinLabel(iBin, Form("%i",fRelevantBitsArray->At(iBin-1)));
1733 // -----------------------------------------------------------------------
1734 void AliAODTrackCutsDiHadronPID::SetXaxisAcceptedFilterBits() {
1736 // Creates the axis for the AcceptedFilterBits histogram.
1737 // See exercise: "FindAllCombinations.C"
1738 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1740 // Step 1: Find the largest bit, in the requested filtermask,
1741 Int_t largestBit = 0;
1742 Int_t baseArraySizeTmp = 0;
1743 Int_t fullArraySizeTmp = 0;
1745 while (((Int_t)fFilterMask) >= (1<<(largestBit)) ) {
1746 if ((((Int_t)fFilterMask)&(1<<largestBit))==(1<<largestBit)) {
1747 fullArraySizeTmp += AliFunctionsDiHadronPID::Power(2, baseArraySizeTmp);
1754 // Step 2: Create and fill base array.
1755 const Int_t baseArraySize = baseArraySizeTmp;
1756 Int_t* baseArray = new Int_t[baseArraySize];
1757 for (Int_t ii = 0; ii < baseArraySize; ++ii) {baseArray[ii] = 0;}
1759 Int_t iBaseArray = 0;
1760 for (Int_t iBit = 0; iBit <= largestBit; ++iBit) {
1762 if ((((Int_t)fFilterMask)&(1<<iBit))==(1<<iBit)) {
1763 baseArray[iBaseArray] = (1<<iBit);
1769 // Step 3: Create and fill full array.
1770 const Int_t fullArraySize = fullArraySizeTmp;
1771 Int_t* fullArray = new Int_t[fullArraySize];
1772 fullArray[0] = baseArray[0];
1773 Int_t iFullArray = 1;
1775 for (Int_t ii = 1; ii < baseArraySize; ++ii) {
1776 Int_t range = (iFullArray + AliFunctionsDiHadronPID::Power(2,ii));
1777 for (Int_t jj = iFullArray; jj < range; ++jj) {
1779 fullArray[jj] = baseArray[ii];
1781 // Add beginning part of the array:
1782 if (jj!=iFullArray) {
1783 fullArray[jj] += fullArray[jj - iFullArray - 1];
1787 iFullArray += AliFunctionsDiHadronPID::Power(2,ii);
1791 // Step 4: Convert to TArrayI object.
1792 fRelevantBitsArray = new TArrayI(fullArraySize, fullArray);
1794 // Delete the temporary arrays.
1800 // -----------------------------------------------------------------------
1801 TH1F* AliAODTrackCutsDiHadronPID::InitializePtSpectrum(const char* name, Int_t histoclass) {
1803 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1805 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1806 Form("p_{T} Spectrum (%s);p_{T} (GeV/c);Count",fHistoName[histoclass].Data()),fNPtBins,fPtAxis);
1808 hout->SetDirectory(0);
1814 // -----------------------------------------------------------------------
1815 TH2F* AliAODTrackCutsDiHadronPID::InitializeRecPtGenPt(const char* name, Int_t histoclass) {
1817 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1819 TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
1820 Form("p_{T} Rec vs Gen (%s);p_{T,gen} (GeV/c); p_{T,rec}",fHistoName[histoclass].Data()),
1821 fNPtBins,fPtAxis,fNPtBins,fPtAxis);
1823 hout->SetDirectory(0);
1829 // -----------------------------------------------------------------------
1830 TH3F* AliAODTrackCutsDiHadronPID::InitializePhiEtaPt(const char* name, Int_t histoclass) {
1832 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1834 TH3F* hout = AliHistToolsDiHadronPID::MakeHist3D(Form("%s%s",name,fHistoName[histoclass].Data()),
1835 Form("Spectrum (%s);#phi;#eta;p_{T} (GeV/c)",fHistoName[histoclass].Data()),
1836 fNPhiBins,0.,2.*TMath::Pi(),
1837 fNEtaBins,-fMaxEta,fMaxEta,
1840 hout->SetDirectory(0);
1846 // -----------------------------------------------------------------------
1847 TH2F* AliAODTrackCutsDiHadronPID::InitializeDCASpectrum(const char* name, Int_t histoclass) {
1849 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1851 TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
1852 Form("DCA_{xy} (%s); p_{T} (GeV/c); DCA_{xy} (cm)",fHistoName[histoclass].Data()),
1856 hout->SetDirectory(0);
1862 // -----------------------------------------------------------------------
1863 TH1F* AliAODTrackCutsDiHadronPID::InitializeNTracksHisto(const char* name, Int_t histoclass) {
1865 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1867 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1868 Form("Number of Accepted Tracks (%s);N%s;N_{event}",fHistoName[histoclass].Data(),fHistoLatex[histoclass].Data()),
1871 hout->SetDirectory(0);
1877 // -----------------------------------------------------------------------
1878 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAxyHisto(const char* name, Int_t histoclass) {
1880 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1882 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1883 Form("DCAxy (%s);DCAxy (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1885 hout->SetDirectory(0);
1891 // -----------------------------------------------------------------------
1892 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAzHisto(const char* name, Int_t histoclass) {
1894 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1896 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1897 Form("DCAz (%s);DCAz (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1899 hout->SetDirectory(0);
1905 // -----------------------------------------------------------------------
1906 TH3F* AliAODTrackCutsDiHadronPID::InitializeAcceptanceHisto(const char* /*name*/, Int_t /*histoclass*/) {
1908 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1913 // -----------------------------------------------------------------------
1914 TH3F* AliAODTrackCutsDiHadronPID::InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1916 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1918 TString PIDaxeslabel;
1919 if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF};n#sigma_{TPC}";}
1920 else {PIDaxeslabel = ";#Delta t (ps);dE/dx (a.u.)";}
1922 TH3F* hout = new TH3F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1923 Form("PID %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1924 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1925 fTOFbins[ptclass][expspecies],
1926 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1927 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]),
1928 fTPCbins[ptclass][expspecies],
1929 (fUseNSigmaOnPIDAxes ? fTPCLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTPCStd : fTPCLowerBound[ptclass][expspecies]),
1930 (fUseNSigmaOnPIDAxes ? fTPCUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTPCStd : fTPCUpperBound[ptclass][expspecies]));
1932 hout->SetDirectory(0);
1938 // -----------------------------------------------------------------------
1939 TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1941 // Is basically the same as InitializeTOFHisto -> CAN BE REMOVED!
1942 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1944 TString PIDaxeslabel;
1945 if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF}";}
1946 else {PIDaxeslabel = ";#Delta t (ps)";}
1948 TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1949 Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1950 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1951 fTOFbins[ptclass][expspecies],
1952 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1953 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]));
1955 hout->SetDirectory(0);
1961 // -----------------------------------------------------------------------
1962 TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1964 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1966 TString PIDaxeslabel;
1967 if (fUseNSigmaOnPIDAxes) {PIDaxeslabel = ";n#sigma_{TOF}";}
1968 else {PIDaxeslabel = ";#Delta t (ps)";}
1970 TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1971 Form("TOF %s (Exp: %s);p_{T} (GeV/c)%s",fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),PIDaxeslabel.Data()),
1972 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1973 fTOFbins[ptclass][expspecies],
1974 (fUseNSigmaOnPIDAxes ? fTOFLowerBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFLowerBound[ptclass][expspecies]),
1975 (fUseNSigmaOnPIDAxes ? fTOFUpperBound[ptclass][expspecies] / AliTrackDiHadronPID::fSigmaTOFStd : fTOFUpperBound[ptclass][expspecies]));
1977 hout->SetDirectory(0);