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