STEER/AliESDpid.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 //           Implementation of the combined PID class
  20 //           For the Event Summary Data Class
  21 //           produced by the reconstruction process
  22 //           and containing information on the particle identification
  23 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  24 //-----------------------------------------------------------------
  25
  26 #include "AliLog.h"
  27 #include "AliPID.h"
  28 #include "AliESDpid.h"
  29 #include "AliESDEvent.h"
  30 #include "AliESDtrack.h"
  31
  32 ClassImp(AliESDpid)
  33
  34 Int_t AliESDpid::MakePID(AliESDEvent *event, Bool_t TPConly, Float_t TimeZeroTOF) const {
  35   //
  36   //  Calculate probabilities for all detectors, except if TPConly==kTRUE
  37   //  and combine PID
  38   //
  39   //   Option TPConly==kTRUE is used during reconstruction,
  40   //  because ITS tracking uses TPC pid
  41   //  HMPID and TRD pid are done in detector reconstructors
  42   //
  43
  44   /*
  45   Float_t TimeZeroTOF = 0;
  46   if (subtractT0)
  47     TimeZeroTOF = event->GetT0();
  48   */
  49   Int_t nTrk=event->GetNumberOfTracks();
  50   for (Int_t iTrk=0; iTrk<nTrk; iTrk++) {
  51     AliESDtrack *track=event->GetTrack(iTrk);
  52     MakeTPCPID(track);
  53     if (!TPConly) {
  54       MakeITSPID(track);
  55       MakeTOFPID(track, TimeZeroTOF);
  56       //MakeHMPIDPID(track);
  57       //MakeTRDPID(track);
  58     }
  59     CombinePID(track);
  60   }
  61   return 0;
  62 }
  63 //_________________________________________________________________________
  64 void AliESDpid::MakeTPCPID(AliESDtrack *track) const
  65 {
  66   //
  67   //  TPC pid using bethe-bloch and gaussian response
  68   //
  69   if ((track->GetStatus()&AliESDtrack::kTPCin )==0)
  70     if ((track->GetStatus()&AliESDtrack::kTPCout)==0) return;
  71
  72     Double_t mom = track->GetP();
  73     const AliExternalTrackParam *in=track->GetInnerParam();
  74     if (in) mom = in->GetP();
  75
  76     Double_t p[AliPID::kSPECIES];
  77     Double_t dedx=track->GetTPCsignal();
  78     Bool_t mismatch=kTRUE, heavy=kTRUE;
  79
  80     for (Int_t j=0; j<AliPID::kSPECIES; j++) {
  81       AliPID::EParticleType type=AliPID::EParticleType(j);
  82       Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
  83       Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
  84       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
  85         p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
  86       } else {
  87         p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
  88         mismatch=kFALSE;
  89       }
  90
  91       // Check for particles heavier than (AliPID::kSPECIES - 1)
  92       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
  93
  94     }
  95
  96     if (mismatch)
  97        for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
  98
  99     track->SetTPCpid(p);
 100
 101     if (heavy) track->ResetStatus(AliESDtrack::kTPCpid);
 102
 103 }
 104 //_________________________________________________________________________
 105 void AliESDpid::MakeITSPID(AliESDtrack *track) const
 106 {
 107   //
 108   // ITS PID
 109   // Two options, depending on fITSPIDmethod:
 110   //  1) Truncated mean method
 111   //  2) Likelihood, using charges measured in all 4 layers and
 112   //     Landau+gaus response functions
 113   //
 114
 115   if ((track->GetStatus()&AliESDtrack::kITSin)==0 &&
 116       (track->GetStatus()&AliESDtrack::kITSout)==0) return;
 117
 118   Double_t mom=track->GetP();
 119   if (fITSPIDmethod == kITSTruncMean) {
 120     Double_t dedx=track->GetITSsignal();
 121     Bool_t isSA=kTRUE;
 122     Double_t momITS=mom;
 123     ULong_t trStatus=track->GetStatus();
 124     if(trStatus&AliESDtrack::kTPCin){
 125       isSA=kFALSE;
 126       const AliExternalTrackParam *in = track->GetInnerParam();
 127       momITS=in->P();
 128     }
 129     UChar_t clumap=track->GetITSClusterMap();
 130     Int_t nPointsForPid=0;
 131     for(Int_t i=2; i<6; i++){
 132       if(clumap&(1<<i)) ++nPointsForPid;
 133     }
 134
 135     if(nPointsForPid<3) { // track not to be used for combined PID purposes
 136       track->ResetStatus(AliESDtrack::kITSpid);
 137       return;
 138     }
 139
 140     Double_t p[10];
 141
 142     Bool_t mismatch=kTRUE, heavy=kTRUE;
 143     for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 144       Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
 145       Double_t bethe=fITSResponse.Bethe(momITS,mass);
 146       Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
 147       if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 148         p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 149       } else {
 150         p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 151         mismatch=kFALSE;
 152       }
 153
 154       // Check for particles heavier than (AliPID::kSPECIES - 1)
 155       if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
 156
 157     }
 158
 159     if (mismatch)
 160        for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 161
 162     track->SetITSpid(p);
 163
 164     if (heavy) track->ResetStatus(AliESDtrack::kITSpid);
 165   }
 166   else {  // Likelihood method
 167     Double_t condprobfun[AliPID::kSPECIES];
 168     Double_t qclu[4];
 169     track->GetITSdEdxSamples(qclu);
 170     fITSResponse.GetITSProbabilities(mom,qclu,condprobfun);
 171     track->SetITSpid(condprobfun);
 172   }
 173
 174 }
 175 //_________________________________________________________________________
 176 void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t TimeZeroTOF) const
 177 {
 178   //
 179   //   TOF PID using gaussian response
 180   //
 181   if ((track->GetStatus()&AliESDtrack::kTOFout)==0) return;
 182   if ((track->GetStatus()&AliESDtrack::kTIME)==0) return;
 183
 184   Double_t time[AliPID::kSPECIESN];
 185   track->GetIntegratedTimes(time);
 186
 187   Double_t sigma[AliPID::kSPECIES];
 188   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
 189     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->GetP(),time[iPart],AliPID::ParticleMass(iPart));
 190   }
 191
 192   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 193            Form("Expected TOF signals [ps]: %f %f %f %f %f",
 194                   time[AliPID::kElectron],
 195                   time[AliPID::kMuon],
 196                   time[AliPID::kPion],
 197                   time[AliPID::kKaon],
 198                   time[AliPID::kProton]));
 199
 200   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 201            Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
 202                   sigma[AliPID::kElectron],
 203                   sigma[AliPID::kMuon],
 204                   sigma[AliPID::kPion],
 205                   sigma[AliPID::kKaon],
 206                   sigma[AliPID::kProton]
 207                   ));
 208
 209   Double_t tof = track->GetTOFsignal() - TimeZeroTOF;
 210
 211   Double_t p[AliPID::kSPECIES];
 212   Bool_t mismatch = kTRUE, heavy = kTRUE;
 213   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 214     Double_t sig = sigma[j];
 215     if (TMath::Abs(tof-time[j]) > fRange*sig) {
 216       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
 217     } else
 218       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
 219
 220     // Check the mismatching
 221     Double_t mass = AliPID::ParticleMass(j);
 222     Double_t pm = fTOFResponse.GetMismatchProbability(track->GetP(),mass);
 223     if (p[j]>pm) mismatch = kFALSE;
 224
 225     // Check for particles heavier than (AliPID::kSPECIES - 1)
 226     if (tof < (time[j] + fRange*sig)) heavy=kFALSE;
 227
 228   }
 229
 230   if (mismatch)
 231     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
 232
 233   track->SetTOFpid(p);
 234
 235   if (heavy) track->ResetStatus(AliESDtrack::kTOFpid);
 236 }
 237 //_________________________________________________________________________
 238 void AliESDpid::MakeTRDPID(AliESDtrack *track) const
 239 {
 240   //
 241   // Method to recalculate the TRD PID probabilities
 242   //
 243   if((track->GetStatus()&AliESDtrack::kTRDout)==0) return;
 244   Double_t prob[AliPID::kSPECIES]; Float_t mom[6];
 245   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
 246   for(Int_t ilayer = 0; ilayer < 6; ilayer++){
 247     mom[ilayer] = track->GetTRDmomentum(ilayer);
 248     for(Int_t islice = 0; islice < track->GetNumberOfTRDslices(); islice++){
 249       dedx[ilayer*track->GetNumberOfTRDslices()+islice] = track->GetTRDslice(ilayer, islice);
 250     }
 251   }
 252   fTRDResponse.GetResponse(track->GetNumberOfTRDslices(), dedx, mom, prob);
 253   track->SetTRDpid(prob);
 254 }
 255 //_________________________________________________________________________
 256 void AliESDpid::CombinePID(AliESDtrack *track) const
 257 {
 258   //
 259   // Combine the information of various detectors
 260   // to determine the Particle Identification
 261   //
 262   Int_t ns=AliPID::kSPECIES;
 263   Double_t p[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
 264
 265   if (track->IsOn(AliESDtrack::kITSpid)) {
 266     Double_t d[10];
 267     track->GetITSpid(d);
 268     for (Int_t j=0; j<ns; j++) p[j]*=d[j];
 269   }
 270
 271   if (track->IsOn(AliESDtrack::kTPCpid)) {
 272     Double_t d[10];
 273     track->GetTPCpid(d);
 274     for (Int_t j=0; j<ns; j++) p[j]*=d[j];
 275   }
 276
 277   if (track->IsOn(AliESDtrack::kTRDpid)) {
 278     Double_t d[10];
 279     track->GetTRDpid(d);
 280     for (Int_t j=0; j<ns; j++) p[j]*=d[j];
 281   }
 282
 283   if (track->IsOn(AliESDtrack::kTOFpid)) {
 284     Double_t d[10];
 285     track->GetTOFpid(d);
 286     for (Int_t j=0; j<ns; j++) p[j]*=d[j];
 287   }
 288
 289   if (track->IsOn(AliESDtrack::kHMPIDpid)) {
 290     Double_t d[10];
 291     track->GetHMPIDpid(d);
 292     for (Int_t j=0; j<ns; j++) p[j]*=d[j];
 293   }
 294
 295   track->SetESDpid(p);
 296 }