Change of PID calculation by Alexandru
authorcblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Jul 2007 11:34:08 +0000 (11:34 +0000)
committercblume <cblume@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 26 Jul 2007 11:34:08 +0000 (11:34 +0000)
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/AliTRDtracker.h

index c29a59c..d53b358 100644 (file)
@@ -4,7 +4,7 @@
  * Author: The ALICE Off-line Project.                                    *
  * Contributors are mentioned in the code where appropriate.              *
  *                                                                        *
- * Permission to use, copy, modify and distribute this software and its   *
+ * Permission to use   , copy, modify and distribute this software and its   *
  * documentation strictly for non-commercial purposes is hereby granted   *
  * without fee, provided that the above copyright notice appears in all   *
  * copies and that both the copyright notice and this permission notice   *
 #include "AliTRDtrack.h"
 #include "AliTRDtracklet.h"
 
+// A. Bercuci - used for PID calculations
+#include "AliTRDcalibDB.h"
+#include "Cal/AliTRDCalPIDLQ.h"
+
 ClassImp(AliTRDtrack)
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -43,6 +47,7 @@ AliTRDtrack::AliTRDtrack()
   ,fSeedLab(-1)
   ,fdEdx(0)
   ,fDE(0)
+       ,fClusterOwner(kFALSE) // A.Bercuci
   ,fStopped(kFALSE)
   ,fLhElectron(0)
   ,fNWrong(0)
@@ -64,13 +69,19 @@ AliTRDtrack::AliTRDtrack()
       fdEdxPlane[i][j] = 0.0;
     }
     fTimBinPlane[i] = -1;
-  }
+               // A.Bercuci additions
+               fMom[i] = -1.;
+               fSnp[i] = 0.;
+               fTgl[i] = 0.;
+       }
 
   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
     fIndex[i]       = 0;
     fIndexBackup[i] = 0;
     fdQdl[i]        = 0;
-  }
+               //A.Bercuci additions
+               fClusters[i]    = 0x0;
+       }
 
   for (Int_t i = 0; i < 3; i++) {
     fBudget[i] = 0;
@@ -79,13 +90,14 @@ AliTRDtrack::AliTRDtrack()
 }
 
 //_____________________________________________________________________________
-AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index 
+AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
                        , const Double_t p[5], const Double_t cov[15] 
                        , Double_t x, Double_t alpha)
   :AliKalmanTrack()
   ,fSeedLab(-1)
   ,fdEdx(0)
   ,fDE(0)
+       ,fClusterOwner(kFALSE) // A.Bercuci
   ,fStopped(kFALSE)
   ,fLhElectron(0)
   ,fNWrong(0)
@@ -96,7 +108,7 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
   ,fNExpectedLast(0)
   ,fNdedx(0)
   ,fChi2Last(1e10)
-  ,fBackupTrack(0x0) 
+  ,fBackupTrack(0x0)
 {
   //
   // The main AliTRDtrack constructor.
@@ -125,12 +137,17 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
   Set(x,alpha,pp,cc);
   SetNumberOfClusters(1);
   fIndex[0] = index;
+       fClusters[0] = c;       // A.Bercuci additions
 
   for (Int_t i = 0; i < kNplane; i++) {
     for (Int_t j = 0; j < kNslice; j++) {
       fdEdxPlane[i][j] = 0.0;
     }
     fTimBinPlane[i] = -1;
+               // A.Bercuci additions
+               fMom[i] = -1.;
+               fSnp[i] = 0.;
+               fTgl[i] = 0.;
   }
 
   Double_t q = TMath::Abs(c->GetQ());
@@ -140,12 +157,13 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
     q *= TMath::Sqrt((1-s*s)/(1+t*t));
   }
 
-  fdQdl[0] = q;  
-
+  fdQdl[0]     = q;    
   for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
     fdQdl[i]        = 0;
     fIndex[i]       = 0;
     fIndexBackup[i] = 0;
+               // A.Bercuci additions
+               fClusters[i]    = 0x0;
   }
 
   for (Int_t i = 0; i < 3;i++) {
@@ -155,11 +173,12 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
 }                              
            
 //_____________________________________________________________________________
-AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
+AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
   :AliKalmanTrack(t) 
   ,fSeedLab(t.GetSeedLabel())
   ,fdEdx(t.fdEdx)
   ,fDE(t.fDE)
+       ,fClusterOwner(kTRUE) // A.Bercuci
   ,fStopped(t.fStopped)
   ,fLhElectron(0)
   ,fNWrong(t.fNWrong)
@@ -170,7 +189,7 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
   ,fNExpectedLast(t.fNExpectedLast)
   ,fNdedx(t.fNdedx)
   ,fChi2Last(t.fChi2Last)
