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