]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDgtuTrack.cxx
Introduce an enumerator for PID methods
[u/mrichter/AliRoot.git] / TRD / AliTRDgtuTrack.cxx
index a4b50bfba43073459801f3119e31f8b1ef387d3d..35390681e8111dce8353b9c5b39320cdad27e03b 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 ///////////////////////////////////////////////////////////////////////////////
 //                                                                           //
 //                                                                           //
 //                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
+#include <TMath.h>
 #include <TMatrixD.h>
 #include <TObjArray.h>
 
+#include "AliLog.h"
+
+#include "AliTRDReconstructor.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDcalibDB.h"
-#include "Cal/AliTRDCalPIDLQ.h"
-
 #include "AliTRDltuTracklet.h"
-
 #include "AliTRDgtuTrack.h"
+#include "Cal/AliTRDCalPID.h"
 
 ClassImp(AliTRDgtuTrack)
 
 //_____________________________________________________________________________
 AliTRDgtuTrack::AliTRDgtuTrack()
+  :TObject()
+  ,fTracklets(new TObjArray(400))
+  ,fYproj(0)
+  ,fZproj(0)
+  ,fSlope(0)
+  ,fDetector(-1)
+  ,fNtracklets(0)
+  ,fNplanes(0)
+  ,fNclusters(0)
+  ,fPt(0)
+  ,fPhi(0)
+  ,fEta(0)
+  ,fLabel(-1)
+  ,fPID(0)
+  ,fIsElectron(kFALSE)
 {
   //
   // Default constructor
   //
 
-  fYproj = 0.0;
-  fZproj = 0.0;
-  fSlope = 0.0;
-  fDetector = -1;
-  fNtracklets = 0;
-  fNplanes = 0;
-  fNclusters = 0;
-  fPt = 0.0;
-  fPhi = 0.0;
-  fEta = 0.0;
-  fLabel = -1;
-  fPID = 0.0;
-  fIsElectron = kFALSE;
-
-  fTracklets = new TObjArray(400);
-
 }
 
 //_____________________________________________________________________________
-AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack& track):
-  TObject(track),
-  fTracklets(NULL),
-  fYproj(track.fYproj),
-  fZproj(track.fZproj),
-  fSlope(track.fSlope),
-  fDetector(track.fDetector),
-  fNtracklets(track.fNtracklets),
-  fNplanes(track.fNplanes),
-  fNclusters(track.fNclusters),
-  fPt(track.fPt),
-  fPhi(track.fPhi),
-  fEta(track.fEta),
-  fLabel(track.fLabel),
-  fPID(track.fPID),
-  fIsElectron(track.fIsElectron)
+AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack &t)
+  :TObject(t)
+  ,fTracklets(NULL)
+  ,fYproj(t.fYproj)
+  ,fZproj(t.fZproj)
+  ,fSlope(t.fSlope)
+  ,fDetector(t.fDetector)
+  ,fNtracklets(t.fNtracklets)
+  ,fNplanes(t.fNplanes)
+  ,fNclusters(t.fNclusters)
+  ,fPt(t.fPt)
+  ,fPhi(t.fPhi)
+  ,fEta(t.fEta)
+  ,fLabel(t.fLabel)
+  ,fPID(t.fPID)
+  ,fIsElectron(t.fIsElectron)
 {
   //
-  // copy contructor
+  // Copy contructor
   //
 
 }
@@ -87,19 +90,20 @@ AliTRDgtuTrack::AliTRDgtuTrack(const AliTRDgtuTrack& track):
 AliTRDgtuTrack &AliTRDgtuTrack::operator=(const AliTRDgtuTrack &t)
 {
   //
-  // assignment operator
+  // Assignment operator
   //
 
   if (this != &t) ((AliTRDgtuTrack &) t).Copy(*this); 
+
   return *this;
 
-}
+} 
 
 //_____________________________________________________________________________
 void AliTRDgtuTrack::Copy(TObject &t) const
 {
   //
-  // copy function
+  // Copy function
   //
 
   ((AliTRDgtuTrack &) t).fTracklets  = NULL;
@@ -123,9 +127,16 @@ void AliTRDgtuTrack::Copy(TObject &t) const
 AliTRDgtuTrack::~AliTRDgtuTrack()
 {
   //
-  // destructor
+  // Destructor
   //
 
+  if (fTracklets) {
+    //fTracklets->Delete();
+    fTracklets->Clear();
+    delete fTracklets;
+    fTracklets = 0;
+  }
+
 }
 
 //_____________________________________________________________________________
@@ -140,31 +151,39 @@ void AliTRDgtuTrack::AddTracklet(AliTRDltuTracklet *trk)
 }
 
 //_____________________________________________________________________________
