STEER/AOD/AliAODpidUtil.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: AliAODpidUtil.cxx 38329 2010-01-17 19:17:24Z hristov $ */
  17
  18 //-----------------------------------------------------------------
  19 //           Implementation of the combined PID class
  20 //           For the AOD Class
  21 //           containing information on the particle identification
  22 //      Origin: Rosa Romita, GSI, r.romita@gsi.de
  23 //-----------------------------------------------------------------
  24
  25 #include "TRandom.h"
  26 #include "AliLog.h"
  27 #include "AliPID.h"
  28 #include "AliAODpidUtil.h"
  29 #include "AliAODEvent.h"
  30 #include "AliAODTrack.h"
  31 #include "AliAODPid.h"
  32 #include "AliTRDPIDResponse.h"
  33 #include "AliESDtrack.h"
  34 #include "AliAODMCHeader.h"
  35 #include "AliAODMCParticle.h"
  36
  37 #include <AliDetectorPID.h>
  38
  39 ClassImp(AliAODpidUtil)
  40
  41   Int_t AliAODpidUtil::MakePID(const AliAODTrack *track,Double_t *p) const {
  42   //
  43   //  Calculate probabilities for all detectors, except if TPConly==kTRUE
  44   //  and combine PID
  45   //
  46   //   Option TPConly==kTRUE is used during reconstruction,
  47   //  because ITS tracking uses TPC pid
  48   //  HMPID and TRD pid are done in detector reconstructors
  49   //
  50
  51   /*
  52     Float_t TimeZeroTOF = 0;
  53     if (subtractT0)
  54     TimeZeroTOF = event->GetT0();
  55   */
  56   Int_t ns=AliPID::kSPECIES;
  57   Double_t tpcPid[AliPID::kSPECIES];
  58   MakeTPCPID(track,tpcPid);
  59   Double_t itsPid[AliPID::kSPECIES];
  60   Double_t tofPid[AliPID::kSPECIES];
  61   Double_t trdPid[AliPID::kSPECIES];
  62   MakeITSPID(track,itsPid);
  63   MakeTOFPID(track,tofPid);
  64   //MakeHMPIDPID(track);
  65   MakeTRDPID(track,trdPid);
  66   for (Int_t j=0; j<ns; j++) {
  67     p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
  68   }
  69
  70   return 0;
  71 }
  72 //_________________________________________________________________________
  73 Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
  74     AliAODTrack *track = (AliAODTrack *) t;
  75     Float_t dedx = track->GetTPCsignalTunedOnData();
  76     if(dedx > 0) return dedx;
  77
  78     Double_t mom = t->GetTPCmomentum();
  79
  80     dedx = t->GetTPCsignal();
  81     track->SetTPCsignalTunedOnData(dedx);
  82
  83     if(dedx < 20) return dedx;
  84
  85
  86     AliPID::EParticleType type = AliPID::kPion;
  87
  88     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(track->GetAODEvent()->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
  89     if (mcHeader) {
  90         TClonesArray *mcArray = (TClonesArray*)track->GetAODEvent()->GetList()->FindObject(AliAODMCParticle::StdBranchName());
  91
  92         Bool_t kGood = kTRUE;
  93
  94         Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(t->GetLabel())))->GetPdgCode());
  95         if(iS==AliPID::ParticleCode(AliPID::kElectron)){
  96             type = AliPID::kElectron;
  97         }
  98         else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
  99             type = AliPID::kMuon;
 100         }
 101         else if(iS==AliPID::ParticleCode(AliPID::kPion)){
 102             type = AliPID::kPion;
 103         }
 104         else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
 105             type = AliPID::kKaon;
 106         }
 107         else if(iS==AliPID::ParticleCode(AliPID::kProton)){
 108             type = AliPID::kProton;
 109         }
 110         else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
 111             type = AliPID::kDeuteron;
 112         }
 113         else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
 114             type = AliPID::kTriton;
 115         }
 116         else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
 117             type = AliPID::kHe3;
 118         }
 119         else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
 120             type = AliPID::kAlpha;
 121         }
 122         else
 123             kGood = kFALSE;
 124
 125         if(kGood){
 126             Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 127             Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
 128             dedx = gRandom->Gaus(bethe,sigma);
 129
 130             if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
 131         }
 132
 133     }
 134
 135     track->SetTPCsignalTunedOnData(dedx);
 136     return dedx;
 137 }
 138 //_________________________________________________________________________
 139 void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const
 140 {
 141   //
 142   //  TPC pid using bethe-bloch and gaussian response
 143   //
 144
 145   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return;
 146
 147   AliAODPid *pidObj = track->GetDetPid();
 148   Double_t mom      = track->P();
 149   Double_t dedx     = 0.;
 150   UShort_t nTPCClus = 0;
 151   if (pidObj) {
 152       nTPCClus = pidObj->GetTPCsignalN();
 153       dedx     = pidObj->GetTPCsignal();
 154       mom      = pidObj->GetTPCmomentum();
 155   }
 156
 157   Bool_t mismatch=kTRUE;
 158
 159   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 160     AliPID::EParticleType type=AliPID::EParticleType(j);
 161     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 162     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
 163     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 164       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 165     } else {
 166       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 167       mismatch=kFALSE;
 168     }
 169
 170   }
 171
 172   if (mismatch)
 173     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 174
 175
 176   return;
 177 }
 178 //_________________________________________________________________________
 179 void AliAODpidUtil::MakeITSPID(const AliAODTrack *track,Double_t *p) const
 180 {
 181   //
 182   // ITS PID
 183   //  1) Truncated mean method
 184   //
 185
 186
 187   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
 188   UChar_t clumap=track->GetITSClusterMap();
 189   Int_t nPointsForPid=0;
 190   for(Int_t i=2; i<6; i++){
 191    if(clumap&(1<<i)) ++nPointsForPid;
 192   }
 193   if(nPointsForPid<3) { // track not to be used for combined PID purposes
 194     for (Int_t j=0; j<AliPID::kSPECIES; j++)
 195       p[j] = 1./AliPID::kSPECIES;
 196     return;
 197   }
 198   Double_t mom=track->P();
 199   AliAODPid *pidObj = track->GetDetPid();
 200
 201   Double_t dedx = 0.;
 202   if (pidObj) {
 203       dedx = pidObj->GetITSsignal();
 204   }
 205
 206   Bool_t mismatch = kTRUE;
 207   Bool_t isSA = kTRUE;
 208   if(track->GetStatus() & AliESDtrack::kTPCin){
 209     isSA = kFALSE;
 210     if (pidObj)
 211       mom = pidObj->GetTPCmomentum();
 212   }
 213   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 214     Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
 215     Double_t bethe = fITSResponse.Bethe(mom,mass);
 216     Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
 217     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 218       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 219     } else {
 220       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 221       mismatch=kFALSE;
 222     }
 223
 224     // Check for particles heavier than (AliPID::kSPECIES - 1)
 225
 226   }
 227
 228   if (mismatch)
 229     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 230
 231   return;
 232
 233 }
 234 //_________________________________________________________________________
 235 void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const
 236 {
 237   //
 238   //   TOF PID using gaussian response
 239   //
 240   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return;
 241   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
 242   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
 243
 244   Double_t time[AliPID::kSPECIES];
 245   Double_t sigma[AliPID::kSPECIES];
 246   AliAODPid *pidObj = track->GetDetPid();
 247   pidObj->GetIntegratedTimes(time);
 248
 249   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
 250     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
 251   }
 252
 253   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 254                   Form("Expected TOF signals [ps]: %f %f %f %f %f",
 255                        time[AliPID::kElectron],
 256                        time[AliPID::kMuon],
 257                        time[AliPID::kPion],
 258                        time[AliPID::kKaon],
 259                        time[AliPID::kProton]));
 260
 261   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 262                   Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
 263                        sigma[AliPID::kElectron],
 264                        sigma[AliPID::kMuon],
 265                        sigma[AliPID::kPion],
 266                        sigma[AliPID::kKaon],
 267                        sigma[AliPID::kProton]
 268                        ));
 269
 270   Double_t tof = pidObj->GetTOFsignal();
 271
 272   Bool_t mismatch = kTRUE;
 273   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 274     Double_t sig = sigma[j];
 275     if (TMath::Abs(tof-time[j]) > fRange*sig) {
 276       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
 277     } else
 278       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
 279
 280     // Check the mismatching
 281     Double_t mass = AliPID::ParticleMass(j);
 282     Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
 283     if (p[j]>pm) mismatch = kFALSE;
 284
 285     // Check for particles heavier than (AliPID::kSPECIES - 1)
 286
 287   }
 288
 289   if (mismatch)
 290     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 291
 292   return;
 293 }
 294 //_________________________________________________________________________
 295 void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const
 296 {
 297   ComputeTRDProbability(track, AliPID::kSPECIES, p);
 298   return;
 299 }
 300
 301 //_________________________________________________________________________
 302 Float_t AliAODpidUtil::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
 303 {
 304   //
 305   // Number of sigma implementation for the TOF
 306   //
 307
 308   AliAODTrack *track=(AliAODTrack*)vtrack;
 309
 310   // look for cached value first
 311   if (track->GetDetectorPID()){
 312     return track->GetDetectorPID()->GetNumberOfSigmas(kTOF, type);
 313   }
 314
 315   Bool_t oldAod=kTRUE;
 316   Double_t sigTOF;
 317   AliAODPid *pidObj = track->GetDetPid();
 318   if (!pidObj) return -999.;
 319   Double_t tofTime=pidObj->GetTOFsignal();
 320   Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type);
 321   Double_t sigmaTOFPid[AliPID::kSPECIES];
 322   pidObj->GetTOFpidResolution(sigmaTOFPid);
 323   AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
 324   if (event) {  // protection if the user didn't call GetTrack, which sets the internal pointer
 325     AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
 326     if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD
 327         sigTOF=fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type)); //fTOFResponse is set in InitialiseEvent
 328         tofTime -= fTOFResponse.GetStartTime(vtrack->P());
 329         oldAod=kFALSE;
 330     }
 331   } else {
 332     AliError("pointer to AliAODEvent not found, please call GetTrack to set it");
 333     return -996.;
 334   }
 335   if (oldAod) { // old AOD
 336     if (type <= AliPID::kProton) {
 337       sigTOF=sigmaTOFPid[type];
 338     } else return -998.;  // light nuclei cannot be supported on old AOD because we don't have timeZero resolution
 339   }
 340   if (sigTOF>0) return (tofTime - expTime)/sigTOF;
 341   else return -997.;
 342 }