-  ,fBackupTrack(0x0) 
+  ,fBackupTrack(0x0)
 {
   //
   // Copy constructor.
@@ -182,6 +201,10 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
     }
     fTimBinPlane[i] = t.fTimBinPlane[i];
     fTracklets[i]   = t.fTracklets[i];
+               // A.Bercuci additions
+               fMom[i] = t.fMom[i];
+               fSnp[i] = t.fSnp[i];
+               fTgl[i] = t.fTgl[i];
   }
 
   Int_t n = t.GetNumberOfClusters(); 
@@ -191,12 +214,17 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
     fIndex[i]       = t.fIndex[i];
     fIndexBackup[i] = t.fIndex[i];
     fdQdl[i]        = t.fdQdl[i];
-  }
+               // A.Bercuci additions
+               if(fClusterOwner && t.fClusters[i]) fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
+               else fClusters[i] = t.fClusters[i];
+       }
 
   for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
     fdQdl[i]        = 0;
     fIndex[i]       = 0;
     fIndexBackup[i] = 0;
+               // A.Bercuci additions
+               fClusters[i]    = 0x0;
   }
 
   for (Int_t i = 0; i < 3;i++) {
@@ -211,6 +239,7 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
   ,fSeedLab(-1)
   ,fdEdx(t.GetPIDsignal())
   ,fDE(0)
+       ,fClusterOwner(kFALSE) // A.Bercuci
   ,fStopped(kFALSE)
   ,fLhElectron(0.0)
   ,fNWrong(0)
@@ -237,12 +266,18 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
       fdEdxPlane[i][j] = 0.0;
     }
     fTimBinPlane[i] = -1;
+               // A.Bercuci additions
+               fMom[i] = -1.;
+               fSnp[i] = 0.;
+               fTgl[i] = 0.;
   }
 
   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
     fdQdl[i]        = 0;
     fIndex[i]       = 0;
     fIndexBackup[i] = 0;
+               // A.Bercuci additions
+               fClusters[i]    = 0x0;
   }
   
   for (Int_t i = 0; i < 3; i++) { 
@@ -257,6 +292,7 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
   ,fSeedLab(-1)
   ,fdEdx(t.GetTRDsignal())
   ,fDE(0)
+       ,fClusterOwner(kFALSE) // A.Bercuci
   ,fStopped(kFALSE)
   ,fLhElectron(0)
   ,fNWrong(0)
@@ -289,6 +325,10 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
       fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
     }
     fTimBinPlane[i] = t.GetTRDTimBin(i);
+               // A.Bercuci additions
+               fMom[i] = -1.;
+               fSnp[i] = 0.;
+               fTgl[i] = 0.;
   }
 
   const AliExternalTrackParam *par = &t;
@@ -303,8 +343,10 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
 
   
   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
-    fdQdl[i]   = 0;
-  }
+    fdQdl[i]     = 0;
+       // A.Bercuci additions
+               fClusters[i] = 0x0;
+       }
 
   for (Int_t i = 0; i < 3; i++) {
     fBudget[i] = 0;
@@ -332,7 +374,15 @@ AliTRDtrack::~AliTRDtrack()
   if (fBackupTrack) {
     delete fBackupTrack;
   }
-  fBackupTrack = 0;
+  fBackupTrack = 0x0;
+  if (fClusterOwner){
+    UInt_t icluster=0;
+               while(icluster<kMAXCLUSTERSPERTRACK && fClusters[icluster]){
+                       delete fClusters[icluster];
+                       fClusters[icluster] = 0x0;
+                       icluster++;
+               }
+  }
 
 }
 
@@ -393,38 +443,213 @@ Int_t AliTRDtrack::Compare(const TObject *o) const
 }                
 
 //_____________________________________________________________________________
