PWGHF/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   , fCutNTracklets(0)
  53   , fRunNumber(0)
  54   , fElectronEfficiency(0.90)
  55   , fTotalChargeInSlice0(kFALSE)
  56   , fTRD2DPID(kFALSE)
  57 {
  58   //
  59   // default  constructor
  60   //
  61   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  62 }
  63
  64 //___________________________________________________________________
  65 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
  66     AliHFEpidBase(name)
  67   , fOADBThresholds(NULL)
  68   , fMinP(0.5)
  69   , fNTracklets(6)
  70   , fCutNTracklets(0)
  71   , fRunNumber(0)
  72   , fElectronEfficiency(0.91)
  73   , fTotalChargeInSlice0(kFALSE)
  74   , fTRD2DPID(kFALSE)
  75 {
  76   //
  77   // default  constructor
  78   //
  79   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  80 }
  81
  82 //___________________________________________________________________
  83 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
  84     AliHFEpidBase("")
  85   , fOADBThresholds(NULL)
  86   , fMinP(ref.fMinP)
  87   , fNTracklets(ref.fNTracklets)
  88   , fCutNTracklets(ref.fCutNTracklets)
  89   , fRunNumber(ref.fRunNumber)
  90   , fElectronEfficiency(ref.fElectronEfficiency)
  91   , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
  92   , fTRD2DPID(ref.fTRD2DPID)
  93 {
  94   //
  95   // Copy constructor
  96   //
  97   memset(fThreshParams, 0, sizeof(Double_t) * kThreshParams);
  98   ref.Copy(*this);
  99 }
 100
 101 //___________________________________________________________________
 102 AliHFEpidTRD &AliHFEpidTRD::operator=(const AliHFEpidTRD &ref){
 103   //
 104   // Assignment operator
 105   //
 106   if(this != &ref){
 107     ref.Copy(*this);
 108   }
 109   return *this;
 110 }
 111
 112 //___________________________________________________________________
 113 void AliHFEpidTRD::Copy(TObject &ref) const {
 114   //
 115   // Performs the copying of the object
 116   //
 117   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 118
 119   target.fMinP = fMinP;
 120   target.fNTracklets = fNTracklets;
 121   target.fCutNTracklets = fCutNTracklets;
 122   target.fRunNumber = fRunNumber;
 123   target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
 124   target.fTRD2DPID = fTRD2DPID;
 125   target.fElectronEfficiency = fElectronEfficiency;
 126   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
 127   AliHFEpidBase::Copy(ref);
 128 }
 129
 130 //___________________________________________________________________
 131 AliHFEpidTRD::~AliHFEpidTRD(){
 132   //
 133   // Destructor
 134   //
 135 }
 136
 137 //______________________________________________________
 138 Bool_t AliHFEpidTRD::InitializePID(Int_t run){
 139   //
 140   // InitializePID call different init function depending on TRD PID method
 141   //
 142   //
 143
 144   //if(fTRD2DPID) return Initialize2D(run);
 145   if(fTRD2DPID) return kTRUE;
 146   else return Initialize1D(run);
 147
 148
 149
 150 }
 151
 152 //______________________________________________________
 153 Bool_t AliHFEpidTRD::Initialize1D(Int_t run){
 154   //
 155   // InitializePID: Load TRD thresholds and create the electron efficiency axis
 156   // to navigate
 157   //
 158   AliDebug(1, Form("Initializing TRD PID for run %d", run));
 159   if(InitParamsFromOADB(run)){
 160     SetBit(kThresholdsInitialized);
 161     return kTRUE;
 162   }
 163   AliDebug(1, Form("Threshold Parameters for %d tracklets and an electron efficiency %f loaded:", fNTracklets, fElectronEfficiency));
 164   AliDebug(1, Form("Params: [%f|%f|%f|%f]", fThreshParams[0], fThreshParams[1], fThreshParams[2], fThreshParams[3]));
 165   fRunNumber = run;
 166   return kFALSE;
 167 }
 168
 169 /*
 170 //______________________________________________________
 171 Bool_t AliHFEpidTRD::Initialize2D(Int_t run){
 172   //
 173   // Initialize2DimPID
 174   //
 175   //
 176
 177     return kTRUE;
 178
 179 }
 180 */
 181
 182 //______________________________________________________
 183 Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 184   //
 185   // Does PID for TRD alone:
 186   // PID algorithm selected according to flag
 187   //
 188   //
 189
 190     if(fTRD2DPID) return IsSelected2D(track, pidqa);
 191     else return IsSelected1D(track, pidqa);
 192
 193
 194 }
 195
 196 //______________________________________________________
 197 Int_t AliHFEpidTRD::IsSelected1D(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 198   //
 199   // Does PID for TRD alone:
 200   // PID thresholds based on 90% Electron Efficiency level approximated by a linear
 201   // step function
 202   //
 203   if(!TestBit(kThresholdsInitialized)) {
 204     AliDebug(1,"Threshold Parameters not available");
 205     return 0;
 206   }
 207   AliDebug(2, "Applying TRD PID");
 208   if(!fkPIDResponse){
 209     AliDebug(2, "Cannot process track");
 210     return 0;
 211   }
 212
 213
 214 /*
 215   const AliESDtrack *esdt = dynamic_cast<const AliESDtrack *>(track->GetRecTrack());
 216   printf("checking IdentifiedAsElectronTRD, number of Tracklets: %d\n", esdt->GetTRDntrackletsPID());
 217   if(fkPIDResponse->IdentifiedAsElectronTRD(dynamic_cast<const AliVTrack *>(track->GetRecTrack()), 0.8)) printf("Track identified as electron\n");
 218   else printf("Track rejected\n");
 219 */
 220   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
 221   Double_t p = GetP(track->GetRecTrack(), anatype);
 222   if(p < fMinP){
 223     AliDebug(2, Form("Track momentum below %f", fMinP));
 224     return 0;
 225   }
 226
 227   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
 228
 229   if(fCutNTracklets > 0){
 230     AliDebug(1, Form("Number of tracklets cut applied: %d\n", fCutNTracklets));
 231     Int_t ntracklets = track->GetRecTrack() ? track->GetRecTrack()->GetTRDntrackletsPID() : 0;
 232     if(TestBit(kExactTrackletCut)){
 233       AliDebug(1, Form("Exact cut applied: %d tracklets found\n", ntracklets));
 234       if(ntracklets != fCutNTracklets) return 0;
 235     } else {
 236       AliDebug(1, Form("Greater Equal cut applied: %d tracklets found\n", ntracklets));
 237       if(ntracklets < fCutNTracklets) return 0;
 238     }
 239   }
 240   AliDebug(1,"Track selected\n");
 241
 242   Double_t electronLike = GetElectronLikelihood(track->GetRecTrack(), anatype);
 243   Double_t threshold;
 244   if(TestBit(kSelectCutOnTheFly)){
 245     threshold = GetTRDthresholds(p, track->GetRecTrack()->GetTRDntrackletsPID());
 246   } else {
 247     threshold = GetTRDthresholds(p);
 248   }
 249   AliDebug(2, Form("Threshold: %f\n", threshold));
 250   if(electronLike > threshold){
 251     if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID);
 252     return 11;
 253   }
 254   return 211;
 255
 256 }
 257
 258 //______________________________________________________
 259 Int_t AliHFEpidTRD::IsSelected2D(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const {
 260   //
 261   // 2D TRD PID
 262   //
 263   //
 264   //
 265   AliDebug(2, "Applying TRD PID");
 266   if(!fkPIDResponse){
 267     AliDebug(2, "Cannot process track");
 268     return 0;
 269   }
 270
 271   AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis;
 272   Double_t p = GetP(track->GetRecTrack(), anatype);
 273   if(p < fMinP){
 274     AliDebug(2, Form("Track momentum below %f", fMinP));
 275     return 0;
 276   }
 277
 278   if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID);
 279
 280   if(fCutNTracklets > 0){
 281     AliDebug(1, Form("Number of tracklets cut applied: %d\n", fCutNTracklets));
 282     Int_t ntracklets = track->GetRecTrack() ? track->GetRecTrack()->GetTRDntrackletsPID() : 0;
 283     if(TestBit(kExactTrackletCut)){
 284       AliDebug(1, Form("Exact cut applied: %d tracklets found\n", ntracklets));
 285       if(ntracklets != fCutNTracklets) return 0;
 286     } else {
 287       AliDebug(1, Form("Greater Equal cut applied: %d tracklets found\n", ntracklets));
 288       if(ntracklets < fCutNTracklets) return 0;
 289     }
 290   }
 291   AliDebug(1,"Track selected\n");
 292
 293   Int_t centralitybin = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
 294   Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
 295   Float_t centrality=fCentralityLimitsdefault[centralitybin]+1;
 296
 297
 298   if(fkPIDResponse->IdentifiedAsElectronTRD(track->GetRecTrack(),fElectronEfficiency,centrality,AliTRDPIDResponse::kLQ2D)){
 299       AliDebug(2, Form("Electron effi: %f\n", fElectronEfficiency));
 300       return 11;
 301   } else return 211;
 302
 303 }
 304
 305 //___________________________________________________________________
 306 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p, UInt_t nTracklets) const {
 307   //
 308   // Return momentum dependent and electron efficiency dependent TRD thresholds
 309   // Determine threshold based on the number of tracklets on the fly, electron efficiency not modified
 310   //
 311   Double_t threshParams[4];
 312   AliDebug(1, Form("Select cut for %d tracklets\n", nTracklets));
 313   // Get threshold paramters for the given number of tracklets from OADB container
 314   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(fRunNumber));
 315   if(!thresholds){
 316     AliDebug(1, Form("Thresholds for run %d not in the OADB", fRunNumber));
 317     return 0.;
 318   }
 319   if(!thresholds->GetThresholdParameters(nTracklets, fElectronEfficiency, threshParams)){
 320     AliDebug(1, "loading thresholds failed\n");
 321     return 0.;
 322   }
 323   Double_t threshold = 1. - threshParams[0] - threshParams[1] * p - threshParams[2] * TMath::Exp(-threshParams[3] * p);
 324   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
 325 }
 326
 327 //___________________________________________________________________
 328 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p) const {
 329   //
 330   // Return momentum dependent and electron efficiency dependent TRD thresholds
 331   //
 332   Double_t threshold = 1. - fThreshParams[0] - fThreshParams[1] * p - fThreshParams[2] * TMath::Exp(-fThreshParams[3] * p);
 333   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
 334 }
 335
 336
 337 //___________________________________________________________________
 338 Bool_t AliHFEpidTRD::InitParamsFromOADB(Int_t run){
 339   //
 340   // The name of the function says it all
 341   //
 342   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(run));
 343   if(!thresholds){
 344     AliDebug(1, Form("Thresholds for run %d not in the OADB", run));
 345     return kFALSE;
 346   }
 347   return thresholds->GetThresholdParameters(fNTracklets, fElectronEfficiency, fThreshParams);
 348 }
 349
 350 //___________________________________________________________________
 351 void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const {
 352   //
 353   // Renormalize likelihoods for electrons and pions neglecting the
 354   // likelihoods for protons, kaons and muons
 355   //
 356   memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
 357   Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
 358   if(norm == 0.) norm = 1.;   // Safety
 359   likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
 360   likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
 361 }
 362
 363 //___________________________________________________________________
 364 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidObject::AnalysisType_t anaType) const {
 365   //
 366   // Get TRD likelihoods for ESD respectively AOD tracks
 367   //
 368   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
 369   if(anaType == AliHFEpidObject::kESDanalysis){
 370     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 371     if(esdtrack) esdtrack->GetTRDpid(pidProbs);
 372   } else {
 373     fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs);
 374   }
 375   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
 376   Double_t probsNew[AliPID::kSPECIES];
 377   RenormalizeElPi(pidProbs, probsNew);
 378   return probsNew[AliPID::kElectron];
 379 }
 380
 381 //___________________________________________________________________
 382 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
 383   //
 384   // Get the Momentum in the TRD
 385   //
 386   Double_t p = 0.;
 387   if(anaType == AliHFEpidObject::kESDanalysis){
 388     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 389     if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
 390   } else {
 391     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 392     if(aodtrack) p = aodtrack->P();
 393   }
 394   return p;
 395 }
 396
 397 //___________________________________________________________________
 398 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
 399   //
 400   // Get the Charge in a single TRD layer
 401   //
 402   if(layer >= 6) return 0.;
 403   Double_t charge = 0.;
 404   if(anaType == AliHFEpidObject::kESDanalysis){
 405     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
 406     if(esdtrack){
 407       // 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
 408       // the neural network.
 409       if(fTotalChargeInSlice0)
 410         charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
 411       else
 412        for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
 413     }
 414   } else {
 415     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
 416     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
 417     if(aoddetpid){
 418       if(fTotalChargeInSlice0)
 419         charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
 420       else
 421        for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
 422     }
 423   }
 424   return charge;
 425 }
 426
 427 //___________________________________________________________________
 428 void AliHFEpidTRD::GetTRDmomenta(const AliVTrack *track, Double_t *mom) const {
 429   //
 430   // Fill Array with momentum information at the TRD tracklet
 431   //
 432   for(Int_t itl = 0; itl < 6; itl++)
 433     mom[itl] = track->GetTRDmomentum(itl);
 434 }
 435
 436 //___________________________________________________________________
 437 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
 438   //
 439   // Calculation of the TRD Signal via truncated mean
 440   // Method 1: Take all Slices available
 441   // cut out 0s
 442   // Order them in increasing order
 443   // Cut out the upper third
 444   // Calculate mean over the last 2/3 slices
 445   //
 446   const Int_t kNSlices = 48;
 447   const Int_t kLastSlice = 6; // Slice 7 is taken out from the truncated mean calculation
 448   const Double_t kVerySmall = 1e-12;
 449   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
 450   // Pions are used as reference for the equalization
 451   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 452   const Double_t kWeightSliceNo0[8] = {1., 1., 1.271, 1.451, 1.531, 1.543, 1.553, 2.163};  // Weighting factors in case slice 0 stores the total charge
 453   const Double_t *kWeightFactor = fTotalChargeInSlice0 ? kWeightSliceNo0 : kWeightSlice;
 454   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
 455   Double_t trdSlices[kNSlices], tmp[kNSlices];
 456   Int_t indices[48];
 457   Int_t icnt = 0;
 458   for(Int_t idet = 0; idet < 6; idet++)
 459     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice <= kLastSlice; islice++){
 460       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
 461       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
 462       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightFactor[islice];
 463     }
 464   AliDebug(1, Form("Number of Slices: %d\n", icnt));
 465   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
 466   TMath::Sort(icnt, trdSlices, indices, kFALSE);
 467   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
 468   for(Int_t ien = 0; ien < icnt; ien++)
 469     trdSlices[ien] = tmp[indices[ien]];
 470   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
 471   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 472   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
 473   return trdSignal;
 474 }
 475
 476 //___________________________________________________________________
 477 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
 478   //
 479   // Calculation of the TRD Signal via truncated mean
 480   // Method 2: Take only first 5 slices per chamber
 481   // Order them in increasing order
 482   // Cut out upper half
 483   // Now do mean with the reamining 3 slices per chamber
 484   //
 485   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
 486   const Int_t kLayers = 6;
 487   const Int_t kSlicesLow = 6;
 488   const Int_t kSlicesHigh = 1;
 489   const Double_t kVerySmall = 1e-12;
 490   Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
 491   Int_t indices[kLayers*kSlicesLow];
 492   Int_t cntLowTime=0, cntRemaining = 0;
 493   for(Int_t idet = 0; idet < 6; idet++)
 494     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
 495       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
 496       if(islice < kSlicesLow){
 497         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 498         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 499       } else{
 500         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
 501         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
 502       }
 503     }
 504   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
 505   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
 506   // Fill the second array with the lower half of the first time bins
 507   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
 508     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
 509   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
 510   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
 511   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
 512   return trdSignal;
 513 }