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(AliESD *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> (Original version)//
  30 // Alex Bercuci (a.bercuci@gsi.de)                                        //
  31 //                                                                        //
  32 ////////////////////////////////////////////////////////////////////////////
  33
  34 #include "AliLog.h"
  35 #include "AliESD.h"
  36 #include "AliESDtrack.h"
  37
  38 #include "AliTRDpidESD.h"
  39 #include "AliTRDgeometry.h"
  40 #include "AliTRDcalibDB.h"
  41 #include "AliRun.h"
  42 #include "AliTRDtrack.h"
  43 #include "Cal/AliTRDCalPIDLQ.h"
  44
  45 ClassImp(AliTRDpidESD)
  46
  47   Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
  48   Bool_t AliTRDpidESD::fCheckKinkStatus  = kFALSE;
  49   Int_t AliTRDpidESD::fMinPlane          = 0;
  50
  51 //_____________________________________________________________________________
  52 AliTRDpidESD::AliTRDpidESD()
  53   :TObject()
  54 {
  55   //
  56   // Default constructor
  57   //
  58
  59 }
  60
  61 //_____________________________________________________________________________
  62 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p)
  63   :TObject(p)
  64 {
  65   //
  66   // AliTRDpidESD copy constructor
  67   //
  68
  69   ((AliTRDpidESD &) p).Copy(*this);
  70
  71 }
  72
  73 //_____________________________________________________________________________
  74 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
  75 {
  76   //
  77   // Assignment operator
  78   //
  79
  80   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
  81   return *this;
  82
  83 }
  84
  85 //_____________________________________________________________________________
  86 void AliTRDpidESD::Copy(TObject &p) const
  87 {
  88   //
  89   // Copy function
  90   //
  91
  92   ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
  93   ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
  94   ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
  95
  96 }
  97
  98 //_____________________________________________________________________________
  99 Int_t AliTRDpidESD::MakePID(AliESD *event)
 100 {
 101   //
 102   // This function calculates the PID probabilities based on TRD signals
 103   //
 104   // The method produces probabilities based on the charge
 105   // and the position of the maximum time bin in each layer.
 106   // The dE/dx information can be used as global charge or 2 to 3
 107   // slices. Check AliTRDCalPIDLQ and AliTRDCalPIDLQRef for the actual
 108   // implementation.
 109   //
 110   // Author
 111   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
 112
 113         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 114         if (!calibration) {
 115                 AliErrorGeneral("AliTRDpidESD::MakePID()"
 116                         ,"No access to calibration data");
 117                 return -1;
 118         }
 119
 120         // Retrieve the CDB container class with the probability distributions
 121         const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
 122         if (!pd) {
 123                 AliErrorGeneral("AliTRDpidESD::MakePID()"
 124                         ,"No access to AliTRDCalPIDLQ");
 125                 return -1;
 126         }
 127
 128
 129         // Loop through all ESD tracks
 130         Double_t p[10];
 131         AliESDtrack *t = 0x0;
 132         Double_t dedx[AliTRDtrack::kNslice], dEdx;
 133         Int_t    timebin;
 134         Float_t mom, length, probTotal;
 135         Int_t nPlanePID;
 136         for (Int_t i=0; i<event->GetNumberOfTracks(); i++) {
 137                 t = event->GetTrack(i);
 138
 139                 // Check track
 140                 if(!CheckTrack(t)) continue;
 141
 142                 // Skip tracks which have no TRD signal at all
 143                 if (t->GetTRDsignal() == 0.) continue;
 144
 145                 // Loop over detector layers
 146                 mom       = 0.; //t->GetP();
 147                 length    = 0.;
 148                 probTotal = 0.;
 149                 nPlanePID    = 0;
 150                 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] = 1.;
 151                 for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
 152                         // read data for track segment
 153                         for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
 154                                 dedx[iSlice] = t->GetTRDsignals(iPlan, iSlice);
 155                         dEdx    = t->GetTRDsignals(iPlan, -1);
 156                         timebin = t->GetTRDTimBin(iPlan);
 157
 158                         // check data
 159                         if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 160
 161                         // retrive kinematic info for this track segment
 162                         if(!GetTrackSegmentKine(t, iPlan, mom, length)) continue;
 163
 164                         // this track segment has fulfilled all requierments
 165                         nPlanePID++;
 166
 167                         // Get the probabilities for the different particle species
 168                         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
 169                                 p[iSpecies] *= pd->GetProbability(iSpecies, mom, dedx, length);
 170                                 p[iSpecies] *= pd->GetProbabilityT(iSpecies, mom, timebin);
 171                                 probTotal   += p[iSpecies];
 172                         }
 173                 }
 174
 175                 // normalize probabilities
 176                 if(probTotal > 0.)
 177                         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
 178                                 if(nPlanePID > fMinPlane) p[iSpecies] /= probTotal;
 179                                 else p[iSpecies] = 1.0;
 180
 181
 182                 // book PID to the track
 183                 t->SetTRDpid(p);
 184         }
 185
 186         return 0;
 187 }
 188
 189 //_____________________________________________________________________________
 190 Bool_t AliTRDpidESD::CheckTrack(AliESDtrack *t)
 191 {
 192   //
 193   // Check if track is eligible for PID calculations
 194   //
 195
 196         // Check the ESD track status
 197         if (fCheckTrackStatus) {
 198                 if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
 199                         ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) return kFALSE;
 200         }
 201
 202         // Check for ESD kink tracks
 203         if (fCheckKinkStatus && (t->GetKinkIndex(0) != 0)) return kFALSE;
 204
 205         return kTRUE;
 206 }
 207
 208 //_____________________________________________________________________________
 209 Bool_t AliTRDpidESD::GetTrackSegmentKine(AliESDtrack *t, Int_t plan, Float_t &mom, Float_t &length)
 210 {
 211   //
 212   // Retrive momentum "mom" and track "length" in TRD chamber from plane
 213   // "plan" according to information stored in AliESDtrack "t".
 214   //
 215
 216         if(!gAlice){
 217                 AliErrorGeneral("AliTRDpidESD::GetTrackSegmentKine()"
 218                 ,"No gAlice object to retrive TRDgeometry and Magnetic fied  - this has to be removed in the future.");
 219                 return kFALSE;
 220         }
 221
 222         // Retrieve TRD geometry -> Maybe there is a better way to do this
 223         AliTRDgeometry trdGeom;
 224         const Float_t kAmHalfWidth = AliTRDgeometry::AmThick() / 2.;
 225         const Float_t kDrWidth     = AliTRDgeometry::DrThick();
 226
 227
 228         // retrive the magnetic field
 229         Double_t xyz0[3] = { 0., 0., 0.}, xyz1[3];
 230         Double_t b[3], alpha;
 231         gAlice->Field(xyz0,b);      // b[] is in kilo Gauss
 232         Float_t field = b[2] * 0.1; // Tesla
 233
 234         // find momentum at chamber entrance and track length in chamber
 235         AliExternalTrackParam *param = (plan<3) ? new AliExternalTrackParam(*t->GetInnerParam()) : new AliExternalTrackParam(*t->GetOuterParam());
 236
 237         param->PropagateTo(trdGeom.GetTime0(plan)+kAmHalfWidth, field);
 238         param->GetXYZ(xyz0);
 239         alpha = param->GetAlpha();
 240         param->PropagateTo(trdGeom.GetTime0(plan)-kAmHalfWidth-kDrWidth, field);
 241         // eliminate track segments which are crossing SM boundaries along chamber
 242         if(TMath::Abs(alpha-param->GetAlpha())>.01){
 243                 delete param;
 244                 return kFALSE;
 245         }
 246         param->GetXYZ(xyz1);
 247         length = sqrt(
 248                 (xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0])+
 249                 (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1])+
 250                 (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2])
 251         );
 252         param->GetPxPyPz(xyz1);
 253         mom = sqrt(xyz1[0]*xyz1[0] + xyz1[1]*xyz1[1] + xyz1[2]*xyz1[2]);
 254         delete param;
 255
 256         return kTRUE;
 257 }
 258