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 ClassImp(AliAODpidUtil)
  38
  39   Int_t AliAODpidUtil::MakePID(const AliAODTrack *track,Double_t *p) const {
  40   //
  41   //  Calculate probabilities for all detectors, except if TPConly==kTRUE
  42   //  and combine PID
  43   //
  44   //   Option TPConly==kTRUE is used during reconstruction,
  45   //  because ITS tracking uses TPC pid
  46   //  HMPID and TRD pid are done in detector reconstructors
  47   //
  48
  49   /*
  50     Float_t TimeZeroTOF = 0;
  51     if (subtractT0)
  52     TimeZeroTOF = event->GetT0();
  53   */
  54   Int_t ns=AliPID::kSPECIES;
  55   Double_t tpcPid[AliPID::kSPECIES];
  56   MakeTPCPID(track,tpcPid);
  57   Double_t itsPid[AliPID::kSPECIES];
  58   Double_t tofPid[AliPID::kSPECIES];
  59   Double_t trdPid[AliPID::kSPECIES];
  60   MakeITSPID(track,itsPid);
  61   MakeTOFPID(track,tofPid);
  62   //MakeHMPIDPID(track);
  63   MakeTRDPID(track,trdPid);
  64   for (Int_t j=0; j<ns; j++) {
  65     p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
  66   }
  67
  68   return 0;
  69 }
  70 //_________________________________________________________________________
  71 Float_t AliAODpidUtil::GetTPCsignalTunedOnData(const AliVTrack *t) const {
  72     AliAODTrack *track = (AliAODTrack *) t;
  73     Float_t dedx = track->GetTPCsignalTunedOnData();
  74     if(dedx > 0) return dedx;
  75
  76     Double_t mom = t->GetTPCmomentum();
  77
  78     dedx = t->GetTPCsignal();
  79     track->SetTPCsignalTunedOnData(dedx);
  80
  81     if(dedx < 20) return dedx;
  82
  83
  84     AliPID::EParticleType type = AliPID::kPion;
  85
  86     AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(track->GetAODEvent()->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
  87     if (mcHeader) {
  88         TClonesArray *mcArray = (TClonesArray*)track->GetAODEvent()->GetList()->FindObject(AliAODMCParticle::StdBranchName());
  89
  90         Bool_t kGood = kTRUE;
  91
  92         Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(t->GetLabel())))->GetPdgCode());
  93         if(iS==AliPID::ParticleCode(AliPID::kElectron)){
  94             type = AliPID::kElectron;
  95         }
  96         else if(iS==AliPID::ParticleCode(AliPID::kMuon)){
  97             type = AliPID::kMuon;
  98         }
  99         else if(iS==AliPID::ParticleCode(AliPID::kPion)){
 100             type = AliPID::kPion;
 101         }
 102         else if(iS==AliPID::ParticleCode(AliPID::kKaon)){
 103             type = AliPID::kKaon;
 104         }
 105         else if(iS==AliPID::ParticleCode(AliPID::kProton)){
 106             type = AliPID::kProton;
 107         }
 108         else if(iS==AliPID::ParticleCode(AliPID::kDeuteron)){ // d
 109             type = AliPID::kDeuteron;
 110         }
 111         else if(iS==AliPID::ParticleCode(AliPID::kTriton)){ // t
 112             type = AliPID::kTriton;
 113         }
 114         else if(iS==AliPID::ParticleCode(AliPID::kHe3)){ // 3He
 115             type = AliPID::kHe3;
 116         }
 117         else if(iS==AliPID::ParticleCode(AliPID::kAlpha)){ // 4He
 118             type = AliPID::kAlpha;
 119         }
 120         else
 121             kGood = kFALSE;
 122
 123         if(kGood){
 124             Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 125             Double_t sigma=fTPCResponse.GetExpectedSigma(mom,t->GetTPCsignalN(),type);
 126             dedx = gRandom->Gaus(bethe,sigma);
 127
 128             if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
 129         }
 130
 131     }
 132
 133     track->SetTPCsignalTunedOnData(dedx);
 134     return dedx;
 135 }
 136 //_________________________________________________________________________
 137 void AliAODpidUtil::MakeTPCPID(const AliAODTrack *track,Double_t *p) const
 138 {
 139   //
 140   //  TPC pid using bethe-bloch and gaussian response
 141   //
 142
 143   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return;
 144
 145   AliAODPid *pidObj = track->GetDetPid();
 146   Double_t mom      = track->P();
 147   Double_t dedx     = 0.;
 148   UShort_t nTPCClus = 0;
 149   if (pidObj) {
 150       nTPCClus = pidObj->GetTPCsignalN();
 151       dedx     = pidObj->GetTPCsignal();
 152       mom      = pidObj->GetTPCmomentum();
 153   }
 154
 155   Bool_t mismatch=kTRUE;
 156
 157   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 158     AliPID::EParticleType type=AliPID::EParticleType(j);
 159     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
 160     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
 161     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 162       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 163     } else {
 164       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 165       mismatch=kFALSE;
 166     }
 167
 168   }
 169
 170   if (mismatch)
 171     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 172
 173
 174   return;
 175 }
 176 //_________________________________________________________________________
 177 void AliAODpidUtil::MakeITSPID(const AliAODTrack *track,Double_t *p) const
 178 {
 179   //
 180   // ITS PID
 181   //  1) Truncated mean method
 182   //
 183
 184
 185   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
 186   UChar_t clumap=track->GetITSClusterMap();
 187   Int_t nPointsForPid=0;
 188   for(Int_t i=2; i<6; i++){
 189    if(clumap&(1<<i)) ++nPointsForPid;
 190   }
 191   if(nPointsForPid<3) { // track not to be used for combined PID purposes
 192     for (Int_t j=0; j<AliPID::kSPECIES; j++)
 193       p[j] = 1./AliPID::kSPECIES;
 194     return;
 195   }
 196   Double_t mom=track->P();
 197   AliAODPid *pidObj = track->GetDetPid();
 198
 199   Double_t dedx = 0.;
 200   if (pidObj) {
 201       dedx = pidObj->GetITSsignal();
 202   }
 203
 204   Bool_t mismatch = kTRUE;
 205   Bool_t isSA = kTRUE;
 206   if(track->GetStatus() & AliESDtrack::kTPCin){
 207     isSA = kFALSE;
 208     if (pidObj)
 209       mom = pidObj->GetTPCmomentum();
 210   }
 211   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 212     Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
 213     Double_t bethe = fITSResponse.Bethe(mom,mass);
 214     Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
 215     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 216       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 217     } else {
 218       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 219       mismatch=kFALSE;
 220     }
 221
 222     // Check for particles heavier than (AliPID::kSPECIES - 1)
 223
 224   }
 225
 226   if (mismatch)
 227     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 228
 229   return;
 230
 231 }
 232 //_________________________________________________________________________
 233 void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const
 234 {
 235   //
 236   //   TOF PID using gaussian response
 237   //
 238   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return;
 239   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
 240   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
 241
 242   Double_t time[AliPID::kSPECIESN];
 243   Double_t sigma[AliPID::kSPECIESN];
 244   AliAODPid *pidObj = track->GetDetPid();
 245   pidObj->GetIntegratedTimes(time);
 246
 247   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
 248     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
 249   }
 250
 251   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 252                   Form("Expected TOF signals [ps]: %f %f %f %f %f",
 253                        time[AliPID::kElectron],
 254                        time[AliPID::kMuon],
 255                        time[AliPID::kPion],
 256                        time[AliPID::kKaon],
 257                        time[AliPID::kProton]));
 258
 259   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 260                   Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
 261                        sigma[AliPID::kElectron],
 262                        sigma[AliPID::kMuon],
 263                        sigma[AliPID::kPion],
 264                        sigma[AliPID::kKaon],
 265                        sigma[AliPID::kProton]
 266                        ));
 267
 268   Double_t tof = pidObj->GetTOFsignal();
 269
 270   Bool_t mismatch = kTRUE;
 271   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 272     Double_t sig = sigma[j];
 273     if (TMath::Abs(tof-time[j]) > fRange*sig) {
 274       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
 275     } else
 276       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
 277
 278     // Check the mismatching
 279     Double_t mass = AliPID::ParticleMass(j);
 280     Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
 281     if (p[j]>pm) mismatch = kFALSE;
 282
 283     // Check for particles heavier than (AliPID::kSPECIES - 1)
 284
 285   }
 286
 287   if (mismatch)
 288     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 289
 290   return;
 291 }
 292 //_________________________________________________________________________
 293 void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const
 294 {
 295   ComputeTRDProbability(track, AliPID::kSPECIES, p);
 296   return;
 297 }
 298