-void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
+void AliTRDtrack::CookdEdx(Double_t low, Double_t up) 
 {
   //
   // Calculates the truncated dE/dx within the "low" and "up" cuts.
   //
 
   // Array to sort the dEdx values according to amplitude
-  Float_t sorted[kMAXCLUSTERSPERTRACK];
-  fdEdx = 0.0;
-   
+
+       Float_t sorted[kMAXCLUSTERSPERTRACK];
+       fdEdx = 0.;
+       
   // Require at least 10 clusters for a dedx measurement
   if (fNdedx < 10) return;
 
+
   // Can fdQdl be negative ????
-  for (Int_t i = 0; i < fNdedx; i++) {
-    sorted[i] = TMath::Abs(fdQdl[i]);
-  }
+  for (Int_t i = 0; i < fNdedx; i++) sorted[i] = TMath::Abs(fdQdl[i]);
   // Sort the dedx values by amplitude
   Int_t *index = new Int_t[fNdedx];
   TMath::Sort(fNdedx, sorted, index, kFALSE);
 
-  // Sum up the truncated charge between lower and upper bounds
+  // Sum up the truncated charge between lower and upper bounds 
   Int_t nl = Int_t(low * fNdedx);
   Int_t nu = Int_t( up * fNdedx);
-  for (Int_t i = nl; i <= nu; i++) {
-    fdEdx += sorted[index[i]];
+  for (Int_t i = nl; i <= nu; i++) fdEdx += sorted[index[i]];
+       fdEdx /= (nu - nl + 1.0);
+
+       delete[] index;
+}                     
+
+//_____________________________________________________________________________
+void AliTRDtrack::CookdEdxTimBin()
+{
+  //
+  // Set fdEdxPlane and fTimBinPlane and also get the 
+  // Time bin for Max. Cluster
+       //
+       // Authors:
+  // Prashant Shukla (shukla@physi.uni-heidelberg.de)
+  // Alexandru Bercuci (A.Bercuci@gsi.de)
+
+
+       Double_t  maxcharge[AliESDtrack::kNPlane]; // max charge in chamber
+       // number of clusters attached to track per chamber and slice
+       Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
+       //number of time bins in chamber
+       Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+       Int_t plane;   // plane of current cluster
+       Int_t tb;      // time bin of current cluster
+       Int_t slice;   // curent slice
+       AliTRDcluster *cluster = 0x0; // pointer to current cluster
+
+       // Reset class and local contors/variables
+  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
+               fTimBinPlane[iPlane] = -1;
+    maxcharge[iPlane] = 0.;
+               for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+      fdEdxPlane[iPlane][iSlice] = 0.;
+      nCluster[iPlane][iSlice]   = 0;
+    }
   }
-  fdEdx /= (nu - nl + 1.0);
 
-  delete[] index;
+       // start looping over clusters attached to this track
+       for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
+    cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
+    if(!cluster) continue;
+
+               // Read info from current cluster
+               plane  = AliTRDgeometry::GetPlane(cluster->GetDetector());
+    if (plane < 0 || plane >= AliESDtrack::kNPlane) {
+      AliError(Form("Wrong plane %d", plane));
+      continue;
+    }
+
+               tb     = cluster->GetLocalTimeBin();
+    if(tb == 0 || tb >= ntb){
+                       AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
+                       AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
+                       continue;
+               }
+       
+    slice   = tb * AliESDtrack::kNSlice / ntb;
+    
+               fdEdxPlane[plane][slice] += fdQdl[iClus];
+    if(fdQdl[iClus] > maxcharge[plane]) {
+      maxcharge[plane] = fdQdl[iClus];
+      fTimBinPlane[plane] = tb;
+    }
+    nCluster[plane][slice]++;
+  } // End of loop over cluster
+
+
+       
+  // Normalize fdEdxPlane to number of clusters and set track segments
+  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
+    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+      if (nCluster[iPlane][iSlice]) fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+    }
+       }
+}
+
+
+//_____________________________________________________________________________
+void   AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
+{
+       if(plane<0 || plane>= kNplane){
+               AliError(Form("Trying to access out of range plane (%d)", plane));
+               return;
+       }
+       
+       fSnp[plane] = GetSnp();
+       fTgl[plane] = GetTgl();
+       Double_t p[3]; GetPxPyPz(p);
+       fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
+}
+
+//_____________________________________________________________________________
+Float_t        AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
+{
+       if(plane < 0 || plane >= kNplane) return 0.;
+  return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1.
+- fSnp[plane]*fSnp[plane]) / (1. + fTgl[plane]*fTgl[plane]));
+
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
+{
+  //
+  // This function calculates the PID probabilities based on TRD signals
+  //
+  // The method produces probabilities based on the charge
+  // and the position of the maximum time bin in each layer.
+  // The dE/dx information can be used as global charge or 2 to 3
+  // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
+  // implementation.
+  //
+  // Author
+  // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
+
+
+       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+       if (!calibration) {
+               AliError("No access to calibration data");
+               return -1;
+       }
+       
+       // Retrieve the CDB container class with the probability distributions
+       const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
+       if (!pd) {
+               AliError("No access to AliTRDCalPIDLQ");
+               return -1;
+       }
+
+
+       Double_t p0 = 1./AliPID::kSPECIES;
+       if(AliPID::kSPECIES != 5){
+               AliError("Probabilities array defined only for 5 values. Please modify !!");
+               return -1;
+       }
+       Double_t p[] = {p0, p0, p0, p0, p0};
+       Float_t length;  // track segment length in chamber
+       Int_t nPlanePID = 0;
+               
+       // Skip tracks which have no TRD signal at all
+       if (fdEdx == 0.) return -1;
+       
+       for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
+
+               length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
+
+               // check data
+               if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <=  0.
+               || fTimBinPlane[iPlane] == -1.) continue;
+
+               
+               // this track segment has fulfilled all requierments
+               nPlanePID++;
+
+               // Get the probabilities for the different particle species
+               for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
+                       p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
+                       //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
+               }
+       }
+       if(nPlanePID == 0) return 0;
+
+       // normalize probabilities
+       Double_t probTotal = 0.;
+       for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
+
+       if(probTotal <= 0.) {
+               AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
+               return -1;
+       }
+       for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
+
+       
+       // book PID to the ESD track
+       esd->SetTRDpid(p);
+       esd->SetTRDpidQuality(nPlanePID);
+       
+       return 0;
+}
+
 
