TRD/AliTRDpidESD.cxx

   1 /**************************************************************************
   2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ////////////////////////////////////////////////////////////////////////////
  19 //                                                                        //
  20 // Implementation of the TRD PID class                                    //
  21 //                                                                        //
  22 // Assigns the electron and pion likelihoods to each ESD track.           //
  23 // The function MakePID(AliESDEvent *event) calculates the probability    //
  24 // of having dedx and a maximum timbin at a given                         //
  25 // momentum (mom) and particle type k                                     //
  26 // from the precalculated distributions.                                  //
  27 //                                                                        //
  28 // Authors :                                                              //
  29 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de> (orig. version) //
  30 //   Alex Bercuci (a.bercuci@gsi.de)                                      //
  31 //                                                                        //
  32 ////////////////////////////////////////////////////////////////////////////
  33
  34 #include "AliESDEvent.h"
  35 #include "AliESDtrack.h"
  36 #include "AliTracker.h"
  37
  38 #include "AliTRDpidESD.h"
  39 #include "AliTRDgeometry.h"
  40 #include "AliTRDcalibDB.h"
  41 #include "Cal/AliTRDCalPID.h"
  42
  43 ClassImp(AliTRDpidESD)
  44
  45 Bool_t  AliTRDpidESD::fgCheckTrackStatus = kTRUE;
  46 Bool_t  AliTRDpidESD::fgCheckKinkStatus  = kFALSE;
  47 Int_t   AliTRDpidESD::fgMinPlane         = 0;
  48
  49 //_____________________________________________________________________________
  50 AliTRDpidESD::AliTRDpidESD()
  51   :TObject()
  52   ,fTrack(NULL)
  53 {
  54   //
  55   // Default constructor
  56   //
  57
  58 }
  59
  60 //_____________________________________________________________________________
  61 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
  62   :TObject(p)
  63   ,fTrack(NULL)
  64 {
  65   //
  66   // AliTRDpidESD copy constructor
  67   //
  68
  69   ((AliTRDpidESD &) p).Copy(*this);
  70
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliTRDpidESD::~AliTRDpidESD()
  75 {
  76   //
  77   // Destructor
  78   //
  79
  80   if(fTrack) delete fTrack;
  81
  82 }
  83
  84 //_____________________________________________________________________________
  85 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
  86 {
  87   //
  88   // Assignment operator
  89   //
  90
  91   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
  92   return *this;
  93
  94 }
  95
  96 //_____________________________________________________________________________
  97 void AliTRDpidESD::Copy(TObject &p) const
  98 {
  99   //
 100   // Copy function
 101   //
 102
 103   ((AliTRDpidESD &) p).fgCheckTrackStatus         = fgCheckTrackStatus;
 104   ((AliTRDpidESD &) p).fgCheckKinkStatus          = fgCheckKinkStatus;
 105   ((AliTRDpidESD &) p).fgMinPlane                 = fgMinPlane;
 106   ((AliTRDpidESD &) p).fTrack                     = NULL;
 107
 108 }
 109
 110 //_____________________________________________________________________________
 111 Int_t AliTRDpidESD::MakePID(AliESDEvent * const event)
 112 {
 113   //
 114   // This function calculates the PID probabilities based on TRD signals
 115   //
 116   // The method produces probabilities based on the charge
 117   // and the position of the maximum time bin in each layer.
 118   // The dE/dx information can be used as global charge or 2 to 3
 119   // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
 120   // implementation.
 121   //
 122   // Author
 123   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
 124   //
 125
 126   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 127   if (!calibration) {
 128     AliErrorGeneral("AliTRDpidESD::MakePID()"
 129       ,"No access to calibration data");
 130     return -1;
 131   }
 132
 133 //   AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
 134 //   if (!rec) {
 135 //     AliErrorGeneral("AliTRDpidESD::MakePID()", "No TRD reco param.");
 136 //     return 0x0;
 137 //   }
 138
 139   // Retrieve the CDB container class with the probability distributions
 140   const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDpidUtil::kLQ);
 141   if (!pd) {
 142     AliErrorGeneral("AliTRDpidESD::MakePID()"
 143       ,"No access to AliTRDCalPID");
 144     return -1;
 145   }
 146
 147   // Loop through all ESD tracks
 148   Double_t p[10];
 149   AliESDtrack *t = NULL;
 150   Float_t dedx[AliTRDCalPID::kNSlicesLQ], dEdx;
 151   Int_t   timebin;
 152   Float_t mom, length;
 153   Int_t   nPlanePID;
 154   for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
 155     t = event->GetTrack(i);
 156
 157     // Check track
 158     if(!CheckTrack(t)) continue;
 159
 160
 161     // Skip tracks which have no TRD signal at all
 162     if (t->GetTRDsignal()<1.e-5) continue;
 163
 164     // Loop over detector layers
 165     mom          = 0.;
 166     length       = 0.;
 167     nPlanePID    = 0;
 168     for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
 169                   p[iSpecies] = 1./AliPID::kSPECIES;
 170
 171     for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) {
 172       // read data for track segment
 173       for(int iSlice=0; iSlice<AliTRDCalPID::kNSlicesLQ; iSlice++)
 174         dedx[iSlice] = t->GetTRDslice(iLayer, iSlice);
 175       dEdx    = t->GetTRDslice(iLayer, -1);
 176       timebin = t->GetTRDTimBin(iLayer);
 177
 178       // check data
 179       if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 180
 181       // retrive kinematic info for this track segment
 182       if(!RecalculateTrackSegmentKine(t, iLayer, mom, length)){
 183         // information is not fully reliable especialy for length
 184         // estimation. To be used in the future.
 185       }
 186
 187       // this track segment has fulfilled all requierments
 188       nPlanePID++;
 189
 190       // Get the probabilities for the different particle species
 191       for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
 192         p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length, iLayer);
 193         //p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
 194       }
 195     }
 196     if(nPlanePID == 0) continue;
 197
 198     // normalize probabilities
 199     Double_t probTotal = 0.;
 200     for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal   += p[iSpecies];
 201     if(probTotal <= 0.){
 202       AliWarning(Form("The total probability (%e) over all species <= 0 in ESD track %d."
 203                                       , probTotal, i));
 204       AliWarning("This may be caused by some error in reference data.");
 205       AliWarning("Calculation continues but results might be corrupted.");
 206       continue;
 207     }
 208     for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
 209
 210     // book PID to the track
 211     t->SetTRDpid(p);
 212     t->SetTRDntracklets(nPlanePID<<3);
 213   }
 214
 215   return 0;
 216
 217 }
 218
 219 //_____________________________________________________________________________
 220 Bool_t AliTRDpidESD::CheckTrack(AliESDtrack * const t)
 221 {
 222   //
 223   // Check if track is eligible for PID calculations
 224   //
 225
 226   // Check the ESD track status
 227   if (fgCheckTrackStatus) {
 228     if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
 229       ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
 230   }
 231
 232   // Check for ESD kink tracks
 233   if (fgCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
 234
 235   return kTRUE;
 236
 237 }
 238
 239 //_____________________________________________________________________________
 240 Bool_t AliTRDpidESD::RecalculateTrackSegmentKine(AliESDtrack * const esd
 241                                               , Int_t layer
 242                                               , Float_t &mom
 243                                               , Float_t &length)
 244 {
 245   //
 246   // Retrive momentum "mom" and track "length" in TRD chamber from plane
 247   // "plan" according to information stored in AliESDtrack "t".
 248   //
 249   // Origin
 250   // Alex Bercuci (A.Bercuci@gsi.de)
 251   //
 252
 253   const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
 254         const Float_t kDrWidth     = AliTRDgeometry::DrThick();
 255   const Float_t kTime0       = AliTRDgeometry::GetTime0(layer);
 256
 257   // set initial length value to chamber height
 258   length = 2 * kAmHalfWidth + kDrWidth;
 259
 260   // retrive track's outer param
 261   const AliExternalTrackParam *op = esd->GetOuterParam();
 262   if(!op){
 263     mom    = esd->GetP();
 264     return kFALSE;
 265   }
 266
 267   AliExternalTrackParam *param = NULL;
 268   if(!fTrack){
 269     fTrack = new AliExternalTrackParam(*op);
 270     param = fTrack;
 271   } else param = new(fTrack) AliExternalTrackParam(*op);
 272
 273   // retrive the magnetic field
 274   Double_t xyz0[3];
 275   op->GetXYZ(xyz0);
 276   Float_t field = AliTracker::GetBz(xyz0); // Bz in kG at point xyz0
 277   Double_t s, t;
 278
 279   // propagate to chamber entrance
 280   if(!param->PropagateTo(kTime0-kAmHalfWidth-kDrWidth, field)){
 281     mom    = op->GetP();
 282     s      = op->GetSnp();
 283     t      = op->GetTgl();
 284     if (s < 1.) length /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
 285     return kFALSE;
 286   }
 287   mom        = param->GetP();
 288   s = param->GetSnp();
 289   t = param->GetTgl();
 290   if (s < 1.) length    /= TMath::Sqrt((1.-s)*(1.+s) / (1. + t*t));
 291
 292   // check if track is crossing tracking sector by propagating to chamber exit- maybe is too much :)
 293   Double_t alpha = param->GetAlpha();
 294   if(!param->PropagateTo(kTime0+kAmHalfWidth, field)) return kFALSE;
 295
 296   // mark track segments which are crossing SM boundaries along chamber
 297   if(TMath::Abs(alpha-param->GetAlpha())>.01) return kFALSE;
 298
 299   return kTRUE;
 300
 301 }