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