1 /**************************************************************************
2 * Copyright(c) 1998-1999, 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 /* $Id: AliESDtrackCuts.cxx 24534 2008-03-16 22:22:11Z fca $ */
18 #include "AliESDtrackCuts.h"
20 #include <AliESDtrack.h>
21 #include <AliESDVertex.h>
22 #include <AliESDEvent.h>
27 #include <TDirectory.h>
31 //____________________________________________________________________
32 ClassImp(AliESDtrackCuts)
35 const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
40 "#Chi^{2}/cluster TPC",
41 "#Chi^{2}/cluster ITS",
57 "trk-to-vtx max dca 2D absolute",
58 "trk-to-vtx max dca xy absolute",
59 "trk-to-vtx max dca z absolute",
60 "trk-to-vtx min dca 2D absolute",
61 "trk-to-vtx min dca xy absolute",
62 "trk-to-vtx min dca z absolute",
63 "SPD cluster requirement",
64 "SDD cluster requirement",
65 "SSD cluster requirement",
66 "require ITS stand-alone"
69 //____________________________________________________________________
70 AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliAnalysisCuts(name,title),
71 fCutMinNClusterTPC(0),
72 fCutMinNClusterITS(0),
73 fCutMaxChi2PerClusterTPC(0),
74 fCutMaxChi2PerClusterITS(0),
80 fCutAcceptKinkDaughters(0),
81 fCutRequireTPCRefit(0),
82 fCutRequireITSRefit(0),
83 fCutRequireITSStandAlone(0),
84 fCutNsigmaToVertex(0),
85 fCutSigmaToVertexRequired(0),
86 fCutMaxDCAToVertexXY(0),
87 fCutMaxDCAToVertexZ(0),
88 fCutMinDCAToVertexXY(0),
89 fCutMinDCAToVertexZ(0),
116 //##############################################################################
117 // setting default cuts
118 SetMinNClustersTPC();
119 SetMinNClustersITS();
120 SetMaxChi2PerClusterTPC();
121 SetMaxChi2PerClusterITS();
122 SetMaxCovDiagonalElements();
123 SetRequireTPCRefit();
124 SetRequireITSRefit();
125 SetRequireITSStandAlone(kFALSE);
126 SetAcceptKinkDaughters();
127 SetMaxNsigmaToVertex();
128 SetMaxDCAToVertexXY();
129 SetMaxDCAToVertexZ();
131 SetMinDCAToVertexXY();
132 SetMinDCAToVertexZ();
140 SetClusterRequirementITS(kSPD);
141 SetClusterRequirementITS(kSDD);
142 SetClusterRequirementITS(kSSD);
147 //_____________________________________________________________________________
148 AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
149 fCutMinNClusterTPC(0),
150 fCutMinNClusterITS(0),
151 fCutMaxChi2PerClusterTPC(0),
152 fCutMaxChi2PerClusterITS(0),
158 fCutAcceptKinkDaughters(0),
159 fCutRequireTPCRefit(0),
160 fCutRequireITSRefit(0),
161 fCutRequireITSStandAlone(0),
162 fCutNsigmaToVertex(0),
163 fCutSigmaToVertexRequired(0),
164 fCutMaxDCAToVertexXY(0),
165 fCutMaxDCAToVertexZ(0),
166 fCutMinDCAToVertexXY(0),
167 fCutMinDCAToVertexZ(0),
168 fCutDCAToVertex2D(0),
192 ((AliESDtrackCuts &) c).Copy(*this);
195 AliESDtrackCuts::~AliESDtrackCuts()
201 for (Int_t i=0; i<2; i++) {
203 if (fhNClustersITS[i])
204 delete fhNClustersITS[i];
205 if (fhNClustersTPC[i])
206 delete fhNClustersTPC[i];
207 if (fhChi2PerClusterITS[i])
208 delete fhChi2PerClusterITS[i];
209 if (fhChi2PerClusterTPC[i])
210 delete fhChi2PerClusterTPC[i];
231 if (fhDXYNormalized[i])
232 delete fhDXYNormalized[i];
233 if (fhDZNormalized[i])
234 delete fhDZNormalized[i];
235 if (fhDXYvsDZNormalized[i])
236 delete fhDXYvsDZNormalized[i];
237 if (fhNSigmaToVertex[i])
238 delete fhNSigmaToVertex[i];
246 delete ffDTheoretical;
249 delete fhCutStatistics;
250 if (fhCutCorrelation)
251 delete fhCutCorrelation;
254 void AliESDtrackCuts::Init()
257 // sets everything to zero
260 fCutMinNClusterTPC = 0;
261 fCutMinNClusterITS = 0;
263 fCutMaxChi2PerClusterTPC = 0;
264 fCutMaxChi2PerClusterITS = 0;
266 for (Int_t i = 0; i < 3; i++)
267 fCutClusterRequirementITS[i] = kOff;
275 fCutAcceptKinkDaughters = 0;
276 fCutRequireTPCRefit = 0;
277 fCutRequireITSRefit = 0;
278 fCutRequireITSStandAlone = 0;
280 fCutNsigmaToVertex = 0;
281 fCutSigmaToVertexRequired = 0;
282 fCutMaxDCAToVertexXY = 0;
283 fCutMaxDCAToVertexZ = 0;
284 fCutDCAToVertex2D = 0;
285 fCutMinDCAToVertexXY = 0;
286 fCutMinDCAToVertexZ = 0;
304 fHistogramsOn = kFALSE;
306 for (Int_t i=0; i<2; ++i)
308 fhNClustersITS[i] = 0;
309 fhNClustersTPC[i] = 0;
311 fhChi2PerClusterITS[i] = 0;
312 fhChi2PerClusterTPC[i] = 0;
325 fhDXYNormalized[i] = 0;
326 fhDZNormalized[i] = 0;
327 fhDXYvsDZNormalized[i] = 0;
328 fhNSigmaToVertex[i] = 0;
336 fhCutCorrelation = 0;
339 //_____________________________________________________________________________
340 AliESDtrackCuts &AliESDtrackCuts::operator=(const AliESDtrackCuts &c)
343 // Assignment operator
346 if (this != &c) ((AliESDtrackCuts &) c).Copy(*this);
350 //_____________________________________________________________________________
351 void AliESDtrackCuts::Copy(TObject &c) const
357 AliESDtrackCuts& target = (AliESDtrackCuts &) c;
361 target.fCutMinNClusterTPC = fCutMinNClusterTPC;
362 target.fCutMinNClusterITS = fCutMinNClusterITS;
364 target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
365 target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
367 for (Int_t i = 0; i < 3; i++)
368 target.fCutClusterRequirementITS[i] = fCutClusterRequirementITS[i];
370 target.fCutMaxC11 = fCutMaxC11;
371 target.fCutMaxC22 = fCutMaxC22;
372 target.fCutMaxC33 = fCutMaxC33;
373 target.fCutMaxC44 = fCutMaxC44;
374 target.fCutMaxC55 = fCutMaxC55;
376 target.fCutAcceptKinkDaughters = fCutAcceptKinkDaughters;
377 target.fCutRequireTPCRefit = fCutRequireTPCRefit;
378 target.fCutRequireITSRefit = fCutRequireITSRefit;
379 target.fCutRequireITSStandAlone = fCutRequireITSStandAlone;
381 target.fCutNsigmaToVertex = fCutNsigmaToVertex;
382 target.fCutSigmaToVertexRequired = fCutSigmaToVertexRequired;
383 target.fCutMaxDCAToVertexXY = fCutMaxDCAToVertexXY;
384 target.fCutMaxDCAToVertexZ = fCutMaxDCAToVertexZ;
385 target.fCutDCAToVertex2D = fCutDCAToVertex2D;
386 target.fCutMinDCAToVertexXY = fCutMinDCAToVertexXY;
387 target.fCutMinDCAToVertexZ = fCutMinDCAToVertexZ;
389 target.fPMin = fPMin;
390 target.fPMax = fPMax;
391 target.fPtMin = fPtMin;
392 target.fPtMax = fPtMax;
393 target.fPxMin = fPxMin;
394 target.fPxMax = fPxMax;
395 target.fPyMin = fPyMin;
396 target.fPyMax = fPyMax;
397 target.fPzMin = fPzMin;
398 target.fPzMax = fPzMax;
399 target.fEtaMin = fEtaMin;
400 target.fEtaMax = fEtaMax;
401 target.fRapMin = fRapMin;
402 target.fRapMax = fRapMax;
404 target.fHistogramsOn = fHistogramsOn;
406 for (Int_t i=0; i<2; ++i)
408 if (fhNClustersITS[i]) target.fhNClustersITS[i] = (TH1F*) fhNClustersITS[i]->Clone();
409 if (fhNClustersTPC[i]) target.fhNClustersTPC[i] = (TH1F*) fhNClustersTPC[i]->Clone();
411 if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
412 if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
414 if (fhC11[i]) target.fhC11[i] = (TH1F*) fhC11[i]->Clone();
415 if (fhC22[i]) target.fhC22[i] = (TH1F*) fhC22[i]->Clone();
416 if (fhC33[i]) target.fhC33[i] = (TH1F*) fhC33[i]->Clone();
417 if (fhC44[i]) target.fhC44[i] = (TH1F*) fhC44[i]->Clone();
418 if (fhC55[i]) target.fhC55[i] = (TH1F*) fhC55[i]->Clone();
420 if (fhDXY[i]) target.fhDXY[i] = (TH1F*) fhDXY[i]->Clone();
421 if (fhDZ[i]) target.fhDZ[i] = (TH1F*) fhDZ[i]->Clone();
422 if (fhDXYDZ[i]) target.fhDXYDZ[i] = (TH1F*) fhDXYDZ[i]->Clone();
423 if (fhDXYvsDZ[i]) target.fhDXYvsDZ[i] = (TH2F*) fhDXYvsDZ[i]->Clone();
425 if (fhDXYNormalized[i]) target.fhDXYNormalized[i] = (TH1F*) fhDXYNormalized[i]->Clone();
426 if (fhDZNormalized[i]) target.fhDZNormalized[i] = (TH1F*) fhDZNormalized[i]->Clone();
427 if (fhDXYvsDZNormalized[i]) target.fhDXYvsDZNormalized[i] = (TH2F*) fhDXYvsDZNormalized[i]->Clone();
428 if (fhNSigmaToVertex[i]) target.fhNSigmaToVertex[i] = (TH1F*) fhNSigmaToVertex[i]->Clone();
430 if (fhPt[i]) target.fhPt[i] = (TH1F*) fhPt[i]->Clone();
431 if (fhEta[i]) target.fhEta[i] = (TH1F*) fhEta[i]->Clone();
433 if (ffDTheoretical) target.ffDTheoretical = (TF1*) ffDTheoretical->Clone();
435 if (fhCutStatistics) target.fhCutStatistics = (TH1F*) fhCutStatistics->Clone();
436 if (fhCutCorrelation) target.fhCutCorrelation = (TH2F*) fhCutCorrelation->Clone();
441 //_____________________________________________________________________________
442 Long64_t AliESDtrackCuts::Merge(TCollection* list) {
443 // Merge a list of AliESDtrackCuts objects with this (needed for PROOF)
444 // Returns the number of merged objects (including this)
451 TIterator* iter = list->MakeIterator();
454 // collection of measured and generated histograms
456 while ((obj = iter->Next())) {
458 AliESDtrackCuts* entry = dynamic_cast<AliESDtrackCuts*>(obj);
462 if (!entry->fHistogramsOn)
465 for (Int_t i=0; i<2; i++) {
467 fhNClustersITS[i] ->Add(entry->fhNClustersITS[i] );
468 fhNClustersTPC[i] ->Add(entry->fhNClustersTPC[i] );
470 fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]);
471 fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]);
473 fhC11[i] ->Add(entry->fhC11[i] );
474 fhC22[i] ->Add(entry->fhC22[i] );
475 fhC33[i] ->Add(entry->fhC33[i] );
476 fhC44[i] ->Add(entry->fhC44[i] );
477 fhC55[i] ->Add(entry->fhC55[i] );
479 fhDXY[i] ->Add(entry->fhDXY[i] );
480 fhDZ[i] ->Add(entry->fhDZ[i] );
481 fhDXYDZ[i] ->Add(entry->fhDXYDZ[i] );
482 fhDXYvsDZ[i] ->Add(entry->fhDXYvsDZ[i] );
484 fhDXYNormalized[i] ->Add(entry->fhDXYNormalized[i] );
485 fhDZNormalized[i] ->Add(entry->fhDZNormalized[i] );
486 fhDXYvsDZNormalized[i] ->Add(entry->fhDXYvsDZNormalized[i]);
487 fhNSigmaToVertex[i] ->Add(entry->fhNSigmaToVertex[i]);
489 fhPt[i] ->Add(entry->fhPt[i]);
490 fhEta[i] ->Add(entry->fhEta[i]);
493 fhCutStatistics ->Add(entry->fhCutStatistics);
494 fhCutCorrelation ->Add(entry->fhCutCorrelation);
502 //____________________________________________________________________
503 Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* esdTrack)
505 // Calculates the number of sigma to the vertex.
510 esdTrack->GetImpactParameters(b,bCov);
512 if (bCov[0]<=0 || bCov[2]<=0) {
513 AliDebugClass(1, "Estimated b resolution lower or equal zero!");
514 bCov[0]=0; bCov[2]=0;
516 bRes[0] = TMath::Sqrt(bCov[0]);
517 bRes[1] = TMath::Sqrt(bCov[2]);
519 // -----------------------------------
520 // How to get to a n-sigma cut?
522 // The accumulated statistics from 0 to d is
524 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
525 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
527 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-d**2)/2)
528 // Can this be expressed in a different way?
530 if (bRes[0] == 0 || bRes[1] ==0)
533 Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
535 // work around precision problem
536 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
537 // 1e-15 corresponds to nsigma ~ 7.7
538 if (TMath::Exp(-d * d / 2) < 1e-15)
541 Float_t nSigma = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
545 void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
547 // enables the branches needed by AcceptTrack, for a list see comment of AcceptTrack
549 tree->SetBranchStatus("fTracks.fFlags", 1);
550 tree->SetBranchStatus("fTracks.fITSncls", 1);
551 tree->SetBranchStatus("fTracks.fTPCncls", 1);
552 tree->SetBranchStatus("fTracks.fITSchi2", 1);
553 tree->SetBranchStatus("fTracks.fTPCchi2", 1);
554 tree->SetBranchStatus("fTracks.fC*", 1);
555 tree->SetBranchStatus("fTracks.fD", 1);
556 tree->SetBranchStatus("fTracks.fZ", 1);
557 tree->SetBranchStatus("fTracks.fCdd", 1);
558 tree->SetBranchStatus("fTracks.fCdz", 1);
559 tree->SetBranchStatus("fTracks.fCzz", 1);
560 tree->SetBranchStatus("fTracks.fP*", 1);
561 tree->SetBranchStatus("fTracks.fR*", 1);
562 tree->SetBranchStatus("fTracks.fKinkIndexes*", 1);
565 //____________________________________________________________________
566 Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
569 // figure out if the tracks survives all the track cuts defined
571 // the different quality parameter and kinematic values are first
572 // retrieved from the track. then it is found out what cuts the
573 // track did not survive and finally the cuts are imposed.
575 // this function needs the following branches:
581 // fTracks.fC //GetExternalCovariance
582 // fTracks.fD //GetImpactParameters
583 // fTracks.fZ //GetImpactParameters
584 // fTracks.fCdd //GetImpactParameters
585 // fTracks.fCdz //GetImpactParameters
586 // fTracks.fCzz //GetImpactParameters
587 // fTracks.fP //GetPxPyPz
588 // fTracks.fR //GetMass
589 // fTracks.fP //GetMass
590 // fTracks.fKinkIndexes
592 UInt_t status = esdTrack->GetStatus();
594 // getting quality parameters from the ESD track
595 Int_t nClustersITS = esdTrack->GetITSclusters(0);
596 Int_t nClustersTPC = esdTrack->GetTPCclusters(0);
598 Float_t chi2PerClusterITS = -1;
599 Float_t chi2PerClusterTPC = -1;
601 chi2PerClusterITS = esdTrack->GetITSchi2()/Float_t(nClustersITS);
603 chi2PerClusterTPC = esdTrack->GetTPCchi2()/Float_t(nClustersTPC);
605 esdTrack->GetExternalCovariance(extCov);
607 // getting the track to vertex parameters
608 Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
612 esdTrack->GetImpactParameters(b,bCov);
613 if (bCov[0]<=0 || bCov[2]<=0) {
614 AliDebug(1, "Estimated b resolution lower or equal zero!");
615 bCov[0]=0; bCov[2]=0;
618 Float_t dcaToVertexXY = b[0];
619 Float_t dcaToVertexZ = b[1];
621 Float_t dcaToVertex = -1;
623 if (fCutDCAToVertex2D)
625 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY/fCutMaxDCAToVertexXY/fCutMaxDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMaxDCAToVertexZ/fCutMaxDCAToVertexZ);
628 dcaToVertex = TMath::Sqrt(dcaToVertexXY*dcaToVertexXY + dcaToVertexZ*dcaToVertexZ);
630 // getting the kinematic variables of the track
631 // (assuming the mass is known)
633 esdTrack->GetPxPyPz(p);
635 Float_t momentum = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2) + TMath::Power(p[2],2));
636 Float_t pt = TMath::Sqrt(TMath::Power(p[0],2) + TMath::Power(p[1],2));
637 Float_t energy = TMath::Sqrt(TMath::Power(esdTrack->GetMass(),2) + TMath::Power(momentum,2));
640 //y-eta related calculations
643 if((momentum != TMath::Abs(p[2]))&&(momentum != 0))
644 eta = 0.5*TMath::Log((momentum + p[2])/(momentum - p[2]));
645 if((energy != TMath::Abs(p[2]))&&(momentum != 0))
646 y = 0.5*TMath::Log((energy + p[2])/(energy - p[2]));
649 //########################################################################
653 for (Int_t i=0; i<kNCuts; i++) cuts[i]=kFALSE;
655 // track quality cuts
656 if (fCutRequireTPCRefit && (status&AliESDtrack::kTPCrefit)==0)
658 if (fCutRequireITSRefit && (status&AliESDtrack::kITSrefit)==0)
660 if (nClustersTPC<fCutMinNClusterTPC)
662 if (nClustersITS<fCutMinNClusterITS)
664 if (chi2PerClusterTPC>fCutMaxChi2PerClusterTPC)
666 if (chi2PerClusterITS>fCutMaxChi2PerClusterITS)
668 if (extCov[0] > fCutMaxC11)
670 if (extCov[2] > fCutMaxC22)
672 if (extCov[5] > fCutMaxC33)
674 if (extCov[9] > fCutMaxC44)
676 if (extCov[14] > fCutMaxC55)
678 if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
680 // if n sigma could not be calculated
681 if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
683 if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
685 // track kinematics cut
686 if((momentum < fPMin) || (momentum > fPMax))
688 if((pt < fPtMin) || (pt > fPtMax))
690 if((p[0] < fPxMin) || (p[0] > fPxMax))
692 if((p[1] < fPyMin) || (p[1] > fPyMax))
694 if((p[2] < fPzMin) || (p[2] > fPzMax))
696 if((eta < fEtaMin) || (eta > fEtaMax))
698 if((y < fRapMin) || (y > fRapMax))
700 if (fCutDCAToVertex2D && dcaToVertex > 1)
702 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) > fCutMaxDCAToVertexXY)
704 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) > fCutMaxDCAToVertexZ)
706 if (fCutDCAToVertex2D && fCutMinDCAToVertexXY > 0 && fCutMinDCAToVertexZ > 0 && dcaToVertexXY*dcaToVertexXY/fCutMinDCAToVertexXY/fCutMinDCAToVertexXY + dcaToVertexZ*dcaToVertexZ/fCutMinDCAToVertexZ/fCutMinDCAToVertexZ < 1)
708 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexXY) < fCutMinDCAToVertexXY)
710 if (!fCutDCAToVertex2D && TMath::Abs(dcaToVertexZ) < fCutMinDCAToVertexZ)
713 for (Int_t i = 0; i < 3; i++)
714 cuts[27+i] = !CheckITSClusterRequirement(fCutClusterRequirementITS[i], esdTrack->HasPointOnITSLayer(i*2), esdTrack->HasPointOnITSLayer(i*2+1));
716 if (fCutRequireITSStandAlone && ((status & AliESDtrack::kITSin) == 0 || (status & AliESDtrack::kTPCin)))
720 for (Int_t i=0; i<kNCuts; i++)
721 if (cuts[i]) {cut = kTRUE;}
725 //########################################################################
726 // filling histograms
728 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n tracks")));
730 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin("n cut tracks")));
732 for (Int_t i=0; i<kNCuts; i++) {
734 fhCutStatistics->Fill(fhCutStatistics->GetBinCenter(fhCutStatistics->GetXaxis()->FindBin(fgkCutNames[i])));
736 for (Int_t j=i; j<kNCuts; j++) {
737 if (cuts[i] && cuts[j]) {
738 Float_t xC = fhCutCorrelation->GetXaxis()->GetBinCenter(fhCutCorrelation->GetXaxis()->FindBin(fgkCutNames[i]));
739 Float_t yC = fhCutCorrelation->GetYaxis()->GetBinCenter(fhCutCorrelation->GetYaxis()->FindBin(fgkCutNames[j]));
740 fhCutCorrelation->Fill(xC, yC);
746 // now we loop over the filling of the histograms twice: once "before" the cut, once "after"
747 // the code is not in a function due to too many local variables that would need to be passed
749 for (Int_t id = 0; id < 2; id++)
751 // id = 0 --> before cut
752 // id = 1 --> after cut
756 fhNClustersITS[id]->Fill(nClustersITS);
757 fhNClustersTPC[id]->Fill(nClustersTPC);
758 fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
759 fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
761 fhC11[id]->Fill(extCov[0]);
762 fhC22[id]->Fill(extCov[2]);
763 fhC33[id]->Fill(extCov[5]);
764 fhC44[id]->Fill(extCov[9]);
765 fhC55[id]->Fill(extCov[14]);
768 fhEta[id]->Fill(eta);
771 bRes[0] = TMath::Sqrt(bCov[0]);
772 bRes[1] = TMath::Sqrt(bCov[2]);
774 fhDZ[id]->Fill(b[1]);
775 fhDXY[id]->Fill(b[0]);
776 fhDXYDZ[id]->Fill(dcaToVertex);
777 fhDXYvsDZ[id]->Fill(b[1],b[0]);
779 if (bRes[0]!=0 && bRes[1]!=0) {
780 fhDZNormalized[id]->Fill(b[1]/bRes[1]);
781 fhDXYNormalized[id]->Fill(b[0]/bRes[0]);
782 fhDXYvsDZNormalized[id]->Fill(b[1]/bRes[1], b[0]/bRes[0]);
783 fhNSigmaToVertex[id]->Fill(nSigmaToVertex);
795 //____________________________________________________________________
796 Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bool_t clusterL1, Bool_t clusterL2)
798 // checks if the cluster requirement is fullfilled (in this case: return kTRUE)
802 case kOff: return kTRUE;
803 case kNone: return !clusterL1 && !clusterL2;
804 case kAny: return clusterL1 || clusterL2;
805 case kFirst: return clusterL1;
806 case kOnlyFirst: return clusterL1 && !clusterL2;
807 case kSecond: return clusterL2;
808 case kOnlySecond: return clusterL2 && !clusterL1;
809 case kBoth: return clusterL1 && clusterL2;
815 //____________________________________________________________________
816 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(AliESDEvent* esd, Int_t iTrack)
818 // creates a TPC only track from the given esd track
819 // the track has to be deleted by the user
821 // NB. most of the functionality to get a TPC only track from an ESD track is in AliESDtrack, where it should be
822 // there are only missing propagations here that are needed for old data
823 // this function will therefore become obsolete
825 // adapted from code provided by CKB
827 if (!esd->GetPrimaryVertexTPC())
828 return 0; // No TPC vertex no TPC tracks
830 if(!esd->GetPrimaryVertexTPC()->GetStatus())
831 return 0; // TPC Vertex is created by default in AliESDEvent, do not use in this case
833 AliESDtrack* track = esd->GetTrack(iTrack);
837 AliESDtrack *tpcTrack = new AliESDtrack();
839 // This should have been done during the reconstruction
840 // fixed by Juri in r26675
841 // but recalculate for older data CKB
843 track->GetImpactParametersTPC(p,cov);
845 track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
848 // only true if we have a tpc track
849 if (!track->FillTPCOnlyTrack(*tpcTrack))
855 // propagate to Vertex
856 // not needed for normal reconstructed ESDs...
857 // Double_t pTPC[2],covTPC[3];
858 // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000, pTPC, covTPC);
863 //____________________________________________________________________
864 TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
867 // returns an array of all tracks that pass the cuts
868 // or an array of TPC only tracks (propagated to the TPC vertex during reco)
869 // tracks that pass the cut
871 TObjArray* acceptedTracks = new TObjArray();
873 // loop over esd tracks
874 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
876 if(!esd->GetPrimaryVertexTPC())return acceptedTracks; // No TPC vertex no TPC tracks
877 if(!esd->GetPrimaryVertexTPC()->GetStatus())return acceptedTracks; // No proper TPC vertex, only the default
879 AliESDtrack *tpcTrack = GetTPCOnlyTrack(esd, iTrack);
883 if (AcceptTrack(tpcTrack)) {
884 acceptedTracks->Add(tpcTrack);
891 AliESDtrack* track = esd->GetTrack(iTrack);
892 if(AcceptTrack(track))
893 acceptedTracks->Add(track);
896 if(bTPC)acceptedTracks->SetOwner(kTRUE);
897 return acceptedTracks;
900 //____________________________________________________________________
901 Int_t AliESDtrackCuts::CountAcceptedTracks(AliESDEvent* esd)
904 // returns an the number of tracks that pass the cuts
909 // loop over esd tracks
910 for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
911 AliESDtrack* track = esd->GetTrack(iTrack);
912 if (AcceptTrack(track))
919 //____________________________________________________________________
920 void AliESDtrackCuts::DefineHistograms(Int_t color) {
922 // diagnostics histograms are defined
927 Bool_t oldStatus = TH1::AddDirectoryStatus();
928 TH1::AddDirectory(kFALSE);
930 //###################################################################################
931 // defining histograms
933 fhCutStatistics = new TH1F("cut_statistics","cut statistics",kNCuts+4,-0.5,kNCuts+3.5);
935 fhCutStatistics->GetXaxis()->SetBinLabel(1,"n tracks");
936 fhCutStatistics->GetXaxis()->SetBinLabel(2,"n cut tracks");
938 fhCutCorrelation = new TH2F("cut_correlation","cut correlation",kNCuts,-0.5,kNCuts-0.5,kNCuts,-0.5,kNCuts-0.5);;
940 for (Int_t i=0; i<kNCuts; i++) {
941 fhCutStatistics->GetXaxis()->SetBinLabel(i+4,fgkCutNames[i]);
942 fhCutCorrelation->GetXaxis()->SetBinLabel(i+1,fgkCutNames[i]);
943 fhCutCorrelation->GetYaxis()->SetBinLabel(i+1,fgkCutNames[i]);
946 fhCutStatistics ->SetLineColor(color);
947 fhCutCorrelation ->SetLineColor(color);
948 fhCutStatistics ->SetLineWidth(2);
949 fhCutCorrelation ->SetLineWidth(2);
951 for (Int_t i=0; i<2; i++) {
952 fhNClustersITS[i] = new TH1F("nClustersITS" ,"",8,-0.5,7.5);
953 fhNClustersTPC[i] = new TH1F("nClustersTPC" ,"",165,-0.5,164.5);
954 fhChi2PerClusterITS[i] = new TH1F("chi2PerClusterITS","",500,0,10);
955 fhChi2PerClusterTPC[i] = new TH1F("chi2PerClusterTPC","",500,0,10);
957 fhC11[i] = new TH1F("covMatrixDiagonal11","",2000,0,20);
958 fhC22[i] = new TH1F("covMatrixDiagonal22","",2000,0,20);
959 fhC33[i] = new TH1F("covMatrixDiagonal33","",1000,0,0.1);
960 fhC44[i] = new TH1F("covMatrixDiagonal44","",1000,0,0.1);
961 fhC55[i] = new TH1F("covMatrixDiagonal55","",1000,0,5);
963 fhDXY[i] = new TH1F("dXY" ,"",500,-10,10);
964 fhDZ[i] = new TH1F("dZ" ,"",500,-10,10);
965 fhDXYDZ[i] = new TH1F("dXYDZ" ,"",500,0,10);
966 fhDXYvsDZ[i] = new TH2F("dXYvsDZ","",200,-10,10,200,-10,10);
968 fhDXYNormalized[i] = new TH1F("dXYNormalized" ,"",500,-10,10);
969 fhDZNormalized[i] = new TH1F("dZNormalized" ,"",500,-10,10);
970 fhDXYvsDZNormalized[i] = new TH2F("dXYvsDZNormalized","",200,-10,10,200,-10,10);
972 fhNSigmaToVertex[i] = new TH1F("nSigmaToVertex","",500,0,10);
974 fhPt[i] = new TH1F("pt" ,"p_{T} distribution;p_{T} (GeV/c)", 800, 0.0, 10.0);
975 fhEta[i] = new TH1F("eta" ,"#eta distribution;#eta",40,-2.0,2.0);
977 fhNClustersITS[i]->SetTitle("n ITS clusters");
978 fhNClustersTPC[i]->SetTitle("n TPC clusters");
979 fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
980 fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
982 fhC11[i]->SetTitle("cov 11 : #sigma_{y}^{2} [cm^{2}]");
983 fhC22[i]->SetTitle("cov 22 : #sigma_{z}^{2} [cm^{2}]");
984 fhC33[i]->SetTitle("cov 33 : #sigma_{sin(#phi)}^{2}");
985 fhC44[i]->SetTitle("cov 44 : #sigma_{tan(#theta_{dip})}^{2}");
986 fhC55[i]->SetTitle("cov 55 : #sigma_{1/p_{T}}^{2} [(c/GeV)^2]");
988 fhDXY[i]->SetXTitle("transverse impact parameter (cm)");
989 fhDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
990 fhDXYDZ[i]->SetTitle("absolute impact parameter;sqrt(dXY**2 + dZ**2) (cm)");
991 fhDXYvsDZ[i]->SetXTitle("longitudinal impact parameter (cm)");
992 fhDXYvsDZ[i]->SetYTitle("transverse impact parameter (cm)");
994 fhDXYNormalized[i]->SetTitle("normalized trans impact par (n#sigma)");
995 fhDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
996 fhDXYvsDZNormalized[i]->SetTitle("normalized long impact par (n#sigma)");
997 fhDXYvsDZNormalized[i]->SetYTitle("normalized trans impact par (n#sigma)");
998 fhNSigmaToVertex[i]->SetTitle("n #sigma to vertex");
1000 fhNClustersITS[i]->SetLineColor(color); fhNClustersITS[i]->SetLineWidth(2);
1001 fhNClustersTPC[i]->SetLineColor(color); fhNClustersTPC[i]->SetLineWidth(2);
1002 fhChi2PerClusterITS[i]->SetLineColor(color); fhChi2PerClusterITS[i]->SetLineWidth(2);
1003 fhChi2PerClusterTPC[i]->SetLineColor(color); fhChi2PerClusterTPC[i]->SetLineWidth(2);
1005 fhC11[i]->SetLineColor(color); fhC11[i]->SetLineWidth(2);
1006 fhC22[i]->SetLineColor(color); fhC22[i]->SetLineWidth(2);
1007 fhC33[i]->SetLineColor(color); fhC33[i]->SetLineWidth(2);
1008 fhC44[i]->SetLineColor(color); fhC44[i]->SetLineWidth(2);
1009 fhC55[i]->SetLineColor(color); fhC55[i]->SetLineWidth(2);
1011 fhDXY[i]->SetLineColor(color); fhDXY[i]->SetLineWidth(2);
1012 fhDZ[i]->SetLineColor(color); fhDZ[i]->SetLineWidth(2);
1013 fhDXYDZ[i]->SetLineColor(color); fhDXYDZ[i]->SetLineWidth(2);
1015 fhDXYNormalized[i]->SetLineColor(color); fhDXYNormalized[i]->SetLineWidth(2);
1016 fhDZNormalized[i]->SetLineColor(color); fhDZNormalized[i]->SetLineWidth(2);
1017 fhNSigmaToVertex[i]->SetLineColor(color); fhNSigmaToVertex[i]->SetLineWidth(2);
1020 // The number of sigmas to the vertex is per definition gaussian
1021 ffDTheoretical = new TF1("nSigmaToVertexTheoretical","([0]/2.506628274)*exp(-(x**2)/2)",0,50);
1022 ffDTheoretical->SetParameter(0,1);
1024 TH1::AddDirectory(oldStatus);
1027 //____________________________________________________________________
1028 Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
1031 // loads the histograms from a file
1032 // if dir is empty a directory with the name of this object is taken (like in SaveHistogram)
1038 if (!gDirectory->cd(dir))
1041 ffDTheoretical = dynamic_cast<TF1*> (gDirectory->Get("nSigmaToVertexTheory"));
1043 fhCutStatistics = dynamic_cast<TH1F*> (gDirectory->Get("cut_statistics"));
1044 fhCutCorrelation = dynamic_cast<TH2F*> (gDirectory->Get("cut_correlation"));
1046 for (Int_t i=0; i<2; i++) {
1049 gDirectory->cd("before_cuts");
1052 gDirectory->cd("after_cuts");
1054 fhNClustersITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("nClustersITS" ));
1055 fhNClustersTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("nClustersTPC" ));
1056 fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
1057 fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
1059 fhC11[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal11"));
1060 fhC22[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal22"));
1061 fhC33[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal33"));
1062 fhC44[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal44"));
1063 fhC55[i] = dynamic_cast<TH1F*> (gDirectory->Get("covMatrixDiagonal55"));
1065 fhDXY[i] = dynamic_cast<TH1F*> (gDirectory->Get("dXY" ));
1066 fhDZ[i] = dynamic_cast<TH1F*> (gDirectory->Get("dZ" ));
1067 fhDXYDZ[i] = dynamic_cast<TH1F*> (gDirectory->Get("dXYDZ"));
1068 fhDXYvsDZ[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZ"));
1070 fhDXYNormalized[i] = dynamic_cast<TH1F*> (gDirectory->Get("dXYNormalized" ));
1071 fhDZNormalized[i] = dynamic_cast<TH1F*> (gDirectory->Get("dZNormalized" ));
1072 fhDXYvsDZNormalized[i] = dynamic_cast<TH2F*> (gDirectory->Get("dXYvsDZNormalized"));
1073 fhNSigmaToVertex[i] = dynamic_cast<TH1F*> (gDirectory->Get("nSigmaToVertex"));
1075 fhPt[i] = dynamic_cast<TH1F*> (gDirectory->Get("pt"));
1076 fhEta[i] = dynamic_cast<TH1F*> (gDirectory->Get("eta"));
1078 gDirectory->cd("../");
1081 gDirectory->cd("..");
1086 //____________________________________________________________________
1087 void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
1089 // saves the histograms in a directory (dir)
1092 if (!fHistogramsOn) {
1093 AliDebug(0, "Histograms not on - cannot save histograms!!!");
1100 gDirectory->mkdir(dir);
1101 gDirectory->cd(dir);
1103 gDirectory->mkdir("before_cuts");
1104 gDirectory->mkdir("after_cuts");
1106 // a factor of 2 is needed since n sigma is positive
1107 ffDTheoretical->SetParameter(0,2*fhNSigmaToVertex[0]->Integral("width"));
1108 ffDTheoretical->Write("nSigmaToVertexTheory");
1110 fhCutStatistics->Write();
1111 fhCutCorrelation->Write();
1113 for (Int_t i=0; i<2; i++) {
1115 gDirectory->cd("before_cuts");
1117 gDirectory->cd("after_cuts");
1119 fhNClustersITS[i] ->Write();
1120 fhNClustersTPC[i] ->Write();
1121 fhChi2PerClusterITS[i] ->Write();
1122 fhChi2PerClusterTPC[i] ->Write();
1132 fhDXYDZ[i] ->Write();
1133 fhDXYvsDZ[i] ->Write();
1135 fhDXYNormalized[i] ->Write();
1136 fhDZNormalized[i] ->Write();
1137 fhDXYvsDZNormalized[i] ->Write();
1138 fhNSigmaToVertex[i] ->Write();
1143 gDirectory->cd("../");
1146 gDirectory->cd("../");
1149 //____________________________________________________________________
1150 void AliESDtrackCuts::DrawHistograms()
1152 // draws some histograms
1154 TCanvas* canvas1 = new TCanvas(Form("%s_1", GetName()), "Track Quality Results1", 800, 800);
1155 canvas1->Divide(2, 2);
1158 fhNClustersTPC[0]->SetStats(kFALSE);
1159 fhNClustersTPC[0]->Draw();
1162 fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1163 fhChi2PerClusterTPC[0]->Draw();
1166 fhNSigmaToVertex[0]->SetStats(kFALSE);
1167 fhNSigmaToVertex[0]->GetXaxis()->SetRangeUser(0, 10);
1168 fhNSigmaToVertex[0]->Draw();
1170 canvas1->SaveAs(Form("%s_%s.gif", GetName(), canvas1->GetName()));
1172 TCanvas* canvas2 = new TCanvas(Form("%s_2", GetName()), "Track Quality Results2", 1200, 800);
1173 canvas2->Divide(3, 2);
1176 fhC11[0]->SetStats(kFALSE);
1181 fhC22[0]->SetStats(kFALSE);
1186 fhC33[0]->SetStats(kFALSE);
1191 fhC44[0]->SetStats(kFALSE);
1196 fhC55[0]->SetStats(kFALSE);
1200 canvas2->SaveAs(Form("%s_%s.gif", GetName(), canvas2->GetName()));
1202 TCanvas* canvas3 = new TCanvas(Form("%s_3", GetName()), "Track Quality Results3", 1200, 800);
1203 canvas3->Divide(3, 2);
1206 fhDXY[0]->SetStats(kFALSE);
1211 fhDZ[0]->SetStats(kFALSE);
1216 fhDXYvsDZ[0]->SetStats(kFALSE);
1218 gPad->SetRightMargin(0.15);
1219 fhDXYvsDZ[0]->Draw("COLZ");
1222 fhDXYNormalized[0]->SetStats(kFALSE);
1224 fhDXYNormalized[0]->Draw();
1227 fhDZNormalized[0]->SetStats(kFALSE);
1229 fhDZNormalized[0]->Draw();
1232 fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1234 gPad->SetRightMargin(0.15);
1235 fhDXYvsDZNormalized[0]->Draw("COLZ");
1237 canvas3->SaveAs(Form("%s_%s.gif", GetName(), canvas3->GetName()));
1239 TCanvas* canvas4 = new TCanvas(Form("%s_4", GetName()), "Track Quality Results4", 800, 500);
1240 canvas4->Divide(2, 1);
1243 fhCutStatistics->SetStats(kFALSE);
1244 fhCutStatistics->LabelsOption("v");
1245 gPad->SetBottomMargin(0.3);
1246 fhCutStatistics->Draw();
1249 fhCutCorrelation->SetStats(kFALSE);
1250 fhCutCorrelation->LabelsOption("v");
1251 gPad->SetBottomMargin(0.3);
1252 gPad->SetLeftMargin(0.3);
1253 fhCutCorrelation->Draw("COLZ");
1255 canvas4->SaveAs(Form("%s_%s.gif", GetName(), canvas4->GetName()));
1258 fhDXYvsDZNormalized[0]->SetStats(kFALSE);
1259 fhDXYvsDZNormalized[0]->DrawCopy("COLZ");
1262 fhNClustersTPC[0]->SetStats(kFALSE);
1263 fhNClustersTPC[0]->DrawCopy();
1266 fhChi2PerClusterITS[0]->SetStats(kFALSE);
1267 fhChi2PerClusterITS[0]->DrawCopy();
1268 fhChi2PerClusterITS[1]->SetLineColor(2);
1269 fhChi2PerClusterITS[1]->DrawCopy("SAME");
1272 fhChi2PerClusterTPC[0]->SetStats(kFALSE);
1273 fhChi2PerClusterTPC[0]->DrawCopy();
1274 fhChi2PerClusterTPC[1]->SetLineColor(2);
1275 fhChi2PerClusterTPC[1]->DrawCopy("SAME");*/
1278 Float_t AliESDtrackCuts::GetMinNsigmaToVertex() const
1280 // deprecated, please use GetMaxNsigmaToVertex
1282 Printf("WARNING: AliESDtrackCuts::GetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use GetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1284 return GetMaxNsigmaToVertex();
1287 void AliESDtrackCuts::SetMinNsigmaToVertex(Float_t sigma)
1289 // deprecated, will be removed in next release
1291 Printf("WARNING: AliESDtrackCuts::SetMinNsigmaToVertex is DEPRECATED and will be removed in the next release. Please use SetMaxNsigmaToVertex instead. Renaming was done to improve code readability.");
1293 SetMaxNsigmaToVertex(sigma);
1296 void AliESDtrackCuts::SetAcceptKingDaughters(Bool_t b)
1298 // deprecated, will be removed in next release
1300 Printf("WARNING: AliESDtrackCuts::SetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use SetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
1302 SetAcceptKinkDaughters(b);
1305 Bool_t AliESDtrackCuts::GetAcceptKingDaughters() const
1307 // deprecated, will be removed in next release
1309 Printf("WARNING: AliESDtrackCuts::GetAcceptKingDaughters is DEPRECATED and will be removed in the next release. Please use GetAcceptKinkDaughters instead. Renaming was done to improve code readability.");
1311 return GetAcceptKinkDaughters();