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 ////////////////////////////////////////////////////////////////////////////
  17 //                                                                        //
  18 // Implementation of the TRD PID class                                    //
  19 //                                                                        //
  20 // Assigns the electron and pion likelihoods to each ESD track.           //
  21 // The function MakePID(AliESD *event) calculates the probability         //
  22 // of having dedx and a maximum timbin at a given                         //
  23 // momentum (mom) and particle type k                                     //
  24 // from the precalculated distributions.                                  //
  25 //                                                                        //
  26 // Original version:                                                      //
  27 // Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>                   //
  28 //                                                                        //
  29 ////////////////////////////////////////////////////////////////////////////
  30
  31 #include "AliLog.h"
  32 #include "AliESD.h"
  33 #include "AliESDtrack.h"
  34
  35 #include "AliTRDpidESD.h"
  36 #include "AliTRDgeometry.h"
  37 #include "AliTRDcalibDB.h"
  38 #include "Cal/AliTRDCalPIDLQ.h"
  39
  40 ClassImp(AliTRDpidESD)
  41
  42   Bool_t AliTRDpidESD::fCheckTrackStatus = kTRUE;
  43   Bool_t AliTRDpidESD::fCheckKinkStatus  = kFALSE;
  44   Int_t AliTRDpidESD::fMinPlane         = 0;
  45
  46 //_____________________________________________________________________________
  47 AliTRDpidESD::AliTRDpidESD():TObject()
  48 {
  49   //
  50   // Default constructor
  51   //
  52
  53
  54 }
  55
  56 //_____________________________________________________________________________
  57 AliTRDpidESD::AliTRDpidESD(const AliTRDpidESD &p):TObject(p)
  58 {
  59   //
  60   // AliTRDpidESD copy constructor
  61   //
  62
  63   ((AliTRDpidESD &) p).Copy(*this);
  64
  65 }
  66
  67 //_____________________________________________________________________________
  68 AliTRDpidESD &AliTRDpidESD::operator=(const AliTRDpidESD &p)
  69 {
  70   //
  71   // Assignment operator
  72   //
  73
  74   if (this != &p) ((AliTRDpidESD &) p).Copy(*this);
  75   return *this;
  76
  77 }
  78
  79 //_____________________________________________________________________________
  80 void AliTRDpidESD::Copy(TObject &p) const
  81 {
  82   //
  83   // Copy function
  84   //
  85
  86  ((AliTRDpidESD &) p).fCheckTrackStatus          = fCheckTrackStatus;
  87  ((AliTRDpidESD &) p).fCheckKinkStatus           = fCheckKinkStatus;
  88  ((AliTRDpidESD &) p).fMinPlane                  = fMinPlane;
  89
  90 }
  91
  92 //_____________________________________________________________________________
  93 Int_t AliTRDpidESD::MakePID(AliESD *event)
  94 {
  95   //
  96   // This function calculates the PID probabilities based on TRD signals
  97   //
  98   // So far this method produces probabilities based on the total charge
  99   // in each layer and the position of the maximum time bin in each layer.
 100   // In a final version this should also exploit the charge measurement in
 101   // the different slices of a given layer.
 102   //
 103
 104   Double_t p[10];
 105   Int_t    nSpecies  = AliPID::kSPECIES;
 106   Int_t    nPlanePID = 0;
 107   Double_t mom       = 0.0;
 108   Double_t probTotal = 0.0;
 109
 110   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 111   if (!calibration) {
 112     //AliError("No access to calibration data\n");
 113     return -1;
 114   }
 115
 116   // Retrieve the CDB container class with the probability distributions
 117   const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
 118   if (!pd) {
 119     //AliError("No access to AliTRDCalPIDLQ\n");
 120     return -1;
 121   }
 122
 123   // Loop through all ESD tracks
 124   Int_t ntrk = event->GetNumberOfTracks();
 125   for (Int_t i = 0; i < ntrk; i++) {
 126
 127     AliESDtrack *t = event->GetTrack(i);
 128
 129     // Check the ESD track status
 130     if (fCheckTrackStatus) {
 131       if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
 132           ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) {
 133         continue;
 134       }
 135     }
 136
 137     // Check for ESD kink tracks
 138     if (fCheckKinkStatus) {
 139       if (t->GetKinkIndex(0) != 0) {
 140         continue;
 141       }
 142     }
 143
 144     // Skip tracks that have no TRD signal at all
 145     if (t->GetTRDsignal() == 0) {
 146       continue;
 147     }
 148
 149     mom       = t->GetP();
 150     probTotal = 0.0;
 151     nPlanePID = 0;
 152     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 153       p[iSpecies] = 1.0;
 154     }
 155
 156     // Check the different detector layers
 157     for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
 158
 159       // Use the total charge in a given plane
 160       Double_t dedx    = t->GetTRDsignals(iPlan,-1);
 161       Int_t    timebin = t->GetTRDTimBin(iPlan);
 162       if ((dedx    >  0.0) &&
 163           (timebin > -1.0)) {
 164
 165         nPlanePID++;
 166
 167         // Get the probabilities for the different particle species
 168         for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 169
 170           p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx);
 171           p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin);
 172           p[iSpecies] *= 100.0; // ??????????????
 173
 174           probTotal   += p[iSpecies];
 175
 176         }
 177
 178       }
 179
 180     }
 181
 182     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 183       if ((probTotal >       0.0) &&
 184           (nPlanePID > fMinPlane)) {
 185         p[iSpecies] /= probTotal;
 186       }
 187       else {
 188         p[iSpecies] = 1.0;
 189       }
 190     }
 191
 192     t->SetTRDpid(p);
 193
 194   }
 195
 196   return 0;
 197
 198 }