-AliTRDltuTrackletAliTRDgtuTrack::GetTracklet(Int_t pos) const
+AliTRDltuTracklet *AliTRDgtuTrack::GetTracklet(Int_t pos) const
 {
   //
   // Return LTU tracklet at position "pos"
   //
 
-  if (fTracklets == 0) return 0;
-  void * trk = fTracklets->UncheckedAt(pos);
-  if (trk == 0) return 0;
+  if (fTracklets == 0) {
+    return 0;
+  }
+  void *trk = fTracklets->UncheckedAt(pos);
+  if (trk == 0) {
+    return 0;
+  }
 
-  return (AliTRDltuTracklet*)trk;
+  return (AliTRDltuTracklet *) trk;
 
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDgtuTrack::Compare(const TObject * o) const
+Int_t AliTRDgtuTrack::Compare(const TObject *o) const
 {
   //
   // Compare function for sorting the tracks
   //
 
-  AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack*)o;
+  AliTRDgtuTrack *gtutrack = (AliTRDgtuTrack *) o;
 
-  if (fYproj <  gtutrack->GetYproj()) return -1;
-  if (fYproj == gtutrack->GetYproj()) return  0;
+  if (fYproj <  gtutrack->GetYproj()) {
+    return -1;
+  }
+  if (fYproj == gtutrack->GetYproj()) {
+    return  0;
+  }
 
   return +1;
 
@@ -177,18 +196,18 @@ void AliTRDgtuTrack::Reset()
   // Reset the track information
   //
 
-  fYproj = 0.0;
-  fZproj = 0.0;
-  fSlope = 0.0;
-  fDetector = -1;
+  fYproj      = 0.0;
+  fZproj      = 0.0;
+  fSlope      = 0.0;
+  fDetector   = -1;
   fNtracklets = 0;
-  fNplanes = 0;
-  fNclusters = 0;
-  fPt = 0.0;
-  fPhi = 0.0;
-  fEta = 0.0;
-  fLabel = -1;
-  fPID = 0.0;
+  fNplanes    = 0;
+  fNclusters  = 0;
+  fPt         = 0.0;
+  fPhi        = 0.0;
+  fEta        = 0.0;
+  fLabel      = -1;
+  fPID        = 0.0;
   fIsElectron = kFALSE;
   
 }
@@ -200,28 +219,34 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
   // Calculate the kinematics of the found track
   //
 
+  Int_t  i    = 0;
+  Int_t  iDet = 0;
+  Int_t  nDet = 0;
+  Bool_t newDetector;
+
   AliTRDltuTracklet *trk;
   Int_t nTracklets = GetNtracklets();
   Float_t fC[kNmaxTrk][3];            // X, Y, Z  coordinates of segments
 
-  fYproj = 0.0;
-  fZproj = 0.0;
-  fSlope = 0.0;
-  fNclusters = 0;
-  fNplanes = 0;
+  fYproj      = 0.0;
+  fZproj      = 0.0;
+  fSlope      = 0.0;
+  fNclusters  = 0;
+  fNplanes    = 0;
   fNtracklets = GetNtracklets();
