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 "AliAODPid.h"
  31 #include "AliAODTrack.h"
  32 #include "AliAODMCParticle.h"
  33 #include "AliESDtrack.h"
  34 #include "AliLog.h"
  35 #include "AliMCParticle.h"
  36 #include "AliOADBContainer.h"
  37 #include "AliPID.h"
  38 #include "AliPIDResponse.h"
  39
  40 #include "AliHFEOADBThresholdsTRD.h"
  41 #include "AliHFEpidQAmanager.h"
  42 #include "AliHFEpidTRD.h"
  43
  44 ClassImp(AliHFEpidTRD)
  45
  46 //___________________________________________________________________
  47 AliHFEpidTRD::AliHFEpidTRD() :
  48     AliHFEpidBase()
  49   , fOADBThresholds(NULL)
  50   , fMinP(0.5)
  51   , fNTracklets(6)
  52   , fRunNumber(0)
  53   , fElectronEfficiency(0.90)
  54   , fTotalChargeInSlice0(kFALSE)
  55 {
  56   //
  57   // default  constructor
  58   //
  59   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  60 }
  61
  62 //___________________________________________________________________
  63 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
  64     AliHFEpidBase(name)
  65   , fOADBThresholds(NULL)
  66   , fMinP(0.5)
  67   , fNTracklets(6)
  68   , fRunNumber(0)
  69   , fElectronEfficiency(0.91)
  70   , fTotalChargeInSlice0(kFALSE)
  71 {
  72   //
  73   // default  constructor
  74   //
  75   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  76 }
  77
  78 //___________________________________________________________________
  79 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
  80     AliHFEpidBase("")
  81   , fOADBThresholds(NULL)
  82   , fMinP(ref.fMinP)
  83   , fNTracklets(ref.fNTracklets)
  84   , fRunNumber(ref.fRunNumber)
  85   , fElectronEfficiency(ref.fElectronEfficiency)
  86   , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
  87 {
  88   //
  89   // Copy constructor
  90   //
  91   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  92   ref.Copy(*this);
  93 }
  94
  95 //___________________________________________________________________
  96 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
  97   //
  98   // Assignment operator
  99   //
 100   if(this != &ref){
 101     ref.Copy(*this);
 102   }
 103   return *this;
 104 }
 105
 106 //___________________________________________________________________
 107 void AliHFEpidTRD::Copy(TObject &ref) const {
 108   //
 109   // Performs the copying of the object
 110   //
 111   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 112
 113   target.fMinP = fMinP;
 114   target.fNTracklets = fNTracklets;
 115   target.fRunNumber = fRunNumber;
 116   target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
 117   target.fElectronEfficiency = fElectronEfficiency;
 118   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
 119   AliHFEpidBase::Copy(ref);
 120 }
 121
 122 //___________________________________________________________________
 123 AliHFEpidTRD::~AliHFEpidTRD(){
 124   //
 125   // Destructor
 126   //
 127 }
 128
 129 //______________________________________________________
 130 Bool_t AliHFEpidTRD::InitializePID(Int_t run){
 131   //
 132   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 133   // to navigate
 134   //
 135   AliDebug(1, Form("Initializing TRD PID for run %d", run));
 136   if(InitParamsFromOADB(run)){
 137     SetBit(kThresholdsInitialized);
 138     return kTRUE;
 139   }
 140   AliDebug(1, Form("Threshold Parameters for %d tracklets and an electron efficiency %f loaded:", fNTracklets, fElectronEfficiency));
 141   AliDebug(1, Form("Params: [%f|%f|%f|%f]", fThreshParams[0], fThreshParams[1], fThreshParams[2], fThreshParams[3]));
 142   fRunNumber = run;
 143   return kFALSE;
 144 }
 145
 146 //______________________________________________________
 147 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 148   //
 149   // Does PID for TRD alone:
 150   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 151   // step function
 152   //
 153   if(!TestBit(kThresholdsInitialized)) {
 154     AliDebug(1,"Threshold Parameters not available");
 155     return 0;
 156   }
 157   AliDebug(2, "Applying TRD PID");
 158   if(!fkPIDResponse){
 159     AliDebug(2, "Cannot process track");
 160     return 0;
 161   }
 162
 163 /*
 164   const AliESDtrack *esdt = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
 165   printf("checking IdentifiedAsElectronTRD, number of Tracklets: %d\n", esdt->GetTRDntrackletsPID());
 166   if(fkPIDResponse->IdentifiedAsElectronTRD(dynamic_cast<const AliVTrack *>(track->GetRecTrack()), 0.8)) printf("Track identified as electron\n");
 167   else printf("Track rejected\n");
 168 */
 169   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
 170   Double_t p = GetP(track->GetRecTrack(), anatype);
 171   if(p < fMinP){
 172     AliDebug(2, Form("Track momentum below %f", fMinP));
 173     return 0;
 174   }
 175
 176   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
 177   Double_t electronLike = GetElectronLikelihood(static_cast<const AliVTrack *>(track->GetRecTrack()), anatype);
 178   Double_t threshold;
 179   if(TestBit(kSelectCutOnTheFly)){
 180     threshold = GetTRDthresholds(p, (static_cast<const AliVTrack *>(track->GetRecTrack()))->GetTRDntrackletsPID());
 181   } else {
 182     threshold = GetTRDthresholds(p);
 183   }
 184   AliDebug(2, Form("Threshold: %f\n", threshold));
 185   if(electronLike > threshold){
 186     if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
 187     return 11;
 188   }
 189   return 211;
 190
 191 }
 192
 193 //___________________________________________________________________
 194 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p, UInt_t nTracklets) const {
 195   //
 196   // Return momentum dependent and electron efficiency dependent TRD thresholds
 197   // Determine threshold based on the number of tracklets on the fly, electron efficiency not modified
 198   //
 199   Double_t threshParams[4];
 200   // Get threshold paramters for the given number of tracklets from OADB container
 201   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(fRunNumber));
 202   if(!thresholds){
 203     AliDebug(1, Form("Thresholds for run %d not in the OADB", fRunNumber));
 204     return 0.;
 205   }
 206   if(!thresholds->GetThresholdParameters(nTracklets, fElectronEfficiency, threshParams)) return 0.;
 207
 208   Double_t threshold = 1. - threshParams[0] - threshParams[1] * p - threshParams[2] * TMath::Exp(-threshParams[3] * p);
 209   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
 210 }
 211
 212 //___________________________________________________________________
 213 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p) const {
 214   //
 215   // Return momentum dependent and electron efficiency dependent TRD thresholds
 216   //
 217   Double_t threshold = 1. - fThreshParams[0] - fThreshParams[1] * p - fThreshParams[2] * TMath::Exp(-fThreshParams[3] * p);
 218   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
 219 }
 220
 221
 222 //___________________________________________________________________
 223 Bool_t AliHFEpidTRD::InitParamsFromOADB(Int_t run){
 224   //
 225   // The name of the function says it all
 226   //
 227   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(run));
 228   if(!thresholds){
 229     AliDebug(1, Form("Thresholds for run %d not in the OADB", run));
 230     return kFALSE;
 231   }
 232   return thresholds->GetThresholdParameters(fNTracklets, fElectronEfficiency, fThreshParams);
 233 }
 234
 235 //___________________________________________________________________
 236 void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const {
 237   //
 238   // Renormalize likelihoods for electrons and pions neglecting the
 239   // likelihoods for protons, kaons and muons
 240   //
 241   memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
 242   Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
 243   if(norm == 0.) norm = 1.;   // Safety
 244   likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
 245   likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
 246 }
 247
 248 //___________________________________________________________________
 249 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidObject::AnalysisType_t anaType) const {
 250   //
 251   // Get TRD likelihoods for ESD respectively AOD tracks
 252   //
 253   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
 254   if(anaType == AliHFEpidObject::kESDanalysis){
 255     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 256     if(esdtrack) esdtrack->GetTRDpid(pidProbs);
 257   } else {
 258     fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs);
 259   }
 260   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
 261   Double_t probsNew[AliPID::kSPECIES];
 262   RenormalizeElPi(pidProbs, probsNew);
 263   return probsNew[AliPID::kElectron];
 264 }
 265
 266 //___________________________________________________________________
 267 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 268   //
 269   // Get the Momentum in the TRD
 270   //
 271   Double_t p = 0.;
 272   if(anaType == AliHFEpidObject::kESDanalysis){
 273     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 274     if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
 275   } else {
 276     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 277     if(aodtrack) p = aodtrack->P();
 278   }
 279   return p;
 280 }
 281
 282 //___________________________________________________________________
 283 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
 284   //
 285   // Get the Charge in a single TRD layer
 286   //
 287   if(layer >= 6) return 0.;
 288   Double_t charge = 0.;
 289   if(anaType == AliHFEpidObject::kESDanalysis){
 290     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 291     if(esdtrack){
 292       // Distinction between old and new reconstruction: in the new reconstruction, the total charge is stored in slice 0, slices 1 to 8 are used for the slices for
 293       // the neural network.
 294       if(fTotalChargeInSlice0)
 295         charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
 296       else
 297        for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
 298     }
 299   } else {
 300     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 301     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
 302     if(aoddetpid){
 303       if(fTotalChargeInSlice0)
 304         charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
 305       else
 306        for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
 307     }
 308   }
 309   return charge;
 310 }
 311
 312 //___________________________________________________________________
 313 void AliHFEpidTRD::GetTRDmomenta(const AliVTrack *track, Double_t *mom) const {
 314   //
 315   // Fill Array with momentum information at the TRD tracklet
 316   //
 317   for(Int_t itl = 0; itl < 6; itl++)
 318     mom[itl] = track->GetTRDmomentum(itl);
 319 }
 320
 321 //___________________________________________________________________
 322 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
 323   //
 324   // Calculation of the TRD Signal via truncated mean
 325   // Method 1: Take all Slices available
 326   // cut out 0s
 327   // Order them in increasing order
 328   // Cut out the upper third
 329   // Calculate mean over the last 2/3 slices
 330   //
 331   const Int_t kNSlices = 48;
 332   const Int_t kSlicePerLayer = 7;
 333   const Double_t kVerySmall = 1e-12;
 334   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
 335   // Pions are used as reference for the equalization
 336   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 337   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
 338   Double_t trdSlices[kNSlices], tmp[kNSlices];
 339   Int_t indices[48];
 340   Int_t icnt = 0;
 341   for(Int_t idet = 0; idet < 6; idet++)
 342     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice < kSlicePerLayer; islice++){
 343       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 344       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
 345       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 346     }
 347   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 348   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 349   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 350   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 351   for(Int_t ien = 0; ien < icnt; ien++)
 352     trdSlices[ien] = tmp[indices[ien]];
 353   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
 354   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 355   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 356   return trdSignal;
 357 }
 358
 359 //___________________________________________________________________
 360 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
 361   //
 362   // Calculation of the TRD Signal via truncated mean
 363   // Method 2: Take only first 5 slices per chamber
 364   // Order them in increasing order
 365   // Cut out upper half
 366   // Now do mean with the reamining 3 slices per chamber
 367   //
 368   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 369   const Int_t kLayers = 6;
 370   const Int_t kSlicesLow = 6;
 371   const Int_t kSlicesHigh = 1;
 372   const Double_t kVerySmall = 1e-12;
 373   Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
 374   Int_t indices[kLayers*kSlicesLow];
 375   Int_t cntLowTime=0, cntRemaining = 0;
 376   for(Int_t idet = 0; idet < 6; idet++)
 377     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
 378       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
 379       if(islice < kSlicesLow){
 380         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 381         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 382       } else{
 383         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 384         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 385       }
 386     }
 387   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
 388   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 389   // Fill the second array with the lower half of the first time bins
 390   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
 391     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 392   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 393   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 394   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 395   return trdSignal;
 396 }