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     AliErrorGeneral("AliTRDpidESD::MakePID"
 113                    ,"No access to calibration data\n");
 114     return -1;
 115   }
 116
 117   // Retrieve the CDB container class with the probability distributions
 118   const AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
 119   if (!pd) {
 120     AliErrorGeneral("AliTRDpidESD::MakePID"
 121                    ,"No access to AliTRDCalPIDLQ\n");
 122     return -1;
 123   }
 124
 125   // Loop through all ESD tracks
 126   Int_t ntrk = event->GetNumberOfTracks();
 127   for (Int_t i = 0; i < ntrk; i++) {
 128
 129     AliESDtrack *t = event->GetTrack(i);
 130
 131     // Check the ESD track status
 132     if (fCheckTrackStatus) {
 133       if (((t->GetStatus() & AliESDtrack::kTRDout  ) == 0) &&
 134           ((t->GetStatus() & AliESDtrack::kTRDrefit) == 0)) {
 135         continue;
 136       }
 137     }
 138
 139     // Check for ESD kink tracks
 140     if (fCheckKinkStatus) {
 141       if (t->GetKinkIndex(0) != 0) {
 142         continue;
 143       }
 144     }
 145
 146     // Skip tracks that have no TRD signal at all
 147     if (t->GetTRDsignal() == 0) {
 148       continue;
 149     }
 150
 151     mom       = t->GetP();
 152     probTotal = 0.0;
 153     nPlanePID = 0;
 154     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 155       p[iSpecies] = 1.0;
 156     }
 157
 158     // Check the different detector layers
 159     for (Int_t iPlan = 0; iPlan < AliTRDgeometry::kNplan; iPlan++) {
 160
 161       // Use the total charge in a given plane
 162       Double_t dedx    = t->GetTRDsignals(iPlan,-1);
 163       Int_t    timebin = t->GetTRDTimBin(iPlan);
 164       if ((dedx    >  0.0) &&
 165           (timebin > -1.0)) {
 166
 167         nPlanePID++;
 168
 169         // Get the probabilities for the different particle species
 170         for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 171
 172           p[iSpecies] *= pd->GetProbability(iSpecies,mom,dedx);
 173           p[iSpecies] *= pd->GetProbabilityT(iSpecies,mom,timebin);
 174           p[iSpecies] *= 100.0; // ??????????????
 175
 176           probTotal   += p[iSpecies];
 177
 178         }
 179
 180       }
 181
 182     }
 183
 184     for (Int_t iSpecies = 0; iSpecies < nSpecies; iSpecies++) {
 185       if ((probTotal >       0.0) &&
 186           (nPlanePID > fMinPlane)) {
 187         p[iSpecies] /= probTotal;
 188       }
 189       else {
 190         p[iSpecies] = 1.0;
 191       }
 192     }
 193
 194     t->SetTRDpid(p);
 195
 196   }
 197
 198   return 0;
 199
 200 }