]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDpidESD.cxx
Make use of new method AliRawReader::GetNumberOfEvents() - goinf to the last event...
[u/mrichter/AliRoot.git] / TRD / AliTRDpidESD.cxx
index aea62d597d2aff268f13df761f1d9a0817b6e9f4..5f0739763aee3953fe304f8a21dac06b344c1059 100644 (file)
 // Implementation of the TRD PID class                                    //
 //                                                                        //
 // Assigns the electron and pion likelihoods to each ESD track.           //
-// The function MakePID(AliESD *event) calculates the probability         //
+// The function MakePID(AliESDEvent *event) calculates the probability    //
 // of having dedx and a maximum timbin at a given                         //
 // momentum (mom) and particle type k                                     //
 // from the precalculated distributions.                                  //
 //                                                                        //
 // Authors :                                                              //
-// Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (Original version)//
-// Alex Bercuci (a.bercuci@gsi.de)                                        //
+//   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (orig. version) //
+//   Alex Bercuci (a.bercuci@gsi.de)                                      //
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
 #include "AliLog.h"
-#include "AliESD.h"
+#include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliTracker.h"
+#include "AliRun.h"
 
+#include "AliTRDReconstructor.h"
 #include "AliTRDpidESD.h"
 #include "AliTRDgeometry.h"
 #include "AliTRDcalibDB.h"
-#include "AliRun.h"
 #include "AliTRDtrack.h"
-#include "Cal/AliTRDCalPIDLQ.h"
-
+#include "Cal/AliTRDCalPID.h"
 
 ClassImp(AliTRDpidESD)
 
-  Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
-  Bool_t AliTRDpidESD::fCheckKinkStatus  = kFALSE;
-  Int_t AliTRDpidESD::fMinPlane          = 0;
+Bool_t  AliTRDpidESD::fgCheckTrackStatus = kTRUE;
+Bool_t  AliTRDpidESD::fgCheckKinkStatus  = kFALSE;
+Int_t   AliTRDpidESD::fgMinPlane         = 0;
 
 //_____________________________________________________________________________
 AliTRDpidESD::AliTRDpidESD()
-  :TObject(), fTrack(0x0)
+  :TObject()
+  ,fTrack(0x0)
 {
   //
   // Default constructor
@@ -62,7 +63,8 @@ AliTRDpidESD::AliTRDpidESD()
 
 //_____________________________________________________________________________
 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
-  :TObject(p), fTrack(0x0)
+  :TObject(p)
+  ,fTrack(0x0)
 {
   //
   // AliTRDpidESD copy constructor
@@ -78,7 +80,9 @@ AliTRDpidESD::~AliTRDpidESD()
   //
   // Destructor
   //
-       if(fTrack) delete fTrack;
+
+  if(fTrack) delete fTrack;
+
 }
 
 //_____________________________________________________________________________
@@ -100,15 +104,15 @@ void AliTRDpidESD::Copy(TObject &p) const
   // Copy function
   //
 
-  ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
-  ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
-  ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
+  ((AliTRDpidESD &) p).fgCheckTrackStatus         = fgCheckTrackStatus;
+  ((AliTRDpidESD &) p).fgCheckKinkStatus          = fgCheckKinkStatus;
+  ((AliTRDpidESD &) p).fgMinPlane                 = fgMinPlane;
   ((AliTRDpidESD &) p).fTrack                     = 0x0;
        
 }
 
 //_____________________________________________________________________________
-Int_t AliTRDpidESD::MakePID(AliESD *event)
+Int_t AliTRDpidESD::MakePID(AliESDEvent *event)
 {
   //
   // This function calculates the PID probabilities based on TRD signals
@@ -116,11 +120,12 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
   // 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
+  // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
   // implementation.
   //
   // Author
   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
+  //
 
        AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
        if (!calibration) {
@@ -129,11 +134,17 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
                return -1;
        }
        
+//   AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
+//   if (!rec) {
+//     AliErrorGeneral("AliTRDpidESD::MakePID()", "No TRD reco param.");
+//     return 0x0;
+//   }
+
        // Retrieve the CDB container class with the probability distributions
-       const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
+       const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDReconstructor::kLQPID/*rec->GetPIDMethod()*/);
        if (!pd) {
                AliErrorGeneral("AliTRDpidESD::MakePID()"
-                       ,"No access to AliTRDCalPIDLQ");
+                       ,"No access to AliTRDCalPID");
                return -1;
        }
 
