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 //___________________________________________________________________
  46 AliHFEpidTRD::AliHFEpidTRD() :
  47     AliHFEpidBase()
  48   , fMinP(1.)
  49   , fElectronEfficiency(0.91)
  50   , fPIDMethod(kNN)
  51 {
  52   //
  53   // default  constructor
  54   //
  55   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  56   SetUseDefaultParameters();
  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   SetUseDefaultParameters();
  71 }
  72
  73 //___________________________________________________________________
  74 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
  75     AliHFEpidBase("")
  76   , fMinP(ref.fMinP)
  77   , fElectronEfficiency(ref.fElectronEfficiency)
  78   , fPIDMethod(ref.fPIDMethod)
  79 {
  80   //
  81   // Copy constructor
  82   //
  83   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  84   ref.Copy(*this);
  85 }
  86
  87 //___________________________________________________________________
  88 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
  89   //
  90   // Assignment operator
  91   //
  92   if(this != &ref){
  93     ref.Copy(*this);
  94   }
  95   return *this;
  96 }
  97
  98 //___________________________________________________________________
  99 void AliHFEpidTRD::Copy(TObject &ref) const {
 100   //
 101   // Performs the copying of the object
 102   //
 103   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 104
 105   Bool_t defaultParameters = UseDefaultParameters();
 106   target.SetUseDefaultParameters(defaultParameters);
 107   target.fMinP = fMinP;
 108   target.fPIDMethod = fPIDMethod;
 109   target.fElectronEfficiency = fElectronEfficiency;
 110   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
 111   AliHFEpidBase::Copy(ref);
 112 }
 113
 114 //___________________________________________________________________
 115 AliHFEpidTRD::~AliHFEpidTRD(){
 116   //
 117   // Destructor
 118   //
 119 }
 120
 121 //______________________________________________________
 122 Bool_t AliHFEpidTRD::InitializePID(){
 123   //
 124   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 125   // to navigate
 126   //
 127   if(UseDefaultParameters()){
 128     if(fPIDMethod == kLQ)
 129       InitParameters1DLQ();
 130     else
 131       InitParameters();
 132   }
 133   return kTRUE;
 134 }
 135
 136 //______________________________________________________
 137 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 138   //
 139   // Does PID for TRD alone:
 140   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 141   // step function
 142   //
 143   AliDebug(1, "Applying TRD PID");
 144   if((!fESDpid && track->IsESDanalysis()) || (!fAODpid && track->IsAODanalysis())){
 145     AliDebug(1, "Cannot process track");
 146     return 0;
 147   }
 148
 149   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
 150   Double_t p = GetP(track->GetRecTrack(), anatype);
 151   if(p < fMinP){
 152     AliDebug(1, Form("Track momentum below %f", fMinP));
 153     return 0;
 154   }
 155
 156   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
 157   Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
 158   Double_t threshold = GetTRDthresholds(fElectronEfficiency, p);
 159   AliDebug(1, Form("Threshold: %f\n", threshold));
 160   if(electronLike > threshold){
 161     if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
 162     return 11;
 163   }
 164   return 211;
 165
 166 }
 167
 168 //___________________________________________________________________
 169 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p) const {
 170   //
 171   // Return momentum dependent and electron efficiency dependent TRD thresholds
 172   //
 173   Double_t params[4];
 174   GetParameters(electronEff, params);
 175   Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p);
 176   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
 177 }
 178
 179 //___________________________________________________________________
 180 void AliHFEpidTRD::SetThresholdParameters(Double_t electronEff, Double_t *params){
 181   //
 182   // Set threshold parameters for the given bin
 183   //
 184   if(electronEff >= 1. || electronEff < 0.7) return;
 185   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 186   memcpy(&fThreshParams[effbin * 4], params, sizeof(Double_t) * 4);
 187   SetUseDefaultParameters(kFALSE);
 188 }
 189
 190 //___________________________________________________________________
 191 void AliHFEpidTRD::InitParameters(){
 192   //
 193   // Fill the Parameters into an array
 194   //
 195
 196   AliDebug(2, "Loading threshold Parameter");
 197   // Parameters for 6 Layers
 198   fThreshParams[0] = -0.001839; // 0.7 electron eff
 199   fThreshParams[1] = 0.000276;
 200   fThreshParams[2] = 0.044902;
 201   fThreshParams[3] = 1.726751;
 202   fThreshParams[4] = -0.002405; // 0.75 electron eff
 203   fThreshParams[5] = 0.000372;
 204   fThreshParams[6] = 0.061775;
 205   fThreshParams[7] = 1.739371;
 206   fThreshParams[8] = -0.003178; // 0.8 electron eff
 207   fThreshParams[9] = 0.000521;
 208   fThreshParams[10] = 0.087585;
 209   fThreshParams[11] = 1.749154;
 210   fThreshParams[12] = -0.004058; // 0.85 electron eff
 211   fThreshParams[13] = 0.000748;
 212   fThreshParams[14] = 0.129583;
 213   fThreshParams[15] = 1.782323;
 214   fThreshParams[16] = -0.004967; // 0.9 electron eff
 215   fThreshParams[17] = 0.001216;
 216   fThreshParams[18] = 0.210128;
 217   fThreshParams[19] = 1.807665;
 218   fThreshParams[20] = -0.000996; // 0.95 electron eff
 219   fThreshParams[21] = 0.002627;
 220   fThreshParams[22] = 0.409099;
 221   fThreshParams[23] = 1.787076;
 222 }
 223
 224 //___________________________________________________________________
 225 void AliHFEpidTRD::InitParameters1DLQ(){
 226   //
 227   // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
 228   //
 229
 230   // Parameters for 6 Layers
 231   AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
 232   fThreshParams[0] = -0.02241; // 0.7 electron eff
 233   fThreshParams[1] = 0.05043;
 234   fThreshParams[2] = 0.7925;
 235   fThreshParams[3] = 2.625;
 236   fThreshParams[4] = 0.07438; // 0.75 electron eff
 237   fThreshParams[5] = 0.05158;
 238   fThreshParams[6] = 2.864;
 239   fThreshParams[7] = 4.356;
 240   fThreshParams[8] = 0.1977; // 0.8 electron eff
 241   fThreshParams[9] = 0.05956;
 242   fThreshParams[10] = 2.853;
 243   fThreshParams[11] = 3.713;
 244   fThreshParams[12] = 0.5206; // 0.85 electron eff
 245   fThreshParams[13] = 0.03077;
 246   fThreshParams[14] = 2.966;
 247   fThreshParams[15] = 4.07;
 248   fThreshParams[16] = 0.8808; // 0.9 electron eff
 249   fThreshParams[17] = 0.002092;
 250   fThreshParams[18] = 1.17;
 251   fThreshParams[19] = 4.506;
 252   fThreshParams[20] = 1.; // 0.95 electron eff
 253   fThreshParams[21] = 0.;
 254   fThreshParams[22] = 0.;
 255   fThreshParams[23] = 0.;
 256
 257 }
 258
 259 //___________________________________________________________________
 260 void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const {
 261   //
 262   // Renormalize likelihoods for electrons and pions neglecting the
 263   // likelihoods for protons, kaons and muons
 264   //
 265   memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
 266   Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
 267   if(norm == 0.) norm = 1.;   // Safety
 268   likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
 269   likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
 270 }
 271
 272 //___________________________________________________________________
 273 void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters) const {
 274   //
 275   // return parameter set for the given efficiency bin
 276   //
 277   Int_t effbin = static_cast<Int_t>((electronEff - 0.7)/0.05);
 278   memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4);
 279 }
 280
 281 //___________________________________________________________________
 282 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 283   //
 284   // Get TRD likelihoods for ESD respectively AOD tracks
 285   //
 286   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
 287   if(anaType == AliHFEpidObject::kESDanalysis){
 288     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 289     if(esdtrack) esdtrack->GetTRDpid(pidProbs);
 290   } else {
 291     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 292     if(aodtrack)fAODpid->MakeTRDPID(const_cast<AliAODTrack *>(aodtrack), pidProbs);
 293   }
 294   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
 295   Double_t probsNew[AliPID::kSPECIES];
 296   RenormalizeElPi(pidProbs, probsNew);
 297   return probsNew[AliPID::kElectron];
 298 }
 299
 300 //___________________________________________________________________
 301 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 302   //
 303   // Get the Momentum in the TRD
 304   //
 305   Double_t p = 0.;
 306   if(anaType == AliHFEpidObject::kESDanalysis){
 307     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 308     if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
 309   } else {
 310     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 311     if(aodtrack) p = aodtrack->P();
 312   }
 313   return p;
 314 }
 315
 316 //___________________________________________________________________
 317 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
 318   //
 319   // Get the Charge in a single TRD layer
 320   //
 321   if(layer >= 6) return 0.;
 322   Double_t charge = 0.;
 323   if(anaType == AliHFEpidObject::kESDanalysis){
 324     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 325     if(esdtrack)
 326       for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
 327   } else {
 328     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 329     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
 330     if(aoddetpid)
 331       for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
 332   }
 333   return charge;
 334 }
 335
 336 //___________________________________________________________________
 337 void AliHFEpidTRD::GetTRDmomenta(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType, Double_t *mom) const {
 338   //
 339   // Fill Array with momentum information at the TRD tracklet
 340   //
 341   if(anaType == AliHFEpidObject::kESDanalysis){
 342     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 343     if(esdtrack)
 344       for(Int_t itl = 0; itl < 6; itl++)
 345         mom[itl] = esdtrack->GetTRDmomentum(itl);
 346   } else {
 347     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 348     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
 349     if(aoddetpid){
 350       Float_t *trdmom = aoddetpid->GetTRDmomentum();
 351       for(Int_t itl = 0; itl < 6; itl++){
 352         mom[itl] = trdmom[itl];
 353       }
 354     }
 355   }
 356 }
 357
 358 //___________________________________________________________________
 359 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
 360   //
 361   // Calculation of the TRD Signal via truncated mean
 362   // Method 1: Take all Slices available
 363   // cut out 0s
 364   // Order them in increasing order
 365   // Cut out the upper third
 366   // Calculate mean over the last 2/3 slices
 367   //
 368   const Int_t kNSlices = 48;
 369   const Int_t kSlicePerLayer = 7;
 370   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
 371   // Pions are used as reference for the equalization
 372   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 373   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
 374   Double_t trdSlices[kNSlices], tmp[kNSlices];
 375   Int_t indices[48];
 376   Int_t icnt = 0;
 377   for(Int_t idet = 0; idet < 6; idet++)
 378     for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
 379       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 380       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 381       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 382     }
 383   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 384   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 385   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 386   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 387   for(Int_t ien = 0; ien < icnt; ien++)
 388     trdSlices[ien] = tmp[indices[ien]];
 389   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
 390   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 391   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 392   return trdSignal;
 393 }
 394
 395 //___________________________________________________________________
 396 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
 397   //
 398   // Calculation of the TRD Signal via truncated mean
 399   // Method 2: Take only first 5 slices per chamber
 400   // Order them in increasing order
 401   // Cut out upper half
 402   // Now do mean with the reamining 3 slices per chamber
 403   //
 404   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 405   const Int_t kLayers = 6;
 406   const Int_t kSlicesLow = 6;
 407   const Int_t kSlicesHigh = 1;
 408   Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
 409   Int_t indices[kLayers*kSlicesLow];
 410   Int_t cntLowTime=0, cntRemaining = 0;
 411   for(Int_t idet = 0; idet < 6; idet++)
 412     for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
 413       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
 414       if(islice < kSlicesLow){
 415         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 416         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 417       } else{
 418         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 419         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 420       }
 421     }
 422   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
 423   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 424   // Fill the second array with the lower half of the first time bins
 425   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
 426     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 427   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 428   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 429   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 430   return trdSignal;
 431 }