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 <TH2I.h>
  29 #include <TList.h>
  30 #include <TMath.h>
  31 //#include <TParticle.h>
  32
  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 //#include "AliVParticle.h"
  42
  43 #include "AliHFEcollection.h"
  44 #include "AliHFEpidTPC.h"
  45
  46
  47
  48 //___________________________________________________________________
  49 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  50   // add a list here
  51   AliHFEpidBase(name)
  52   , fLineCrossingType(0)
  53   , fLineCrossingsEnabled(0)
  54   , fNsigmaTPC(3)
  55   , fRejectionEnabled(0)
  56   , fPID(NULL)
  57   , fESDpid(NULL)
  58   , fQAList(NULL)
  59 {
  60   //
  61   // default  constructor
  62   //
  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   fPID = new AliPID;
  67   fESDpid = new AliESDpid;
  68 }
  69
  70 //___________________________________________________________________
  71 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  72   AliHFEpidBase("")
  73   , fLineCrossingType(0)
  74   , fLineCrossingsEnabled(0)
  75   , fNsigmaTPC(2)
  76   , fRejectionEnabled(0)
  77   , fPID(NULL)
  78   , fESDpid(NULL)
  79   , fQAList(NULL)
  80 {
  81   //
  82   // Copy constructor
  83   //
  84   ref.Copy(*this);
  85 }
  86
  87 //___________________________________________________________________
  88 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
  89   //
  90   // Assignment operator
  91   //
  92   if(this != &ref){
  93     ref.Copy(*this);
  94   }
  95   return *this;
  96 }
  97
  98 //___________________________________________________________________
  99 void AliHFEpidTPC::Copy(TObject &o) const{
 100   //
 101   // Copy function
 102   // called in copy constructor and assigment operator
 103   //
 104   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 105
 106   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
 107   target.fNsigmaTPC = fNsigmaTPC;
 108   target.fRejectionEnabled = fRejectionEnabled;
 109   target.fPID = new AliPID(*fPID);
 110   target.fESDpid = new AliESDpid(*fESDpid);
 111   target.fQAList = new AliHFEcollection(*fQAList);
 112   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
 113   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
 114   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
 115
 116   AliHFEpidBase::Copy(target);
 117 }
 118
 119 //___________________________________________________________________
 120 AliHFEpidTPC::~AliHFEpidTPC(){
 121   //
 122   // Destructor
 123   //
 124   if(fPID) delete fPID;
 125   if(fESDpid) delete fESDpid;
 126   if(fQAList){
 127     delete fQAList;
 128   }
 129 }
 130
 131 //___________________________________________________________________
 132 Bool_t AliHFEpidTPC::InitializePID(){
 133   //
 134   // Add TPC dE/dx Line crossings
 135   //
 136   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
 137   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
 138   return kTRUE;
 139 }
 140
 141 //___________________________________________________________________
 142 Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
 143 {
 144   //
 145   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
 146   // for electrons
 147   // exclusion of the crossing points
 148   //
 149   if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
 150     AliESDtrack *esdTrack = NULL;
 151     if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
 152     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track->fMCtrack);
 153     return MakePIDesd(esdTrack, mctrack);
 154   }else{
 155     AliAODTrack *aodtrack = NULL;
 156     if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
 157     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
 158     return MakePIDaod(aodtrack, aodmctrack);
 159   }
 160 }
 161 //___________________________________________________________________
 162 Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
 163   //
 164   //  Doing TPC PID as explained in IsSelected for ESD tracks
 165   //
 166   Float_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron);
 167   if(IsQAon()){
 168     FillTPChistograms(esdTrack, mctrack);
 169     fQAList->Fill("fHistSigmaElectronAll",  esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), nsigma);
 170   }
 171   // exclude crossing points:
 172   // Determine the bethe values for each particle species
 173   Bool_t isLineCrossing = kFALSE;
 174   fLineCrossingType = 0;  // default value
 175   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 176     if(ispecies == AliPID::kElectron) continue;
 177     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 178     if(TMath::Abs(fESDpid->NumberOfSigmasTPC(esdTrack, (AliPID::EParticleType)ispecies)) < fLineCrossingSigma[ispecies] && TMath::Abs(nsigma) < fNsigmaTPC){
 179       // Point in a line crossing region, no PID possible, but !PID still possible ;-)
 180       isLineCrossing = kTRUE;
 181       fLineCrossingType = ispecies;
 182       break;
 183     }
 184   }
 185   if(isLineCrossing) return 0;
 186
 187   // Check particle rejection
 188   if(HasParticleRejection()){
 189     Int_t reject = Reject(esdTrack);
 190     if(reject != 0) return reject;
 191   }
 192   // Check whether distance from the electron line is smaller than n-sigma
 193
 194   // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
 195   Float_t p = 0.;
 196   Int_t pdg = 0;
 197   if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
 198     if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
 199   } else {
 200     if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
 201   }
 202   if(IsQAon() && pdg != 0){
 203     fQAList->Fill("fHistTPCselected", esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
 204     fQAList->Fill("fHistSigmaElectronSelected",  esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), nsigma);
 205   }
 206
 207   return pdg;
 208 }
 209
 210 //___________________________________________________________________
 211 Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mctrack*/){
 212   AliError("AOD PID not yet implemented");
 213   return 0;
 214 }
 215
 216 //___________________________________________________________________
 217 Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
 218   //
 219   // reject particles based on asymmetric sigma cut
 220   //
 221   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
 222   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 223   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 224     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
 225     // Particle rejection enabled
 226     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
 227     Double_t sigma = fESDpid->NumberOfSigmasTPC(track, static_cast<AliPID::EParticleType>(ispec));
 228     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
 229   }
 230   return 0;
 231 }
 232
 233 //___________________________________________________________________
 234 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
 235   //
 236   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 237   // Stores line center and line sigma
 238   //
 239   if(species >= AliPID::kSPECIES){
 240     AliError("Species doesn't exist");
 241     return;
 242   }
 243   fLineCrossingsEnabled |= 1 << species;
 244   fLineCrossingSigma[species] = sigma;
 245 }
 246
 247 //___________________________________________________________________
 248 Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
 249 {
 250   //gives probability for track to come from a certain particle species;
 251   //no priors considered -> probabilities for equal abundances of all species!
 252   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
 253
 254   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
 255
 256   if(!track) return -1.;
 257   Bool_t outlier = kTRUE;
 258   // Check whether distance from the respective particle line is smaller than r sigma
 259   for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
 260     if(TMath::Abs(fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)hypo)) > rsig)
 261       outlier = kTRUE;
 262     else {
 263             outlier = kFALSE;
 264             break;
 265           }
 266   }
 267   if(outlier)
 268     return -2.;
 269
 270   Double_t tpcProb[5];
 271
 272   track->GetTPCpid(tpcProb);
 273
 274   return tpcProb[species];
 275 }
 276
 277 //___________________________________________________________________
 278 Double_t  AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
 279 {
 280   //ratio of likelihoods to be whatever species/to be an electron;
 281   //as a cross-check for possible particle type suppression compared to electrons
 282   const Double_t kVerySmall = 1e-10;
 283   if(!track) return -20;
 284   if((TMath::Abs(Likelihood(track,species) + 2) < kVerySmall)||(TMath::Abs(Likelihood(track,0) + 2 ) < kVerySmall))
 285     return -30;
 286   if(TMath::Abs(Likelihood(track,species)) < kVerySmall)
 287     return -10;
 288   if (TMath::Abs(Likelihood(track,0)) < kVerySmall)
 289     return 10.;
 290   else
 291     return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
 292
 293
 294 }
 295
 296 //___________________________________________________________________
 297 void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack){
 298   //
 299   // Fill the QA histogtrams
 300   //
 301   if(!track)
 302     return;
 303
 304   Double_t tpcSignal = track->GetTPCsignal();
 305   Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
 306   if(HasMCData()){
 307     switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
 308       case 11:    fQAList->Fill("fHistTPCelectron", p, tpcSignal);
 309                         fQAList->Fill("fHistTPCprobEl", p, Likelihood(track, 0));
 310                         //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
 311                         fQAList->Fill("fHistTPCenhanceElPi", p, -Suppression(track, 2));
 312                         fQAList->Fill("fHistTPCenhanceElMu", p, -Suppression(track, 1));
 313                         fQAList->Fill("fHistTPCenhanceElKa", p, -Suppression(track, 3));
 314                         fQAList->Fill("fHistTPCenhanceElPro", p, -Suppression(track, 4));
 315                         //___________________________________________________________________________________________
 316                         //Likelihoods for electrons to be other particle species
 317                         fQAList->Fill("fHistTPCElprobPi", p, Likelihood(track, 2));
 318                         fQAList->Fill("fHistTPCElprobMu", p, Likelihood(track, 1));
 319                         fQAList->Fill("fHistTPCElprobKa", p, Likelihood(track, 3));
 320                         fQAList->Fill("fHistTPCElprobPro", p, Likelihood(track, 4));
 321                         break;
 322             //___________________________________________________________________________________________
 323       case 13:    fQAList->Fill("fHistTPCmuon", p, tpcSignal);
 324                   //Likelihood of muon to be an electron
 325                   fQAList->Fill("fHistTPCprobMu", p, Likelihood(track, 0));
 326                   //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
 327                   //below functions are the same for other species
 328                         fQAList->Fill("fHistTPCsuppressMu", p, Suppression(track, 1));
 329                   break;
 330       case 211:   fQAList->Fill("fHistTPCpion", p, tpcSignal);
 331                         fQAList->Fill("fHistTPCprobPi", p, Likelihood(track, 0));
 332                         fQAList->Fill("fHistTPCsuppressPi", p, Suppression(track, 2));
 333                   break;
 334       case 321:   fQAList->Fill("fHistTPCkaon", p, tpcSignal);
 335                         fQAList->Fill("fHistTPCprobKa", p, Likelihood(track, 0));
 336                         fQAList->Fill("fHistTPCsuppressKa", p, Suppression(track, 3));
 337                   break;
 338       case 2212:  fQAList->Fill("fHistTPCproton", p, tpcSignal);
 339                         fQAList->Fill("fHistTPCprobPro", p, Likelihood(track, 0));
 340                         fQAList->Fill("fHistTPCsuppressPro", p, Suppression(track, 4));
 341                   break;
 342       default:    fQAList->Fill("fHistTPCothers", p, tpcSignal);
 343                         fQAList->Fill("fHistTPCprobOth", p, Likelihood(track, 0));
 344
 345                         break;
 346     }
 347   }
 348   //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
 349   fQAList->Fill("fHistTPCall", p, tpcSignal);
 350   fQAList->Fill("kHistTPCprobAll", p, Likelihood(track, 0));
 351 }
 352
 353 //___________________________________________________________________
 354 void AliHFEpidTPC::AddQAhistograms(TList *qaList){
 355   //
 356   // Create QA histograms for TPC PID
 357   //
 358   fQAList = new AliHFEcollection;
 359
 360   fQAList->CreateTH2F("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600);
 361   fQAList->CreateTH2F("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600);
 362   fQAList->CreateTH2F("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600);
 363   fQAList->CreateTH2F("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600);
 364   fQAList->CreateTH2F("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600);
 365   fQAList->CreateTH2F("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600);
 366   fQAList->CreateTH2F("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600);
 367   fQAList->CreateTH2F("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 60, 0, 600);
 368
 369   fQAList->CreateTH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.);
 370   fQAList->CreateTH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 371   fQAList->CreateTH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 372   fQAList->CreateTH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 373   fQAList->CreateTH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 374   fQAList->CreateTH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 375   fQAList->CreateTH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.);
 376
 377   fQAList->CreateTH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8);
 378   fQAList->CreateTH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8);
 379   fQAList->CreateTH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8);
 380   fQAList->CreateTH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8);
 381
 382   fQAList->CreateTH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8);
 383   fQAList->CreateTH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8);
 384   fQAList->CreateTH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8);
 385   fQAList->CreateTH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8);
 386
 387   fQAList->CreateTH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.);
 388   fQAList->CreateTH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.);
 389   fQAList->CreateTH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.);
 390   fQAList->CreateTH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.);
 391
 392   fQAList->CreateTH2F("fHistSigmaElectronAll", "TPC NSigma around the Electron Line", 200, 0, 20, 40, -10, 10);
 393   fQAList->CreateTH2F("fHistSigmaElectronSelected", "TPC NSigma around the Electron Line for selected Tracks", 200, 0, 20, 40, -10, 10);
 394
 395   qaList->AddLast(fQAList->GetList());
 396 }
 397
 398 //___________________________________________________________________
 399 void AliHFEpidTPC::SetBetheBlochParameters(Double_t *pars){
 400   //
 401   // Set non-default Bethe-Bloch Parameters
 402   //
 403   fESDpid->GetTPCResponse().SetBetheBlochParameters(pars[0], pars[1], pars[2], pars[3], pars[4]);
 404 }