-}                    
 
 //_____________________________________________________________________________
 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
@@ -544,7 +769,8 @@ Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
   Int_t n = GetNumberOfClusters();
   fIndex[n] = index;
   SetNumberOfClusters(n+1);
-
+       
+       
   SetChi2(GetChi2()+chisq);
 
   return kTRUE;     
@@ -552,8 +778,7 @@ Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 }        
 
 //_____________________________________________________________________________
-Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index
-                          , Double_t h01, Int_t /*plane*/) 
+Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/)
 {
   //
   // Assignes found cluster to the track and updates track information
@@ -621,8 +846,10 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index
     return kFALSE;
   }
 
+       // Register cluster to track
   Int_t n = GetNumberOfClusters();
   fIndex[n] = index;
+       fClusters[n] = c; // A.Bercuci 25.07.07
   SetNumberOfClusters(n+1);
   SetChi2(GetChi2() + chisq);
 
index 058b96d..236fea3 100644 (file)
@@ -20,24 +20,22 @@ class AliESDtrack;
 class AliTrackReference;
 class AliTRDcluster;
 
-const unsigned kMAXCLUSTERSPERTRACK = 210; 
+const unsigned kMAXCLUSTERSPERTRACK = 210;
 
 class AliTRDtrack : public AliKalmanTrack {
 
  public:
    
-  enum { kNdet      = 540
+  enum { kNdet     = 540
        , kNstacks   =  90
        , kNplane    =   6
        , kNcham     =   5
        , kNsect     =  18
-       , kNslice    =   3
-       , kMidTimeBin=  14};
+       , kNslice    =   3};
 
    AliTRDtrack();
-   AliTRDtrack(const AliTRDcluster *c, Int_t index, const Double_t xx[5]
-              ,const Double_t cc[15], Double_t xr, Double_t alpha);  
-   AliTRDtrack(const AliTRDtrack &t);    
+   AliTRDtrack(/*const */AliTRDcluster *c, Int_t index, const Double_t xx[5], const Double_t cc[15], Double_t xr, Double_t alpha);
+   AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner = kTRUE*/);    
    AliTRDtrack(const AliESDtrack &t);
    ~AliTRDtrack();
    AliTRDtrack(const AliKalmanTrack &t, Double_t alpha); 
@@ -91,17 +89,29 @@ class AliTRDtrack : public AliKalmanTrack {
            void     SetNExpectedLast(Int_t nexp)                            { fNExpectedLast             = nexp;   }
            void     SetChi2Last(Float_t chi2)                               { fChi2Last                  = chi2;   }
            void     SetTracklets(Int_t i, AliTRDtracklet t)                 { fTracklets[i]              = t;      }
-           void     SetBudget(Int_t i, Float_t budget)                      { fBudget[i]                 = budget; }
-
-           Int_t    PropagateToX(Double_t xr, Double_t step);
-           Int_t    PropagateToR(Double_t xr, Double_t step);
-           Int_t    UpdateMI(const AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane);
-           void     MakeBackupTrack();
+                                        void     SetBudget(Int_t i, Float_t budget)                      { fBudget[i]                 = budget; }
+
+// A. Bercuci new functions
+                                       void            SetTrackSegmentDirMom(const Int_t plane);
+                                       void            CookdEdx(Double_t low = 0.05, Double_t up = 0.7);
+                                       void            CookdEdxTimBin();
+                                       Int_t           CookPID(AliESDtrack *p);
+                                       void            SetCluster(AliTRDcluster* cl, Int_t index = -1) {fClusters[(index == -1) ? GetNumberOfClusters()-1 : index] = cl;}
+       inline  AliTRDcluster* GetCluster(Int_t layer){ return (layer >= 0 && layer < GetNumberOfClusters()) ?  fClusters[layer] : 0x0;}
+       inline  Float_t GetMomentumPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fMom[plane] : 0.;}
+       inline  Float_t GetSnpPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fSnp[plane] : 0.;}
+       inline  Float_t GetTglPlane(Int_t plane) const {return (plane >= 0 && plane < kNplane) ? fTgl[plane] : 0.;}
+       inline  Float_t GetTrackLengthPlane(Int_t plane) const;
+//end A.Bercuci
+
+                                        void     MakeBackupTrack();
            Bool_t   PropagateTo(Double_t xr, Double_t x0 = 8.72, Double_t rho = 5.86e-3);
