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