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