PWG3/hfe/AliHFEpidTOF.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 TOF PID
  17 // Implements the abstract base class AliHFEpidBase
  18 // IsInitialized() does the PID decision
  19 //
  20 // Authors:
  21 //   Markus Fasel  <M.Fasel@gsi.de>
  22 //   Matus Kalisky <matus.kalisky@cern.ch>  (contact)
  23 //
  24
  25 #include <TH2F.h>
  26 #include <TList.h>
  27 #include <TMath.h>
  28
  29 #include "AliAODTrack.h"
  30 #include "AliAODMCParticle.h"
  31 #include "AliESDtrack.h"
  32 #include "AliLog.h"
  33 #include "AliMCParticle.h"
  34 #include "AliPID.h"
  35 #include "AliESDpid.h"
  36
  37 #include "AliHFEpidTOF.h"
  38 #include "AliHFEpidBase.h"
  39
  40
  41 ClassImp(AliHFEpidTOF)
  42
  43 //___________________________________________________________________
  44 AliHFEpidTOF::AliHFEpidTOF(const Char_t *name):
  45   AliHFEpidBase(name)
  46   , fPID(0x0)
  47   , fQAList(0x0)
  48   , fESDpid(NULL)
  49   , fNsigmaTOF(3)
  50 {
  51   //
  52   // Constructor
  53   //
  54
  55   fESDpid = new AliESDpid;
  56
  57 }
  58 //___________________________________________________________________
  59 AliHFEpidTOF::AliHFEpidTOF(const AliHFEpidTOF &c):
  60   AliHFEpidBase("")
  61   , fPID(0x0)
  62   , fQAList(0x0)
  63   , fESDpid(NULL)
  64   , fNsigmaTOF(3)
  65 {
  66   //
  67   // Copy operator
  68   //
  69
  70   c.Copy(*this);
  71 }
  72 //___________________________________________________________________
  73 AliHFEpidTOF &AliHFEpidTOF::operator=(const AliHFEpidTOF &ref){
  74   //
  75   // Assignment operator
  76   //
  77
  78   if(this != &ref){
  79     ref.Copy(*this);
  80   }
  81
  82   return *this;
  83 }
  84 //___________________________________________________________________
  85 AliHFEpidTOF::~AliHFEpidTOF(){
  86   //
  87   // Destructor
  88   //
  89   if(fPID) delete fPID;
  90   if(fESDpid) delete fESDpid;
  91   if(fQAList){
  92     fQAList->Delete();
  93     delete fQAList;
  94   }
  95 }
  96 //___________________________________________________________________
  97 void AliHFEpidTOF::Copy(TObject &ref) const {
  98   //
  99   // Performs the copying of the object
 100   //
 101   AliHFEpidTOF &target = dynamic_cast<AliHFEpidTOF &>(ref);
 102
 103   target.fPID = fPID;
 104   target.fQAList = fQAList;
 105   target.fESDpid = new AliESDpid(*fESDpid);
 106
 107   AliHFEpidBase::Copy(ref);
 108 }
 109 //___________________________________________________________________
 110 Bool_t AliHFEpidTOF::InitializePID(){
 111   //
 112   // InitializePID: TOF experts have to implement code here
 113   //
 114   return kTRUE;
 115 }
 116
 117 //___________________________________________________________________
 118 Int_t AliHFEpidTOF::IsSelected(AliHFEpidObject *vtrack)
 119 {
 120
 121   //
 122   // as of 22/05/2006 :
 123   // returns AliPID based on the ESD TOF PID decision
 124   // the ESD PID will be checked and if necessary improved
 125   // in the mean time. Best of luck
 126   //
 127   // returns 10 (== kUnknown)if PID can not be assigned for whatever reason
 128   //
 129
 130   if(vtrack->fAnalysisType == AliHFEpidObject::kESDanalysis){
 131     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(vtrack->fRecTrack);
 132     if(!esdTrack) return 0;
 133     AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(vtrack->fMCtrack);
 134     return MakePIDesd(esdTrack, mcTrack);
 135   } else {
 136     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(vtrack->fRecTrack);
 137     if(!aodTrack) return 0;
 138     AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(vtrack->fMCtrack);
 139     return MakePIDaod(aodTrack, aodmc);
 140   }
 141 }
 142
 143 //___________________________________________________________________
 144 Int_t AliHFEpidTOF::MakePIDesd(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
 145   //
 146   // Does particle identification as discribed in IsSelected
 147   //
 148   Long_t status = 0;
 149   status = track->GetStatus();
 150
 151   if(!(status & AliESDtrack::kTOFout)) return 0;
 152
 153   if(IsQAon())(dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(0.);
 154
 155   Double_t tItrackL = track->GetIntegratedLength();
 156   Double_t tTOFsignal = track->GetTOFsignal();
 157
 158   if(IsQAon()){
 159     if(tItrackL > 0)
 160       (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(1.);
 161
 162     if(tTOFsignal > 0)
 163       (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(2.);
 164   }
 165
 166
 167   if(tItrackL <=0 || tTOFsignal <=0) return 0;
 168
 169   if(IsQAon()){
 170     (dynamic_cast<TH1F *>(fQAList->At(kHistTOFpidFlags)))->Fill(3.);
 171     (dynamic_cast<TH1F *>(fQAList->At(kHistTOFsignal)))->Fill(tTOFsignal/1000.);
 172     (dynamic_cast<TH1F *>(fQAList->At(kHistTOFlength)))->Fill(tItrackL);
 173   }
 174   // get the TOF pid probabilities
 175   Double_t tESDpid[5] = {0., 0., 0., 0., 0.};
 176   Float_t tTOFpidSum = 0.;
 177   // find the largest PID probability
 178   track->GetTOFpid(tESDpid);
 179   Double_t tMAXpid = 0.;
 180   Int_t tMAXindex = -1;
 181   for(Int_t i=0; i<5; ++i){
 182     tTOFpidSum += tESDpid[i];
 183     if(tESDpid[i] > tMAXpid){
 184       tMAXpid = tESDpid[i];
 185       tMAXindex = i;
 186     }
 187   }
 188
 189   Int_t pdg = 0;
 190
 191   switch(tMAXindex){
 192     case 0:    pdg = 11; break;
 193     case 1:    pdg = 13; break;
 194     case 2:    pdg = 211; break;
 195     case 3:    pdg = 321; break;
 196     case 4:    pdg = 2212; break;
 197     default:   pdg = 0;
 198   };
 199
 200
 201   Double_t p = track->GetOuterParam()->P();
 202   Double_t beta = (tItrackL/100.)/(TMath::C()*(tTOFsignal/1e12));
 203
 204   // sanity check, should not be necessary
 205   if(TMath::Abs(tTOFpidSum - 1) > 0.01) return 0;
 206
 207   Double_t nSigmas = fESDpid->NumberOfSigmasTOF(track, (AliPID::EParticleType)tMAXindex, 0.);
 208   if(TMath::Abs(nSigmas) > fNsigmaTOF) return 0;
 209
 210
 211   // should be the same as AliPID flags
 212
 213   if(IsQAon()){
 214     (dynamic_cast<TH2F *>(fQAList->At(kHistTOFpid0+tMAXindex)))->Fill(beta, p);
 215     (dynamic_cast<TH2F *>(fQAList->At(kHistTOFpidBetavP)))->Fill(beta, p);
 216   }
 217   //return tMAXindex;
 218   return pdg;
 219
 220 }
 221 //___________________________________________________________________
 222 Double_t AliHFEpidTOF::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig){
 223
 224   //gives probability for track to come from a certain particle species;
 225   //no priors considered -> probabilities for equal abundances of all species!
 226   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
 227
 228   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
 229
 230   if(!track) return -1.;
 231   Bool_t outlier = kTRUE;
 232   // Check whether distance from the respective particle line is smaller than r sigma
 233   for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
 234     if(TMath::Abs(fESDpid->NumberOfSigmasTOF(track, (AliPID::EParticleType)hypo, 0.)) > rsig)
 235       outlier = kTRUE;
 236     else {
 237       outlier = kFALSE;
 238       break;
 239     }
 240   }
 241   if(outlier)
 242     return -2.;
 243
 244   Double_t tofProb[5];
 245
 246   track->GetTOFpid(tofProb);
 247
 248   return tofProb[species];
 249 }
 250 //___________________________________________________________________
 251 Int_t AliHFEpidTOF::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mcTrack*/){
 252   AliError("AOD PID not yet implemented");
 253   return 0;
 254 }
 255
 256 //___________________________________________________________________
 257 void AliHFEpidTOF::AddQAhistograms(TList *qaList){
 258   //
 259   // Create QA histograms for TOF PID
 260   //
 261
 262   fQAList = new TList;
 263   fQAList->SetName("fTOFqaHistos");
 264   fQAList->AddAt(new TH1F("hTOF_flags", "TOF flags;flags (see code for info);counts", 10, -0.25, 4.75), kHistTOFpidFlags);
 265   fQAList->AddAt(new TH2F("fTOFbeta_v_P_no","beta -v- P; beta;momentum [GeV/c]", 120, 0, 1.2, 200, 0, 20), kHistTOFpidBetavP);
 266   fQAList->AddAt(new TH1F("hTOF_signal", "TOF signal; TOF signal [ns];counts", 1000, 12, 50), kHistTOFsignal);
 267   fQAList->AddAt(new TH1F("hTOF_length", "TOF track length; length [cm];counts", 400, 300, 700), kHistTOFlength);
 268   fQAList->AddAt(new TH2F("hTOFpid_electron", "TOF reco electron; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5), kHistTOFpid0);
 269   fQAList->AddAt(new TH2F("hTOFpid_muon", "TOF reco muon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5), kHistTOFpid1);
 270   fQAList->AddAt(new TH2F("hTOFpid_pion", "TOF reco pion; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5), kHistTOFpid2);
 271   fQAList->AddAt(new TH2F("hTOFpid_kaon", "TOF reco kaon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5), kHistTOFpid3);
 272   fQAList->AddAt(new TH2F("hTOFpid_proton", "TOF reco proton; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5), kHistTOFpid4);
 273
 274   qaList->AddLast(fQAList);
 275 }
 276
 277