-  Int_t inDetector[kNplan];
-  for (Int_t i = 0; i < kNplan; i++) inDetector[i] = -1;
-  Int_t iDet, nDet = 0;
-  Bool_t newDetector;
-  for (Int_t i = 0; i < nTracklets; i++) {
+  Int_t inDetector[kNlayer];
+  for (i = 0; i < kNlayer; i++) {
+    inDetector[i] = -1;
+  }
 
-    trk = GetTracklet(i);
-    fYproj += trk->GetYproj(xpl);
-    fZproj += trk->GetZproj(xpl);
-    fSlope += trk->GetSlope();
+  for (i = 0; i < nTracklets; i++) {
+
+    trk         = GetTracklet(i);
+    fYproj     += trk->GetYproj(xpl);
+    fZproj     += trk->GetZproj(xpl);
+    fSlope     += trk->GetSlope();
     fNclusters += trk->GetNclusters();
-    iDet = trk->GetDetector();
+    iDet        = trk->GetDetector();
 
     newDetector = kTRUE;
     for (Int_t id = 0; id < nDet; id++) {
@@ -240,26 +265,32 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
     fC[i][2] = trk->GetRowz();
 
   }
-  fYproj /= (Float_t)nTracklets;
-  fZproj /= (Float_t)nTracklets;
-  fSlope /= (Float_t)nTracklets;
-
-  Float_t x[kNmaxTrk+1], y[kNmaxTrk+1], z[kNmaxTrk+1];
-  Bool_t count[kNmaxTrk];
-  for (Int_t i = 0; i < kNmaxTrk; i++) count[i] = kFALSE;
+  fYproj /= (Float_t) nTracklets;
+  fZproj /= (Float_t) nTracklets;
+  fSlope /= (Float_t) nTracklets;
+
+  Float_t x[kNmaxTrk+1];
+  Float_t y[kNmaxTrk+1];
+  Float_t z[kNmaxTrk+1];
+  Bool_t  count[kNmaxTrk];
+  for (i = 0; i < kNmaxTrk; i++) {
+    count[i] = kFALSE;
+  }
 
   Int_t iXmin = -1;
-  Int_t j = 0;
+  Int_t j     =  0;
   x[0] = y[0] = z[0] = 0.0;
   while (j < nTracklets) {
     iXmin = -1;
-    for (Int_t i = 0; i < nTracklets; i++) {
+    for (i = 0; i < nTracklets; i++) {
       if (count[i]) continue;
       if (iXmin == -1) {
        iXmin = i;
        continue;
       }
-      if (fC[i][0] < fC[iXmin][0]) iXmin = i;
+      if (fC[i][0] < fC[iXmin][0]) {
+        iXmin = i;
+      }
     }
     x[j+1] = fC[iXmin][0];
     y[j+1] = fC[iXmin][1];
@@ -276,9 +307,9 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
 
   smatrix.Zero();
   sums.Zero();
-  for (Int_t i = 0; i < nTracklets; i++) {
-    xv = (Double_t)x[i+1];
-    yv = (Double_t)y[i+1];
+  for (i = 0; i < nTracklets; i++) {
+    xv = (Double_t) x[i+1];
+    yv = (Double_t) y[i+1];
     smatrix(0,0) += 1.0;
     smatrix(1,1) += xv*xv;
     smatrix(0,1) += xv;
@@ -287,31 +318,31 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
     sums(1,0)    += xv*yv;
   }
   res = smatrix.Invert() * sums;
-  a = res(0,0);
-  b = res(1,0);
-
-  Float_t dist = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
-
-  Float_t fx1 = x[1]          + dist * (Float_t)(nTracklets-1)/6.0;
-  Float_t fy1 = a + b * fx1;
-  Float_t fx2 = x[nTracklets] - dist * (Float_t)(nTracklets-1)/6.0;
-  Float_t fy2 = a + b * fx2;
-  Float_t d12 = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
+  a   = res(0,0);
+  b   = res(1,0);
+
+  Float_t dist  = AliTRDgeometry::GetTime0(1) - AliTRDgeometry::GetTime0(0);
+  Float_t fx1   = x[1]          + dist * (Float_t) (nTracklets-1) / 6.0;
+  Float_t fy1   = a + b * fx1;
+  Float_t fx2   = x[nTracklets] - dist * (Float_t) (nTracklets-1) / 6.0;
+  Float_t fy2   = a + b * fx2;
+  Float_t d12   = TMath::Sqrt((fx2-fx1)*(fx2-fx1)+(fy2-fy1)*(fy2-fy1));
   Float_t alpha = TMath::ATan(fy2/fx2) - TMath::ATan(fy1/fx1);
-  Float_t r = (d12/2.0)/TMath::Sin(alpha);
+  Float_t r     = (d12 / 2.0) / TMath::Sin(alpha);
 
   fPt = 0.3 * field * 0.01 * r;
 
-  Float_t d1 = fx1*fx1+fy1*fy1;
-  Float_t d2 = fx2*fx2+fy2*fy2;
-  Float_t d  = fx1*fy2-fx2*fy1;
+  Float_t d1 = fx1*fx1 + fy1*fy1;
+  Float_t d2 = fx2*fx2 + fy2*fy2;
+  Float_t d  = fx1*fy2 - fx2*fy1;
   
-  Float_t xc = (d1*fy2-d2*fy1)/(2*d);
-  Float_t yc = (d2*fx1-d1*fx2)/(2*d);
+  Float_t xc = (d1*fy2 - d2*fy1) / (2.0*d);
+  Float_t yc = (d2*fx1 - d1*fx2) / (2.0*d);
 
   if (yc != 0.0) {
     fPhi = TMath::ATan(xc/yc);
-  } else {
+  } 
+  else {
     fPhi = TMath::PiOver2();
   }
 
@@ -319,9 +350,9 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
 
   smatrix.Zero();
   sums.Zero();
-  for (Int_t i = 0; i < nTracklets+1; i++) {
-    xv = (Double_t)z[i];
-    yv = (Double_t)x[i];
+  for (i = 0; i < nTracklets+1; i++) {
+    xv = (Double_t) z[i];
+    yv = (Double_t) x[i];
     smatrix(0,0) += 1.0;
     smatrix(1,1) += xv*xv;
     smatrix(0,1) += xv;
@@ -330,15 +361,17 @@ void AliTRDgtuTrack::Track(Float_t xpl, Float_t field)
     sums(1,0)    += xv*yv;
   }
   res = smatrix.Invert() * sums;
-  a = res(0,0);
-  b = res(1,0);
+  a   = res(0,0);
+  b   = res(1,0);
   Float_t theta = TMath::ATan(b);
   
-  if (theta < 0.0) theta = TMath::Pi() + theta;
-  
+  if (theta < 0.0) {
+    theta = TMath::Pi() + theta;
+  }
   if (theta == 0.0) {
     fEta = 0.0;
-  } else {
+  } 
+  else {
     fEta = -TMath::Log(TMath::Tan(theta/2.0));
   }
   
@@ -351,36 +384,49 @@ void AliTRDgtuTrack::MakePID()
   // Electron likelihood signal
   //
 
-  AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
-  if (!calibration)
-  {
-    Error("MakePID","No instance of AliTRDcalibDB.");
+  Int_t i = 0;
+
+  AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+  if (!calibration) {
+    AliError("No instance of AliTRDcalibDB.");
     return;  
   }
-  const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
+
+  AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
+  if (!rec) {
+    AliError("No TRD reco param.");
+    return;
+  }
+
+  const AliTRDCalPID *pd = calibration->GetPIDObject(rec->GetPIDMethod());
   
   AliTRDltuTracklet *trk;
-  Int_t nTracklets = GetNtracklets();
-  Int_t det, pla;
-  Float_t sl, th, q, probPio = 1.0, probEle = 1.0;
-  for (Int_t i = 0; i < nTracklets; i++) {
+  Int_t   nTracklets = GetNtracklets();
+  Int_t   det;
+  Int_t   pla;
+  Float_t sl;
+  Float_t th;
+  Float_t q;
+  Float_t probPio = 1.0;
+  Float_t probEle = 1.0;
+
+  for (i = 0; i < nTracklets; i++) {
 
     trk = GetTracklet(i);
 
-    sl = TMath::Abs(trk->GetSlope());     // tracklet inclination in X-y plane
-    th = trk->GetRowz()/trk->GetTime0();  // tracklet inclination in X-z plane
-    th = TMath::ATan(TMath::Abs(th));
+    sl  = TMath::Abs(trk->GetSlope());     // Tracklet inclination in X-y plane
+    th  = trk->GetRowz()/trk->GetTime0();  // Tracklet inclination in X-z plane
+    th  = TMath::ATan(TMath::Abs(th));
 
-    q = trk->GetQ() 
-      * TMath::Cos(sl/180.0*TMath::Pi()) 
-      * TMath::Cos(th/180.0*TMath::Pi());
+    q   = trk->GetQ() 
+        * TMath::Cos(sl/180.0*TMath::Pi()) 
+        * TMath::Cos(th/180.0*TMath::Pi());
 
     det = trk->GetDetector();
     pla = trk->GetPlane(det);
 
-    // unclear yet factor to match the PS distributions = 5.8
+    // Unclear yet factor to match the PS distributions = 5.8
     // not explained only by the tail filter ...
-
     // AliRoot v4-03-07 , v4-03-Release
     //q = q * 5.8;
 
@@ -410,14 +456,22 @@ void AliTRDgtuTrack::MakePID()
 
     // HEAD28Mar06 new distributions: factor = 6.0
 
-    probEle *= pd->GetProbability(0,TMath::Abs(fPt),q*6.0);
-    probPio *= pd->GetProbability(2,TMath::Abs(fPt),q*6.0);
+               // A.Bercuci on 2nd May 2007
+               // Temporary modification to deal with more energy slices. The values
+               // attached to charge slices and track length are dummy
+                Float_t dedx[3];
+                dedx[0] = dedx[1] = q*3.; dedx[2] = 0.;
+                Float_t length = 3.7;
+
+    probEle *= pd->GetProbability(0, TMath::Abs(fPt), dedx, length, 0);
+    probPio *= pd->GetProbability(2, TMath::Abs(fPt), dedx, length, 0);
 
   }
 
-  if ((probEle+probPio) > 0.0) {
-    fPID = probEle/(probEle+probPio);
-  } else {
+  if ((probEle + probPio) > 0.0) {
+    fPID = probEle/(probEle + probPio);
+  } 
+  else {
     fPID = 0.0;
   }
 
@@ -438,7 +492,7 @@ void AliTRDgtuTrack::MakePID()
   //Float_t fPIDcut = 0.829 - 0.032 * TMath::Abs(fPt);
 
   // HEAD28Mar06 with new generated distributions (pol2 fit)
-  Float_t absPt = TMath::Abs(fPt);
+  Float_t absPt   = TMath::Abs(fPt);
   //Float_t fPIDcut = 0.9575 - 0.0832 * absPt + 0.0037 * absPt*absPt;  // 800 bins in dEdx
   Float_t fPIDcut = 0.9482 - 0.0795 * absPt + 0.0035 * absPt*absPt;    // 200 bins in dEdx
 
@@ -463,16 +517,19 @@ void AliTRDgtuTrack::CookLabel()
   Int_t trackCount[kMaxTracks];
   for (Int_t it = 0; it < kMaxTracks; it++) {
     trackLabel[it] = -1;
-    trackCount[it] = 0;
+    trackCount[it] =  0;
   }
 
   Bool_t counted;
-  Int_t label, nTracks = 0;
+  Int_t  label;
+  Int_t  nTracks = 0;
   for (Int_t itrk = 0; itrk < fNtracklets; itrk++) {
 
     trk = GetTracklet(itrk);
 
-    if (trk->GetLabel() == -1) continue;
+    if (trk->GetLabel() == -1) {
+      continue;
+    }
 
     label = trk->GetLabel();
 
@@ -497,9 +554,9 @@ void AliTRDgtuTrack::CookLabel()
     
   }
 
-  Float_t frac = 4.0/5.0;
+  Float_t frac = 4.0 / 5.0;
   for (Int_t it = 0; it < kMaxTracks; it++) {
-    if (trackCount[it] >= (Int_t)(frac*fNtracklets)) {
+    if (trackCount[it] >= (Int_t) (frac*fNtracklets)) {
       fLabel = trackLabel[it];
       break;
     }
@@ -507,3 +564,46 @@ void AliTRDgtuTrack::CookLabel()
 
 }
 
+//_____________________________________________________________________________
+void AliTRDgtuTrack::ResetTracklets() 
+{ 
+  //
+  // Resets the list of tracklets
+  //
+
+  if (fTracklets) {
+    //fTracklets->Delete(); 
+    fTracklets->Clear();
+  }
+
+}
+
+//_____________________________________________________________________________
+TObjArray *AliTRDgtuTrack::Tracklets() 
+{ 
+  //
+  // Returns the list of tracklets
+  //
+
+  if (!fTracklets) {
+    fTracklets = new TObjArray(400); 
+  }
+
+  return fTracklets; 
+
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDgtuTrack::GetNtracklets() const 
+{
+  //
+  // Returns the number of tracklets
+  //
+
+  if (fTracklets) {
+    return fTracklets->GetEntriesFast();
+  }
+
+  return 0;
+
+}