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 <TList.h>
  26 #include <TMath.h>
  27 #include <THnSparse.h>
  28 #include <TDatabasePDG.h>
  29
  30 #include "AliAODTrack.h"
  31 #include "AliAODMCParticle.h"
  32 #include "AliESDtrack.h"
  33 #include "AliMCParticle.h"
  34 #include "AliPID.h"
  35 #include "AliESDpid.h"
  36
  37 #include "AliHFEcollection.h"
  38 #include "AliHFEpidTOF.h"
  39 #include "AliHFEpidBase.h"
  40
  41
  42 ClassImp(AliHFEpidTOF)
  43
  44 //___________________________________________________________________
  45 AliHFEpidTOF::AliHFEpidTOF(const Char_t *name):
  46   AliHFEpidBase(name)
  47   , fPID(0x0)
  48   , fQAList(0x0)
  49   , fNsigmaTOF(3)
  50 {
  51   //
  52   // Constructor
  53   //
  54 }
  55 //___________________________________________________________________
  56 AliHFEpidTOF::AliHFEpidTOF(const AliHFEpidTOF &c):
  57   AliHFEpidBase("")
  58   , fPID(0x0)
  59   , fQAList(0x0)
  60   , fNsigmaTOF(3)
  61 {
  62   //
  63   // Copy operator
  64   //
  65
  66   c.Copy(*this);
  67 }
  68 //___________________________________________________________________
  69 AliHFEpidTOF &AliHFEpidTOF::operator=(const AliHFEpidTOF &ref){
  70   //
  71   // Assignment operator
  72   //
  73
  74   if(this != &ref){
  75     ref.Copy(*this);
  76   }
  77
  78   return *this;
  79 }
  80 //___________________________________________________________________
  81 AliHFEpidTOF::~AliHFEpidTOF(){
  82   //
  83   // Destructor
  84   //
  85   if(fPID) delete fPID;
  86   if(fQAList){
  87     fQAList->Delete();
  88     delete fQAList;
  89   }
  90 }
  91 //___________________________________________________________________
  92 void AliHFEpidTOF::Copy(TObject &ref) const {
  93   //
  94   // Performs the copying of the object
  95   //
  96   AliHFEpidTOF &target = dynamic_cast<AliHFEpidTOF &>(ref);
  97
  98   target.fPID = fPID;
  99   target.fQAList = fQAList;
 100   target.fNsigmaTOF = fNsigmaTOF;
 101
 102   AliHFEpidBase::Copy(ref);
 103 }
 104 //___________________________________________________________________
 105 Bool_t AliHFEpidTOF::InitializePID(){
 106   //
 107   // InitializePID: TOF experts have to implement code here
 108   //
 109   return kTRUE;
 110 }
 111
 112 //___________________________________________________________________
 113 Int_t AliHFEpidTOF::IsSelected(AliHFEpidObject *vtrack)
 114 {
 115
 116   //
 117   // as of 22/05/2006 :
 118   // returns AliPID based on the ESD TOF PID decision
 119   // the ESD PID will be checked and if necessary improved
 120   // in the mean time. Best of luck
 121   //
 122   // returns 10 (== kUnknown)if PID can not be assigned for whatever reason
 123   //
 124
 125   if(vtrack->fAnalysisType == AliHFEpidObject::kESDanalysis){
 126     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(vtrack->fRecTrack);
 127     if(!esdTrack) return 0;
 128     AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(vtrack->fMCtrack);
 129     return MakePIDesdV3(esdTrack, mcTrack);
 130   } else {
 131     AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(vtrack->fRecTrack);
 132     if(!aodTrack) return 0;
 133     AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(vtrack->fMCtrack);
 134     return MakePIDaod(aodTrack, aodmc);
 135   }
 136 }
 137
 138 //___________________________________________________________________
 139 Int_t AliHFEpidTOF::MakePIDesd(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
 140   //
 141   // Does particle identification as discribed in IsSelected
 142   //
 143   if(!fESDpid){
 144     AliError("No ESD PID object available");
 145     return kFALSE;
 146   }
 147   Long_t status = 0;
 148   status = track->GetStatus();
 149
 150   if(!(status & AliESDtrack::kTOFout)) return 0;
 151
 152   if(IsQAon()) fQAList->Fill("hTOF_flags", 0.);
 153
 154   Double_t tItrackL = track->GetIntegratedLength();
 155   Double_t tTOFsignal = track->GetTOFsignal();
 156
 157   if(IsQAon()){
 158     if(tItrackL > 0)
 159       fQAList->Fill("hTOF_flags", 1.);
 160
 161     if(tTOFsignal > 0)
 162       fQAList->Fill("hTOF_flags", 2.);
 163   }
 164
 165
 166   if(tItrackL <=0 || tTOFsignal <=0) return 0;
 167
 168   if(IsQAon()){
 169     fQAList->Fill("hTOF_flags", 3.);
 170     fQAList->Fill("hTOF_signal", tTOFsignal/1000.);
 171     fQAList->Fill("hTOF_length", tItrackL);
 172   }
 173   // get the TOF pid probabilities
 174   Double_t tESDpid[5] = {0., 0., 0., 0., 0.};
 175   Float_t tTOFpidSum = 0.;
 176   // find the largest PID probability
 177   track->GetTOFpid(tESDpid);
 178   Double_t tMAXpid = 0.;
 179   Int_t tMAXindex = -1;
 180   for(Int_t i=0; i<5; ++i){
 181     tTOFpidSum += tESDpid[i];
 182     if(tESDpid[i] > tMAXpid){
 183       tMAXpid = tESDpid[i];
 184       tMAXindex = i;
 185     }
 186   }
 187
 188   Int_t pdg = 0;
 189
 190   TString specname;
 191   switch(tMAXindex){
 192     case 0:    pdg = 11; specname = "electron";  break;
 193     case 1:    pdg = 13; specname = "muon"; break;
 194     case 2:    pdg = 211; specname = "pion"; break;
 195     case 3:    pdg = 321; specname = "kaon"; break;
 196     case 4:    pdg = 2212; specname = "proton"; 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     TString histname = "hTOFpid_" + specname;
 215     fQAList->Fill(histname.Data(), beta, p);
 216     fQAList->Fill("fTOFbeta_v_P_no", beta, p);
 217   }
 218   //return tMAXindex;
 219   return pdg;
 220
 221 }
 222 //__________________________________________________________________
 223 Int_t AliHFEpidTOF::MakePIDesdV2(AliESDtrack *track, AliMCParticle * /*mcTrack*/){
 224   //
 225   // Computes the PID response based on TOF & T0 signal
 226   //
 227
 228   if(!fESDpid){
 229     AliError("No ESD PID object available");
 230     return kFALSE;
 231   }
 232   Long_t status = 0;
 233   status = track->GetStatus();
 234   if(!(status & AliESDtrack::kTOFpid)) return 0;
 235
 236   Double_t p = track->GetOuterParam()->P();
 237   // track integrated times for 5 different hypothesis and T0 time
 238   Double_t times[5];
 239   track->GetIntegratedTimes(times);
 240   Double_t tItrackL = track->GetIntegratedLength();
 241   Double_t tTOFsignal = track->GetTOFsignal();
 242   Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
 243   //printf("-D: tof: %f, T0: %f\n", tTOFsignal, t0);
 244   // suppress missing or wrong T0 information
 245   if(t0 > 999990.0) return 0;
 246   Double_t tof = tTOFsignal - t0;
 247   Double_t beta = (tItrackL/100.)/((TMath::C()*tof)/1e12);
 248   //if(IsQAon())fQAList->Fill("hTOFbetaV2all", p, beta);
 249
 250   const Int_t pdg[5] = {11, 13, 211, 321, 2212};
 251   const Double_t invMass[5] = {TDatabasePDG::Instance()->GetParticle(11)->Mass(),
 252                                TDatabasePDG::Instance()->GetParticle(13)->Mass(),
 253                                TDatabasePDG::Instance()->GetParticle(211)->Mass(),
 254                                TDatabasePDG::Instance()->GetParticle(321)->Mass(),
 255                                TDatabasePDG::Instance()->GetParticle(2212)->Mass()};
 256
 257
 258   // accepted beta bands as function of momentum - parameters
 259   // line: par[0]/p + par[1] + expected_tof
 260   const Double_t bMin[5][2] = {{0., -0.03}, {-0.005, -0.02}, {-0.005, -0.02}, {-0.02, -0.006}, {-0.03, -0.005}};
 261   const Double_t bMax[5][2] = {{0., 0.03}, {0.005, 0.02}, {0.005, 0.02}, {0.02, 0.006}, {0.03, 0.005}};
 262
 263   Int_t index = -1;
 264   Double_t rdiff = 1.;
 265   for(Int_t i=0; i<5; ++i){
 266     Double_t d = (TMath::Abs(times[i] - tof))/times[i];
 267     if(d < rdiff){
 268       rdiff = d;
 269       index = i;
 270     }
 271   }
 272
 273   // stupid and unnecessary complicated - to be improved soon
 274   Double_t a = p/(invMass[index]);
 275   a *= a;
 276   Double_t betaMatch = TMath::Sqrt(a/(1+a));
 277
 278   // check wheter the most probable match is within allowed region of beta for given momentum and species
 279   Double_t min = bMin[index][0]/p + bMin[index][1] + betaMatch;
 280   Double_t max = bMax[index][0]/p + bMax[index][1] + betaMatch;
 281
 282   // debug
 283   //printf("-D: p: %f, beta: %f, pdg: %i, min: %f, max: %f\n", p, beta, pdg[index], min, max);
 284
 285   //
 286   // PID decision - can be simpler than the QA histograms above could indicate !
 287   //
 288
 289   // suppress nonsense
 290   if(beta < 0.2) return 0;
 291
 292   // 1) Simple version - protect electrons
 293   if(beta > (1+bMin[0][1]) && beta < (1+bMax[0][1])){
 294     //if(IsQAon())fQAList->Fill("hTOFbetaV2electron", p, beta);
 295     return 11;
 296   }
 297   else return 0;
 298
 299
 300   // NOT ACTIVE when (1) activated
 301   // 2) more complex version - return true PID of the particle based on the best TOF estimate
 302   // under development - still keep protecting electrons
 303   if(beta > (1+bMin[0][1]) && beta < (1+bMax[0][1])){
 304     if(IsQAon())fQAList->Fill("hTOFbetaV2_electron", p, beta);
 305     return 11;
 306   }
 307   // above 3 GeV/c the supression gets weak
 308   if(p > 3.0) return 0;
 309   if(beta > min && beta < max) {
 310     if(IsQAon())fQAList->Fill("hTOFbetaV2selected", p, beta);
 311     return pdg[index];
 312   }
 313
 314
 315   return 0;
 316 }
 317
 318 //___________________________________________________________________
 319 Int_t AliHFEpidTOF::MakePIDesdV3(AliESDtrack *track, AliMCParticle * /*mctrack*/){
 320   //
 321   // TOF PID based on n-Sigma cut
 322   // Selects Protons and Kaons via n-sigma cut up to 3 GeV/c
 323   // In addition histos for n-sigma before (all species) and after (only closest species) are filled
 324   //
 325   if(!fESDpid){
 326     AliError("No ESD pid Object available. Return");
 327     return 0;
 328   }
 329   if(!(track->GetStatus() & AliESDtrack::kTOFpid)) return 0;
 330
 331   Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
 332   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 333
 334   // Fill before selection
 335   Double_t sigEle = fESDpid->NumberOfSigmasTOF(track, AliPID::kElectron, t0);
 336   //printf("-D: p: %f, t0: %f, nSigma: %f\n", p, t0, sigEle);
 337   Int_t pdg = 0;
 338   if(TMath::Abs(sigEle) < fNsigmaTOF)
 339     pdg = 11;
 340   if(IsQAon()){
 341     // Draw TPC cleanup
 342     Double_t nsigTPC = fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron);
 343     if(nsigTPC > -3 && nsigTPC < 5) fQAList->Fill("hTOFsigmaTPCcleanup", p, sigEle);
 344     // Draw TOF signal
 345     Double_t hcontent[3] = {p, sigEle, 0};
 346     THnSparseF * hptr = dynamic_cast<THnSparseF *>(fQAList->Get("hTOFsigmaElectron"));
 347     hptr->Fill(hcontent);
 348     if(pdg == 11){
 349       hcontent[2] = 1;
 350       hptr->Fill(hcontent);
 351     }
 352   }
 353
 354   return pdg;
 355 }
 356 //___________________________________________________________________
 357 Double_t AliHFEpidTOF::Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig){
 358
 359   //gives probability for track to come from a certain particle species;
 360   //no priors considered -> probabilities for equal abundances of all species!
 361   //optional restriction to r-sigma bands of different particle species; default: rsig = 2. (see below)
 362
 363   //IMPORTANT: Tracks which are judged to be outliers get negative likelihoods -> unusable for combination with further detectors!
 364
 365   if(!track) return -1.;
 366   Bool_t outlier = kTRUE;
 367   // Check whether distance from the respective particle line is smaller than r sigma
 368   for(Int_t hypo = 0; hypo < AliPID::kSPECIES; hypo++){
 369     if(TMath::Abs(fESDpid->NumberOfSigmasTOF(track, (AliPID::EParticleType)hypo, 0.)) > rsig)
 370       outlier = kTRUE;
 371     else {
 372       outlier = kFALSE;
 373       break;
 374     }
 375   }
 376   if(outlier)
 377     return -2.;
 378
 379   Double_t tofProb[5];
 380
 381   track->GetTOFpid(tofProb);
 382
 383   return tofProb[species];
 384 }
 385 //___________________________________________________________________
 386 Int_t AliHFEpidTOF::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*mcTrack*/){
 387   AliError("AOD PID not yet implemented");
 388   return 0;
 389 }
 390
 391 //___________________________________________________________________
 392 void AliHFEpidTOF::AddQAhistograms(TList *qaList){
 393   //
 394   // Create QA histograms for TOF PID
 395   //
 396
 397   fQAList = new AliHFEcollection("TOFqaHistos", "Collection for TOF PID histograms");
 398   //fQAList->SetName("fTOFqaHistos");
 399   fQAList->CreateTH1F("hTOF_flags", "TOF flags;flags (see code for info);counts", 10, -0.25, 4.75);
 400   fQAList->CreateTH2F("fTOFbeta_v_P_no","beta -v- P; beta;momentum [GeV/c]", 120, 0, 1.2, 200, 0, 20);
 401   fQAList->CreateTH1F("hTOF_signal", "TOF signal; TOF signal [ns];counts", 1000, 12, 50);
 402   fQAList->CreateTH1F("hTOF_length", "TOF track length; length [cm];counts", 400, 300, 700);
 403   fQAList->CreateTH2F("hTOFpid_electron", "TOF reco electron; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
 404   fQAList->CreateTH2F("hTOFpid_muon", "TOF reco muon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
 405   fQAList->CreateTH2F("hTOFpid_pion", "TOF reco pion; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
 406   fQAList->CreateTH2F("hTOFpid_kaon", "TOF reco kaon; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
 407   fQAList->CreateTH2F("hTOFpid_proton", "TOF reco proton; beta ; momentum [GeV/c]", 120, 0, 1.2, 200, 0, 5);
 408   //fQAList->CreateTH2F("hTOFbetaV2all", "TOF #beta vs p for all tracks; momentum [GeV/c]; beta ", 400, 0.1, 10., 1200, 0, 1.2, 0);
 409   //fQAList->CreateTH2F("hTOFbetaV2electron", "TOF #beta vs p for selected electron tracks; momentum [GeV/c]; beta", 400, 0.1, 10., 1200, 0, 1.2, 0);
 410
 411   // histograms for sigma cut
 412   const Int_t kNdim= 3;
 413   Int_t nBins[kNdim] = {1000, 1400, 2};
 414   Double_t binMin[kNdim] = {0.1, -12., 0};
 415   Double_t binMax[kNdim] = {20, 12., 2};
 416   fQAList->CreateTHnSparse("hTOFsigmaElectron", "TOF N#sigma around the Electron Line for all tracks; p (GeV/c); N sigma; Selection Step", kNdim, nBins, binMin, binMax);
 417   fQAList->CreateTH2F("hTOFsigmaTPCcleanup", "TOF N#sigma around the Electron Line for all tracks after TPC cleanup; p (GeV/c); N sigma;", 100, 0.1, 20, 140, -12, 12);
 418
 419   qaList->AddLast(fQAList->GetList());
 420 }
 421
 422