-           Bool_t   Update(const AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01);
+           Int_t    PropagateToR(Double_t xr, Double_t step);
+           Int_t    PropagateToX(Double_t xr, Double_t step);
            Bool_t   Rotate(Double_t angle, Bool_t absolute = kFALSE);
-           void     CookdEdx(Double_t low = 0.05, Double_t up = 0.7);   
            Float_t  StatusForTOF();
+           Bool_t   Update(const AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01);
+           Int_t    UpdateMI(/*const */AliTRDcluster *c, Double_t chi2, Int_t i, Double_t h01, Int_t plane);
  
  protected:
 
@@ -117,12 +127,20 @@ class AliTRDtrack : public AliKalmanTrack {
            Float_t  fdEdxPlane[kNplane][kNslice];       //  dE/dx from all 6 planes in 3 slices each
            Int_t    fTimBinPlane[kNplane];              //  Time bin of Max cluster from all 6 planes
 
-           Bool_t   fStopped;                           //  Track stop indication
+                                       // A.Bercuci
+                                        Float_t  fMom[kNplane];                      //  Track momentum at chamber entrance
+                                        Float_t  fSnp[kNplane];                   //  track direction
+                                        Float_t  fTgl[kNplane];                   //  track direction
+           AliTRDcluster*   fClusters[kMAXCLUSTERSPERTRACK];
+           Bool_t   fClusterOwner;         // indicates the track is owner of cluster
+           // end A.Bercuci
+                                        
+                                        Bool_t   fStopped;                           //  Track stop indication
            Int_t    fIndex[kMAXCLUSTERSPERTRACK];       //  Global indexes of clusters  
            Int_t    fIndexBackup[kMAXCLUSTERSPERTRACK]; //  Backup indexes of clusters - used in iterations
            Float_t  fdQdl[kMAXCLUSTERSPERTRACK];        //  Cluster amplitudes corrected for track angles    
-                           
-           Float_t  fLhElectron;                        //  Likelihood to be an electron    
+           
+                                        Float_t  fLhElectron;                        //  Likelihood to be an electron
            Int_t    fNWrong;                            //  Number of wrong clusters
            Int_t    fNRotate;                           //  Number of rotation
            Int_t    fNCross;                            //  Number of the cross materials
index 120b345..7585475 100644 (file)
@@ -716,7 +716,7 @@ Int_t AliTRDtracker::PropagateBack(AliESDEvent *event)
       Int_t foundClr = track->GetNumberOfClusters();
       if (foundClr >= foundMin) {
        track->CookdEdx(); 
-       CookdEdxTimBin(*track);
+       track->CookdEdxTimBin(); // A.Bercuci 25.07.07
        CookLabel(track,1 - fgkLabelFraction);
        if (track->GetBackupTrack()) {
           UseClusters(track->GetBackupTrack());
@@ -977,18 +977,21 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
     }  
         
     //AliTRDtrack *pt = seed2;
-    AliTRDtrack &t = *pt; 
-    FollowProlongation(t); 
-    if (t.GetNumberOfClusters() >= foundMin) {
+               // A.Bercuci 25.07.07
+               //AliTRDtrack &t = *pt; 
+    FollowProlongation(*pt); 
+    if (pt->GetNumberOfClusters() >= foundMin) {
       //UseClusters(&t);
       //CookLabel(pt, 1-fgkLabelFraction);
-      t.CookdEdx();
-      CookdEdxTimBin(t);
+      pt->CookdEdx();
+      pt->CookdEdxTimBin();
+                       pt->CookPID(seed);
+                       //pt->Calibrate(); // slot for calibration
     }
     found++;
 
     Double_t xTPC = 250.0;
-    if (PropagateToX(t,xTPC,fgkMaxStep)) {
+    if (PropagateToX(*pt,xTPC,fgkMaxStep)) {
 
       seed->UpdateTrackParams(pt,AliESDtrack::kTRDrefit);
       fHRefit->Fill(5);
@@ -1009,21 +1012,27 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
       delete seed2;
       if (PropagateToX(*pt2,xTPC,fgkMaxStep)) { 
         pt2->CookdEdx( ); 
-        CookdEdxTimBin(*pt2);
-       seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
-       fHRefit->Fill(6);
+        pt2->CookdEdxTimBin(); // A.Bercuci 25.07.07
+                               seed->UpdateTrackParams(pt2,AliESDtrack::kTRDrefit);
+                               fHRefit->Fill(6);
 
         for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) {
           for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) {
             seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
-         }
+                               }
           seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
         }
       }
+                       // A.Bercuci 25.07.07
+                       // Add TRD track to ESDfriendTrack - maybe this tracks are not useful for post-processing - TODO make decission
+                       if (AliTRDReconstructor::StreamLevel()>0)  seed->AddCalibObject(new AliTRDtrack(*pt2/*, kTRUE*/));
       delete pt2;
     }
 
     delete pt;
+               // A.Bercuci 25.07.07
+               // Add TRD track to ESDfriendTrack
+               if (AliTRDReconstructor::StreamLevel()>0)  seed->AddCalibObject(new AliTRDtrack(*pt/*, kTRUE*/));
 
   }   
 
