PWG3/hfe/AliHFEpidTRD.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 TRD PID
  17 // Implements the abstract base class AliHFEpidbase
  18 // Make PID does the PID decision
  19 // Class further contains TRD specific cuts and QA histograms
  20 //
  21 // Authors:
  22 //   Markus Fasel <M.Fasel@gsi.de>
  23 //
  24 #include <TH1F.h>
  25 #include <TH2F.h>
  26 #include <TList.h>
  27 #include <TString.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
  36 #include "AliHFEpidTRD.h"
  37
  38 ClassImp(AliHFEpidTRD)
  39
  40 const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
  41
  42 //___________________________________________________________________
  43 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
  44     AliHFEpidBase(name)
  45   , fPIDMethod(kNN)
  46   , fContainer(0x0)
  47 {
  48   //
  49   // default  constructor
  50   //
  51   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  52 }
  53
  54 //___________________________________________________________________
  55 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
  56     AliHFEpidBase("")
  57   , fPIDMethod(kLQ)
  58   , fContainer(0x0)
  59 {
  60   //
  61   // Copy constructor
  62   //
  63   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  64   ref.Copy(*this);
  65 }
  66
  67 //___________________________________________________________________
  68 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
  69   //
  70   // Assignment operator
  71   //
  72   if(this != &ref){
  73     ref.Copy(*this);
  74   }
  75   return *this;
  76 }
  77
  78 //___________________________________________________________________
  79 void AliHFEpidTRD::Copy(TObject &ref) const {
  80   //
  81   // Performs the copying of the object
  82   //
  83   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
  84
  85   target.fPIDMethod = fPIDMethod;
  86   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
  87   if(fContainer) target.fContainer = dynamic_cast<TList *>(fContainer->Clone());
  88   AliHFEpidBase::Copy(ref);
  89 }
  90
  91 //___________________________________________________________________
  92 AliHFEpidTRD::~AliHFEpidTRD(){
  93   //
  94   // Destructor
  95   //
  96   if(fContainer){
  97     fContainer->Clear();
  98     delete fContainer;
  99   }
 100 }
 101
 102 //______________________________________________________
 103 Bool_t AliHFEpidTRD::InitializePID(){
 104   //
 105   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 106   // to navigate
 107   //
 108   InitParameters();
 109   return kTRUE;
 110 }
 111
 112 //______________________________________________________
 113 Int_t AliHFEpidTRD::IsSelected(AliHFEpidObject *track){
 114   //
 115   // Does PID for TRD alone:
 116   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 117   // step function
 118   //
 119   if(track->fAnalysisType == AliHFEpidObject::kESDanalysis){
 120       AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
 121       if(!esdTrack) return 0;
 122       AliMCParticle *esdmc = dynamic_cast<AliMCParticle *>(track->fMCtrack);
 123       return MakePIDesd(esdTrack, esdmc);
 124    } else {
 125      AliAODTrack *aodTrack = dynamic_cast<AliAODTrack *>(track->fRecTrack);
 126      if(!aodTrack) return 0;
 127      AliAODMCParticle *aodmc = dynamic_cast<AliAODMCParticle *>(track->fMCtrack);
 128      return MakePIDaod(aodTrack, aodmc);
 129    }
 130 }
 131
 132 //______________________________________________________
 133 Int_t AliHFEpidTRD::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /*aodmc*/){
 134   //
 135   // Does PID decision for AOD tracks as discribed above
 136   //
 137   AliError("AOD PID not yet implemented");
 138   return 0;
 139 }
 140
 141 //______________________________________________________
 142 Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle * /*mcTrack*/){
 143   //
 144   // Does PID decision for ESD tracks as discribed above
 145   //
 146   Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
 147   if(p < 2.0) return 0;
 148
 149   Double_t pidProbs[AliPID::kSPECIES];
 150   esdTrack->GetTRDpid(pidProbs);
 151   if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]);
 152   Double_t threshold = GetTRDthresholds(0.91, p);
 153   AliDebug(1, Form("Threshold: %f\n", threshold));
 154   if(IsQAon()) (dynamic_cast<TH2F *>(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold);
 155   if(pidProbs[AliPID::kElectron] > threshold){
 156     if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]);
 157     return 11;
 158   }
 159   return 211;
 160 }
 161
 162 //___________________________________________________________________
 163 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){
 164   //
 165   // Return momentum dependent and electron efficiency dependent TRD thresholds
 166   //
 167   Double_t params[4];
 168   GetParameters(electronEff, params);
 169   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
 170   return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values
 171 }
 172
 173 //___________________________________________________________________
 174 void AliHFEpidTRD::InitParameters(){
 175   //
 176   // Fill the Parameters into an array
 177   //
 178
 179   // Parameters for 6 Layers
 180   fThreshParams[0] = -0.001839; // 0.7 electron eff
 181   fThreshParams[1] = 0.000276;
 182   fThreshParams[2] = 0.044902;
 183   fThreshParams[3] = 1.726751;
 184   fThreshParams[4] = -0.002405; // 0.75 electron eff
 185   fThreshParams[5] = 0.000372;
 186   fThreshParams[6] = 0.061775;
 187   fThreshParams[7] = 1.739371;
 188   fThreshParams[8] = -0.003178; // 0.8 electron eff
 189   fThreshParams[9] = 0.000521;
 190   fThreshParams[10] = 0.087585;
 191   fThreshParams[11] = 1.749154;
 192   fThreshParams[12] = -0.004058; // 0.85 electron eff
 193   fThreshParams[13] = 0.000748;
 194   fThreshParams[14] = 0.129583;
 195   fThreshParams[15] = 1.782323;
 196   fThreshParams[16] = -0.004967; // 0.9 electron eff
 197   fThreshParams[17] = 0.001216;
 198   fThreshParams[18] = 0.210128;
 199   fThreshParams[19] = 1.807665;
 200   fThreshParams[20] = -0.000996; // 0.95 electron eff
 201   fThreshParams[21] = 0.002627;
 202   fThreshParams[22] = 0.409099;
 203   fThreshParams[23] = 1.787076;
 204 }
 205
 206 //___________________________________________________________________
 207 void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
 208   //
 209   // return parameter set for the given efficiency bin
 210   //
 211   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 212   memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
 213 }
 214
 215 //___________________________________________________________________
 216 Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
 217   //
 218   // Calculation of the TRD Signal via truncated mean
 219   // Method 1: Take all Slices available
 220   // cut out 0s
 221   // Order them in increasing order
 222   // Cut out the upper third
 223   // Calculate mean over the last 2/3 slices
 224   //
 225   const Int_t kNSlices = 48;
 226   AliDebug(2, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
 227   Double_t trdSlices[kNSlices], tmp[kNSlices];
 228   Int_t indices[48];
 229   Int_t icnt = 0;
 230   for(Int_t idet = 0; idet < 6; idet++)
 231     for(Int_t islice = 0; islice < 8; islice++){
 232       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 233       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 234       trdSlices[icnt++] = track->GetTRDslice(idet, islice);
 235     }
 236   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 237   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 238   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 239   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 240   for(Int_t ien = 0; ien < icnt; ien++)
 241     trdSlices[ien] = tmp[indices[ien]];
 242   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
 243   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 244   AliDebug(1, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 245   if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV1(trdSignal, mom, mcPID);
 246   return trdSignal;
 247 }
 248
 249 //___________________________________________________________________
 250 Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
 251   //
 252   // Calculation of the TRD Signal via truncated mean
 253   // Method 2: Take only first 5 slices per chamber
 254   // Order them in increasing order
 255   // Cut out upper half
 256   // Now do mean with the reamining 3 slices per chamber
 257   //
 258   Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
 259   Int_t indices[30];
 260   Int_t cntLowTime=0, cntRemaining = 0;
 261   for(Int_t idet = 0; idet < 6; idet++)
 262     for(Int_t islice = 0; islice < 8; islice++){
 263       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 264       if(islice < 5){
 265         AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 266         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
 267       } else{
 268         AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 269         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
 270       }
 271     }
 272   if(cntLowTime < 4 || cntRemaining < 2) return 0.; // Min. Number of Slices at high time is 2 (matches with 1 layer), for the truncated mean we need at least 4 Slices
 273   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 274   // Fill the second array with the lower half of the first time bins
 275   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
 276     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 277   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 278   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 279   AliDebug(1, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 280   if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV2(trdSignal, mom, mcPID);
 281   return trdSignal;
 282 }
 283
 284 //___________________________________________________________________
 285 void AliHFEpidTRD::FillHistogramsTRDSignalV1(Double_t signal, Double_t mom, Int_t species){
 286   //
 287   // Fill histograms for TRD Signal from Method 1 vs. p for different particle species
 288   //
 289   (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV1)))->Fill(mom, signal);
 290   if(species < 0 || species >= AliPID::kSPECIES) return;
 291   (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + species)))->Fill(mom, signal);
 292 }
 293
 294 //___________________________________________________________________
 295 void AliHFEpidTRD::FillHistogramsTRDSignalV2(Double_t signal, Double_t mom, Int_t species){
 296   //
 297   // Fill histograms for TRD Signal from Method 2 vs. p for different particle species
 298   //
 299   (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV2)))->Fill(mom, signal);
 300   if(species < 0 || species >= AliPID::kSPECIES) return;
 301   (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + AliPID::kSPECIES + species)))->Fill(mom, signal);
 302 }
 303
 304 //___________________________________________________________________
 305 void AliHFEpidTRD::AddQAhistograms(TList *l){
 306   //
 307   // Adding QA histograms for the TRD PID
 308   // QA histograms are:
 309   // + TRD Signal from Meth. 1 vs p for all species
 310   // + TRD Signal from Meth. 2 vs p for all species
 311   // + For each species
 312   //    - TRD Signal from Meth. 1 vs p
 313   //    - TRD Signal from Meth. 2 vs p
 314   //
 315   const Int_t kMomentumBins = 41;
 316   const Double_t kPtMin = 0.1;
 317   const Double_t kPtMax = 10.;
 318   const Int_t kSigBinsMeth1 = 100;
 319   const Int_t kSigBinsMeth2 = 100;
 320   const Double_t kMinSig = 0.;
 321   const Double_t kMaxSigMeth1 = 10000.;
 322   const Double_t kMaxSigMeth2 = 10000.;
 323
 324   if(!fContainer) fContainer = new TList;
 325   fContainer->SetName("fTRDqaHistograms");
 326
 327   Double_t momentumBins[kMomentumBins];
 328   for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
 329     momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
 330   // Likelihood Histograms
 331   fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore);
 332   fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter);
 333   fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds);
 334   // Signal Histograms
 335   fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1);
 336   fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2);
 337   for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 338     fContainer->AddAt(new TH2F(Form("fTRDSigV1%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 1)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistOverallSpecies + ispec);
 339     fContainer->AddAt(new TH2F(Form("fTRDSigV2%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 2)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth2), kHistOverallSpecies + AliPID::kSPECIES + ispec);
 340   }
 341   l->AddLast(fContainer);
 342 }
 343
 344 //___________________________________________________________________
 345 void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){
 346   //
 347   // Fill Likelihood Histogram before respectively after decision
 348   //
 349   TH2F *histo = NULL;
 350   if(whenFilled)
 351     histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeAfter));
 352   else
 353     histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeBefore));
 354   if(!histo){
 355     AliError("QA histograms not found");
 356     return;
 357   }
 358   histo->Fill(p, elProb);
 359 }