@@ -158,19 +169,21 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
                mom          = 0.;
                length       = 0.;
                nPlanePID    = 0;
-               for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1./AliPID::kSPECIES;
-               for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
+               for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) 
+                  p[iSpecies] = 1./AliPID::kSPECIES;
+
+               for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) {
                        // read data for track segment
                        for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
-                               dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice);
-                       dEdx    = t->GetTRDsignals(iPlan, -1);
-                       timebin = t->GetTRDTimBin(iPlan);
+                               dedx[iSlice] = t->GetTRDslice(iLayer, iSlice);
+                       dEdx    = t->GetTRDslice(iLayer, -1);
+                       timebin = t->GetTRDTimBin(iLayer);
 
                        // check data
                        if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 
                        // retrive kinematic info for this track segment
-                       if(!RecalculateTrackSegmentKine(t, iPlan, mom, length)){
+                       if(!RecalculateTrackSegmentKine(t, iLayer, mom, length)){
                                // information is not fully reliable especialy for length
                                // estimation. To be used in the future. 
                        }
@@ -180,7 +193,7 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
 
                        // Get the probabilities for the different particle species
                        for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
-                               p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
+                               p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer);
                                //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
                        }
                }
@@ -190,7 +203,10 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
                Double_t probTotal = 0.;
                for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal   += p[iSpecies];
                if(probTotal <= 0.){
-                       AliWarningGeneral("AliTRDpidESD::MakePID()", Form("The total probability (%e) over all species <= 0 in ESD track %d. This may be caused by some error in reference data. Calculation continue but results might be corrupted.", probTotal, i));
+                       AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d."
+                                       , probTotal, i));
+                       AliWarning("This may be caused by some error in reference data.");
+                       AliWarning("Calculation continues but results might be corrupted.");
                        continue;
                }
                for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
@@ -201,6 +217,7 @@ Int_t AliTRDpidESD::MakePID(AliESD *event)
        }
        
        return 0;
+
 }
 
 //_____________________________________________________________________________
@@ -211,19 +228,23 @@ Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t)
   //
        
        // Check the ESD track status
-       if (fCheckTrackStatus) {
+       if (fgCheckTrackStatus) {
                if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
                        ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
        }
 
        // Check for ESD kink tracks
-       if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
+       if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
 
        return kTRUE;
+
 }
 
 //_____________________________________________________________________________
-Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd, Int_t plan, Float_t &mom, Float_t &length)
+Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd
+                                               , Int_t layer
+                                               , Float_t &mom
+                                               , Float_t &length)
 {
   //
   // Retrive momentum "mom" and track "length" in TRD chamber from plane
@@ -231,10 +252,11 @@ Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd, Int_t plan, F
   //
   // Origin
   // Alex Bercuci (A.Bercuci@gsi.de)   
+  //
 
        const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
         const Float_t kDrWidth     = AliTRDgeometry::DrThick();
-       const Float_t kTime0       = AliTRDgeometry::GetTime0(plan);
+       const Float_t kTime0       = AliTRDgeometry::GetTime0(layer);
 
        // set initial length value to chamber height 
        length = 2 * kAmHalfWidth + kDrWidth;
@@ -263,13 +285,13 @@ Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd, Int_t plan, F
                mom    = op->GetP();
                s      = op->GetSnp();
                t      = op->GetTgl();
-         length /= TMath::Sqrt((1. - s*s) / (1. - t*t));
+               if (s < 1.) length /= TMath::Sqrt((1. - s*s) / (1. + t*t));
                return kFALSE;
        }
        mom        = param->GetP();
        s = param->GetSnp();
        t = param->GetTgl();
-       length    /= TMath::Sqrt((1. - s*s) / (1. - t*t));
+       if (s < 1.) length    /= TMath::Sqrt((1. - s*s) / (1. + t*t));
 
        // check if track is crossing tracking sector by propagating to chamber exit- maybe is too much :)
        Double_t alpha = param->GetAlpha();
@@ -279,5 +301,5 @@ Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack *esd, Int_t plan, F
        if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE;
        
        return kTRUE;
-}
 
+}