Updates for pA and new flags:
[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),
60fOADBPath(),
00a38d07 61fCustomTPCpidResponse(),
4ec8e76d 62fBeamType("PP"),
63fLHCperiod(),
64fMCperiodTPC(),
fd21ec8d 65fMCperiodUser(),
ea235c90 66fCurrentFile(),
4ec8e76d 67fRecoPass(0),
fd21ec8d 68fRecoPassUser(-1),
4ec8e76d 69fRun(0),
70fOldRun(0),
71fArrPidResponseMaster(0x0),
72fResolutionCorrection(0x0),
db0e2c5f 73fTRDPIDResponseObject(0x0),
0b39f221 74fTOFtail(1.1),
b79db598 75fTOFPIDParams(0x0),
b2138b40 76fEMCALPIDParams(0x0),
80f28562 77fCurrentEvent(0x0),
539a5a59 78fCurrCentrality(0.0),
79fTuneMConData(kFALSE)
4ec8e76d 80{
81 //
82 // default ctor
83 //
a635821f 84 AliLog::SetClassDebugLevel("AliPIDResponse",0);
85 AliLog::SetClassDebugLevel("AliESDpid",0);
86 AliLog::SetClassDebugLevel("AliAODpidUtil",0);
ea235c90 87
88 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 89}
90
91//______________________________________________________________________________
92AliPIDResponse::~AliPIDResponse()
93{
94 //
95 // dtor
96 //
00a38d07 97 delete fArrPidResponseMaster;
98 delete fTRDPIDResponseObject;
99 delete fTOFPIDParams;
4ec8e76d 100}
101
102//______________________________________________________________________________
103AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
104TNamed(other),
105fITSResponse(other.fITSResponse),
106fTPCResponse(other.fTPCResponse),
107fTRDResponse(other.fTRDResponse),
108fTOFResponse(other.fTOFResponse),
e96b9916 109fEMCALResponse(other.fEMCALResponse),
fd21ec8d 110fRange(other.fRange),
111fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 112fIsMC(other.fIsMC),
113fOADBPath(other.fOADBPath),
00a38d07 114fCustomTPCpidResponse(other.fCustomTPCpidResponse),
4ec8e76d 115fBeamType("PP"),
116fLHCperiod(),
117fMCperiodTPC(),
fd21ec8d 118fMCperiodUser(other.fMCperiodUser),
ea235c90 119fCurrentFile(),
4ec8e76d 120fRecoPass(0),
fd21ec8d 121fRecoPassUser(other.fRecoPassUser),
4ec8e76d 122fRun(0),
123fOldRun(0),
124fArrPidResponseMaster(0x0),
125fResolutionCorrection(0x0),
db0e2c5f 126fTRDPIDResponseObject(0x0),
0b39f221 127fTOFtail(1.1),
b79db598 128fTOFPIDParams(0x0),
b2138b40 129fEMCALPIDParams(0x0),
80f28562 130fCurrentEvent(0x0),
539a5a59 131fCurrCentrality(0.0),
132fTuneMConData(kFALSE)
4ec8e76d 133{
134 //
135 // copy ctor
136 //
ea235c90 137 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 138}
139
140//______________________________________________________________________________
141AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
142{
143 //
144 // copy ctor
145 //
146 if(this!=&other) {
147 delete fArrPidResponseMaster;
148 TNamed::operator=(other);
149 fITSResponse=other.fITSResponse;
150 fTPCResponse=other.fTPCResponse;
151 fTRDResponse=other.fTRDResponse;
152 fTOFResponse=other.fTOFResponse;
e96b9916 153 fEMCALResponse=other.fEMCALResponse;
fd21ec8d 154 fRange=other.fRange;
155 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 156 fOADBPath=other.fOADBPath;
00a38d07 157 fCustomTPCpidResponse=other.fCustomTPCpidResponse;
4ec8e76d 158 fIsMC=other.fIsMC;
159 fBeamType="PP";
160 fLHCperiod="";
161 fMCperiodTPC="";
fd21ec8d 162 fMCperiodUser=other.fMCperiodUser;
ea235c90 163 fCurrentFile="";
4ec8e76d 164 fRecoPass=0;
fd21ec8d 165 fRecoPassUser=other.fRecoPassUser;
4ec8e76d 166 fRun=0;
167 fOldRun=0;
168 fArrPidResponseMaster=0x0;
169 fResolutionCorrection=0x0;
db0e2c5f 170 fTRDPIDResponseObject=0x0;
b2138b40 171 fEMCALPIDParams=0x0;
ea235c90 172 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
0b39f221 173 fTOFtail=1.1;
b79db598 174 fTOFPIDParams=0x0;
e96b9916 175 fCurrentEvent=other.fCurrentEvent;
4ec8e76d 176 }
177 return *this;
178}
179
180//______________________________________________________________________________
fd21ec8d 181Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
182{
183 //
184 // NumberOfSigmas for 'detCode'
185 //
186
187 switch (detCode){
188 case kDetITS: return NumberOfSigmasITS(track, type); break;
189 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
190 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
00a38d07 191 case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
fd21ec8d 192 default: return -999.;
193 }
194
195}
196
197//______________________________________________________________________________
00a38d07 198Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
199{
200 //
201 // NumberOfSigmas for 'detCode'
202 //
203 return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
204}
205
206//______________________________________________________________________________
207Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
208{
209 //
210 // Calculate the number of sigmas in the ITS
211 //
212
213 AliVTrack *track=(AliVTrack*)vtrack;
214
215 // look for cached value first
216 // only the non SA tracks are cached
217 if ( track->GetDetectorPID() ){
218 return track->GetDetectorPID()->GetNumberOfSigmas(kITS, type);
219 }
220
221 Float_t dEdx=track->GetITSsignal();
222 if (dEdx<=0) return -999.;
223
224 UChar_t clumap=track->GetITSClusterMap();
225 Int_t nPointsForPid=0;
226 for(Int_t i=2; i<6; i++){
227 if(clumap&(1<<i)) ++nPointsForPid;
228 }
229 Float_t mom=track->P();
230
231 //check for ITS standalone tracks
232 Bool_t isSA=kTRUE;
233 if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
234
235 //TODO: in case of the electron, use the SA parametrisation,
236 // this needs to be changed if ITS provides a parametrisation
237 // for electrons also for ITS+TPC tracks
238 return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
239}
240
241//______________________________________________________________________________
242Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
243{
244 //
245 // Calculate the number of sigmas in the TPC
246 //
247
248 AliVTrack *track=(AliVTrack*)vtrack;
249
250 // look for cached value first
251 if (track->GetDetectorPID()){
252 return track->GetDetectorPID()->GetNumberOfSigmas(kTPC, type);
253 }
254
255 Double_t mom = track->GetTPCmomentum();
256 Double_t sig = track->GetTPCsignal();
257 UInt_t sigN = track->GetTPCsignalN();
258
259 Double_t nSigma = -999.;
260 if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
261
262 return nSigma;
263}
264
265//______________________________________________________________________________
266Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
267{
268 //
269 // Calculate the number of sigmas in the EMCAL
270 //
271
272 AliVTrack *track=(AliVTrack*)vtrack;
e96b9916 273
00a38d07 274 // look for cached value first
275 if (track->GetDetectorPID()){
276 return track->GetDetectorPID()->GetNumberOfSigmas(kEMCAL, type);
277 }
278
e96b9916 279 AliVCluster *matchedClus = NULL;
280
281 Double_t mom = -1.;
282 Double_t pt = -1.;
283 Double_t EovP = -1.;
284 Double_t fClsE = -1.;
285
286 Int_t nMatchClus = -1;
287 Int_t charge = 0;
288
289 // Track matching
290 nMatchClus = track->GetEMCALcluster();
291 if(nMatchClus > -1){
6d0064aa 292
e96b9916 293 mom = track->P();
294 pt = track->Pt();
295 charge = track->Charge();
296
297 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
298
299 if(matchedClus){
300
6d0064aa 301 // matched cluster is EMCAL
302 if(matchedClus->IsEMCAL()){
303
304 fClsE = matchedClus->E();
305 EovP = fClsE/mom;
306
307
308 // NSigma value really meaningful only for electrons!
309 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
310 }
e96b9916 311 }
312 }
6d0064aa 313
e96b9916 314 return -999;
315
316}
317
318//______________________________________________________________________________
00a38d07 319Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
6d0064aa 320
00a38d07 321 AliVTrack *track=(AliVTrack*)vtrack;
322
6d0064aa 323 AliVCluster *matchedClus = NULL;
324
325 Double_t mom = -1.;
326 Double_t pt = -1.;
327 Double_t EovP = -1.;
328 Double_t fClsE = -1.;
32fa24d6 329
330 // initialize eop and shower shape parameters
331 eop = -1.;
332 for(Int_t i = 0; i < 4; i++){
333 showershape[i] = -1.;
334 }
6d0064aa 335
336 Int_t nMatchClus = -1;
337 Int_t charge = 0;
338
339 // Track matching
340 nMatchClus = track->GetEMCALcluster();
341 if(nMatchClus > -1){
342
343 mom = track->P();
344 pt = track->Pt();
345 charge = track->Charge();
346
347 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
348
349 if(matchedClus){
350
351 // matched cluster is EMCAL
352 if(matchedClus->IsEMCAL()){
353
354 fClsE = matchedClus->E();
355 EovP = fClsE/mom;
356
357 // fill used EMCAL variables here
358 eop = EovP; // E/p
359 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
360 showershape[1] = matchedClus->GetM02(); // long axis
361 showershape[2] = matchedClus->GetM20(); // short axis
362 showershape[3] = matchedClus->GetDispersion(); // dispersion
363
364 // NSigma value really meaningful only for electrons!
365 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
366 }
367 }
368 }
369 return -999;
370}
371
372
373//______________________________________________________________________________
fd21ec8d 374AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
375{
376 //
377 // Compute PID response of 'detCode'
378 //
379
380 switch (detCode){
381 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
382 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
fd21ec8d 383 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
00a38d07 384 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
fd21ec8d 385 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
386 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
387 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
388 default: return kDetNoSignal;
389 }
390}
391
392//______________________________________________________________________________
00a38d07 393AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
394{
395 //
396 // Compute PID response of 'detCode'
397 //
398
399 return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
400}
401
402//______________________________________________________________________________
fd21ec8d 403AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
404{
405 //
406 // Compute PID response for the ITS
407 //
408
00a38d07 409 // look for cached value first
410 // only the non SA tracks are cached
411 if (track->GetDetectorPID()){
412 return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
413 }
414
fd21ec8d 415 // set flat distribution (no decision)
416 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
417
418 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
419 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
00a38d07 420
421 //check for ITS standalone tracks
422 Bool_t isSA=kTRUE;
423 if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
fd21ec8d 424
425 Double_t mom=track->P();
426 Double_t dedx=track->GetITSsignal();
fd21ec8d 427 Double_t momITS=mom;
fd21ec8d 428 UChar_t clumap=track->GetITSClusterMap();
429 Int_t nPointsForPid=0;
430 for(Int_t i=2; i<6; i++){
431 if(clumap&(1<<i)) ++nPointsForPid;
432 }
433
434 if(nPointsForPid<3) { // track not to be used for combined PID purposes
435 // track->ResetStatus(AliVTrack::kITSpid);
436 return kDetNoSignal;
437 }
438
439 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
440 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
00a38d07 441 Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
442 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
443 Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
444 //TODO: in case of the electron, use the SA parametrisation,
445 // this needs to be changed if ITS provides a parametrisation
446 // for electrons also for ITS+TPC tracks
447 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
fd21ec8d 448 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
449 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
450 } else {
451 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
452 mismatch=kFALSE;
453 }
454
455 // Check for particles heavier than (AliPID::kSPECIES - 1)
456 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
457
458 }
459
460 if (mismatch){
461 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
462 return kDetNoSignal;
463 }
464
465
466 return kDetPidOk;
467}
468//______________________________________________________________________________
469AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
470{
471 //
472 // Compute PID response for the TPC
473 //
00a38d07 474
475 // look for cached value first
476 if (track->GetDetectorPID()){
477 return track->GetDetectorPID()->GetRawProbability(kTPC, p, nSpecies);
478 }
479
fd21ec8d 480 // set flat distribution (no decision)
481 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
482
483 // check quality of the track
484 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
485
486 Double_t mom = track->GetTPCmomentum();
487
488 Double_t dedx=track->GetTPCsignal();
489 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
490
539a5a59 491 if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
492
00a38d07 493 for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
fd21ec8d 494 AliPID::EParticleType type=AliPID::EParticleType(j);
495 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
496 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
497 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
498 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
499 } else {
500 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
501 mismatch=kFALSE;
502 }
fd21ec8d 503 }
504
505 if (mismatch){
506 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
507 return kDetNoSignal;
508 }
509
510 return kDetPidOk;
511}
512//______________________________________________________________________________
513AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
514{
515 //
516 // Compute PID response for the
517 //
00a38d07 518
519 // look for cached value first
520 if (track->GetDetectorPID()){
521 return track->GetDetectorPID()->GetRawProbability(kTOF, p, nSpecies);
522 }
523
0b39f221 524 Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
525
fd21ec8d 526 // set flat distribution (no decision)
527 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
528
529 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
530 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
531
fd21ec8d 532 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
00a38d07 533 for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
fd21ec8d 534 AliPID::EParticleType type=AliPID::EParticleType(j);
0b39f221 535 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
536
67376d1d 537 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
538 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
fd21ec8d 539 if (TMath::Abs(nsigmas) > (fRange+2)) {
0b39f221 540 if(nsigmas < fTOFtail)
541 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
542 else
543 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
544 } else{
545 if(nsigmas < fTOFtail)
546 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
547 else
548 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
549 }
550
fd21ec8d 551 if (TMath::Abs(nsigmas)<5.){
552 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
553 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
554 }
555 }
556
557 if (mismatch){
558 return kDetMismatch;
559 }
560
561 // TODO: Light nuclei
562
563 return kDetPidOk;
564}
565//______________________________________________________________________________
ea235c90 566AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 567{
568 //
569 // Compute PID response for the
570 //
00a38d07 571
572 // look for cached value first
573 if (track->GetDetectorPID()){
574 return track->GetDetectorPID()->GetRawProbability(kTRD, p, nSpecies);
575 }
576
fd21ec8d 577 // set flat distribution (no decision)
578 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 579 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
580
00a38d07 581 Float_t mom[6]={0.};
582 Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
ea235c90 583 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
584 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
585 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
586 mom[ilayer] = track->GetTRDmomentum(ilayer);
587 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
588 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
589 }
590 }
591 fTRDResponse.GetResponse(nslices, dedx, mom, p);
592 return kDetPidOk;
fd21ec8d 593}
594//______________________________________________________________________________
e96b9916 595AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 596{
597 //
598 // Compute PID response for the EMCAL
599 //
00a38d07 600
601 // look for cached value first
602 if (track->GetDetectorPID()){
603 return track->GetDetectorPID()->GetRawProbability(kEMCAL, p, nSpecies);
604 }
605
606 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
fd21ec8d 607
e96b9916 608 AliVCluster *matchedClus = NULL;
609
610 Double_t mom = -1.;
611 Double_t pt = -1.;
612 Double_t EovP = -1.;
613 Double_t fClsE = -1.;
614
615 Int_t nMatchClus = -1;
616 Int_t charge = 0;
617
618 // Track matching
619 nMatchClus = track->GetEMCALcluster();
620
621 if(nMatchClus > -1){
00a38d07 622
e96b9916 623 mom = track->P();
624 pt = track->Pt();
625 charge = track->Charge();
626
627 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
628
629 if(matchedClus){
00a38d07 630
631 // matched cluster is EMCAL
632 if(matchedClus->IsEMCAL()){
633
e96b9916 634 fClsE = matchedClus->E();
635 EovP = fClsE/mom;
636
637
638 // compute the probabilities
00a38d07 639 if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){
640
641 // in case everything is OK
642 return kDetPidOk;
643 }
e96b9916 644 }
645 }
646 }
e96b9916 647
648 // in all other cases set flat distribution (no decision)
00a38d07 649 for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
fd21ec8d 650 return kDetNoSignal;
e96b9916 651
fd21ec8d 652}
653//______________________________________________________________________________
654AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
655{
656 //
657 // Compute PID response for the PHOS
658 //
00a38d07 659
660 // look for cached value first
661// if (track->GetDetectorPID()){
662// return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
663// }
664
fd21ec8d 665 // set flat distribution (no decision)
666 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
667 return kDetNoSignal;
668}
669//______________________________________________________________________________
ea235c90 670AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 671{
672 //
673 // Compute PID response for the HMPID
674 //
675
00a38d07 676
677 // look for cached value first
678 if (track->GetDetectorPID()){
679 return track->GetDetectorPID()->GetRawProbability(kHMPID, p, nSpecies);
680 }
681
fd21ec8d 682 // set flat distribution (no decision)
683 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 684 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
685
686 track->GetHMPIDpid(p);
687
688 return kDetPidOk;
fd21ec8d 689}
690
691//______________________________________________________________________________
00a38d07 692void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
4ec8e76d 693{
694 //
695 // Apply settings for the current event
696 //
697 fRecoPass=pass;
e96b9916 698
699 fCurrentEvent=0x0;
4ec8e76d 700 if (!event) return;
e96b9916 701 fCurrentEvent=event;
00a38d07 702 if (run>0) fRun=run;
703 else fRun=event->GetRunNumber();
4ec8e76d 704
705 if (fRun!=fOldRun){
706 ExecNewRun();
707 fOldRun=fRun;
708 }
709
710 //TPC resolution parametrisation PbPb
711 if ( fResolutionCorrection ){
712 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
713 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
714 }
715
716 //TOF resolution
b79db598 717 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
718
80f28562 719
720 // Get and set centrality
721 AliCentrality *centrality = event->GetCentrality();
722 if(centrality){
723 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
724 }
725 else{
726 fCurrCentrality = -1;
727 }
4ec8e76d 728}
729
730//______________________________________________________________________________
731void AliPIDResponse::ExecNewRun()
732{
733 //
734 // Things to Execute upon a new run
735 //
736 SetRecoInfo();
737
738 SetITSParametrisation();
739
740 SetTPCPidResponseMaster();
741 SetTPCParametrisation();
53d016dc 742
743 SetTRDPidResponseMaster();
744 InitializeTRDResponse();
b2138b40 745
746 SetEMCALPidResponseMaster();
747 InitializeEMCALResponse();
4ec8e76d 748
b79db598 749 SetTOFPidResponseMaster();
750 InitializeTOFResponse();
4ec8e76d 751}
752
753//_____________________________________________________
754Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
755{
756 //
757 // Get TPC multiplicity in bins of 150
758 //
759
760 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
761 Double_t tpcMulti=0.;
762 if(vertexTPC){
763 Double_t vertexContribTPC=vertexTPC->GetNContributors();
764 tpcMulti=vertexContribTPC/150.;
765 if (tpcMulti>20.) tpcMulti=20.;
766 }
767
768 return tpcMulti;
769}
770
771//______________________________________________________________________________
772void AliPIDResponse::SetRecoInfo()
773{
774 //
775 // Set reconstruction information
776 //
777
778 //reset information
779 fLHCperiod="";
780 fMCperiodTPC="";
781
782 fBeamType="";
783
784 fBeamType="PP";
785
1436d6bb 786
787 TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
a78fd045 788 TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
1436d6bb 789
4ec8e76d 790 //find the period by run number (UGLY, but not stored in ESD and AOD... )
791 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
792 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
793 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
99e9d5ec 794 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 795 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
796 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
ea235c90 797 else if (fRun>=136851&&fRun<=139517) {
798 fLHCperiod="LHC10H";
799 fMCperiodTPC="LHC10H8";
800 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
801 fBeamType="PBPB";
802 }
3077a03d 803 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
2e97a211 804 //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
805 else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
806 else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
00a38d07 807 // also for 11e,f use 11d
808 else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
809 else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
810
3077a03d 811 else if (fRun>=166529) {
812 fLHCperiod="LHC11H";
813 fMCperiodTPC="LHC11A10";
814 fBeamType="PBPB";
a78fd045 815 if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
3077a03d 816 }
817
80ab5635 818
819 //exception new pp MC productions from 2011
820 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
4a527e08 821 // exception for 11f1
00a38d07 822 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
4ec8e76d 823}
824
825//______________________________________________________________________________
826void AliPIDResponse::SetITSParametrisation()
827{
828 //
829 // Set the ITS parametrisation
830 //
831}
832
833//______________________________________________________________________________
834void AliPIDResponse::SetTPCPidResponseMaster()
835{
836 //
837 // Load the TPC pid response functions from the OADB
838 //
09b50a42 839 //don't load twice for the moment
840 if (fArrPidResponseMaster) return;
841
842
4ec8e76d 843 //reset the PID response functions
844 delete fArrPidResponseMaster;
845 fArrPidResponseMaster=0x0;
846
847 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
00a38d07 848 if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
4ec8e76d 849
ea235c90 850 TFile *f=TFile::Open(fileName.Data());
851 if (f && f->IsOpen() && !f->IsZombie()){
852 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 853 }
ea235c90 854 delete f;
4ec8e76d 855
856 if (!fArrPidResponseMaster){
857 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
858 return;
859 }
860 fArrPidResponseMaster->SetOwner();
861}
862
863//______________________________________________________________________________
864void AliPIDResponse::SetTPCParametrisation()
865{
866 //
867 // Change BB parametrisation for current run
868 //
869
870 if (fLHCperiod.IsNull()) {
871 AliFatal("No period set, not changing parametrisation");
872 return;
873 }
874
875 //
876 // Set default parametrisations for data and MC
877 //
878
879 //data type
880 TString datatype="DATA";
881 //in case of mc fRecoPass is per default 1
882 if (fIsMC) {
539a5a59 883 if(!fTuneMConData) datatype="MC";
884 fRecoPass=1;
4ec8e76d 885 }
886
887 //
888 //reset old splines
889 //
00a38d07 890 for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
4ec8e76d 891 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
892 }
4a527e08 893
894 // period
895 TString period=fLHCperiod;
539a5a59 896 if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
4a527e08 897
898 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
899 Bool_t found=kFALSE;
4ec8e76d 900 //
901 //set the new PID splines
902 //
4ec8e76d 903 if (fArrPidResponseMaster){
539a5a59 904 Int_t recopass = fRecoPass;
905 if(fTuneMConData) recopass = fRecoPassUser;
4ec8e76d 906 TObject *grAll=0x0;
907 //for MC don't use period information
908// if (fIsMC) period="[A-Z0-9]*";
909 //for MC use MC period information
4ec8e76d 910//pattern for the default entry (valid for all particles)
539a5a59 911 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
912
4ec8e76d 913 //loop over entries and filter them
914 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
915 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
916 if (responseFunction==0x0) continue;
917 TString responseName=responseFunction->GetName();
918
919 if (!reg.MatchB(responseName)) continue;
920
921 TObjArray *arr=reg.MatchS(responseName);
922 TString particleName=arr->At(1)->GetName();
923 delete arr;
924 if (particleName.IsNull()) continue;
00a38d07 925 if (!grAll&&particleName=="ALL") grAll=responseFunction;
4ec8e76d 926 else {
927 //find particle id
00a38d07 928 for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
4ec8e76d 929 TString particle=AliPID::ParticleName(ispec);
930 particle.ToUpper();
931 if ( particle == particleName ){
932 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
933 fTPCResponse.SetUseDatabase(kTRUE);
934 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
00a38d07 935 // overwrite default with proton spline (for light nuclei)
936 if (ispec==AliPID::kProton) grAll=responseFunction;
4a527e08 937 found=kTRUE;
4ec8e76d 938 break;
939 }
940 }
941 }
942 }
943
944 //set default response function to all particles which don't have a specific one
945 if (grAll){
00a38d07 946 for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
4ec8e76d 947 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
948 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
949 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
4a527e08 950 found=kTRUE;
4ec8e76d 951 }
952 }
953 }
954 }
4a527e08 955
956 if (!found){
957 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
958 }
4ec8e76d 959 //
960 // Setup resolution parametrisation
961 //
962
963 //default
964 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
965
966 if (fRun>=122195){
967 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
968 }
23425eb2 969 if (fArrPidResponseMaster)
4ec8e76d 970 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
971
972 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
973}
974
ea235c90 975//______________________________________________________________________________
976void AliPIDResponse::SetTRDPidResponseMaster()
977{
978 //
979 // Load the TRD pid params and references from the OADB
980 //
db0e2c5f 981 if(fTRDPIDResponseObject) return;
53d016dc 982 AliOADBContainer contParams("contParams");
983
db0e2c5f 984 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
985 if(statusResponse){
986 AliError("Failed initializing PID Response Object from OADB");
59a8e853 987 } else {
db0e2c5f 988 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
989 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
990 if(!fTRDPIDResponseObject){
991 AliError(Form("TRD Response not found in run %d", fRun));
59a8e853 992 }
993 }
db0e2c5f 994 /*
53d016dc 995 AliOADBContainer contRefs("contRefs");
59a8e853 996 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
997 if(statusRefs){
998 AliInfo("Failed Loading References for TRD");
999 } else {
1000 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
1001 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
1002 if(!fTRDPIDReference){
1003 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
1004 }
db0e2c5f 1005 }
1006 */
ea235c90 1007}
1008
1009//______________________________________________________________________________
1010void AliPIDResponse::InitializeTRDResponse(){
1011 //
1012 // Set PID Params and references to the TRD PID response
1013 //
db0e2c5f 1014 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
f2762b1c 1015 SetTRDPIDmethod();
1016}
1017
1018void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
1019
1020 fTRDResponse.SetPIDmethod(method);
1021 if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1022 // backward compatibility for setting with 8 slices
1023 fTRDslicesForPID[0] = 0;
1024 fTRDslicesForPID[1] = 7;
1025 }
1026 else{
1027 if(method==AliTRDPIDResponse::kLQ1D){
1028 fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1029 fTRDslicesForPID[1] = 0;
1030 }
1031 if(method==AliTRDPIDResponse::kLQ2D){
1032 fTRDslicesForPID[0] = 1;
1033 fTRDslicesForPID[1] = 7;
db0e2c5f 1034 }
f2762b1c 1035 }
1036 AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
ea235c90 1037}
1038
b79db598 1039//______________________________________________________________________________
1040void AliPIDResponse::SetTOFPidResponseMaster()
1041{
1042 //
1043 // Load the TOF pid params from the OADB
1044 //
00a38d07 1045
1046 if (fTOFPIDParams) delete fTOFPIDParams;
1047 fTOFPIDParams=0x0;
1048
b79db598 1049 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
00a38d07 1050 if (oadbf && oadbf->IsOpen()) {
b79db598 1051 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1052 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
00a38d07 1053 if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 1054 oadbf->Close();
1055 delete oadbc;
b79db598 1056 }
1057 delete oadbf;
1058
00a38d07 1059 if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1060}
b79db598 1061
1062//______________________________________________________________________________
1063void AliPIDResponse::InitializeTOFResponse(){
1064 //
1065 // Set PID Params to the TOF PID response
00a38d07 1066 //
1067
1068 AliInfo("TOF PID Params loaded from OADB");
1069 AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1070 AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1071 AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1072 fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1073
b79db598 1074 for (Int_t i=0;i<4;i++) {
1075 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1076 }
1077 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1078
1079}
1080
1081
ea235c90 1082//_________________________________________________________________________
1083Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
1084 //
1085 // Check whether track is identified as electron under a given electron efficiency hypothesis
1086 //
1087 Double_t probs[AliPID::kSPECIES];
1088 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
1089
99e9d5ec 1090 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
1091 // Take mean of the TRD momenta in the given tracklets
1092 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1093 Int_t nmomenta = 0;
ea235c90 1094 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1095 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 1096 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 1097 }
1098 }
99e9d5ec 1099 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 1100
1101 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
1102}
1103
b2138b40 1104//______________________________________________________________________________
1105void AliPIDResponse::SetEMCALPidResponseMaster()
1106{
1107 //
1108 // Load the EMCAL pid response functions from the OADB
1109 //
1110 TObjArray* fEMCALPIDParamsRun = NULL;
1111 TObjArray* fEMCALPIDParamsPass = NULL;
1112
1113 if(fEMCALPIDParams) return;
1114 AliOADBContainer contParams("contParams");
1115
1116 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1117 if(statusPars){
1118 AliError("Failed initializing PID Params from OADB");
1119 }
1120 else {
1121 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1122
1123 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1124 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1125 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1126
1127 if(!fEMCALPIDParams){
f8d39067 1128 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 1129 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 1130
1f631618 1131 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1132 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 1133 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1134
1135 if(!fEMCALPIDParams){
1f631618 1136 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 1137 }
1138 }
1139 }
1140}
1141
1142//______________________________________________________________________________
1143void AliPIDResponse::InitializeEMCALResponse(){
1144 //
1145 // Set PID Params to the EMCAL PID response
1146 //
1147 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1148
1149}
00a38d07 1150
1151//_____________________________________________________
1152void AliPIDResponse::FillTrackDetectorPID()
1153{
1154 //
1155 // create detector PID information and setup the transient pointer in the track
1156 //
1157
1158 if (!fCurrentEvent) return;
1159
1160 //TODO: which particles to include? See also the loops below...
1161 Double_t values[AliPID::kSPECIESC]={0};
1162
1163 for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1164 AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1165 if (!track) continue;
1166
1167 AliDetectorPID *detPID=new AliDetectorPID;
1168 for (Int_t idet=0; idet<kNdetectors; ++idet){
1169
1170 //nsigmas
1171 for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1172 values[ipart]=NumberOfSigmas((EDetector)idet,track,(AliPID::EParticleType)ipart);
1173 detPID->SetNumberOfSigmas((EDetector)idet, values, (Int_t)AliPID::kSPECIESC);
1174
1175 //probabilities
1176 EDetPidStatus status=ComputePIDProbability((EDetector)idet,track,AliPID::kSPECIESC,values);
1177 detPID->SetRawProbability((EDetector)idet, values, (Int_t)AliPID::kSPECIESC, status);
1178 }
1179
1180 track->SetDetectorPID(detPID);
1181 }
1182}
1183
5f8db5fe 1184//_________________________________________________________________________
1185void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1186 //
1187 // Set TOF response function
1188 // Input option for event_time used
1189 //
1190
1191 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1192 if(t0spread < 10) t0spread = 80;
1193
1194 // T0 from TOF algorithm
1195
1196 Bool_t flagT0TOF=kFALSE;
1197 Bool_t flagT0T0=kFALSE;
1198 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1199 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1200 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1201
1202 // T0-TOF arrays
1203 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1204 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1205 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1206 estimatedT0event[i]=0.0;
1207 estimatedT0resolution[i]=0.0;
1208 startTimeMask[i] = 0;
1209 }
1210
1211 Float_t resT0A=75,resT0C=65,resT0AC=55;
1212 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1213 flagT0T0=kTRUE;
1214 }
1215
1216
1217 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1218
1219 if (tofHeader) { // read global info and T0-TOF
1220 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1221 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1222 if(t0spread < 10) t0spread = 80;
1223
1224 flagT0TOF=kTRUE;
1225 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1226 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1227 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1228 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1229 }
1230
1231 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1232 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1233 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1234 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1235 Int_t icurrent = (Int_t)ibin->GetAt(j);
1236 startTime[icurrent]=t0Bin->GetAt(j);
1237 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1238 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1239 }
1240 }
1241
1242 // for cut of 3 sigma on t0 spread
1243 Float_t t0cut = 3 * t0spread;
1244 if(t0cut < 500) t0cut = 500;
1245
1246 if(option == kFILL_T0){ // T0-FILL is used
1247 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1248 estimatedT0event[i]=0.0;
1249 estimatedT0resolution[i]=t0spread;
1250 }
1251 fTOFResponse.SetT0event(estimatedT0event);
1252 fTOFResponse.SetT0resolution(estimatedT0resolution);
1253 }
1254
1255 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1256 if(flagT0TOF){
1257 fTOFResponse.SetT0event(startTime);
1258 fTOFResponse.SetT0resolution(startTimeRes);
1259 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1260 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1261 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1262 }
1263 }
1264 else{
1265 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1266 estimatedT0event[i]=0.0;
1267 estimatedT0resolution[i]=t0spread;
1268 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1269 }
1270 fTOFResponse.SetT0event(estimatedT0event);
1271 fTOFResponse.SetT0resolution(estimatedT0resolution);
1272 }
1273 }
1274 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1275 Float_t t0AC=-10000;
1276 Float_t t0A=-10000;
1277 Float_t t0C=-10000;
1278 if(flagT0T0){
1279 t0AC= vevent->GetT0TOF()[0];
1280 t0A= vevent->GetT0TOF()[1];
1281 t0C= vevent->GetT0TOF()[2];
1282 }
1283
1284 Float_t t0t0Best = 0;
1285 Float_t t0t0BestRes = 9999;
1286 Int_t t0used=0;
1287 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1288 t0t0Best = t0AC;
1289 t0t0BestRes = resT0AC;
1290 t0used=6;
1291 }
1292 else if(TMath::Abs(t0C) < t0cut){
1293 t0t0Best = t0C;
1294 t0t0BestRes = resT0C;
1295 t0used=4;
1296 }
1297 else if(TMath::Abs(t0A) < t0cut){
1298 t0t0Best = t0A;
1299 t0t0BestRes = resT0A;
1300 t0used=2;
1301 }
1302
1303 if(flagT0TOF){ // if T0-TOF info is available
1304 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1305 if(t0t0BestRes < 999){
1306 if(startTimeRes[i] < t0spread){
1307 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1308 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1309 estimatedT0event[i]=t0best / wtot;
1310 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1311 startTimeMask[i] = t0used+1;
1312 }
1313 else {
1314 estimatedT0event[i]=t0t0Best;
1315 estimatedT0resolution[i]=t0t0BestRes;
1316 startTimeMask[i] = t0used;
1317 }
1318 }
1319 else{
1320 estimatedT0event[i]=startTime[i];
1321 estimatedT0resolution[i]=startTimeRes[i];
1322 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1323 }
1324 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1325 }
1326 fTOFResponse.SetT0event(estimatedT0event);
1327 fTOFResponse.SetT0resolution(estimatedT0resolution);
1328 }
1329 else{ // if no T0-TOF info is available
1330 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1331 fTOFResponse.SetT0binMask(i,t0used);
1332 if(t0t0BestRes < 999){
1333 estimatedT0event[i]=t0t0Best;
1334 estimatedT0resolution[i]=t0t0BestRes;
1335 }
1336 else{
1337 estimatedT0event[i]=0.0;
1338 estimatedT0resolution[i]=t0spread;
1339 }
1340 }
1341 fTOFResponse.SetT0event(estimatedT0event);
1342 fTOFResponse.SetT0resolution(estimatedT0resolution);
1343 }
1344 }
1345
1346 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1347 Float_t t0AC=-10000;
1348 Float_t t0A=-10000;
1349 Float_t t0C=-10000;
1350 if(flagT0T0){
1351 t0AC= vevent->GetT0TOF()[0];
1352 t0A= vevent->GetT0TOF()[1];
1353 t0C= vevent->GetT0TOF()[2];
1354 }
1355
1356 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1357 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1358 estimatedT0event[i]=t0AC;
1359 estimatedT0resolution[i]=resT0AC;
1360 fTOFResponse.SetT0binMask(i,6);
1361 }
1362 }
1363 else if(TMath::Abs(t0C) < t0cut){
1364 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1365 estimatedT0event[i]=t0C;
1366 estimatedT0resolution[i]=resT0C;
1367 fTOFResponse.SetT0binMask(i,4);
1368 }
1369 }
1370 else if(TMath::Abs(t0A) < t0cut){
1371 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1372 estimatedT0event[i]=t0A;
1373 estimatedT0resolution[i]=resT0A;
1374 fTOFResponse.SetT0binMask(i,2);
1375 }
1376 }
1377 else{
1378 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1379 estimatedT0event[i]=0.0;
1380 estimatedT0resolution[i]=t0spread;
1381 fTOFResponse.SetT0binMask(i,0);
1382 }
1383 }
1384 fTOFResponse.SetT0event(estimatedT0event);
1385 fTOFResponse.SetT0resolution(estimatedT0resolution);
1386 }
1387 delete [] startTime;
1388 delete [] startTimeRes;
1389 delete [] startTimeMask;
1390 delete [] estimatedT0event;
1391 delete [] estimatedT0resolution;
1392}