1. Update to AliESDtrack
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Oct 2011 10:20:15 +0000 (10:20 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 18 Oct 2011 10:20:15 +0000 (10:20 +0000)
- Adding "golden TPC cut" on chi2 between tpc only track constrained
  to the vertex and the global track. The code is an adapted version
  from the implemention provided by the RAA group (Jacek et al).
- The calculated chi2 value is cached to avoid the computing-intensive
  operation several times in train running.
- New function that returns the number of crossed rows and caches it
  (suggestion by Jens W).

2. AliESDtrackCuts has been optimized
- Golden cut added. Therefore the prototype of AcceptTrack needed to
  be changed (backward compatible). One needs to pass the event as well
  when the new cut is to be used.
- CPU intensive cuts are only evaluated when other cuts have
  passed (Nsigma, golden TPC cut).
  NOTE: Control histograms are only filled for those quantities when the other cuts have passed

ANALYSIS/AliESDtrackCuts.cxx
ANALYSIS/AliESDtrackCuts.h
STEER/ESD/AliESDtrack.cxx
STEER/ESD/AliESDtrack.h
STEER/STEERBase/AliExternalTrackParam.cxx
STEER/STEERBase/AliExternalTrackParam.h

index 5a9b84a..4c458e4 100644 (file)
@@ -72,7 +72,8 @@ const Char_t* AliESDtrackCuts::fgkCutNames[kNCuts] = {
  "require ITS Pid",
  "n crossed rows TPC",
  "n crossed rows / n findable clusters",
- "missing ITS points"
+ "missing ITS points",
+ "#Chi^{2} TPC constrained vs. global"
 };
 
 AliESDtrackCuts* AliESDtrackCuts::fgMultEstTrackCuts[AliESDtrackCuts::kNMultEstTrackCuts] = { 0, 0, 0, 0 };
@@ -85,6 +86,8 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
   fCutMaxChi2PerClusterTPC(0),
   fCutMaxChi2PerClusterITS(0),
+  fCutMaxChi2TPCConstrainedVsGlobal(0),
+  fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
   fCutMaxMissingITSPoints(0),
   fCutMaxC11(0),
   fCutMaxC22(0),
@@ -148,7 +151,9 @@ AliESDtrackCuts::AliESDtrackCuts(const Char_t* name, const Char_t* title) : AliA
   SetMinNCrossedRowsTPC();
   SetMinRatioCrossedRowsOverFindableClustersTPC();
   SetMaxChi2PerClusterTPC();
-  SetMaxChi2PerClusterITS();                               
+  SetMaxChi2PerClusterITS();
+  SetMaxChi2TPCConstrainedGlobal();
+  SetMaxChi2TPCConstrainedGlobalVertexType();
   SetMaxNOfMissingITSPoints();
   SetMaxCovDiagonalElements();
   SetMaxRel1PtUncertainty();
@@ -189,6 +194,8 @@ AliESDtrackCuts::AliESDtrackCuts(const AliESDtrackCuts &c) : AliAnalysisCuts(c),
   fCutMinRatioCrossedRowsOverFindableClustersTPC(0),
   fCutMaxChi2PerClusterTPC(0),
   fCutMaxChi2PerClusterITS(0),
+  fCutMaxChi2TPCConstrainedVsGlobal(0),
+  fCutMaxChi2TPCConstrainedVsGlobalVertexType(kVertexTracks | kVertexSPD),
   fCutMaxMissingITSPoints(0),
   fCutMaxC11(0),
   fCutMaxC22(0),
@@ -268,6 +275,8 @@ AliESDtrackCuts::~AliESDtrackCuts()
       delete fhChi2PerClusterITS[i];       
     if (fhChi2PerClusterTPC[i])
       delete fhChi2PerClusterTPC[i];    
+    if (fhChi2TPCConstrainedVsGlobal[i])
+      delete fhChi2TPCConstrainedVsGlobal[i];
     if(fhNClustersForITSPID[i])
       delete fhNClustersForITSPID[i];
     if(fhNMissingITSPoints[i])
@@ -339,6 +348,8 @@ void AliESDtrackCuts::Init()
 
   fCutMaxChi2PerClusterTPC = 0;
   fCutMaxChi2PerClusterITS = 0;
+  fCutMaxChi2TPCConstrainedVsGlobal = 0;
+  fCutMaxChi2TPCConstrainedVsGlobalVertexType = kVertexTracks | kVertexSPD;
   fCutMaxMissingITSPoints  = 0;
   
   for (Int_t i = 0; i < 3; i++)
@@ -413,6 +424,7 @@ void AliESDtrackCuts::Init()
 
     fhChi2PerClusterITS[i] = 0;
     fhChi2PerClusterTPC[i] = 0;
+    fhChi2TPCConstrainedVsGlobal[i] = 0;
     fhNClustersForITSPID[i] = 0;
     fhNMissingITSPoints[i] = 0;
 
@@ -473,6 +485,8 @@ void AliESDtrackCuts::Copy(TObject &c) const
 
   target.fCutMaxChi2PerClusterTPC = fCutMaxChi2PerClusterTPC;
   target.fCutMaxChi2PerClusterITS = fCutMaxChi2PerClusterITS;
+  target.fCutMaxChi2TPCConstrainedVsGlobal = fCutMaxChi2TPCConstrainedVsGlobal;
+  target.fCutMaxChi2TPCConstrainedVsGlobalVertexType = fCutMaxChi2TPCConstrainedVsGlobalVertexType;
   target.fCutMaxMissingITSPoints = fCutMaxMissingITSPoints;
 
   for (Int_t i = 0; i < 3; i++)
@@ -543,6 +557,7 @@ void AliESDtrackCuts::Copy(TObject &c) const
 
     if (fhChi2PerClusterITS[i]) target.fhChi2PerClusterITS[i] = (TH1F*) fhChi2PerClusterITS[i]->Clone();
     if (fhChi2PerClusterTPC[i]) target.fhChi2PerClusterTPC[i] = (TH1F*) fhChi2PerClusterTPC[i]->Clone();
+    if (fhChi2TPCConstrainedVsGlobal[i]) target.fhChi2TPCConstrainedVsGlobal[i] = (TH1F*) fhChi2TPCConstrainedVsGlobal[i]->Clone();
     if (fhNClustersForITSPID[i]) target.fhNClustersForITSPID[i] = (TH1F*) fhNClustersForITSPID[i]->Clone();
     if (fhNMissingITSPoints[i]) target.fhNMissingITSPoints[i] = (TH1F*) fhNMissingITSPoints[i]->Clone();
 
@@ -612,6 +627,8 @@ Long64_t AliESDtrackCuts::Merge(TCollection* list) {
                                                                    
       fhChi2PerClusterITS[i] ->Add(entry->fhChi2PerClusterITS[i]); 
       fhChi2PerClusterTPC[i] ->Add(entry->fhChi2PerClusterTPC[i]); 
+      if (fhChi2TPCConstrainedVsGlobal[i])
+       fhChi2TPCConstrainedVsGlobal[i]->Add(entry->fhChi2TPCConstrainedVsGlobal[i]);
       if (fhNClustersForITSPID[i])
        fhNClustersForITSPID[i]->Add(entry->fhNClustersForITSPID[i]);
       if (fhNMissingITSPoints[i])
@@ -887,7 +904,7 @@ Int_t AliESDtrackCuts::GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t t
 }
 
 //____________________________________________________________________
-Float_t AliESDtrackCuts::GetSigmaToVertex(AliESDtrack* const esdTrack)
+Float_t AliESDtrackCuts::GetSigmaToVertex(const AliESDtrack* const esdTrack)
 {
   // Calculates the number of sigma to the vertex.
 
@@ -950,7 +967,7 @@ void AliESDtrackCuts::EnableNeededBranches(TTree* tree)
 }
 
 //____________________________________________________________________
-Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack) 
+Bool_t AliESDtrackCuts::AcceptTrack(const AliESDtrack* esdTrack, const AliESDEvent* esdEvent) 
 {
   // 
   // figure out if the tracks survives all the track cuts defined
@@ -975,6 +992,8 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
   // fTracks.fR   //GetMass
   // fTracks.fP   //GetMass
   // fTracks.fKinkIndexes
+  //
+  // esdEvent is only required for the MaxChi2TPCConstrainedVsGlobal
 
   UInt_t status = esdTrack->GetStatus();
 
@@ -987,7 +1006,7 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
   else {
     nClustersTPC = esdTrack->GetTPCclusters(0);
   }
-  Float_t nCrossedRowsTPC = esdTrack->GetTPCClusterInfo(2,1); 
+  Float_t nCrossedRowsTPC = esdTrack->GetTPCCrossedRows();
   Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
   if (esdTrack->GetTPCNclsF()>0) {
     ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC / esdTrack->GetTPCNclsF();
@@ -1012,9 +1031,6 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
   Double_t extCov[15];
   esdTrack->GetExternalCovariance(extCov);
 
-  // getting the track to vertex parameters
-  Float_t nSigmaToVertex = GetSigmaToVertex(esdTrack);
-      
   Float_t b[2];
   Float_t bCov[3];
   esdTrack->GetImpactParameters(b,bCov);
@@ -1095,11 +1111,9 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
     cuts[10]=kTRUE;  
   if (extCov[14]  > fCutMaxC55) 
     cuts[11]=kTRUE;  
-  if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
-    cuts[12] = kTRUE;
-  // if n sigma could not be calculated
-  if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
-    cuts[13]=kTRUE;
+  
+  // cut 12 and 13 see below
+  
   if (!fCutAcceptKinkDaughters && esdTrack->GetKinkIndex(0)>0)
     cuts[14]=kTRUE;
   // track kinematics cut
@@ -1177,12 +1191,61 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
     if(retc && statusLay==5) ++nMissITSpts;
   }
   if(nMissITSpts>fCutMaxMissingITSPoints) cuts[38] = kTRUE;
-
+  
   Bool_t cut=kFALSE;
   for (Int_t i=0; i<kNCuts; i++) 
     if (cuts[i]) {cut = kTRUE;}
 
-
+  // for performance evaluate the CPU intensive cuts only when the others have passed, and when they are requested
+  Double_t chi2TPCConstrainedVsGlobal = -2;
+  Float_t nSigmaToVertex = -2;
+  if (!cut)
+  {
+    // getting the track to vertex parameters
+    if (fCutSigmaToVertexRequired)
+    {
+      nSigmaToVertex = GetSigmaToVertex(esdTrack);
+      if (nSigmaToVertex > fCutNsigmaToVertex && fCutSigmaToVertexRequired)
+      {
+       cuts[12] = kTRUE;
+       cut = kTRUE;
+      }
+      // if n sigma could not be calculated
+      if (nSigmaToVertex<0 && fCutSigmaToVertexRequired)
+      {
+       cuts[13] = kTRUE;
+       cut = kTRUE;
+      }
+    }
+      
+    // max chi2 TPC constrained vs global track only if track passed the other cut
+    if (fCutMaxChi2TPCConstrainedVsGlobal < 1e9)
+    {
+      if (!esdEvent)
+       AliFatal("fCutMaxChi2TPCConstrainedVsGlobal set but ESD event not provided.");
+      
+      // get vertex
+      const AliESDVertex* vertex = 0;
+      if (fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTracks)
+       vertex = esdEvent->GetPrimaryVertexTracks();
+      
+      if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexSPD)
+       vertex = esdEvent->GetPrimaryVertexSPD();
+       
+      if ((!vertex || !vertex->GetStatus()) && fCutMaxChi2TPCConstrainedVsGlobalVertexType & kVertexTPC)
+       vertex = esdEvent->GetPrimaryVertexTPC();
+
+      if (vertex->GetStatus())
+       chi2TPCConstrainedVsGlobal = esdTrack->GetChi2TPCConstrainedVsGlobal(vertex);
+      
+      if (chi2TPCConstrainedVsGlobal < 0 || chi2TPCConstrainedVsGlobal > fCutMaxChi2TPCConstrainedVsGlobal)
+      {
+       cuts[39] = kTRUE;
+       cut = kTRUE;
+      }
+    }
+  }
+  
   //########################################################################
   // filling histograms
   if (fHistogramsOn) {
@@ -1224,6 +1287,7 @@ Bool_t AliESDtrackCuts::AcceptTrack(AliESDtrack* esdTrack)
       fhRatioCrossedRowsOverFindableClustersTPC[id]->Fill(ratioCrossedRowsOverFindableClustersTPC);
       fhChi2PerClusterITS[id]->Fill(chi2PerClusterITS);
       fhChi2PerClusterTPC[id]->Fill(chi2PerClusterTPC);
+      fhChi2TPCConstrainedVsGlobal[id]->Fill(chi2TPCConstrainedVsGlobal);
       fhNClustersForITSPID[id]->Fill(nITSPointsForPid);
       fhNMissingITSPoints[id]->Fill(nMissITSpts);
 
@@ -1286,9 +1350,7 @@ Bool_t AliESDtrackCuts::CheckITSClusterRequirement(ITSClusterRequirement req, Bo
 //____________________________________________________________________
 AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTrack)
 {
-  
-  // Utility function to 
-  // create a TPC only track from the given esd track
+  // Utility function to create a TPC only track from the given esd track
   // 
   // IMPORTANT: The track has to be deleted by the user
   //
@@ -1317,16 +1379,11 @@ AliESDtrack* AliESDtrackCuts::GetTPCOnlyTrack(const AliESDEvent* esd, Int_t iTra
     return 0;
   }
 
-  // propagate to Vertex
-  // not needed for normal reconstructed ESDs...
-  // Double_t pTPC[2],covTPC[3];
-  // tpcTrack->PropagateToDCA(esd->GetPrimaryVertexTPC(), esd->GetMagneticField(), 10000,  pTPC, covTPC);
-
   return tpcTrack;
 }
 
 //____________________________________________________________________
-TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
+TObjArray* AliESDtrackCuts::GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC)
 {
   //
   // returns an array of all tracks that pass the cuts
@@ -1347,7 +1404,7 @@ TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
       if (!tpcTrack)
         continue;
 
-      if (AcceptTrack(tpcTrack)) {
+      if (AcceptTrack(tpcTrack, esd)) {
         acceptedTracks->Add(tpcTrack);
       }
       else
@@ -1356,7 +1413,7 @@ TObjArray* AliESDtrackCuts::GetAcceptedTracks(AliESDEvent* esd,Bool_t bTPC)
     else
     {
       AliESDtrack* track = esd->GetTrack(iTrack);
-      if(AcceptTrack(track))
+      if(AcceptTrack(track, esd))
         acceptedTracks->Add(track);
     }
   } 
@@ -1376,7 +1433,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
   // loop over esd tracks
   for (Int_t iTrack = 0; iTrack < esd->GetNumberOfTracks(); iTrack++) {
     AliESDtrack* track = esd->GetTrack(iTrack);
-    if (AcceptTrack(track))
+    if (AcceptTrack(track, esd))
       count++;
   }
 
@@ -1423,6 +1480,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
     fhRatioCrossedRowsOverFindableClustersTPC[i]     = new TH1F("ratioCrossedRowsOverFindableClustersTPC"  ,"",60,0,1.5);
     fhChi2PerClusterITS[i]   = new TH1F("chi2PerClusterITS","",500,0,10);
     fhChi2PerClusterTPC[i]   = new TH1F("chi2PerClusterTPC","",500,0,10);
+    fhChi2TPCConstrainedVsGlobal[i]   = new TH1F("chi2TPCConstrainedVsGlobal","",600,-2,50);
     fhNClustersForITSPID[i]  = new TH1F("nPointsForITSpid","",5,-0.5,4.5);
     fhNMissingITSPoints[i]   = new TH1F("nMissingITSClusters","",7,-0.5,6.5);
 
@@ -1453,6 +1511,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
     fhNSharedClustersTPC[i]->SetTitle("n TPC shared clusters");
     fhChi2PerClusterITS[i]->SetTitle("#Chi^{2} per ITS cluster");
     fhChi2PerClusterTPC[i]->SetTitle("#Chi^{2} per TPC cluster");
+    fhChi2TPCConstrainedVsGlobal[i]->SetTitle("#Chi^{2} TPC constrained track vs global track");
     fhNClustersForITSPID[i]->SetTitle("n ITS points for PID");
     fhNMissingITSPoints[i]->SetTitle("n ITS layers with missing cluster");
 
@@ -1481,6 +1540,7 @@ Int_t AliESDtrackCuts::CountAcceptedTracks(const AliESDEvent* const esd)
     fhNSharedClustersTPC[i]->SetLineColor(color);  fhNSharedClustersTPC[i]->SetLineWidth(2);
     fhChi2PerClusterITS[i]->SetLineColor(color);   fhChi2PerClusterITS[i]->SetLineWidth(2);
     fhChi2PerClusterTPC[i]->SetLineColor(color);   fhChi2PerClusterTPC[i]->SetLineWidth(2);
+    fhChi2TPCConstrainedVsGlobal[i]->SetLineColor(color);   fhChi2TPCConstrainedVsGlobal[i]->SetLineWidth(2);
     fhNClustersForITSPID[i]->SetLineColor(color);  fhNClustersForITSPID[i]->SetLineWidth(2);
     fhNMissingITSPoints[i]->SetLineColor(color);   fhNMissingITSPoints[i]->SetLineWidth(2);
 
@@ -1542,6 +1602,7 @@ Bool_t AliESDtrackCuts::LoadHistograms(const Char_t* dir)
     fhRatioCrossedRowsOverFindableClustersTPC[i]   = dynamic_cast<TH1F*> (gDirectory->Get("ratioCrossedRowsOverFindableClustersTPC"  ));
     fhChi2PerClusterITS[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterITS"));
     fhChi2PerClusterTPC[i] = dynamic_cast<TH1F*> (gDirectory->Get("chi2PerClusterTPC"));
+    fhChi2TPCConstrainedVsGlobal[i] = dynamic_cast<TH1F*> (gDirectory->Get("fhChi2TPCConstrainedVsGlobal"));
     fhNClustersForITSPID[i] = dynamic_cast<TH1F*> (gDirectory->Get("nPointsForITSpid"));
     fhNMissingITSPoints[i] = dynamic_cast<TH1F*> (gDirectory->Get("nMissingITSClusters"));
 
@@ -1614,6 +1675,7 @@ void AliESDtrackCuts::SaveHistograms(const Char_t* dir) {
     fhRatioCrossedRowsOverFindableClustersTPC[i]     ->Write();
     fhChi2PerClusterITS[i]   ->Write();
     fhChi2PerClusterTPC[i]   ->Write();
+    fhChi2TPCConstrainedVsGlobal[i]   ->Write();
     fhNClustersForITSPID[i]  ->Write();
     fhNMissingITSPoints[i]   ->Write();
 
index c6cce5b..8e8ab9e 100644 (file)
@@ -43,6 +43,7 @@ public:
   enum Detector { kSPD = 0, kSDD, kSSD };
   enum MultEstTrackCuts { kMultEstTrackCutGlobal = 0, kMultEstTrackCutITSSA, kMultEstTrackCutDCAwSPD, kMultEstTrackCutDCAwoSPD, kNMultEstTrackCuts /* this must always be the last */};
   enum MultEstTrackType { kTrackletsITSTPC = 0, kTrackletsITSSA, kTracklets };
+  enum VertexType { kVertexTracks = 0x1, kVertexSPD = 0x2, kVertexTPC = 0x4 };
   
   AliESDtrackCuts(const Char_t* name = "AliESDtrackCuts", const Char_t* title = "");
   virtual ~AliESDtrackCuts();
@@ -51,8 +52,8 @@ public:
        {return AcceptTrack((AliESDtrack*)obj);}
   virtual Bool_t IsSelected(TList* /*list*/) {return kTRUE;}
 
-  Bool_t AcceptTrack(AliESDtrack* esdTrack);
-  TObjArray* GetAcceptedTracks(AliESDEvent* esd, Bool_t bTPC = kFALSE);
+  Bool_t AcceptTrack(const AliESDtrack* esdTrack, const AliESDEvent* esdEvent = 0);
+  TObjArray* GetAcceptedTracks(const AliESDEvent* esd, Bool_t bTPC = kFALSE);
   Int_t CountAcceptedTracks(const AliESDEvent* const esd);
   
   static Int_t GetReferenceMultiplicity(const AliESDEvent* esd, Bool_t tpcOnly);
@@ -87,6 +88,8 @@ public:
   void SetClusterRequirementITS(Detector det, ITSClusterRequirement req = kOff) { fCutClusterRequirementITS[det] = req; }
   void SetMaxChi2PerClusterTPC(Float_t max=1e10) {fCutMaxChi2PerClusterTPC=max;}
   void SetMaxChi2PerClusterITS(Float_t max=1e10) {fCutMaxChi2PerClusterITS=max;}
+  void SetMaxChi2TPCConstrainedGlobal(Float_t max=1e10) {fCutMaxChi2TPCConstrainedVsGlobal = max; }
+  void SetMaxChi2TPCConstrainedGlobalVertexType(Int_t vertexType = kVertexTracks | kVertexSPD) { fCutMaxChi2TPCConstrainedVsGlobalVertexType = vertexType; }
   void SetMaxNOfMissingITSPoints(Int_t max=6)    {fCutMaxMissingITSPoints=max;}
   void SetRequireTPCRefit(Bool_t b=kFALSE)       {fCutRequireTPCRefit=b;}
   void SetRequireTPCStandAlone(Bool_t b=kFALSE)  {fCutRequireTPCStandAlone=b;}
@@ -125,6 +128,8 @@ public:
   ITSClusterRequirement GetClusterRequirementITS(Detector det) const { return fCutClusterRequirementITS[det]; }
   Float_t GetMaxChi2PerClusterTPC()  const   { return fCutMaxChi2PerClusterTPC;}
   Float_t GetMaxChi2PerClusterITS()  const   { return fCutMaxChi2PerClusterITS;}
+  Float_t GetMaxChi2TPCConstrainedGlobal() const { return fCutMaxChi2TPCConstrainedVsGlobal; }
+  Int_t GetMaxChi2TPCConstrainedGlobalVertexType() const { return fCutMaxChi2TPCConstrainedVsGlobalVertexType; }
   Int_t   GetMaxNOfMissingITSPoints() const   { return  fCutMaxMissingITSPoints;}
   Bool_t  GetRequireTPCRefit()       const   { return fCutRequireTPCRefit;}
   Bool_t  GetRequireTPCStandAlone()  const   { return fCutRequireTPCStandAlone;}
@@ -172,7 +177,7 @@ public:
   void SaveHistograms(const Char_t* dir = 0);
   void DrawHistograms();
 
-  static Float_t GetSigmaToVertex(AliESDtrack* const esdTrack);
+  static Float_t GetSigmaToVertex(const AliESDtrack* const esdTrack);
   
   static void EnableNeededBranches(TTree* tree);
 
@@ -189,7 +194,7 @@ protected:
   Bool_t CheckPtDepDCA(TString dist,Bool_t print=kFALSE) const;
   void SetPtDepDCACuts(Double_t pt);
 
-  enum { kNCuts = 39 }; 
+  enum { kNCuts = 40 }; 
 
   //######################################################
   // esd track quality cuts
@@ -205,6 +210,8 @@ protected:
 
   Float_t fCutMaxChi2PerClusterTPC;   // max tpc fit chi2 per tpc cluster
   Float_t fCutMaxChi2PerClusterITS;   // max its fit chi2 per its cluster
+  Float_t fCutMaxChi2TPCConstrainedVsGlobal; // max chi2 TPC track constrained with vtx vs. global track
+  Int_t fCutMaxChi2TPCConstrainedVsGlobalVertexType; // vertex type for max chi2 TPC track constrained with vtx vs. global track (can be configured to accept several vertex types)
   Int_t   fCutMaxMissingITSPoints;    // max n. of missing ITS points
 
   Float_t fCutMaxC11;                 // max cov. matrix diag. elements (res. y^2)
@@ -268,6 +275,7 @@ protected:
 
   TH1F* fhChi2PerClusterITS[2];       //->
   TH1F* fhChi2PerClusterTPC[2];       //->
+  TH1F* fhChi2TPCConstrainedVsGlobal[2];       //->
   TH1F* fhNClustersForITSPID[2];      //-> number of points in SDD+SSD (ITS PID selection)
   TH1F* fhNMissingITSPoints[2];       //-> number of missing ITS points
 
@@ -297,7 +305,7 @@ protected:
   TH1F*  fhCutStatistics;             //-> statistics of what cuts the tracks did not survive
   TH2F*  fhCutCorrelation;            //-> 2d statistics plot
 
-  ClassDef(AliESDtrackCuts, 16)
+  ClassDef(AliESDtrackCuts, 17)
 };
 
 
index dceb3e9..4f0ce25 100644 (file)
 #include <TMath.h>
 #include <TParticle.h>
 #include <TDatabasePDG.h>
+#include <TMatrixD.h>
 
 #include "AliESDVertex.h"
 #include "AliESDtrack.h"
 #include "AliLog.h"
 #include "AliTrackPointArray.h"
 #include "TPolyMarker3D.h"
+#include "AliTrackerBase.h"
 
 ClassImp(AliESDtrack)
 
@@ -227,7 +229,10 @@ AliESDtrack::AliESDtrack() :
   fTRDnSlices(0),
   fTRDslices(0x0),
   fVertexID(-2),// -2 means an orphan track 
-  fESDEvent(0)
+  fESDEvent(0),
+  fCacheNCrossedRows(-10),
+  fCacheChi2TPCConstrainedVsGlobal(-10),
+  fCacheChi2TPCConstrainedVsGlobalVertex(0)
 {
   //
   // The default ESD constructor 
@@ -332,7 +337,10 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fTRDnSlices(track.fTRDnSlices),
   fTRDslices(0x0),
   fVertexID(track.fVertexID),
-  fESDEvent(track.fESDEvent)
+  fESDEvent(track.fESDEvent),
+  fCacheNCrossedRows(track.fCacheNCrossedRows),
+  fCacheChi2TPCConstrainedVsGlobal(track.fCacheChi2TPCConstrainedVsGlobal),
+  fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex)
 {
   //
   //copy constructor
@@ -446,7 +454,10 @@ AliESDtrack::AliESDtrack(const AliVTrack *track) :
   fTRDnSlices(0),
   fTRDslices(0x0),
   fVertexID(-2),  // -2 means an orphan track
-  fESDEvent(0)  
+  fESDEvent(0),
+  fCacheNCrossedRows(-10),
+  fCacheChi2TPCConstrainedVsGlobal(-10),
+  fCacheChi2TPCConstrainedVsGlobalVertex(0)
 {
   //
   // ESD track from AliVTrack.
@@ -583,7 +594,10 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fTRDnSlices(0),
   fTRDslices(0x0),
   fVertexID(-2),  // -2 means an orphan track
-  fESDEvent(0)
+  fESDEvent(0),
+  fCacheNCrossedRows(-10),
+  fCacheChi2TPCConstrainedVsGlobal(-10),
+  fCacheChi2TPCConstrainedVsGlobalVertex(0)  
 {
   //
   // ESD track from TParticle
@@ -1709,6 +1723,19 @@ UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 }
 
 //_______________________________________________________________________
+Float_t AliESDtrack::GetTPCCrossedRows() const
+{
+  // This function calls GetTPCClusterInfo with some default parameters which are used in the track selection and caches the outcome
+  // because GetTPCClusterInfo is quite time-consuming
+  
+  if (fCacheNCrossedRows > -1)
+    return fCacheNCrossedRows;
+  
+  fCacheNCrossedRows = GetTPCClusterInfo(2, 1);
+  return fCacheNCrossedRows;
+}
+
+//_______________________________________________________________________
 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
 {
   //
@@ -2386,3 +2413,90 @@ UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
   //
   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
 }
+
+//____________________________________________________________________
+Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const
+{
+  // Calculates the chi2 between the TPC track (TPCinner) constrained to the primary vertex and the global track
+  //
+  // Returns -1 in case the calculation failed
+  //
+  // Value is cached as a non-persistent member.
+  //
+  // Code adapted from original code by GSI group (Jacek, Marian, Michael)
+  
+  // cache, ignoring that a different vertex might be passed
+  if (fCacheChi2TPCConstrainedVsGlobalVertex == vtx)
+    return fCacheChi2TPCConstrainedVsGlobal;
+  
+  fCacheChi2TPCConstrainedVsGlobal = -1;
+  fCacheChi2TPCConstrainedVsGlobalVertex = vtx;
+  
+  Double_t x[3];
+  GetXYZ(x);
+  Double_t b[3];
+  AliTrackerBase::GetBxByBz(x,b);
+
+  if (!fTPCInner)  { 
+    AliWarning("Could not get TPC Inner Param.");
+    return fCacheChi2TPCConstrainedVsGlobal;
+  }
+  
+  // clone for constraining
+  AliExternalTrackParam* tpcInnerC = new AliExternalTrackParam(*fTPCInner);
+  if (!tpcInnerC) { 
+    AliWarning("Clone of TPCInnerParam failed.");
+    return fCacheChi2TPCConstrainedVsGlobal;  
+  }
+  
+  // transform to the track reference frame 
+  Bool_t isOK = tpcInnerC->Rotate(GetAlpha());
+  isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
+  if (!isOK) { 
+    delete tpcInnerC;
+    tpcInnerC = 0; 
+    AliWarning("Rotation/Propagation of track failed.") ; 
+    return fCacheChi2TPCConstrainedVsGlobal;    
+  }  
+
+  // constrain TPCinner 
+  isOK = tpcInnerC->ConstrainToVertex(vtx, b);
+  
+  // transform to the track reference frame 
+  isOK &= tpcInnerC->Rotate(GetAlpha());
+  isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
+
+  if (!isOK) {
+    AliWarning("ConstrainTPCInner failed.") ;
+    delete tpcInnerC;
+    tpcInnerC = 0; 
+    return fCacheChi2TPCConstrainedVsGlobal;  
+  }
+  
+  // calculate chi2 between vi and vj vectors
+  // with covi and covj covariance matrices
+  // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
+  TMatrixD deltaT(5,1);
+  TMatrixD delta(1,5);
+  TMatrixD covarM(5,5);
+
+  for (Int_t ipar=0; ipar<5; ipar++) {
+    deltaT(ipar,0) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
+    delta(0,ipar) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
+
+    for (Int_t jpar=0; jpar<5; jpar++) {
+      Int_t index = GetIndex(ipar,jpar);
+      covarM(ipar,jpar) = GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
+    }
+  }
+  // chi2 distance TPC constrained and TPC+ITS
+  TMatrixD covarMInv = covarM.Invert();
+  TMatrixD mat2 = covarMInv*deltaT;
+  TMatrixD chi2 = delta*mat2; 
+  
+  delete tpcInnerC; 
+  tpcInnerC = 0;
+  
+  fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
+  return fCacheChi2TPCConstrainedVsGlobal;
+}
index d815b34..784e8e7 100644 (file)
@@ -102,6 +102,7 @@ public:
               (Double_t &alpha, Double_t &x, Double_t p[5]) const;
   Bool_t GetConstrainedExternalCovariance(Double_t cov[15]) const;
   Double_t GetConstrainedChi2() const {return fCchi2;}
+  Double_t GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const;
   //
   
   // global track chi2
@@ -231,6 +232,7 @@ public:
   void    SetTPCClusterMap(const TBits amap) {fTPCClusterMap = amap;}
   void    SetTPCSharedMap(const TBits amap) {fTPCSharedMap = amap;}
   Float_t GetTPCClusterInfo(Int_t nNeighbours=3, Int_t type=0, Int_t row0=0, Int_t row1=159) const;
+  Float_t GetTPCCrossedRows() const;
   
   void    SetTRDpid(const Double_t *p);
   void    SetTRDsignal(Double_t sig) {fTRDsignal = sig;}
@@ -499,11 +501,15 @@ protected:
   Char_t  fVertexID; // ID of the primary vertex this track belongs to
   AliESDEvent*   fESDEvent; //!Pointer back to event to which the track belongs
   
+  mutable Float_t fCacheNCrossedRows; //! Cache for the number of crossed rows
+  mutable Float_t fCacheChi2TPCConstrainedVsGlobal; //! Cache for the chi2 of constrained TPC vs global track
+  mutable const AliESDVertex* fCacheChi2TPCConstrainedVsGlobalVertex; //! Vertex for which the cache is valid
+  
  private:
   static bool fgkOnlineMode; //! indicate the online mode to skip some of the functionality
 
   AliESDtrack & operator=(const AliESDtrack & );
-  ClassDef(AliESDtrack,60)  //ESDtrack 
+  ClassDef(AliESDtrack,61)  //ESDtrack 
 };
 
 
index 77c6d3d..d1f41f7 100644 (file)
@@ -2248,3 +2248,30 @@ void AliExternalTrackParam::CheckCovariance() {
 //       eig.Print();
 //     }
 }
+
+Bool_t AliExternalTrackParam::ConstrainToVertex(const AliVVertex* vtx, Double_t b[3])
+{
+  // Constrain TPC inner params constrained
+  //
+  if (!vtx) 
+    return kFALSE;
+
+  Double_t dz[2], cov[3];
+  if (!PropagateToDCABxByBz(vtx, b, 3, dz, cov)) 
+    return kFALSE; 
+
+  Double_t covar[6]; 
+  vtx->GetCovarianceMatrix(covar);
+  
+  Double_t p[2]= { fP[0] - dz[0], fP[1] - dz[1] };
+  Double_t c[3]= { covar[2], 0., covar[5] };
+  
+  Double_t chi2C = GetPredictedChi2(p,c);
+  if (chi2C>kVeryBig) 
+    return kFALSE; 
+
+  if (!Update(p,c)) 
+    return kFALSE; 
+
+  return kTRUE;
+}
index d6a7f19..ce4b033 100644 (file)
@@ -200,6 +200,7 @@ class AliExternalTrackParam: public AliVTrack {
                         Double_t dz[2]=0, Double_t cov[3]=0);
   Bool_t PropagateToDCABxByBz(const AliVVertex *vtx, Double_t b[3], 
          Double_t maxd, Double_t dz[2]=0, Double_t cov[3]=0);
+  Bool_t ConstrainToVertex(const AliVVertex* vtx, Double_t b[3]);
   
   void GetDirection(Double_t d[3]) const;
   Bool_t GetPxPyPz(Double_t *p) const;