@@ -1037,142 +1046,149 @@ Int_t AliTRDtracker::RefitInward(AliESDEvent *event)
 //_____________________________________________________________________________
 Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
 {
-  //
-  // Starting from current position on track=t this function tries
-  // to extrapolate the track up to timeBin=0 and to confirm prolongation
-  // if a close cluster is found. Returns the number of clusters
-  // expected to be found in sensitive layers
-  // GeoManager used to estimate mean density
-  //
-
-  Int_t    sector;
-  Int_t    lastplane = GetLastPlane(&t);
-  Double_t radLength = 0.0;
-  Double_t rho       = 0.0;
-  Int_t    expectedNumberOfClusters = 0;
-
-  for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
-
-    Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
-    Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
-
-    //
-    // Propagate track close to the plane if neccessary
-    //
-    Double_t currentx  = fTrSec[0]->GetLayer(rowlast)->GetX();
-    if (currentx < (-fgkMaxStep + t.GetX())) {
-      // Propagate closer to chamber - safety space fgkMaxStep      
-      if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
-        break;
-      }
-    }
-
-    if (!AdjustSector(&t)) {
-      break;
-    }
+       //
+       // Starting from current position on track=t this function tries
+       // to extrapolate the track up to timeBin=0 and to confirm prolongation
+       // if a close cluster is found. Returns the number of clusters
+       // expected to be found in sensitive layers
+       // GeoManager used to estimate mean density
+       //
 
-    //
-    // Get material budget
-    //
-    Double_t xyz0[3];
-    Double_t xyz1[3];
-    Double_t param[7];
-    Double_t x;
-    Double_t y;
-    Double_t z;
+       Int_t    sector;
+       Int_t    lastplane = GetLastPlane(&t);
+       Double_t radLength = 0.0;
+       Double_t rho       = 0.0;
+       Int_t    expectedNumberOfClusters = 0;
 
-    // Starting global position
-    t.GetXYZ(xyz0);   
-    // End global position
-    x = fTrSec[0]->GetLayer(row0)->GetX();
-    if (!t.GetProlongation(x,y,z)) {
-      break;
-    }
-    xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); 
-    xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
-    xyz1[2] =  z;
-    AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);       
-    rho       = param[0];
-    radLength = param[1]; // Get mean propagation parameters
+       for (Int_t iplane = lastplane; iplane >= 0; iplane--) {
 
-    //
-    // Propagate and update
-    //
-    sector = t.GetSector();
-    //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
-    for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
+               Int_t row0    = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1);
+               Int_t rowlast = GetGlobalTimeBin(0,iplane,0);
 
-      Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
-      expectedNumberOfClusters++;       
-      t.SetNExpected(t.GetNExpected() + 1);
-      if (t.GetX() > 345.0) {
-        t.SetNExpectedLast(t.GetNExpectedLast() + 1);
-      }
-      AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
-      AliTRDcluster *cl = 0;
-      UInt_t   index   = 0;
-      Double_t maxChi2 = fgkMaxChi2;
-      x = timeBin.GetX();
+               //
+               // Propagate track close to the plane if neccessary
+               //
+               Double_t currentx  = fTrSec[0]->GetLayer(rowlast)->GetX();
+               if (currentx < (-fgkMaxStep + t.GetX())) {
+                       // Propagate closer to chamber - safety space fgkMaxStep
+                       if (!PropagateToX(t,currentx+fgkMaxStep,fgkMaxStep)) {
+                               break;
+                       }
+               }
 
