PWG3/vertexingHF/AliAODPidHF.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  *                                                                        *
   5  * Author: The ALICE Off-line Project.                                    *
   6  * Contributors are mentioned in the code where appropriate.              *
   7  *
   8  * Permission to use, copy, modify and distribute this software and its   *
   9  * documentation strictly for non-commercial purposes is hereby granted   *
  10  * without fee, provided that the above copyright notice appears in all   *
  11  * copies and that both the copyright notice and this permission notice   *
  12  * appear in the supporting documentation. The authors make no claims     *
  13  * about the suitability of this software for any purpose. It is          *
  14  * provided "as is" without express or implied warranty.                  *
  15  * *************************************************************************/
  16
  17 //***********************************************************
  18 // Class AliAODPidHF
  19 // class for PID with AliAODRecoDecayHF
  20 // Authors: D. Caffarri caffarri@bo.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de
  21 //***********************************************************
  22 #include "AliAODPidHF.h"
  23 #include "AliAODPid.h"
  24 #include "AliAODTrack.h"
  25 #include "AliPID.h"
  26 #include "AliTPCPIDResponse.h"
  27 #include "AliITSPIDResponse.h"
  28 #include "AliTOFPIDResponse.h"
  29 #include "AliAODpidUtil.h"
  30
  31
  32 ClassImp(AliAODPidHF)
  33
  34 //------------------------------
  35 AliAODPidHF::AliAODPidHF():
  36   fSigma(3.),
  37   fPriors()
  38 {
  39  //
  40  // Default constructor
  41  //
  42
  43 }
  44 //----------------------
  45 AliAODPidHF::~AliAODPidHF()
  46 {
  47       // destructor
  48 }
  49 //------------------------
  50 AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
  51   AliAODPid(pid),
  52   fSigma(pid.fSigma),
  53   fPriors(pid.fPriors)
  54   {
  55
  56   for(Int_t i=0;i<5;i++){
  57     fPriors[i]=pid.fPriors[i];
  58   }
  59
  60   }
  61
  62 //----------------------
  63 Int_t AliAODPidHF::RawSignalPID(AliAODTrack *track, TString detector){
  64
  65    if(detector.Contains("ITS")) return ApplyPidITSRaw(track,0);
  66    if(detector.Contains("TPC")) return ApplyPidTPCRaw(track,0);
  67    if(detector.Contains("TOF")) return ApplyPidTOFRaw(track,0);
  68
  69   return 0;
  70
  71 }
  72 //---------------------------
  73 Bool_t AliAODPidHF::IsKaonRaw (AliAODTrack *track, TString detector){
  74
  75  Int_t specie=0;
  76
  77  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,3);
  78  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,3);
  79  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,3);
  80
  81  if(specie==3) return kTRUE;
  82  return kFALSE;
  83 }
  84 //---------------------------
  85 Bool_t AliAODPidHF::IsPionRaw (AliAODTrack *track, TString detector){
  86
  87  Int_t specie=0;
  88
  89  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,2);
  90  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,2);
  91  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,2);
  92
  93  if(specie==2) return kTRUE;
  94  return kFALSE;
  95 }
  96 //---------------------------
  97 Bool_t AliAODPidHF::IsProtonRaw (AliAODTrack *track, TString detector){
  98
  99  Int_t specie=0;
 100  if(detector.Contains("ITS")) specie=ApplyPidITSRaw(track,4);
 101  if(detector.Contains("TPC")) specie=ApplyPidTPCRaw(track,4);
 102  if(detector.Contains("TOF")) specie=ApplyPidTOFRaw(track,4);
 103
 104  if(specie==4) return kTRUE;
 105
 106  return kFALSE;
 107 }
 108 //---------------------------
 109 Int_t AliAODPidHF::ApplyPidTPCRaw(AliAODTrack *track,Int_t specie){
 110
 111   if(track->P()>2.) printf("Warning: TPC resolution not good in this region");
 112   AliAODPid *pidObj = track->GetDetPid();
 113
 114   Double_t dedx=pidObj->GetTPCsignal();
 115   Double_t mom = pidObj->GetTPCmomentum();
 116   AliTPCPIDResponse tpcResponse;
 117   Int_t pid=0;
 118   if(specie==0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
 119    Double_t nsigmaMax=fSigma;
 120    for(Int_t ipart=2;ipart<5;ipart++){
 121     AliPID::EParticleType type=AliPID::EParticleType(ipart);
 122     Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,50,type));
 123     if((nsigma<nsigmaMax) && (nsigma<fSigma)) {
 124      pid=ipart;
 125      nsigmaMax=nsigma;
 126     }
 127    }
 128   }else{ // asks only for one particle specie
 129    AliPID::EParticleType type=AliPID::EParticleType(specie);
 130     Double_t nsigma = TMath::Abs(tpcResponse.GetNumberOfSigmas(mom,dedx,50,type));
 131    if (nsigma>fSigma) {
 132     pid=0;
 133    }else{
 134     pid=specie;
 135    }
 136   }
 137
 138  return pid;
 139
 140 }
 141 //----------------------------
 142 Int_t AliAODPidHF::ApplyPidITSRaw(AliAODTrack *track,Int_t specie){
 143   Double_t mom=track->P();
 144   AliAODPid *pidObj = track->GetDetPid();
 145
 146   Double_t dedx=pidObj->GetITSsignal();
 147   AliITSPIDResponse itsResponse;
 148   Int_t pid=0;
 149   if(specie==0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
 150    Double_t nsigmaMax=fSigma;
 151    for(Int_t ipart=2;ipart<5;ipart++){
 152     AliPID::EParticleType type=AliPID::EParticleType(ipart);
 153     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
 154     if((nsigma<nsigmaMax) && (nsigma<fSigma)) {
 155      pid=ipart;
 156      nsigmaMax=nsigma;
 157     }
 158    }
 159   }else{ // asks only for one particle specie
 160    AliPID::EParticleType type=AliPID::EParticleType(specie);
 161     Double_t nsigma = TMath::Abs(itsResponse.GetNumberOfSigmas(mom,dedx,type));
 162    if (nsigma>fSigma) {
 163     pid=0;
 164    }else{
 165     pid=specie;
 166    }
 167   }
 168  return pid;
 169 }
 170 //----------------------------
 171 Int_t AliAODPidHF::ApplyPidTOFRaw(AliAODTrack *track,Int_t specie){
 172  Double_t time[AliPID::kSPECIESN];
 173  AliAODPid *pidObj = track->GetDetPid();
 174  pidObj->GetIntegratedTimes(time);
 175  AliTOFPIDResponse tofResponse;
 176  Int_t pid=0;
 177
 178   if(specie==0){  // from RawSignalPID : should return the particle specie to wich the de/dx is closer to the bethe-block curve -> performance to be checked
 179    Double_t nsigmaMax=fSigma;
 180    for(Int_t ipart=2;ipart<5;ipart++){
 181     AliPID::EParticleType type=AliPID::EParticleType(ipart);
 182     Double_t nsigma = tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type));
 183     if((nsigma<nsigmaMax) && (nsigma<fSigma)) {
 184      pid=ipart;
 185      nsigmaMax=nsigma;
 186     }
 187    }
 188   }else{ // asks only for one particle specie
 189    AliPID::EParticleType type=AliPID::EParticleType(specie);
 190    Double_t nsigma = TMath::Abs(tofResponse.GetExpectedSigma(track->P(),time[type],AliPID::ParticleMass(type)));
 191    if (nsigma>fSigma) {
 192     pid=0;
 193    }else{
 194     pid=specie;
 195    }
 196   }
 197  return pid;
 198
 199 }
 200 //------------------------------
 201 void AliAODPidHF::CombinedProbability(AliAODTrack *track,Bool_t *type){
 202
 203  const Double_t *pid=track->PID();
 204  Float_t max=0.;
 205  Int_t k=-1;
 206  for (Int_t i=0; i<10; i++) {
 207    if (pid[i]>max) {k=i; max=pid[i];}
 208  }
 209
 210  if(k==2) type[0]=kTRUE;
 211  if(k==3) type[1]=kTRUE;
 212  if(k==4) type[2]=kTRUE;
 213
 214  return;
 215 }
 216 //--------------------
 217 void AliAODPidHF::BayesianProbability(AliAODTrack *track,TString detectors,Double_t *pid){
 218
 219   if(detectors.Contains("ITS")) {BayesianProbabilityITS(track,pid);return;}
 220   if(detectors.Contains("TPC")) {BayesianProbabilityTPC(track,pid);return;}
 221   if(detectors.Contains("TOF")) {BayesianProbabilityTOF(track,pid);return;}
 222
 223   if(detectors.Contains("All")) {
 224     Double_t probITS[5]={0.,0.,0.,0.,0.};
 225     Double_t probTPC[5]={0.,0.,0.,0.,0.};
 226     Double_t probTOF[5]={0.,0.,0.,0.,0.};
 227     BayesianProbabilityITS(track,probITS);
 228     BayesianProbabilityTPC(track,probTPC);
 229     BayesianProbabilityTOF(track,probTOF);
 230     Double_t probTot[5]={0.,0.,0.,0.,0.};
 231     for(Int_t i=0;i<5;i++){
 232      probTot[i]=probITS[i]*probTPC[i]*probTOF[i];
 233     }
 234     for(Int_t i2=0;i2<5;i2++){
 235      pid[i2]=probTot[i2]*fPriors[i2]/(probTot[0]*fPriors[0]+probTot[1]*fPriors[1]+probTot[2]*fPriors[2]+probTot[3]*fPriors[3]+probTot[4]*fPriors[4]);
 236     }
 237    }
 238
 239  return;
 240
 241 }
 242 //------------------------------------
 243 void AliAODPidHF::BayesianProbabilityITS(AliAODTrack *track,Double_t *prob){
 244
 245  AliAODpidUtil pid;
 246  Double_t itspid[AliPID::kSPECIES];
 247  pid.MakeITSPID(track,itspid);
 248  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 249   prob[ind]=itspid[ind]*fPriors[ind]/(itspid[0]*fPriors[0]+itspid[1]*fPriors[1]+itspid[2]*fPriors[2]+itspid[3]*fPriors[3]+itspid[4]*fPriors[4]);
 250  }
 251  return;
 252
 253 }
 254 //------------------------------------
 255 void AliAODPidHF::BayesianProbabilityTPC(AliAODTrack *track,Double_t *prob){
 256
 257  AliAODpidUtil pid;
 258  Double_t tpcpid[AliPID::kSPECIES];
 259  pid.MakeTPCPID(track,tpcpid);
 260  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 261   if(tpcpid[ind]>0.) {
 262    prob[ind]=tpcpid[ind]*fPriors[ind]/(tpcpid[0]*fPriors[0]+tpcpid[1]*fPriors[1]+tpcpid[2]*fPriors[2]+tpcpid[3]*fPriors[3]+tpcpid[4]*fPriors[4]);
 263   }else{
 264    prob[ind]=0.;
 265   }
 266  }
 267  return;
 268
 269 }
 270 //------------------------------------
 271 void AliAODPidHF::BayesianProbabilityTOF(AliAODTrack *track,Double_t *prob){
 272
 273  AliAODpidUtil pid;
 274  Double_t tofpid[AliPID::kSPECIES];
 275  Float_t time=0.1; //needs to be changed!
 276  pid.MakeTOFPID(track,time,tofpid);
 277  for(Int_t ind=0;ind<AliPID::kSPECIES;ind++){
 278   prob[ind]=tofpid[ind]*fPriors[ind]/(tofpid[0]*fPriors[0]+tofpid[1]*fPriors[1]+tofpid[2]*fPriors[2]+tofpid[3]*fPriors[3]+tofpid[4]*fPriors[4]);
 279  }
 280  return;
 281
 282 }