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