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