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 "AliTPCdEdxInfo.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 "AliPIDResponse.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   , fkEtaCorrection(NULL)
  53   , fkCentralityCorrection(NULL)
  54   , fHasCutModel(kFALSE)
  55   , fUseOnlyOROC(kFALSE)
  56   , fNsigmaTPC(3)
  57   , fRejectionEnabled(0)
  58 {
  59   //
  60   // default  constructor
  61   //
  62
  63   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
  64   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
  65
  66   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  67   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  68   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  69   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  70
  71 }
  72
  73 //___________________________________________________________________
  74 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  75   // add a list here
  76   AliHFEpidBase(name)
  77   , fLineCrossingsEnabled(0)
  78   , fkEtaCorrection(NULL)
  79   , fkCentralityCorrection(NULL)
  80   , fHasCutModel(kFALSE)
  81   , fUseOnlyOROC(kFALSE)
  82   , fNsigmaTPC(3)
  83   , fRejectionEnabled(0)
  84 {
  85   //
  86   // default  constructor
  87   //
  88   //
  89   memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
  90   memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
  91
  92   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
  93   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  94   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
  95   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
  96 }
  97
  98 //___________________________________________________________________
  99 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
 100   AliHFEpidBase("")
 101   , fLineCrossingsEnabled(0)
 102   , fkEtaCorrection(NULL)
 103   , fkCentralityCorrection(NULL)
 104   , fHasCutModel(ref.fHasCutModel)
 105   , fUseOnlyOROC(ref.fUseOnlyOROC)
 106   , fNsigmaTPC(2)
 107   , fRejectionEnabled(0)
 108 {
 109   //
 110   // Copy constructor
 111   //
 112   ref.Copy(*this);
 113 }
 114
 115 //___________________________________________________________________
 116 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
 117   //
 118   // Assignment operator
 119   //
 120   if(this != &ref){
 121     ref.Copy(*this);
 122   }
 123   return *this;
 124 }
 125 //___________________________________________________________________
 126 void AliHFEpidTPC::Copy(TObject &o) const{
 127   //
 128   // Copy function
 129   // called in copy constructor and assigment operator
 130   //
 131   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 132
 133   target.fkEtaCorrection = fkEtaCorrection;
 134   target.fkCentralityCorrection = fkCentralityCorrection;
 135   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
 136   target.fHasCutModel = fHasCutModel;
 137   target.fUseOnlyOROC = fUseOnlyOROC;
 138   target.fNsigmaTPC = fNsigmaTPC;
 139   target.fRejectionEnabled = fRejectionEnabled;
 140
 141   memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
 142   memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
 143
 144   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
 145   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
 146   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
 147
 148   AliHFEpidBase::Copy(target);
 149 }
 150
 151 //___________________________________________________________________
 152 AliHFEpidTPC::~AliHFEpidTPC(){
 153   //
 154   // Destructor
 155   //
 156 }
 157
 158 //___________________________________________________________________
 159 Bool_t AliHFEpidTPC::InitializePID(Int_t /*run*/){
 160   //
 161   // Add TPC dE/dx Line crossings
 162   //
 163   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
 164   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
 165   return kTRUE;
 166 }
 167
 168 //___________________________________________________________________
 169 Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
 170 {
 171   //
 172   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
 173   // for electrons
 174   // exclusion of the crossing points
 175   //
 176
 177   if(!fkPIDResponse) return 0;
 178
 179   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 180
 181   // Make clone track in order to be able to apply correction
 182   const AliVTrack *rectrack;
 183   AliESDtrack esdtrack;
 184   AliAODTrack aodtrack;
 185   /*if(fUseOnlyOROC && !(fkEtaCorrection || fkCentralityCorrection)) {
 186     if(track->IsESDanalysis()){
 187       esdtrack.~AliESDtrack();
 188       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
 189       UseOROC(&esdtrack, anatype);
 190       rectrack = &esdtrack;
 191     } else {
 192       aodtrack.~AliAODTrack();
 193       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
 194       UseOROC(&aodtrack, anatype);
 195       rectrack = &aodtrack;
 196     }
 197   }
 198   else if(fkEtaCorrection || fkCentralityCorrection){*/
 199   if(fkEtaCorrection || fkCentralityCorrection){
 200     // Correction available
 201     // apply it on copy
 202     if(track->IsESDanalysis()){
 203       esdtrack.~AliESDtrack();
 204       new(&esdtrack) AliESDtrack(*(static_cast<const AliESDtrack *>(track->GetRecTrack())));
 205       if(track->IsPbPb() && HasCentralityCorrection())
 206         ApplyCentralityCorrection(&esdtrack, track->GetMultiplicity(), anatype);
 207       if(HasEtaCorrection())
 208         ApplyEtaCorrection(&esdtrack, anatype);
 209       rectrack = &esdtrack;
 210     } else {
 211       aodtrack.~AliAODTrack();
 212       new(&aodtrack) AliAODTrack(*(static_cast<const AliAODTrack *>(track->GetRecTrack())));
 213       if(track->IsPbPb() && HasCentralityCorrection())
 214         ApplyCentralityCorrection(&aodtrack, track->GetMultiplicity(), anatype);
 215       if(HasEtaCorrection())
 216         ApplyEtaCorrection(&aodtrack, anatype);
 217       rectrack = &aodtrack;
 218     }
 219   } else {
 220     // Correction not available - no need to copy
 221     rectrack = track->GetRecTrack();
 222   }
 223   AliHFEpidObject tpctrack(*track);
 224   tpctrack.SetRecTrack(rectrack);
 225
 226   // QA before selection (after correction)
 227   if(pidqa) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kBeforePID);
 228   AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
 229
 230   // make copy of the track in order to allow for applying the correction
 231   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(rectrack, AliPID::kElectron);
 232   AliDebug(1, Form("TPC NSigma: %f", nsigma));
 233   // exclude crossing points:
 234   // Determine the bethe values for each particle species
 235   Bool_t isLineCrossing = kFALSE;
 236   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 237     if(ispecies == AliPID::kElectron) continue;
 238     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 239     if(TMath::Abs(fkPIDResponse->NumberOfSigmasTPC(rectrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
 240       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
 241       isLineCrossing = kTRUE;
 242       break;
 243     }
 244   }
 245   if(isLineCrossing) return 0;
 246
 247   // Check particle rejection
 248   if(HasParticleRejection()){
 249     Int_t reject = Reject(rectrack, anatype);
 250     if(reject != 0) return reject;
 251   }
 252
 253   // Check if we have an asymmetric sigma model set
 254   Int_t pdg = 0;
 255   if(fHasCutModel){
 256     pdg = CutSigmaModel(&tpctrack) ? 11 : 0;
 257   } else {
 258     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
 259     Float_t p = rectrack->P();
 260     if(HasAsymmetricSigmaCut()) {
 261
 262       //printf("p %f, fPAsigCut[0] %f, fPAsigCut[1] %f\n",p,fPAsigCut[0],fPAsigCut[1]);
 263       if(p >= fPAsigCut[0] && p <= fPAsigCut[1]) {
 264         if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
 265         else pdg = 0;
 266       }
 267       else pdg = 0;
 268
 269     }
 270     else {
 271       if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
 272     }
 273   }
 274   if(pidqa && pdg != 0) pidqa->ProcessTrack(&tpctrack, AliHFEpid::kTPCpid, AliHFEdetPIDqa::kAfterPID);
 275   return pdg;
 276
 277 }
 278
 279 //___________________________________________________________________
 280 Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
 281   //
 282   // N SigmaCut using parametrization of the cuts
 283   //
 284   Bool_t isSelected = kTRUE;
 285   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
 286   Float_t nsigma = fkPIDResponse->NumberOfSigmasTPC(track->GetRecTrack(), AliPID::kElectron);
 287   Double_t p = GetP(track->GetRecTrack(), anatype);
 288   Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
 289   AliDebug(2, Form("Centrality: %d\n", centrality));
 290   if(centrality > 11) return kFALSE;
 291   const TF1 *cutfunction;
 292   if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
 293   if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
 294   return isSelected;
 295 }
 296
 297 //___________________________________________________________________
 298 Int_t AliHFEpidTPC::Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const{
 299   //
 300   // reject particles based on asymmetric sigma cut
 301   //
 302   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
 303   Double_t p = GetP(track, anaType);
 304   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 305     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
 306     // Particle rejection enabled
 307     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
 308     Double_t sigma = fkPIDResponse->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
 309     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
 310   }
 311   return 0;
 312 }
 313
 314 //___________________________________________________________________
 315 void AliHFEpidTPC::ApplyEtaCorrection(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
 316   //
 317   // Apply correction for the eta dependence
 318   // N.B. This correction has to be applied on a copy track
 319   //
 320   AliDebug(1, Form("Applying correction function %s\n", fkEtaCorrection->GetName()));
 321   Double_t original = track->GetTPCsignal(),
 322            eta = track->Eta();
 323   if(anatype == AliHFEpidObject::kESDanalysis){
 324     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 325     if(fkEtaCorrection->Eval(eta)>0.0) esdtrack->SetTPCsignal(original/fkEtaCorrection->Eval(eta), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
 326   } else {
 327     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 328     AliAODPid *pid = aodtrack->GetDetPid();
 329     if(pid && fkEtaCorrection->Eval(eta)>0.0) pid->SetTPCsignal(original/fkEtaCorrection->Eval(eta));
 330   }
 331 }
 332
 333 //___________________________________________________________________
 334 void AliHFEpidTPC::ApplyCentralityCorrection(AliVTrack *track, Double_t centralityEstimator, AliHFEpidObject::AnalysisType_t anatype) const{
 335   //
 336   // Apply correction for the eta dependence
 337   // N.B. This correction has to be applied on a copy track
 338   //
 339   AliDebug(1, Form("Applying correction function %s\n", fkCentralityCorrection->GetName()));
 340   Double_t original = track->GetTPCsignal();
 341   if(anatype == AliHFEpidObject::kESDanalysis){
 342     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 343     if(fkCentralityCorrection->Eval(centralityEstimator)>0.0) esdtrack->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator), esdtrack->GetTPCsignalSigma(), esdtrack->GetTPCsignalN());
 344   } else {
 345     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 346     AliAODPid *pid = aodtrack->GetDetPid();
 347     if(pid && fkCentralityCorrection->Eval(centralityEstimator)>0.0) pid->SetTPCsignal(original/fkCentralityCorrection->Eval(centralityEstimator));
 348   }
 349 }
 350
 351 //___________________________________________________________________
 352 void AliHFEpidTPC::UseOROC(AliVTrack *track, AliHFEpidObject::AnalysisType_t anatype) const{
 353   //
 354   // Use TPC signal from the OROC
 355   // N.B. This correction has to be applied on a copy track
 356   //
 357   //Double_t original = track->GetTPCsignal();
 358
 359   if(anatype == AliHFEpidObject::kESDanalysis){
 360     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
 361     AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
 362     Double32_t  TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
 363     Char_t      TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
 364     Char_t      TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
 365     dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
 366     esdtrack->SetTPCsignal(TPCsignalRegion[3],esdtrack->GetTPCsignalSigma(),(TPCsignalNRegion[1]+TPCsignalNRegion[2])); // the two last are not ok
 367   } else {
 368     AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
 369     AliTPCdEdxInfo *dEdxInfo = track->GetTPCdEdxInfo();
 370     Double32_t  TPCsignalRegion[4]; // TPC dEdx signal in 4 different regions - 0 - IROC, 1- OROC medium, 2 - OROC long, 3- OROC all, (default truncation used)
 371     Char_t      TPCsignalNRegion[3]; // number of clusters above threshold used in the dEdx calculation
 372     Char_t      TPCsignalNRowRegion[3]; // number of crosed rows used in the dEdx calculation - signal below threshold included
 373     dEdxInfo->GetTPCSignalRegionInfo(TPCsignalRegion,TPCsignalNRegion,TPCsignalNRowRegion);
 374     AliAODPid *pid = aodtrack->GetDetPid();
 375     if(pid) pid->SetTPCsignal(TPCsignalRegion[3]);
 376     if(pid) pid->SetTPCsignalN((TPCsignalNRegion[1]+TPCsignalNRegion[2]));
 377   }
 378
 379 }
 380 //___________________________________________________________________
 381 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
 382   //
 383   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 384   // Stores line center and line sigma
 385   //
 386   if(species >= AliPID::kSPECIES){
 387     AliError("Species doesn't exist");
 388     return;
 389   }
 390   fLineCrossingsEnabled |= 1 << species;
 391   fLineCrossingSigma[species] = sigma;
 392 }
 393
 394 //___________________________________________________________________
 395 Double_t AliHFEpidTPC::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype) const {
 396   //
 397   // Get the momentum at the inner wall of the TPC
 398   //
 399   Double_t p = -1;
 400   if(anatype == AliHFEpidObject::kESDanalysis){
 401     // ESD analysis: Use Inner Params for the momentum estimate
 402     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 403     if(esdtrack) p = esdtrack->GetInnerParam() ? esdtrack->GetInnerParam()->GetP() : esdtrack->P();
 404   } else {
 405     // AOD analysis: Use TPC momentum stored in the AliAODpid object
 406     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 407     if(aodtrack) p = aodtrack->GetDetPid() ? aodtrack->GetDetPid()->GetTPCmomentum() : aodtrack->P();
 408   }
 409   return p;
 410 }