]>
Commit | Line | Data |
---|---|---|
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 | ||
44 | ClassImp(AliHFEpidTRD) | |
45 | ||
faee3b18 | 46 | //___________________________________________________________________ |
47 | AliHFEpidTRD::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 | //___________________________________________________________________ |
66 | AliHFEpidTRD::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 | //___________________________________________________________________ | |
85 | AliHFEpidTRD::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 | //___________________________________________________________________ | |
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); | |
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 | //___________________________________________________________________ | |
135 | AliHFEpidTRD::~AliHFEpidTRD(){ | |
136 | // | |
137 | // Destructor | |
138 | // | |
809a4336 | 139 | } |
140 | ||
141 | //______________________________________________________ | |
8c1c76e9 | 142 | Bool_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 | ||
b593c849 | 152 | } |
153 | ||
154 | //______________________________________________________ | |
155 | Bool_t AliHFEpidTRD::Initialize1D(Int_t run){ | |
809a4336 | 156 | // |
157 | // InitializePID: Load TRD thresholds and create the electron efficiency axis | |
158 | // to navigate | |
159 | // | |
8c1c76e9 | 160 | AliDebug(1, Form("Initializing TRD PID for run %d", run)); |
161 | if(InitParamsFromOADB(run)){ | |
162 | SetBit(kThresholdsInitialized); | |
163 | return kTRUE; | |
bf892a6a | 164 | } |
8c1c76e9 | 165 | AliDebug(1, Form("Threshold Parameters for %d tracklets and an electron efficiency %f loaded:", fNTracklets, fElectronEfficiency)); |
166 | AliDebug(1, Form("Params: [%f|%f|%f|%f]", fThreshParams[0], fThreshParams[1], fThreshParams[2], fThreshParams[3])); | |
167 | fRunNumber = run; | |
168 | return kFALSE; | |
809a4336 | 169 | } |
170 | ||
b593c849 | 171 | /* |
172 | //______________________________________________________ | |
173 | Bool_t AliHFEpidTRD::Initialize2D(Int_t run){ | |
174 | // | |
175 | // Initialize2DimPID | |
176 | // | |
177 | // | |
178 | ||
179 | return kTRUE; | |
180 | ||
181 | } | |
182 | */ | |
183 | ||
809a4336 | 184 | //______________________________________________________ |
6555e2ad | 185 | Int_t AliHFEpidTRD::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const { |
b593c849 | 186 | // |
187 | // Does PID for TRD alone: | |
188 | // PID algorithm selected according to flag | |
189 | // | |
190 | // | |
191 | ||
7bdde22f | 192 | if(fTRDOldPIDMethod) return IsSelected1D(track, pidqa); |
193 | else return IsSelectedTRDPID(track, pidqa); | |
b593c849 | 194 | } |
195 | ||
196 | //______________________________________________________ | |
197 | Int_t AliHFEpidTRD::IsSelected1D(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const { | |
809a4336 | 198 | // |
199 | // Does PID for TRD alone: | |
200 | // PID thresholds based on 90% Electron Efficiency level approximated by a linear | |
201 | // step function | |
202 | // | |
8c1c76e9 | 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"); | |
3a72645a | 210 | return 0; |
211 | } | |
722347d8 | 212 | |
e17c1f86 | 213 | |
8c1c76e9 | 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 | */ | |
3a72645a | 220 | AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis: AliHFEpidObject::kAODanalysis; |
221 | Double_t p = GetP(track->GetRecTrack(), anatype); | |
222 | if(p < fMinP){ | |
8c1c76e9 | 223 | AliDebug(2, Form("Track momentum below %f", fMinP)); |
3a72645a | 224 | return 0; |
225 | } | |
809a4336 | 226 | |
3a72645a | 227 | if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID); |
e17c1f86 | 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); | |
8c1c76e9 | 243 | Double_t threshold; |
244 | if(TestBit(kSelectCutOnTheFly)){ | |
e17c1f86 | 245 | threshold = GetTRDthresholds(p, track->GetRecTrack()->GetTRDntrackletsPID()); |
8c1c76e9 | 246 | } else { |
247 | threshold = GetTRDthresholds(p); | |
248 | } | |
249 | AliDebug(2, Form("Threshold: %f\n", threshold)); | |
3a72645a | 250 | if(electronLike > threshold){ |
251 | if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID); | |
0792aa82 | 252 | return 11; |
253 | } | |
254 | return 211; | |
3a72645a | 255 | |
809a4336 | 256 | } |
257 | ||
b593c849 | 258 | //______________________________________________________ |
93f99c26 | 259 | Int_t AliHFEpidTRD::IsSelectedTRDPID(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const { |
b593c849 | 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){ | |
7bdde22f | 274 | AliDebug(2, Form("Track momentum below %f", fMinP)); |
b593c849 | 275 | return 0; |
276 | } | |
b593c849 | 277 | |
7bdde22f | 278 | if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kBeforePID); |
279 | ||
b593c849 | 280 | if(fCutNTracklets > 0){ |
281 | AliDebug(1, Form("Number of tracklets cut applied: %d\n", fCutNTracklets)); | |
7bdde22f | 282 | Int_t ntracklets = track->GetRecTrack() ? track->GetRecTrack()->GetTRDntrackletsPID() : 0; |
b593c849 | 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 | ||
84e1e407 | 293 | Int_t centralitybin = track->IsPbPb() ? track->GetCentrality() : -2; |
b593c849 | 294 | Float_t fCentralityLimitsdefault[12]= {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.}; |
d44e3c84 | 295 | Float_t centrality=-1; |
296 | if(centralitybin>=0) centrality=fCentralityLimitsdefault[centralitybin]+1; | |
b593c849 | 297 | |
93f99c26 R |
298 | AliTRDPIDResponse::ETRDPIDMethod fTRDPIDMethod = AliTRDPIDResponse::kLQ1D; |
299 | if(fTRD2DPID) fTRDPIDMethod = AliTRDPIDResponse::kLQ2D; | |
300 | ||
7bdde22f | 301 | // printf("TRD PID Method in use: %i \n",fTRDPIDMethod); |
302 | // if(fkPIDResponse->IdentifiedAsElectronTRD(track->GetRecTrack(),fElectronEfficiency,centrality,fTRDPIDMethod)){ | |
bad10286 | 303 | Int_t ntrackletsPID=0; |
304 | Bool_t iselectron=kFALSE; | |
305 | iselectron=fkPIDResponse->IdentifiedAsElectronTRD(track->GetRecTrack(),ntrackletsPID,fElectronEfficiency,centrality,fTRDPIDMethod); | |
306 | if((ntrackletsPID==fCutNTracklets) && iselectron){ | |
93f99c26 | 307 | AliDebug(2, Form("Electron effi: %f %i %i %f %i\n", fElectronEfficiency,track->GetCentrality(),centralitybin,centrality,fTRDPIDMethod)); |
3cd4d968 | 308 | if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kTRDpid, AliHFEdetPIDqa::kAfterPID); |
b593c849 | 309 | return 11; |
310 | } else return 211; | |
311 | ||
d44e3c84 | 312 | |
313 | ||
b593c849 | 314 | } |
315 | ||
809a4336 | 316 | //___________________________________________________________________ |
8c1c76e9 | 317 | Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p, UInt_t nTracklets) const { |
dbe3abbe | 318 | // |
319 | // Return momentum dependent and electron efficiency dependent TRD thresholds | |
8c1c76e9 | 320 | // Determine threshold based on the number of tracklets on the fly, electron efficiency not modified |
dbe3abbe | 321 | // |
8c1c76e9 | 322 | Double_t threshParams[4]; |
e17c1f86 | 323 | AliDebug(1, Form("Select cut for %d tracklets\n", nTracklets)); |
8c1c76e9 | 324 | // Get threshold paramters for the given number of tracklets from OADB container |
325 | AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(fRunNumber)); | |
326 | if(!thresholds){ | |
327 | AliDebug(1, Form("Thresholds for run %d not in the OADB", fRunNumber)); | |
328 | return 0.; | |
329 | } | |
e17c1f86 | 330 | if(!thresholds->GetThresholdParameters(nTracklets, fElectronEfficiency, threshParams)){ |
331 | AliDebug(1, "loading thresholds failed\n"); | |
332 | return 0.; | |
333 | } | |
8c1c76e9 | 334 | Double_t threshold = 1. - threshParams[0] - threshParams[1] * p - threshParams[2] * TMath::Exp(-threshParams[3] * p); |
0792aa82 | 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 |
809a4336 | 336 | } |
337 | ||
bf892a6a | 338 | //___________________________________________________________________ |
8c1c76e9 | 339 | Double_t AliHFEpidTRD::GetTRDthresholds(Double_t p) const { |
bf892a6a | 340 | // |
8c1c76e9 | 341 | // Return momentum dependent and electron efficiency dependent TRD thresholds |
342 | // | |
343 | Double_t threshold = 1. - fThreshParams[0] - fThreshParams[1] * p - fThreshParams[2] * TMath::Exp(-fThreshParams[3] * p); | |
344 | 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 | 345 | } |
346 | ||
809a4336 | 347 | |
e3fc062d | 348 | //___________________________________________________________________ |
8c1c76e9 | 349 | Bool_t AliHFEpidTRD::InitParamsFromOADB(Int_t run){ |
e3fc062d | 350 | // |
8c1c76e9 | 351 | // The name of the function says it all |
e3fc062d | 352 | // |
8c1c76e9 | 353 | AliHFEOADBThresholdsTRD *thresholds = dynamic_cast<AliHFEOADBThresholdsTRD *>(fOADBThresholds->GetObject(run)); |
354 | if(!thresholds){ | |
355 | AliDebug(1, Form("Thresholds for run %d not in the OADB", run)); | |
356 | return kFALSE; | |
357 | } | |
358 | return thresholds->GetThresholdParameters(fNTracklets, fElectronEfficiency, fThreshParams); | |
e3fc062d | 359 | } |
360 | ||
c2690925 | 361 | //___________________________________________________________________ |
362 | void AliHFEpidTRD::RenormalizeElPi(const Double_t * const likein, Double_t * const likeout) const { | |
363 | // | |
364 | // Renormalize likelihoods for electrons and pions neglecting the | |
365 | // likelihoods for protons, kaons and muons | |
366 | // | |
367 | memset(likeout, 0, sizeof(Double_t) * AliPID::kSPECIES); | |
368 | Double_t norm = likein[AliPID::kElectron] + likein[AliPID::kPion]; | |
369 | if(norm == 0.) norm = 1.; // Safety | |
370 | likeout[AliPID::kElectron] = likein[AliPID::kElectron] / norm; | |
371 | likeout[AliPID::kPion] = likein[AliPID::kPion] / norm; | |
372 | } | |
373 | ||
809a4336 | 374 | //___________________________________________________________________ |
8c1c76e9 | 375 | Double_t AliHFEpidTRD::GetElectronLikelihood(const AliVTrack *track, AliHFEpidObject::AnalysisType_t anaType) const { |
3a72645a | 376 | // |
377 | // Get TRD likelihoods for ESD respectively AOD tracks | |
378 | // | |
379 | Double_t pidProbs[AliPID::kSPECIES]; memset(pidProbs, 0, sizeof(Double_t) * AliPID::kSPECIES); | |
380 | if(anaType == AliHFEpidObject::kESDanalysis){ | |
381 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
bf892a6a | 382 | if(esdtrack) esdtrack->GetTRDpid(pidProbs); |
3a72645a | 383 | } else { |
7bdde22f | 384 | if(fTRD2DPID) fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ2D); |
7a35b7aa | 385 | else fkPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, pidProbs,AliTRDPIDResponse::kLQ1D); |
3a72645a | 386 | } |
c2690925 | 387 | if(!IsRenormalizeElPi()) return pidProbs[AliPID::kElectron]; |
388 | Double_t probsNew[AliPID::kSPECIES]; | |
389 | RenormalizeElPi(pidProbs, probsNew); | |
390 | return probsNew[AliPID::kElectron]; | |
3a72645a | 391 | } |
392 | ||
393 | //___________________________________________________________________ | |
6555e2ad | 394 | Double_t AliHFEpidTRD::GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const { |
3a72645a | 395 | // |
396 | // Get the Momentum in the TRD | |
397 | // | |
398 | Double_t p = 0.; | |
399 | if(anaType == AliHFEpidObject::kESDanalysis){ | |
400 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
bf892a6a | 401 | if(esdtrack) p = esdtrack->GetOuterParam() ? esdtrack->GetOuterParam()->P() : esdtrack->P(); |
3a72645a | 402 | } else { |
403 | const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track); | |
bf892a6a | 404 | if(aodtrack) p = aodtrack->P(); |
3a72645a | 405 | } |
406 | return p; | |
407 | } | |
408 | ||
409 | //___________________________________________________________________ | |
6555e2ad | 410 | Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, AliHFEpidObject::AnalysisType_t anaType) const { |
411 | // | |
412 | // Get the Charge in a single TRD layer | |
413 | // | |
414 | if(layer >= 6) return 0.; | |
415 | Double_t charge = 0.; | |
416 | if(anaType == AliHFEpidObject::kESDanalysis){ | |
417 | const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track); | |
e156c3bb | 418 | if(esdtrack){ |
419 | // 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 | |
420 | // the neural network. | |
421 | if(fTotalChargeInSlice0) | |
422 | charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0); | |
423 | else | |
424 | for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice); | |
425 | } | |
6555e2ad | 426 | } else { |
427 | const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track); | |
bf892a6a | 428 | AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL; |
e156c3bb | 429 | if(aoddetpid){ |
430 | if(fTotalChargeInSlice0) | |
6736efd5 | 431 | charge = aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices()]; |
e156c3bb | 432 | else |
6736efd5 | 433 | for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDslices()[layer * aoddetpid->GetTRDnSlices() + islice]; |
e156c3bb | 434 | } |
6555e2ad | 435 | } |
436 | return charge; | |
437 | } | |
438 | ||
439 | //___________________________________________________________________ | |
8c1c76e9 | 440 | void AliHFEpidTRD::GetTRDmomenta(const AliVTrack *track, Double_t *mom) const { |
bf892a6a | 441 | // |
442 | // Fill Array with momentum information at the TRD tracklet | |
443 | // | |
8c1c76e9 | 444 | for(Int_t itl = 0; itl < 6; itl++) |
445 | mom[itl] = track->GetTRDmomentum(itl); | |
bf892a6a | 446 | } |
447 | ||
448 | //___________________________________________________________________ | |
449 | Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncation) const { | |
75d81601 | 450 | // |
451 | // Calculation of the TRD Signal via truncated mean | |
452 | // Method 1: Take all Slices available | |
453 | // cut out 0s | |
454 | // Order them in increasing order | |
455 | // Cut out the upper third | |
456 | // Calculate mean over the last 2/3 slices | |
457 | // | |
458 | const Int_t kNSlices = 48; | |
11ff28c5 | 459 | const Int_t kLastSlice = 6; // Slice 7 is taken out from the truncated mean calculation |
3ccf8f4c | 460 | const Double_t kVerySmall = 1e-12; |
bf892a6a | 461 | // Weight the slice to equalize the MPV of the dQ/dl-distribution per slice to the one in the first slice |
462 | // Pions are used as reference for the equalization | |
463 | const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0}; | |
11ff28c5 | 464 | 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 |
465 | const Double_t *kWeightFactor = fTotalChargeInSlice0 ? kWeightSliceNo0 : kWeightSlice; | |
bf892a6a | 466 | AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDntrackletsPID())); |
75d81601 | 467 | Double_t trdSlices[kNSlices], tmp[kNSlices]; |
468 | Int_t indices[48]; | |
469 | Int_t icnt = 0; | |
470 | for(Int_t idet = 0; idet < 6; idet++) | |
11ff28c5 | 471 | for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice <= kLastSlice; islice++){ |
75d81601 | 472 | AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice))); |
3ccf8f4c | 473 | if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;; |
11ff28c5 | 474 | trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightFactor[islice]; |
75d81601 | 475 | } |
476 | AliDebug(1, Form("Number of Slices: %d\n", icnt)); | |
477 | if(icnt < 6) return 0.; // We need at least 6 Slices for the truncated mean | |
478 | TMath::Sort(icnt, trdSlices, indices, kFALSE); | |
479 | memcpy(tmp, trdSlices, sizeof(Double_t) * icnt); | |
480 | for(Int_t ien = 0; ien < icnt; ien++) | |
481 | trdSlices[ien] = tmp[indices[ien]]; | |
bf892a6a | 482 | Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Float_t>(icnt) * truncation), trdSlices); |
75d81601 | 483 | Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1; |
e3fc062d | 484 | AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal)); |
75d81601 | 485 | return trdSignal; |
486 | } | |
487 | ||
488 | //___________________________________________________________________ | |
bf892a6a | 489 | Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncation) const { |
75d81601 | 490 | // |
491 | // Calculation of the TRD Signal via truncated mean | |
492 | // Method 2: Take only first 5 slices per chamber | |
493 | // Order them in increasing order | |
494 | // Cut out upper half | |
495 | // Now do mean with the reamining 3 slices per chamber | |
496 | // | |
bf892a6a | 497 | const Double_t kWeightSlice[8] = {1., 2.122, 1.8, 1.635, 1.595, 1.614, 1.16, 7.0}; |
498 | const Int_t kLayers = 6; | |
499 | const Int_t kSlicesLow = 6; | |
500 | const Int_t kSlicesHigh = 1; | |
3ccf8f4c | 501 | const Double_t kVerySmall = 1e-12; |
bf892a6a | 502 | Double_t trdSlicesLowTime[kLayers*kSlicesLow], trdSlicesRemaining[kLayers*(kSlicesHigh + kSlicesLow)]; |
503 | Int_t indices[kLayers*kSlicesLow]; | |
75d81601 | 504 | Int_t cntLowTime=0, cntRemaining = 0; |
505 | for(Int_t idet = 0; idet < 6; idet++) | |
e156c3bb | 506 | for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){ |
3ccf8f4c | 507 | if(TMath::Abs(track->GetTRDslice(idet, islice)) < kVerySmall) continue;; |
bf892a6a | 508 | if(islice < kSlicesLow){ |
e3fc062d | 509 | AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice))); |
bf892a6a | 510 | trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice]; |
75d81601 | 511 | } else{ |
e3fc062d | 512 | AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice))); |
bf892a6a | 513 | trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice]; |
75d81601 | 514 | } |
515 | } | |
516 | 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 | |
517 | TMath::Sort(cntLowTime, trdSlicesLowTime, indices, kFALSE); | |
518 | // Fill the second array with the lower half of the first time bins | |
bf892a6a | 519 | for(Int_t ien = 0; ien < static_cast<Int_t>(static_cast<Float_t>(cntLowTime) * truncation); ien++) |
75d81601 | 520 | trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]]; |
521 | Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining); | |
522 | Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1; | |
e3fc062d | 523 | AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal)); |
75d81601 | 524 | return trdSignal; |
525 | } |