]>
Commit | Line | Data |
---|---|---|
29bf19f2 | 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 | ||
4ec8e76d | 16 | /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */ |
29bf19f2 | 17 | |
18 | //----------------------------------------------------------------- | |
4ec8e76d | 19 | // Base class for handling the pid response // |
20 | // functions of all detectors // | |
21 | // and give access to the nsigmas // | |
22 | // // | |
23 | // Origin: Jens Wiechula, Uni Tuebingen, jens.wiechula@cern.ch // | |
29bf19f2 | 24 | //----------------------------------------------------------------- |
25 | ||
4ec8e76d | 26 | #include <TList.h> |
27 | #include <TObjArray.h> | |
28 | #include <TPRegexp.h> | |
29 | #include <TF1.h> | |
30 | #include <TSpline.h> | |
31 | #include <TFile.h> | |
00a38d07 | 32 | #include <TArrayI.h> |
db0e2c5f | 33 | #include <TArrayF.h> |
4ec8e76d | 34 | |
35 | #include <AliVEvent.h> | |
fd21ec8d | 36 | #include <AliVTrack.h> |
4ec8e76d | 37 | #include <AliLog.h> |
38 | #include <AliPID.h> | |
ea235c90 | 39 | #include <AliOADBContainer.h> |
db0e2c5f | 40 | #include <AliTRDPIDResponseObject.h> |
b79db598 | 41 | #include <AliTOFPIDParams.h> |
29bf19f2 | 42 | |
43 | #include "AliPIDResponse.h" | |
00a38d07 | 44 | #include "AliDetectorPID.h" |
29bf19f2 | 45 | |
80f28562 | 46 | #include "AliCentrality.h" |
47 | ||
29bf19f2 | 48 | ClassImp(AliPIDResponse); |
49 | ||
4ec8e76d | 50 | AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) : |
51 | TNamed("PIDResponse","PIDResponse"), | |
52 | fITSResponse(isMC), | |
53 | fTPCResponse(), | |
54 | fTRDResponse(), | |
55 | fTOFResponse(), | |
e96b9916 | 56 | fEMCALResponse(), |
fd21ec8d | 57 | fRange(5.), |
58 | fITSPIDmethod(kITSTruncMean), | |
4ec8e76d | 59 | fIsMC(isMC), |
1c9d11be | 60 | fCachePID(kTRUE), |
4ec8e76d | 61 | fOADBPath(), |
00a38d07 | 62 | fCustomTPCpidResponse(), |
4ec8e76d | 63 | fBeamType("PP"), |
64 | fLHCperiod(), | |
65 | fMCperiodTPC(), | |
fd21ec8d | 66 | fMCperiodUser(), |
ea235c90 | 67 | fCurrentFile(), |
4ec8e76d | 68 | fRecoPass(0), |
fd21ec8d | 69 | fRecoPassUser(-1), |
4ec8e76d | 70 | fRun(0), |
71 | fOldRun(0), | |
78cbd205 | 72 | fResT0A(75.), |
73 | fResT0C(65.), | |
74 | fResT0AC(55.), | |
644666df | 75 | fArrPidResponseMaster(NULL), |
76 | fResolutionCorrection(NULL), | |
77 | fOADBvoltageMaps(NULL), | |
78 | fTRDPIDResponseObject(NULL), | |
0b39f221 | 79 | fTOFtail(1.1), |
644666df | 80 | fTOFPIDParams(NULL), |
81 | fEMCALPIDParams(NULL), | |
82 | fCurrentEvent(NULL), | |
539a5a59 | 83 | fCurrCentrality(0.0), |
84 | fTuneMConData(kFALSE) | |
4ec8e76d | 85 | { |
86 | // | |
87 | // default ctor | |
88 | // | |
a635821f | 89 | AliLog::SetClassDebugLevel("AliPIDResponse",0); |
90 | AliLog::SetClassDebugLevel("AliESDpid",0); | |
91 | AliLog::SetClassDebugLevel("AliAODpidUtil",0); | |
ea235c90 | 92 | |
4ec8e76d | 93 | } |
94 | ||
95 | //______________________________________________________________________________ | |
96 | AliPIDResponse::~AliPIDResponse() | |
97 | { | |
98 | // | |
99 | // dtor | |
100 | // | |
00a38d07 | 101 | delete fArrPidResponseMaster; |
102 | delete fTRDPIDResponseObject; | |
103 | delete fTOFPIDParams; | |
4ec8e76d | 104 | } |
105 | ||
106 | //______________________________________________________________________________ | |
107 | AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) : | |
108 | TNamed(other), | |
109 | fITSResponse(other.fITSResponse), | |
110 | fTPCResponse(other.fTPCResponse), | |
111 | fTRDResponse(other.fTRDResponse), | |
112 | fTOFResponse(other.fTOFResponse), | |
e96b9916 | 113 | fEMCALResponse(other.fEMCALResponse), |
fd21ec8d | 114 | fRange(other.fRange), |
115 | fITSPIDmethod(other.fITSPIDmethod), | |
4ec8e76d | 116 | fIsMC(other.fIsMC), |
1c9d11be | 117 | fCachePID(other.fCachePID), |
4ec8e76d | 118 | fOADBPath(other.fOADBPath), |
00a38d07 | 119 | fCustomTPCpidResponse(other.fCustomTPCpidResponse), |
4ec8e76d | 120 | fBeamType("PP"), |
121 | fLHCperiod(), | |
122 | fMCperiodTPC(), | |
fd21ec8d | 123 | fMCperiodUser(other.fMCperiodUser), |
ea235c90 | 124 | fCurrentFile(), |
4ec8e76d | 125 | fRecoPass(0), |
fd21ec8d | 126 | fRecoPassUser(other.fRecoPassUser), |
4ec8e76d | 127 | fRun(0), |
128 | fOldRun(0), | |
78cbd205 | 129 | fResT0A(75.), |
130 | fResT0C(65.), | |
131 | fResT0AC(55.), | |
644666df | 132 | fArrPidResponseMaster(NULL), |
133 | fResolutionCorrection(NULL), | |
134 | fOADBvoltageMaps(NULL), | |
135 | fTRDPIDResponseObject(NULL), | |
0b39f221 | 136 | fTOFtail(1.1), |
644666df | 137 | fTOFPIDParams(NULL), |
138 | fEMCALPIDParams(NULL), | |
139 | fCurrentEvent(NULL), | |
539a5a59 | 140 | fCurrCentrality(0.0), |
141 | fTuneMConData(kFALSE) | |
4ec8e76d | 142 | { |
143 | // | |
144 | // copy ctor | |
145 | // | |
146 | } | |
147 | ||
148 | //______________________________________________________________________________ | |
149 | AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other) | |
150 | { | |
151 | // | |
152 | // copy ctor | |
153 | // | |
154 | if(this!=&other) { | |
155 | delete fArrPidResponseMaster; | |
156 | TNamed::operator=(other); | |
157 | fITSResponse=other.fITSResponse; | |
158 | fTPCResponse=other.fTPCResponse; | |
159 | fTRDResponse=other.fTRDResponse; | |
160 | fTOFResponse=other.fTOFResponse; | |
e96b9916 | 161 | fEMCALResponse=other.fEMCALResponse; |
fd21ec8d | 162 | fRange=other.fRange; |
163 | fITSPIDmethod=other.fITSPIDmethod; | |
4ec8e76d | 164 | fOADBPath=other.fOADBPath; |
00a38d07 | 165 | fCustomTPCpidResponse=other.fCustomTPCpidResponse; |
4ec8e76d | 166 | fIsMC=other.fIsMC; |
1c9d11be | 167 | fCachePID=other.fCachePID; |
4ec8e76d | 168 | fBeamType="PP"; |
169 | fLHCperiod=""; | |
170 | fMCperiodTPC=""; | |
fd21ec8d | 171 | fMCperiodUser=other.fMCperiodUser; |
ea235c90 | 172 | fCurrentFile=""; |
4ec8e76d | 173 | fRecoPass=0; |
fd21ec8d | 174 | fRecoPassUser=other.fRecoPassUser; |
4ec8e76d | 175 | fRun=0; |
176 | fOldRun=0; | |
78cbd205 | 177 | fResT0A=75.; |
178 | fResT0C=65.; | |
179 | fResT0AC=55.; | |
644666df | 180 | fArrPidResponseMaster=NULL; |
181 | fResolutionCorrection=NULL; | |
182 | fOADBvoltageMaps=NULL; | |
183 | fTRDPIDResponseObject=NULL; | |
184 | fEMCALPIDParams=NULL; | |
0b39f221 | 185 | fTOFtail=1.1; |
644666df | 186 | fTOFPIDParams=NULL; |
e96b9916 | 187 | fCurrentEvent=other.fCurrentEvent; |
bd58d4b9 | 188 | |
4ec8e76d | 189 | } |
190 | return *this; | |
191 | } | |
192 | ||
fd21ec8d | 193 | //______________________________________________________________________________ |
194 | Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const | |
195 | { | |
196 | // | |
197 | // NumberOfSigmas for 'detCode' | |
198 | // | |
199 | ||
200 | switch (detCode){ | |
201 | case kDetITS: return NumberOfSigmasITS(track, type); break; | |
202 | case kDetTPC: return NumberOfSigmasTPC(track, type); break; | |
203 | case kDetTOF: return NumberOfSigmasTOF(track, type); break; | |
00a38d07 | 204 | case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break; |
fd21ec8d | 205 | default: return -999.; |
206 | } | |
207 | ||
208 | } | |
209 | ||
e96b9916 | 210 | //______________________________________________________________________________ |
00a38d07 | 211 | Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const |
212 | { | |
213 | // | |
214 | // NumberOfSigmas for 'detCode' | |
215 | // | |
216 | return NumberOfSigmas((EDetCode)(1<<detCode), track, type); | |
217 | } | |
218 | ||
1c9d11be | 219 | //______________________________________________________________________________ |
220 | // public buffered versions of the PID calculation | |
221 | // | |
222 | ||
00a38d07 | 223 | //______________________________________________________________________________ |
224 | Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
225 | { | |
226 | // | |
227 | // Calculate the number of sigmas in the ITS | |
228 | // | |
229 | ||
230 | AliVTrack *track=(AliVTrack*)vtrack; | |
231 | ||
232 | // look for cached value first | |
1c9d11be | 233 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
234 | const EDetector detector=kITS; | |
235 | ||
236 | if ( detPID && detPID->HasNumberOfSigmas(detector)){ | |
237 | return detPID->GetNumberOfSigmas(detector, type); | |
238 | } else if (fCachePID) { | |
239 | FillTrackDetectorPID(track, detector); | |
240 | detPID=track->GetDetectorPID(); | |
241 | return detPID->GetNumberOfSigmas(detector, type); | |
00a38d07 | 242 | } |
00a38d07 | 243 | |
1c9d11be | 244 | return GetNumberOfSigmasITS(track, type); |
00a38d07 | 245 | } |
246 | ||
247 | //______________________________________________________________________________ | |
248 | Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
249 | { | |
250 | // | |
251 | // Calculate the number of sigmas in the TPC | |
252 | // | |
253 | ||
254 | AliVTrack *track=(AliVTrack*)vtrack; | |
255 | ||
256 | // look for cached value first | |
1c9d11be | 257 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
258 | const EDetector detector=kTPC; | |
259 | ||
260 | if ( detPID && detPID->HasNumberOfSigmas(detector)){ | |
261 | return detPID->GetNumberOfSigmas(detector, type); | |
262 | } else if (fCachePID) { | |
263 | FillTrackDetectorPID(track, detector); | |
264 | detPID=track->GetDetectorPID(); | |
265 | return detPID->GetNumberOfSigmas(detector, type); | |
00a38d07 | 266 | } |
267 | ||
1c9d11be | 268 | return GetNumberOfSigmasTPC(track, type); |
00a38d07 | 269 | } |
270 | ||
644666df | 271 | //______________________________________________________________________________ |
272 | Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack, | |
273 | AliPID::EParticleType type, | |
274 | AliTPCPIDResponse::ETPCdEdxSource dedxSource) | |
275 | { | |
276 | //get number of sigmas according the selected TPC gain configuration scenario | |
277 | const AliVTrack *track=static_cast<const AliVTrack*>(vtrack); | |
278 | ||
279 | Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource); | |
280 | ||
281 | return nSigma; | |
282 | } | |
283 | ||
00a38d07 | 284 | //______________________________________________________________________________ |
1c9d11be | 285 | Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const |
00a38d07 | 286 | { |
287 | // | |
1c9d11be | 288 | // Calculate the number of sigmas in the TOF |
00a38d07 | 289 | // |
290 | ||
291 | AliVTrack *track=(AliVTrack*)vtrack; | |
1c9d11be | 292 | |
00a38d07 | 293 | // look for cached value first |
1c9d11be | 294 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
295 | const EDetector detector=kTOF; | |
296 | ||
297 | if ( detPID && detPID->HasNumberOfSigmas(detector)){ | |
298 | return detPID->GetNumberOfSigmas(detector, type); | |
299 | } else if (fCachePID) { | |
300 | FillTrackDetectorPID(track, detector); | |
301 | detPID=track->GetDetectorPID(); | |
302 | return detPID->GetNumberOfSigmas(detector, type); | |
00a38d07 | 303 | } |
304 | ||
1c9d11be | 305 | return GetNumberOfSigmasTOF(track, type); |
306 | } | |
e96b9916 | 307 | |
1c9d11be | 308 | //______________________________________________________________________________ |
309 | Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
310 | { | |
311 | // | |
312 | // Calculate the number of sigmas in the EMCAL | |
313 | // | |
e96b9916 | 314 | |
1c9d11be | 315 | AliVTrack *track=(AliVTrack*)vtrack; |
316 | ||
317 | // look for cached value first | |
318 | const AliDetectorPID *detPID=track->GetDetectorPID(); | |
319 | const EDetector detector=kEMCAL; | |
320 | ||
321 | if ( detPID && detPID->HasNumberOfSigmas(detector)){ | |
322 | return detPID->GetNumberOfSigmas(detector, type); | |
323 | } else if (fCachePID) { | |
324 | FillTrackDetectorPID(track, detector); | |
325 | detPID=track->GetDetectorPID(); | |
326 | return detPID->GetNumberOfSigmas(detector, type); | |
e96b9916 | 327 | } |
6d0064aa | 328 | |
1c9d11be | 329 | return GetNumberOfSigmasEMCAL(track, type); |
e96b9916 | 330 | } |
331 | ||
6d0064aa | 332 | //______________________________________________________________________________ |
1c9d11be | 333 | Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const |
334 | { | |
335 | // | |
336 | // emcal nsigma with eop and showershape | |
337 | // | |
00a38d07 | 338 | AliVTrack *track=(AliVTrack*)vtrack; |
339 | ||
6d0064aa | 340 | AliVCluster *matchedClus = NULL; |
341 | ||
342 | Double_t mom = -1.; | |
343 | Double_t pt = -1.; | |
344 | Double_t EovP = -1.; | |
345 | Double_t fClsE = -1.; | |
32fa24d6 | 346 | |
347 | // initialize eop and shower shape parameters | |
348 | eop = -1.; | |
349 | for(Int_t i = 0; i < 4; i++){ | |
350 | showershape[i] = -1.; | |
351 | } | |
6d0064aa | 352 | |
353 | Int_t nMatchClus = -1; | |
354 | Int_t charge = 0; | |
355 | ||
356 | // Track matching | |
357 | nMatchClus = track->GetEMCALcluster(); | |
358 | if(nMatchClus > -1){ | |
359 | ||
360 | mom = track->P(); | |
361 | pt = track->Pt(); | |
362 | charge = track->Charge(); | |
363 | ||
364 | matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); | |
365 | ||
366 | if(matchedClus){ | |
367 | ||
368 | // matched cluster is EMCAL | |
369 | if(matchedClus->IsEMCAL()){ | |
370 | ||
371 | fClsE = matchedClus->E(); | |
372 | EovP = fClsE/mom; | |
373 | ||
374 | // fill used EMCAL variables here | |
375 | eop = EovP; // E/p | |
376 | showershape[0] = matchedClus->GetNCells(); // number of cells in cluster | |
377 | showershape[1] = matchedClus->GetM02(); // long axis | |
378 | showershape[2] = matchedClus->GetM20(); // short axis | |
379 | showershape[3] = matchedClus->GetDispersion(); // dispersion | |
1c9d11be | 380 | |
381 | // look for cached value first | |
382 | const AliDetectorPID *detPID=track->GetDetectorPID(); | |
383 | const EDetector detector=kEMCAL; | |
384 | ||
385 | if ( detPID && detPID->HasNumberOfSigmas(detector)){ | |
386 | return detPID->GetNumberOfSigmas(detector, type); | |
387 | } else if (fCachePID) { | |
388 | FillTrackDetectorPID(track, detector); | |
389 | detPID=track->GetDetectorPID(); | |
390 | return detPID->GetNumberOfSigmas(detector, type); | |
391 | } | |
392 | ||
393 | // NSigma value really meaningful only for electrons! | |
394 | return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); | |
6d0064aa | 395 | } |
396 | } | |
397 | } | |
398 | return -999; | |
399 | } | |
400 | ||
401 | ||
fd21ec8d | 402 | //______________________________________________________________________________ |
1c9d11be | 403 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
fd21ec8d | 404 | { |
405 | // | |
406 | // Compute PID response of 'detCode' | |
407 | // | |
408 | ||
409 | switch (detCode){ | |
410 | case kDetITS: return ComputeITSProbability(track, nSpecies, p); break; | |
411 | case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break; | |
fd21ec8d | 412 | case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break; |
00a38d07 | 413 | case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break; |
fd21ec8d | 414 | case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break; |
415 | case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break; | |
416 | case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break; | |
417 | default: return kDetNoSignal; | |
418 | } | |
419 | } | |
420 | ||
00a38d07 | 421 | //______________________________________________________________________________ |
422 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
423 | { | |
424 | // | |
425 | // Compute PID response of 'detCode' | |
426 | // | |
427 | ||
428 | return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p); | |
429 | } | |
430 | ||
fd21ec8d | 431 | //______________________________________________________________________________ |
432 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
433 | { | |
434 | // | |
435 | // Compute PID response for the ITS | |
436 | // | |
437 | ||
00a38d07 | 438 | // look for cached value first |
1c9d11be | 439 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
440 | const EDetector detector=kITS; | |
441 | ||
442 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ | |
443 | return detPID->GetRawProbability(detector, p, nSpecies); | |
444 | } else if (fCachePID) { | |
445 | FillTrackDetectorPID(track, detector); | |
446 | detPID=track->GetDetectorPID(); | |
447 | return detPID->GetRawProbability(detector, p, nSpecies); | |
fd21ec8d | 448 | } |
449 | ||
1c9d11be | 450 | return GetComputeITSProbability(track, nSpecies, p); |
fd21ec8d | 451 | } |
452 | //______________________________________________________________________________ | |
453 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
454 | { | |
455 | // | |
456 | // Compute PID response for the TPC | |
457 | // | |
00a38d07 | 458 | |
459 | // look for cached value first | |
1c9d11be | 460 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
461 | const EDetector detector=kTPC; | |
462 | ||
463 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ | |
464 | return detPID->GetRawProbability(detector, p, nSpecies); | |
465 | } else if (fCachePID) { | |
466 | FillTrackDetectorPID(track, detector); | |
467 | detPID=track->GetDetectorPID(); | |
468 | return detPID->GetRawProbability(detector, p, nSpecies); | |
00a38d07 | 469 | } |
470 | ||
1c9d11be | 471 | return GetComputeTPCProbability(track, nSpecies, p); |
fd21ec8d | 472 | } |
473 | //______________________________________________________________________________ | |
474 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
475 | { | |
476 | // | |
477 | // Compute PID response for the | |
478 | // | |
00a38d07 | 479 | |
1c9d11be | 480 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
481 | const EDetector detector=kTOF; | |
fd21ec8d | 482 | |
1c9d11be | 483 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ |
484 | return detPID->GetRawProbability(detector, p, nSpecies); | |
485 | } else if (fCachePID) { | |
486 | FillTrackDetectorPID(track, detector); | |
487 | detPID=track->GetDetectorPID(); | |
488 | return detPID->GetRawProbability(detector, p, nSpecies); | |
fd21ec8d | 489 | } |
1c9d11be | 490 | |
491 | return GetComputeTOFProbability(track, nSpecies, p); | |
fd21ec8d | 492 | } |
493 | //______________________________________________________________________________ | |
bd58d4b9 | 494 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const |
fd21ec8d | 495 | { |
496 | // | |
497 | // Compute PID response for the | |
1c9d11be | 498 | // |
00a38d07 | 499 | |
1c9d11be | 500 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
501 | const EDetector detector=kTRD; | |
ea235c90 | 502 | |
1c9d11be | 503 | // chacke only for the default method (1d at the moment) |
504 | if (PIDmethod!=AliTRDPIDResponse::kLQ1D) return GetComputeTRDProbability(track, nSpecies, p, PIDmethod); | |
505 | ||
506 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ | |
507 | return detPID->GetRawProbability(detector, p, nSpecies); | |
508 | } else if (fCachePID) { | |
509 | FillTrackDetectorPID(track, detector); | |
510 | detPID=track->GetDetectorPID(); | |
511 | return detPID->GetRawProbability(detector, p, nSpecies); | |
ea235c90 | 512 | } |
1c9d11be | 513 | |
514 | return GetComputeTRDProbability(track, nSpecies, p); | |
fd21ec8d | 515 | } |
516 | //______________________________________________________________________________ | |
e96b9916 | 517 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const |
fd21ec8d | 518 | { |
519 | // | |
520 | // Compute PID response for the EMCAL | |
521 | // | |
00a38d07 | 522 | |
1c9d11be | 523 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
524 | const EDetector detector=kEMCAL; | |
525 | ||
526 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ | |
527 | return detPID->GetRawProbability(detector, p, nSpecies); | |
528 | } else if (fCachePID) { | |
529 | FillTrackDetectorPID(track, detector); | |
530 | detPID=track->GetDetectorPID(); | |
531 | return detPID->GetRawProbability(detector, p, nSpecies); | |
00a38d07 | 532 | } |
1c9d11be | 533 | |
534 | return GetComputeEMCALProbability(track, nSpecies, p); | |
535 | } | |
536 | //______________________________________________________________________________ | |
537 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const | |
538 | { | |
539 | // | |
540 | // Compute PID response for the PHOS | |
541 | // | |
542 | ||
543 | // look for cached value first | |
544 | // if (track->GetDetectorPID()){ | |
545 | // return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies); | |
546 | // } | |
547 | ||
548 | // set flat distribution (no decision) | |
00a38d07 | 549 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; |
1c9d11be | 550 | return kDetNoSignal; |
551 | } | |
552 | //______________________________________________________________________________ | |
553 | AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
554 | { | |
555 | // | |
556 | // Compute PID response for the HMPID | |
557 | // | |
fd21ec8d | 558 | |
1c9d11be | 559 | const AliDetectorPID *detPID=track->GetDetectorPID(); |
560 | const EDetector detector=kHMPID; | |
e96b9916 | 561 | |
1c9d11be | 562 | if ( detPID && detPID->HasRawProbabilitiy(detector)){ |
563 | return detPID->GetRawProbability(detector, p, nSpecies); | |
564 | } else if (fCachePID) { | |
565 | FillTrackDetectorPID(track, detector); | |
566 | detPID=track->GetDetectorPID(); | |
567 | return detPID->GetRawProbability(detector, p, nSpecies); | |
e96b9916 | 568 | } |
00a38d07 | 569 | |
1c9d11be | 570 | return GetComputeHMPIDProbability(track, nSpecies, p); |
fd21ec8d | 571 | } |
572 | ||
4ec8e76d | 573 | //______________________________________________________________________________ |
00a38d07 | 574 | void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run) |
4ec8e76d | 575 | { |
576 | // | |
577 | // Apply settings for the current event | |
578 | // | |
579 | fRecoPass=pass; | |
e96b9916 | 580 | |
78cbd205 | 581 | |
644666df | 582 | fCurrentEvent=NULL; |
4ec8e76d | 583 | if (!event) return; |
e96b9916 | 584 | fCurrentEvent=event; |
00a38d07 | 585 | if (run>0) fRun=run; |
586 | else fRun=event->GetRunNumber(); | |
4ec8e76d | 587 | |
588 | if (fRun!=fOldRun){ | |
589 | ExecNewRun(); | |
590 | fOldRun=fRun; | |
591 | } | |
592 | ||
593 | //TPC resolution parametrisation PbPb | |
594 | if ( fResolutionCorrection ){ | |
595 | Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event)); | |
596 | fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04); | |
597 | } | |
598 | ||
599 | //TOF resolution | |
b79db598 | 600 | SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod()); |
601 | ||
80f28562 | 602 | |
603 | // Get and set centrality | |
604 | AliCentrality *centrality = event->GetCentrality(); | |
605 | if(centrality){ | |
606 | fCurrCentrality = centrality->GetCentralityPercentile("V0M"); | |
607 | } | |
608 | else{ | |
609 | fCurrCentrality = -1; | |
610 | } | |
4ec8e76d | 611 | } |
612 | ||
613 | //______________________________________________________________________________ | |
614 | void AliPIDResponse::ExecNewRun() | |
615 | { | |
616 | // | |
617 | // Things to Execute upon a new run | |
618 | // | |
619 | SetRecoInfo(); | |
620 | ||
621 | SetITSParametrisation(); | |
622 | ||
623 | SetTPCPidResponseMaster(); | |
624 | SetTPCParametrisation(); | |
53d016dc | 625 | |
626 | SetTRDPidResponseMaster(); | |
627 | InitializeTRDResponse(); | |
b2138b40 | 628 | |
629 | SetEMCALPidResponseMaster(); | |
630 | InitializeEMCALResponse(); | |
4ec8e76d | 631 | |
b79db598 | 632 | SetTOFPidResponseMaster(); |
633 | InitializeTOFResponse(); | |
644666df | 634 | |
635 | if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField()); | |
4ec8e76d | 636 | } |
637 | ||
1c9d11be | 638 | //______________________________________________________________________________ |
4ec8e76d | 639 | Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event) |
640 | { | |
641 | // | |
642 | // Get TPC multiplicity in bins of 150 | |
643 | // | |
644 | ||
645 | const AliVVertex* vertexTPC = event->GetPrimaryVertex(); | |
646 | Double_t tpcMulti=0.; | |
647 | if(vertexTPC){ | |
648 | Double_t vertexContribTPC=vertexTPC->GetNContributors(); | |
649 | tpcMulti=vertexContribTPC/150.; | |
650 | if (tpcMulti>20.) tpcMulti=20.; | |
651 | } | |
652 | ||
653 | return tpcMulti; | |
654 | } | |
655 | ||
656 | //______________________________________________________________________________ | |
657 | void AliPIDResponse::SetRecoInfo() | |
658 | { | |
659 | // | |
660 | // Set reconstruction information | |
661 | // | |
662 | ||
663 | //reset information | |
664 | fLHCperiod=""; | |
665 | fMCperiodTPC=""; | |
666 | ||
667 | fBeamType=""; | |
668 | ||
669 | fBeamType="PP"; | |
670 | ||
1436d6bb | 671 | TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*"); |
a78fd045 | 672 | TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*"); |
1436d6bb | 673 | |
4ec8e76d | 674 | //find the period by run number (UGLY, but not stored in ESD and AOD... ) |
675 | if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; } | |
676 | else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; } | |
677 | else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; } | |
99e9d5ec | 678 | else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; } |
4ec8e76d | 679 | else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; } |
680 | else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; } | |
ea235c90 | 681 | else if (fRun>=136851&&fRun<=139517) { |
682 | fLHCperiod="LHC10H"; | |
683 | fMCperiodTPC="LHC10H8"; | |
684 | if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10"; | |
685 | fBeamType="PBPB"; | |
686 | } | |
3077a03d | 687 | else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; } |
2e97a211 | 688 | //TODO: periods 11B, 11C are not yet treated assume 11d for the moment |
689 | else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } | |
690 | else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } | |
00a38d07 | 691 | // also for 11e,f use 11d |
692 | else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } | |
693 | else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; } | |
694 | ||
8af51a65 | 695 | else if (fRun>=166529 && fRun<=170718) { |
3077a03d | 696 | fLHCperiod="LHC11H"; |
697 | fMCperiodTPC="LHC11A10"; | |
698 | fBeamType="PBPB"; | |
a78fd045 | 699 | if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17"; |
3077a03d | 700 | } |
8af51a65 | 701 | if (fRun>=170719 && fRun<=177311) { fLHCperiod="LHC12A"; fBeamType="PP"; /*fMCperiodTPC="";*/ } |
702 | // for the moment use LHC12b parameters up to LHC12e | |
703 | if (fRun>=177312 /*&& fRun<=179356*/) { fLHCperiod="LHC12B"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
704 | // if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
705 | // if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
706 | // if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
707 | ||
708 | // if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
709 | // if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; /*fMCperiodTPC="";*/ } | |
710 | // if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ } | |
711 | // for the moment use 12g parametrisation for all full gain runs (LHC12f+) | |
712 | if (fRun >= 186636 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; /*fMCperiodTPC="";*/ } | |
80ab5635 | 713 | |
714 | //exception new pp MC productions from 2011 | |
715 | if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2"; | |
4a527e08 | 716 | // exception for 11f1 |
00a38d07 | 717 | if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1"; |
4ec8e76d | 718 | } |
719 | ||
720 | //______________________________________________________________________________ | |
721 | void AliPIDResponse::SetITSParametrisation() | |
722 | { | |
723 | // | |
724 | // Set the ITS parametrisation | |
725 | // | |
726 | } | |
727 | ||
728 | //______________________________________________________________________________ | |
729 | void AliPIDResponse::SetTPCPidResponseMaster() | |
730 | { | |
731 | // | |
732 | // Load the TPC pid response functions from the OADB | |
644666df | 733 | // Load the TPC voltage maps from OADB |
4ec8e76d | 734 | // |
09b50a42 | 735 | //don't load twice for the moment |
736 | if (fArrPidResponseMaster) return; | |
737 | ||
738 | ||
4ec8e76d | 739 | //reset the PID response functions |
740 | delete fArrPidResponseMaster; | |
644666df | 741 | fArrPidResponseMaster=NULL; |
4ec8e76d | 742 | |
743 | TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data())); | |
644666df | 744 | TFile *f=NULL; |
00a38d07 | 745 | if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse; |
4ec8e76d | 746 | |
644666df | 747 | TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data())); |
748 | f=TFile::Open(fileNamePIDresponse.Data()); | |
ea235c90 | 749 | if (f && f->IsOpen() && !f->IsZombie()){ |
750 | fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse")); | |
4ec8e76d | 751 | } |
ea235c90 | 752 | delete f; |
644666df | 753 | |
754 | TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data())); | |
755 | f=TFile::Open(fileNameVoltageMaps.Data()); | |
756 | if (f && f->IsOpen() && !f->IsZombie()){ | |
757 | fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings")); | |
758 | } | |
759 | delete f; | |
4ec8e76d | 760 | |
761 | if (!fArrPidResponseMaster){ | |
644666df | 762 | AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data())); |
4ec8e76d | 763 | return; |
764 | } | |
765 | fArrPidResponseMaster->SetOwner(); | |
644666df | 766 | |
767 | if (!fOADBvoltageMaps) | |
768 | { | |
769 | AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data())); | |
770 | } | |
771 | fArrPidResponseMaster->SetOwner(); | |
4ec8e76d | 772 | } |
773 | ||
774 | //______________________________________________________________________________ | |
775 | void AliPIDResponse::SetTPCParametrisation() | |
776 | { | |
777 | // | |
778 | // Change BB parametrisation for current run | |
779 | // | |
780 | ||
781 | if (fLHCperiod.IsNull()) { | |
782 | AliFatal("No period set, not changing parametrisation"); | |
783 | return; | |
784 | } | |
785 | ||
786 | // | |
787 | // Set default parametrisations for data and MC | |
788 | // | |
789 | ||
790 | //data type | |
791 | TString datatype="DATA"; | |
792 | //in case of mc fRecoPass is per default 1 | |
793 | if (fIsMC) { | |
539a5a59 | 794 | if(!fTuneMConData) datatype="MC"; |
795 | fRecoPass=1; | |
4ec8e76d | 796 | } |
797 | ||
798 | // | |
799 | //reset old splines | |
800 | // | |
644666df | 801 | fTPCResponse.ResetSplines(); |
4a527e08 | 802 | |
803 | // period | |
804 | TString period=fLHCperiod; | |
539a5a59 | 805 | if (fIsMC && !fTuneMConData) period=fMCperiodTPC; |
4a527e08 | 806 | |
807 | AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); | |
808 | Bool_t found=kFALSE; | |
4ec8e76d | 809 | // |
810 | //set the new PID splines | |
811 | // | |
4ec8e76d | 812 | if (fArrPidResponseMaster){ |
539a5a59 | 813 | Int_t recopass = fRecoPass; |
814 | if(fTuneMConData) recopass = fRecoPassUser; | |
4ec8e76d | 815 | //for MC don't use period information |
644666df | 816 | //if (fIsMC) period="[A-Z0-9]*"; |
4ec8e76d | 817 | //for MC use MC period information |
644666df | 818 | //pattern for the default entry (valid for all particles) |
de678885 | 819 | TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data())); |
644666df | 820 | |
821 | //find particle id ang gain scenario | |
822 | for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++) | |
823 | { | |
824 | TObject *grAll=NULL; | |
825 | TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario); | |
826 | gainScenario.ToUpper(); | |
827 | //loop over entries and filter them | |
828 | for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp) | |
829 | { | |
830 | TObject *responseFunction=fArrPidResponseMaster->At(iresp); | |
831 | if (responseFunction==NULL) continue; | |
832 | TString responseName=responseFunction->GetName(); | |
833 | ||
834 | if (!reg.MatchB(responseName)) continue; | |
835 | ||
836 | TObjArray *arr=reg.MatchS(responseName); if (!arr) continue; | |
837 | TObject* tmp=NULL; | |
838 | tmp=arr->At(1); if (!tmp) continue; | |
839 | TString particleName=tmp->GetName(); | |
840 | tmp=arr->At(3); if (!tmp) continue; | |
841 | TString gainScenarioName=tmp->GetName(); | |
842 | delete arr; | |
843 | if (particleName.IsNull()) continue; | |
844 | if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction; | |
845 | else | |
846 | { | |
847 | for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec) | |
848 | { | |
849 | TString particle=AliPID::ParticleName(ispec); | |
850 | particle.ToUpper(); | |
851 | //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl; | |
852 | if ( particle == particleName && gainScenario == gainScenarioName ) | |
853 | { | |
854 | fTPCResponse.SetResponseFunction( responseFunction, | |
855 | (AliPID::EParticleType)ispec, | |
856 | (AliTPCPIDResponse::ETPCgainScenario)igainScenario ); | |
857 | fTPCResponse.SetUseDatabase(kTRUE); | |
858 | AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,responseFunction->GetName())); | |
859 | found=kTRUE; | |
860 | // overwrite default with proton spline (for light nuclei) | |
861 | if (ispec==AliPID::kProton) grAll=responseFunction; | |
862 | break; | |
863 | } | |
4ec8e76d | 864 | } |
865 | } | |
866 | } | |
644666df | 867 | if (grAll) |
868 | { | |
869 | for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec) | |
870 | { | |
871 | if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec, | |
872 | (AliTPCPIDResponse::ETPCgainScenario)igainScenario)) | |
873 | { | |
874 | fTPCResponse.SetResponseFunction( grAll, | |
875 | (AliPID::EParticleType)ispec, | |
876 | (AliTPCPIDResponse::ETPCgainScenario)igainScenario ); | |
877 | fTPCResponse.SetUseDatabase(kTRUE); | |
878 | AliInfo(Form("Adding graph: %d %d - %s",ispec,igainScenario,grAll->GetName())); | |
879 | found=kTRUE; | |
880 | } | |
4ec8e76d | 881 | } |
882 | } | |
883 | } | |
884 | } | |
644666df | 885 | else AliInfo("no fArrPidResponseMaster"); |
4a527e08 | 886 | |
887 | if (!found){ | |
888 | AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); | |
889 | } | |
644666df | 890 | |
4ec8e76d | 891 | // |
892 | // Setup resolution parametrisation | |
893 | // | |
894 | ||
895 | //default | |
896 | fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04); | |
897 | ||
898 | if (fRun>=122195){ | |
899 | fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02); | |
900 | } | |
8af51a65 | 901 | |
902 | if (fRun>=186636){ | |
903 | // if (fRun>=188356){ | |
723c4874 | 904 | fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05); |
905 | } | |
906 | ||
23425eb2 | 907 | if (fArrPidResponseMaster) |
4ec8e76d | 908 | fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data())); |
909 | ||
910 | if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName())); | |
644666df | 911 | |
912 | //read in the voltage map | |
913 | TVectorF* gsm = dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun)); | |
914 | if (gsm) | |
915 | { | |
916 | fTPCResponse.SetVoltageMap(*gsm); | |
917 | TString vals; | |
918 | AliInfo(Form("Reading the voltage map for run %d\n",fRun)); | |
919 | vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);} | |
920 | AliInfo(vals.Data()); | |
921 | vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);} | |
922 | AliInfo(vals.Data()); | |
923 | vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);} | |
924 | AliInfo(vals.Data()); | |
925 | vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);} | |
926 | AliInfo(vals.Data()); | |
927 | } | |
928 | else AliInfo("no voltage map, ideal default assumed"); | |
4ec8e76d | 929 | } |
930 | ||
ea235c90 | 931 | //______________________________________________________________________________ |
932 | void AliPIDResponse::SetTRDPidResponseMaster() | |
933 | { | |
934 | // | |
935 | // Load the TRD pid params and references from the OADB | |
936 | // | |
db0e2c5f | 937 | if(fTRDPIDResponseObject) return; |
53d016dc | 938 | AliOADBContainer contParams("contParams"); |
939 | ||
db0e2c5f | 940 | Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject"); |
941 | if(statusResponse){ | |
942 | AliError("Failed initializing PID Response Object from OADB"); | |
59a8e853 | 943 | } else { |
db0e2c5f | 944 | AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data())); |
945 | fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun)); | |
946 | if(!fTRDPIDResponseObject){ | |
947 | AliError(Form("TRD Response not found in run %d", fRun)); | |
59a8e853 | 948 | } |
949 | } | |
ea235c90 | 950 | } |
951 | ||
952 | //______________________________________________________________________________ | |
953 | void AliPIDResponse::InitializeTRDResponse(){ | |
954 | // | |
955 | // Set PID Params and references to the TRD PID response | |
956 | // | |
db0e2c5f | 957 | fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject); |
f2762b1c | 958 | } |
959 | ||
bd58d4b9 | 960 | //______________________________________________________________________________ |
961 | void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{ | |
962 | ||
963 | if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){ | |
964 | // backward compatibility for setting with 8 slices | |
965 | TRDslicesForPID[0] = 0; | |
966 | TRDslicesForPID[1] = 7; | |
f2762b1c | 967 | } |
bd58d4b9 | 968 | else{ |
969 | if(method==AliTRDPIDResponse::kLQ1D){ | |
970 | TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx | |
971 | TRDslicesForPID[1] = 0; | |
972 | } | |
973 | if(method==AliTRDPIDResponse::kLQ2D){ | |
974 | TRDslicesForPID[0] = 1; | |
975 | TRDslicesForPID[1] = 7; | |
976 | } | |
db0e2c5f | 977 | } |
bd58d4b9 | 978 | AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1])); |
ea235c90 | 979 | } |
980 | ||
b79db598 | 981 | //______________________________________________________________________________ |
982 | void AliPIDResponse::SetTOFPidResponseMaster() | |
983 | { | |
984 | // | |
985 | // Load the TOF pid params from the OADB | |
986 | // | |
00a38d07 | 987 | |
988 | if (fTOFPIDParams) delete fTOFPIDParams; | |
644666df | 989 | fTOFPIDParams=NULL; |
00a38d07 | 990 | |
b79db598 | 991 | TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data())); |
00a38d07 | 992 | if (oadbf && oadbf->IsOpen()) { |
b79db598 | 993 | AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data())); |
994 | AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb"); | |
00a38d07 | 995 | if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams")); |
b79db598 | 996 | oadbf->Close(); |
997 | delete oadbc; | |
b79db598 | 998 | } |
999 | delete oadbf; | |
1000 | ||
00a38d07 | 1001 | if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved"); |
1002 | } | |
b79db598 | 1003 | |
1004 | //______________________________________________________________________________ | |
1005 | void AliPIDResponse::InitializeTOFResponse(){ | |
1006 | // | |
1007 | // Set PID Params to the TOF PID response | |
00a38d07 | 1008 | // |
1009 | ||
1010 | AliInfo("TOF PID Params loaded from OADB"); | |
1011 | AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution())); | |
1012 | AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod())); | |
1013 | AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f", | |
1014 | fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3))); | |
1015 | ||
b79db598 | 1016 | for (Int_t i=0;i<4;i++) { |
1017 | fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i)); | |
1018 | } | |
1019 | fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution()); | |
1020 | ||
78cbd205 | 1021 | AliInfo("TZERO resolution loaded from ESDrun/AODheader"); |
1022 | Float_t t0Spread[4]; | |
1023 | for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i); | |
1024 | AliInfo(Form(" TZERO spreads from data: (A+C)/2 %f A %f C %f (A'-C')/2: %f",t0Spread[0],t0Spread[1],t0Spread[2],t0Spread[3])); | |
1025 | Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3]; | |
1026 | Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3]; | |
1027 | if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) { | |
1028 | fResT0AC=t0Spread[3]; | |
1029 | fResT0A=TMath::Sqrt(a); | |
1030 | fResT0C=TMath::Sqrt(c); | |
1031 | } else { | |
1032 | AliInfo(" TZERO spreads not present or inconsistent, loading default"); | |
1033 | fResT0A=75.; | |
1034 | fResT0C=65.; | |
1035 | fResT0AC=55.; | |
1036 | } | |
1037 | AliInfo(Form(" TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC)); | |
1038 | ||
b79db598 | 1039 | } |
1040 | ||
1041 | ||
1c9d11be | 1042 | //______________________________________________________________________________ |
bd58d4b9 | 1043 | Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const { |
ea235c90 | 1044 | // |
1045 | // Check whether track is identified as electron under a given electron efficiency hypothesis | |
bd58d4b9 | 1046 | // |
1047 | ||
ea235c90 | 1048 | Double_t probs[AliPID::kSPECIES]; |
bd58d4b9 | 1049 | ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs,PIDmethod); |
ea235c90 | 1050 | |
99e9d5ec | 1051 | Int_t ntracklets = vtrack->GetTRDntrackletsPID(); |
1052 | // Take mean of the TRD momenta in the given tracklets | |
1053 | Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes]; | |
1054 | Int_t nmomenta = 0; | |
ea235c90 | 1055 | for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){ |
1056 | if(vtrack->GetTRDmomentum(iPl) > 0.){ | |
99e9d5ec | 1057 | trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl); |
ea235c90 | 1058 | } |
1059 | } | |
99e9d5ec | 1060 | p = TMath::Mean(nmomenta, trdmomenta); |
ea235c90 | 1061 | |
bd58d4b9 | 1062 | return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod); |
ea235c90 | 1063 | } |
1064 | ||
b2138b40 | 1065 | //______________________________________________________________________________ |
1066 | void AliPIDResponse::SetEMCALPidResponseMaster() | |
1067 | { | |
1068 | // | |
1069 | // Load the EMCAL pid response functions from the OADB | |
1070 | // | |
1071 | TObjArray* fEMCALPIDParamsRun = NULL; | |
1072 | TObjArray* fEMCALPIDParamsPass = NULL; | |
1073 | ||
1074 | if(fEMCALPIDParams) return; | |
1075 | AliOADBContainer contParams("contParams"); | |
1076 | ||
1077 | Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams"); | |
1078 | if(statusPars){ | |
1079 | AliError("Failed initializing PID Params from OADB"); | |
1080 | } | |
1081 | else { | |
1082 | AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data())); | |
1083 | ||
1084 | fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun)); | |
1085 | if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass))); | |
1086 | if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles"))); | |
1087 | ||
1088 | if(!fEMCALPIDParams){ | |
f8d39067 | 1089 | AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass)); |
1f631618 | 1090 | AliInfo("Will take the standard LHC11d instead ..."); |
b2138b40 | 1091 | |
1f631618 | 1092 | fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477)); |
1093 | if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1))); | |
b2138b40 | 1094 | if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles"))); |
1095 | ||
1096 | if(!fEMCALPIDParams){ | |
1f631618 | 1097 | AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data())); |
b2138b40 | 1098 | } |
1099 | } | |
1100 | } | |
1101 | } | |
1102 | ||
1103 | //______________________________________________________________________________ | |
1104 | void AliPIDResponse::InitializeEMCALResponse(){ | |
1105 | // | |
1106 | // Set PID Params to the EMCAL PID response | |
1107 | // | |
1108 | fEMCALResponse.SetPIDParams(fEMCALPIDParams); | |
1109 | ||
1110 | } | |
00a38d07 | 1111 | |
1c9d11be | 1112 | //______________________________________________________________________________ |
1113 | void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const | |
00a38d07 | 1114 | { |
1115 | // | |
1116 | // create detector PID information and setup the transient pointer in the track | |
1117 | // | |
1c9d11be | 1118 | |
1119 | // check if detector number is inside accepted range | |
1120 | if (detector == kNdetectors) return; | |
1121 | ||
1122 | // get detector pid | |
1123 | AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID()); | |
1124 | if (!detPID) { | |
1125 | detPID=new AliDetectorPID; | |
1126 | (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID); | |
1127 | } | |
1128 | ||
1129 | //check if values exist | |
1130 | if (detPID->HasRawProbabilitiy(detector) && detPID->HasNumberOfSigmas(detector)) return; | |
00a38d07 | 1131 | |
1132 | //TODO: which particles to include? See also the loops below... | |
1133 | Double_t values[AliPID::kSPECIESC]={0}; | |
1c9d11be | 1134 | |
1135 | //nsigmas | |
1136 | for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart) | |
1137 | values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart); | |
1138 | detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC); | |
1139 | ||
1140 | //probabilities | |
1141 | EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values); | |
1142 | detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status); | |
1143 | } | |
1144 | ||
1145 | //______________________________________________________________________________ | |
1146 | void AliPIDResponse::FillTrackDetectorPID() | |
1147 | { | |
1148 | // | |
1149 | // create detector PID information and setup the transient pointer in the track | |
1150 | // | |
1151 | ||
1152 | if (!fCurrentEvent) return; | |
00a38d07 | 1153 | |
1154 | for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){ | |
1155 | AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack)); | |
1156 | if (!track) continue; | |
1157 | ||
00a38d07 | 1158 | for (Int_t idet=0; idet<kNdetectors; ++idet){ |
1c9d11be | 1159 | FillTrackDetectorPID(track, (EDetector)idet); |
00a38d07 | 1160 | } |
00a38d07 | 1161 | } |
1162 | } | |
1163 | ||
1c9d11be | 1164 | //______________________________________________________________________________ |
5f8db5fe | 1165 | void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){ |
1166 | // | |
1167 | // Set TOF response function | |
1168 | // Input option for event_time used | |
1169 | // | |
1170 | ||
1171 | Float_t t0spread = 0.; //vevent->GetEventTimeSpread(); | |
1172 | if(t0spread < 10) t0spread = 80; | |
1173 | ||
1174 | // T0 from TOF algorithm | |
1175 | ||
1176 | Bool_t flagT0TOF=kFALSE; | |
1177 | Bool_t flagT0T0=kFALSE; | |
1178 | Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()]; | |
1179 | Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()]; | |
1180 | Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()]; | |
1181 | ||
1182 | // T0-TOF arrays | |
1183 | Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()]; | |
1184 | Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()]; | |
1185 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1186 | estimatedT0event[i]=0.0; | |
1187 | estimatedT0resolution[i]=0.0; | |
1188 | startTimeMask[i] = 0; | |
1189 | } | |
1190 | ||
78cbd205 | 1191 | Float_t resT0A=fResT0A; |
1192 | Float_t resT0C=fResT0C; | |
1193 | Float_t resT0AC=fResT0AC; | |
5f8db5fe | 1194 | if(vevent->GetT0TOF()){ // check if T0 detector information is available |
1195 | flagT0T0=kTRUE; | |
1196 | } | |
1197 | ||
1198 | ||
1199 | AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader(); | |
1200 | ||
1201 | if (tofHeader) { // read global info and T0-TOF | |
1202 | fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution()); | |
1203 | t0spread = tofHeader->GetT0spread(); // read t0 sprad | |
1204 | if(t0spread < 10) t0spread = 80; | |
1205 | ||
1206 | flagT0TOF=kTRUE; | |
1207 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value | |
1208 | startTime[i]=tofHeader->GetDefaultEventTimeVal(); | |
1209 | startTimeRes[i]=tofHeader->GetDefaultEventTimeRes(); | |
1210 | if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread; | |
1211 | } | |
1212 | ||
1213 | TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues(); | |
1214 | TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues(); | |
1215 | TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes(); | |
1216 | for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins | |
1217 | Int_t icurrent = (Int_t)ibin->GetAt(j); | |
1218 | startTime[icurrent]=t0Bin->GetAt(j); | |
1219 | startTimeRes[icurrent]=t0ResBin->GetAt(j); | |
1220 | if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread; | |
1221 | } | |
1222 | } | |
1223 | ||
1224 | // for cut of 3 sigma on t0 spread | |
1225 | Float_t t0cut = 3 * t0spread; | |
1226 | if(t0cut < 500) t0cut = 500; | |
1227 | ||
1228 | if(option == kFILL_T0){ // T0-FILL is used | |
1229 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1230 | estimatedT0event[i]=0.0; | |
1231 | estimatedT0resolution[i]=t0spread; | |
1232 | } | |
1233 | fTOFResponse.SetT0event(estimatedT0event); | |
1234 | fTOFResponse.SetT0resolution(estimatedT0resolution); | |
1235 | } | |
1236 | ||
1237 | if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD | |
1238 | if(flagT0TOF){ | |
1239 | fTOFResponse.SetT0event(startTime); | |
1240 | fTOFResponse.SetT0resolution(startTimeRes); | |
1241 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1242 | if(startTimeRes[i]<t0spread) startTimeMask[i]=1; | |
1243 | fTOFResponse.SetT0binMask(i,startTimeMask[i]); | |
1244 | } | |
1245 | } | |
1246 | else{ | |
1247 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1248 | estimatedT0event[i]=0.0; | |
1249 | estimatedT0resolution[i]=t0spread; | |
1250 | fTOFResponse.SetT0binMask(i,startTimeMask[i]); | |
1251 | } | |
1252 | fTOFResponse.SetT0event(estimatedT0event); | |
1253 | fTOFResponse.SetT0resolution(estimatedT0resolution); | |
1254 | } | |
1255 | } | |
1256 | else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD | |
1257 | Float_t t0AC=-10000; | |
1258 | Float_t t0A=-10000; | |
1259 | Float_t t0C=-10000; | |
1260 | if(flagT0T0){ | |
1261 | t0AC= vevent->GetT0TOF()[0]; | |
1262 | t0A= vevent->GetT0TOF()[1]; | |
1263 | t0C= vevent->GetT0TOF()[2]; | |
1264 | } | |
1265 | ||
1266 | Float_t t0t0Best = 0; | |
1267 | Float_t t0t0BestRes = 9999; | |
1268 | Int_t t0used=0; | |
1269 | if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){ | |
1270 | t0t0Best = t0AC; | |
1271 | t0t0BestRes = resT0AC; | |
1272 | t0used=6; | |
1273 | } | |
1274 | else if(TMath::Abs(t0C) < t0cut){ | |
1275 | t0t0Best = t0C; | |
1276 | t0t0BestRes = resT0C; | |
1277 | t0used=4; | |
1278 | } | |
1279 | else if(TMath::Abs(t0A) < t0cut){ | |
1280 | t0t0Best = t0A; | |
1281 | t0t0BestRes = resT0A; | |
1282 | t0used=2; | |
1283 | } | |
1284 | ||
1285 | if(flagT0TOF){ // if T0-TOF info is available | |
1286 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1287 | if(t0t0BestRes < 999){ | |
1288 | if(startTimeRes[i] < t0spread){ | |
1289 | Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes; | |
1290 | Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes; | |
1291 | estimatedT0event[i]=t0best / wtot; | |
1292 | estimatedT0resolution[i]=1./TMath::Sqrt(wtot); | |
1293 | startTimeMask[i] = t0used+1; | |
1294 | } | |
1295 | else { | |
1296 | estimatedT0event[i]=t0t0Best; | |
1297 | estimatedT0resolution[i]=t0t0BestRes; | |
1298 | startTimeMask[i] = t0used; | |
1299 | } | |
1300 | } | |
1301 | else{ | |
1302 | estimatedT0event[i]=startTime[i]; | |
1303 | estimatedT0resolution[i]=startTimeRes[i]; | |
1304 | if(startTimeRes[i]<t0spread) startTimeMask[i]=1; | |
1305 | } | |
1306 | fTOFResponse.SetT0binMask(i,startTimeMask[i]); | |
1307 | } | |
1308 | fTOFResponse.SetT0event(estimatedT0event); | |
1309 | fTOFResponse.SetT0resolution(estimatedT0resolution); | |
1310 | } | |
1311 | else{ // if no T0-TOF info is available | |
1312 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1313 | fTOFResponse.SetT0binMask(i,t0used); | |
1314 | if(t0t0BestRes < 999){ | |
1315 | estimatedT0event[i]=t0t0Best; | |
1316 | estimatedT0resolution[i]=t0t0BestRes; | |
1317 | } | |
1318 | else{ | |
1319 | estimatedT0event[i]=0.0; | |
1320 | estimatedT0resolution[i]=t0spread; | |
1321 | } | |
1322 | } | |
1323 | fTOFResponse.SetT0event(estimatedT0event); | |
1324 | fTOFResponse.SetT0resolution(estimatedT0resolution); | |
1325 | } | |
1326 | } | |
1327 | ||
1328 | else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise) | |
1329 | Float_t t0AC=-10000; | |
1330 | Float_t t0A=-10000; | |
1331 | Float_t t0C=-10000; | |
1332 | if(flagT0T0){ | |
1333 | t0AC= vevent->GetT0TOF()[0]; | |
1334 | t0A= vevent->GetT0TOF()[1]; | |
1335 | t0C= vevent->GetT0TOF()[2]; | |
1336 | } | |
1337 | ||
1338 | if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){ | |
1339 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1340 | estimatedT0event[i]=t0AC; | |
1341 | estimatedT0resolution[i]=resT0AC; | |
1342 | fTOFResponse.SetT0binMask(i,6); | |
1343 | } | |
1344 | } | |
1345 | else if(TMath::Abs(t0C) < t0cut){ | |
1346 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1347 | estimatedT0event[i]=t0C; | |
1348 | estimatedT0resolution[i]=resT0C; | |
1349 | fTOFResponse.SetT0binMask(i,4); | |
1350 | } | |
1351 | } | |
1352 | else if(TMath::Abs(t0A) < t0cut){ | |
1353 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1354 | estimatedT0event[i]=t0A; | |
1355 | estimatedT0resolution[i]=resT0A; | |
1356 | fTOFResponse.SetT0binMask(i,2); | |
1357 | } | |
1358 | } | |
1359 | else{ | |
1360 | for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ | |
1361 | estimatedT0event[i]=0.0; | |
1362 | estimatedT0resolution[i]=t0spread; | |
1363 | fTOFResponse.SetT0binMask(i,0); | |
1364 | } | |
1365 | } | |
1366 | fTOFResponse.SetT0event(estimatedT0event); | |
1367 | fTOFResponse.SetT0resolution(estimatedT0resolution); | |
1368 | } | |
1369 | delete [] startTime; | |
1370 | delete [] startTimeRes; | |
1371 | delete [] startTimeMask; | |
1372 | delete [] estimatedT0event; | |
1373 | delete [] estimatedT0resolution; | |
1374 | } | |
1c9d11be | 1375 | |
1376 | //______________________________________________________________________________ | |
1377 | // private non cached versions of the PID calculation | |
1378 | // | |
1379 | ||
1380 | ||
1381 | //______________________________________________________________________________ | |
1382 | Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const | |
1383 | { | |
1384 | // | |
1385 | // NumberOfSigmas for 'detCode' | |
1386 | // | |
1387 | ||
1388 | switch (detCode){ | |
1389 | case kITS: return GetNumberOfSigmasITS(track, type); break; | |
1390 | case kTPC: return GetNumberOfSigmasTPC(track, type); break; | |
1391 | case kTOF: return GetNumberOfSigmasTOF(track, type); break; | |
1392 | case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break; | |
1393 | default: return -999.; | |
1394 | } | |
1395 | ||
1396 | } | |
1397 | ||
1398 | ||
1399 | ||
1400 | //______________________________________________________________________________ | |
1401 | Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
1402 | { | |
1403 | // | |
1404 | // Calculate the number of sigmas in the ITS | |
1405 | // | |
1406 | ||
1407 | AliVTrack *track=(AliVTrack*)vtrack; | |
1408 | ||
1409 | Float_t dEdx=track->GetITSsignal(); | |
1410 | if (dEdx<=0) return -999.; | |
1411 | ||
1412 | UChar_t clumap=track->GetITSClusterMap(); | |
1413 | Int_t nPointsForPid=0; | |
1414 | for(Int_t i=2; i<6; i++){ | |
1415 | if(clumap&(1<<i)) ++nPointsForPid; | |
1416 | } | |
1417 | Float_t mom=track->P(); | |
1418 | ||
1419 | //check for ITS standalone tracks | |
1420 | Bool_t isSA=kTRUE; | |
1421 | if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE; | |
1422 | ||
1423 | //TODO: in case of the electron, use the SA parametrisation, | |
1424 | // this needs to be changed if ITS provides a parametrisation | |
1425 | // for electrons also for ITS+TPC tracks | |
1426 | return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron)); | |
1427 | } | |
1428 | ||
1429 | //______________________________________________________________________________ | |
1430 | Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
1431 | { | |
1432 | // | |
1433 | // Calculate the number of sigmas in the TPC | |
1434 | // | |
1435 | ||
1436 | AliVTrack *track=(AliVTrack*)vtrack; | |
1437 | ||
1438 | Double_t mom = track->GetTPCmomentum(); | |
1439 | Double_t sig = track->GetTPCsignal(); | |
1440 | UInt_t sigN = track->GetTPCsignalN(); | |
1441 | ||
1442 | Double_t nSigma = -999.; | |
1443 | if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type); | |
1444 | ||
1445 | return nSigma; | |
1446 | } | |
1447 | ||
1448 | //______________________________________________________________________________ | |
1449 | Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const | |
1450 | { | |
1451 | // | |
1452 | // Calculate the number of sigmas in the EMCAL | |
1453 | // | |
1454 | ||
1455 | AliVTrack *track=(AliVTrack*)vtrack; | |
1456 | ||
1457 | AliVCluster *matchedClus = NULL; | |
1458 | ||
1459 | Double_t mom = -1.; | |
1460 | Double_t pt = -1.; | |
1461 | Double_t EovP = -1.; | |
1462 | Double_t fClsE = -1.; | |
1463 | ||
1464 | Int_t nMatchClus = -1; | |
1465 | Int_t charge = 0; | |
1466 | ||
1467 | // Track matching | |
1468 | nMatchClus = track->GetEMCALcluster(); | |
1469 | if(nMatchClus > -1){ | |
1470 | ||
1471 | mom = track->P(); | |
1472 | pt = track->Pt(); | |
1473 | charge = track->Charge(); | |
1474 | ||
1475 | matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); | |
1476 | ||
1477 | if(matchedClus){ | |
1478 | ||
1479 | // matched cluster is EMCAL | |
1480 | if(matchedClus->IsEMCAL()){ | |
1481 | ||
1482 | fClsE = matchedClus->E(); | |
1483 | EovP = fClsE/mom; | |
1484 | ||
1485 | ||
1486 | // NSigma value really meaningful only for electrons! | |
1487 | return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); | |
1488 | } | |
1489 | } | |
1490 | } | |
1491 | ||
1492 | return -999; | |
1493 | ||
1494 | } | |
1495 | ||
1496 | ||
1497 | //______________________________________________________________________________ | |
1498 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1499 | { | |
1500 | // | |
1501 | // Compute PID response of 'detCode' | |
1502 | // | |
1503 | ||
1504 | switch (detCode){ | |
1505 | case kITS: return GetComputeITSProbability(track, nSpecies, p); break; | |
1506 | case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break; | |
1507 | case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break; | |
1508 | case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break; | |
1509 | case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break; | |
1510 | case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break; | |
1511 | case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break; | |
1512 | default: return kDetNoSignal; | |
1513 | } | |
1514 | } | |
1515 | ||
1516 | //______________________________________________________________________________ | |
1517 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1518 | { | |
1519 | // | |
1520 | // Compute PID response for the ITS | |
1521 | // | |
1522 | ||
1523 | // look for cached value first | |
1524 | // only the non SA tracks are cached | |
1525 | if (track->GetDetectorPID()){ | |
1526 | return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies); | |
1527 | } | |
1528 | ||
1529 | // set flat distribution (no decision) | |
1530 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1531 | ||
1532 | if ((track->GetStatus()&AliVTrack::kITSin)==0 && | |
1533 | (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal; | |
1534 | ||
1535 | //check for ITS standalone tracks | |
1536 | Bool_t isSA=kTRUE; | |
1537 | if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE; | |
1538 | ||
1539 | Double_t mom=track->P(); | |
1540 | Double_t dedx=track->GetITSsignal(); | |
1541 | Double_t momITS=mom; | |
1542 | UChar_t clumap=track->GetITSClusterMap(); | |
1543 | Int_t nPointsForPid=0; | |
1544 | for(Int_t i=2; i<6; i++){ | |
1545 | if(clumap&(1<<i)) ++nPointsForPid; | |
1546 | } | |
1547 | ||
1548 | if(nPointsForPid<3) { // track not to be used for combined PID purposes | |
1549 | // track->ResetStatus(AliVTrack::kITSpid); | |
1550 | return kDetNoSignal; | |
1551 | } | |
1552 | ||
1553 | Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; | |
1554 | for (Int_t j=0; j<AliPID::kSPECIES; j++) { | |
1555 | Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2 | |
1556 | const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.); | |
1557 | Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor; | |
1558 | //TODO: in case of the electron, use the SA parametrisation, | |
1559 | // this needs to be changed if ITS provides a parametrisation | |
1560 | // for electrons also for ITS+TPC tracks | |
1561 | Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron)); | |
1562 | if (TMath::Abs(dedx-bethe) > fRange*sigma) { | |
1563 | p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; | |
1564 | } else { | |
1565 | p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; | |
1566 | mismatch=kFALSE; | |
1567 | } | |
1568 | ||
1569 | // Check for particles heavier than (AliPID::kSPECIES - 1) | |
1570 | // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE; | |
1571 | ||
1572 | } | |
1573 | ||
1574 | if (mismatch){ | |
1575 | for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES; | |
1576 | return kDetNoSignal; | |
1577 | } | |
1578 | ||
1579 | ||
1580 | return kDetPidOk; | |
1581 | } | |
1582 | //______________________________________________________________________________ | |
1583 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1584 | { | |
1585 | // | |
1586 | // Compute PID response for the TPC | |
1587 | // | |
1588 | ||
1589 | // set flat distribution (no decision) | |
1590 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1591 | ||
1592 | // check quality of the track | |
1593 | if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal; | |
1594 | ||
1595 | Double_t mom = track->GetTPCmomentum(); | |
1596 | ||
1597 | Double_t dedx=track->GetTPCsignal(); | |
1598 | Bool_t mismatch=kTRUE/*, heavy=kTRUE*/; | |
1599 | ||
1600 | if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track); | |
1601 | ||
1602 | for (Int_t j=0; j<AliPID::kSPECIESC; j++) { | |
1603 | AliPID::EParticleType type=AliPID::EParticleType(j); | |
1604 | Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type); | |
1605 | Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type); | |
1606 | if (TMath::Abs(dedx-bethe) > fRange*sigma) { | |
1607 | p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma; | |
1608 | } else { | |
1609 | p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma; | |
1610 | mismatch=kFALSE; | |
1611 | } | |
1612 | } | |
1613 | ||
1614 | if (mismatch){ | |
1615 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1616 | return kDetNoSignal; | |
1617 | } | |
1618 | ||
1619 | return kDetPidOk; | |
1620 | } | |
1621 | //______________________________________________________________________________ | |
1622 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1623 | { | |
1624 | // | |
1625 | // Compute PID response for the | |
1626 | // | |
1627 | ||
1628 | Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1) | |
1629 | ||
1630 | // set flat distribution (no decision) | |
1631 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1632 | ||
1633 | if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal; | |
1634 | if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal; | |
1635 | ||
1636 | Bool_t mismatch = kTRUE/*, heavy = kTRUE*/; | |
1637 | for (Int_t j=0; j<AliPID::kSPECIESC; j++) { | |
1638 | AliPID::EParticleType type=AliPID::EParticleType(j); | |
1639 | Double_t nsigmas=GetNumberOfSigmasTOF(track,type) + meanCorrFactor; | |
1640 | ||
1641 | Double_t expTime = fTOFResponse.GetExpectedSignal(track,type); | |
1642 | Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type)); | |
1643 | if (TMath::Abs(nsigmas) > (fRange+2)) { | |
1644 | if(nsigmas < fTOFtail) | |
1645 | p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig; | |
1646 | else | |
1647 | p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig; | |
1648 | } else{ | |
1649 | if(nsigmas < fTOFtail) | |
1650 | p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig; | |
1651 | else | |
1652 | p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig; | |
1653 | } | |
1654 | ||
1655 | if (TMath::Abs(nsigmas)<5.){ | |
1656 | Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type); | |
1657 | if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE; | |
1658 | } | |
1659 | } | |
1660 | ||
1661 | if (mismatch){ | |
1662 | return kDetMismatch; | |
1663 | } | |
1664 | ||
1665 | // TODO: Light nuclei | |
1666 | ||
1667 | return kDetPidOk; | |
1668 | } | |
1669 | //______________________________________________________________________________ | |
1670 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod/*=AliTRDPIDResponse::kLQ1D*/) const | |
1671 | { | |
1672 | // | |
1673 | // Compute PID response for the | |
1674 | // | |
1675 | ||
1676 | UInt_t TRDslicesForPID[2]; | |
1677 | SetTRDSlices(TRDslicesForPID,PIDmethod); | |
1678 | // set flat distribution (no decision) | |
1679 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1680 | if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal; | |
1681 | ||
1682 | Float_t mom[6]={0.}; | |
1683 | Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices | |
1684 | Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1; | |
1685 | AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", TRDslicesForPID[0], TRDslicesForPID[1], nslices)); | |
1686 | for(UInt_t ilayer = 0; ilayer < 6; ilayer++){ | |
1687 | mom[ilayer] = track->GetTRDmomentum(ilayer); | |
1688 | for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){ | |
1689 | dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice); | |
1690 | } | |
1691 | } | |
1692 | fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod); | |
1693 | return kDetPidOk; | |
1694 | ||
1695 | } | |
1696 | //______________________________________________________________________________ | |
1697 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1698 | { | |
1699 | // | |
1700 | // Compute PID response for the EMCAL | |
1701 | // | |
1702 | ||
1703 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1704 | ||
1705 | AliVCluster *matchedClus = NULL; | |
1706 | ||
1707 | Double_t mom = -1.; | |
1708 | Double_t pt = -1.; | |
1709 | Double_t EovP = -1.; | |
1710 | Double_t fClsE = -1.; | |
1711 | ||
1712 | Int_t nMatchClus = -1; | |
1713 | Int_t charge = 0; | |
1714 | ||
1715 | // Track matching | |
1716 | nMatchClus = track->GetEMCALcluster(); | |
1717 | ||
1718 | if(nMatchClus > -1){ | |
1719 | ||
1720 | mom = track->P(); | |
1721 | pt = track->Pt(); | |
1722 | charge = track->Charge(); | |
1723 | ||
1724 | matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus); | |
1725 | ||
1726 | if(matchedClus){ | |
1727 | ||
1728 | // matched cluster is EMCAL | |
1729 | if(matchedClus->IsEMCAL()){ | |
1730 | ||
1731 | fClsE = matchedClus->E(); | |
1732 | EovP = fClsE/mom; | |
1733 | ||
1734 | ||
1735 | // compute the probabilities | |
1736 | if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){ | |
1737 | ||
1738 | // in case everything is OK | |
1739 | return kDetPidOk; | |
1740 | } | |
1741 | } | |
1742 | } | |
1743 | } | |
1744 | ||
1745 | // in all other cases set flat distribution (no decision) | |
1746 | for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies; | |
1747 | return kDetNoSignal; | |
1748 | ||
1749 | } | |
1750 | //______________________________________________________________________________ | |
1751 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const | |
1752 | { | |
1753 | // | |
1754 | // Compute PID response for the PHOS | |
1755 | // | |
1756 | ||
1757 | // set flat distribution (no decision) | |
1758 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1759 | return kDetNoSignal; | |
1760 | } | |
1761 | //______________________________________________________________________________ | |
1762 | AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const | |
1763 | { | |
1764 | // | |
1765 | // Compute PID response for the HMPID | |
1766 | // | |
1767 | // set flat distribution (no decision) | |
1768 | for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies; | |
1769 | if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal; | |
1770 | ||
1771 | track->GetHMPIDpid(p); | |
1772 | ||
1773 | return kDetPidOk; | |
1774 | } |