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