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