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