PWGHF/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 "AliAODPid.h"
  32 #include "AliAODTrack.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliESDtrack.h"
  35 #include "AliExternalTrackParam.h"
  36 #include "AliLog.h"
  37 #include "AliMCParticle.h"
  38 #include "AliPID.h"
  39 #include "AliPIDResponse.h"
  40
  41 #include "AliHFEpidTPC.h"
  42 #include "AliHFEpidQAmanager.h"
  43
  44 ClassImp(AliHFEpidTPC)
  45
  46 //___________________________________________________________________
  47 AliHFEpidTPC::AliHFEpidTPC() :
  48   // add a list here
  49   AliHFEpidBase()
  50   , fLineCrossingsEnabled(0)
  51   , fkEtaCorrection(NULL)
  52   , fHasCutModel(kFALSE)
  53   , fNsigmaTPC(3)
  54   , fRejectionEnabled(0)
  55 {
  56   //
  57   // default  constructor
  58   //
  59
  60   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
  61   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
  62
  63   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  64   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  65   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  66   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  67
  68 }
  69
  70 //___________________________________________________________________
  71 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  72   // add a list here
  73   AliHFEpidBase(name)
  74   , fLineCrossingsEnabled(0)
  75   , fkEtaCorrection(NULL)
  76   , fHasCutModel(kFALSE)
  77   , fNsigmaTPC(3)
  78   , fRejectionEnabled(0)
  79 {
  80   //
  81   // default  constructor
  82   //
  83   //
  84   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
  85   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
  86
  87   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  88   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  89   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  90   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  91 }
  92
  93 //___________________________________________________________________
  94 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  95   AliHFEpidBase("")
  96   , fLineCrossingsEnabled(0)
  97   , fkEtaCorrection(NULL)
  98   , fHasCutModel(ref.fHasCutModel)
  99   , fNsigmaTPC(2)
 100   , fRejectionEnabled(0)
 101 {
 102   //
 103   // Copy constructor
 104   //
 105   ref.Copy(*this);
 106 }
 107
 108 //___________________________________________________________________
 109 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
 110   //
 111   // Assignment operator
 112   //
 113   if(this != &ref){
 114     ref.Copy(*this);
 115   }
 116   return *this;
 117 }
 118 //___________________________________________________________________
 119 void AliHFEpidTPC::Copy(TObject &o) const{
 120   //
 121   // Copy function
 122   // called in copy constructor and assigment operator
 123   //
 124   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 125
 126   target.fkEtaCorrection =fkEtaCorrection;
 127   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
 128   target.fHasCutModel = fHasCutModel;
 129   target.fNsigmaTPC = fNsigmaTPC;
 130   target.fRejectionEnabled = fRejectionEnabled;
 131
 132   memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
 133   memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
 134
 135   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
 136   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
 137   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
 138
 139   AliHFEpidBase::Copy(target);
 140 }
 141
 142 //___________________________________________________________________
 143 AliHFEpidTPC::~AliHFEpidTPC(){
 144   //
 145   // Destructor
 146   //
 147 }
 148
 149 //___________________________________________________________________
 150 Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){
 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(!fkPIDResponse) return 0;
 169
 170   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 171
 172   // Make clone track in order to be able to apply correction
 173   const AliVTrack *rectrack;
 174   AliESDtrack esdtrack;
 175   AliAODTrack aodtrack;
 176   if(fkEtaCorrection){
 177     // Correction available
 178     // apply it on copy
 179     if(track->IsESDanalysis()){
 180       esdtrack.~AliESDtrack();
 181       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
 182       ApplyEtaCorrection(&esdtrack, anatype);
 183       rectrack = &esdtrack;
 184     } else {
 185       aodtrack.~AliAODTrack();
 186       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
 187       ApplyEtaCorrection(&aodtrack, anatype);
 188       rectrack = &aodtrack;
 189     }
 190   } else {
 191     // Correction not available - no need to copy
 192     rectrack = track->GetRecTrack();
 193   }
 194   AliHFEpidObject tpctrack(*track);
 195   tpctrack.SetRecTrack(rectrack);
 196
 197   // QA before selection (after correction)
 198   if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
 199   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
 200
 201   // make copy of the track in order to allow for applying the correction
 202   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
 203   AliDebug(1, Form("TPC NSigma: %f", nsigma));
 204   // exclude crossing points:
 205   // Determine the bethe values for each particle species
 206   Bool_t isLineCrossing = kFALSE;
 207   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 208     if(ispecies == AliPID::kElectron) continue;
 209     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 210     if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
 211       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
 212       isLineCrossing = kTRUE;
 213       break;
 214     }
 215   }
 216   if(isLineCrossing) return 0;
 217
 218   // Check particle rejection
 219   if(HasParticleRejection()){
 220     Int_t reject = Reject(rectrack, anatype);
 221     if(reject != 0) return reject;
 222   }
 223
 224   // Check if we have an asymmetric sigma model set
 225   Int_t pdg = 0;
 226   if(fHasCutModel){
 227     pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
 228   } else {
 229     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
 230     Float_t p = rectrack->P();
 231     if(HasAsymmetricSigmaCut()) {
 232
 233       //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]);
 234       if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) {
 235         if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
 236         else pdg = 0;
 237       }
 238       else pdg = 0;
 239
 240     }
 241     else {
 242       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
 243     }
 244   }
 245   if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
 246   return pdg;
 247
 248 }
 249
 250 //___________________________________________________________________
 251 Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
 252   //
 253   // N SigmaCut using parametrization of the cuts
 254   //
 255   Bool_t isSelected = kTRUE;
 256   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 257   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
 258   Double_t p = GetP(track->GetRecTrack(), anatype);
 259   Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
 260   AliDebug(2, Form("Centrality: %d\n", centrality));
 261   if(centrality > 11) return kFALSE;
 262   const TF1 *cutfunction;
 263   if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
 264   if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
 265   return isSelected;
 266 }
 267
 268 //___________________________________________________________________
 269 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
 270   //
 271   // reject particles based on asymmetric sigma cut
 272   //
 273   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
 274   Double_t p = GetP(track, anaType);
 275   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 276     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
 277     // Particle rejection enabled
 278     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
 279     Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
 280     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
 281   }
 282   return 0;
 283 }
 284
 285 //___________________________________________________________________
 286 void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
 287   //
 288   // Apply correction for the eta dependence
 289   // N.B. This correction has to be applied on a copy track
 290   //
 291   AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
 292   Double_t original = track->GetTPCsignal(),
 293            eta = track->Eta();
 294   if(anatype == AliHFEpidObject::kESDanalysis){
 295     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 296     esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
 297   } else {
 298     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 299     AliAODPid *pid = aodtrack->GetDetPid();
 300     if(pid) pid->SetTPCsignal(original);
 301   }
 302 }
 303
 304 //___________________________________________________________________
 305 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
 306   //
 307   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 308   // Stores line center and line sigma
 309   //
 310   if(species >= AliPID::kSPECIES){
 311     AliError("Species doesn't exist");
 312     return;
 313   }
 314   fLineCrossingsEnabled |= 1 << species;
 315   fLineCrossingSigma[species] = sigma;
 316 }
 317
 318 //___________________________________________________________________
 319 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
 320   //
 321   // Get the momentum at the inner wall of the TPC
 322   //
 323   Double_t p = -1;
 324   if(anatype == AliHFEpidObject::kESDanalysis){
 325     // ESD analysis: Use Inner Params for the momentum estimate
 326     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 327     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
 328   } else {
 329     // AOD analysis: Use TPC momentum stored in the AliAODpid object
 330     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 331     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
 332   }
 333   return p;
 334 }