STEER/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 "AliLog.h"
  26 #include "AliPID.h"
  27 #include "AliAODpidUtil.h"
  28 #include "AliAODEvent.h"
  29 #include "AliAODTrack.h"
  30 #include "AliAODPid.h"
  31 #include "AliTRDPIDResponse.h"
  32 #include "AliESDtrack.h"
  33
  34 ClassImp(AliAODpidUtil)
  35
  36   Int_t AliAODpidUtil::MakePID(AliAODTrack *track,Double_t *p) const {
  37   //
  38   //  Calculate probabilities for all detectors, except if TPConly==kTRUE
  39   //  and combine PID
  40   //
  41   //   Option TPConly==kTRUE is used during reconstruction,
  42   //  because ITS tracking uses TPC pid
  43   //  HMPID and TRD pid are done in detector reconstructors
  44   //
  45
  46   /*
  47     Float_t TimeZeroTOF = 0;
  48     if (subtractT0)
  49     TimeZeroTOF = event->GetT0();
  50   */
  51   Int_t ns=AliPID::kSPECIES;
  52   Double_t tpcPid[AliPID::kSPECIES];
  53   MakeTPCPID(track,tpcPid);
  54   Double_t itsPid[AliPID::kSPECIES];
  55   Double_t tofPid[AliPID::kSPECIES];
  56   Double_t trdPid[AliPID::kSPECIES];
  57   MakeITSPID(track,itsPid);
  58   MakeTOFPID(track,tofPid);
  59   //MakeHMPIDPID(track);
  60   MakeTRDPID(track,trdPid);
  61   for (Int_t j=0; j<ns; j++) {
  62     p[j]=tpcPid[j]*itsPid[j]*tofPid[j]*trdPid[j];
  63   }
  64
  65   return 0;
  66 }
  67 //_________________________________________________________________________
  68 void AliAODpidUtil::MakeTPCPID(AliAODTrack *track,Double_t *p) const
  69 {
  70   //
  71   //  TPC pid using bethe-bloch and gaussian response
  72   //
  73
  74   if ((track->GetStatus()&AliESDtrack::kTPCin )==0) return;
  75
  76   AliAODPid *pidObj = track->GetDetPid();
  77   Double_t mom      = track->P();
  78   Double_t dedx     = 0.;
  79   UShort_t nTPCClus = 0;
  80   if (pidObj) {
  81       nTPCClus = pidObj->GetTPCsignalN();
  82       dedx     = pidObj->GetTPCsignal();
  83       mom      = pidObj->GetTPCmomentum();
  84   }
  85
  86   Bool_t mismatch=kTRUE;
  87
  88   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
  89     AliPID::EParticleType type=AliPID::EParticleType(j);
  90     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
  91     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,nTPCClus,type);
  92     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
  93       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
  94     } else {
  95       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
  96       mismatch=kFALSE;
  97     }
  98
  99   }
 100
 101   if (mismatch)
 102     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
 103
 104
 105   return;
 106 }
 107 //_________________________________________________________________________
 108 void AliAODpidUtil::MakeITSPID(AliAODTrack *track,Double_t *p) const
 109 {
 110   //
 111   // ITS PID
 112   //  1) Truncated mean method
 113   //
 114
 115
 116   if ((track->GetStatus()&AliESDtrack::kITSin)==0) return;
 117   UChar_t clumap=track->GetITSClusterMap();
 118   Int_t nPointsForPid=0;
 119   for(Int_t i=2; i<6; i++){
 120    if(clumap&(1<<i)) ++nPointsForPid;
 121   }
 122   if(nPointsForPid<3) { // track not to be used for combined PID purposes
 123     for (Int_t j=0; j<AliPID::kSPECIES; j++)
 124       p[j] = 1./AliPID::kSPECIES;
 125     return;
 126   }
 127   Double_t mom=track->P();
 128   AliAODPid *pidObj = track->GetDetPid();
 129
 130   Double_t dedx = 0.;
 131   if (pidObj) {
 132       dedx = pidObj->GetITSsignal();
 133   }
 134
 135   Bool_t mismatch = kTRUE;
 136   Bool_t isSA = kTRUE;
 137   if(track->GetStatus() & AliESDtrack::kTPCin){
 138     isSA = kFALSE;
 139     if (pidObj)
 140       mom = pidObj->GetTPCmomentum();
 141   }
 142   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 143     Double_t mass = AliPID::ParticleMass(j);//GeV/c^2
 144     Double_t bethe = fITSResponse.Bethe(mom,mass);
 145     Double_t sigma = fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
 146     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
 147       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
 148     } else {
 149       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
 150       mismatch=kFALSE;
 151     }
 152
 153     // Check for particles heavier than (AliPID::kSPECIES - 1)
 154
 155   }
 156
 157   if (mismatch)
 158     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
 159
 160   return;
 161
 162 }
 163 //_________________________________________________________________________
 164 void AliAODpidUtil::MakeTOFPID(AliAODTrack *track, Double_t *p) const
 165 {
 166   //
 167   //   TOF PID using gaussian response
 168   //
 169   if ((track->GetStatus()&AliESDtrack::kTOFout )==0)    return;
 170   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
 171   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
 172
 173   Double_t time[AliPID::kSPECIESN];
 174   Double_t sigma[AliPID::kSPECIESN];
 175   AliAODPid *pidObj = track->GetDetPid();
 176   pidObj->GetIntegratedTimes(time);
 177
 178   for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
 179     sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
 180   }
 181
 182   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 183                   Form("Expected TOF signals [ps]: %f %f %f %f %f",
 184                        time[AliPID::kElectron],
 185                        time[AliPID::kMuon],
 186                        time[AliPID::kPion],
 187                        time[AliPID::kKaon],
 188                        time[AliPID::kProton]));
 189
 190   AliDebugGeneral("AliESDpid::MakeTOFPID",2,
 191                   Form("Expected TOF std deviations [ps]: %f %f %f %f %f",
 192                        sigma[AliPID::kElectron],
 193                        sigma[AliPID::kMuon],
 194                        sigma[AliPID::kPion],
 195                        sigma[AliPID::kKaon],
 196                        sigma[AliPID::kProton]
 197                        ));
 198
 199   Double_t tof = pidObj->GetTOFsignal();
 200
 201   Bool_t mismatch = kTRUE;
 202   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
 203     Double_t sig = sigma[j];
 204     if (TMath::Abs(tof-time[j]) > fRange*sig) {
 205       p[j] = TMath::Exp(-0.5*fRange*fRange)/sig;
 206     } else
 207       p[j] = TMath::Exp(-0.5*(tof-time[j])*(tof-time[j])/(sig*sig))/sig;
 208
 209     // Check the mismatching
 210     Double_t mass = AliPID::ParticleMass(j);
 211     Double_t pm = fTOFResponse.GetMismatchProbability(track->P(),mass);
 212     if (p[j]>pm) mismatch = kFALSE;
 213
 214     // Check for particles heavier than (AliPID::kSPECIES - 1)
 215
 216   }
 217
 218   if (mismatch)
 219     for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1/AliPID::kSPECIES;
 220
 221   return;
 222 }
 223 //_________________________________________________________________________
 224 void AliAODpidUtil::MakeTRDPID(AliAODTrack *track,Double_t *p) const
 225 {
 226
 227   // Method to recalculate the TRD PID probabilities
 228   if ((track->GetStatus()&AliESDtrack::kTRDout )==0)   return;
 229
 230   AliAODPid *pidObj = track->GetDetPid();
 231   Float_t *mom=pidObj->GetTRDmomentum();
 232   Int_t ntracklets=0;
 233   for(Int_t iPl=0;iPl<6;iPl++){
 234    if(mom[iPl]>0.) ntracklets++;
 235   }
 236    if(ntracklets<4) return;
 237
 238   Double_t* dedx=pidObj->GetTRDsignal();
 239   Bool_t norm=kTRUE;
 240   fTRDResponse.GetResponse(pidObj->GetTRDnSlices(),dedx,mom,p,norm);
 241   return;
 242 }