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  *                                                                      *
  17  * Class for TPC PID                                                    *
  18  * Implements the abstract base class AliHFEpidBase                     *
  19  *                                                                      *
  20  * Class contains TPC specific cuts and QA histograms                   *
  21  * Two selection strategies are offered: Selection of certain value     *
  22  * regions in the TPC dE/dx (by IsSelected), and likelihoods            *
  23  *                                                                      *
  24  * Authors:                                                             *
  25  *                                                                      *
  26  *   Markus Fasel <M.Fasel@gsi.de>                                      *
  27  *   Markus Heide <mheide@uni-muenster.de>                              *
  28  *                                                                      *
  29  *                                                                      *
  30  ************************************************************************/
  31 #include <TH2I.h>
  32 #include <TList.h>
  33 #include <TMath.h>
  34 #include <TParticle.h>
  35
  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 "AliTPCpidESD.h"
  44 #include "AliVParticle.h"
  45
  46 #include "AliHFEpidTPC.h"
  47
  48
  49
  50 //___________________________________________________________________
  51 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  52   // add a list here
  53     AliHFEpidBase(name)
  54   , fLineCrossingsEnabled(0)
  55   , fNsigmaTPC(3)
  56   , fRejectionEnabled(0)
  57   , fPID(NULL)
  58   , fPIDtpcESD(NULL)
  59   , fQAList(NULL)
  60 {
  61   //
  62   // default  constructor
  63   //
  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   fPID = new AliPID;
  68   fPIDtpcESD = new AliTPCpidESD;
  69 }
  70
  71 //___________________________________________________________________
  72 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  73     AliHFEpidBase("")
  74   , fLineCrossingsEnabled(0)
  75   , fNsigmaTPC(2)
  76   , fRejectionEnabled(0)
  77   , fPID(NULL)
  78   , fPIDtpcESD(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.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD);
 111   target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
 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(fPIDtpcESD) delete fPIDtpcESD;
 126   if(fQAList){
 127     fQAList->Delete();
 128     delete fQAList;
 129   }
 130 }
 131
 132 //___________________________________________________________________
 133 Bool_t AliHFEpidTPC::InitializePID(){
 134   //
 135   // Add TPC dE/dx Line crossings
 136   //
 137   //AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
 138   //AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
 139   return kTRUE;
 140 }
 141
 142 //___________________________________________________________________
 143 Int_t AliHFEpidTPC::IsSelected(AliHFEpidObject *track)
 144 {
 145   //
 146   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
 147   // for electrons
 148   // exclusion of the crossing points
 149   //
 150   if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
 151     AliESDtrack *esdTrack = NULL;
 152     if(!(esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack))) return 0;
 153     AliMCParticle *mctrack = dynamic_cast<AliMCParticle *>(track->fMCtrack);
 154     return MakePIDesd(esdTrack, mctrack);
 155   }else{
 156     AliAODTrack *aodtrack = NULL;
 157     if(!(aodtrack = dynamic_cast<AliAODTrack *>(track->fRecTrack))) return 0;
 158     AliAODMCParticle *aodmctrack = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
 159     return MakePIDaod(aodtrack, aodmctrack);
 160   }
 161 }
 162
 163 //___________________________________________________________________
 164 Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
 165   //
 166   //  Doing TPC PID as explained in IsSelected for ESD tracks
 167   //
 168   if(IsQAon()) FillTPChistograms(esdTrack, mctrack);
 169   // exclude crossing points:
 170   // Determine the bethe values for each particle species
 171   Bool_t isLineCrossing = kFALSE;
 172   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 173     if(ispecies == AliPID::kElectron) continue;
 174     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 175     if(fPIDtpcESD->GetNumberOfSigmas(esdTrack, (AliPID::EParticleType)ispecies) < fLineCrossingSigma[ispecies]){
 176       // Point in a line crossing region, no PID possible
 177       isLineCrossing = kTRUE;
 178       break;
 179     }
 180   }
 181   if(isLineCrossing) return 0;
 182
 183   // Check particle rejection
 184   if(HasParticleRejection()){
 185     Int_t reject = Reject(esdTrack);
 186     AliDebug(1, Form("PID code from Rejection: %d", reject));
 187     if(reject != 0) return reject;
 188   }
 189   // Check whether distance from the electron line is smaller than n-sigma
 190
 191   // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
 192   Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron);
 193   Float_t p = 0.;
 194   Int_t pdg = 0;
 195   if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){
 196     if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11;
 197   } else {
 198     if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11;
 199   }
 200   if(IsQAon() && pdg != 0) (dynamic_cast<TH2I *>(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal());
 201   return pdg;
 202 }
 203
 204 //___________________________________________________________________
 205 Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mctrack*/){
 206   AliError("AOD PID not yet implemented");
 207   return 0;
 208 }
 209
 210 //___________________________________________________________________
 211 Int_t AliHFEpidTPC::Reject(AliESDtrack *track){
 212   //
 213   // reject particles based on asymmetric sigma cut
 214   //
 215   Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212};
 216   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 217   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 218     if(!TESTBIT(fRejectionEnabled, ispec)) continue;
 219     AliDebug(1, Form("Particle Rejection enabled for species %d", ispec));
 220     // Particle rejection enabled
 221     if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue;
 222     Double_t sigma = fPIDtpcESD->GetNumberOfSigmas(track, static_cast<AliPID::EParticleType>(ispec));
 223     AliDebug(1, Form("Sigma %f, min %f, max %f", sigma, fRejection[4*ispec + 1], fRejection[4*ispec+3]));
 224     if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge();
 225   }
 226   return 0;
 227 }
 228
 229 //___________________________________________________________________
 230 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){
 231   //
 232   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 233   // Stores line center and line sigma
 234   //
 235   if(species >= AliPID::kSPECIES){
 236     AliError("Species doesn't exist");
 237     return;
 238   }
 239   fLineCrossingsEnabled |= 1 << species;
 240   fLineCrossingSigma[species] = sigma;
 241 }
 242
 243 //___________________________________________________________________
 244 Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
 245 {
 246   //gives probability for track to come from a certain particle species;
 247   //no priors considered -> probabilities for equal abundances of all species!
 248   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
 249
 250   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
 251
 252   if(!track) return -1.;
 253   Bool_t outlier = kTRUE;
 254   // Check whether distance from the respective particle line is smaller than r sigma
 255   for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
 256     if(TMath::Abs(fPIDtpcESD->GetNumberOfSigmas(track, (AliPID::EParticleType)hypo)) > rsig)
 257       outlier = kTRUE;
 258     else {
 259             outlier = kFALSE;
 260             break;
 261           }
 262   }
 263   if(outlier)
 264     return -2.;
 265
 266   Double_t tpcProb[5];
 267
 268   track->GetTPCpid(tpcProb);
 269
 270   return tpcProb[species];
 271 }
 272
 273 //___________________________________________________________________
 274 Double_t  AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
 275 {
 276   //ratio of likelihoods to be whatever species/to be an electron;
 277   //as a cross-check for possible particle type suppression compared to electrons
 278   if(!track) return -20;
 279   if((Likelihood(track,species) == -2.)||(Likelihood(track,0)== -2.))
 280     return -30;
 281   if(Likelihood(track,species) == 0.)
 282     return -10;
 283   if (Likelihood(track,0) == 0.)
 284     return 10.;
 285   else
 286     return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
 287
 288
 289 }
 290
 291 //___________________________________________________________________
 292 void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack){
 293   //
 294   if(!track)
 295     return;
 296
 297   Double_t tpcSignal = track->GetTPCsignal();
 298   Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
 299   if(HasMCData()){
 300     switch(TMath::Abs(mctrack->Particle()->GetPdgCode())){
 301       case 11:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpcSignal);
 302                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
 303                         //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
 304                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
 305                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
 306                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
 307                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
 308                         //___________________________________________________________________________________________
 309                         //Likelihoods for electrons to be other particle species
 310                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
 311                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
 312                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
 313                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
 314                         break;
 315             //___________________________________________________________________________________________
 316       case 13:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpcSignal);
 317                   //Likelihood of muon to be an electron
 318                   (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
 319                   //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
 320                   //below functions are the same for other species
 321                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
 322                   break;
 323       case 211:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpcSignal);
 324                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
 325                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
 326                   break;
 327       case 321:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpcSignal);
 328                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
 329                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
 330                   break;
 331       case 2212:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpcSignal);
 332                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
 333                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
 334                   break;
 335       default:    (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpcSignal);
 336                         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
 337
 338                         break;
 339     }
 340   }
 341   //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
 342   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpcSignal);
 343   (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
 344 }
 345
 346 //___________________________________________________________________
 347 void AliHFEpidTPC::AddQAhistograms(TList *qaList){
 348   fQAList = new TList;
 349   fQAList->SetName("fTPCqaHistos");
 350
 351   fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 200, 0, 200), kHistTPCelectron);
 352   fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 200, 0, 200), kHistTPCmuon);
 353   fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 200, 0, 200), kHistTPCpion);
 354   fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 200, 0, 200), kHistTPCkaon);
 355   fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 200, 0, 200), kHistTPCproton);
 356   fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 200, 0, 200), kHistTPCothers);
 357   fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 200, 0, 200), kHistTPCall);
 358   fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 200, 0, 200), kHistTPCselected);
 359
 360   fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
 361   fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
 362   fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
 363   fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
 364   fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
 365   fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
 366   fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobAll);
 367
 368
 369  fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 370  fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 371  fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 372  fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
 373
 374  fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 375  fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 376  fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 377  fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
 378
 379  fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
 380  fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
 381  fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
 382  fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);
 383
 384
 385   qaList->AddLast(fQAList);
 386 }