Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /* $Id: AliPIDResponse.cxx 46193 2010-12-21 09:00:14Z wiechula $ */
17 :
18 : //-----------------------------------------------------------------
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 //
24 : //-----------------------------------------------------------------
25 :
26 : #include <TList.h>
27 : #include <TObjArray.h>
28 : #include <TPRegexp.h>
29 : #include <TF1.h>
30 : #include <TH2D.h>
31 : #include <TSpline.h>
32 : #include <TFile.h>
33 : #include <TArrayI.h>
34 : #include <TArrayF.h>
35 : #include <TLinearFitter.h>
36 : #include <TSystem.h>
37 : #include <TMD5.h>
38 : #include "TRandom.h"
39 :
40 : #include <AliVEvent.h>
41 : #include <AliVTrack.h>
42 : #include <AliMCEvent.h>
43 : #include <AliLog.h>
44 : #include <AliPID.h>
45 : #include <AliOADBContainer.h>
46 : #include <AliTRDPIDResponseObject.h>
47 : #include <AliTRDdEdxParams.h>
48 : #include <AliTOFPIDParams.h>
49 : #include <AliHMPIDPIDParams.h>
50 :
51 : #include "AliPIDResponse.h"
52 : #include "AliDetectorPID.h"
53 :
54 : #include "AliCentrality.h"
55 :
56 176 : ClassImp(AliPIDResponse);
57 :
58 : Float_t AliPIDResponse::fgTOFmismatchProb = 0.0;
59 :
60 : AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
61 10 : TNamed("PIDResponse","PIDResponse"),
62 10 : fITSResponse(isMC),
63 10 : fTPCResponse(),
64 10 : fTRDResponse(),
65 10 : fTOFResponse(),
66 10 : fHMPIDResponse(),
67 10 : fEMCALResponse(),
68 10 : fRange(5.),
69 10 : fITSPIDmethod(kITSTruncMean),
70 10 : fTuneMConData(kFALSE),
71 10 : fTuneMConDataMask(kDetTOF|kDetTPC),
72 10 : fIsMC(isMC),
73 10 : fCachePID(kFALSE),
74 10 : fOADBPath(),
75 10 : fCustomTPCpidResponse(),
76 10 : fCustomTPCpidResponseOADBFile(),
77 10 : fCustomTPCetaMaps(),
78 10 : fBeamType("PP"),
79 10 : fLHCperiod(),
80 10 : fMCperiodTPC(),
81 10 : fMCperiodUser(),
82 10 : fCurrentFile(),
83 10 : fCurrentAliRootRev(-1),
84 10 : fRecoPass(0),
85 10 : fRecoPassUser(-1),
86 10 : fRun(-1),
87 10 : fOldRun(-1),
88 10 : fResT0A(75.),
89 10 : fResT0C(65.),
90 10 : fResT0AC(55.),
91 10 : fTPCPIDResponseArray(NULL),
92 10 : fArrPidResponseMaster(NULL),
93 10 : fResolutionCorrection(NULL),
94 10 : fOADBvoltageMaps(NULL),
95 10 : fUseTPCEtaCorrection(kFALSE),
96 10 : fUseTPCMultiplicityCorrection(kFALSE),
97 10 : fUseTPCNewResponse(kTRUE),
98 10 : fTRDPIDResponseObject(NULL),
99 10 : fTRDdEdxParams(NULL),
100 10 : fUseTRDEtaCorrection(kFALSE),
101 10 : fTOFtail(0.9),
102 10 : fTOFPIDParams(NULL),
103 10 : fHMPIDPIDParams(NULL),
104 10 : fEMCALPIDParams(NULL),
105 10 : fCurrentEvent(NULL),
106 10 : fCurrentMCEvent(NULL),
107 10 : fCurrCentrality(0.0),
108 10 : fBeamTypeNum(kPP),
109 10 : fNoTOFmism(kFALSE)
110 20 : {
111 : //
112 : // default ctor
113 : //
114 10 : AliLog::SetClassDebugLevel("AliPIDResponse",0);
115 10 : AliLog::SetClassDebugLevel("AliESDpid",0);
116 10 : AliLog::SetClassDebugLevel("AliAODpidUtil",0);
117 :
118 10 : }
119 :
120 : //______________________________________________________________________________
121 : AliPIDResponse::~AliPIDResponse()
122 20 : {
123 : //
124 : // dtor
125 : //
126 10 : delete fTPCPIDResponseArray;
127 10 : delete fArrPidResponseMaster;
128 10 : delete fTRDPIDResponseObject;
129 10 : delete fTRDdEdxParams;
130 10 : delete fTOFPIDParams;
131 10 : }
132 :
133 : //______________________________________________________________________________
134 : AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
135 0 : TNamed(other),
136 0 : fITSResponse(other.fITSResponse),
137 0 : fTPCResponse(other.fTPCResponse),
138 0 : fTRDResponse(other.fTRDResponse),
139 0 : fTOFResponse(other.fTOFResponse),
140 0 : fHMPIDResponse(other.fHMPIDResponse),
141 0 : fEMCALResponse(other.fEMCALResponse),
142 0 : fRange(other.fRange),
143 0 : fITSPIDmethod(other.fITSPIDmethod),
144 0 : fTuneMConData(other.fTuneMConData),
145 0 : fTuneMConDataMask(other.fTuneMConDataMask),
146 0 : fIsMC(other.fIsMC),
147 0 : fCachePID(other.fCachePID),
148 0 : fOADBPath(other.fOADBPath),
149 0 : fCustomTPCpidResponse(other.fCustomTPCpidResponse),
150 0 : fCustomTPCpidResponseOADBFile(other.fCustomTPCpidResponseOADBFile),
151 0 : fCustomTPCetaMaps(other.fCustomTPCetaMaps),
152 0 : fBeamType("PP"),
153 0 : fLHCperiod(),
154 0 : fMCperiodTPC(),
155 0 : fMCperiodUser(other.fMCperiodUser),
156 0 : fCurrentFile(),
157 0 : fCurrentAliRootRev(other.fCurrentAliRootRev),
158 0 : fRecoPass(0),
159 0 : fRecoPassUser(other.fRecoPassUser),
160 0 : fRun(-1),
161 0 : fOldRun(-1),
162 0 : fResT0A(75.),
163 0 : fResT0C(65.),
164 0 : fResT0AC(55.),
165 0 : fTPCPIDResponseArray(NULL),
166 0 : fArrPidResponseMaster(NULL),
167 0 : fResolutionCorrection(NULL),
168 0 : fOADBvoltageMaps(NULL),
169 0 : fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
170 0 : fUseTPCMultiplicityCorrection(other.fUseTPCMultiplicityCorrection),
171 0 : fUseTPCNewResponse(other.fUseTPCNewResponse),
172 0 : fTRDPIDResponseObject(NULL),
173 0 : fTRDdEdxParams(NULL),
174 0 : fUseTRDEtaCorrection(other.fUseTRDEtaCorrection),
175 0 : fTOFtail(0.9),
176 0 : fTOFPIDParams(NULL),
177 0 : fHMPIDPIDParams(NULL),
178 0 : fEMCALPIDParams(NULL),
179 0 : fCurrentEvent(NULL),
180 0 : fCurrentMCEvent(NULL),
181 0 : fCurrCentrality(0.0),
182 0 : fBeamTypeNum(kPP),
183 0 : fNoTOFmism(other.fNoTOFmism)
184 0 : {
185 : //
186 : // copy ctor
187 : //
188 0 : }
189 :
190 : //______________________________________________________________________________
191 : AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
192 : {
193 : //
194 : // copy ctor
195 : //
196 0 : if(this!=&other) {
197 0 : delete fArrPidResponseMaster;
198 0 : TNamed::operator=(other);
199 0 : fITSResponse=other.fITSResponse;
200 0 : fTPCResponse=other.fTPCResponse;
201 0 : fTRDResponse=other.fTRDResponse;
202 0 : fTOFResponse=other.fTOFResponse;
203 0 : fHMPIDResponse=other.fHMPIDResponse;
204 0 : fEMCALResponse=other.fEMCALResponse;
205 0 : fRange=other.fRange;
206 0 : fITSPIDmethod=other.fITSPIDmethod;
207 0 : fOADBPath=other.fOADBPath;
208 0 : fCustomTPCpidResponse=other.fCustomTPCpidResponse;
209 0 : fCustomTPCpidResponseOADBFile=other.fCustomTPCpidResponseOADBFile;
210 0 : fCustomTPCetaMaps=other.fCustomTPCetaMaps;
211 0 : fTuneMConData=other.fTuneMConData;
212 0 : fTuneMConDataMask=other.fTuneMConDataMask;
213 0 : fIsMC=other.fIsMC;
214 0 : fCachePID=other.fCachePID;
215 0 : fBeamType="PP";
216 0 : fBeamTypeNum=kPP;
217 0 : fLHCperiod="";
218 0 : fMCperiodTPC="";
219 0 : fMCperiodUser=other.fMCperiodUser;
220 0 : fCurrentFile="";
221 0 : fCurrentAliRootRev=other.fCurrentAliRootRev;
222 0 : fRecoPass=0;
223 0 : fRecoPassUser=other.fRecoPassUser;
224 0 : fRun=-1;
225 0 : fOldRun=-1;
226 0 : fResT0A=75.;
227 0 : fResT0C=65.;
228 0 : fResT0AC=55.;
229 0 : fTPCPIDResponseArray=NULL;
230 0 : fArrPidResponseMaster=NULL;
231 0 : fResolutionCorrection=NULL;
232 0 : fOADBvoltageMaps=NULL;
233 0 : fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
234 0 : fUseTPCMultiplicityCorrection=other.fUseTPCMultiplicityCorrection;
235 0 : fUseTPCNewResponse=other.fUseTPCNewResponse;
236 0 : fTRDPIDResponseObject=NULL;
237 0 : fTRDdEdxParams=NULL;
238 0 : fUseTRDEtaCorrection=other.fUseTRDEtaCorrection;
239 0 : fEMCALPIDParams=NULL;
240 0 : fTOFtail=0.9;
241 0 : fTOFPIDParams=NULL;
242 0 : fHMPIDPIDParams=NULL;
243 0 : fCurrentEvent=other.fCurrentEvent;
244 0 : fCurrentMCEvent=other.fCurrentMCEvent;
245 0 : fNoTOFmism = other.fNoTOFmism;
246 :
247 0 : }
248 0 : return *this;
249 : }
250 :
251 : //______________________________________________________________________________
252 : Float_t AliPIDResponse::NumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
253 : {
254 : //
255 : // NumberOfSigmas for 'detCode'
256 : //
257 :
258 0 : const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
259 : // look for cached value first
260 0 : const AliDetectorPID *detPID=track->GetDetectorPID();
261 :
262 0 : if ( detPID && detPID->HasNumberOfSigmas(detector)){
263 0 : return detPID->GetNumberOfSigmas(detector, type);
264 0 : } else if (fCachePID) {
265 0 : FillTrackDetectorPID(track, detector);
266 0 : detPID=track->GetDetectorPID();
267 0 : return detPID->GetNumberOfSigmas(detector, type);
268 : }
269 :
270 0 : return GetNumberOfSigmas(detector, track, type);
271 0 : }
272 :
273 : //______________________________________________________________________________
274 : AliPIDResponse::EDetPidStatus AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track,
275 : AliPID::EParticleType type, Double_t &val) const
276 : {
277 : //
278 : // NumberOfSigmas with detector status as return value
279 : //
280 :
281 0 : val=NumberOfSigmas(detCode, track, type);
282 0 : return CheckPIDStatus(detCode, (AliVTrack*)track);
283 : }
284 :
285 : //______________________________________________________________________________
286 : // public buffered versions of the PID calculation
287 : //
288 :
289 : //______________________________________________________________________________
290 : Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
291 : {
292 : //
293 : // Calculate the number of sigmas in the ITS
294 : //
295 :
296 0 : return NumberOfSigmas(kITS, vtrack, type);
297 : }
298 :
299 : //______________________________________________________________________________
300 : Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
301 : {
302 : //
303 : // Calculate the number of sigmas in the TPC
304 : //
305 :
306 0 : return NumberOfSigmas(kTPC, vtrack, type);
307 : }
308 :
309 : //______________________________________________________________________________
310 : Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
311 : AliPID::EParticleType type,
312 : AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
313 : {
314 : //get number of sigmas according the selected TPC gain configuration scenario
315 0 : const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
316 :
317 0 : Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
318 :
319 0 : return nSigma;
320 : }
321 :
322 : //______________________________________________________________________________
323 : Float_t AliPIDResponse::NumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
324 : {
325 : //
326 : // Calculate the number of sigmas in the TRD
327 : //
328 0 : return NumberOfSigmas(kTRD, vtrack, type);
329 : }
330 :
331 : //______________________________________________________________________________
332 : Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
333 : {
334 : //
335 : // Calculate the number of sigmas in the TOF
336 : //
337 :
338 0 : return NumberOfSigmas(kTOF, vtrack, type);
339 : }
340 :
341 : //______________________________________________________________________________
342 : Float_t AliPIDResponse::NumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
343 : {
344 : //
345 : // Calculate the number of sigmas in the EMCAL
346 : //
347 :
348 0 : return NumberOfSigmas(kHMPID, vtrack, type);
349 : }
350 :
351 : //______________________________________________________________________________
352 : Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
353 : {
354 : //
355 : // Calculate the number of sigmas in the EMCAL
356 : //
357 :
358 0 : return NumberOfSigmas(kEMCAL, vtrack, type);
359 : }
360 :
361 : //______________________________________________________________________________
362 : Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const
363 : {
364 : //
365 : // emcal nsigma with eop and showershape
366 : //
367 0 : AliVTrack *track=(AliVTrack*)vtrack;
368 :
369 : AliVCluster *matchedClus = NULL;
370 :
371 : Double_t mom = -1.;
372 : Double_t pt = -1.;
373 : Double_t EovP = -1.;
374 : Double_t fClsE = -1.;
375 :
376 : // initialize eop and shower shape parameters
377 0 : eop = -1.;
378 0 : for(Int_t i = 0; i < 4; i++){
379 0 : showershape[i] = -1.;
380 : }
381 :
382 : Int_t nMatchClus = -1;
383 : Int_t charge = 0;
384 :
385 : // Track matching
386 0 : nMatchClus = track->GetEMCALcluster();
387 0 : if(nMatchClus > -1){
388 :
389 0 : mom = track->P();
390 0 : pt = track->Pt();
391 0 : charge = track->Charge();
392 :
393 0 : matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
394 :
395 0 : if(matchedClus){
396 :
397 : // matched cluster is EMCAL
398 0 : if(matchedClus->IsEMCAL()){
399 :
400 0 : fClsE = matchedClus->E();
401 0 : EovP = fClsE/mom;
402 :
403 : // fill used EMCAL variables here
404 0 : eop = EovP; // E/p
405 0 : showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
406 0 : showershape[1] = matchedClus->GetM02(); // long axis
407 0 : showershape[2] = matchedClus->GetM20(); // short axis
408 0 : showershape[3] = matchedClus->GetDispersion(); // dispersion
409 :
410 : // look for cached value first
411 0 : const AliDetectorPID *detPID=track->GetDetectorPID();
412 : const EDetector detector=kEMCAL;
413 :
414 0 : if ( detPID && detPID->HasNumberOfSigmas(detector)){
415 0 : return detPID->GetNumberOfSigmas(detector, type);
416 0 : } else if (fCachePID) {
417 0 : FillTrackDetectorPID(track, detector);
418 0 : detPID=track->GetDetectorPID();
419 0 : return detPID->GetNumberOfSigmas(detector, type);
420 : }
421 :
422 : // NSigma value really meaningful only for electrons!
423 0 : return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
424 : }
425 : }
426 : }
427 0 : return -999;
428 0 : }
429 :
430 : //______________________________________________________________________________
431 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDelta(EDetector detector, const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
432 : {
433 : //
434 : //
435 : //
436 0 : val=-9999.;
437 0 : switch (detector){
438 0 : case kITS: return GetSignalDeltaITS(track,type,val,ratio); break;
439 0 : case kTPC: return GetSignalDeltaTPC(track,type,val,ratio); break;
440 0 : case kTRD: return GetSignalDeltaTRD(track,type,val,ratio); break;
441 0 : case kTOF: return GetSignalDeltaTOF(track,type,val,ratio); break;
442 0 : case kHMPID: return GetSignalDeltaHMPID(track,type,val,ratio); break;
443 0 : default: return kDetNoSignal;
444 : }
445 : return kDetNoSignal;
446 0 : }
447 :
448 : //______________________________________________________________________________
449 : Double_t AliPIDResponse::GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
450 : {
451 : //
452 : //
453 : //
454 0 : Double_t val=-9999.;
455 0 : EDetPidStatus stat=GetSignalDelta(detCode, track, type, val, ratio);
456 0 : if ( stat==kDetNoSignal ) val=-9999.;
457 0 : return val;
458 0 : }
459 :
460 : //______________________________________________________________________________
461 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
462 : {
463 : // Compute PID response of 'detCode'
464 :
465 : // find detector code from detector bit mask
466 : Int_t detector=-1;
467 0 : for (Int_t idet=0; idet<kNdetectors; ++idet) if ( (detCode&(1<<idet)) ) { detector=idet; break; }
468 0 : if (detector==-1) return kDetNoSignal;
469 :
470 0 : return ComputePIDProbability((EDetector)detector, track, nSpecies, p);
471 0 : }
472 :
473 : //______________________________________________________________________________
474 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detector, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
475 : {
476 : //
477 : // Compute PID response of 'detector'
478 : //
479 :
480 576 : const AliDetectorPID *detPID=track->GetDetectorPID();
481 :
482 288 : if ( detPID && detPID->HasRawProbability(detector)){
483 0 : return detPID->GetRawProbability(detector, p, nSpecies);
484 288 : } else if (fCachePID) {
485 0 : FillTrackDetectorPID(track, detector);
486 0 : detPID=track->GetDetectorPID();
487 0 : return detPID->GetRawProbability(detector, p, nSpecies);
488 : }
489 :
490 : //if no caching return values calculated from scratch
491 288 : return GetComputePIDProbability(detector, track, nSpecies, p);
492 288 : }
493 :
494 : //______________________________________________________________________________
495 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
496 : {
497 : // Compute PID response for the ITS
498 0 : return ComputePIDProbability(kITS, track, nSpecies, p);
499 : }
500 :
501 : //______________________________________________________________________________
502 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
503 : {
504 : // Compute PID response for the TPC
505 0 : return ComputePIDProbability(kTPC, track, nSpecies, p);
506 : }
507 :
508 : //______________________________________________________________________________
509 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
510 : {
511 : // Compute PID response for the
512 0 : return ComputePIDProbability(kTOF, track, nSpecies, p);
513 : }
514 :
515 : //______________________________________________________________________________
516 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
517 : {
518 : // Compute PID response for the
519 0 : return ComputePIDProbability(kTRD, track, nSpecies, p);
520 : }
521 :
522 : //______________________________________________________________________________
523 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
524 : {
525 : // Compute PID response for the EMCAL
526 0 : return ComputePIDProbability(kEMCAL, track, nSpecies, p);
527 : }
528 : //______________________________________________________________________________
529 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
530 : {
531 : // Compute PID response for the PHOS
532 :
533 : // set flat distribution (no decision)
534 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
535 0 : return kDetNoSignal;
536 : }
537 :
538 : //______________________________________________________________________________
539 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
540 : {
541 : // Compute PID response for the HMPID
542 0 : return ComputePIDProbability(kHMPID, track, nSpecies, p);
543 : }
544 :
545 : //______________________________________________________________________________
546 : AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
547 : {
548 : // Compute PID response for the
549 0 : return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
550 : }
551 :
552 : //______________________________________________________________________________
553 : AliPIDResponse::EDetPidStatus AliPIDResponse::CheckPIDStatus(EDetector detector, const AliVTrack *track) const
554 : {
555 : // calculate detector pid status
556 :
557 : const Int_t iDetCode=(Int_t)detector;
558 0 : if (iDetCode<0||iDetCode>=kNdetectors) return kDetNoSignal;
559 0 : const AliDetectorPID *detPID=track->GetDetectorPID();
560 :
561 0 : if ( detPID ){
562 0 : return detPID->GetPIDStatus(detector);
563 0 : } else if (fCachePID) {
564 0 : FillTrackDetectorPID(track, detector);
565 0 : detPID=track->GetDetectorPID();
566 0 : return detPID->GetPIDStatus(detector);
567 : }
568 :
569 : // if not buffered and no buffering is requested
570 0 : return GetPIDStatus(detector, track);
571 0 : }
572 :
573 : //______________________________________________________________________________
574 : void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
575 : {
576 : //
577 : // Apply settings for the current event
578 : //
579 0 : fRecoPass=pass;
580 :
581 :
582 0 : fCurrentEvent=NULL;
583 0 : if (!event) return;
584 0 : fCurrentEvent=event;
585 0 : if (run>0) fRun=run;
586 0 : else fRun=event->GetRunNumber();
587 :
588 0 : if (fRun!=fOldRun){
589 0 : ExecNewRun();
590 0 : fOldRun=fRun;
591 0 : }
592 :
593 : //TPC resolution parametrisation PbPb
594 0 : if ( fResolutionCorrection ){
595 0 : Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
596 0 : fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
597 0 : }
598 :
599 : // Set up TPC multiplicity for PbPb
600 0 : if (fUseTPCMultiplicityCorrection) {
601 0 : Int_t numESDtracks = event->GetNumberOfESDTracks();
602 0 : if (numESDtracks < 0) {
603 0 : AliError("Cannot obtain event multiplicity (number of ESD tracks < 0). If you are using AODs, this might be a too old production. Please disable the multiplicity correction to get a reliable PID result!");
604 : numESDtracks = 0;
605 0 : }
606 0 : fTPCResponse.SetCurrentEventMultiplicity(numESDtracks);
607 0 : }
608 : else
609 0 : fTPCResponse.SetCurrentEventMultiplicity(0);
610 :
611 : //TOF resolution
612 0 : SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
613 :
614 :
615 : // Get and set centrality
616 0 : AliCentrality *centrality = event->GetCentrality();
617 0 : if(centrality){
618 0 : fCurrCentrality = centrality->GetCentralityPercentile("V0M");
619 0 : }
620 : else{
621 0 : fCurrCentrality = -1;
622 : }
623 :
624 : // Set centrality percentile for EMCAL
625 0 : fEMCALResponse.SetCentrality(fCurrCentrality);
626 :
627 : // switch off some TOF channel according to OADB to match data TOF matching eff
628 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) && fTOFPIDParams->GetTOFmatchingLossMC() > 0.01){
629 0 : Int_t ntrk = event->GetNumberOfTracks();
630 0 : for(Int_t i=0;i < ntrk;i++){
631 0 : AliVParticle *trk = event->GetTrack(i);
632 0 : Int_t channel = GetTOFResponse().GetTOFchannel(trk);
633 0 : Int_t swoffEachOfThem = Int_t(100./fTOFPIDParams->GetTOFmatchingLossMC() + 0.5);
634 0 : if(!(channel%swoffEachOfThem)) ((AliVTrack *) trk)->ResetStatus(AliVTrack::kTOFout);
635 : }
636 0 : }
637 :
638 0 : }
639 :
640 : //______________________________________________________________________________
641 : void AliPIDResponse::ExecNewRun()
642 : {
643 : //
644 : // Things to Execute upon a new run
645 : //
646 0 : SetRecoInfo();
647 :
648 : // ===| ITS part |============================================================
649 0 : SetITSParametrisation();
650 :
651 : // ===| TPC part |============================================================
652 : // new treatment for loading the TPC PID response if requested
653 : // for the moment fall back to old method if no PID response array is found
654 : // by the new method for backward compatibility
655 : Bool_t doOldTPCPID=kTRUE;
656 0 : if (fUseTPCNewResponse) {
657 0 : doOldTPCPID = !InitializeTPCResponse();
658 0 : if (doOldTPCPID) {
659 0 : AliWarning("No TPC response parametrisations found using the new method. Falling back to the old method.");
660 0 : }
661 : }
662 :
663 0 : if (doOldTPCPID) {
664 0 : SetTPCPidResponseMaster();
665 0 : SetTPCParametrisation();
666 0 : }
667 0 : SetTPCEtaMaps();
668 :
669 : // ===| TRD part |============================================================
670 0 : SetTRDPidResponseMaster();
671 : //has to precede InitializeTRDResponse(), otherwise the read-out fTRDdEdxParams is not pased in TRDResponse!
672 0 : SetTRDdEdxParams();
673 0 : SetTRDEtaMaps();
674 0 : InitializeTRDResponse();
675 :
676 : // ===| TOF part |============================================================
677 0 : SetTOFPidResponseMaster();
678 0 : InitializeTOFResponse();
679 :
680 : // ===| EMCAL part |==========================================================
681 0 : SetEMCALPidResponseMaster();
682 0 : InitializeEMCALResponse();
683 :
684 : // ===| HMPID part |==========================================================
685 0 : SetHMPIDPidResponseMaster();
686 0 : InitializeHMPIDResponse();
687 :
688 0 : if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
689 0 : if (fCurrentEvent) fTRDResponse.SetMagField(fCurrentEvent->GetMagneticField());
690 0 : }
691 :
692 : //______________________________________________________________________________
693 : Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
694 : {
695 : //
696 : // Get TPC multiplicity in bins of 150
697 : //
698 :
699 0 : const AliVVertex* vertexTPC = event->GetPrimaryVertex();
700 : Double_t tpcMulti=0.;
701 0 : if(vertexTPC){
702 0 : Double_t vertexContribTPC=vertexTPC->GetNContributors();
703 0 : tpcMulti=vertexContribTPC/150.;
704 0 : if (tpcMulti>20.) tpcMulti=20.;
705 0 : }
706 :
707 0 : return tpcMulti;
708 : }
709 :
710 : //______________________________________________________________________________
711 : void AliPIDResponse::SetRecoInfo()
712 : {
713 : //
714 : // Set reconstruction information
715 : //
716 :
717 : //reset information
718 0 : fLHCperiod="";
719 0 : fMCperiodTPC="";
720 :
721 0 : fBeamType="";
722 :
723 0 : fBeamType="PP";
724 0 : fBeamTypeNum=kPP;
725 :
726 0 : Bool_t hasProdInfo=(fCurrentFile.BeginsWith("LHC"));
727 :
728 0 : TPRegexp reg(".*(LHC1[1-3][a-z]+[0-9]+[a-z_]*)[/_].*");
729 0 : if (hasProdInfo) reg=TPRegexp("LHC1[1-2][a-z]+[0-9]+[a-z_]*");
730 0 : TPRegexp reg12a17("LHC1[2-4][a-z]");
731 :
732 : //find the period by run number (UGLY, but not stored in ESD and AOD... )
733 0 : if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
734 0 : else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
735 0 : else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
736 0 : else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
737 0 : else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
738 0 : else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
739 0 : else if (fRun>=136851&&fRun<=139846) {
740 0 : fLHCperiod="LHC10H";
741 0 : fMCperiodTPC="LHC10H8";
742 0 : if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
743 : // exception for 13d2 and later
744 0 : if (fCurrentAliRootRev >= 62714) fMCperiodTPC="LHC13D2";
745 0 : fBeamType="PBPB";
746 0 : fBeamTypeNum=kPBPB;
747 0 : }
748 0 : else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
749 : //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
750 0 : else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
751 0 : else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
752 : // also for 11e (159650-162750),f(162751-165771) use 11d
753 0 : else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
754 0 : else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
755 :
756 0 : else if (fRun>=165772 && fRun<=170718) {
757 0 : fLHCperiod="LHC11H";
758 0 : fMCperiodTPC="LHC11A10";
759 0 : fBeamType="PBPB";
760 0 : fBeamTypeNum=kPBPB;
761 0 : if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
762 : }
763 0 : if (fRun>=170719 && fRun<=177311) {
764 0 : fLHCperiod="LHC12A";
765 0 : fBeamType="PP";
766 0 : fBeamTypeNum=kPP;
767 0 : fMCperiodTPC="LHC10F6A";
768 0 : if (fCurrentAliRootRev >= 62714)
769 0 : fMCperiodTPC="LHC14E2";
770 : }
771 : // for the moment use LHC12b parameters up to LHC12d
772 0 : if (fRun>=177312 /*&& fRun<=179356*/) {
773 0 : fLHCperiod="LHC12B";
774 0 : fBeamType="PP";
775 0 : fBeamTypeNum=kPP;
776 0 : fMCperiodTPC="LHC10F6A";
777 0 : if (fCurrentAliRootRev >= 62714)
778 0 : fMCperiodTPC="LHC14E2";
779 : }
780 : // if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
781 : // if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
782 : // if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
783 :
784 : // if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
785 : // if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
786 : // if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fBeamTypeNum=kPPB;/*fMCperiodTPC="";*/ }
787 : // for the moment use 12g parametrisation for all full gain runs (LHC12e+)
788 :
789 : // Dedicated splines for periods 12g and 12i(j) (and use more appropriate MC)
790 0 : if (fRun >= 188720 && fRun <= 192738) {
791 0 : fLHCperiod="LHC12H";
792 0 : fBeamType="PP";
793 0 : fBeamTypeNum=kPP;
794 0 : fMCperiodTPC="LHC10F6A";
795 0 : if (fCurrentAliRootRev >= 62714)
796 0 : fMCperiodTPC="LHC13B2_FIXn1";
797 : }
798 0 : if (fRun >= 192739 && fRun <= 194479) {
799 0 : fLHCperiod="LHC12I";
800 0 : fBeamType="PP";
801 0 : fBeamTypeNum=kPP;
802 0 : fMCperiodTPC="LHC10F6A";
803 0 : if (fCurrentAliRootRev >= 62714)
804 0 : fMCperiodTPC="LHC13B2_FIXn1";
805 : }
806 :
807 : // Use for pp and pPb 12E-G pass 1 the PPB runs
808 0 : if (fRecoPass==1 && fRun >= 186346 && fRun < 188719) { fLHCperiod="LHC12G"; fBeamType="PPB";fBeamTypeNum=kPPB; fMCperiodTPC="LHC12G"; }
809 :
810 : // settings for pass2 of 2012 data
811 0 : if (fRecoPass>=2 && fRun>=170719 && fRun<=194479) {
812 0 : fBeamType = "PP";
813 0 : fBeamTypeNum = kPP;
814 0 : fMCperiodTPC = "";
815 0 : if (fRun >= 170719 && fRun <= 177311) { fLHCperiod="LHC12A"; }
816 0 : if (fRun >= 177312 && fRun <= 179356) { fLHCperiod="LHC12B"; }
817 0 : if (fRun >= 179357 && fRun <= 183173) { fLHCperiod="LHC12C"; }
818 0 : if (fRun >= 183174 && fRun <= 186345) { fLHCperiod="LHC12D"; }
819 0 : if (fRun >= 186346 && fRun <= 186635) { fLHCperiod="LHC12E"; }
820 0 : if (fRun >= 186636 && fRun <= 188166) { fLHCperiod="LHC12F"; }
821 0 : if (fRun >= 188167 && fRun <= 188719) { fLHCperiod="LHC12G"; }
822 0 : if (fRun >= 188720 && fRun <= 192738) { fLHCperiod="LHC12H"; }
823 0 : if (fRun >= 192739 && fRun <= 193766) { fLHCperiod="LHC12I"; }
824 : // no special parametrisations for 12J, use 12I instead
825 0 : if (fRun >= 193767 && fRun <= 194479) { /*fLHCperiod="LHC12J";*/ fLHCperiod="LHC12I"; }
826 :
827 : // overwriting for the PPB period
828 0 : if (fRun >= 188167 && fRun <= 188418) { fLHCperiod="LHC12G"; fBeamType="PPB";fBeamTypeNum=kPPB; fMCperiodTPC="LHC12G"; }
829 : }
830 :
831 : // New parametrisation for 2013 pPb runs
832 0 : if (fRun >= 194480) {
833 0 : fLHCperiod="LHC13B";
834 0 : fBeamType="PPB";
835 0 : fBeamTypeNum=kPPB;
836 0 : fMCperiodTPC="LHC12G";
837 :
838 0 : if (fCurrentAliRootRev >= 61605)
839 0 : fMCperiodTPC="LHC13B2_FIX";
840 0 : if (fCurrentAliRootRev >= 62714)
841 0 : fMCperiodTPC="LHC13B2_FIXn1";
842 :
843 : // High luminosity pPb runs require different parametrisations
844 0 : if (fRun >= 195875 && fRun <= 197411) {
845 0 : fLHCperiod="LHC13F";
846 : }
847 : }
848 :
849 : // New parametrisation for the first 2015 pp runs
850 0 : if (fRun >= 208505) { // <<< This is the first run in 15a
851 0 : fLHCperiod="LHC15F";
852 0 : fBeamType="PP";
853 0 : fBeamTypeNum=kPP;
854 0 : fMCperiodTPC="LHC15G3";
855 : }
856 :
857 : //exception new pp MC productions from 2011 (11a periods have 10f6a splines!)
858 0 : if (fBeamType=="PP" && reg.MatchB(fCurrentFile) && !fCurrentFile.Contains("LHC11a")) { fMCperiodTPC="LHC11B2"; fBeamType="PP";fBeamTypeNum=kPP; }
859 : // exception for 11f1
860 0 : if (fCurrentFile.Contains("LHC11f1")) fMCperiodTPC="LHC11F1";
861 : // exception for 12f1a, 12f1b and 12i3
862 0 : if (fCurrentFile.Contains("LHC12f1") || fCurrentFile.Contains("LHC12i3")) fMCperiodTPC="LHC12F1";
863 : // exception for 12c4
864 0 : if (fCurrentFile.Contains("LHC12c4")) fMCperiodTPC="LHC12C4";
865 : // exception for 13d1 11d anchored prod
866 0 : if (fLHCperiod=="LHC11D" && fCurrentFile.Contains("LHC13d1")) fMCperiodTPC="LHC13D1";
867 0 : }
868 :
869 : //______________________________________________________________________________
870 : void AliPIDResponse::SetITSParametrisation()
871 : {
872 : //
873 : // Set the ITS parametrisation
874 : //
875 0 : }
876 :
877 :
878 : //______________________________________________________________________________
879 : void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
880 : {
881 0 : if (h->GetBinContent(binX, binY) <= 1e-4)
882 : return; // Reject bins without content (within some numerical precision) or with strange content
883 :
884 0 : Double_t coord[2] = {0, 0};
885 0 : coord[0] = h->GetXaxis()->GetBinCenter(binX);
886 0 : coord[1] = h->GetYaxis()->GetBinCenter(binY);
887 0 : Double_t binError = h->GetBinError(binX, binY);
888 0 : if (binError <= 0) {
889 : binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
890 0 : printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
891 0 : }
892 0 : linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
893 0 : }
894 :
895 :
896 : //______________________________________________________________________________
897 : TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
898 : {
899 0 : if (!h)
900 0 : return 0x0;
901 :
902 : // Interpolate to finer map
903 0 : TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
904 :
905 0 : Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
906 0 : Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
907 : Int_t nBinsX = 30;
908 : // Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
909 : // scale the number of bins correspondingly
910 0 : Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
911 0 : Int_t nBinsXrefined = nBinsX * refineFactorX;
912 0 : Int_t nBinsYrefined = nBinsY * refineFactorY;
913 :
914 0 : TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()), Form("%s (refined)", h->GetTitle()),
915 0 : nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
916 : nBinsYrefined, lowerMapBoundY, upperMapBoundY);
917 :
918 0 : for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
919 0 : for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
920 :
921 0 : hRefined->SetBinContent(binX, binY, 1); // Default value is 1
922 :
923 0 : Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
924 0 : Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
925 :
926 : /*OLD
927 : linExtrapolation->ClearPoints();
928 :
929 : // For interpolation: Just take the corresponding bin from the old histo.
930 : // For extrapolation: take the last available bin from the old histo.
931 : // If the boundaries are to be skipped, also skip the corresponding bins
932 : Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
933 : if (oldBinX < 1)
934 : oldBinX = 1;
935 : if (oldBinX > nBinsX)
936 : oldBinX = nBinsX;
937 :
938 : Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
939 : if (oldBinY < 1)
940 : oldBinY = 1;
941 : if (oldBinY > nBinsY)
942 : oldBinY = nBinsY;
943 :
944 : // Neighbours left column
945 : if (oldBinX >= 2) {
946 : if (oldBinY >= 2) {
947 : AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
948 : }
949 :
950 : AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
951 :
952 : if (oldBinY < nBinsY) {
953 : AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
954 : }
955 : }
956 :
957 : // Neighbours (and point itself) same column
958 : if (oldBinY >= 2) {
959 : AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
960 : }
961 :
962 : AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
963 :
964 : if (oldBinY < nBinsY) {
965 : AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
966 : }
967 :
968 : // Neighbours right column
969 : if (oldBinX < nBinsX) {
970 : if (oldBinY >= 2) {
971 : AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
972 : }
973 :
974 : AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
975 :
976 : if (oldBinY < nBinsY) {
977 : AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
978 : }
979 : }
980 :
981 :
982 : // Fit 2D-hyperplane
983 : if (linExtrapolation->GetNpoints() <= 0)
984 : continue;
985 :
986 : if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
987 : continue;
988 :
989 : // Fill the bin of the refined histogram with the extrapolated value
990 : Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
991 : + linExtrapolation->GetParameter(2) * centerY;
992 : */
993 0 : Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
994 0 : hRefined->SetBinContent(binX, binY, interpolatedValue);
995 : }
996 : }
997 :
998 :
999 : // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
1000 : // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
1001 : // Assume line through these points and extropolate to last bin of refined map
1002 0 : const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
1003 0 : const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
1004 :
1005 0 : const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
1006 :
1007 0 : const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
1008 0 : const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
1009 :
1010 0 : for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
1011 0 : Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
1012 :
1013 0 : const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
1014 0 : const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
1015 :
1016 0 : const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
1017 : const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
1018 :
1019 :
1020 0 : const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
1021 0 : const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
1022 :
1023 0 : const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
1024 : const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
1025 :
1026 0 : for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
1027 0 : Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
1028 :
1029 0 : if (centerX < firstOldXbinCenter) {
1030 0 : Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
1031 0 : hRefined->SetBinContent(binX, binY, extrapolatedValue);
1032 0 : }
1033 0 : else if (centerX <= lastOldXbinCenter) {
1034 0 : continue;
1035 : }
1036 : else {
1037 0 : Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
1038 0 : hRefined->SetBinContent(binX, binY, extrapolatedValue);
1039 : }
1040 0 : }
1041 : }
1042 :
1043 0 : delete linExtrapolation;
1044 :
1045 : return hRefined;
1046 0 : }
1047 :
1048 : //______________________________________________________________________________
1049 : void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
1050 : Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
1051 : {
1052 : //
1053 : // Load the TPC eta correction maps from the OADB
1054 : //
1055 :
1056 0 : if (fUseTPCEtaCorrection == kFALSE) {
1057 : // Disable eta correction via setting no maps
1058 0 : if (!fTPCResponse.SetEtaCorrMap(0x0))
1059 0 : AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
1060 : else
1061 0 : AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
1062 :
1063 0 : if (!fTPCResponse.SetSigmaParams(0x0, 0))
1064 0 : AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
1065 : else
1066 0 : AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
1067 :
1068 : return;
1069 : }
1070 :
1071 0 : TString dataType = "DATA";
1072 0 : TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
1073 :
1074 0 : if (fIsMC) {
1075 0 : if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
1076 0 : period=fMCperiodTPC;
1077 0 : dataType="MC";
1078 : }
1079 0 : fRecoPass = 1;
1080 :
1081 0 : if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) && fMCperiodTPC.IsNull()) {
1082 0 : AliError("***** Risk for unreliable TPC PID detected: ********");
1083 0 : AliError(" MC detected, but no MC period set -> Not changing eta maps!");
1084 0 : return;
1085 : }
1086 : }
1087 :
1088 0 : Int_t recopass = fRecoPass;
1089 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) )
1090 0 : recopass = fRecoPassUser;
1091 :
1092 0 : TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
1093 :
1094 0 : AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
1095 :
1096 : // Invalidate old maps
1097 0 : fTPCResponse.SetEtaCorrMap(0x0);
1098 0 : fTPCResponse.SetSigmaParams(0x0, 0);
1099 :
1100 :
1101 0 : TString fileNameMaps(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
1102 0 : if (!fCustomTPCetaMaps.IsNull()) fileNameMaps=fCustomTPCetaMaps;
1103 :
1104 : // Load the eta correction maps
1105 0 : AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
1106 :
1107 0 : Int_t statusCont = etaMapsCont.InitFromFile(fileNameMaps.Data(), Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
1108 0 : if (statusCont) {
1109 0 : AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
1110 0 : fUseTPCEtaCorrection = kFALSE;
1111 0 : }
1112 : else {
1113 0 : AliInfo(Form("Loading TPC eta correction map from %s", fileNameMaps.Data()));
1114 :
1115 : TH2D* etaMap = 0x0;
1116 :
1117 0 : if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
1118 0 : TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
1119 0 : etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
1120 0 : if (!etaMap) {
1121 : // Try default object
1122 0 : etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
1123 0 : }
1124 0 : }
1125 : else {
1126 0 : etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
1127 : }
1128 :
1129 :
1130 0 : if (!etaMap) {
1131 0 : AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
1132 0 : fUseTPCEtaCorrection = kFALSE;
1133 0 : }
1134 : else {
1135 0 : TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
1136 :
1137 0 : if (etaMapRefined) {
1138 0 : if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
1139 0 : AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
1140 0 : fTPCResponse.SetEtaCorrMap(0x0);
1141 0 : fUseTPCEtaCorrection = kFALSE;
1142 0 : }
1143 : else {
1144 0 : AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s: %s (MD5(map) = %s)",
1145 : refineFactorMapX, refineFactorMapY, fileNameMaps.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle(),
1146 : AliTPCPIDResponse::GetChecksum(fTPCResponse.GetEtaCorrMap()).Data()));
1147 : }
1148 :
1149 0 : delete etaMapRefined;
1150 : }
1151 : else {
1152 0 : AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
1153 0 : fUseTPCEtaCorrection = kFALSE;
1154 : }
1155 : }
1156 : }
1157 :
1158 : // If there was some problem loading the eta maps, it makes no sense to load the sigma maps (that require eta corrected data)
1159 0 : if (fUseTPCEtaCorrection == kFALSE) {
1160 0 : AliError("Failed to load TPC eta correction map required by sigma maps -> Using old parametrisation for sigma");
1161 0 : return;
1162 : }
1163 :
1164 : // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
1165 0 : AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
1166 :
1167 0 : statusCont = etaSigmaMapsCont.InitFromFile(fileNameMaps.Data(), Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
1168 0 : if (statusCont) {
1169 0 : AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
1170 : }
1171 : else {
1172 0 : AliInfo(Form("Loading TPC eta sigma map from %s", fileNameMaps.Data()));
1173 :
1174 : TObjArray* etaSigmaPars = 0x0;
1175 :
1176 0 : if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
1177 0 : TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
1178 0 : etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
1179 0 : if (!etaSigmaPars) {
1180 : // Try default object
1181 0 : etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
1182 0 : }
1183 0 : }
1184 : else {
1185 0 : etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
1186 : }
1187 :
1188 0 : if (!etaSigmaPars) {
1189 0 : AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
1190 : }
1191 : else {
1192 0 : TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
1193 0 : TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
1194 : Double_t sigmaPar0 = 0.0;
1195 :
1196 0 : if (sigmaPar0Info) {
1197 0 : TString sigmaPar0String = sigmaPar0Info->GetTitle();
1198 0 : sigmaPar0 = sigmaPar0String.Atof();
1199 0 : }
1200 : else {
1201 : // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
1202 : etaSigmaPar1Map = 0x0;
1203 : }
1204 :
1205 0 : TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
1206 :
1207 :
1208 0 : if (etaSigmaPar1MapRefined) {
1209 0 : if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
1210 0 : AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
1211 0 : fTPCResponse.SetSigmaParams(0x0, 0);
1212 : }
1213 : else {
1214 0 : AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s: %s (MD5(map) = %s, sigmaPar0 = %f)",
1215 : refineFactorSigmaMapX, refineFactorSigmaMapY, fileNameMaps.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle(),
1216 : AliTPCPIDResponse::GetChecksum(fTPCResponse.GetSigmaPar1Map()).Data(), sigmaPar0));
1217 : }
1218 :
1219 0 : delete etaSigmaPar1MapRefined;
1220 : }
1221 : else {
1222 0 : AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
1223 : fRun));
1224 : }
1225 : }
1226 : }
1227 0 : }
1228 :
1229 :
1230 : //______________________________________________________________________________
1231 : Bool_t AliPIDResponse::InitializeTPCResponse()
1232 : {
1233 : // Load the Array with TPC PID response information
1234 : // This is the new method which will completely replace the old one at some point
1235 :
1236 : //
1237 : // Setup old resolution parametrisation
1238 : // TODO: This should be moved to the initialisation and vanish completely at some point
1239 :
1240 : //default
1241 0 : fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1242 :
1243 0 : if (fRun>=122195){ //LHC10d
1244 0 : fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1245 0 : }
1246 :
1247 0 : if (fRun>=170719){ // LHC12a
1248 0 : fTPCResponse.SetSigma(2.95714e-03, 1.01953e+05);
1249 0 : }
1250 :
1251 0 : if (fRun>=177312){ // LHC12b
1252 0 : fTPCResponse.SetSigma(3.74633e-03, 7.11829e+04 );
1253 0 : }
1254 :
1255 0 : if (fRun>=186346){ // LHC12e
1256 0 : fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1257 0 : }
1258 :
1259 :
1260 0 : AliInfo("---------------------------- TPC Response Configuration (New) ----------------------------");
1261 : // ===| load TPC response array from OADB |===================================
1262 0 : TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponseOADB.root", fOADBPath.Data()));
1263 0 : if (!fCustomTPCpidResponseOADBFile.IsNull()) fileNamePIDresponse=fCustomTPCpidResponseOADBFile;
1264 :
1265 :
1266 : // ---| In case of MC and NO tune on data fall back to old method |-----------
1267 0 : if (fIsMC) {
1268 0 : if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) return kFALSE;
1269 : }
1270 :
1271 : // ---| set reco pass |-------------------------------------------------------
1272 0 : Int_t recopass = fRecoPass;
1273 0 : if(fIsMC && fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser;
1274 :
1275 0 : const Bool_t returnValue = fTPCResponse.InitFromOADB(fRun, TString::Format("%d", recopass), fileNamePIDresponse, fUseTPCMultiplicityCorrection);
1276 0 : AliInfo("------------------------------------------------------------------------------------------");
1277 :
1278 0 : return returnValue;
1279 0 : }
1280 :
1281 : //______________________________________________________________________________
1282 : void AliPIDResponse::SetTPCPidResponseMaster()
1283 : {
1284 : //
1285 : // Load the TPC pid response functions from the OADB
1286 : // Load the TPC voltage maps from OADB
1287 : //
1288 : //don't load twice for the moment
1289 0 : if (fArrPidResponseMaster) return;
1290 :
1291 :
1292 : //reset the PID response functions
1293 0 : delete fArrPidResponseMaster;
1294 0 : fArrPidResponseMaster=NULL;
1295 :
1296 : TFile *f=NULL;
1297 :
1298 0 : TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
1299 0 : if (!fCustomTPCpidResponse.IsNull()) fileNamePIDresponse=fCustomTPCpidResponse;
1300 :
1301 0 : f=TFile::Open(fileNamePIDresponse.Data());
1302 0 : if (f && f->IsOpen() && !f->IsZombie()){
1303 0 : fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
1304 0 : }
1305 0 : delete f;
1306 :
1307 0 : TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
1308 0 : f=TFile::Open(fileNameVoltageMaps.Data());
1309 0 : if (f && f->IsOpen() && !f->IsZombie()){
1310 0 : fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
1311 0 : }
1312 0 : delete f;
1313 :
1314 0 : if (!fArrPidResponseMaster){
1315 0 : AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
1316 0 : return;
1317 : }
1318 0 : fArrPidResponseMaster->SetOwner();
1319 :
1320 0 : if (!fOADBvoltageMaps)
1321 : {
1322 0 : AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
1323 : }
1324 0 : fArrPidResponseMaster->SetOwner();
1325 0 : }
1326 :
1327 : //______________________________________________________________________________
1328 : void AliPIDResponse::SetTPCParametrisation()
1329 : {
1330 : //
1331 : // Change BB parametrisation for current run
1332 : //
1333 :
1334 0 : AliInfo("---------------------------- TPC Response Configuration (Old) ----------------------------");
1335 :
1336 : //
1337 : //reset old splines
1338 : //
1339 0 : fTPCResponse.ResetSplines();
1340 :
1341 0 : if (fLHCperiod.IsNull()) {
1342 0 : AliError("No period set, not changing parametrisation");
1343 0 : AliInfo("------------------------------------------------------------------------------------------");
1344 0 : return;
1345 : }
1346 :
1347 : //
1348 : // Set default parametrisations for data and MC
1349 : //
1350 :
1351 : //data type
1352 0 : TString datatype="DATA";
1353 : //in case of mc fRecoPass is per default 1
1354 0 : if (fIsMC) {
1355 0 : if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) datatype="MC";
1356 0 : fRecoPass=1;
1357 0 : } else {
1358 0 : if (fRecoPass<=0) {
1359 0 : fTPCResponse.SetUseDatabase(kFALSE);
1360 0 : AliError("******** Risk for unreliable TPC PID detected **********");
1361 0 : AliError(" no proper reco pass was set, no splines can be set");
1362 0 : AliError(" an outdate Bethe Bloch parametrisation will be used");
1363 0 : AliInfo("------------------------------------------------------------------------------------------");
1364 0 : return;
1365 : }
1366 : }
1367 :
1368 : // period
1369 0 : TString period=fLHCperiod;
1370 0 : if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) period=fMCperiodTPC;
1371 :
1372 0 : Int_t recopass = fRecoPass;
1373 0 : if(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser;
1374 :
1375 0 : AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1376 : Bool_t found=kFALSE;
1377 : //
1378 : //set the new PID splines
1379 : //
1380 0 : if (fArrPidResponseMaster){
1381 : //for MC don't use period information
1382 : //if (fIsMC) period="[A-Z0-9]*";
1383 : //for MC use MC period information
1384 : //pattern for the default entry (valid for all particles)
1385 0 : TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1386 :
1387 : //find particle id and gain scenario
1388 0 : for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1389 : {
1390 : TObject *grAll=NULL;
1391 0 : TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1392 0 : gainScenario.ToUpper();
1393 : //loop over entries and filter them
1394 0 : for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1395 : {
1396 0 : TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1397 0 : if (responseFunction==NULL) continue;
1398 0 : TString responseName=responseFunction->GetName();
1399 :
1400 0 : if (!reg.MatchB(responseName)) continue;
1401 :
1402 0 : TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1403 : TObject* tmp=NULL;
1404 0 : tmp=arr->At(1); if (!tmp) continue;
1405 0 : TString particleName=tmp->GetName();
1406 0 : tmp=arr->At(3); if (!tmp) continue;
1407 0 : TString gainScenarioName=tmp->GetName();
1408 0 : delete arr;
1409 0 : if (particleName.IsNull()) continue;
1410 0 : if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1411 : else
1412 : {
1413 0 : for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1414 : {
1415 0 : TString particle=AliPID::ParticleName(ispec);
1416 0 : particle.ToUpper();
1417 : //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1418 0 : if ( particle == particleName && gainScenario == gainScenarioName )
1419 : {
1420 0 : fTPCResponse.SetResponseFunction( responseFunction,
1421 : (AliPID::EParticleType)ispec,
1422 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1423 0 : fTPCResponse.SetUseDatabase(kTRUE);
1424 0 : AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunction->GetName(),
1425 : AliTPCPIDResponse::GetChecksum((TSpline3*)responseFunction).Data()));
1426 : found=kTRUE;
1427 0 : break;
1428 : }
1429 0 : }
1430 : }
1431 0 : }
1432 :
1433 : // Retrieve responsefunction for pions - will (if available) be used for muons if there are no dedicated muon splines.
1434 : // For light nuclei, try to set the proton spline, if no dedicated splines are available.
1435 : // In both cases: Use default splines, if no dedicated splines and no pion/proton splines are available.
1436 0 : TObject* responseFunctionPion = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kPion,
1437 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1438 0 : TObject* responseFunctionProton = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kProton,
1439 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1440 :
1441 0 : for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1442 : {
1443 0 : if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1444 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
1445 : {
1446 0 : if (ispec == AliPID::kMuon) { // Muons
1447 0 : if (responseFunctionPion) {
1448 0 : fTPCResponse.SetResponseFunction( responseFunctionPion,
1449 : (AliPID::EParticleType)ispec,
1450 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1451 0 : fTPCResponse.SetUseDatabase(kTRUE);
1452 0 : AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionPion->GetName(),
1453 : AliTPCPIDResponse::GetChecksum((TSpline3*)responseFunctionPion).Data()));
1454 : found=kTRUE;
1455 0 : }
1456 0 : else if (grAll) {
1457 0 : fTPCResponse.SetResponseFunction( grAll,
1458 : (AliPID::EParticleType)ispec,
1459 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1460 0 : fTPCResponse.SetUseDatabase(kTRUE);
1461 0 : AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
1462 : AliTPCPIDResponse::GetChecksum((TSpline3*)grAll).Data()));
1463 : found=kTRUE;
1464 0 : }
1465 : //else
1466 : // AliError(Form("No splines found for muons (also no pion splines and no default splines) for gain scenario %d!", igainScenario));
1467 : }
1468 0 : else if (ispec >= AliPID::kSPECIES) { // Light nuclei
1469 0 : if (responseFunctionProton) {
1470 0 : fTPCResponse.SetResponseFunction( responseFunctionProton,
1471 : (AliPID::EParticleType)ispec,
1472 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1473 0 : fTPCResponse.SetUseDatabase(kTRUE);
1474 0 : AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionProton->GetName(),
1475 : AliTPCPIDResponse::GetChecksum((TSpline3*)responseFunctionProton).Data()));
1476 : found=kTRUE;
1477 0 : }
1478 0 : else if (grAll) {
1479 0 : fTPCResponse.SetResponseFunction( grAll,
1480 : (AliPID::EParticleType)ispec,
1481 : (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1482 0 : fTPCResponse.SetUseDatabase(kTRUE);
1483 0 : AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
1484 : AliTPCPIDResponse::GetChecksum((TSpline3*)grAll).Data()));
1485 : found=kTRUE;
1486 0 : }
1487 : //else
1488 : // AliError(Form("No splines found for species %d (also no proton splines and no default splines) for gain scenario %d!",
1489 : // ispec, igainScenario));
1490 : }
1491 : }
1492 : }
1493 0 : }
1494 0 : }
1495 0 : else AliInfo("no fArrPidResponseMaster");
1496 :
1497 0 : if (!found){
1498 0 : AliError("***** Risk for unreliable TPC PID detected: ********");
1499 0 : AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1500 : }
1501 :
1502 :
1503 : //
1504 : // Setup multiplicity correction (only used for non-pp collisions)
1505 : //
1506 :
1507 0 : const Bool_t isPP = (fBeamType.CompareTo("PP") == 0);
1508 :
1509 : // 2013 pPb data taking at low luminosity
1510 0 : const Bool_t isPPb2013LowLuminosity = period.Contains("LHC13B") || period.Contains("LHC13C") || period.Contains("LHC13D");
1511 : // PbPb 2010, period 10h.pass2
1512 : //TODO Needs further development const Bool_t is10hpass2 = period.Contains("LHC10H") && recopass == 2;
1513 :
1514 :
1515 : // In case of MC without(!) tune on data activated for the TPC, don't use the multiplicity correction for the moment
1516 0 : Bool_t isMCandNotTPCtuneOnData = fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC));
1517 :
1518 : // If correction is available, but disabled (highly NOT recommended!), print warning
1519 0 : if (!fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
1520 : //TODO: Needs further development if (is10hpass2 || isPPb2013LowLuminosity) {
1521 0 : if (isPPb2013LowLuminosity) {
1522 0 : AliWarning("Mulitplicity correction disabled, but correction parameters for this period exist. It is highly recommended to use enable the correction. Otherwise the splines might be off!");
1523 : }
1524 : }
1525 :
1526 0 : if (fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
1527 0 : AliInfo("Multiplicity correction enabled!");
1528 :
1529 : //TODO After testing, load parameters from outside
1530 : /*TODO no correction for MC
1531 : if (period.Contains("LHC11A10")) {//LHC11A10A
1532 : AliInfo("Using multiplicity correction parameters for 11a10!");
1533 : fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06);
1534 : fTPCResponse.SetParameterMultiplicityCorrection(1, -1.22123e-03);
1535 : fTPCResponse.SetParameterMultiplicityCorrection(2, 1.80220e-02);
1536 : fTPCResponse.SetParameterMultiplicityCorrection(3, 0.1);
1537 : fTPCResponse.SetParameterMultiplicityCorrection(4, 6.45306e-03);
1538 :
1539 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -2.85505e-07);
1540 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, -1.31911e-06);
1541 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1542 :
1543 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -4.29665e-05);
1544 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 1.37023e-02);
1545 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01);
1546 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02);
1547 : }
1548 0 : else*/ if (isPPb2013LowLuminosity) {// 2013 pPb data taking at low luminosity
1549 0 : AliInfo("Using multiplicity correction parameters for 13b.pass2 (at least also valid for 13{c,d} and pass 3)!");
1550 :
1551 0 : fTPCResponse.SetParameterMultiplicityCorrection(0, -5.906e-06);
1552 0 : fTPCResponse.SetParameterMultiplicityCorrection(1, -5.064e-04);
1553 0 : fTPCResponse.SetParameterMultiplicityCorrection(2, -3.521e-02);
1554 0 : fTPCResponse.SetParameterMultiplicityCorrection(3, 2.469e-02);
1555 0 : fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1556 :
1557 0 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.32e-06);
1558 0 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.177e-05);
1559 0 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1560 :
1561 0 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 0.);
1562 0 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 0.);
1563 0 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 0.);
1564 0 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 0.);
1565 :
1566 : /* Not too bad, but far from perfect in the details
1567 : fTPCResponse.SetParameterMultiplicityCorrection(0, -6.27187e-06);
1568 : fTPCResponse.SetParameterMultiplicityCorrection(1, -4.60649e-04);
1569 : fTPCResponse.SetParameterMultiplicityCorrection(2, -4.26450e-02);
1570 : fTPCResponse.SetParameterMultiplicityCorrection(3, 2.40590e-02);
1571 : fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1572 :
1573 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.338e-06);
1574 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.220e-05);
1575 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1576 :
1577 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 7.89237e-05);
1578 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, -1.30662e-02);
1579 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 8.91548e-01);
1580 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.47931e-02);
1581 : */
1582 : }
1583 : /*TODO: Needs further development
1584 : else if (is10hpass2) {
1585 : AliInfo("Using multiplicity correction parameters for 10h.pass2!");
1586 : fTPCResponse.SetParameterMultiplicityCorrection(0, 3.21636e-07);
1587 : fTPCResponse.SetParameterMultiplicityCorrection(1, -6.65876e-04);
1588 : fTPCResponse.SetParameterMultiplicityCorrection(2, 1.28786e-03);
1589 : fTPCResponse.SetParameterMultiplicityCorrection(3, 1.47677e-02);
1590 : fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1591 :
1592 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, 7.23591e-08);
1593 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 2.7469e-06);
1594 : fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1595 :
1596 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -1.22590e-05);
1597 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 6.88888e-03);
1598 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -3.20788e-01);
1599 : fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.07345e-02);
1600 : }
1601 : */
1602 : else {
1603 0 : AliError(Form("Multiplicity correction is enabled, but no multiplicity correction parameters have been found for period %s.pass%d -> Mulitplicity correction DISABLED!", period.Data(), recopass));
1604 0 : fUseTPCMultiplicityCorrection = kFALSE;
1605 0 : fTPCResponse.ResetMultiplicityCorrectionFunctions();
1606 : }
1607 : }
1608 : else {
1609 : // Just set parameters such that overall correction factor is 1, i.e. no correction.
1610 : // This is just a reasonable choice for the parameters for safety reasons. Disabling
1611 : // the multiplicity correction will anyhow skip the calculation of the corresponding
1612 : // correction factor inside THIS class. Nevertheless, experts can access the TPCPIDResponse
1613 : // directly and use it for calculations - which should still give valid results, even if
1614 : // the multiplicity correction is explicitely enabled in such expert calls.
1615 :
1616 0 : TString reasonForDisabling = "requested by user";
1617 0 : if (fUseTPCMultiplicityCorrection) {
1618 0 : if (isPP)
1619 0 : reasonForDisabling = "pp collisions";
1620 : else
1621 0 : reasonForDisabling = "MC w/o tune on data";
1622 : }
1623 :
1624 0 : AliInfo(Form("Multiplicity correction %sdisabled (%s)!", fUseTPCMultiplicityCorrection ? "automatically " : "",
1625 : reasonForDisabling.Data()));
1626 :
1627 0 : fUseTPCMultiplicityCorrection = kFALSE;
1628 0 : fTPCResponse.ResetMultiplicityCorrectionFunctions();
1629 0 : }
1630 :
1631 0 : if (fUseTPCMultiplicityCorrection) {
1632 0 : for (Int_t i = 0; i <= 4 + 1; i++) {
1633 0 : AliInfo(Form("parMultCorr: %d, %e", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i)));
1634 : }
1635 0 : for (Int_t j = 0; j <= 2 + 1; j++) {
1636 0 : AliInfo(Form("parMultCorrTanTheta: %d, %e", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j)));
1637 : }
1638 0 : for (Int_t j = 0; j <= 3 + 1; j++) {
1639 0 : AliInfo(Form("parMultSigmaCorr: %d, %e", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j)));
1640 : }
1641 0 : }
1642 :
1643 : //
1644 : // Setup old resolution parametrisation
1645 : //
1646 :
1647 : //default
1648 0 : fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1649 :
1650 0 : if (fRun>=122195){ //LHC10d
1651 0 : fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1652 : }
1653 :
1654 0 : if (fRun>=170719){ // LHC12a
1655 0 : fTPCResponse.SetSigma(2.95714e-03, 1.01953e+05);
1656 : }
1657 :
1658 0 : if (fRun>=177312){ // LHC12b
1659 0 : fTPCResponse.SetSigma(3.74633e-03, 7.11829e+04 );
1660 : }
1661 :
1662 0 : if (fRun>=186346){ // LHC12e
1663 0 : fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1664 : }
1665 :
1666 0 : if (fArrPidResponseMaster)
1667 0 : fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
1668 :
1669 0 : if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s (MD5(corr function) = %s)",
1670 : fResolutionCorrection->GetName(), AliTPCPIDResponse::GetChecksum(fResolutionCorrection).Data()));
1671 :
1672 : //read in the voltage map
1673 : TVectorF* gsm = 0x0;
1674 0 : if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
1675 0 : if (gsm)
1676 : {
1677 0 : fTPCResponse.SetVoltageMap(*gsm);
1678 0 : TString vals;
1679 0 : AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1680 0 : vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1681 0 : AliInfo(vals.Data());
1682 0 : vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1683 0 : AliInfo(vals.Data());
1684 0 : vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1685 0 : AliInfo(vals.Data());
1686 0 : vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1687 0 : AliInfo(vals.Data());
1688 0 : }
1689 0 : else AliInfo("no voltage map, ideal default assumed");
1690 0 : AliInfo("------------------------------------------------------------------------------------------");
1691 0 : }
1692 :
1693 : //______________________________________________________________________________
1694 : void AliPIDResponse::SetTRDPidResponseMaster()
1695 : {
1696 : //
1697 : // Load the TRD pid params and references from the OADB
1698 : //
1699 0 : if(fTRDPIDResponseObject) return;
1700 0 : AliOADBContainer contParams("contParams");
1701 :
1702 0 : Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1703 0 : if(statusResponse){
1704 0 : AliError("Failed initializing PID Response Object from OADB");
1705 : } else {
1706 0 : AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1707 0 : fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1708 0 : if(!fTRDPIDResponseObject){
1709 0 : AliError(Form("TRD Response not found in run %d", fRun));
1710 : }
1711 : }
1712 0 : }
1713 :
1714 : //______________________________________________________________________________
1715 : void AliPIDResponse::InitializeTRDResponse(){
1716 : //
1717 : // Set PID Params and references to the TRD PID response
1718 : //
1719 0 : fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
1720 0 : fTRDResponse.SetdEdxParams(fTRDdEdxParams);
1721 0 : }
1722 :
1723 : //______________________________________________________________________________
1724 : void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1725 :
1726 0 : if(fLHCperiod.Contains("LHC10D") || fLHCperiod.Contains("LHC10E")){
1727 : // backward compatibility for setting with 8 slices
1728 0 : TRDslicesForPID[0] = 0;
1729 0 : TRDslicesForPID[1] = 7;
1730 0 : }
1731 : else{
1732 0 : if(method==AliTRDPIDResponse::kLQ1D){
1733 0 : TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1734 0 : TRDslicesForPID[1] = 0;
1735 0 : }
1736 0 : if((method==AliTRDPIDResponse::kLQ2D)||(method==AliTRDPIDResponse::kLQ3D)||(method==AliTRDPIDResponse::kLQ7D)){
1737 0 : TRDslicesForPID[0] = 1;
1738 0 : TRDslicesForPID[1] = 7;
1739 0 : }
1740 : }
1741 0 : AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
1742 0 : }
1743 : //______________________________________________________________________________
1744 : void AliPIDResponse::SetTRDdEdxParams()
1745 : {
1746 0 : if(fTRDdEdxParams) return;
1747 :
1748 0 : const TString containerName = "TRDdEdxParamsContainer";
1749 0 : AliOADBContainer cont(containerName.Data());
1750 :
1751 0 : const TString filePathNamePackage=Form("%s/COMMON/PID/data/TRDdEdxParams.root", fOADBPath.Data());
1752 :
1753 0 : const Int_t statusCont = cont.InitFromFile(filePathNamePackage.Data(), cont.GetName());
1754 0 : if (statusCont){
1755 0 : AliFatal("Failed initializing settings from OADB");
1756 : }
1757 : else{
1758 0 : AliInfo(Form("Loading %s from %s\n", cont.GetName(), filePathNamePackage.Data()));
1759 :
1760 0 : fTRDdEdxParams = (AliTRDdEdxParams*)(cont.GetObject(fRun, "default"));
1761 :
1762 0 : if(!fTRDdEdxParams){
1763 0 : AliError(Form("TRD dEdx Params default not found"));
1764 : }
1765 : }
1766 0 : }
1767 :
1768 : //______________________________________________________________________________
1769 : void AliPIDResponse::SetTRDEtaMaps()
1770 : {
1771 : //
1772 : // Load the TRD eta correction map from the OADB
1773 : //
1774 :
1775 0 : if (fIsMC) fUseTRDEtaCorrection = kFALSE;
1776 0 : if (fUseTRDEtaCorrection == kFALSE) {
1777 : // fTRDResponse.SetEtaCorrMap(0,0x0);
1778 0 : AliInfo("Request to disable TRD eta correction -> Eta correction has been disabled");
1779 0 : return;
1780 : }
1781 : TH2D* etaMap[1];
1782 : etaMap[0] = 0x0;
1783 :
1784 :
1785 0 : const TString containerName = "TRDEtaCorrectionMap";
1786 0 : AliOADBContainer cont(containerName.Data());
1787 :
1788 0 : const TString filePathNamePackage=Form("%s/COMMON/PID/data/TRDdEdxEtaCorrectionParams.root", fOADBPath.Data());
1789 :
1790 0 : const Int_t statusCont = cont.InitFromFile(filePathNamePackage.Data(), cont.GetName());
1791 0 : if (statusCont){
1792 0 : AliFatal("Failed initializing TRD Eta Correction settings from OADB");
1793 0 : return;
1794 : }
1795 : else{
1796 0 : AliInfo(Form("Loading %s from %s\n", cont.GetName(), filePathNamePackage.Data()));
1797 :
1798 0 : TObject* etaarray=(TObject*)cont.GetObject(fRun);
1799 :
1800 0 : if(etaarray){
1801 0 : etaMap[0] = (TH2D *)etaarray->FindObject("TRDEtaMap");
1802 0 : fTRDResponse.SetEtaCorrMap(0,etaMap[0]);
1803 : }
1804 : else{
1805 0 : AliError(Form("TRD Eta Correction Params not found"));
1806 0 : fUseTRDEtaCorrection = kFALSE;
1807 0 : return;
1808 : //fTRDResponse.SetEtaCorrMap(0,0x0);
1809 : }
1810 :
1811 :
1812 :
1813 0 : if (!etaMap[0]) {
1814 0 : AliError(Form("TRD Eta Correction Params not found"));
1815 0 : fUseTRDEtaCorrection = kFALSE;
1816 0 : return;
1817 : //fTRDResponse.SetEtaCorrMap(0,0x0);
1818 : }
1819 :
1820 :
1821 0 : }
1822 :
1823 :
1824 0 : }
1825 :
1826 :
1827 : //______________________________________________________________________________
1828 : void AliPIDResponse::SetTOFPidResponseMaster()
1829 : {
1830 : //
1831 : // Load the TOF pid params from the OADB
1832 : //
1833 :
1834 0 : if (fTOFPIDParams) delete fTOFPIDParams;
1835 0 : fTOFPIDParams=NULL;
1836 :
1837 0 : TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1838 0 : if (oadbf && oadbf->IsOpen()) {
1839 0 : Int_t recoPass = fRecoPass;
1840 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) recoPass = fRecoPassUser;
1841 0 : TString passName = Form("pass%d",recoPass);
1842 0 : AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root for run %d (pass: %s)", fOADBPath.Data(),fRun,passName.Data()));
1843 0 : AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
1844 0 : if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams",passName));
1845 0 : oadbf->Close();
1846 0 : delete oadbc;
1847 0 : } else {
1848 0 : AliError(Form("TOF OADB file not found!!! %s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
1849 : }
1850 0 : delete oadbf;
1851 :
1852 0 : if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1853 :
1854 0 : if (TString(fTOFPIDParams->GetOADBentryTag()) == "default") {
1855 0 : AliWarning("******* Risk for unreliable TOF PID detected *********");
1856 0 : if (!fIsMC && fRecoPass<=0) {
1857 0 : AliWarningF(" Invalid reco pass for data (%d) was detected", fRecoPass);
1858 0 : }
1859 0 : AliWarning(" The default object was loaded");
1860 0 : }
1861 0 : }
1862 :
1863 : //______________________________________________________________________________
1864 : void AliPIDResponse::InitializeTOFResponse()
1865 : {
1866 : //
1867 : // Set PID Params to the TOF PID response
1868 : //
1869 0 : AliInfo("---------------------------- TOF Response Configuration ----------------------------");
1870 0 : AliInfo(Form("TOF PID Params loaded from OADB [entryTag: %s]",fTOFPIDParams->GetOADBentryTag()));
1871 0 : AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1872 0 : AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1873 0 : AliInfo(Form(" TOF res. mom. params: %6.3f %6.3f %6.3f %6.3f",
1874 : fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1875 0 : AliInfo(Form(" Start Time Offset %6.2f ps",fTOFPIDParams->GetTOFtimeOffset()));
1876 0 : AliInfo(Form(" Fraction of tracks within gaussian behaviour: %6.4f",fTOFPIDParams->GetTOFtail()));
1877 :
1878 0 : if (fIsMC) {
1879 0 : AliInfo("MC data:");
1880 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) ) {
1881 0 : AliInfo(Form(" MC: TuneOnData option enabled for TOF data (target data reco pass is %d)",fRecoPassUser));
1882 0 : } else AliInfo(" MC: TuneOnData option NOT enabled for TOF data");
1883 : }
1884 0 : AliInfo(Form(" MC: Fraction of tracks (percentage) to cut to fit matching in data: %6.2f%%",fTOFPIDParams->GetTOFmatchingLossMC()));
1885 0 : AliInfo(Form(" MC: Fraction of random hits (percentage) to add to fit mismatch in data: %6.2f%%",fTOFPIDParams->GetTOFadditionalMismForMC()));
1886 :
1887 0 : for (Int_t i=0;i<4;i++) {
1888 0 : fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1889 : }
1890 0 : fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1891 :
1892 0 : AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1893 0 : Float_t t0Spread[4];
1894 0 : for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1895 0 : 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]));
1896 0 : Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1897 0 : Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1898 0 : if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1899 0 : fResT0AC=t0Spread[3];
1900 0 : fResT0A=TMath::Sqrt(a);
1901 0 : fResT0C=TMath::Sqrt(c);
1902 0 : } else {
1903 0 : AliInfo(" TZERO spreads not present or inconsistent, loading default");
1904 0 : fResT0A=75.;
1905 0 : fResT0C=65.;
1906 0 : fResT0AC=55.;
1907 : }
1908 0 : AliInfo(Form(" TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1909 0 : AliInfo("------------------------------------------------------------------------------------");
1910 :
1911 0 : }
1912 :
1913 : //______________________________________________________________________________
1914 : void AliPIDResponse::SetHMPIDPidResponseMaster()
1915 : {
1916 : //
1917 : // Load the HMPID pid params from the OADB
1918 : //
1919 :
1920 0 : if (fHMPIDPIDParams) delete fHMPIDPIDParams;
1921 0 : fHMPIDPIDParams=NULL;
1922 :
1923 : TFile *oadbf;
1924 0 : if(!fIsMC) oadbf = new TFile(Form("%s/COMMON/PID/data/HMPIDPIDParams.root",fOADBPath.Data()));
1925 0 : else oadbf = new TFile(Form("%s/COMMON/PID/MC/HMPIDPIDParams.root",fOADBPath.Data()));
1926 0 : if (oadbf && oadbf->IsOpen()) {
1927 0 : AliInfo(Form("Loading HMPID Params from %s/COMMON/PID/data/HMPIDPIDParams.root", fOADBPath.Data()));
1928 0 : AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("HMPoadb");
1929 0 : if (oadbc) fHMPIDPIDParams = dynamic_cast<AliHMPIDPIDParams *>(oadbc->GetObject(fRun,"HMPparams"));
1930 0 : oadbf->Close();
1931 0 : delete oadbc;
1932 0 : }
1933 0 : delete oadbf;
1934 :
1935 0 : if (!fHMPIDPIDParams) AliFatal("HMPIDPIDParams could not be retrieved");
1936 0 : }
1937 :
1938 : //______________________________________________________________________________
1939 : void AliPIDResponse::InitializeHMPIDResponse(){
1940 : //
1941 : // Set PID Params to the HMPID PID response
1942 : //
1943 :
1944 0 : fHMPIDResponse.SetRefIndexArray(fHMPIDPIDParams->GetHMPIDrefIndex());
1945 0 : }
1946 :
1947 : //______________________________________________________________________________
1948 : Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1949 : // old function for compatibility
1950 0 : Int_t ntracklets=0;
1951 0 : return IdentifiedAsElectronTRD(vtrack,ntracklets,efficiencyLevel,centrality,PIDmethod);
1952 0 : }
1953 :
1954 : //______________________________________________________________________________
1955 : Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Int_t &ntracklets,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1956 : //
1957 : // Check whether track is identified as electron under a given electron efficiency hypothesis
1958 : //
1959 : // ntracklets is the number of tracklets that has been used to calculate the PID signal
1960 :
1961 0 : Double_t probs[AliPID::kSPECIES];
1962 :
1963 0 : ntracklets =CalculateTRDResponse(vtrack,probs,PIDmethod);
1964 :
1965 : // Take mean of the TRD momenta in the given tracklets
1966 0 : Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1967 : Int_t nmomenta = 0;
1968 0 : for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1969 0 : if(vtrack->GetTRDmomentum(iPl) > 0.){
1970 0 : trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
1971 0 : }
1972 : }
1973 0 : p = TMath::Mean(nmomenta, trdmomenta);
1974 :
1975 0 : return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
1976 0 : }
1977 :
1978 : //______________________________________________________________________________
1979 : void AliPIDResponse::SetEMCALPidResponseMaster()
1980 : {
1981 : //
1982 : // Load the EMCAL pid response functions from the OADB
1983 : //
1984 : TObjArray* fEMCALPIDParamsRun = NULL;
1985 : TObjArray* fEMCALPIDParamsPass = NULL;
1986 :
1987 0 : if(fEMCALPIDParams) return;
1988 0 : AliOADBContainer contParams("contParams");
1989 :
1990 0 : Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1991 0 : if(statusPars){
1992 0 : AliError("Failed initializing PID Params from OADB");
1993 : }
1994 : else {
1995 0 : AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1996 :
1997 0 : fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1998 0 : if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1999 0 : if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
2000 :
2001 0 : if(!fEMCALPIDParams){
2002 0 : AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
2003 0 : AliInfo("Will take the standard LHC11d instead ...");
2004 :
2005 0 : fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
2006 0 : if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
2007 0 : if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
2008 :
2009 0 : if(!fEMCALPIDParams){
2010 0 : AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
2011 : }
2012 : }
2013 : }
2014 0 : }
2015 :
2016 : //______________________________________________________________________________
2017 : void AliPIDResponse::InitializeEMCALResponse(){
2018 : //
2019 : // Set PID Params to the EMCAL PID response
2020 : //
2021 0 : fEMCALResponse.SetPIDParams(fEMCALPIDParams);
2022 :
2023 0 : }
2024 :
2025 : //______________________________________________________________________________
2026 : void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
2027 : {
2028 : //
2029 : // create detector PID information and setup the transient pointer in the track
2030 : //
2031 :
2032 : // check if detector number is inside accepted range
2033 0 : if (detector == kNdetectors) return;
2034 :
2035 : // get detector pid
2036 0 : AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
2037 0 : if (!detPID) {
2038 0 : detPID=new AliDetectorPID;
2039 0 : (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
2040 0 : }
2041 :
2042 : //check if values exist
2043 0 : if (detPID->HasRawProbability(detector) && detPID->HasNumberOfSigmas(detector)) return;
2044 :
2045 : //TODO: which particles to include? See also the loops below...
2046 0 : Double_t values[AliPID::kSPECIESC]={0};
2047 :
2048 : //probabilities
2049 0 : EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
2050 0 : detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
2051 :
2052 : //nsigmas
2053 0 : for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
2054 0 : values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
2055 : // the pid status is the same for probabilities and nSigmas, so it is
2056 : // fine to use the one from the probabilities also here
2057 0 : detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC, status);
2058 :
2059 0 : }
2060 :
2061 : //______________________________________________________________________________
2062 : void AliPIDResponse::FillTrackDetectorPID()
2063 : {
2064 : //
2065 : // create detector PID information and setup the transient pointer in the track
2066 : //
2067 :
2068 0 : if (!fCurrentEvent) return;
2069 :
2070 0 : for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
2071 0 : AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
2072 0 : if (!track) continue;
2073 :
2074 0 : for (Int_t idet=0; idet<kNdetectors; ++idet){
2075 0 : FillTrackDetectorPID(track, (EDetector)idet);
2076 : }
2077 0 : }
2078 0 : }
2079 :
2080 : //______________________________________________________________________________
2081 : void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
2082 : //
2083 : // Set TOF response function
2084 : // Input option for event_time used
2085 : //
2086 :
2087 : Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
2088 24 : if(t0spread < 10) t0spread = 80;
2089 :
2090 : // T0-FILL and T0-TO offset (because of TOF misallignment
2091 : Float_t starttimeoffset = 0;
2092 8 : if(fTOFPIDParams && !(fIsMC)) starttimeoffset=fTOFPIDParams->GetTOFtimeOffset();
2093 8 : if(fTOFPIDParams){
2094 0 : fTOFtail = fTOFPIDParams->GetTOFtail();
2095 0 : GetTOFResponse().SetTOFtail(fTOFtail);
2096 0 : }
2097 :
2098 : // T0 from TOF algorithm
2099 : Bool_t flagT0TOF=kFALSE;
2100 : Bool_t flagT0T0=kFALSE;
2101 8 : Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
2102 8 : Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
2103 8 : Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
2104 :
2105 : // T0-TOF arrays
2106 8 : Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
2107 8 : Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
2108 176 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2109 80 : estimatedT0event[i]=0.0;
2110 80 : estimatedT0resolution[i]=0.0;
2111 80 : startTimeMask[i] = 0;
2112 : }
2113 :
2114 8 : Float_t resT0A=fResT0A;
2115 8 : Float_t resT0C=fResT0C;
2116 8 : Float_t resT0AC=fResT0AC;
2117 8 : if(vevent->GetT0TOF()){ // check if T0 detector information is available
2118 : flagT0T0=kTRUE;
2119 8 : }
2120 :
2121 :
2122 8 : AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
2123 :
2124 8 : if (tofHeader) { // read global info and T0-TOF
2125 : // fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution()); // read from OADB in the initialization
2126 8 : t0spread = tofHeader->GetT0spread(); // read t0 sprad
2127 8 : if(t0spread < 10) t0spread = 80;
2128 :
2129 : flagT0TOF=kTRUE;
2130 176 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
2131 80 : startTime[i]=tofHeader->GetDefaultEventTimeVal();
2132 80 : startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
2133 80 : if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
2134 :
2135 80 : if(startTimeRes[i] > t0spread - 10 && TMath::Abs(startTime[i]) < 0.001) startTime[i] = -starttimeoffset; // apply offset for T0-fill
2136 : }
2137 :
2138 8 : TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
2139 8 : TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
2140 8 : TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
2141 102 : for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
2142 43 : Int_t icurrent = (Int_t)ibin->GetAt(j);
2143 43 : startTime[icurrent]=t0Bin->GetAt(j);
2144 43 : startTimeRes[icurrent]=t0ResBin->GetAt(j);
2145 43 : if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
2146 43 : if(startTimeRes[icurrent] > t0spread - 10 && TMath::Abs(startTime[icurrent]) < 0.001) startTime[icurrent] = -starttimeoffset; // apply offset for T0-fill
2147 : }
2148 8 : }
2149 :
2150 : // for cut of 3 sigma on t0 spread
2151 8 : Float_t t0cut = 3 * t0spread;
2152 8 : if(t0cut < 500) t0cut = 500;
2153 :
2154 8 : if(option == kFILL_T0){ // T0-FILL is used
2155 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2156 0 : estimatedT0event[i]=0.0-starttimeoffset;
2157 0 : estimatedT0resolution[i]=t0spread;
2158 : }
2159 0 : fTOFResponse.SetT0event(estimatedT0event);
2160 0 : fTOFResponse.SetT0resolution(estimatedT0resolution);
2161 0 : }
2162 :
2163 8 : if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
2164 8 : if(flagT0TOF){
2165 8 : fTOFResponse.SetT0event(startTime);
2166 8 : fTOFResponse.SetT0resolution(startTimeRes);
2167 176 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2168 160 : if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
2169 80 : fTOFResponse.SetT0binMask(i,startTimeMask[i]);
2170 : }
2171 8 : }
2172 : else{
2173 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2174 0 : estimatedT0event[i]=0.0-starttimeoffset;
2175 0 : estimatedT0resolution[i]=t0spread;
2176 0 : fTOFResponse.SetT0binMask(i,startTimeMask[i]);
2177 : }
2178 0 : fTOFResponse.SetT0event(estimatedT0event);
2179 0 : fTOFResponse.SetT0resolution(estimatedT0resolution);
2180 : }
2181 : }
2182 0 : else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
2183 : Float_t t0AC=-10000;
2184 : Float_t t0A=-10000;
2185 : Float_t t0C=-10000;
2186 0 : if(flagT0T0){
2187 0 : t0A= vevent->GetT0TOF()[1] - starttimeoffset;
2188 0 : t0C= vevent->GetT0TOF()[2] - starttimeoffset;
2189 0 : t0AC= vevent->GetT0TOF()[0] - starttimeoffset;
2190 : //t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
2191 : // resT0AC= 1./TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
2192 : // t0AC *= resT0AC*resT0AC;
2193 0 : }
2194 :
2195 : Float_t t0t0Best = 0;
2196 : Float_t t0t0BestRes = 9999;
2197 : Int_t t0used=0;
2198 0 : if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
2199 : t0t0Best = t0AC;
2200 : t0t0BestRes = resT0AC;
2201 : t0used=6;
2202 0 : }
2203 0 : else if(TMath::Abs(t0C) < t0cut){
2204 : t0t0Best = t0C;
2205 : t0t0BestRes = resT0C;
2206 : t0used=4;
2207 0 : }
2208 0 : else if(TMath::Abs(t0A) < t0cut){
2209 : t0t0Best = t0A;
2210 : t0t0BestRes = resT0A;
2211 : t0used=2;
2212 0 : }
2213 :
2214 0 : if(flagT0TOF){ // if T0-TOF info is available
2215 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2216 0 : if(t0t0BestRes < 999){
2217 0 : if(startTimeRes[i] < t0spread){
2218 0 : Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
2219 0 : Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
2220 0 : estimatedT0event[i]=t0best / wtot;
2221 0 : estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
2222 0 : startTimeMask[i] = t0used+1;
2223 0 : }
2224 : else {
2225 0 : estimatedT0event[i]=t0t0Best;
2226 0 : estimatedT0resolution[i]=t0t0BestRes;
2227 0 : startTimeMask[i] = t0used;
2228 : }
2229 : }
2230 : else{
2231 0 : estimatedT0event[i]=startTime[i];
2232 0 : estimatedT0resolution[i]=startTimeRes[i];
2233 0 : if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
2234 : }
2235 0 : fTOFResponse.SetT0binMask(i,startTimeMask[i]);
2236 : }
2237 0 : fTOFResponse.SetT0event(estimatedT0event);
2238 0 : fTOFResponse.SetT0resolution(estimatedT0resolution);
2239 0 : }
2240 : else{ // if no T0-TOF info is available
2241 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2242 0 : fTOFResponse.SetT0binMask(i,t0used);
2243 0 : if(t0t0BestRes < 999){
2244 0 : estimatedT0event[i]=t0t0Best;
2245 0 : estimatedT0resolution[i]=t0t0BestRes;
2246 0 : }
2247 : else{
2248 0 : estimatedT0event[i]=0.0-starttimeoffset;
2249 0 : estimatedT0resolution[i]=t0spread;
2250 : }
2251 : }
2252 0 : fTOFResponse.SetT0event(estimatedT0event);
2253 0 : fTOFResponse.SetT0resolution(estimatedT0resolution);
2254 : }
2255 0 : }
2256 :
2257 0 : else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
2258 : Float_t t0AC=-10000;
2259 : Float_t t0A=-10000;
2260 : Float_t t0C=-10000;
2261 0 : if(flagT0T0){
2262 0 : t0A= vevent->GetT0TOF()[1] - starttimeoffset;
2263 0 : t0C= vevent->GetT0TOF()[2] - starttimeoffset;
2264 0 : t0AC= vevent->GetT0TOF()[0] - starttimeoffset;
2265 : // t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
2266 : // resT0AC= 1./TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
2267 : // t0AC *= resT0AC*resT0AC;
2268 0 : }
2269 :
2270 0 : if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
2271 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2272 0 : estimatedT0event[i]=t0AC;
2273 0 : estimatedT0resolution[i]=resT0AC;
2274 0 : fTOFResponse.SetT0binMask(i,6);
2275 : }
2276 0 : }
2277 0 : else if(TMath::Abs(t0C) < t0cut){
2278 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2279 0 : estimatedT0event[i]=t0C;
2280 0 : estimatedT0resolution[i]=resT0C;
2281 0 : fTOFResponse.SetT0binMask(i,4);
2282 : }
2283 0 : }
2284 0 : else if(TMath::Abs(t0A) < t0cut){
2285 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2286 0 : estimatedT0event[i]=t0A;
2287 0 : estimatedT0resolution[i]=resT0A;
2288 0 : fTOFResponse.SetT0binMask(i,2);
2289 : }
2290 0 : }
2291 : else{
2292 0 : for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2293 0 : estimatedT0event[i]= 0.0 - starttimeoffset;
2294 0 : estimatedT0resolution[i]=t0spread;
2295 0 : fTOFResponse.SetT0binMask(i,0);
2296 : }
2297 : }
2298 0 : fTOFResponse.SetT0event(estimatedT0event);
2299 0 : fTOFResponse.SetT0resolution(estimatedT0resolution);
2300 0 : }
2301 :
2302 16 : delete [] startTime;
2303 16 : delete [] startTimeRes;
2304 16 : delete [] startTimeMask;
2305 16 : delete [] estimatedT0event;
2306 16 : delete [] estimatedT0resolution;
2307 8 : }
2308 :
2309 : //______________________________________________________________________________
2310 : // private non cached versions of the PID calculation
2311 : //
2312 :
2313 :
2314 : //______________________________________________________________________________
2315 : Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
2316 : {
2317 : //
2318 : // NumberOfSigmas for 'detCode'
2319 : //
2320 :
2321 0 : const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
2322 :
2323 0 : switch (detector){
2324 0 : case kITS: return GetNumberOfSigmasITS(track, type); break;
2325 0 : case kTPC: return GetNumberOfSigmasTPC(track, type); break;
2326 0 : case kTRD: return GetNumberOfSigmasTRD(track, type); break;
2327 0 : case kTOF: return GetNumberOfSigmasTOF(track, type); break;
2328 0 : case kHMPID: return GetNumberOfSigmasHMPID(track, type); break;
2329 0 : case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
2330 0 : default: return -999.;
2331 : }
2332 :
2333 : return -999.;
2334 0 : }
2335 :
2336 : //______________________________________________________________________________
2337 : Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
2338 : {
2339 : //
2340 : // Calculate the number of sigmas in the ITS
2341 : //
2342 :
2343 0 : AliVTrack *track=(AliVTrack*)vtrack;
2344 :
2345 0 : const EDetPidStatus pidStatus=GetITSPIDStatus(track);
2346 0 : if (pidStatus!=kDetPidOk) return -999.;
2347 :
2348 : // the following call is needed in order to fill the transient data member
2349 : // fITSsignalTuned which is used in the ITSPIDResponse to judge
2350 : // if using tuned on data
2351 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetITS) == kDetITS))
2352 0 : GetITSsignalTunedOnData(track);
2353 :
2354 0 : return fITSResponse.GetNumberOfSigmas(track,type);
2355 0 : }
2356 :
2357 : //______________________________________________________________________________
2358 : Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
2359 : {
2360 : //
2361 : // Calculate the number of sigmas in the TPC
2362 : //
2363 :
2364 0 : AliVTrack *track=(AliVTrack*)vtrack;
2365 :
2366 0 : const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
2367 0 : if (pidStatus==kDetNoSignal) return -999.;
2368 :
2369 : // the following call is needed in order to fill the transient data member
2370 : // fTPCsignalTuned which is used in the TPCPIDResponse to judge
2371 : // if using tuned on data
2372 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
2373 0 : GetTPCsignalTunedOnData(track);
2374 :
2375 0 : return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
2376 0 : }
2377 :
2378 : //______________________________________________________________________________
2379 : Float_t AliPIDResponse::GetNumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
2380 : {
2381 : //
2382 : // Calculate the number of sigmas in the TRD
2383 : //
2384 :
2385 0 : AliVTrack *track=(AliVTrack*)vtrack;
2386 :
2387 0 : const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
2388 0 : if (pidStatus!=kDetPidOk) return -999.;
2389 :
2390 0 : return fTRDResponse.GetNumberOfSigmas(track,type, fUseTRDEtaCorrection);
2391 0 : }
2392 :
2393 : //______________________________________________________________________________
2394 : Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
2395 : {
2396 : //
2397 : // Calculate the number of sigmas in the TOF
2398 : //
2399 :
2400 0 : AliVTrack *track=(AliVTrack*)vtrack;
2401 :
2402 0 : const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
2403 0 : if (pidStatus!=kDetPidOk) return -999.;
2404 :
2405 0 : return GetNumberOfSigmasTOFold(vtrack, type);
2406 0 : }
2407 : //______________________________________________________________________________
2408 :
2409 : Float_t AliPIDResponse::GetNumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
2410 : {
2411 : //
2412 : // Calculate the number of sigmas in the HMPID
2413 : //
2414 0 : AliVTrack *track=(AliVTrack*)vtrack;
2415 :
2416 0 : const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2417 0 : if (pidStatus!=kDetPidOk) return -999.;
2418 :
2419 0 : return fHMPIDResponse.GetNumberOfSigmas(track, type);
2420 0 : }
2421 :
2422 : //______________________________________________________________________________
2423 : Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
2424 : {
2425 : //
2426 : // Calculate the number of sigmas in the EMCAL
2427 : //
2428 :
2429 0 : AliVTrack *track=(AliVTrack*)vtrack;
2430 :
2431 0 : const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
2432 0 : if (pidStatus!=kDetPidOk) return -999.;
2433 :
2434 0 : const Int_t nMatchClus = track->GetEMCALcluster();
2435 0 : AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2436 :
2437 0 : const Double_t mom = track->P();
2438 0 : const Double_t pt = track->Pt();
2439 0 : const Int_t charge = track->Charge();
2440 0 : const Double_t fClsE = matchedClus->E();
2441 0 : const Double_t EovP = fClsE/mom;
2442 :
2443 0 : return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
2444 0 : }
2445 :
2446 : //______________________________________________________________________________
2447 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaITS(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2448 : {
2449 : //
2450 : // Signal minus expected Signal for ITS
2451 : //
2452 0 : AliVTrack *track=(AliVTrack*)vtrack;
2453 :
2454 : // the following call is needed in order to fill the transient data member
2455 : // fITSsignalTuned which is used in the ITSPIDResponse to judge
2456 : // if using tuned on data
2457 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetITS) == kDetITS))
2458 0 : GetITSsignalTunedOnData(track);
2459 :
2460 0 : val=fITSResponse.GetSignalDelta(track,type,ratio);
2461 :
2462 0 : return GetITSPIDStatus(track);
2463 : }
2464 :
2465 : //______________________________________________________________________________
2466 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2467 : {
2468 : //
2469 : // Signal minus expected Signal for TPC
2470 : //
2471 0 : AliVTrack *track=(AliVTrack*)vtrack;
2472 :
2473 : // the following call is needed in order to fill the transient data member
2474 : // fTPCsignalTuned which is used in the TPCPIDResponse to judge
2475 : // if using tuned on data
2476 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
2477 0 : GetTPCsignalTunedOnData(track);
2478 :
2479 0 : val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, ratio);
2480 :
2481 0 : return GetTPCPIDStatus(track);
2482 : }
2483 :
2484 : //______________________________________________________________________________
2485 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTRD(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2486 : {
2487 : //
2488 : // Signal minus expected Signal for TRD
2489 : //
2490 0 : AliVTrack *track=(AliVTrack*)vtrack;
2491 0 : val=fTRDResponse.GetSignalDelta(track,type,ratio,fUseTRDEtaCorrection);
2492 :
2493 0 : return GetTRDPIDStatus(track);
2494 : }
2495 :
2496 : //______________________________________________________________________________
2497 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2498 : {
2499 : //
2500 : // Signal minus expected Signal for TOF
2501 : //
2502 0 : AliVTrack *track=(AliVTrack*)vtrack;
2503 0 : val=GetSignalDeltaTOFold(track, type, ratio);
2504 :
2505 0 : return GetTOFPIDStatus(track);
2506 : }
2507 :
2508 : //______________________________________________________________________________
2509 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2510 : {
2511 : //
2512 : // Signal minus expected Signal for HMPID
2513 : //
2514 0 : AliVTrack *track=(AliVTrack*)vtrack;
2515 0 : val=fHMPIDResponse.GetSignalDelta(track, type, ratio);
2516 :
2517 0 : return GetHMPIDPIDStatus(track);
2518 : }
2519 :
2520 : //______________________________________________________________________________
2521 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2522 : {
2523 : //
2524 : // Compute PID response of 'detCode'
2525 : //
2526 :
2527 576 : switch (detCode){
2528 0 : case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
2529 288 : case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
2530 0 : case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
2531 0 : case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
2532 0 : case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
2533 0 : case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
2534 0 : case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
2535 0 : default: return kDetNoSignal;
2536 : }
2537 288 : }
2538 :
2539 : //______________________________________________________________________________
2540 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2541 : {
2542 : //
2543 : // Compute PID response for the ITS
2544 : //
2545 :
2546 : // set flat distribution (no decision)
2547 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2548 :
2549 0 : const EDetPidStatus pidStatus=GetITSPIDStatus(track);
2550 0 : if (pidStatus!=kDetPidOk) return pidStatus;
2551 :
2552 0 : if (track->GetDetectorPID()){
2553 0 : return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
2554 : }
2555 :
2556 : //check for ITS standalone tracks
2557 : Bool_t isSA=kTRUE;
2558 0 : if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
2559 :
2560 0 : Double_t mom=track->P();
2561 0 : Double_t dedx=track->GetITSsignal();
2562 0 : if (fTuneMConData && ((fTuneMConDataMask & kDetITS) == kDetITS)) dedx = GetITSsignalTunedOnData(track);
2563 :
2564 : Double_t momITS=mom;
2565 0 : UChar_t clumap=track->GetITSClusterMap();
2566 : Int_t nPointsForPid=0;
2567 0 : for(Int_t i=2; i<6; i++){
2568 0 : if(clumap&(1<<i)) ++nPointsForPid;
2569 : }
2570 :
2571 : Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
2572 0 : for (Int_t j=0; j<nSpecies; j++) {
2573 0 : const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
2574 : //TODO: in case of the electron, use the SA parametrisation,
2575 : // this needs to be changed if ITS provides a parametrisation
2576 : // for electrons also for ITS+TPC tracks
2577 0 : Double_t bethe=fITSResponse.Bethe(momITS,(AliPID::EParticleType)j,isSA || (j==(Int_t)AliPID::kElectron))*chargeFactor;
2578 0 : Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
2579 0 : Double_t nSigma=fITSResponse.GetNumberOfSigmas(track, (AliPID::EParticleType)j);
2580 0 : if (TMath::Abs(nSigma) > fRange) {
2581 0 : p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
2582 0 : } else {
2583 0 : p[j]=TMath::Exp(-0.5*nSigma*nSigma)/sigma;
2584 : mismatch=kFALSE;
2585 : }
2586 : }
2587 :
2588 0 : if (mismatch){
2589 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2590 0 : }
2591 :
2592 : return kDetPidOk;
2593 0 : }
2594 : //______________________________________________________________________________
2595 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2596 : {
2597 : //
2598 : // Compute PID response for the TPC
2599 : //
2600 :
2601 : // set flat distribution (no decision)
2602 6048 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2603 :
2604 288 : const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
2605 442 : if (pidStatus==kDetNoSignal) return pidStatus;
2606 :
2607 134 : Double_t dedx=track->GetTPCsignal();
2608 : Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
2609 :
2610 134 : if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) dedx = GetTPCsignalTunedOnData(track);
2611 :
2612 : Double_t bethe = 0.;
2613 : Double_t sigma = 0.;
2614 :
2615 2680 : for (Int_t j=0; j<nSpecies; j++) {
2616 : AliPID::EParticleType type=AliPID::EParticleType(j);
2617 :
2618 1206 : bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
2619 1206 : sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
2620 :
2621 1206 : if (TMath::Abs(dedx-bethe) > fRange*sigma) {
2622 710 : p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
2623 710 : } else {
2624 496 : p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
2625 : mismatch=kFALSE;
2626 : }
2627 : }
2628 :
2629 134 : if (mismatch){
2630 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2631 0 : }
2632 :
2633 : return pidStatus;
2634 288 : }
2635 : //______________________________________________________________________________
2636 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],Bool_t kNoMism) const
2637 : {
2638 : //
2639 : // Compute PID probabilities for TOF
2640 : //
2641 :
2642 0 : fgTOFmismatchProb = 1E-8;
2643 :
2644 : // centrality --> fCurrCentrality
2645 : // Beam type --> fBeamTypeNum
2646 : // N TOF cluster --> TOF header --> to get the TOF header we need to add a virtual method in AliVTrack extended to ESD and AOD tracks
2647 : // isMC --> fIsMC
2648 0 : Float_t pt = track->Pt();
2649 0 : Float_t mismPropagationFactor[10] = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
2650 0 : if(! (kNoMism | fNoTOFmism)){ // this flag allows to disable mismatch for iterative procedure to get priors
2651 0 : mismPropagationFactor[3] = 1 + TMath::Exp(1 - 1.12*pt);// it has to be alligned with the one in AliPIDCombined
2652 0 : mismPropagationFactor[4] = 1 + 1./(4.71114 - 5.72372*pt + 2.94715*pt*pt);// it has to be alligned with the one in AliPIDCombined
2653 :
2654 : Int_t nTOFcluster = 0;
2655 0 : if(track->GetTOFHeader() && track->GetTOFHeader()->GetTriggerMask() && track->GetTOFHeader()->GetNumberOfTOFclusters() > -1){ // N TOF clusters available
2656 0 : nTOFcluster = track->GetTOFHeader()->GetNumberOfTOFclusters();
2657 0 : if(fIsMC) nTOFcluster = Int_t(nTOFcluster * 1.5); // +50% in MC
2658 : }
2659 : else{
2660 0 : switch(fBeamTypeNum){
2661 : case kPP: // pp
2662 : nTOFcluster = 80;
2663 0 : break;
2664 : case kPPB: // pPb 5.05 ATeV
2665 0 : nTOFcluster = Int_t(308 - 2.12*fCurrCentrality + TMath::Exp(4.917 -0.1604*fCurrCentrality));
2666 0 : break;
2667 : case kPBPB: // PbPb 2.76 ATeV
2668 0 : nTOFcluster = Int_t(TMath::Exp(9.4 - 0.022*fCurrCentrality));
2669 0 : break;
2670 : }
2671 : }
2672 :
2673 0 : switch(fBeamTypeNum){ // matching window factors for 3 cm and 10 cm (about (10/3)^2)
2674 : case kPP: // pp 7 TeV
2675 0 : nTOFcluster *= 10;
2676 0 : break;
2677 : case kPPB: // pPb 5.05 ATeV
2678 0 : nTOFcluster *= 10;
2679 0 : break;
2680 : case kPBPB: // pPb 5.05 ATeV
2681 : // nTOFcluster *= 1;
2682 : break;
2683 : }
2684 :
2685 0 : if(nTOFcluster < 0) nTOFcluster = 10;
2686 :
2687 :
2688 0 : fgTOFmismatchProb=fTOFResponse.GetMismatchProbability(track->GetTOFsignal(),track->Eta()) * nTOFcluster *6E-6 * (1 + 2.90505e-01/pt/pt); // mism weight * tof occupancy (including matching window factor) * pt dependence
2689 :
2690 0 : }
2691 :
2692 : // set flat distribution (no decision)
2693 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2694 :
2695 0 : const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
2696 0 : if (pidStatus!=kDetPidOk) return pidStatus;
2697 :
2698 0 : const Double_t meanCorrFactor = 0.07/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2699 :
2700 0 : for (Int_t j=0; j<nSpecies; j++) {
2701 : AliPID::EParticleType type=AliPID::EParticleType(j);
2702 0 : const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2703 :
2704 0 : const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
2705 0 : const Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
2706 :
2707 0 : if(nsigmas < fTOFtail)
2708 0 : p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
2709 : else
2710 0 : p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
2711 :
2712 0 : p[j] += fgTOFmismatchProb*mismPropagationFactor[j];
2713 : }
2714 :
2715 : return kDetPidOk;
2716 0 : }
2717 :
2718 : Int_t AliPIDResponse::CalculateTRDResponse(const AliVTrack *track,Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
2719 : {
2720 : // new function for backward compatibility
2721 : // returns number of tracklets PID
2722 :
2723 0 : UInt_t TRDslicesForPID[2];
2724 0 : SetTRDSlices(TRDslicesForPID,PIDmethod);
2725 :
2726 0 : Float_t mom[6]={0.};
2727 0 : Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
2728 0 : Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2729 0 : for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2730 0 : mom[ilayer] = track->GetTRDmomentum(ilayer);
2731 0 : for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2732 : // do not consider tracklets with no momentum (should be non functioning chambers)
2733 0 : if(mom[ilayer]>0) dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2734 : }
2735 : }
2736 :
2737 0 : return fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2738 :
2739 0 : }
2740 : //______________________________________________________________________________
2741 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
2742 : {
2743 : //
2744 : // Compute PID probabilities for the TRD
2745 : //
2746 :
2747 : // set flat distribution (no decision)
2748 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2749 0 : const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
2750 0 : if (pidStatus!=kDetPidOk) return pidStatus;
2751 :
2752 0 : CalculateTRDResponse(track,p,PIDmethod);
2753 :
2754 0 : return kDetPidOk;
2755 0 : }
2756 :
2757 : //______________________________________________________________________________
2758 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2759 : {
2760 : //
2761 : // Compute PID response for the EMCAL
2762 : //
2763 :
2764 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2765 :
2766 0 : const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
2767 0 : if (pidStatus!=kDetPidOk) return pidStatus;
2768 :
2769 0 : const Int_t nMatchClus = track->GetEMCALcluster();
2770 0 : AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2771 :
2772 0 : const Double_t mom = track->P();
2773 0 : const Double_t pt = track->Pt();
2774 0 : const Int_t charge = track->Charge();
2775 0 : const Double_t fClsE = matchedClus->E();
2776 0 : const Double_t EovP = fClsE/mom;
2777 :
2778 : // compute the probabilities
2779 0 : fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p);
2780 : return kDetPidOk;
2781 0 : }
2782 :
2783 : //______________________________________________________________________________
2784 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2785 : {
2786 : //
2787 : // Compute PID response for the PHOS
2788 : //
2789 :
2790 : // set flat distribution (no decision)
2791 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2792 0 : return kDetNoSignal;
2793 : }
2794 :
2795 : //______________________________________________________________________________
2796 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2797 : {
2798 : //
2799 : // Compute PID response for the HMPID
2800 : //
2801 :
2802 : // set flat distribution (no decision)
2803 0 : for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2804 :
2805 0 : const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2806 0 : if (pidStatus!=kDetPidOk) return pidStatus;
2807 :
2808 0 : fHMPIDResponse.GetProbability(track,nSpecies,p);
2809 :
2810 0 : return kDetPidOk;
2811 0 : }
2812 :
2813 : //______________________________________________________________________________
2814 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetITSPIDStatus(const AliVTrack *track) const
2815 : {
2816 : // compute ITS pid status
2817 :
2818 : // check status bits
2819 0 : if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
2820 0 : (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
2821 :
2822 0 : const Float_t dEdx=track->GetITSsignal();
2823 0 : if (dEdx<=0) return kDetNoSignal;
2824 :
2825 : // requite at least 3 pid clusters
2826 0 : const UChar_t clumap=track->GetITSClusterMap();
2827 : Int_t nPointsForPid=0;
2828 0 : for(Int_t i=2; i<6; i++){
2829 0 : if(clumap&(1<<i)) ++nPointsForPid;
2830 : }
2831 :
2832 0 : if(nPointsForPid<3) {
2833 0 : return kDetNoSignal;
2834 : }
2835 :
2836 0 : return kDetPidOk;
2837 0 : }
2838 :
2839 : //______________________________________________________________________________
2840 : AliPIDResponse::EDetPidStatus AliPIDResponse:: GetTPCPIDStatus(const AliVTrack *track) const
2841 : {
2842 : // compute TPC pid status
2843 :
2844 : // check quality of the track
2845 608 : if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
2846 :
2847 : // check pid values
2848 272 : const Double_t dedx=track->GetTPCsignal();
2849 272 : const UShort_t signalN=track->GetTPCsignalN();
2850 410 : if (signalN<10 || dedx<10) return kDetNoSignal;
2851 :
2852 268 : if (!fTPCResponse.GetResponseFunction(AliPID::kPion)) return kDetNoParams;
2853 :
2854 0 : return kDetPidOk;
2855 288 : }
2856 :
2857 : //______________________________________________________________________________
2858 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetTRDPIDStatus(const AliVTrack *track) const
2859 : {
2860 : // compute TRD pid status
2861 :
2862 0 : if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2863 0 : return kDetPidOk;
2864 0 : }
2865 :
2866 : //______________________________________________________________________________
2867 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetTOFPIDStatus(const AliVTrack *track) const
2868 : {
2869 : // compute TOF pid status
2870 :
2871 0 : if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2872 0 : if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2873 :
2874 0 : return kDetPidOk;
2875 0 : }
2876 :
2877 : //______________________________________________________________________________
2878 : Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
2879 : {
2880 : // compute mismatch probability cross-checking at 5 sigmas with TPC
2881 : // currently just implemented as a 5 sigma compatibility cut
2882 :
2883 0 : if(!track) return fgTOFmismatchProb;
2884 :
2885 : // check pid status
2886 0 : const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
2887 0 : if (tofStatus!=kDetPidOk) return 0.;
2888 :
2889 : //mismatch
2890 0 : const EDetPidStatus tpcStatus=GetTPCPIDStatus(track);
2891 0 : if (tpcStatus==kDetNoSignal) return 0.;
2892 :
2893 0 : const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2894 : Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2895 0 : for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
2896 : AliPID::EParticleType type=AliPID::EParticleType(j);
2897 0 : const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2898 :
2899 0 : if (TMath::Abs(nsigmas)<5.){
2900 0 : const Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2901 0 : if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2902 0 : }
2903 : }
2904 :
2905 0 : if (mismatch){
2906 0 : return 1.;
2907 : }
2908 :
2909 0 : return 0.;
2910 0 : }
2911 :
2912 : //______________________________________________________________________________
2913 : AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
2914 : {
2915 : // compute HMPID pid status
2916 :
2917 0 : Int_t ch = track->GetHMPIDcluIdx()/1000000;
2918 0 : Double_t HMPIDsignal = track->GetHMPIDsignal();
2919 :
2920 0 : if((track->GetStatus()&AliVTrack::kHMPIDpid)==0 || ch<0 || ch>6 || HMPIDsignal<0) return kDetNoSignal;
2921 :
2922 0 : return kDetPidOk;
2923 0 : }
2924 :
2925 : //______________________________________________________________________________
2926 : AliPIDResponse::EDetPidStatus AliPIDResponse:: GetPHOSPIDStatus(const AliVTrack */*track*/) const
2927 : {
2928 : // compute PHOS pid status
2929 0 : return kDetNoSignal;
2930 : }
2931 :
2932 : //______________________________________________________________________________
2933 : AliPIDResponse::EDetPidStatus AliPIDResponse:: GetEMCALPIDStatus(const AliVTrack *track) const
2934 : {
2935 : // compute EMCAL pid status
2936 :
2937 :
2938 : // Track matching
2939 0 : const Int_t nMatchClus = track->GetEMCALcluster();
2940 0 : if (nMatchClus<0) return kDetNoSignal;
2941 :
2942 0 : AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2943 :
2944 0 : if (!(matchedClus && matchedClus->IsEMCAL())) return kDetNoSignal;
2945 :
2946 0 : const Int_t charge = track->Charge();
2947 0 : if (TMath::Abs(charge)!=1) return kDetNoSignal;
2948 :
2949 0 : if (!(fEMCALPIDParams && fEMCALPIDParams->At(AliPID::kElectron))) return kDetNoParams;
2950 :
2951 0 : return kDetPidOk;
2952 :
2953 0 : }
2954 :
2955 : //______________________________________________________________________________
2956 : AliPIDResponse::EDetPidStatus AliPIDResponse::GetPIDStatus(EDetector detector, const AliVTrack *track) const
2957 : {
2958 : //
2959 : // check pid status for a track
2960 : //
2961 :
2962 0 : switch (detector){
2963 0 : case kITS: return GetITSPIDStatus(track); break;
2964 0 : case kTPC: return GetTPCPIDStatus(track); break;
2965 0 : case kTRD: return GetTRDPIDStatus(track); break;
2966 0 : case kTOF: return GetTOFPIDStatus(track); break;
2967 0 : case kPHOS: return GetPHOSPIDStatus(track); break;
2968 0 : case kEMCAL: return GetEMCALPIDStatus(track); break;
2969 0 : case kHMPID: return GetHMPIDPIDStatus(track); break;
2970 0 : default: return kDetNoSignal;
2971 : }
2972 : return kDetNoSignal;
2973 :
2974 0 : }
2975 :
2976 : //_________________________________________________________________________
2977 : Float_t AliPIDResponse::GetITSsignalTunedOnData(const AliVTrack *t) const
2978 : {
2979 : /// Create gaussian signal response based on the dE/dx response observed in data
2980 : /// Currently only for deuterons and triton. The other particles are fine in MC
2981 :
2982 0 : Float_t dedx = t->GetITSsignalTunedOnData();
2983 0 : if(dedx > 0) return dedx;
2984 :
2985 0 : dedx = t->GetITSsignal();
2986 0 : ((AliVTrack*)t)->SetITSsignalTunedOnData(dedx);
2987 0 : if(dedx < 20) return dedx;
2988 :
2989 : Int_t nITSpid = 0;
2990 0 : for(Int_t i=2; i<6; i++){
2991 0 : if(t->HasPointOnITSLayer(i)) nITSpid++;
2992 : }
2993 :
2994 : Bool_t isSA=kTRUE;
2995 0 : if( t->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
2996 :
2997 : // check if we have MC information
2998 0 : if (!fCurrentMCEvent) {
2999 0 : AliError("Tune On Data requested, but MC event not set. Call 'SetCurrentMCEvent' before!");
3000 0 : return dedx;
3001 : }
3002 :
3003 : // get MC particle
3004 0 : AliVParticle *mcPart = fCurrentMCEvent->GetTrack(TMath::Abs(t->GetLabel()));
3005 :
3006 0 : if (mcPart != NULL) { // protect against label-0 track (initial proton in Pythia events)
3007 : AliPID::EParticleType type = AliPID::kPion;
3008 : Bool_t kGood = kFALSE;
3009 0 : Int_t iS = TMath::Abs(mcPart->PdgCode());
3010 :
3011 0 : for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart) {
3012 0 : if (iS == AliPID::ParticleCode(ipart)) {
3013 : type = static_cast<AliPID::EParticleType>(ipart);
3014 : kGood=kTRUE;
3015 0 : break;
3016 : }
3017 : }
3018 :
3019 : // NOTE: currently only tune for deuterons and triton
3020 0 : if(kGood && (iS == AliPID::ParticleCode(AliPID::kDeuteron) || (iS == AliPID::ParticleCode(AliPID::kTriton))) ) {
3021 0 : Double_t bethe = fITSResponse.Bethe(t->P(),type,isSA);
3022 0 : Double_t sigma = fITSResponse.GetResolution(bethe,nITSpid,isSA,t->P(),type);
3023 0 : dedx = gRandom->Gaus(bethe,sigma);
3024 0 : }
3025 0 : }
3026 :
3027 0 : const_cast<AliVTrack*>(t)->SetITSsignalTunedOnData(dedx);
3028 : return dedx;
3029 0 : }
3030 :
3031 : //_________________________________________________________________________
3032 : Float_t AliPIDResponse::GetTPCsignalTunedOnData(const AliVTrack *t) const
3033 : {
3034 : /// Create gaussian signal response based on the dE/dx response observed in data
3035 :
3036 0 : Float_t dedx = t->GetTPCsignalTunedOnData();
3037 :
3038 0 : if(dedx > 0) return dedx;
3039 :
3040 0 : dedx = t->GetTPCsignal();
3041 0 : ((AliVTrack*)t)->SetTPCsignalTunedOnData(dedx);
3042 0 : if(dedx < 20) return dedx;
3043 :
3044 : // check if we have MC information
3045 0 : if (!fCurrentMCEvent) {
3046 0 : AliError("Tune On Data requested, but MC event not set. Call 'SetCurrentMCEvent' before!");
3047 0 : return dedx;
3048 : }
3049 :
3050 : // get MC particle
3051 0 : AliVParticle *mcPart = fCurrentMCEvent->GetTrack(TMath::Abs(t->GetLabel()));
3052 :
3053 0 : if (mcPart != NULL) { // protect against label-0 track (initial proton in Pythia events)
3054 : AliPID::EParticleType type = AliPID::kPion;
3055 : Bool_t kGood = kFALSE;
3056 0 : Int_t iS = TMath::Abs(mcPart->PdgCode());
3057 :
3058 0 : for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart) {
3059 0 : if (iS == AliPID::ParticleCode(ipart)) {
3060 : type = static_cast<AliPID::EParticleType>(ipart);
3061 : kGood=kTRUE;
3062 0 : break;
3063 : }
3064 : }
3065 :
3066 0 : if(kGood){
3067 : //TODO maybe introduce different dEdxSources?
3068 0 : Double_t bethe = fTPCResponse.GetExpectedSignal(t, type, AliTPCPIDResponse::kdEdxDefault, UseTPCEtaCorrection(),
3069 0 : UseTPCMultiplicityCorrection());
3070 0 : Double_t sigma = fTPCResponse.GetExpectedSigma(t, type, AliTPCPIDResponse::kdEdxDefault, UseTPCEtaCorrection(),
3071 0 : UseTPCMultiplicityCorrection());
3072 0 : dedx = gRandom->Gaus(bethe,sigma);
3073 : // if(iS == AliPID::ParticleCode(AliPID::kHe3) || iS == AliPID::ParticleCode(AliPID::kAlpha)) dedx *= 5;
3074 0 : }
3075 0 : }
3076 :
3077 0 : const_cast<AliVTrack*>(t)->SetTPCsignalTunedOnData(dedx);
3078 : return dedx;
3079 0 : }
3080 :
3081 : //_________________________________________________________________________
3082 : Float_t AliPIDResponse::GetTOFsignalTunedOnData(const AliVTrack *t) const
3083 : {
3084 : /// Calculate the TOF signal tuned on data by adding a tail
3085 0 : Double_t tofSignal = t->GetTOFsignalTunedOnData();
3086 :
3087 0 : if(tofSignal < 99999) return (Float_t)tofSignal; // it has been already set
3088 :
3089 : // read additional mismatch fraction
3090 0 : Float_t addmism = GetTOFPIDParams()->GetTOFadditionalMismForMC();
3091 0 : if(addmism > 1.){
3092 0 : Float_t centr = GetCurrentCentrality();
3093 0 : if(centr > 50) addmism *= 0.1667;
3094 0 : else if(centr > 20) addmism *= 0.33;
3095 0 : }
3096 :
3097 0 : tofSignal = t->GetTOFsignal() + fTOFResponse.GetTailRandomValue(t->Pt(),t->Eta(),t->GetTOFsignal(),addmism);
3098 0 : const_cast<AliVTrack*>(t)->SetTOFsignalTunedOnData(tofSignal);
3099 0 : return (Float_t)tofSignal;
3100 0 : }
|