-      if (timeBin) {
+               if (!AdjustSector(&t)) {
+                       break;
+               }
 
-       AliTRDcluster *cl0 = timeBin[0];
-       if (!cl0) {
-          // No clusters in given time bin
-          continue;         
-       }
+               //
+               // Get material budget
+               //
+               Double_t xyz0[3];
+               Double_t xyz1[3];
+               Double_t param[7];
+               Double_t x;
+               Double_t y;
+               Double_t z;
+
+               // Starting global position
+               t.GetXYZ(xyz0);
+               // End global position
+               x = fTrSec[0]->GetLayer(row0)->GetX();
+               if (!t.GetProlongation(x,y,z)) {
+                       break;
+               }
+               xyz1[0] =  x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
+               xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
+               xyz1[2] =  z;
+               AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+               rho       = param[0];
+               radLength = param[1]; // Get mean propagation parameters
+               // A.Bercuci 25.07.07
+               // Flags for marking the track momentum and direction. The track is
+               // marked only if it has at least 1 cluster picked up in the current
+               // chamber.
+               Bool_t kUPDATED = kFALSE, kMARKED = kFALSE;
 
-       Int_t plane   = fGeom->GetPlane(cl0->GetDetector());
-       if (plane > lastplane) {
-          continue;
+               //
+               // Propagate and update
+               //
+               sector = t.GetSector();
+               //for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
+               for (Int_t itime = 0 ; itime < GetTimeBinsPerPlane(); itime++) {
+                       // A.Bercuci 25.07.07
+                       // Mark track kinematics
+                       if(itime > 10 && kUPDATED && !kMARKED){
+                               t.SetTrackSegmentDirMom(iplane);
+                               kMARKED = kTRUE;
+                       }
+
+                       Int_t ilayer = GetGlobalTimeBin(0,iplane,itime);
+                       expectedNumberOfClusters++;
+                       t.SetNExpected(t.GetNExpected() + 1);
+                       if (t.GetX() > 345.0) {
+                               t.SetNExpectedLast(t.GetNExpectedLast() + 1);
+                       }
+                       AliTRDpropagationLayer &timeBin = *(fTrSec[sector]->GetLayer(ilayer));
+                       AliTRDcluster *cl = 0;
+                       UInt_t   index   = 0;
+                       Double_t maxChi2 = fgkMaxChi2;
+                       x = timeBin.GetX();
+
+                       if (timeBin) {
+
+                               AliTRDcluster *cl0 = timeBin[0];
+                               if (!cl0) {
+                                       // No clusters in given time bin
+                                       continue;
+                               }
+
+                               Int_t plane   = fGeom->GetPlane(cl0->GetDetector());
+                               if (plane > lastplane) {
+                                                               continue;
+                               }
+
+                               Int_t timebin = cl0->GetLocalTimeBin();
+                               AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
+
+                               if (cl2) {
+                                       cl = cl2;
+                                       //Double_t h01 = GetTiltFactor(cl);    //I.B's fix
+                                       //maxChi2=t.GetPredictedChi2(cl,h01);
+                               }
+                               if (cl) {
+                                       //if (cl->GetNPads()<5)
+                                       Double_t dxsample = timeBin.GetdX();
+                                       t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample));
+                                                               Double_t h01      = GetTiltFactor(cl);
+                                       Int_t    det      = cl->GetDetector();
+                                       Int_t    plane    = fGeom->GetPlane(det);
+                                       if (t.GetX() > 345.0) {
+                                               t.SetNLast(t.GetNLast() + 1);
+                                               t.SetChi2Last(t.GetChi2Last() + maxChi2);
+                                       }
+
+                                       Double_t xcluster = cl->GetX();
+                                       t.PropagateTo(xcluster,radLength,rho);
+                                       
+                                       if (!AdjustSector(&t)) {
+                                               break;           //I.B's fix
+                                       }
+                                       maxChi2 = t.GetPredictedChi2(cl,h01);
+                                       
+                                       if (maxChi2<1e+10)
+                                               if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
+                                                       // ????
+                                               } else {
+                                                       // A.Bercuci 25.07.07
+                                                       //SetCluster(cl, GetNumberOfClusters()-1);
+                                                       kUPDATED = kTRUE;
+                                               }
+                               }
+                       }
+               }
        }
 
-       Int_t timebin = cl0->GetLocalTimeBin();
-       AliTRDcluster *cl2 = GetCluster(&t,plane,timebin,index);
-
-       if (cl2) {
-
-         cl = cl2;     
-         //Double_t h01 = GetTiltFactor(cl);    //I.B's fix
-         //maxChi2=t.GetPredictedChi2(cl,h01);
-
-       }       
-        if (cl) {
-
-         //if (cl->GetNPads()<5) 
-         Double_t dxsample = timeBin.GetdX();
-         t.SetSampledEdx(TMath::Abs(cl->GetQ()/dxsample)); 
-          Double_t h01      = GetTiltFactor(cl);
-         Int_t    det      = cl->GetDetector();    
-         Int_t    plane    = fGeom->GetPlane(det);
-         if (t.GetX() > 345.0) {
-           t.SetNLast(t.GetNLast() + 1);
-           t.SetChi2Last(t.GetChi2Last() + maxChi2);
-         }
-
-         Double_t xcluster = cl->GetX();
-         t.PropagateTo(xcluster,radLength,rho);
-
-          if (!AdjustSector(&t)) {
-            break;           //I.B's fix
-         }
-         maxChi2 = t.GetPredictedChi2(cl,h01);
-          
-         if (maxChi2<1e+10)
-           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
-             // ????
-           }
-
-       }                       
-
-      }
-
-    } 
-
-  }
-
-  return expectedNumberOfClusters;  
-
+       return expectedNumberOfClusters;
 }                
 
 //_____________________________________________________________________________
