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