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)
29 #include "TIterator.h"
32 #include "AliAODTrack.h"
33 #include "AliAODEvent.h"
34 #include "AliAODVertex.h"
35 #include "AliAODMCParticle.h"
39 #include "AliPIDResponse.h"
40 #include "AliTPCPIDResponse.h"
42 #include "AliTrackDiHadronPID.h"
43 #include "AliHistToolsDiHadronPID.h"
45 #include "AliAODTrackCutsDiHadronPID.h"
49 ClassImp(AliAODTrackCutsDiHadronPID);
51 // -----------------------------------------------------------------------
52 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID():
60 fMinSPDHitsForPtDeptDCAcut(0),
61 fPtDeptDCAxyCutFormula(0x0),
65 fTestFilterMask(kFALSE),
67 fTestMaxRapidity(kFALSE),
69 fTestTOFmismatch(kFALSE),
70 fTestPtDeptDCAcut(kFALSE),
71 fDataTrackQAHistos(0x0),
72 fPrimRecMCTrackQAHistos(0x0),
73 fPrimGenMCTrackQAHistos(0x0),
74 fSecRecMCTrackQAHistos(0x0),
75 fSecGenMCTrackQAHistos(0x0),
84 // Default constructor
87 cout<<"AliAODTrackCutsDiHadronPID Default Constructor Called."<<endl;
88 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
90 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
92 // Initialize the data histograms.
94 fHistDataPt[iHistoName] = 0x0;
95 fHistDataPhiEtaPt[iHistoName] = 0x0;
96 fHistDataNTracks[iHistoName] = 0x0;
99 fHistDataDCAxy[iHistoName] = 0x0;
100 fHistDataDCAz[iHistoName] = 0x0;
102 // Initialize PID histograms.
103 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
104 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
105 fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
106 fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
107 fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
112 fNTracks[iHistoName] = 0;
114 // Initialize data DCA for one sigma.
115 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
117 // Initialize MC histograms.
118 fHistPrimGenMCPt[iHistoName] = 0x0;
119 fHistPrimRecMCPt[iHistoName] = 0x0;
120 fHistPrimRecNTracks[iHistoName] = 0x0;
122 fHistSecGenMCPt[iHistoName] = 0x0;
123 fHistSecRecMCPt[iHistoName] = 0x0;
125 // Initialze MC DCA histograms
126 fHistPrimRecMCDCA[iHistoName] = 0x0;
127 fHistSecRecMCDCAMat[iHistoName] = 0x0;
128 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
130 // Initialize other stuff.
131 fHistRequested[iHistoName] = kFALSE;
132 f3DSpectraEnabeled[iHistoName] = kFALSE;
133 fPIDHistosEnabeled[iHistoName] = kFALSE;
136 // TODO: It would be a bit neater to initialize all values to zero
137 // instead of the default values...
138 InitializeDefaultHistoNamesAndAxes();
142 // -----------------------------------------------------------------------
143 AliAODTrackCutsDiHadronPID::AliAODTrackCutsDiHadronPID(const char* name):
144 TNamed(name,"AOD Track Cuts"),
151 fMinSPDHitsForPtDeptDCAcut(0),
152 fPtDeptDCAxyCutFormula(0x0),
156 fTestFilterMask(kFALSE),
158 fTestMaxRapidity(kFALSE),
160 fTestTOFmismatch(kFALSE),
161 fTestPtDeptDCAcut(kFALSE),
162 fDataTrackQAHistos(0x0),
163 fPrimRecMCTrackQAHistos(0x0),
164 fPrimGenMCTrackQAHistos(0x0),
165 fSecRecMCTrackQAHistos(0x0),
166 fSecGenMCTrackQAHistos(0x0),
178 cout<<"AliAODTrackCutsDiHadronPID Named Constructor Called."<<endl;
179 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
181 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
183 // Initialize the data histograms.
184 if (iHistoName < 3) {
185 fHistDataPt[iHistoName] = 0x0;
186 fHistDataPhiEtaPt[iHistoName] = 0x0;
187 fHistDataNTracks[iHistoName] = 0x0;
190 fHistDataDCAxy[iHistoName] = 0x0;
191 fHistDataDCAz[iHistoName] = 0x0;
193 // Initialize PID histograms.
194 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
195 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
196 fHistDataPID[iHistoName][iSpecies][iPtClass] = 0x0;
197 fHistTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
198 fHistTPCTOFMismatch[iHistoName][iSpecies][iPtClass] = 0x0;
204 fNTracks[iHistoName] = 0;
206 // Initialize data DCA for one sigma.
207 fHistDataDCAxyOneSigma[iHistoName] = 0x0;
209 // Initialize MC histograms.
210 fHistPrimGenMCPt[iHistoName] = 0x0;
211 fHistPrimRecMCPt[iHistoName] = 0x0;
212 fHistPrimRecNTracks[iHistoName] = 0x0;
214 fHistSecGenMCPt[iHistoName] = 0x0;
215 fHistSecRecMCPt[iHistoName] = 0x0;
217 // Initialze MC DCA histograms
218 fHistPrimRecMCDCA[iHistoName] = 0x0;
219 fHistSecRecMCDCAMat[iHistoName] = 0x0;
220 fHistSecRecMCDCAWeak[iHistoName] = 0x0;
222 // Initialize other stuff.
223 fHistRequested[iHistoName] = kFALSE;
224 f3DSpectraEnabeled[iHistoName] = kFALSE;
225 fPIDHistosEnabeled[iHistoName] = kFALSE;
228 InitializeDefaultHistoNamesAndAxes();
232 // -----------------------------------------------------------------------
233 void AliAODTrackCutsDiHadronPID::InitializeDefaultHistoNamesAndAxes() {
235 // Initializes the histogram name conventions and the ranges of all the histograms
236 // by their default values. This method should only be called by the
237 // (default) constructor.
239 // TODO: User should be able to change these standard values at initialization.
240 // TODO: User should be able to retrieve all these variables with appropriate Getters.
242 cout<<"AliAODTrackCutsDiHadronPID - Initializing Default Histogram Names and axes..."<<endl;
243 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
245 // Setting the Pt axis for all histograms except the PID and Mismatch histograms.
246 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,
247 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,
248 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};
249 for (Int_t iPtBins = 0; iPtBins < 57; iPtBins++) {fPtAxis[iPtBins] = ptaxis[iPtBins];}
252 // Setting the pt range of the five PID histogram pt classes.
253 Double_t ptboundarypid[6] = {0.5,0.7,1.0,1.7,3.0,5.0};
254 for (Int_t iPtBoundary = 0; iPtBoundary < 6; iPtBoundary++) {fPtBoundaryPID[iPtBoundary] = ptboundarypid[iPtBoundary];}
256 // Setting the number of bins in pt for the five PID histogram pt classes.
257 Int_t nptbinspid[5] = {8,6,14,13,10};
258 for (Int_t iPtBins = 0; iPtBins < 5; iPtBins++) {fNPtBinsPID[iPtBins] = nptbinspid[iPtBins];}
260 // Setting the TOF axes for the PID histograms.
261 Double_t toflowerbound[5][3] = {{-2000.,-6000.,-10000.},{-2000.,-4000.,-10000.},{-1000.,-2000.,-5000.},{-1000.,-1000.,-2500.},{-500.,-500.,-1000.}};
262 Double_t tofupperbound[5][3] = {{10000.,10000.,6000.},{10000.,8000.,6000.},{6000.,6000.,6000.},{6000.,6000.,6000.},{4000.,4000.,6000.}};
263 Int_t tofbins[5][3] = {{120,160,160},{120,120,160},{140,140,165},{140,140,170},{90,90,140}};
264 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
265 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
266 fTOFLowerBound[iPtClass][iSpecies] = toflowerbound[iPtClass][iSpecies];
267 fTOFUpperBound[iPtClass][iSpecies] = tofupperbound[iPtClass][iSpecies];
268 fTOFbins[iPtClass][iSpecies] = tofbins[iPtClass][iSpecies];
272 // Setting the TPC axes for the PID histograms.
273 Double_t tpclowerbound[5][3] = {{-20.,-50.,-100.},{-20.,-30.,-80.},{-25.,-25.,-45.},{-25.,-25.,-45.},{-25.,-20.,-20.}};
274 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
275 Int_t tpcbins[5][3] = {{80,80,120},{80,70,100},{75,75,70},{70,70,70},{50,50,50}};
276 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
277 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
278 fTPCLowerBound[iPtClass][iSpecies] = tpclowerbound[iPtClass][iSpecies];
279 fTPCUpperBound[iPtClass][iSpecies] = tpcupperbound[iPtClass][iSpecies];
280 fTPCbins[iPtClass][iSpecies] = tpcbins[iPtClass][iSpecies];
284 // Names for the 12 (species,charge) combinations considered.
285 const char* histoname[12] = {"AllCharged","Pos","Neg","AllPion","PosPion","NegPion","AllKaon","PosKaon","NegKaon","AllProton","PosProton","NegProton"};
286 const char* histolatex[12] = {"ch","ch^{+}","ch^{-}","#pi^{+/-}","#pi^{+}","#pi^{-}","K^{+/-}","K^{+}","K^{-}","p/#bar{p}","p","#bar{p}"};
287 for (Int_t iHistoName = 0; iHistoName < 12; iHistoName++) {
288 fHistoName[iHistoName] = histoname[iHistoName];
289 fHistoLatex[iHistoName] = histolatex[iHistoName];
292 // Names of the 3 species considered.
293 const char* particlename[3] = {"Pion","Kaon","Proton"};
294 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {fParticleName[iSpecies] = particlename[iSpecies];}
296 // Names of the 5 pt classes considered for the PID histograms.
297 const char* ptclassname[5] = {"VeryLowPt","LowPt","MedPt","MedHighPt","HighPt"};
298 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {fPtClassName[iPtClass] = ptclassname[iPtClass];}
302 // -----------------------------------------------------------------------
303 AliAODTrackCutsDiHadronPID::~AliAODTrackCutsDiHadronPID() {
309 cout<<"AliAODTrackCutsDiHadronPID Destructor Called."<<endl;
310 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
312 if (fDataTrackQAHistos) delete fDataTrackQAHistos;
313 fDataTrackQAHistos = 0x0;
314 if (fPrimRecMCTrackQAHistos) delete fPrimRecMCTrackQAHistos;
315 fPrimRecMCTrackQAHistos = 0x0;
316 if (fPrimGenMCTrackQAHistos) delete fPrimGenMCTrackQAHistos;
317 fPrimGenMCTrackQAHistos = 0x0;
318 if (fSecRecMCTrackQAHistos) delete fSecRecMCTrackQAHistos;
319 fSecRecMCTrackQAHistos = 0x0;
320 if (fSecGenMCTrackQAHistos) delete fSecGenMCTrackQAHistos;
321 fSecGenMCTrackQAHistos = 0x0;
325 // -----------------------------------------------------------------------
326 Long64_t AliAODTrackCutsDiHadronPID::Merge(TCollection* list) {
332 cout<<"AliAODTrackCutsDiHadronPID Merger Called."<<endl;
333 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
335 Bool_t HistosOK = kTRUE;
338 if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos||!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {HistosOK = kFALSE;}
340 if (!fDataTrackQAHistos) {HistosOK = kFALSE;}
344 cout<<"AliAODTrackCutsDiHadronPID::Merge() - Warning, current object's histograms are missing... Generating."<<endl;
349 //cout<<"No list available..."<<endl;
352 if (list->IsEmpty()) {
353 //cout<<"List is empty..."<<endl;
357 //Int_t NEntries = list->GetEntries();
358 //cout<<"Supplied list has "<<NEntries<<" entries."<<endl;
360 TIterator* iter = list->MakeIterator();
363 // List of collections
364 TList collection_fDataTrackQAHistos;
365 TList collection_fPrimRecMCTrackQAHistos;
366 TList collection_fPrimGenMCTrackQAHistos;
367 TList collection_fSecRecMCTrackQAHistos;
368 TList collection_fSecGenMCTrackQAHistos;
372 while ((obj = iter->Next())) {
373 AliAODTrackCutsDiHadronPID* entry = dynamic_cast<AliAODTrackCutsDiHadronPID*> (obj);
374 if (entry == 0) continue;
376 // Check if the object to be merged really has the same name! (FIXME!)
378 // Getting the lists from obj.
379 TList* list_fDataTrackQAHistos = entry->GetListOfDataQAHistos();
380 TList* list_fPrimRecMCTrackQAHistos = entry->GetListOfPrimRecMCTrackQAHistos();
381 TList* list_fPrimGenMCTrackQAHistos = entry->GetListOfPrimGenMCTrackQAHistos();
382 TList* list_fSecRecMCTrackQAHistos = entry->GetListOfSecRecMCTrackQAHistos();
383 TList* list_fSecGenMCTrackQAHistos = entry->GetListOfSecGenMCTrackQAHistos();
385 // Adding the retrieved lists to the collection.
386 if (list_fDataTrackQAHistos) collection_fDataTrackQAHistos.Add(list_fDataTrackQAHistos);
387 if (list_fPrimRecMCTrackQAHistos) collection_fPrimRecMCTrackQAHistos.Add(list_fPrimRecMCTrackQAHistos);
388 if (list_fPrimGenMCTrackQAHistos) collection_fPrimGenMCTrackQAHistos.Add(list_fPrimGenMCTrackQAHistos);
389 if (list_fSecRecMCTrackQAHistos) collection_fSecRecMCTrackQAHistos.Add(list_fSecRecMCTrackQAHistos);
390 if (list_fSecGenMCTrackQAHistos) collection_fSecGenMCTrackQAHistos.Add(list_fSecGenMCTrackQAHistos);
392 //cout<<"Entries intermediate list Data: "<<collection_fDataTrackQAHistos.GetEntries()<<endl;
393 //cout<<"Entries intermediate list RecMC: "<<collection_fPrimRecMCTrackQAHistos.GetEntries()<<endl;
398 // Merging. Note that we require the original list to exist.
399 // * Assume that if the collection happens to be empty, then nothing will happen.
400 // |-> This of course leads to some problems, since if the first file does not
401 // | have such a list, then the merged file will have no results...
402 // | IDEA: if fDataTrackQAHistos does not exist, then create a new one. (TO DO)
403 // * All other variables are taken from the original object.
404 if (fDataTrackQAHistos) fDataTrackQAHistos->Merge(&collection_fDataTrackQAHistos);
405 if (fPrimRecMCTrackQAHistos) fPrimRecMCTrackQAHistos->Merge(&collection_fPrimRecMCTrackQAHistos);
406 if (fPrimGenMCTrackQAHistos) fPrimGenMCTrackQAHistos->Merge(&collection_fPrimGenMCTrackQAHistos);
407 if (fSecRecMCTrackQAHistos) fSecRecMCTrackQAHistos->Merge(&collection_fSecRecMCTrackQAHistos);
408 if (fSecGenMCTrackQAHistos) fSecGenMCTrackQAHistos->Merge(&collection_fSecGenMCTrackQAHistos);
416 // -----------------------------------------------------------------------
417 void AliAODTrackCutsDiHadronPID::PrintCuts() const { /* NOT IMPLEMENTED */ }
419 // -----------------------------------------------------------------------
420 TList* AliAODTrackCutsDiHadronPID::GetListOfDataQAHistos() const {
422 // Returns the list of data histograms.
423 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
425 if (fDataTrackQAHistos) {
426 return fDataTrackQAHistos;
433 // -----------------------------------------------------------------------
434 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimRecMCTrackQAHistos() const {
436 // Returns the list of histograms of reconstructed primary tracks.
437 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
439 if (fPrimRecMCTrackQAHistos) {
440 return fPrimRecMCTrackQAHistos;
447 // -----------------------------------------------------------------------
448 TList* AliAODTrackCutsDiHadronPID::GetListOfPrimGenMCTrackQAHistos() const {
450 // Returns the list of histograms of generator level primary particles.
451 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
453 if (fPrimGenMCTrackQAHistos) {
454 return fPrimGenMCTrackQAHistos;
461 // -----------------------------------------------------------------------
462 TList* AliAODTrackCutsDiHadronPID::GetListOfSecRecMCTrackQAHistos() const {
464 // Returns the list of histograms of reconstructed secondary tracks.
465 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
467 if (fSecRecMCTrackQAHistos) {
468 return fSecRecMCTrackQAHistos;
475 // -----------------------------------------------------------------------
476 TList* AliAODTrackCutsDiHadronPID::GetListOfSecGenMCTrackQAHistos() const {
478 // Returns the list of histograms of generator level secondary particles.
479 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
481 if (fSecGenMCTrackQAHistos) {
482 return fSecGenMCTrackQAHistos;
489 // -----------------------------------------------------------------------
490 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
492 // Returns a projection in TOF of the data histogram.
493 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
495 Int_t ptclass = GetPtClass(ptbin);
496 if (ptclass == -1) return 0x0;
497 Int_t bininptclass = GetBinInPtClass(ptbin);
499 // Retrieve original the 3D histogram (we don't own it).
500 TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
502 // Make the projection.
503 TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojection_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
505 // Some settings of the output histogram.
506 htmp_proj->SetDirectory(0);
507 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
514 // -----------------------------------------------------------------------
515 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFProjection(Int_t charge, Int_t species) {
517 // Returns a TObjArray with all TOF histograms (as needed by AliSpectrumDiHadronPID)
518 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
520 TObjArray* aout = new TObjArray(GetNPtBinsPID());
521 aout->SetOwner(kTRUE);
522 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
523 aout->AddLast((TH1F*)GetHistDataTOFProjection(charge, species, iPtBin));
530 // -----------------------------------------------------------------------
531 TH1F* AliAODTrackCutsDiHadronPID::GetHistDataTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
533 // Returns a projection in TOF of the mismatch histogram.
534 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
536 Int_t ptclass = GetPtClass(ptbin);
537 if (ptclass == -1) return 0x0;
538 Int_t bininptclass = GetBinInPtClass(ptbin);
540 // Retrieve original the 2D histogram (we don't own it).
541 TH2F* htmp = (TH2F*)GetHistData(Form("fHistTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
543 // Make the projection.
544 TH1F* htmp_proj = (TH1F*)htmp->ProjectionY(Form("TOFprojectionMismatch_%i_%i_%i",charge,species,ptbin),bininptclass,bininptclass);
546 // Some settings of the output histogram.
547 htmp_proj->SetDirectory(0);
548 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
555 // -----------------------------------------------------------------------
556 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTOFMismatch(Int_t charge, Int_t species) {
558 // Returns a TObjArray with all TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
559 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
561 TObjArray* aout = new TObjArray(GetNPtBinsPID());
562 aout->SetOwner(kTRUE);
563 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
564 aout->AddLast((TH1F*)GetHistDataTOFMismatch(charge, species, iPtBin));
571 // -----------------------------------------------------------------------
572 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFProjection(Int_t charge, Int_t species, Int_t ptbin) {
574 // Returns a projection in TOF of the data histogram.
575 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
577 Int_t ptclass = GetPtClass(ptbin);
578 if (ptclass == -1) return 0x0;
579 Int_t bininptclass = GetBinInPtClass(ptbin);
581 // Retrieve original the 3D histogram (we don't own it).
582 TH3F* htmp = (TH3F*)GetHistData(Form("fHistDataPID%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
584 // Make the projection.
585 TAxis* ptaxis = htmp->GetXaxis();
586 Int_t NbinsPt = ptaxis->GetNbins();
587 //cout<<"bin in pt class: "<<bininptclass<<endl;
588 ptaxis->SetRange(bininptclass, bininptclass);
589 TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
590 htmp_proj->SetName(Form("TPCTOFprojection_%i_%i_%i",charge,species,ptbin));
592 //cout<<"projection: "<<htmp_proj<<endl;
594 // Some settings of the output histogram.
595 htmp_proj->SetDirectory(0);
596 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
599 // Putting back the range on the p_T axis.
600 ptaxis->SetRange(1, NbinsPt);
606 // -----------------------------------------------------------------------
607 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFProjection(Int_t charge, Int_t species) {
609 // TODO: This is basically a copy of GetDataTOFProjection -> Make into one function.
611 // Returns a TObjArray with all TPC/TOF histograms (as needed by AliSpectrumDiHadronPID)
612 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
614 TObjArray* aout = new TObjArray(GetNPtBinsPID());
615 aout->SetOwner(kTRUE);
616 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
618 aout->AddLast((TH2F*)GetHistDataTPCTOFProjection(charge, species, iPtBin));
619 //cout<<"aout entries: "<<aout->GetEntriesFast()<<endl;
626 // -----------------------------------------------------------------------
627 TH2F* AliAODTrackCutsDiHadronPID::GetHistDataTPCTOFMismatch(Int_t charge, Int_t species, Int_t ptbin) {
629 // Returns a projection in TPC/TOF of the mismatch histogram.
630 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
632 Int_t ptclass = GetPtClass(ptbin);
633 if (ptclass == -1) return 0x0;
634 Int_t bininptclass = GetBinInPtClass(ptbin);
636 // Retrieve original the 3D histogram (we don't own it).
637 TH3F* htmp = (TH3F*)GetHistData(Form("fHistTPCTOFMismatch%s%s%s",fHistoName[charge].Data(),fParticleName[species].Data(),fPtClassName[ptclass].Data()));
639 // Make the projection.
640 TAxis* ptaxis = htmp->GetXaxis();
641 Int_t NbinsPt = ptaxis->GetNbins();
642 ptaxis->SetRange(bininptclass, bininptclass);
643 TH2F* htmp_proj = (TH2F*)htmp->Project3D("zy");
644 htmp_proj->SetName(Form("TPCTOFprojectionMismatch_%i_%i_%i",charge,species,ptbin));
646 // Some settings of the output histogram.
647 htmp_proj->SetDirectory(0);
648 htmp_proj->SetTitle(Form("%5.3f < p_{T} < %5.3f",GetPtMinPID(ptbin),GetPtMaxPID(ptbin)));
651 // Putting back the range on the p_T axis.
652 ptaxis->SetRange(1, NbinsPt);
658 // -----------------------------------------------------------------------
659 TObjArray* AliAODTrackCutsDiHadronPID::GetDataTPCTOFMismatch(Int_t charge, Int_t species) {
661 // Returns a TObjArray with all TPC/TOF mismatch histograms (as needed by AliSpectrumDiHadronPID)
662 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
664 TObjArray* aout = new TObjArray(GetNPtBinsPID());
665 aout->SetOwner(kTRUE);
666 for (Int_t iPtBin = 1; iPtBin < (GetNPtBinsPID() + 1); iPtBin++) {
667 aout->AddLast((TH2F*)GetHistDataTPCTOFMismatch(charge, species, iPtBin));
674 // -----------------------------------------------------------------------
675 Double_t AliAODTrackCutsDiHadronPID::GetPtMin(const Int_t bin) const {
677 // Same as: TAxis::GetBinLowEdge()
679 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
681 if ((bin < 1) || (bin > fNPtBins + 1)) {
682 cout<<"Bin is out of range..."<<endl;
685 return fPtAxis[bin - 1];
690 // -----------------------------------------------------------------------
691 Double_t AliAODTrackCutsDiHadronPID::GetPtMax(const Int_t bin) const {
693 // Same as: TAxis::GetBinUpEdge()
695 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
697 if ((bin < 1) || (bin > fNPtBins + 1)) {
698 cout<<"Bin is out of range..."<<endl;
706 // -----------------------------------------------------------------------
707 Int_t AliAODTrackCutsDiHadronPID::GetNPtBinsPID(const Int_t ptclass) const {
709 // Returns the number of pt bins that are used in every "pt class",
710 // where a "pt class" is a range in pT which have the same range in
711 // TOF/ TPC. If class == -1, then the retrun value is the total number
712 // of bins in all classes together.
714 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
717 Int_t nptbinspid = 0;
718 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
719 nptbinspid += fNPtBinsPID[iPtClass];
722 } else if (ptclass >= 0 && ptclass < 5) {
723 return fNPtBinsPID[ptclass];
730 // -----------------------------------------------------------------------
731 Double_t* AliAODTrackCutsDiHadronPID::GetPtAxisPID() const {
733 // Returns an array representing the pT axis for the PID
736 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
738 const Int_t nptbinspid = GetNPtBinsPID();
739 Double_t* ptaxis = new Double_t[nptbinspid + 1];
741 for (Int_t iPtBin = 0; iPtBin < nptbinspid; iPtBin++) {
742 ptaxis[iPtBin] = GetPtMinPID(iPtBin + 1);
745 ptaxis[nptbinspid] = GetPtMaxPID(nptbinspid);
750 // -----------------------------------------------------------------------
751 Double_t AliAODTrackCutsDiHadronPID::GetPtMinPID(const Int_t bin) const {
753 // Same as: TAxis::GetBinLowEdge()
755 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
757 Int_t ptclass = GetPtClass(bin);
759 if (ptclass == -1) {return -999.;}
761 Int_t bininptclass = GetBinInPtClass(bin);
763 Double_t minpt = fPtBoundaryPID[ptclass];
764 Double_t maxpt = fPtBoundaryPID[ptclass+1];
765 Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
767 return (minpt + ptres * ((Double_t)(bininptclass - 1)) );
771 // -----------------------------------------------------------------------
772 Double_t AliAODTrackCutsDiHadronPID::GetPtMaxPID(const Int_t bin) const {
774 // Same as: TAxis::GetBinUpEdge()
776 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
778 Int_t ptclass = GetPtClass(bin);
780 if (ptclass == -1) {return -999.;}
782 Int_t bininptclass = GetBinInPtClass(bin);
784 Double_t minpt = fPtBoundaryPID[ptclass];
785 Double_t maxpt = fPtBoundaryPID[ptclass+1];
786 Double_t ptres = (maxpt - minpt)/((Double_t)fNPtBinsPID[ptclass]);
788 return (minpt + ptres * ((Double_t)(bininptclass)) );
792 // -----------------------------------------------------------------------
793 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMin(const Int_t ptclass) const {
795 // Returns the minimum p_T of a p_T class.
797 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
799 if (ptclass >= 0 && ptclass < 5) {
800 return fPtBoundaryPID[ptclass];
807 // -----------------------------------------------------------------------
808 Double_t AliAODTrackCutsDiHadronPID::GetPtClassMax(const Int_t ptclass) const {
810 // Returns the maximum p_T of a p_T class.
812 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
814 if (ptclass >= 0 && ptclass < 5) {
815 return fPtBoundaryPID[ptclass+1];
822 // -----------------------------------------------------------------------
823 Bool_t AliAODTrackCutsDiHadronPID::RequestQAHistos(const Int_t histoclass, const Bool_t Enable3DSpectra, const Bool_t EnablePIDHistos) {
825 // Request certain histograms to be filled.
827 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
829 if ((histoclass > -1) && (histoclass < 12)) {
830 fHistRequested[histoclass] = kTRUE;
831 f3DSpectraEnabeled[histoclass] = Enable3DSpectra;
832 fPIDHistosEnabeled[histoclass] = EnablePIDHistos;
833 //cout<<"histoclass: "<<histoclass<<" requested: "<<fHistRequested[histoclass]<<endl;
840 // -----------------------------------------------------------------------
841 void AliAODTrackCutsDiHadronPID::StartNewEvent() {
843 // FIXME: This method is now only suited for Data (3 histo classes only.)
844 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
846 // Resetting the counters.
847 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
848 fNTracks[iHistoClass] = 0;
852 // -----------------------------------------------------------------------
853 void AliAODTrackCutsDiHadronPID::EventIsDone(Bool_t IsMC) {
855 // FIXME: This method is now only suited for Data (3 histo classes only.)
856 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
858 // Fill NTracks histos.
860 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
862 if (fHistRequested[iHistoClass]) {
865 // THE FOLLOWING SHOULD NEVER HAPPEN.
866 //if (!fHistPrimRecMCPt[iHistoClass]) InitializeRecMCHistos(iHistoClass);
867 fHistPrimRecNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
869 // THE FOLLOWING SHOULD NEVER HAPPEN.
870 //if (!fHistDataPt[iHistoClass]) InitializeDataHistos(iHistoClass);
871 fHistDataNTracks[iHistoClass]->Fill(fNTracks[iHistoClass]);
877 // -----------------------------------------------------------------------
878 void AliAODTrackCutsDiHadronPID::CreateHistos() {
880 // Function should be called by the UserCreateOutput() function of the
881 // analysis task. This function then generates all the histograms that
882 // were requested locally on the workernode. Even if case that the
883 // histograms are not filled, it is imperative that they are still
884 // created, and in the same order, otherwise problems with merging will
887 // TODO: In principle it should never happen that if this method is called,
888 // that the lists of QA histograms already exist, so this may become a fatal error
889 // instead of a warning.
891 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
895 if (!fPrimGenMCTrackQAHistos) {
896 cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Gen. MC Track QA TList..."<<endl;
897 fPrimGenMCTrackQAHistos = new TList();
898 fPrimGenMCTrackQAHistos->SetName("PrimGenMCTrackQAHistos");
899 fPrimGenMCTrackQAHistos->SetOwner(kTRUE);
900 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Gen. MC Track QA TList was already created..."<<endl;}
902 if (!fSecGenMCTrackQAHistos) {
903 cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Gen. MC Track QA TList..."<<endl;
904 fSecGenMCTrackQAHistos = new TList();
905 fSecGenMCTrackQAHistos->SetName("SecGenMCTrackQAHistos");
906 fSecGenMCTrackQAHistos->SetOwner(kTRUE);
907 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {if (fHistRequested[iHistoClass]) {InitializeGenMCHistos(iHistoClass);}}
908 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Gen. MC Track QA TList was already created..."<<endl;}
910 if (!fPrimRecMCTrackQAHistos) {
911 cout<<"AliAODTrackCutsDiHadronPID - Creating Prim. Rec. MC Track QA TList..."<<endl;
912 fPrimRecMCTrackQAHistos = new TList();
913 fPrimRecMCTrackQAHistos->SetName("PrimRecMCTrackQAHistos");
914 fPrimRecMCTrackQAHistos->SetOwner(kTRUE);
915 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Prim. Rec. MC Track QA TList was already created..."<<endl;}
917 if (!fSecRecMCTrackQAHistos) {
918 cout<<"AliAODTrackCutsDiHadronPID - Creating Sec. Rec. MC Track QA TList..."<<endl;
919 fSecRecMCTrackQAHistos = new TList();
920 fSecRecMCTrackQAHistos->SetName("SecRecMCTrackQAHistos");
921 fSecRecMCTrackQAHistos->SetOwner(kTRUE);
922 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
923 if (fHistRequested[iHistoClass]) {
924 InitializeRecMCHistos(iHistoClass);
927 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Sec. Rec. MC Track QA TList was already created..."<<endl;}
931 if (!fDataTrackQAHistos) {
932 cout<<"AliAODTrackCutsDiHadronPID - Creating Data Track QA TList..."<<endl;
933 fDataTrackQAHistos = new TList();
934 fDataTrackQAHistos->SetName("DataTrackQAHistos");
935 fDataTrackQAHistos->SetOwner(kTRUE);
936 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
937 fHistDataDCAxyOneSigma[iHistoClass] = InitializeDCASpectrum("fHistDataDCAxyOneSigma",iHistoClass);
938 fDataTrackQAHistos->Add(fHistDataDCAxyOneSigma[iHistoClass]);
939 if (fHistRequested[iHistoClass]) {
940 if (iHistoClass < 3) {InitializeDataHistos(iHistoClass);}
943 } else {cout<<"AliAODTrackCutsDiHadronPID - Warning, Data QA TList was already created..."<<endl;}
948 // -----------------------------------------------------------------------
949 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedData(AliTrackDiHadronPID* track, Double_t randomhittime) {
952 // Checks performed on a data track.
955 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
957 if (!track) return kFALSE;
959 if (!fDataTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
961 // Check the track cuts (NB. check function return kTRUE if track is accepted.)
962 if (!CheckPt(track->Pt())) return kFALSE;
963 if (!CheckMaxEta(track->Eta())) return kFALSE;
964 if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
965 if (!CheckFlags(track->GetFlags())) return kFALSE;
966 if (!CheckTOFmismatch(track->IsTOFmismatch())) return kFALSE;
969 if (track->HasPointOnITSLayer(0)) NSPDhits++;
970 if (track->HasPointOnITSLayer(1)) NSPDhits++;
971 if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
973 // Track has passed the cuts, fill QA histograms.
974 for (Int_t iHistoClass = 0; iHistoClass < 3; iHistoClass++) {
976 if (fHistRequested[iHistoClass]) {
978 // Check the charge (could be neater).
979 if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
980 if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
981 if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
983 FillDataHistos(iHistoClass, track);
985 // Ignore if random hit is < -1.e20.
986 if (randomhittime > -1.e20) FillTOFMismatchHistos(iHistoClass,track,randomhittime);
988 fNTracks[iHistoClass]++;
998 // -----------------------------------------------------------------------
999 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedGeneratedMC(AliAODMCParticle* particle) {
1002 // Checks performed on a generated MC particle.
1005 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1007 // PDG codes for particles:
1008 // pi+ 211, K+ 321, p 2212
1010 if (!particle) return kFALSE;
1012 if (!fPrimGenMCTrackQAHistos||!fSecGenMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1014 // Check the track cuts.
1015 if (!CheckPt(particle->Pt())) return kFALSE;
1016 if (!CheckMaxEta(particle->Eta())) return kFALSE;
1017 if (!CheckRapidity(particle->Y())) return kFALSE; // NEW: rapidity cut. (not done on data)
1019 // NOT YET IMPLEMENTED! - Histoclasses for specific particles.
1020 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1022 if (fHistRequested[iHistoClass]) {
1024 // Check the charge (could be neater).
1025 if ((iHistoClass == 0) && (particle->Charge() == 0)) continue;
1026 if ((iHistoClass == 1) && (particle->Charge() <= 0)) continue;
1027 if ((iHistoClass == 2) && (particle->Charge() >= 0)) continue;
1028 if ((iHistoClass == 3) && (TMath::Abs(particle->GetPdgCode()) != 211)) continue;
1029 if ((iHistoClass == 4) && (particle->GetPdgCode()) != 211) continue;
1030 if ((iHistoClass == 5) && (particle->GetPdgCode()) != -211) continue;
1031 if ((iHistoClass == 6) && (TMath::Abs(particle->GetPdgCode()) != 321)) continue;
1032 if ((iHistoClass == 7) && (particle->GetPdgCode()) != 321) continue;
1033 if ((iHistoClass == 8) && (particle->GetPdgCode()) != -321) continue;
1034 if ((iHistoClass == 9) && (TMath::Abs(particle->GetPdgCode()) != 2212)) continue;
1035 if ((iHistoClass == 10) && (particle->GetPdgCode()) != 2212) continue;
1036 if ((iHistoClass == 11) && (particle->GetPdgCode()) != -2212) continue;
1038 // Secondary specification not set.
1039 //cout<<"Charge: "<<particle->Charge()<<" PDG code: "<<particle->GetPdgCode()<<endl;
1040 //cout<<particle->IsPhysicalPrimary()<<" "<<particle->IsSecondaryFromWeakDecay()<<" "<<particle->IsSecondaryFromMaterial()<<endl;
1042 // These two functions are not implemented...
1043 if (particle->IsSecondaryFromWeakDecay()) cout<<"Secondary From Weak Decay!"<<endl;
1044 if (particle->IsSecondaryFromMaterial()) cout<<"Secondary From Material!"<<endl;
1046 FillGenMCHistos(iHistoClass, particle);
1051 // Particle has Passed.
1056 // -----------------------------------------------------------------------
1057 Bool_t AliAODTrackCutsDiHadronPID::IsSelectedReconstructedMC(AliTrackDiHadronPID* track) {
1059 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1061 if (!track) return kFALSE;
1063 if (!fPrimRecMCTrackQAHistos||!fSecRecMCTrackQAHistos) {cout<<"AliAODTrackCutsDiHadronPID - Histograms were not created, you should have called CreateHistos()..."<<endl;}
1065 // Check the track cuts.
1066 if (!CheckPt(track->MCPt())) return kFALSE; // Kinematic cuts are done on the MC particles
1067 if (!CheckMaxEta(track->MCEta())) return kFALSE; // and not on the data.
1068 if (!CheckRapidity(track->MCY())) return kFALSE; // NEW: rapidity cut.
1069 if (!CheckFilterMask(track->GetFilterMap())) return kFALSE;
1070 if (!CheckFlags(track->GetFlags())) return kFALSE;
1073 if (track->HasPointOnITSLayer(0)) NSPDhits++;
1074 if (track->HasPointOnITSLayer(1)) NSPDhits++;
1075 if (!CheckPtDeptDCACut(track->GetZAtDCA(),track->GetXYAtDCA(),track->Pt(),NSPDhits)) return kFALSE;
1077 // Track has passed the cuts, fill QA histograms.
1078 for (Int_t iHistoClass = 0; iHistoClass < 12; iHistoClass++) {
1080 if (fHistRequested[iHistoClass]) {
1082 // Check the charge (could be neater).
1083 if ((iHistoClass == 0) && (track->Charge() == 0)) continue;
1084 if ((iHistoClass == 1) && (track->Charge() <= 0)) continue;
1085 if ((iHistoClass == 2) && (track->Charge() >= 0)) continue;
1086 if ((iHistoClass == 3) && (TMath::Abs(track->GetPdgCode()) != 211)) continue;
1087 if ((iHistoClass == 4) && (track->GetPdgCode()) != 211) continue;
1088 if ((iHistoClass == 5) && (track->GetPdgCode()) != -211) continue;
1089 if ((iHistoClass == 6) && (TMath::Abs(track->GetPdgCode()) != 321)) continue;
1090 if ((iHistoClass == 7) && (track->GetPdgCode()) != 321) continue;
1091 if ((iHistoClass == 8) && (track->GetPdgCode()) != -321) continue;
1092 if ((iHistoClass == 9) && (TMath::Abs(track->GetPdgCode()) != 2212)) continue;
1093 if ((iHistoClass == 10) && (track->GetPdgCode()) != 2212) continue;
1094 if ((iHistoClass == 11) && (track->GetPdgCode()) != -2212) continue;
1096 FillRecMCHistos(iHistoClass, track);
1098 fNTracks[iHistoClass]++;
1102 // Track Has Passed.
1107 // -----------------------------------------------------------------------
1108 Int_t AliAODTrackCutsDiHadronPID::GetPtClass(const Int_t ptbin) const {
1110 // Returns a p_T class as a function of bin (PID histos)
1111 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1113 Int_t currentptclass = 0;
1114 Int_t currentptbin = fNPtBinsPID[0];
1116 while (currentptbin < ptbin) {
1118 if (currentptclass == 5) {break;}
1119 currentptbin += fNPtBinsPID[currentptclass];
1122 if (currentptclass == 5) {
1123 cout<<"GetPtClass -> ptbin out of range!"<<endl;
1127 return currentptclass;
1130 // -----------------------------------------------------------------------
1131 Int_t AliAODTrackCutsDiHadronPID::GetBinInPtClass(const Int_t ptbin) const {
1133 // Returns the bin withing the p_T class (i.e., 1..Nbins)
1134 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1136 Int_t ptclass = GetPtClass(ptbin);
1137 if (ptclass == -1) {return -1;}
1139 Int_t ptbinout = ptbin;
1140 for (Int_t iPtClass = 0; iPtClass < ptclass; iPtClass++) {ptbinout -= fNPtBinsPID[iPtClass];}
1146 // -----------------------------------------------------------------------
1147 Bool_t AliAODTrackCutsDiHadronPID::FillDataHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1149 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1151 // Fill the histograms.
1152 fHistDataPt[histoclass]->Fill(track->Pt());
1153 fHistDataPhiEtaPt[histoclass]->Fill(track->Phi(),track->Eta(),track->Pt());
1155 // Fill DCA histograms.
1156 fHistDataDCAxy[histoclass]->Fill(track->GetXYAtDCA());
1157 fHistDataDCAz[histoclass]->Fill(track->GetZAtDCA());
1159 // Fill DCA_{xy} one sigma histograms.
1161 // histoclass = 0: all charges; histoclass = 1: positive; histoclass = 2: negative;
1162 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(0) * track->GetNumberOfSigmasTOF(0) +
1163 track->GetNumberOfSigmasTPC(0) * track->GetNumberOfSigmasTPC(0)) < 1.) {
1164 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1165 fHistDataDCAxyOneSigma[3 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Pions.
1166 // checkSum++; cout<<"Pion found: nSigTOF: "<<track->GetNumberOfSigmasTOF(0)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(0)<<endl;
1168 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(1) * track->GetNumberOfSigmasTOF(1) +
1169 track->GetNumberOfSigmasTPC(1) * track->GetNumberOfSigmasTPC(1)) < 1.) {
1170 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1171 fHistDataDCAxyOneSigma[6 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Kaons.
1172 // checkSum++; cout<<"Kaon found: nSigTOF: "<<track->GetNumberOfSigmasTOF(1)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(1)<<endl;
1174 if (TMath::Sqrt(track->GetNumberOfSigmasTOF(2) * track->GetNumberOfSigmasTOF(2) +
1175 track->GetNumberOfSigmasTPC(2) * track->GetNumberOfSigmasTPC(2)) < 1.) {
1176 fHistDataDCAxyOneSigma[0 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // All species.
1177 fHistDataDCAxyOneSigma[9 + histoclass]->Fill(track->Pt(),track->GetXYAtDCA()); // Protons.
1178 // checkSum++; cout<<"Proton found: nSigTOF: "<<track->GetNumberOfSigmasTOF(2)<<"; nSigTPC: "<<track->GetNumberOfSigmasTPC(2)<<endl;
1181 // check for double identification (or triple..)
1182 if(checkSum > 1) {AliError("More than one particle identified for the same track!"); }
1185 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1187 // Note that a possible Y cut is only done on the PID cuts!
1188 if (!CheckRapidity(track->Y(iSpecies))) continue;
1190 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1191 fHistDataPID[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),track->GetTOFsignalMinusExpected(iSpecies),track->GetTPCsignalMinusExpected(iSpecies));
1199 // -----------------------------------------------------------------------
1200 Bool_t AliAODTrackCutsDiHadronPID::FillTOFMismatchHistos(Int_t histoclass, AliTrackDiHadronPID* track, Double_t randomhittime) {
1202 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1204 // Fill TOF mismatch.
1205 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1207 // Note that a possible Y cut is only done on the PID cuts!
1208 if (!CheckRapidity(track->Y(iSpecies))) continue;
1210 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1211 fHistTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),randomhittime - track->GetTOFsignalExpected(iSpecies));
1212 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->Fill(track->Pt(),randomhittime - track->GetTOFsignalExpected(iSpecies),track->GetTPCsignalMinusExpected(iSpecies));
1220 // -----------------------------------------------------------------------
1221 Bool_t AliAODTrackCutsDiHadronPID::FillGenMCHistos(Int_t histoclass, AliAODMCParticle* particle) {
1223 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1225 // Fill the histograms.
1226 if (particle->IsPhysicalPrimary()) {
1227 fHistPrimGenMCPt[histoclass]->Fill(particle->Pt());
1229 fHistSecGenMCPt[histoclass]->Fill(particle->Pt());
1236 // -----------------------------------------------------------------------
1237 Bool_t AliAODTrackCutsDiHadronPID::FillRecMCHistos(Int_t histoclass, AliTrackDiHadronPID* track) {
1239 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1241 // Fill the Pt histograms.
1242 if (track->IsPhysicalPrimary()) {
1243 fHistPrimRecMCPt[histoclass]->Fill(track->Pt());
1245 fHistSecRecMCPt[histoclass]->Fill(track->Pt());
1248 // Fill the DCA histograms.
1249 if (track->IsPhysicalPrimary()) {
1250 fHistPrimRecMCDCA[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1252 if (track->IsSecondaryFromMaterial()) {
1253 fHistSecRecMCDCAMat[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1255 if (track->IsSecondaryFromWeakDecay()) {
1256 fHistSecRecMCDCAWeak[histoclass]->Fill(track->Pt(), track->GetXYAtDCA());
1263 // -----------------------------------------------------------------------
1264 Bool_t AliAODTrackCutsDiHadronPID::InitializeDataHistos(Int_t histoclass) {
1266 cout<<"AliAODTrackCutsDiHadronPID - Creating Data Histograms of Class "<<histoclass<<endl;
1267 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1269 // Add the Pt spectra.
1270 fHistDataPt[histoclass] = InitializePtSpectrum("fHistDataPt",histoclass);
1271 fDataTrackQAHistos->Add(fHistDataPt[histoclass]);
1272 fHistDataPhiEtaPt[histoclass] = InitializePhiEtaPt("fHistDataPhiEtaPt",histoclass);
1273 fDataTrackQAHistos->Add(fHistDataPhiEtaPt[histoclass]);
1275 // Add the NTrack histograms.
1276 fHistDataNTracks[histoclass] = InitializeNTracksHisto("fHistDataNTracks",histoclass);
1277 fDataTrackQAHistos->Add(fHistDataNTracks[histoclass]);
1279 // Add the DCA histograms.
1280 fHistDataDCAxy[histoclass] = InitializeDCAxyHisto("fHistDCAxy",histoclass);
1281 fDataTrackQAHistos->Add(fHistDataDCAxy[histoclass]);
1282 fHistDataDCAz[histoclass] = InitializeDCAzHisto("fHistDCAz",histoclass);
1283 fDataTrackQAHistos->Add(fHistDataDCAz[histoclass]);
1285 // Add the PID histograms. (FIXME shoudl be able to turn the creation of these histos off.)
1286 for (Int_t iSpecies = 0; iSpecies < 3; iSpecies++) {
1287 for (Int_t iPtClass = 0; iPtClass < 5; iPtClass++) {
1289 fHistDataPID[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistDataPID",histoclass,iSpecies,iPtClass);
1290 fDataTrackQAHistos->Add(fHistDataPID[histoclass][iSpecies][iPtClass]);
1292 fHistTOFMismatch[histoclass][iSpecies][iPtClass] = InitializeTOFMismatchHisto("fHistTOFMismatch",histoclass,iSpecies,iPtClass);
1293 fDataTrackQAHistos->Add(fHistTOFMismatch[histoclass][iSpecies][iPtClass]);
1295 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass] = InitializePIDHisto("fHistTPCTOFMismatch",histoclass,iSpecies,iPtClass);
1296 fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]->SetTitle(Form("PID %s (Exp: %s)",fHistoName[histoclass].Data(),fParticleName[iSpecies].Data()));
1297 fDataTrackQAHistos->Add(fHistTPCTOFMismatch[histoclass][iSpecies][iPtClass]);
1306 // -----------------------------------------------------------------------
1307 Bool_t AliAODTrackCutsDiHadronPID::InitializeGenMCHistos(Int_t histoclass) {
1309 cout<<"AliAODTrackCutsDiHadronPID - Creating Generated MC Histograms of Class "<<histoclass<<endl;
1310 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1312 // Primary Particles.
1313 fHistPrimGenMCPt[histoclass] = InitializePtSpectrum("fHistPrimGenMCPt",histoclass);
1314 fPrimGenMCTrackQAHistos->Add(fHistPrimGenMCPt[histoclass]);
1316 // Secondary Particles.
1317 fHistSecGenMCPt[histoclass] = InitializePtSpectrum("fHistSecGenMCPt",histoclass);
1318 fSecGenMCTrackQAHistos->Add(fHistSecGenMCPt[histoclass]);
1324 // -----------------------------------------------------------------------
1325 Bool_t AliAODTrackCutsDiHadronPID::InitializeRecMCHistos(Int_t histoclass) {
1327 cout<<"AliAODTrackCutsDiHadronPID - Creating Reconstructed MC Histograms of Class "<<histoclass<<endl;
1328 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1330 // Primary Particles.
1331 fHistPrimRecMCPt[histoclass] = InitializePtSpectrum("fHistPrimRecMCPt",histoclass);
1332 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCPt[histoclass]);
1334 fHistPrimRecNTracks[histoclass] = InitializeNTracksHisto("fHistPrimRecNTracks",histoclass);
1335 fPrimRecMCTrackQAHistos->Add(fHistPrimRecNTracks[histoclass]);
1337 fHistPrimRecMCDCA[histoclass] = InitializeDCASpectrum("fHistPrimRecDCA",histoclass);
1338 fPrimRecMCTrackQAHistos->Add(fHistPrimRecMCDCA[histoclass]);
1340 // Secondary Particles.
1341 fHistSecRecMCPt[histoclass] = InitializePtSpectrum("fHistSecRecMCPt",histoclass);
1342 fSecRecMCTrackQAHistos->Add(fHistSecRecMCPt[histoclass]);
1344 fHistSecRecMCDCAMat[histoclass] = InitializeDCASpectrum("fHistSecRecDCAMat",histoclass);
1345 fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAMat[histoclass]);
1347 fHistSecRecMCDCAWeak[histoclass] = InitializeDCASpectrum("fHistSecRecDCAWeak",histoclass);
1348 fSecRecMCTrackQAHistos->Add(fHistSecRecMCDCAWeak[histoclass]);
1354 // -----------------------------------------------------------------------
1355 TH1F* AliAODTrackCutsDiHadronPID::InitializePtSpectrum(const char* name, Int_t histoclass) {
1357 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1359 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1360 Form("p_{T} Spectrum (%s);p_{T} (GeV/c);Count",fHistoName[histoclass].Data()),fNPtBins,fPtAxis);
1362 hout->SetDirectory(0);
1368 // -----------------------------------------------------------------------
1369 TH3F* AliAODTrackCutsDiHadronPID::InitializePhiEtaPt(const char* name, Int_t histoclass) {
1371 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1373 TH3F* hout = AliHistToolsDiHadronPID::MakeHist3D(Form("%s%s",name,fHistoName[histoclass].Data()),
1374 Form("Spectrum (%s);#phi;#eta;p_{T} (GeV/c)",fHistoName[histoclass].Data()),
1375 fNPhiBins,0.,2.*TMath::Pi(),
1376 fNEtaBins,-fMaxEta,fMaxEta,
1379 hout->SetDirectory(0);
1385 // -----------------------------------------------------------------------
1386 TH2F* AliAODTrackCutsDiHadronPID::InitializeDCASpectrum(const char* name, Int_t histoclass) {
1388 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1390 TH2F* hout = new TH2F(Form("%s%s",name,fHistoName[histoclass].Data()),
1391 Form("DCA_{xy} (%s); p_{T} (GeV/c); DCA_{xy} (cm)",fHistoName[histoclass].Data()),
1395 hout->SetDirectory(0);
1401 // -----------------------------------------------------------------------
1402 TH1F* AliAODTrackCutsDiHadronPID::InitializeNTracksHisto(const char* name, Int_t histoclass) {
1404 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1406 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1407 Form("Number of Accepted Tracks (%s);N%s;Count",fHistoName[histoclass].Data(),fHistoLatex[histoclass].Data()),
1410 hout->SetDirectory(0);
1416 // -----------------------------------------------------------------------
1417 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAxyHisto(const char* name, Int_t histoclass) {
1419 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1421 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1422 Form("DCAxy (%s);DCAxy (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1424 hout->SetDirectory(0);
1430 // -----------------------------------------------------------------------
1431 TH1F* AliAODTrackCutsDiHadronPID::InitializeDCAzHisto(const char* name, Int_t histoclass) {
1433 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1435 TH1F* hout = new TH1F(Form("%s%s",name,fHistoName[histoclass].Data()),
1436 Form("DCAz (%s);DCAz (cm);Count",fHistoName[histoclass].Data()),300,-15.,15.);
1438 hout->SetDirectory(0);
1444 // -----------------------------------------------------------------------
1445 TH3F* AliAODTrackCutsDiHadronPID::InitializeAcceptanceHisto(const char* /*name*/, Int_t /*histoclass*/) {
1447 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1452 // -----------------------------------------------------------------------
1453 TH3F* AliAODTrackCutsDiHadronPID::InitializePIDHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1455 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1457 TH3F* hout = new TH3F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1458 Form("PID %s (Exp: %s);p_{T} (GeV/c);#Delta t (ps);dE/dx (a.u.)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
1459 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1460 fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies],
1461 fTPCbins[ptclass][expspecies],fTPCLowerBound[ptclass][expspecies],fTPCUpperBound[ptclass][expspecies]);
1463 hout->SetDirectory(0);
1469 // -----------------------------------------------------------------------
1470 TH2F* AliAODTrackCutsDiHadronPID::InitializeTOFMismatchHisto(const char* name, Int_t histoclass, Int_t expspecies, Int_t ptclass) {
1472 if (fDebug > 1) {cout << Form("File: %s, Line: %i, Function: %s",__FILE__,__LINE__,__func__) << endl;}
1474 TH2F* hout = new TH2F(Form("%s%s%s%s",name,fHistoName[histoclass].Data(),fParticleName[expspecies].Data(),fPtClassName[ptclass].Data()),
1475 Form("TOF Mismatch %s (Exp: %s);p_{T} (GeV/c);#Delta t (ps)",fHistoName[histoclass].Data(),fParticleName[expspecies].Data()),
1476 fNPtBinsPID[ptclass],fPtBoundaryPID[ptclass],fPtBoundaryPID[ptclass+1],
1477 fTOFbins[ptclass][expspecies],fTOFLowerBound[ptclass][expspecies],fTOFUpperBound[ptclass][expspecies]);
1479 hout->SetDirectory(0);