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
  35 #include "AliESDtrack.h"
  36 #include "AliExternalTrackParam.h"
  37 #include "AliLog.h"
  38 #include "AliPID.h"
  39 #include "AliTPCpidESD.h"
  40 #include "AliVParticle.h"
  41
  42 #include "AliHFEpidTPC.h"
  43
  44
  45
  46 //___________________________________________________________________
  47 AliHFEpidTPC::AliHFEpidTPC(const char* name) :
  48   // add a list here
  49     AliHFEpidBase(name)
  50   , fLineCrossingsEnabled(0)
  51   , fNsigmaTPC(2)
  52   , fPID(0x0)
  53   , fPIDtpcESD(0x0)
  54   , fQAList(0x0)
  55 {
  56   //
  57   // default  constructor
  58   //
  59   memset(fLineCrossingCenter, 0, sizeof(Double_t) * AliPID::kSPECIES);
  60   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
  61   fPID = new AliPID;
  62   fPIDtpcESD = new AliTPCpidESD;
  63 }
  64
  65 //___________________________________________________________________
  66 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
  67     AliHFEpidBase("")
  68   , fLineCrossingsEnabled(0)
  69   , fNsigmaTPC(2)
  70   , fPID(0x0)
  71   , fPIDtpcESD(0x0)
  72   , fQAList(0x0)
  73 {
  74   //
  75   // Copy constructor
  76   //
  77   ref.Copy(*this);
  78 }
  79
  80 //___________________________________________________________________
  81 AliHFEpidTPC &AliHFEpidTPC::operator=(const AliHFEpidTPC &ref){
  82   //
  83   // Assignment operator
  84   //
  85   if(this != &ref){
  86     ref.Copy(*this);
  87   }
  88   return *this;
  89 }
  90
  91 void AliHFEpidTPC::Copy(TObject &o) const{
  92   //
  93   // Copy function
  94   // called in copy constructor and assigment operator
  95   //
  96   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
  97
  98   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
  99   target.fNsigmaTPC = fNsigmaTPC;
 100   target.fPID = new AliPID(*fPID);
 101   target.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD);
 102   target.fQAList = dynamic_cast<TList *>(fQAList->Clone());
 103
 104   AliHFEpidBase::Copy(target);
 105 }
 106
 107 //___________________________________________________________________
 108 AliHFEpidTPC::~AliHFEpidTPC(){
 109   //
 110   // Destructor
 111   //
 112   if(fPID) delete fPID;
 113   if(fPIDtpcESD) delete fPIDtpcESD;
 114   if(fQAList){
 115     fQAList->Delete();
 116     delete fQAList;
 117   }
 118 }
 119
 120 //___________________________________________________________________
 121 Bool_t AliHFEpidTPC::InitializePID(){
 122   //
 123   // Add TPC dE/dx Line crossings
 124   //
 125   AddTPCdEdxLineCrossing(AliPID::kKaon, 0.3, 0.018);
 126   AddTPCdEdxLineCrossing(AliPID::kProton, 0.9, 0.054);
 127   return kTRUE;
 128 }
 129
 130 //___________________________________________________________________
 131 Int_t AliHFEpidTPC::IsSelected(AliVParticle *track)
 132 {
 133   //
 134   // For the TPC pid we use the 2-sigma band around the bethe bloch curve
 135   // for electrons
 136   // exclusion of the crossing points
 137   //
 138   AliESDtrack *esdTrack = 0x0;
 139   if(!(esdTrack = dynamic_cast<AliESDtrack *>(track))) return kFALSE;
 140   if(IsQAon()) FillTPChistograms(esdTrack);
 141   Double_t TPCsignal = esdTrack->GetTPCsignal();
 142   // exclude crossing points:
 143   // Determine the bethe values for each particle species
 144   Double_t p = esdTrack->GetInnerParam()->P();
 145   Bool_t isLineCrossing = kFALSE;
 146   for(Int_t ispecies = 0; ispecies < AliPID::kSPECIES; ispecies++){
 147     if(!(fLineCrossingsEnabled & 1 << ispecies)) continue;
 148     if(TMath::Abs(p - fLineCrossingCenter[ispecies]) < fLineCrossingSigma[ispecies]){
 149       // Point in a line crossing region, no PID possible
 150       isLineCrossing = kTRUE;
 151       break;
 152     }
 153   }
 154   if(isLineCrossing) return 0;
 155   // Check whether distance from the electron line is smaller than n-sigma
 156   Double_t beta = p/fPID->ParticleMass(AliPID::kElectron);
 157   if(TMath::Abs(TPCsignal - 50*fPIDtpcESD->Bethe(beta)) < GetTPCsigma(p,0)) return 11;
 158   return 0;
 159 }
 160
 161 //___________________________________________________________________
 162 void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t p, Double_t sigma_p){
 163   //
 164   // Add exclusion point for the TPC PID where a dEdx line crosses the electron line
 165   // Stores line center and line sigma
 166   //
 167   if(species >= AliPID::kSPECIES){
 168     AliError("Species doesn't exist");
 169     return;
 170   }
 171   fLineCrossingsEnabled |= 1 << species;
 172   fLineCrossingCenter[species] = p;
 173   fLineCrossingSigma[species] = sigma_p;
 174 }
 175
 176 //___________________________________________________________________
 177 Double_t AliHFEpidTPC::GetTPCsigma(Double_t p, Int_t species){
 178   //
 179   // return the TPC sigma, momentum dependent
 180   //
 181   if(p < 0.1 || p > 20.) return 0.;
 182   Double_t beta = p/fPID->ParticleMass(species);
 183
 184
 185   return 50*fPIDtpcESD->Bethe(beta) * 0.06;
 186 }
 187
 188 //___________________________________________________________________
 189 Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig)
 190 {
 191   //gives probability for track to come from a certain particle species;
 192   //no priors considered -> probabilities for equal abundances of all species!
 193   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
 194
 195   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
 196
 197   if(!track) return -1.;
 198   Int_t hypo; //marks particle hypotheses for 2-sigma bands
 199   Double_t beta;//well o.k., it corresponds to gamma * beta
 200   Double_t p = track->GetInnerParam()->P();
 201   Double_t TPCsignal = track->GetTPCsignal();
 202   Bool_t outlier = kTRUE;
 203   // Check whether distance from the respective particle line is smaller than r sigma
 204   for(hypo = 0; hypo < 5; hypo++)
 205     {
 206       beta = p/fPID->ParticleMass(hypo);
 207       if(TMath::Abs(TPCsignal - (GetTPCsigma(p, hypo))/0.06) > (rsig * GetTPCsigma(p,hypo)))
 208         outlier = kTRUE;
 209       else
 210         {
 211           outlier = kFALSE;
 212           break;
 213         }
 214     }
 215   if(outlier)
 216     return -2.;
 217
 218   Double_t TPCprob[5];
 219
 220   track->GetTPCpid(TPCprob);
 221
 222   return TPCprob[species];
 223 }
 224 //___________________________________________________________________
 225 Double_t AliHFEpidTPC::Likelihood(const AliESDtrack *track, Int_t species)
 226 {
 227   //default: rsig = 2.
 228   // for everything else, see above!
 229
 230   if(!track) return -1.;
 231   Int_t hypo; //marks particle hypotheses for 2-sigma bands
 232   Double_t beta;
 233   Double_t p = track->GetInnerParam()->P();
 234   Double_t TPCsignal = track->GetTPCsignal();
 235   Bool_t outlier = kTRUE;
 236   // Check whether distance from the respective particle line is smaller than 2 sigma
 237   for(hypo = 0; hypo < 5; hypo++)
 238     {
 239       beta = p/fPID->ParticleMass(hypo);
 240
 241       if(TMath::Abs(TPCsignal - (GetTPCsigma(p, hypo))/0.06) > (2. * GetTPCsigma(p,hypo)))
 242         outlier = kTRUE;
 243       else
 244         {
 245           outlier = kFALSE;
 246           break;
 247         }
 248     }
 249   if(outlier == kTRUE)
 250     return -2.;
 251
 252   Double_t TPCprob[5];
 253
 254   track->GetTPCpid(TPCprob);
 255
 256   return TPCprob[species];
 257 }
 258 //___________________________________________________________________
 259 Double_t  AliHFEpidTPC::Suppression(const AliESDtrack *track, Int_t species)
 260 {
 261   //ratio of likelihoods to be whatever species/to be an electron;
 262   //as a cross-check for possible particle type suppression compared to electrons
 263   if(!track) return -20;
 264   if((Likelihood(track,species) == -2.)||(Likelihood(track,0)== -2.))
 265     return -30;
 266   if(Likelihood(track,species) == 0.)
 267     return -10;
 268   if (Likelihood(track,0) == 0.)
 269     return 10.;
 270   else
 271     return TMath::Log10(Likelihood(track,species)/(Likelihood(track,0)));
 272
 273
 274 }
 275 //___________________________________________________________________
 276 void AliHFEpidTPC::FillTPChistograms(const AliESDtrack *track){
 277   //
 278   if(!track)
 279     return;
 280
 281   Double_t tpc_signal = track->GetTPCsignal();
 282   Double_t p = track->GetInnerParam() ? track->GetInnerParam()->P() : track->P();
 283   if(HasMCData()){
 284     switch(TMath::Abs(GetPdgCode(const_cast<AliESDtrack *>(track)))){
 285       case 11:   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCelectron)))->Fill(p, tpc_signal);
 286         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobEl)))->Fill(p, Likelihood(track, 0));
 287         //histograms with ratio of likelihood to be electron/to be other species (a check for quality of likelihood PID);
 288         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPi)))->Fill(p, -Suppression(track, 2));
 289         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElMu)))->Fill(p, -Suppression(track, 1));
 290         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElKa)))->Fill(p, -Suppression(track, 3));
 291         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCenhanceElPro)))->Fill(p, -Suppression(track, 4));
 292         //___________________________________________________________________________________________
 293         //Likelihoods for electrons to be other particle species
 294         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPi)))->Fill(p, Likelihood(track, 2));
 295         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobMu)))->Fill(p, Likelihood(track, 1));
 296         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobKa)))->Fill(p, Likelihood(track, 3));
 297         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCElprobPro)))->Fill(p, Likelihood(track, 4));
 298         break;
 299         //___________________________________________________________________________________________
 300     case 13: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCmuon)))->Fill(p, tpc_signal);
 301       //Likelihood of muon to be an electron
 302       (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobMu)))->Fill(p, Likelihood(track, 0));
 303       //ratio of likelihood for muon to be a muon/an electron -> indicator for quality of muon suppression
 304       //below functions are the same for other species
 305         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressMu)))->Fill(p, Suppression(track, 1));
 306       break;
 307       case 211:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCpion)))->Fill(p, tpc_signal);
 308         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPi)))->Fill(p, Likelihood(track, 0));
 309         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPi)))->Fill(p, Suppression(track, 2));
 310           break;
 311       case 321:  (dynamic_cast<TH2I *>(fQAList->At(kHistTPCkaon)))->Fill(p, tpc_signal);
 312         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobKa)))->Fill(p, Likelihood(track, 0));
 313         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressKa)))->Fill(p, Suppression(track, 3));
 314           break;
 315       case 2212: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCproton)))->Fill(p, tpc_signal);
 316         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobPro)))->Fill(p, Likelihood(track, 0));
 317         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCsuppressPro)))->Fill(p, Suppression(track, 4));
 318           break;
 319       default: (dynamic_cast<TH2I *>(fQAList->At(kHistTPCothers)))->Fill(p, tpc_signal);
 320         (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobOth)))->Fill(p, Likelihood(track, 0));
 321
 322         break;
 323     }
 324   }
 325   //TPC signal and Likelihood to be electron for all tracks (independent of MC information)
 326   (dynamic_cast<TH2I *>(fQAList->At(kHistTPCall)))->Fill(p, tpc_signal);
 327  (dynamic_cast<TH2F *>(fQAList->At(kHistTPCprobAll)))->Fill(p, Likelihood(track, 0));
 328
 329
 330
 331 }
 332
 333 //___________________________________________________________________
 334 void AliHFEpidTPC::AddQAhistograms(TList *qaList){
 335   fQAList = new TList;
 336   fQAList->SetName("fTPCqaHistos");
 337
 338   fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron);
 339   fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon);
 340   fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion);
 341   fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon);
 342   fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton);
 343   fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers);
 344   fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall);
 345
 346   fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl);
 347   fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPi);
 348   fQAList->AddAt(new TH2F("fHistTPCprobMu","TPC likelihood for muons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobMu);
 349   fQAList->AddAt(new TH2F("fHistTPCprobKa","TPC likelihood for kaons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobKa);
 350   fQAList->AddAt(new TH2F("fHistTPCprobPro","TPC likelihood for protons to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobPro);
 351   fQAList->AddAt(new TH2F("fHistTPCprobOth","TPC likelihood for other particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobOth);
 352   fQAList->AddAt(new TH2F("fHistTPCprobAll","TPC likelihood for all particles to be an electron vs. p",  200, 0.,20.,200, 0.,1.), kHistTPCprobAll);
 353
 354
 355  fQAList->AddAt(new TH2F("fHistTPCsuppressPi","log10 of TPC Likelihood(pion)/Likelihood(elec) for pions vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 356  fQAList->AddAt(new TH2F("fHistTPCsuppressMu","log10 of TPC Likelihood(muon)/Likelihood(elec) for muons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 357  fQAList->AddAt(new TH2F("fHistTPCsuppressKa","log10 of TPC Likelihood(kaon)/Likelihood(elec) for kaons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 358  fQAList->AddAt(new TH2F("fHistTPCsuppressPro","log10 of TPC Likelihood(proton)/Likelihood(elec)for protons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
 359
 360  fQAList->AddAt(new TH2F("fHistTPCenhanceElPi","log10 of TPC Likelihood(elec)/Likelihood(pion) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPi);
 361  fQAList->AddAt(new TH2F("fHistTPCenhanceElMu","log10 of TPC Likelihood(elec)/Likelihood(muon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressMu);
 362  fQAList->AddAt(new TH2F("fHistTPCenhanceElKa","log10 of TPC Likelihood(elec)/Likelihood(kaon) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressKa);
 363  fQAList->AddAt(new TH2F("fHistTPCenhanceElPro","log10 of TPC Likelihood(elec)/Likelihood(proton) for electrons vs. p", 200, 0.,20.,200,-1.,5.8), kHistTPCsuppressPro);
 364
 365  fQAList->AddAt(new TH2F("fHistTPCElprobPi","TPC likelihood for electrons to be a pion vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPi);
 366  fQAList->AddAt(new TH2F("fHistTPCElprobMu","TPC likelihood for electrons to be a muon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobMu);
 367  fQAList->AddAt(new TH2F("fHistTPCElprobKa","TPC likelihood for electrons to be a kaon vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobKa);
 368  fQAList->AddAt(new TH2F("fHistTPCElprobPro","TPC likelihood for electrons to be a proton vs. p", 200, 0.,20.,200,0.,1.), kHistTPCElprobPro);
 369
 370
 371   qaList->AddLast(fQAList);
 372 }