bca6a60433ce74f09fc89a98ad9247ac5813e821
[u/mrichter/AliRoot.git] / 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() : -2;
294   Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
295   Float_t centrality=-1;
296   if(centralitybin>=0) centrality=fCentralityLimitsdefault[centralitybin]+1;
297
298   if(fkPIDResponse->IdentifiedAsElectronTRD(track->GetRecTrack(),fElectronEfficiency,centrality,AliTRDPIDResponse::kLQ2D)){
299       AliDebug(2, Form("Electron effi: %f %i %i %f\n", fElectronEfficiency,track->GetCentrality(),centralitybin,centrality));
300       return 11;
301   } else return 211;
302
303
304
305 }
306
307 //___________________________________________________________________
308 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p, UInt_t nTracklets) const { 
309   //
310   // Return momentum dependent and electron efficiency dependent TRD thresholds
311   // Determine threshold based on the number of tracklets on the fly, electron efficiency not modified
312   // 
313   Double_t threshParams[4];
314   AliDebug(1, Form("Select cut for %d tracklets\n", nTracklets));
315   // Get threshold paramters for the given number of tracklets from OADB container
316   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(fRunNumber));
317   if(!thresholds){
318     AliDebug(1, Form("Thresholds for run %d not in the OADB", fRunNumber));
319     return 0.;
320   }
321   if(!thresholds->GetThresholdParameters(nTracklets, fElectronEfficiency, threshParams)){
322     AliDebug(1, "loading thresholds failed\n");
323     return 0.;
324   }
325   Double_t threshold = 1. - threshParams[0] - threshParams[1] * p - threshParams[2] * TMath::Exp(-threshParams[3] * p);
326   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
327 }
328
329 //___________________________________________________________________
330 Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p) const { 
331   //
332   // Return momentum dependent and electron efficiency dependent TRD thresholds
333   // 
334   Double_t threshold = 1. - fThreshParams[0] - fThreshParams[1] * p - fThreshParams[2] * TMath::Exp(-fThreshParams[3] * p);
335   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
336 }
337
338
339 //___________________________________________________________________
340 Bool_t AliHFEpidTRD::InitParamsFromOADB(Int_t run){
341   //
342   // The name of the function says it all
343   //
344   AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(run));
345   if(!thresholds){
346     AliDebug(1, Form("Thresholds for run %d not in the OADB", run));
347     return kFALSE;
348   }
349   return thresholds->GetThresholdParameters(fNTracklets, fElectronEfficiency, fThreshParams);
350 }
351
352 //___________________________________________________________________
353 void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const {
354   //
355   // Renormalize likelihoods for electrons and pions neglecting the 
356   // likelihoods for protons, kaons and muons
357   //
358   memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES);
359   Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion];
360   if(norm == 0.) norm = 1.;   // Safety
361   likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm;
362   likeout[AliPID::kPion] = likein[AliPID::kPion] / norm;
363 }
364
365 //___________________________________________________________________
366 Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidObject::AnalysisType_t anaType) const {
367   //
368   // Get TRD likelihoods for ESD respectively AOD tracks
369   //
370   Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES);
371   if(anaType == AliHFEpidObject::kESDanalysis){
372     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
373     if(esdtrack) esdtrack->GetTRDpid(pidProbs);
374   } else {
375       if(fTRD2DPID) fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ2D); 
376       else fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ1D);
377   }
378   if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron];
379   Double_t probsNew[AliPID::kSPECIES];
380   RenormalizeElPi(pidProbs, probsNew);
381   return probsNew[AliPID::kElectron];
382 }
383
384 //___________________________________________________________________
385 Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
386   //
387   // Get the Momentum in the TRD
388   //
389   Double_t p = 0.;
390   if(anaType == AliHFEpidObject::kESDanalysis){
391     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
392     if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P();
393   } else {
394     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
395     if(aodtrack) p = aodtrack->P();
396   }
397   return p;
398 }
399
400 //___________________________________________________________________
401 Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const {
402   //
403   // Get the Charge in a single TRD layer
404   //
405   if(layer >= 6) return 0.;
406   Double_t charge = 0.;
407   if(anaType == AliHFEpidObject::kESDanalysis){
408     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
409     if(esdtrack){
410       // 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 
411       // the neural network. 
412       if(fTotalChargeInSlice0)
413         charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
414       else
415        for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
416     }
417   } else {
418     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
419     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
420     if(aoddetpid){
421       if(fTotalChargeInSlice0)
422         charge = aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices()];
423       else
424        for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices() + islice];
425     }
426   }
427   return charge;
428 }
429
430 //___________________________________________________________________
431 void AliHFEpidTRD::GetTRDmomenta(const AliVTrack *track, Double_t *mom) const {
432   //
433   // Fill Array with momentum information at the TRD tracklet
434   //
435   for(Int_t itl = 0; itl < 6; itl++) 
436     mom[itl] = track->GetTRDmomentum(itl);
437 }
438
439 //___________________________________________________________________
440 Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const {
441   //
442   // Calculation of the TRD Signal via truncated mean
443   // Method 1: Take all Slices available
444   // cut out 0s
445   // Order them in increasing order
446   // Cut out the upper third
447   // Calculate mean over the last 2/3 slices
448   //
449   const Int_t kNSlices = 48;
450   const Int_t kLastSlice = 6; // Slice 7 is taken out from the truncated mean calculation
451   const Double_t kVerySmall = 1e-12;
452   // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice
453   // Pions are used as reference for the equalization
454   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
455   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
456   const Double_t *kWeightFactor = fTotalChargeInSlice0 ? kWeightSliceNo0 : kWeightSlice;
457   AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID()));
458   Double_t trdSlices[kNSlices], tmp[kNSlices];
459   Int_t indices[48];
460   Int_t icnt = 0;
461   for(Int_t idet = 0; idet < 6; idet++)
462     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice <= kLastSlice; islice++){
463       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
464       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
465       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightFactor[islice];
466     }
467   AliDebug(1, Form("Number of Slices: %d\n", icnt));
468   if(icnt < 6) return 0.;   // We need at least 6 Slices for the truncated mean
469   TMath::Sort(icnt, trdSlices, indices, kFALSE);
470   memcpy(tmp, trdSlices, sizeof(Double_t) * icnt);
471   for(Int_t ien = 0; ien < icnt; ien++)
472     trdSlices[ien] = tmp[indices[ien]];
473   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices);
474   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
475   AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
476   return trdSignal;
477 }
478
479 //___________________________________________________________________
480 Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const {
481   //
482   // Calculation of the TRD Signal via truncated mean
483   // Method 2: Take only first 5 slices per chamber
484   // Order them in increasing order
485   // Cut out upper half 
486   // Now do mean with the reamining 3 slices per chamber
487   //
488   const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0};
489   const Int_t kLayers = 6;
490   const Int_t kSlicesLow = 6;
491   const Int_t kSlicesHigh = 1;
492   const Double_t kVerySmall = 1e-12;
493   Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)];
494   Int_t indices[kLayers*kSlicesLow];
495   Int_t cntLowTime=0, cntRemaining = 0;
496   for(Int_t idet = 0; idet < 6; idet++)
497     for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
498       if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;;
499       if(islice < kSlicesLow){
500         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
501         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
502       } else{
503         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
504         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
505       }
506     }
507   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
508   TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE);
509   // Fill the second array with the lower half of the first time bins
510   for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++)
511     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
512   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
513   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
514   AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
515   return trdSignal;
516 }