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 <THnSparse.h>
  27 #include <TList.h>
  28 #include <TString.h>
  29
  30 #include "AliAODpidUtil.h"
  31 #include "AliAODPid.h"
  32 #include "AliAODTrack.h"
  33 #include "AliAODMCParticle.h"
  34 #include "AliESDtrack.h"
  35 #include "AliESDpid.h"
  36 #include "AliLog.h"
  37 #include "AliMCParticle.h"
  38 #include "AliPID.h"
  39
  40 #include "AliHFEpidQAmanager.h"
  41 #include "AliHFEpidTRD.h"
  42
  43 ClassImp(AliHFEpidTRD)
  44
  45 const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
  46
  47 //___________________________________________________________________
  48 AliHFEpidTRD::AliHFEpidTRD() :
  49     AliHFEpidBase()
  50   , fMinP(1.)
  51   , fElectronEfficiency(0.91)
  52   , fPIDMethod(kNN)
  53 {
  54   //
  55   // default  constructor
  56   //
  57 }
  58
  59 //___________________________________________________________________
  60 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
  61     AliHFEpidBase(name)
  62   , fMinP(1.)
  63   , fElectronEfficiency(0.91)
  64   , fPIDMethod(kNN)
  65 {
  66   //
  67   // default  constructor
  68   //
  69   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  70 }
  71
  72 //___________________________________________________________________
  73 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
  74     AliHFEpidBase("")
  75   , fMinP(ref.fMinP)
  76   , fElectronEfficiency(ref.fElectronEfficiency)
  77   , fPIDMethod(ref.fPIDMethod)
  78 {
  79   //
  80   // Copy constructor
  81   //
  82   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  83   ref.Copy(*this);
  84 }
  85
  86 //___________________________________________________________________
  87 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
  88   //
  89   // Assignment operator
  90   //
  91   if(this != &ref){
  92     ref.Copy(*this);
  93   }
  94   return *this;
  95 }
  96
  97 //___________________________________________________________________
  98 void AliHFEpidTRD::Copy(TObject &ref) const {
  99   //
 100   // Performs the copying of the object
 101   //
 102   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 103
 104   target.fMinP = fMinP;
 105   target.fPIDMethod = fPIDMethod;
 106   target.fElectronEfficiency = fElectronEfficiency;
 107   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
 108   AliHFEpidBase::Copy(ref);
 109 }
 110
 111 //___________________________________________________________________
 112 AliHFEpidTRD::~AliHFEpidTRD(){
 113   //
 114   // Destructor
 115   //
 116 }
 117
 118 //______________________________________________________
 119 Bool_t AliHFEpidTRD::InitializePID(){
 120   //
 121   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 122   // to navigate
 123   //
 124   if(fPIDMethod == kLQ)
 125     InitParameters1DLQ();
 126   else
 127     InitParameters();
 128   return kTRUE;
 129 }
 130
 131 //______________________________________________________
 132 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 133   //
 134   // Does PID for TRD alone:
 135   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 136   // step function
 137   //
 138   AliDebug(1, "Applying TRD PID");
 139   if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())){
 140     AliDebug(1, "Cannot process track");
 141     return 0;
 142   }
 143
 144   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
 145   Double_t p = GetP(track->GetRecTrack(), anatype);
 146   if(p < fMinP){
 147     AliDebug(1, Form("Track momentum below %f", fMinP));
 148     return 0;
 149   }
 150
 151   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
 152   Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
 153   Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
 154   AliDebug(1, Form("Threshold: %f\n", threshold));
 155   if(electronLike > threshold){
 156     if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
 157     return 11;
 158   }
 159   return 211;
 160
 161 }
 162
 163 //___________________________________________________________________
 164 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p) const {
 165   //
 166   // Return momentum dependent and electron efficiency dependent TRD thresholds
 167   //
 168   Double_t params[4];
 169   GetParameters(electronEff, params);
 170   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
 171   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
 172 }
 173
 174 //___________________________________________________________________
 175 void AliHFEpidTRD::InitParameters(){
 176   //
 177   // Fill the Parameters into an array
 178   //
 179
 180   AliDebug(2, "Loading threshold Parameter");
 181   // Parameters for 6 Layers
 182   fThreshParams[0] = -0.001839; // 0.7 electron eff
 183   fThreshParams[1] = 0.000276;
 184   fThreshParams[2] = 0.044902;
 185   fThreshParams[3] = 1.726751;
 186   fThreshParams[4] = -0.002405; // 0.75 electron eff
 187   fThreshParams[5] = 0.000372;
 188   fThreshParams[6] = 0.061775;
 189   fThreshParams[7] = 1.739371;
 190   fThreshParams[8] = -0.003178; // 0.8 electron eff
 191   fThreshParams[9] = 0.000521;
 192   fThreshParams[10] = 0.087585;
 193   fThreshParams[11] = 1.749154;
 194   fThreshParams[12] = -0.004058; // 0.85 electron eff
 195   fThreshParams[13] = 0.000748;
 196   fThreshParams[14] = 0.129583;
 197   fThreshParams[15] = 1.782323;
 198   fThreshParams[16] = -0.004967; // 0.9 electron eff
 199   fThreshParams[17] = 0.001216;
 200   fThreshParams[18] = 0.210128;
 201   fThreshParams[19] = 1.807665;
 202   fThreshParams[20] = -0.000996; // 0.95 electron eff
 203   fThreshParams[21] = 0.002627;
 204   fThreshParams[22] = 0.409099;
 205   fThreshParams[23] = 1.787076;
 206 }
 207
 208 //___________________________________________________________________
 209 void AliHFEpidTRD::InitParameters1DLQ(){
 210   //
 211   // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
 212   //
 213
 214   // Parameters for 6 Layers
 215   AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
 216   fThreshParams[0] = -0.02241; // 0.7 electron eff
 217   fThreshParams[1] = 0.05043;
 218   fThreshParams[2] = 0.7925;
 219   fThreshParams[3] = 2.625;
 220   fThreshParams[4] = 0.07438; // 0.75 electron eff
 221   fThreshParams[5] = 0.05158;
 222   fThreshParams[6] = 2.864;
 223   fThreshParams[7] = 4.356;
 224   fThreshParams[8] = 0.1977; // 0.8 electron eff
 225   fThreshParams[9] = 0.05956;
 226   fThreshParams[10] = 2.853;
 227   fThreshParams[11] = 3.713;
 228   fThreshParams[12] = 0.5206; // 0.85 electron eff
 229   fThreshParams[13] = 0.03077;
 230   fThreshParams[14] = 2.966;
 231   fThreshParams[15] = 4.07;
 232   fThreshParams[16] = 0.8808; // 0.9 electron eff
 233   fThreshParams[17] = 0.002092;
 234   fThreshParams[18] = 1.17;
 235   fThreshParams[19] = 4.506;
 236   fThreshParams[20] = 1.; // 0.95 electron eff
 237   fThreshParams[21] = 0.;
 238   fThreshParams[22] = 0.;
 239   fThreshParams[23] = 0.;
 240
 241 }
 242
 243 //___________________________________________________________________
 244 void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters) const {
 245   //
 246   // return parameter set for the given efficiency bin
 247   //
 248   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 249   memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
 250 }
 251
 252 //___________________________________________________________________
 253 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 254   //
 255   // Get TRD likelihoods for ESD respectively AOD tracks
 256   //
 257   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
 258   if(anaType == AliHFEpidObject::kESDanalysis){
 259     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 260     esdtrack->GetTRDpid(pidProbs);
 261   } else {
 262     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 263     fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
 264   }
 265   return pidProbs[AliPID::kElectron];
 266 }
 267
 268 //___________________________________________________________________
 269 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 270   //
 271   // Get the Momentum in the TRD
 272   //
 273   Double_t p = 0.;
 274   if(anaType == AliHFEpidObject::kESDanalysis){
 275     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 276     p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
 277   } else {
 278     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 279     p = aodtrack->P();
 280   }
 281   return p;
 282 }
 283
 284 //___________________________________________________________________
 285 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
 286   //
 287   // Get the Charge in a single TRD layer
 288   //
 289   if(layer >= 6) return 0.;
 290   Double_t charge = 0.;
 291   if(anaType == AliHFEpidObject::kESDanalysis){
 292     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 293     for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
 294   } else {
 295     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 296     AliAODPid *aoddetpid = aodtrack->GetDetPid();
 297     for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
 298   }
 299   return charge;
 300 }
 301
 302 //___________________________________________________________________
 303 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track) const {
 304   //
 305   // Calculation of the TRD Signal via truncated mean
 306   // Method 1: Take all Slices available
 307   // cut out 0s
 308   // Order them in increasing order
 309   // Cut out the upper third
 310   // Calculate mean over the last 2/3 slices
 311   //
 312   const Int_t kNSlices = 48;
 313   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
 314   Double_t trdSlices[kNSlices], tmp[kNSlices];
 315   Int_t indices[48];
 316   Int_t icnt = 0;
 317   for(Int_t idet = 0; idet < 6; idet++)
 318     for(Int_t islice = 0; islice < 8; islice++){
 319       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 320       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 321       trdSlices[icnt++] = track->GetTRDslice(idet, islice);
 322     }
 323   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 324   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 325   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 326   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 327   for(Int_t ien = 0; ien < icnt; ien++)
 328     trdSlices[ien] = tmp[indices[ien]];
 329   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
 330   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 331   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 332   return trdSignal;
 333 }
 334
 335 //___________________________________________________________________
 336 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track) const {
 337   //
 338   // Calculation of the TRD Signal via truncated mean
 339   // Method 2: Take only first 5 slices per chamber
 340   // Order them in increasing order
 341   // Cut out upper half
 342   // Now do mean with the reamining 3 slices per chamber
 343   //
 344   Double_t trdSlicesLowTime[30], trdSlicesRemaining[32];
 345   Int_t indices[30];
 346   Int_t cntLowTime=0, cntRemaining = 0;
 347   for(Int_t idet = 0; idet < 6; idet++)
 348     for(Int_t islice = 0; islice < 8; islice++){
 349       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 350       if(islice < 5){
 351         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 352         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
 353       } else{
 354         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 355         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
 356       }
 357     }
 358   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
 359   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 360   // Fill the second array with the lower half of the first time bins
 361   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Double_t>(cntLowTime) * 0.5); ien++)
 362     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 363   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 364   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 365   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 366   return trdSignal;
 367 }