]> 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
correct also single distributions N+ and N-
[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
89 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 90}
91
92//______________________________________________________________________________
93AliPIDResponse::~AliPIDResponse()
94{
95 //
96 // dtor
97 //
00a38d07 98 delete fArrPidResponseMaster;
99 delete fTRDPIDResponseObject;
100 delete fTOFPIDParams;
4ec8e76d 101}
102
103//______________________________________________________________________________
104AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
105TNamed(other),
106fITSResponse(other.fITSResponse),
107fTPCResponse(other.fTPCResponse),
108fTRDResponse(other.fTRDResponse),
109fTOFResponse(other.fTOFResponse),
e96b9916 110fEMCALResponse(other.fEMCALResponse),
fd21ec8d 111fRange(other.fRange),
112fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 113fIsMC(other.fIsMC),
114fOADBPath(other.fOADBPath),
00a38d07 115fCustomTPCpidResponse(other.fCustomTPCpidResponse),
4ec8e76d 116fBeamType("PP"),
117fLHCperiod(),
118fMCperiodTPC(),
fd21ec8d 119fMCperiodUser(other.fMCperiodUser),
ea235c90 120fCurrentFile(),
4ec8e76d 121fRecoPass(0),
fd21ec8d 122fRecoPassUser(other.fRecoPassUser),
4ec8e76d 123fRun(0),
124fOldRun(0),
644666df 125fArrPidResponseMaster(NULL),
126fResolutionCorrection(NULL),
127fOADBvoltageMaps(NULL),
128fTRDPIDResponseObject(NULL),
0b39f221 129fTOFtail(1.1),
644666df 130fTOFPIDParams(NULL),
131fEMCALPIDParams(NULL),
132fCurrentEvent(NULL),
539a5a59 133fCurrCentrality(0.0),
134fTuneMConData(kFALSE)
4ec8e76d 135{
136 //
137 // copy ctor
138 //
ea235c90 139 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 140}
141
142//______________________________________________________________________________
143AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
144{
145 //
146 // copy ctor
147 //
148 if(this!=&other) {
149 delete fArrPidResponseMaster;
150 TNamed::operator=(other);
151 fITSResponse=other.fITSResponse;
152 fTPCResponse=other.fTPCResponse;
153 fTRDResponse=other.fTRDResponse;
154 fTOFResponse=other.fTOFResponse;
e96b9916 155 fEMCALResponse=other.fEMCALResponse;
fd21ec8d 156 fRange=other.fRange;
157 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 158 fOADBPath=other.fOADBPath;
00a38d07 159 fCustomTPCpidResponse=other.fCustomTPCpidResponse;
4ec8e76d 160 fIsMC=other.fIsMC;
161 fBeamType="PP";
162 fLHCperiod="";
163 fMCperiodTPC="";
fd21ec8d 164 fMCperiodUser=other.fMCperiodUser;
ea235c90 165 fCurrentFile="";
4ec8e76d 166 fRecoPass=0;
fd21ec8d 167 fRecoPassUser=other.fRecoPassUser;
4ec8e76d 168 fRun=0;
169 fOldRun=0;
644666df 170 fArrPidResponseMaster=NULL;
171 fResolutionCorrection=NULL;
172 fOADBvoltageMaps=NULL;
173 fTRDPIDResponseObject=NULL;
174 fEMCALPIDParams=NULL;
ea235c90 175 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
0b39f221 176 fTOFtail=1.1;
644666df 177 fTOFPIDParams=NULL;
e96b9916 178 fCurrentEvent=other.fCurrentEvent;
4ec8e76d 179 }
180 return *this;
181}
182
fd21ec8d 183//______________________________________________________________________________
184Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
185{
186 //
187 // NumberOfSigmas for 'detCode'
188 //
189
190 switch (detCode){
191 case kDetITS: return NumberOfSigmasITS(track, type); break;
192 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
193 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
00a38d07 194 case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
fd21ec8d 195 default: return -999.;
196 }
197
198}
199
e96b9916 200//______________________________________________________________________________
00a38d07 201Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
202{
203 //
204 // NumberOfSigmas for 'detCode'
205 //
206 return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
207}
208
209//______________________________________________________________________________
210Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
211{
212 //
213 // Calculate the number of sigmas in the ITS
214 //
215
216 AliVTrack *track=(AliVTrack*)vtrack;
217
218 // look for cached value first
219 // only the non SA tracks are cached
220 if ( track->GetDetectorPID() ){
221 return track->GetDetectorPID()->GetNumberOfSigmas(kITS, type);
222 }
223
224 Float_t dEdx=track->GetITSsignal();
225 if (dEdx<=0) return -999.;
226
227 UChar_t clumap=track->GetITSClusterMap();
228 Int_t nPointsForPid=0;
229 for(Int_t i=2; i<6; i++){
230 if(clumap&(1<<i)) ++nPointsForPid;
231 }
232 Float_t mom=track->P();
233
234 //check for ITS standalone tracks
235 Bool_t isSA=kTRUE;
236 if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
237
238 //TODO: in case of the electron, use the SA parametrisation,
239 // this needs to be changed if ITS provides a parametrisation
240 // for electrons also for ITS+TPC tracks
241 return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
242}
243
244//______________________________________________________________________________
245Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
246{
247 //
248 // Calculate the number of sigmas in the TPC
249 //
250
251 AliVTrack *track=(AliVTrack*)vtrack;
252
253 // look for cached value first
254 if (track->GetDetectorPID()){
255 return track->GetDetectorPID()->GetNumberOfSigmas(kTPC, type);
256 }
257
258 Double_t mom = track->GetTPCmomentum();
259 Double_t sig = track->GetTPCsignal();
260 UInt_t sigN = track->GetTPCsignalN();
261
262 Double_t nSigma = -999.;
263 if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
264
265 return nSigma;
266}
267
644666df 268//______________________________________________________________________________
269Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
270 AliPID::EParticleType type,
271 AliTPCPIDResponse::ETPCdEdxSource dedxSource)
272{
273 //get number of sigmas according the selected TPC gain configuration scenario
274 const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
275
276 Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource);
277
278 return nSigma;
279}
280
00a38d07 281//______________________________________________________________________________
282Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
283{
284 //
285 // Calculate the number of sigmas in the EMCAL
286 //
287
288 AliVTrack *track=(AliVTrack*)vtrack;
e96b9916 289
00a38d07 290 // look for cached value first
291 if (track->GetDetectorPID()){
292 return track->GetDetectorPID()->GetNumberOfSigmas(kEMCAL, type);
293 }
294
e96b9916 295 AliVCluster *matchedClus = NULL;
296
297 Double_t mom = -1.;
298 Double_t pt = -1.;
299 Double_t EovP = -1.;
300 Double_t fClsE = -1.;
301
302 Int_t nMatchClus = -1;
303 Int_t charge = 0;
304
305 // Track matching
306 nMatchClus = track->GetEMCALcluster();
307 if(nMatchClus > -1){
6d0064aa 308
e96b9916 309 mom = track->P();
310 pt = track->Pt();
311 charge = track->Charge();
312
313 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
314
315 if(matchedClus){
316
6d0064aa 317 // matched cluster is EMCAL
318 if(matchedClus->IsEMCAL()){
319
320 fClsE = matchedClus->E();
321 EovP = fClsE/mom;
322
323
324 // NSigma value really meaningful only for electrons!
325 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
326 }
e96b9916 327 }
328 }
6d0064aa 329
e96b9916 330 return -999;
331
332}
333
6d0064aa 334//______________________________________________________________________________
00a38d07 335Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
6d0064aa 336
00a38d07 337 AliVTrack *track=(AliVTrack*)vtrack;
338
6d0064aa 339 AliVCluster *matchedClus = NULL;
340
341 Double_t mom = -1.;
342 Double_t pt = -1.;
343 Double_t EovP = -1.;
344 Double_t fClsE = -1.;
32fa24d6 345
346 // initialize eop and shower shape parameters
347 eop = -1.;
348 for(Int_t i = 0; i < 4; i++){
349 showershape[i] = -1.;
350 }
6d0064aa 351
352 Int_t nMatchClus = -1;
353 Int_t charge = 0;
354
355 // Track matching
356 nMatchClus = track->GetEMCALcluster();
357 if(nMatchClus > -1){
358
359 mom = track->P();
360 pt = track->Pt();
361 charge = track->Charge();
362
363 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
364
365 if(matchedClus){
366
367 // matched cluster is EMCAL
368 if(matchedClus->IsEMCAL()){
369
370 fClsE = matchedClus->E();
371 EovP = fClsE/mom;
372
373 // fill used EMCAL variables here
374 eop = EovP; // E/p
375 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
376 showershape[1] = matchedClus->GetM02(); // long axis
377 showershape[2] = matchedClus->GetM20(); // short axis
378 showershape[3] = matchedClus->GetDispersion(); // dispersion
379
380 // NSigma value really meaningful only for electrons!
381 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
382 }
383 }
384 }
385 return -999;
386}
387
388
fd21ec8d 389//______________________________________________________________________________
390AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
391{
392 //
393 // Compute PID response of 'detCode'
394 //
395
396 switch (detCode){
397 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
398 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
fd21ec8d 399 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
00a38d07 400 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
fd21ec8d 401 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
402 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
403 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
404 default: return kDetNoSignal;
405 }
406}
407
00a38d07 408//______________________________________________________________________________
409AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
410{
411 //
412 // Compute PID response of 'detCode'
413 //
414
415 return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
416}
417
fd21ec8d 418//______________________________________________________________________________
419AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
420{
421 //
422 // Compute PID response for the ITS
423 //
424
00a38d07 425 // look for cached value first
426 // only the non SA tracks are cached
427 if (track->GetDetectorPID()){
428 return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
429 }
430
fd21ec8d 431 // set flat distribution (no decision)
432 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
433
434 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
435 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
00a38d07 436
437 //check for ITS standalone tracks
438 Bool_t isSA=kTRUE;
439 if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
fd21ec8d 440
441 Double_t mom=track->P();
442 Double_t dedx=track->GetITSsignal();
fd21ec8d 443 Double_t momITS=mom;
fd21ec8d 444 UChar_t clumap=track->GetITSClusterMap();
445 Int_t nPointsForPid=0;
446 for(Int_t i=2; i<6; i++){
447 if(clumap&(1<<i)) ++nPointsForPid;
448 }
449
450 if(nPointsForPid<3) { // track not to be used for combined PID purposes
451 // track->ResetStatus(AliVTrack::kITSpid);
452 return kDetNoSignal;
453 }
454
455 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
456 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
00a38d07 457 Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
458 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
459 Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
460 //TODO: in case of the electron, use the SA parametrisation,
461 // this needs to be changed if ITS provides a parametrisation
462 // for electrons also for ITS+TPC tracks
463 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
fd21ec8d 464 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
465 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
466 } else {
467 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
468 mismatch=kFALSE;
469 }
470
471 // Check for particles heavier than (AliPID::kSPECIES - 1)
472 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
473
474 }
475
476 if (mismatch){
477 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
478 return kDetNoSignal;
479 }
480
481
482 return kDetPidOk;
483}
484//______________________________________________________________________________
485AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
486{
487 //
488 // Compute PID response for the TPC
489 //
00a38d07 490
491 // look for cached value first
492 if (track->GetDetectorPID()){
493 return track->GetDetectorPID()->GetRawProbability(kTPC, p, nSpecies);
494 }
495
fd21ec8d 496 // set flat distribution (no decision)
497 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
498
499 // check quality of the track
500 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
501
502 Double_t mom = track->GetTPCmomentum();
503
504 Double_t dedx=track->GetTPCsignal();
505 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
506
539a5a59 507 if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
508
00a38d07 509 for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
fd21ec8d 510 AliPID::EParticleType type=AliPID::EParticleType(j);
511 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
512 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
513 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
514 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
515 } else {
516 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
517 mismatch=kFALSE;
518 }
fd21ec8d 519 }
520
521 if (mismatch){
522 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
523 return kDetNoSignal;
524 }
525
526 return kDetPidOk;
527}
528//______________________________________________________________________________
529AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
530{
531 //
532 // Compute PID response for the
533 //
00a38d07 534
535 // look for cached value first
536 if (track->GetDetectorPID()){
537 return track->GetDetectorPID()->GetRawProbability(kTOF, p, nSpecies);
538 }
539
0b39f221 540 Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
541
fd21ec8d 542 // set flat distribution (no decision)
543 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
544
545 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
546 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
547
fd21ec8d 548 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
00a38d07 549 for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
fd21ec8d 550 AliPID::EParticleType type=AliPID::EParticleType(j);
0b39f221 551 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
552
67376d1d 553 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
554 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
0b39f221 555 if (TMath::Abs(nsigmas) > (fRange+2)) {
556 if(nsigmas < fTOFtail)
557 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
558 else
559 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
560 } else{
561 if(nsigmas < fTOFtail)
562 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
563 else
564 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
565 }
566
fd21ec8d 567 if (TMath::Abs(nsigmas)<5.){
568 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
569 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
570 }
571 }
572
573 if (mismatch){
574 return kDetMismatch;
575 }
576
577 // TODO: Light nuclei
578
579 return kDetPidOk;
580}
581//______________________________________________________________________________
ea235c90 582AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 583{
584 //
585 // Compute PID response for the
586 //
00a38d07 587
588 // look for cached value first
589 if (track->GetDetectorPID()){
590 return track->GetDetectorPID()->GetRawProbability(kTRD, p, nSpecies);
591 }
592
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
ea235c90 599 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
600 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
601 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
602 mom[ilayer] = track->GetTRDmomentum(ilayer);
603 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
604 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
605 }
606 }
607 fTRDResponse.GetResponse(nslices, dedx, mom, p);
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 SetTRDPIDmethod();
1088}
1089
1090void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
1091
1092 fTRDResponse.SetPIDmethod(method);
1093 if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
1094 // backward compatibility for setting with 8 slices
1095 fTRDslicesForPID[0] = 0;
1096 fTRDslicesForPID[1] = 7;
1097 }
1098 else{
1099 if(method==AliTRDPIDResponse::kLQ1D){
1100 fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1101 fTRDslicesForPID[1] = 0;
1102 }
1103 if(method==AliTRDPIDResponse::kLQ2D){
1104 fTRDslicesForPID[0] = 1;
1105 fTRDslicesForPID[1] = 7;
db0e2c5f 1106 }
f2762b1c 1107 }
1108 AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
ea235c90 1109}
1110
b79db598 1111//______________________________________________________________________________
1112void AliPIDResponse::SetTOFPidResponseMaster()
1113{
1114 //
1115 // Load the TOF pid params from the OADB
1116 //
00a38d07 1117
1118 if (fTOFPIDParams) delete fTOFPIDParams;
644666df 1119 fTOFPIDParams=NULL;
00a38d07 1120
b79db598 1121 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
00a38d07 1122 if (oadbf && oadbf->IsOpen()) {
b79db598 1123 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1124 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
00a38d07 1125 if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 1126 oadbf->Close();
1127 delete oadbc;
b79db598 1128 }
1129 delete oadbf;
1130
00a38d07 1131 if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1132}
b79db598 1133
1134//______________________________________________________________________________
1135void AliPIDResponse::InitializeTOFResponse(){
1136 //
1137 // Set PID Params to the TOF PID response
00a38d07 1138 //
1139
1140 AliInfo("TOF PID Params loaded from OADB");
1141 AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1142 AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1143 AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1144 fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
1145
b79db598 1146 for (Int_t i=0;i<4;i++) {
1147 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1148 }
1149 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1150
1151}
1152
1153
ea235c90 1154//_________________________________________________________________________
1155Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
1156 //
1157 // Check whether track is identified as electron under a given electron efficiency hypothesis
1158 //
1159 Double_t probs[AliPID::kSPECIES];
1160 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
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
1173 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
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}