PWG3/hfe/AliHFEpidTPC.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 // Class for TPC PID
  17 // Implements the abstract base class AliHFEpidBase
  18 //
  19 // Class contains TPC specific cuts and QA histograms
  20 // Two selection strategies are offered: Selection of certain value
  21 // regions in the TPC dE/dx (by IsSelected), and likelihoods
  22 //
  23 // Authors:
  24 //
  25 //   Markus Fasel <M.Fasel@gsi.de>
  26 //   Markus Heide <mheide@uni-muenster.de>
  27 //
  28 #include <TF1.h>
  29 #include <TMath.h>
  30
  31 #include "AliAODpidUtil.h"
  32 #include "AliAODPid.h"
  33 #include "AliAODTrack.h"
  34 #include "AliAODMCParticle.h"
  35 #include "AliESDtrack.h"
  36 #include "AliExternalTrackParam.h"
  37 #include "AliLog.h"
  38 #include "AliMCParticle.h"
  39 #include "AliPID.h"
  40 #include "AliESDpid.h"
  41
  42 #include "AliHFEpidTPC.h"
  43 #include "AliHFEpidQAmanager.h"
  44
  45 ClassImp(AliHFEpidTPC)
  46
  47 //___________________________________________________________________
  48 AliHFEpidTPC::AliHFEpidTPC() :
  49   // add a list here
  50   AliHFEpidBase()
  51   , fLineCrossingsEnabled(0)
  52   , fUpperSigmaCut(NULL)
  53   , fLowerSigmaCut(NULL)
  54   , fElectronMeanCorrection(NULL)
  55   , fNsigmaTPC(3)
  56   , fRejectionEnabled(0)
  57   , fPID(NULL)
  58 {
  59   //
  60   // default  constructor
  61   //
  62   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  63   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  64   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  65   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  66
  67 }
  68
  69 //___________________________________________________________________
  70 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  71   // add a list here
  72   AliHFEpidBase(name)
  73   , fLineCrossingsEnabled(0)
  74   , fUpperSigmaCut(NULL)
  75   , fLowerSigmaCut(NULL)
  76   , fElectronMeanCorrection(NULL)
  77   , fNsigmaTPC(3)
  78   , fRejectionEnabled(0)
  79   , fPID(NULL)
  80 {
  81   //
  82   // default  constructor
  83   //
  84   //
  85   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  86   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  87   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  88   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  89   fPID = new AliPID;
  90 }
  91
  92 //___________________________________________________________________
  93 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  94   AliHFEpidBase("")
  95   , fLineCrossingsEnabled(0)
  96   , fUpperSigmaCut(NULL)
  97   , fLowerSigmaCut(NULL)
  98   , fElectronMeanCorrection(NULL)
  99   , fNsigmaTPC(2)
 100   , fRejectionEnabled(0)
 101   , fPID(NULL)
 102 {
 103   //
 104   // Copy constructor
 105   //
 106   ref.Copy(*this);
 107 }
 108
 109 //___________________________________________________________________
 110 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
 111   //
 112   // Assignment operator
 113   //
 114   if(this != &ref){
 115     ref.Copy(*this);
 116   }
 117   return *this;
 118 }
 119 //___________________________________________________________________
 120 void AliHFEpidTPC::Copy(TObject &o) const{
 121   //
 122   // Copy function
 123   // called in copy constructor and assigment operator
 124   //
 125   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 126
 127   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
 128   target.fUpperSigmaCut = fUpperSigmaCut;
 129   target.fLowerSigmaCut = fLowerSigmaCut;
 130   target.fElectronMeanCorrection = fElectronMeanCorrection;
 131   target.fNsigmaTPC = fNsigmaTPC;
 132   target.fRejectionEnabled = fRejectionEnabled;
 133   target.fPID = new AliPID(*fPID);
 134   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
 135   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
 136   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
 137
 138   AliHFEpidBase::Copy(target);
 139 }
 140
 141 //___________________________________________________________________
 142 AliHFEpidTPC::~AliHFEpidTPC(){
 143   //
 144   // Destructor
 145   //
 146   if(fPID) delete fPID;
 147 }
 148
 149 //___________________________________________________________________
 150 Bool_t AliHFEpidTPC::InitializePID(){
 151   //
 152   // Add TPC dE/dx Line crossings
 153   //
 154   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
 155   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
 156   return kTRUE;
 157 }
 158
 159 //___________________________________________________________________
 160 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
 161 {
 162   //
 163   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
 164   // for electrons
 165   // exclusion of the crossing points
 166   //
 167
 168   if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())) return 0;
 169
 170   // QA before selection
 171   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
 172   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
 173   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 174   Float_t nsigma = NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
 175   AliDebug(1, Form("TPC NSigma: %f", nsigma));
 176   // exclude crossing points:
 177   // Determine the bethe values for each particle species
 178   Bool_t isLineCrossing = kFALSE;
 179   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 180     if(ispecies == AliPID::kElectron) continue;
 181     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 182     if(TMath::Abs(NumberOfSigmas(track->GetRecTrack(), (AliPID::EParticleType)ispecies, anatype)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
 183       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
 184       isLineCrossing = kTRUE;
 185       break;
 186     }
 187   }
 188   if(isLineCrossing) return 0;
 189
 190   // Check particle rejection
 191   if(HasParticleRejection()){
 192     Int_t reject = Reject(track->GetRecTrack(), anatype);
 193     if(reject != 0) return reject;
 194   }
 195
 196   // Check if we have an asymmetric sigma model set
 197   Int_t pdg = 0;
 198   if(fUpperSigmaCut || fLowerSigmaCut){
 199     pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
 200   } else {
 201     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
 202     Float_t p = 0.;
 203     if(HasAsymmetricSigmaCut() && (p = track->GetRecTrack()->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
 204       if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
 205     } else {
 206       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
 207     }
 208   }
 209   if(pidqa && pdg != 0) pidqa->ProcessTrack(track, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
 210   return pdg;
 211
 212 }
 213
 214 //___________________________________________________________________
 215 Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 216   //
 217   // N SigmaCut using parametrization of the cuts
 218   //
 219   Bool_t isSelected = kTRUE;
 220   Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
 221   Double_t p = GetP(track, anaType);
 222   if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
 223   if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
 224   return isSelected;
 225 }
 226
 227 //___________________________________________________________________
 228 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
 229   //
 230   // reject particles based on asymmetric sigma cut
 231   //
 232   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
 233   Double_t p = GetP(track, anaType);
 234   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 235     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
 236     // Particle rejection enabled
 237     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
 238     Double_t sigma = NumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec), anaType);
 239     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
 240   }
 241   return 0;
 242 }
 243
 244 //___________________________________________________________________
 245 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
 246   //
 247   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 248   // Stores line center and line sigma
 249   //
 250   if(species >= AliPID::kSPECIES){
 251     AliError("Species doesn't exist");
 252     return;
 253   }
 254   fLineCrossingsEnabled |= 1 << species;
 255   fLineCrossingSigma[species] = sigma;
 256 }
 257
 258 //___________________________________________________________________
 259 Double_t AliHFEpidTPC::NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anaType) const {
 260   //
 261   // Get the number of sigmas
 262   //
 263   Double_t nSigmas = 100;
 264   if(anaType == AliHFEpidObject::kESDanalysis){
 265     // ESD analysis
 266     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 267     if(esdtrack && fESDpid) nSigmas = fESDpid->NumberOfSigmasTPC(esdtrack, species);
 268   } else {
 269     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 270     if(aodtrack && fAODpid) nSigmas = fAODpid->NumberOfSigmasTPC(aodtrack, species);
 271   }
 272   // Correct for the mean o
 273   if(fElectronMeanCorrection)
 274     nSigmas -= fElectronMeanCorrection->Eval(GetP(track, anaType));
 275   return nSigmas;
 276 }
 277
 278 //___________________________________________________________________
 279 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
 280   //
 281   // Get the momentum at the inner wall of the TPC
 282   //
 283   Double_t p = -1;
 284   if(anatype == AliHFEpidObject::kESDanalysis){
 285     // ESD analysis: Use Inner Params for the momentum estimate
 286     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 287     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
 288   } else {
 289     // AOD analysis: Use TPC momentum stored in the AliAODpid object
 290     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 291     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
 292   }
 293   return p;
 294 }