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