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