@@ -1324,7 +1340,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
              if (!t.Update(cl,maxChi2,index,h01)) {
                // ????
              }
-           }  
+           } // else SetCluster(cl, GetNumberOfClusters()-1); // A.Bercuci 25.07.07
+  
 
           if (calibra->GetMITracking()) {
             calibra->UpdateHistograms(cl,&t);
@@ -3437,94 +3454,6 @@ Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster *c)
 
 }
 
-//_____________________________________________________________________________
-void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
-{
-  //
-  // This is setting fdEdxPlane and fTimBinPlane
-  // Sums up the charge in each plane for track TRDtrack and also get the 
-  // Time bin for Max. Cluster
-  // Prashant Shukla (shukla@physi.uni-heidelberg.de)
-  //
-
-  Double_t  clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
-  Double_t  maxclscharge[AliESDtrack::kNPlane];
-  Int_t     nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
-  Int_t     timebin[AliESDtrack::kNPlane];
-
-  // Initialization of cluster charge per plane.  
-  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
-    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      clscharge[iPlane][iSlice] = 0.0;
-      nCluster[iPlane][iSlice]  = 0;
-    }
-  }
-
-  // Initialization of cluster charge per plane.  
-  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
-    timebin[iPlane]      =  -1;
-    maxclscharge[iPlane] = 0.0;
-  }
-
-  // Loop through all clusters associated to track TRDtrack
-  Int_t nClus = TRDtrack.GetNumberOfClusters();  // from Kalmantrack
-  for (Int_t iClus = 0; iClus < nClus; iClus++) {
-    Double_t charge = TRDtrack.GetClusterdQdl(iClus);
-    Int_t    index  = TRDtrack.GetClusterIndex(iClus);
-    AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index); 
-    if (!pTRDcluster) {
-      continue;
-    }
-    Int_t    tb     = pTRDcluster->GetLocalTimeBin();
-    if (!tb) {
-      continue;
-    }
-    Int_t detector  = pTRDcluster->GetDetector();
-    Int_t iPlane    = fGeom->GetPlane(detector);
-    if (iPlane >= AliESDtrack::kNPlane) {
-      AliError(Form("Wrong plane %d",iPlane));
-      continue;
-    }
-    Int_t iSlice    = tb * AliESDtrack::kNSlice 
-                         / AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
-    if (iSlice >= AliESDtrack::kNSlice) {
-      AliError(Form("Wrong slice %d",iSlice));
-      continue;
-    }
-    clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
-    if (charge > maxclscharge[iPlane]) {
-      maxclscharge[iPlane] = charge;
-      timebin[iPlane]      = tb;
-    }
-    nCluster[iPlane][iSlice]++;
-  } // End of loop over cluster
-
-  // Setting the fdEdxPlane and fTimBinPlane variabales 
-  Double_t totalCharge = 0.0;
-
-  for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
-    for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
-      if (nCluster[iPlane][iSlice]) {
-        clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
-      }
-      TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice],iPlane,iSlice);
-      totalCharge = totalCharge+clscharge[iPlane][iSlice];
-    }
-    TRDtrack.SetPIDTimBin(timebin[iPlane],iPlane);     
-  }
-
-  // Still needed ????
-  //  Int_t i;
-  //  Int_t nc=TRDtrack.GetNumberOfClusters(); 
-  //  Float_t dedx=0;
-  //  for (i=0; i<nc; i++) dedx += TRDtrack.GetClusterdQdl(i);
-  //  dedx /= nc;
-  //  for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
-  //    TRDtrack.SetPIDsignals(dedx, iPlane);
-  //    TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
-  //  }
-
-}
 
 //_____________________________________________________________________________
 Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
@@ -4291,7 +4220,8 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
                                       ,c
                                       ,params[0]
                                       ,params[6]*alpha+shift);
-  track->PropagateTo(params[0]-5.0);
+       // SetCluster(cl, 0); // A. Bercuci
+       track->PropagateTo(params[0]-5.0);
   track->ResetCovariance(1);
 
   Int_t rc = FollowBackProlongation(*track);
@@ -4301,7 +4231,7 @@ AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed *seeds, Double_t *params)
   }
   else {
     track->CookdEdx();
-    CookdEdxTimBin(*track);
+    track->CookdEdxTimBin();
     CookLabel(track,0.9);
   }
 
index 164e907..89a0fba 100644 (file)
@@ -265,7 +265,6 @@ class AliTRDtracker : public AliTracker {
   
   Int_t    FollowBackProlongation(AliTRDtrack &t);
   Int_t    FollowProlongation(AliTRDtrack &t);
-  void     CookdEdxTimBin(AliTRDtrack &t);  
   Int_t    PropagateToX(AliTRDtrack &t, Double_t xToGo, Double_t maxStep);
   Double_t ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const;
   Double_t ExpectedSigmaZ2(Double_t r, Double_t tgl) const;