]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliPIDResponse.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update timestamp for new data points simulation
[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>
f84b18dd 30#include <TH2D.h>
4ec8e76d 31#include <TSpline.h>
32#include <TFile.h>
00a38d07 33#include <TArrayI.h>
db0e2c5f 34#include <TArrayF.h>
f84b18dd 35#include <TLinearFitter.h>
5a9dc560 36#include <TSystem.h>
37#include <TMD5.h>
4ec8e76d 38
39#include <AliVEvent.h>
fd21ec8d 40#include <AliVTrack.h>
4ec8e76d 41#include <AliLog.h>
42#include <AliPID.h>
ea235c90 43#include <AliOADBContainer.h>
db0e2c5f 44#include <AliTRDPIDResponseObject.h>
9c499471 45#include <AliTRDdEdxParams.h>
b79db598 46#include <AliTOFPIDParams.h>
567624b5 47#include <AliHMPIDPIDParams.h>
29bf19f2 48
49#include "AliPIDResponse.h"
00a38d07 50#include "AliDetectorPID.h"
29bf19f2 51
80f28562 52#include "AliCentrality.h"
53
29bf19f2 54ClassImp(AliPIDResponse);
55
42fcc729 56Float_t AliPIDResponse::fgTOFmismatchProb = 0.0;
57
4ec8e76d 58AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
59TNamed("PIDResponse","PIDResponse"),
60fITSResponse(isMC),
61fTPCResponse(),
62fTRDResponse(),
63fTOFResponse(),
567624b5 64fHMPIDResponse(),
e96b9916 65fEMCALResponse(),
fd21ec8d 66fRange(5.),
67fITSPIDmethod(kITSTruncMean),
1ceae0ac 68fTuneMConData(kFALSE),
69fTuneMConDataMask(kDetTOF|kDetTPC),
4ec8e76d 70fIsMC(isMC),
a017c06a 71fCachePID(kFALSE),
4ec8e76d 72fOADBPath(),
00a38d07 73fCustomTPCpidResponse(),
d5113358 74fCustomTPCetaMaps(),
4ec8e76d 75fBeamType("PP"),
76fLHCperiod(),
77fMCperiodTPC(),
fd21ec8d 78fMCperiodUser(),
ea235c90 79fCurrentFile(),
87da0205 80fCurrentAliRootRev(-1),
4ec8e76d 81fRecoPass(0),
fd21ec8d 82fRecoPassUser(-1),
1ceae0ac 83fRun(-1),
1b9e31a6 84fOldRun(-1),
78cbd205 85fResT0A(75.),
86fResT0C(65.),
87fResT0AC(55.),
644666df 88fArrPidResponseMaster(NULL),
89fResolutionCorrection(NULL),
90fOADBvoltageMaps(NULL),
87da0205 91fUseTPCEtaCorrection(kFALSE),
92fUseTPCMultiplicityCorrection(kFALSE),
644666df 93fTRDPIDResponseObject(NULL),
9c499471 94fTRDdEdxParams(NULL),
c5fb644a 95fTOFtail(0.9),
644666df 96fTOFPIDParams(NULL),
567624b5 97fHMPIDPIDParams(NULL),
644666df 98fEMCALPIDParams(NULL),
99fCurrentEvent(NULL),
42fcc729 100fCurrCentrality(0.0),
6839b659 101fBeamTypeNum(kPP),
102fNoTOFmism(kFALSE)
4ec8e76d 103{
104 //
105 // default ctor
106 //
a635821f 107 AliLog::SetClassDebugLevel("AliPIDResponse",0);
108 AliLog::SetClassDebugLevel("AliESDpid",0);
109 AliLog::SetClassDebugLevel("AliAODpidUtil",0);
ea235c90 110
4ec8e76d 111}
112
113//______________________________________________________________________________
114AliPIDResponse::~AliPIDResponse()
115{
116 //
117 // dtor
118 //
00a38d07 119 delete fArrPidResponseMaster;
120 delete fTRDPIDResponseObject;
9c499471 121 delete fTRDdEdxParams;
00a38d07 122 delete fTOFPIDParams;
4ec8e76d 123}
124
125//______________________________________________________________________________
126AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
127TNamed(other),
128fITSResponse(other.fITSResponse),
129fTPCResponse(other.fTPCResponse),
130fTRDResponse(other.fTRDResponse),
131fTOFResponse(other.fTOFResponse),
567624b5 132fHMPIDResponse(other.fHMPIDResponse),
e96b9916 133fEMCALResponse(other.fEMCALResponse),
fd21ec8d 134fRange(other.fRange),
135fITSPIDmethod(other.fITSPIDmethod),
1ceae0ac 136fTuneMConData(other.fTuneMConData),
137fTuneMConDataMask(other.fTuneMConDataMask),
4ec8e76d 138fIsMC(other.fIsMC),
1c9d11be 139fCachePID(other.fCachePID),
4ec8e76d 140fOADBPath(other.fOADBPath),
00a38d07 141fCustomTPCpidResponse(other.fCustomTPCpidResponse),
d5113358 142fCustomTPCetaMaps(other.fCustomTPCetaMaps),
4ec8e76d 143fBeamType("PP"),
144fLHCperiod(),
145fMCperiodTPC(),
fd21ec8d 146fMCperiodUser(other.fMCperiodUser),
ea235c90 147fCurrentFile(),
87da0205 148fCurrentAliRootRev(other.fCurrentAliRootRev),
4ec8e76d 149fRecoPass(0),
fd21ec8d 150fRecoPassUser(other.fRecoPassUser),
1ceae0ac 151fRun(-1),
1b9e31a6 152fOldRun(-1),
78cbd205 153fResT0A(75.),
154fResT0C(65.),
155fResT0AC(55.),
644666df 156fArrPidResponseMaster(NULL),
157fResolutionCorrection(NULL),
158fOADBvoltageMaps(NULL),
f84b18dd 159fUseTPCEtaCorrection(other.fUseTPCEtaCorrection),
87da0205 160fUseTPCMultiplicityCorrection(other.fUseTPCMultiplicityCorrection),
644666df 161fTRDPIDResponseObject(NULL),
9c499471 162fTRDdEdxParams(NULL),
c5fb644a 163fTOFtail(0.9),
644666df 164fTOFPIDParams(NULL),
567624b5 165fHMPIDPIDParams(NULL),
644666df 166fEMCALPIDParams(NULL),
167fCurrentEvent(NULL),
42fcc729 168fCurrCentrality(0.0),
6839b659 169fBeamTypeNum(kPP),
170fNoTOFmism(other.fNoTOFmism)
4ec8e76d 171{
172 //
173 // copy ctor
174 //
175}
176
177//______________________________________________________________________________
178AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
179{
180 //
181 // copy ctor
182 //
183 if(this!=&other) {
184 delete fArrPidResponseMaster;
185 TNamed::operator=(other);
186 fITSResponse=other.fITSResponse;
187 fTPCResponse=other.fTPCResponse;
188 fTRDResponse=other.fTRDResponse;
189 fTOFResponse=other.fTOFResponse;
567624b5 190 fHMPIDResponse=other.fHMPIDResponse;
e96b9916 191 fEMCALResponse=other.fEMCALResponse;
fd21ec8d 192 fRange=other.fRange;
193 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 194 fOADBPath=other.fOADBPath;
00a38d07 195 fCustomTPCpidResponse=other.fCustomTPCpidResponse;
d5113358 196 fCustomTPCetaMaps=other.fCustomTPCetaMaps;
1ceae0ac 197 fTuneMConData=other.fTuneMConData;
198 fTuneMConDataMask=other.fTuneMConDataMask;
4ec8e76d 199 fIsMC=other.fIsMC;
1c9d11be 200 fCachePID=other.fCachePID;
4ec8e76d 201 fBeamType="PP";
42fcc729 202 fBeamTypeNum=kPP;
4ec8e76d 203 fLHCperiod="";
204 fMCperiodTPC="";
fd21ec8d 205 fMCperiodUser=other.fMCperiodUser;
ea235c90 206 fCurrentFile="";
87da0205 207 fCurrentAliRootRev=other.fCurrentAliRootRev;
4ec8e76d 208 fRecoPass=0;
fd21ec8d 209 fRecoPassUser=other.fRecoPassUser;
1ceae0ac 210 fRun=-1;
1b9e31a6 211 fOldRun=-1;
78cbd205 212 fResT0A=75.;
213 fResT0C=65.;
214 fResT0AC=55.;
644666df 215 fArrPidResponseMaster=NULL;
216 fResolutionCorrection=NULL;
217 fOADBvoltageMaps=NULL;
1ceae0ac 218 fUseTPCEtaCorrection=other.fUseTPCEtaCorrection;
87da0205 219 fUseTPCMultiplicityCorrection=other.fUseTPCMultiplicityCorrection;
644666df 220 fTRDPIDResponseObject=NULL;
9c499471 221 fTRDdEdxParams=NULL;
644666df 222 fEMCALPIDParams=NULL;
c5fb644a 223 fTOFtail=0.9;
644666df 224 fTOFPIDParams=NULL;
567624b5 225 fHMPIDPIDParams=NULL;
e96b9916 226 fCurrentEvent=other.fCurrentEvent;
6839b659 227 fNoTOFmism = other.fNoTOFmism;
87da0205 228
4ec8e76d 229 }
230 return *this;
231}
232
fd21ec8d 233//______________________________________________________________________________
355b831b 234Float_t AliPIDResponse::NumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
fd21ec8d 235{
236 //
237 // NumberOfSigmas for 'detCode'
238 //
355b831b 239
240 const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
241 // look for cached value first
242 const AliDetectorPID *detPID=track->GetDetectorPID();
243
244 if ( detPID && detPID->HasNumberOfSigmas(detector)){
245 return detPID->GetNumberOfSigmas(detector, type);
246 } else if (fCachePID) {
247 FillTrackDetectorPID(track, detector);
248 detPID=track->GetDetectorPID();
249 return detPID->GetNumberOfSigmas(detector, type);
fd21ec8d 250 }
355b831b 251
252 return GetNumberOfSigmas(detector, track, type);
fd21ec8d 253}
254
e96b9916 255//______________________________________________________________________________
355b831b 256AliPIDResponse::EDetPidStatus AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track,
257 AliPID::EParticleType type, Double_t &val) const
00a38d07 258{
259 //
355b831b 260 // NumberOfSigmas with detector status as return value
00a38d07 261 //
355b831b 262
263 val=NumberOfSigmas(detCode, track, type);
264 return CheckPIDStatus(detCode, (AliVTrack*)track);
00a38d07 265}
266
1c9d11be 267//______________________________________________________________________________
268// public buffered versions of the PID calculation
269//
270
00a38d07 271//______________________________________________________________________________
272Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
273{
274 //
275 // Calculate the number of sigmas in the ITS
276 //
277
355b831b 278 return NumberOfSigmas(kITS, vtrack, type);
00a38d07 279}
280
281//______________________________________________________________________________
282Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
283{
284 //
285 // Calculate the number of sigmas in the TPC
286 //
287
355b831b 288 return NumberOfSigmas(kTPC, vtrack, type);
00a38d07 289}
290
644666df 291//______________________________________________________________________________
292Float_t AliPIDResponse::NumberOfSigmasTPC( const AliVParticle *vtrack,
293 AliPID::EParticleType type,
f84b18dd 294 AliTPCPIDResponse::ETPCdEdxSource dedxSource) const
644666df 295{
296 //get number of sigmas according the selected TPC gain configuration scenario
297 const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
298
87da0205 299 Float_t nSigma=fTPCResponse.GetNumberOfSigmas(track, type, dedxSource, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
644666df 300
301 return nSigma;
302}
303
5cd0300d 304//______________________________________________________________________________
305Float_t AliPIDResponse::NumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
306{
307 //
308 // Calculate the number of sigmas in the TRD
309 //
310 return NumberOfSigmas(kTRD, vtrack, type);
311}
312
00a38d07 313//______________________________________________________________________________
1c9d11be 314Float_t AliPIDResponse::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
00a38d07 315{
316 //
1c9d11be 317 // Calculate the number of sigmas in the TOF
00a38d07 318 //
319
355b831b 320 return NumberOfSigmas(kTOF, vtrack, type);
1c9d11be 321}
e96b9916 322
567624b5 323//______________________________________________________________________________
324Float_t AliPIDResponse::NumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
325{
326 //
327 // Calculate the number of sigmas in the EMCAL
328 //
329
330 return NumberOfSigmas(kHMPID, vtrack, type);
331}
332
1c9d11be 333//______________________________________________________________________________
334Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
335{
336 //
337 // Calculate the number of sigmas in the EMCAL
338 //
e96b9916 339
355b831b 340 return NumberOfSigmas(kEMCAL, vtrack, type);
e96b9916 341}
342
6d0064aa 343//______________________________________________________________________________
1c9d11be 344Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const
345{
346 //
347 // emcal nsigma with eop and showershape
348 //
00a38d07 349 AliVTrack *track=(AliVTrack*)vtrack;
350
6d0064aa 351 AliVCluster *matchedClus = NULL;
352
353 Double_t mom = -1.;
354 Double_t pt = -1.;
355 Double_t EovP = -1.;
356 Double_t fClsE = -1.;
32fa24d6 357
358 // initialize eop and shower shape parameters
359 eop = -1.;
360 for(Int_t i = 0; i < 4; i++){
361 showershape[i] = -1.;
362 }
6d0064aa 363
364 Int_t nMatchClus = -1;
365 Int_t charge = 0;
366
367 // Track matching
368 nMatchClus = track->GetEMCALcluster();
369 if(nMatchClus > -1){
370
371 mom = track->P();
372 pt = track->Pt();
373 charge = track->Charge();
374
375 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
376
377 if(matchedClus){
378
379 // matched cluster is EMCAL
380 if(matchedClus->IsEMCAL()){
381
382 fClsE = matchedClus->E();
383 EovP = fClsE/mom;
384
385 // fill used EMCAL variables here
386 eop = EovP; // E/p
387 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
388 showershape[1] = matchedClus->GetM02(); // long axis
389 showershape[2] = matchedClus->GetM20(); // short axis
390 showershape[3] = matchedClus->GetDispersion(); // dispersion
1c9d11be 391
392 // look for cached value first
393 const AliDetectorPID *detPID=track->GetDetectorPID();
394 const EDetector detector=kEMCAL;
395
396 if ( detPID && detPID->HasNumberOfSigmas(detector)){
397 return detPID->GetNumberOfSigmas(detector, type);
398 } else if (fCachePID) {
399 FillTrackDetectorPID(track, detector);
400 detPID=track->GetDetectorPID();
401 return detPID->GetNumberOfSigmas(detector, type);
402 }
403
404 // NSigma value really meaningful only for electrons!
405 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
6d0064aa 406 }
407 }
408 }
409 return -999;
410}
411
567624b5 412//______________________________________________________________________________
1d59271b 413AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDelta(EDetector detector, const AliVParticle *track, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
567624b5 414{
415 //
416 //
417 //
418 val=-9999.;
419 switch (detector){
1d59271b 420 case kITS: return GetSignalDeltaITS(track,type,val,ratio); break;
421 case kTPC: return GetSignalDeltaTPC(track,type,val,ratio); break;
5cd0300d 422 case kTRD: return GetSignalDeltaTRD(track,type,val,ratio); break;
1d59271b 423 case kTOF: return GetSignalDeltaTOF(track,type,val,ratio); break;
424 case kHMPID: return GetSignalDeltaHMPID(track,type,val,ratio); break;
567624b5 425 default: return kDetNoSignal;
426 }
427 return kDetNoSignal;
428}
429
430//______________________________________________________________________________
1d59271b 431Double_t AliPIDResponse::GetSignalDelta(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type, Bool_t ratio/*=kFALSE*/) const
567624b5 432{
433 //
434 //
435 //
436 Double_t val=-9999.;
1d59271b 437 EDetPidStatus stat=GetSignalDelta(detCode, track, type, val, ratio);
567624b5 438 if ( stat==kDetNoSignal ) val=-9999.;
439 return val;
440}
441
00a38d07 442//______________________________________________________________________________
355b831b 443AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
00a38d07 444{
00a38d07 445 // Compute PID response of 'detCode'
355b831b 446
447 // find detector code from detector bit mask
448 Int_t detector=-1;
449 for (Int_t idet=0; idet<kNdetectors; ++idet) if ( (detCode&(1<<idet)) ) { detector=idet; break; }
450 if (detector==-1) return kDetNoSignal;
00a38d07 451
355b831b 452 return ComputePIDProbability((EDetector)detector, track, nSpecies, p);
00a38d07 453}
454
fd21ec8d 455//______________________________________________________________________________
355b831b 456AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetector detector, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 457{
458 //
355b831b 459 // Compute PID response of 'detector'
fd21ec8d 460 //
461
1c9d11be 462 const AliDetectorPID *detPID=track->GetDetectorPID();
355b831b 463
464 if ( detPID && detPID->HasRawProbability(detector)){
1c9d11be 465 return detPID->GetRawProbability(detector, p, nSpecies);
466 } else if (fCachePID) {
467 FillTrackDetectorPID(track, detector);
468 detPID=track->GetDetectorPID();
469 return detPID->GetRawProbability(detector, p, nSpecies);
fd21ec8d 470 }
471
355b831b 472 //if no caching return values calculated from scratch
473 return GetComputePIDProbability(detector, track, nSpecies, p);
474}
475
476//______________________________________________________________________________
477AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
478{
479 // Compute PID response for the ITS
480 return ComputePIDProbability(kITS, track, nSpecies, p);
fd21ec8d 481}
355b831b 482
fd21ec8d 483//______________________________________________________________________________
484AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
485{
fd21ec8d 486 // Compute PID response for the TPC
355b831b 487 return ComputePIDProbability(kTPC, track, nSpecies, p);
fd21ec8d 488}
355b831b 489
fd21ec8d 490//______________________________________________________________________________
491AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
492{
fd21ec8d 493 // Compute PID response for the
355b831b 494 return ComputePIDProbability(kTOF, track, nSpecies, p);
fd21ec8d 495}
355b831b 496
fd21ec8d 497//______________________________________________________________________________
355b831b 498AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 499{
fd21ec8d 500 // Compute PID response for the
355b831b 501 return ComputePIDProbability(kTRD, track, nSpecies, p);
fd21ec8d 502}
355b831b 503
fd21ec8d 504//______________________________________________________________________________
e96b9916 505AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 506{
fd21ec8d 507 // Compute PID response for the EMCAL
355b831b 508 return ComputePIDProbability(kEMCAL, track, nSpecies, p);
1c9d11be 509}
510//______________________________________________________________________________
511AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
512{
1c9d11be 513 // Compute PID response for the PHOS
1c9d11be 514
515 // set flat distribution (no decision)
00a38d07 516 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1c9d11be 517 return kDetNoSignal;
518}
355b831b 519
1c9d11be 520//______________________________________________________________________________
521AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
522{
1c9d11be 523 // Compute PID response for the HMPID
355b831b 524 return ComputePIDProbability(kHMPID, track, nSpecies, p);
525}
fd21ec8d 526
355b831b 527//______________________________________________________________________________
528AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
529{
530 // Compute PID response for the
531 return GetComputeTRDProbability(track, nSpecies, p, PIDmethod);
532}
e96b9916 533
355b831b 534//______________________________________________________________________________
535AliPIDResponse::EDetPidStatus AliPIDResponse::CheckPIDStatus(EDetector detector, const AliVTrack *track) const
536{
537 // calculate detector pid status
538
539 const Int_t iDetCode=(Int_t)detector;
540 if (iDetCode<0||iDetCode>=kNdetectors) return kDetNoSignal;
541 const AliDetectorPID *detPID=track->GetDetectorPID();
542
543 if ( detPID ){
544 return detPID->GetPIDStatus(detector);
1c9d11be 545 } else if (fCachePID) {
546 FillTrackDetectorPID(track, detector);
547 detPID=track->GetDetectorPID();
355b831b 548 return detPID->GetPIDStatus(detector);
e96b9916 549 }
355b831b 550
551 // if not buffered and no buffering is requested
552 return GetPIDStatus(detector, track);
fd21ec8d 553}
554
4ec8e76d 555//______________________________________________________________________________
00a38d07 556void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
4ec8e76d 557{
558 //
559 // Apply settings for the current event
560 //
561 fRecoPass=pass;
e96b9916 562
78cbd205 563
644666df 564 fCurrentEvent=NULL;
4ec8e76d 565 if (!event) return;
e96b9916 566 fCurrentEvent=event;
00a38d07 567 if (run>0) fRun=run;
568 else fRun=event->GetRunNumber();
4ec8e76d 569
570 if (fRun!=fOldRun){
571 ExecNewRun();
572 fOldRun=fRun;
573 }
574
575 //TPC resolution parametrisation PbPb
576 if ( fResolutionCorrection ){
577 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
578 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
579 }
580
87da0205 581 // Set up TPC multiplicity for PbPb
22158469
JW		 582 if (fUseTPCMultiplicityCorrection) {
583 Int_t numESDtracks = event->GetNumberOfESDTracks();
584 if (numESDtracks < 0) {
585 AliError("Cannot obtain event multiplicity (number of ESD tracks < 0). If you are using AODs, this might be a too old production. Please disable the multiplicity correction to get a reliable PID result!");
586 numESDtracks = 0;
587 }
588 fTPCResponse.SetCurrentEventMultiplicity(numESDtracks);
87da0205 589 }
87da0205 590 else
591 fTPCResponse.SetCurrentEventMultiplicity(0);
592
4ec8e76d 593 //TOF resolution
b79db598 594 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
595
80f28562 596
597 // Get and set centrality
598 AliCentrality *centrality = event->GetCentrality();
599 if(centrality){
600 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
601 }
602 else{
603 fCurrCentrality = -1;
604 }
87da0205 605
606 // Set centrality percentile for EMCAL
607 fEMCALResponse.SetCentrality(fCurrCentrality);
608
c53e310b 609 // switch off some TOF channel according to OADB to match data TOF matching eff
610 if (fTuneMConData && ((fTuneMConDataMask & kDetTOF) == kDetTOF) && fTOFPIDParams->GetTOFmatchingLossMC() > 0.01){
611 Int_t ntrk = event->GetNumberOfTracks();
612 for(Int_t i=0;i < ntrk;i++){
613 AliVParticle *trk = event->GetTrack(i);
614 Int_t channel = GetTOFResponse().GetTOFchannel(trk);
615 Int_t swoffEachOfThem = Int_t(100./fTOFPIDParams->GetTOFmatchingLossMC() + 0.5);
616 if(!(channel%swoffEachOfThem)) ((AliVTrack *) trk)->ResetStatus(AliVTrack::kTOFout);
617 }
618 }
619
4ec8e76d 620}
621
622//______________________________________________________________________________
623void AliPIDResponse::ExecNewRun()
624{
625 //
626 // Things to Execute upon a new run
627 //
628 SetRecoInfo();
629
630 SetITSParametrisation();
631
632 SetTPCPidResponseMaster();
633 SetTPCParametrisation();
f84b18dd 634 SetTPCEtaMaps();
53d016dc 635
636 SetTRDPidResponseMaster();
5eb3e944 637 //has to precede InitializeTRDResponse(), otherwise the read-out fTRDdEdxParams is not pased in TRDResponse!
9c499471 638 SetTRDdEdxParams();
53d016dc 639 InitializeTRDResponse();
b2138b40 640
641 SetEMCALPidResponseMaster();
642 InitializeEMCALResponse();
4ec8e76d 643
b79db598 644 SetTOFPidResponseMaster();
645 InitializeTOFResponse();
644666df 646
567624b5 647 SetHMPIDPidResponseMaster();
648 InitializeHMPIDResponse();
649
644666df 650 if (fCurrentEvent) fTPCResponse.SetMagField(fCurrentEvent->GetMagneticField());
4ec8e76d 651}
652
1c9d11be 653//______________________________________________________________________________
4ec8e76d 654Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
655{
656 //
657 // Get TPC multiplicity in bins of 150
658 //
659
660 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
661 Double_t tpcMulti=0.;
662 if(vertexTPC){
663 Double_t vertexContribTPC=vertexTPC->GetNContributors();
664 tpcMulti=vertexContribTPC/150.;
665 if (tpcMulti>20.) tpcMulti=20.;
666 }
667
668 return tpcMulti;
669}
670
671//______________________________________________________________________________
672void AliPIDResponse::SetRecoInfo()
673{
674 //
675 // Set reconstruction information
676 //
677
678 //reset information
679 fLHCperiod="";
680 fMCperiodTPC="";
681
682 fBeamType="";
683
684 fBeamType="PP";
42fcc729 685 fBeamTypeNum=kPP;
bbce5a64 686
687 Bool_t hasProdInfo=(fCurrentFile.BeginsWith("LHC"));
4ec8e76d 688
3e95957e 689 TPRegexp reg(".*(LHC1[1-3][a-z]+[0-9]+[a-z_]*)[/_].*");
bbce5a64 690 if (hasProdInfo) reg=TPRegexp("LHC1[1-2][a-z]+[0-9]+[a-z_]*");
ea6ba565 691 TPRegexp reg12a17("LHC1[2-4][a-z]");
1436d6bb 692
4ec8e76d 693 //find the period by run number (UGLY, but not stored in ESD and AOD... )
694 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
695 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
696 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
99e9d5ec 697 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 698 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
699 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
12d3abbc 700 else if (fRun>=136851&&fRun<=139846) {
ea235c90 701 fLHCperiod="LHC10H";
702 fMCperiodTPC="LHC10H8";
703 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
ef7661fd 704 // exception for 13d2 and later
705 if (fCurrentAliRootRev >= 62714) fMCperiodTPC="LHC13D2";
ea235c90 706 fBeamType="PBPB";
42fcc729 707 fBeamTypeNum=kPBPB;
ea235c90 708 }
12d3abbc 709 else if (fRun>=139847&&fRun<=146974) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
710 //TODO: periods 11B (146975-150721), 11C (150722-155837) are not yet treated assume 11d for the moment
711 else if (fRun>=146975&&fRun<=155837) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
712 else if (fRun>=155838&&fRun<=159649) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
713 // also for 11e (159650-162750),f(162751-165771) use 11d
714 else if (fRun>=159650&&fRun<=162750) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
715 else if (fRun>=162751&&fRun<=165771) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
00a38d07 716
12d3abbc 717 else if (fRun>=165772 && fRun<=170718) {
3077a03d 718 fLHCperiod="LHC11H";
719 fMCperiodTPC="LHC11A10";
720 fBeamType="PBPB";
42fcc729 721 fBeamTypeNum=kPBPB;
a78fd045 722 if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
3077a03d 723 }
03696add 724 if (fRun>=170719 && fRun<=177311) {
725 fLHCperiod="LHC12A";
726 fBeamType="PP";
727 fBeamTypeNum=kPP;
728 fMCperiodTPC="LHC10F6A";
729 if (fCurrentAliRootRev >= 62714)
dbeb0132 730 fMCperiodTPC="LHC14E2";
03696add 731 }
3b3bf053 732 // for the moment use LHC12b parameters up to LHC12d
03696add 733 if (fRun>=177312 /*&& fRun<=179356*/) {
734 fLHCperiod="LHC12B";
735 fBeamType="PP";
736 fBeamTypeNum=kPP;
737 fMCperiodTPC="LHC10F6A";
738 if (fCurrentAliRootRev >= 62714)
dbeb0132 739 fMCperiodTPC="LHC14E2";
03696add 740 }
42fcc729 741// if (fRun>=179357 && fRun<=183173) { fLHCperiod="LHC12C"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
742// if (fRun>=183174 && fRun<=186345) { fLHCperiod="LHC12D"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
743// if (fRun>=186346 && fRun<=186635) { fLHCperiod="LHC12E"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
744
745// if (fRun>=186636 && fRun<=188166) { fLHCperiod="LHC12F"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
746// if (fRun >= 188167 && fRun <= 188355 ) { fLHCperiod="LHC12G"; fBeamType="PP"; fBeamTypeNum=kPP;/*fMCperiodTPC="";*/ }
747// if (fRun >= 188356 && fRun <= 188503 ) { fLHCperiod="LHC12G"; fBeamType="PPB"; fBeamTypeNum=kPPB;/*fMCperiodTPC="";*/ }
3b3bf053 748// for the moment use 12g parametrisation for all full gain runs (LHC12e+)
5b7afd75 749 if (fRun >= 186346 && fRun < 188719) { fLHCperiod="LHC12G"; fBeamType="PPB";fBeamTypeNum=kPPB; fMCperiodTPC="LHC12G"; }
87da0205 750
5b7afd75 751 // Dedicated splines for periods 12g and 12i(j) (and use more appropriate MC)
752 if (fRun >= 188720 && fRun <= 192738) {
753 fLHCperiod="LHC12H";
754 fBeamType="PP";
755 fBeamTypeNum=kPP;
756 fMCperiodTPC="LHC10F6A";
757 if (fCurrentAliRootRev >= 62714)
758 fMCperiodTPC="LHC13B2_FIXn1";
759 }
760 if (fRun >= 192739 && fRun <= 194479) {
761 fLHCperiod="LHC12I";
762 fBeamType="PP";
763 fBeamTypeNum=kPP;
764 fMCperiodTPC="LHC10F6A";
765 if (fCurrentAliRootRev >= 62714)
766 fMCperiodTPC="LHC13B2_FIXn1";
767 }
768
87da0205 769 // New parametrisation for 2013 pPb runs
5b7afd75 770 if (fRun >= 194480) {
87da0205 771 fLHCperiod="LHC13B";
772 fBeamType="PPB";
42fcc729 773 fBeamTypeNum=kPPB;
bbce5a64 774 fMCperiodTPC="LHC12G";
87da0205 775
776 if (fCurrentAliRootRev >= 61605)
777 fMCperiodTPC="LHC13B2_FIX";
bbce5a64 778 if (fCurrentAliRootRev >= 62714)
779 fMCperiodTPC="LHC13B2_FIXn1";
c4bec231 780
781 // High luminosity pPb runs require different parametrisations
782 if (fRun >= 195875 && fRun <= 197411) {
783 fLHCperiod="LHC13F";
784 }
87da0205 785 }
80ab5635 786
66de625c 787 //exception new pp MC productions from 2011 (11a periods have 10f6a splines!)
42fcc729 788 if (fBeamType=="PP" && reg.MatchB(fCurrentFile) && !fCurrentFile.Contains("LHC11a")) { fMCperiodTPC="LHC11B2"; fBeamType="PP";fBeamTypeNum=kPP; }
4a527e08 789 // exception for 11f1
bbce5a64 790 if (fCurrentFile.Contains("LHC11f1")) fMCperiodTPC="LHC11F1";
bf26ce58 791 // exception for 12f1a, 12f1b and 12i3
bbce5a64 792 if (fCurrentFile.Contains("LHC12f1") || fCurrentFile.Contains("LHC12i3")) fMCperiodTPC="LHC12F1";
c3ee524d 793 // exception for 12c4
bbce5a64 794 if (fCurrentFile.Contains("LHC12c4")) fMCperiodTPC="LHC12C4";
361813a2 795 // exception for 13d1 11d anchored prod
796 if (fLHCperiod=="LHC11D" && fCurrentFile.Contains("LHC13d1")) fMCperiodTPC="LHC13D1";
4ec8e76d 797}
798
799//______________________________________________________________________________
800void AliPIDResponse::SetITSParametrisation()
801{
802 //
803 // Set the ITS parametrisation
804 //
805}
806
f84b18dd 807
808//______________________________________________________________________________
809void AliPIDResponse::AddPointToHyperplane(TH2D* h, TLinearFitter* linExtrapolation, Int_t binX, Int_t binY)
810{
811 if (h->GetBinContent(binX, binY) <= 1e-4)
812 return; // Reject bins without content (within some numerical precision) or with strange content
813
814 Double_t coord[2] = {0, 0};
815 coord[0] = h->GetXaxis()->GetBinCenter(binX);
816 coord[1] = h->GetYaxis()->GetBinCenter(binY);
817 Double_t binError = h->GetBinError(binX, binY);
818 if (binError <= 0) {
819 binError = 1000; // Should not happen because bins without content are rejected for the map (TH2D* h)
820 printf("ERROR: This should never happen: Trying to add bin in addPointToHyperplane with error not set....\n");
821 }
822 linExtrapolation->AddPoint(coord, h->GetBinContent(binX, binY, binError));
823}
824
825
826//______________________________________________________________________________
827TH2D* AliPIDResponse::RefineHistoViaLinearInterpolation(TH2D* h, Double_t refineFactorX, Double_t refineFactorY)
828{
829 if (!h)
830 return 0x0;
831
832 // Interpolate to finer map
833 TLinearFitter* linExtrapolation = new TLinearFitter(2, "hyp2", "");
834
835 Double_t upperMapBoundY = h->GetYaxis()->GetBinUpEdge(h->GetYaxis()->GetNbins());
836 Double_t lowerMapBoundY = h->GetYaxis()->GetBinLowEdge(1);
1b45e564 837 Int_t nBinsX = 30;
f84b18dd 838 // Binning was find to yield good results, if 40 bins are chosen for the range 0.0016 to 0.02. For the new variable range,
839 // scale the number of bins correspondingly
1b45e564 840 Int_t nBinsY = TMath::Nint((upperMapBoundY - lowerMapBoundY) / (0.02 - 0.0016) * 40);
f84b18dd 841 Int_t nBinsXrefined = nBinsX * refineFactorX;
842 Int_t nBinsYrefined = nBinsY * refineFactorY;
843
844 TH2D* hRefined = new TH2D(Form("%s_refined", h->GetName()), Form("%s (refined)", h->GetTitle()),
845 nBinsXrefined, h->GetXaxis()->GetBinLowEdge(1), h->GetXaxis()->GetBinUpEdge(h->GetXaxis()->GetNbins()),
846 nBinsYrefined, lowerMapBoundY, upperMapBoundY);
847
848 for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
849 for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
850
851 hRefined->SetBinContent(binX, binY, 1); // Default value is 1
852
853 Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
854 Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
855
1b45e564 856 /*OLD
f84b18dd 857 linExtrapolation->ClearPoints();
858
859 // For interpolation: Just take the corresponding bin from the old histo.
860 // For extrapolation: take the last available bin from the old histo.
861 // If the boundaries are to be skipped, also skip the corresponding bins
862 Int_t oldBinX = h->GetXaxis()->FindBin(centerX);
863 if (oldBinX < 1)
864 oldBinX = 1;
865 if (oldBinX > nBinsX)
866 oldBinX = nBinsX;
867
868 Int_t oldBinY = h->GetYaxis()->FindBin(centerY);
869 if (oldBinY < 1)
870 oldBinY = 1;
871 if (oldBinY > nBinsY)
872 oldBinY = nBinsY;
873
874 // Neighbours left column
875 if (oldBinX >= 2) {
876 if (oldBinY >= 2) {
877 AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY - 1);
878 }
879
880 AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY);
881
882 if (oldBinY < nBinsY) {
883 AddPointToHyperplane(h, linExtrapolation, oldBinX - 1, oldBinY + 1);
884 }
885 }
886
887 // Neighbours (and point itself) same column
888 if (oldBinY >= 2) {
889 AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY - 1);
890 }
891
892 AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY);
893
894 if (oldBinY < nBinsY) {
895 AddPointToHyperplane(h, linExtrapolation, oldBinX, oldBinY + 1);
896 }
897
898 // Neighbours right column
899 if (oldBinX < nBinsX) {
900 if (oldBinY >= 2) {
901 AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY - 1);
902 }
903
904 AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY);
905
906 if (oldBinY < nBinsY) {
907 AddPointToHyperplane(h, linExtrapolation, oldBinX + 1, oldBinY + 1);
908 }
909 }
910
911
912 // Fit 2D-hyperplane
913 if (linExtrapolation->GetNpoints() <= 0)
914 continue;
915
916 if (linExtrapolation->Eval() != 0)// EvalRobust -> Takes much, much, [...], much more time (~hours instead of seconds)
917 continue;
918
919 // Fill the bin of the refined histogram with the extrapolated value
920 Double_t interpolatedValue = linExtrapolation->GetParameter(0) + linExtrapolation->GetParameter(1) * centerX
921 + linExtrapolation->GetParameter(2) * centerY;
922 */
f85a3764 923 Double_t interpolatedValue = h->Interpolate(centerX, centerY) ;
f84b18dd 924 hRefined->SetBinContent(binX, binY, interpolatedValue);
925 }
926 }
927
1b45e564 928
929 // Problem: Interpolation does not work before/beyond center of first/last bin (as the name suggests).
930 // Therefore, for each row in dEdx: Take last bin from old map and interpolate values from center and edge.
931 // Assume line through these points and extropolate to last bin of refined map
932 const Double_t firstOldXbinUpEdge = h->GetXaxis()->GetBinUpEdge(1);
933 const Double_t firstOldXbinCenter = h->GetXaxis()->GetBinCenter(1);
934
935 const Double_t oldXbinHalfWidth = firstOldXbinUpEdge - firstOldXbinCenter;
936
937 const Double_t lastOldXbinLowEdge = h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
938 const Double_t lastOldXbinCenter = h->GetXaxis()->GetBinCenter(h->GetNbinsX());
939
940 for (Int_t binY = 1; binY <= nBinsYrefined; binY++) {
941 Double_t centerY = hRefined->GetYaxis()->GetBinCenter(binY);
942
943 const Double_t interpolatedCenterFirstXbin = h->Interpolate(firstOldXbinCenter, centerY);
944 const Double_t interpolatedUpEdgeFirstXbin = h->Interpolate(firstOldXbinUpEdge, centerY);
945
946 const Double_t extrapolationSlopeFirstXbin = (interpolatedUpEdgeFirstXbin - interpolatedCenterFirstXbin) / oldXbinHalfWidth;
947 const Double_t extrapolationOffsetFirstXbin = interpolatedCenterFirstXbin;
948
949
950 const Double_t interpolatedCenterLastXbin = h->Interpolate(lastOldXbinCenter, centerY);
951 const Double_t interpolatedLowEdgeLastXbin = h->Interpolate(lastOldXbinLowEdge, centerY);
952
953 const Double_t extrapolationSlopeLastXbin = (interpolatedCenterLastXbin - interpolatedLowEdgeLastXbin) / oldXbinHalfWidth;
954 const Double_t extrapolationOffsetLastXbin = interpolatedCenterLastXbin;
955
956 for (Int_t binX = 1; binX <= nBinsXrefined; binX++) {
957 Double_t centerX = hRefined->GetXaxis()->GetBinCenter(binX);
958
959 if (centerX < firstOldXbinCenter) {
960 Double_t extrapolatedValue = extrapolationOffsetFirstXbin + (centerX - firstOldXbinCenter) * extrapolationSlopeFirstXbin;
961 hRefined->SetBinContent(binX, binY, extrapolatedValue);
962 }
963 else if (centerX <= lastOldXbinCenter) {
964 continue;
965 }
966 else {
967 Double_t extrapolatedValue = extrapolationOffsetLastXbin + (centerX - lastOldXbinCenter) * extrapolationSlopeLastXbin;
968 hRefined->SetBinContent(binX, binY, extrapolatedValue);
969 }
970 }
971 }
972
f84b18dd 973 delete linExtrapolation;
974
975 return hRefined;
976}
977
978//______________________________________________________________________________
979void AliPIDResponse::SetTPCEtaMaps(Double_t refineFactorMapX, Double_t refineFactorMapY,
980 Double_t refineFactorSigmaMapX, Double_t refineFactorSigmaMapY)
981{
982 //
983 // Load the TPC eta correction maps from the OADB
984 //
985
f85a3764 986 if (fUseTPCEtaCorrection == kFALSE) {
987 // Disable eta correction via setting no maps
988 if (!fTPCResponse.SetEtaCorrMap(0x0))
1b45e564 989 AliInfo("Request to disable TPC eta correction -> Eta correction has been disabled");
f85a3764 990 else
991 AliError("Request to disable TPC eta correction -> Some error occured when unloading the correction maps");
992
993 if (!fTPCResponse.SetSigmaParams(0x0, 0))
1b45e564 994 AliInfo("Request to disable TPC eta correction -> Using old parametrisation for sigma");
995 else
f85a3764 996 AliError("Request to disable TPC eta correction -> Some error occured when unloading the sigma maps");
997
998 return;
999 }
1b45e564 1000
f84b18dd 1001 TString dataType = "DATA";
1002 TString period = fLHCperiod.IsNull() ? "No period information" : fLHCperiod;
1003
1004 if (fIsMC) {
87da0205 1005 if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
f85a3764 1006 period=fMCperiodTPC;
1007 dataType="MC";
1008 }
f84b18dd 1009 fRecoPass = 1;
1010
87da0205 1011 if (!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) && fMCperiodTPC.IsNull()) {
f84b18dd 1012 AliFatal("MC detected, but no MC period set -> Not changing eta maps!");
1013 return;
1014 }
f84b18dd 1015 }
f85a3764 1016
1017 Int_t recopass = fRecoPass;
a2c30af1 1018 if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC) )
f85a3764 1019 recopass = fRecoPassUser;
f84b18dd 1020
f85a3764 1021 TString defaultObj = Form("Default_%s_pass%d", dataType.Data(), recopass);
f84b18dd 1022
f85a3764 1023 AliInfo(Form("Current period and reco pass: %s.pass%d", period.Data(), recopass));
f84b18dd 1024
1025 // Invalidate old maps
1026 fTPCResponse.SetEtaCorrMap(0x0);
1027 fTPCResponse.SetSigmaParams(0x0, 0);
1028
d5113358 1029
1030 TString fileNameMaps(Form("%s/COMMON/PID/data/TPCetaMaps.root", fOADBPath.Data()));
1031 if (!fCustomTPCetaMaps.IsNull()) fileNameMaps=fCustomTPCetaMaps;
1032
f84b18dd 1033 // Load the eta correction maps
f85a3764 1034 AliOADBContainer etaMapsCont(Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
f84b18dd 1035
d5113358 1036 Int_t statusCont = etaMapsCont.InitFromFile(fileNameMaps.Data(), Form("TPCetaMaps_%s_pass%d", dataType.Data(), recopass));
f84b18dd 1037 if (statusCont) {
1038 AliError("Failed initializing TPC eta correction maps from OADB -> Disabled eta correction");
87da0205 1039 fUseTPCEtaCorrection = kFALSE;
f84b18dd 1040 }
1041 else {
d5113358 1042 AliInfo(Form("Loading TPC eta correction map from %s", fileNameMaps.Data()));
f84b18dd 1043
1044 TH2D* etaMap = 0x0;
1045
87da0205 1046 if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
f85a3764 1047 TString searchMap = Form("TPCetaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
f84b18dd 1048 etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(searchMap.Data()));
1049 if (!etaMap) {
1050 // Try default object
1051 etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetDefaultObject(defaultObj.Data()));
1052 }
1053 }
1054 else {
1055 etaMap = dynamic_cast<TH2D *>(etaMapsCont.GetObject(fRun, defaultObj.Data()));
1056 }
1057
1058
1059 if (!etaMap) {
1060 AliError(Form("TPC eta correction map not found for run %d and also no default map found -> Disabled eta correction!!!", fRun));
87da0205 1061 fUseTPCEtaCorrection = kFALSE;
f84b18dd 1062 }
1063 else {
1064 TH2D* etaMapRefined = RefineHistoViaLinearInterpolation(etaMap, refineFactorMapX, refineFactorMapY);
1065
1066 if (etaMapRefined) {
1067 if (!fTPCResponse.SetEtaCorrMap(etaMapRefined)) {
1068 AliError(Form("Failed to set TPC eta correction map for run %d -> Disabled eta correction!!!", fRun));
1069 fTPCResponse.SetEtaCorrMap(0x0);
87da0205 1070 fUseTPCEtaCorrection = kFALSE;
f84b18dd 1071 }
1072 else {
d5113358 1073 AliInfo(Form("Loaded TPC eta correction map (refine factors %.2f/%.2f) from %s: %s (MD5(map) = %s)",
1074 refineFactorMapX, refineFactorMapY, fileNameMaps.Data(), fTPCResponse.GetEtaCorrMap()->GetTitle(),
5a9dc560 1075 GetChecksum(fTPCResponse.GetEtaCorrMap()).Data()));
f84b18dd 1076 }
1077
1078 delete etaMapRefined;
1079 }
1080 else {
1081 AliError(Form("Failed to set TPC eta correction map for run %d (map was loaded, but couldn't be refined) -> Disabled eta correction!!!", fRun));
87da0205 1082 fUseTPCEtaCorrection = kFALSE;
f84b18dd 1083 }
1084 }
1085 }
1086
87da0205 1087 // If there was some problem loading the eta maps, it makes no sense to load the sigma maps (that require eta corrected data)
1088 if (fUseTPCEtaCorrection == kFALSE) {
1089 AliError("Failed to load TPC eta correction map required by sigma maps -> Using old parametrisation for sigma");
1090 return;
1091 }
1092
f84b18dd 1093 // Load the sigma parametrisation (1/dEdx vs tanTheta_local (~eta))
f85a3764 1094 AliOADBContainer etaSigmaMapsCont(Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
f84b18dd 1095
d5113358 1096 statusCont = etaSigmaMapsCont.InitFromFile(fileNameMaps.Data(), Form("TPCetaSigmaMaps_%s_pass%d", dataType.Data(), recopass));
f84b18dd 1097 if (statusCont) {
1098 AliError("Failed initializing TPC eta sigma maps from OADB -> Using old sigma parametrisation");
1099 }
1100 else {
d5113358 1101 AliInfo(Form("Loading TPC eta sigma map from %s", fileNameMaps.Data()));
f84b18dd 1102
1103 TObjArray* etaSigmaPars = 0x0;
1104
87da0205 1105 if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) {
f85a3764 1106 TString searchMap = Form("TPCetaSigmaMaps_%s_%s_pass%d", dataType.Data(), period.Data(), recopass);
f84b18dd 1107 etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(searchMap.Data()));
1108 if (!etaSigmaPars) {
1109 // Try default object
1110 etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetDefaultObject(defaultObj.Data()));
1111 }
1112 }
1113 else {
1114 etaSigmaPars = dynamic_cast<TObjArray *>(etaSigmaMapsCont.GetObject(fRun, defaultObj.Data()));
1115 }
1116
1117 if (!etaSigmaPars) {
1118 AliError(Form("TPC eta sigma parametrisation not found for run %d -> Using old sigma parametrisation!!!", fRun));
1119 }
1120 else {
1121 TH2D* etaSigmaPar1Map = dynamic_cast<TH2D *>(etaSigmaPars->FindObject("sigmaPar1Map"));
1122 TNamed* sigmaPar0Info = dynamic_cast<TNamed *>(etaSigmaPars->FindObject("sigmaPar0"));
1123 Double_t sigmaPar0 = 0.0;
1124
1125 if (sigmaPar0Info) {
1126 TString sigmaPar0String = sigmaPar0Info->GetTitle();
1127 sigmaPar0 = sigmaPar0String.Atof();
1128 }
1129 else {
1130 // Something is weired because the object for parameter 0 could not be loaded -> New sigma parametrisation can not be used!
1131 etaSigmaPar1Map = 0x0;
1132 }
1133
1134 TH2D* etaSigmaPar1MapRefined = RefineHistoViaLinearInterpolation(etaSigmaPar1Map, refineFactorSigmaMapX, refineFactorSigmaMapY);
1135
1136
1137 if (etaSigmaPar1MapRefined) {
1138 if (!fTPCResponse.SetSigmaParams(etaSigmaPar1MapRefined, sigmaPar0)) {
1139 AliError(Form("Failed to set TPC eta sigma map for run %d -> Using old sigma parametrisation!!!", fRun));
1140 fTPCResponse.SetSigmaParams(0x0, 0);
1141 }
1142 else {
d5113358 1143 AliInfo(Form("Loaded TPC sigma correction map (refine factors %.2f/%.2f) from %s: %s (MD5(map) = %s, sigmaPar0 = %f)",
1144 refineFactorSigmaMapX, refineFactorSigmaMapY, fileNameMaps.Data(), fTPCResponse.GetSigmaPar1Map()->GetTitle(),
5a9dc560 1145 GetChecksum(fTPCResponse.GetSigmaPar1Map()).Data(), sigmaPar0));
f84b18dd 1146 }
1147
1148 delete etaSigmaPar1MapRefined;
1149 }
1150 else {
1151 AliError(Form("Failed to set TPC eta sigma map for run %d (map was loaded, but couldn't be refined) -> Using old sigma parametrisation!!!",
1152 fRun));
1153 }
1154 }
1155 }
1156}
1157
4ec8e76d 1158//______________________________________________________________________________
1159void AliPIDResponse::SetTPCPidResponseMaster()
1160{
1161 //
1162 // Load the TPC pid response functions from the OADB
644666df 1163 // Load the TPC voltage maps from OADB
4ec8e76d 1164 //
09b50a42 1165 //don't load twice for the moment
1166 if (fArrPidResponseMaster) return;
1167
1168
4ec8e76d 1169 //reset the PID response functions
1170 delete fArrPidResponseMaster;
644666df 1171 fArrPidResponseMaster=NULL;
4ec8e76d 1172
644666df 1173 TFile *f=NULL;
4ec8e76d 1174
644666df 1175 TString fileNamePIDresponse(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
d5113358 1176 if (!fCustomTPCpidResponse.IsNull()) fileNamePIDresponse=fCustomTPCpidResponse;
1177
644666df 1178 f=TFile::Open(fileNamePIDresponse.Data());
ea235c90 1179 if (f && f->IsOpen() && !f->IsZombie()){
1180 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 1181 }
ea235c90 1182 delete f;
644666df 1183
1184 TString fileNameVoltageMaps(Form("%s/COMMON/PID/data/TPCvoltageSettings.root", fOADBPath.Data()));
1185 f=TFile::Open(fileNameVoltageMaps.Data());
1186 if (f && f->IsOpen() && !f->IsZombie()){
1187 fOADBvoltageMaps=dynamic_cast<AliOADBContainer*>(f->Get("TPCvoltageSettings"));
1188 }
1189 delete f;
4ec8e76d 1190
1191 if (!fArrPidResponseMaster){
644666df 1192 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileNamePIDresponse.Data()));
4ec8e76d 1193 return;
1194 }
1195 fArrPidResponseMaster->SetOwner();
644666df 1196
1197 if (!fOADBvoltageMaps)
1198 {
1199 AliFatal(Form("Could not retrieve the TPC voltage maps from: %s",fileNameVoltageMaps.Data()));
1200 }
1201 fArrPidResponseMaster->SetOwner();
4ec8e76d 1202}
1203
1204//______________________________________________________________________________
1205void AliPIDResponse::SetTPCParametrisation()
1206{
1207 //
1208 // Change BB parametrisation for current run
1209 //
1210
12d3abbc 1211 //
1212 //reset old splines
1213 //
1214 fTPCResponse.ResetSplines();
1215
4ec8e76d 1216 if (fLHCperiod.IsNull()) {
12d3abbc 1217 AliError("No period set, not changing parametrisation");
4ec8e76d 1218 return;
1219 }
1220
1221 //
1222 // Set default parametrisations for data and MC
1223 //
1224
1225 //data type
1226 TString datatype="DATA";
1227 //in case of mc fRecoPass is per default 1
1228 if (fIsMC) {
87da0205 1229 if(!(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) datatype="MC";
539a5a59 1230 fRecoPass=1;
4ec8e76d 1231 }
f84b18dd 1232
4a527e08 1233 // period
1234 TString period=fLHCperiod;
87da0205 1235 if (fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))) period=fMCperiodTPC;
4a527e08 1236
f85a3764 1237 Int_t recopass = fRecoPass;
87da0205 1238 if(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) recopass = fRecoPassUser;
f85a3764 1239
1240 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
4a527e08 1241 Bool_t found=kFALSE;
4ec8e76d 1242 //
1243 //set the new PID splines
1244 //
4ec8e76d 1245 if (fArrPidResponseMaster){
4ec8e76d 1246 //for MC don't use period information
644666df 1247 //if (fIsMC) period="[A-Z0-9]*";
4ec8e76d 1248 //for MC use MC period information
644666df 1249 //pattern for the default entry (valid for all particles)
de678885 1250 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN(_*)([A-Z1-9]*)",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
644666df 1251
f85a3764 1252 //find particle id and gain scenario
644666df 1253 for (Int_t igainScenario=0; igainScenario<AliTPCPIDResponse::fgkNumberOfGainScenarios; igainScenario++)
1254 {
1255 TObject *grAll=NULL;
1256 TString gainScenario = AliTPCPIDResponse::GainScenarioName(igainScenario);
1257 gainScenario.ToUpper();
1258 //loop over entries and filter them
1259 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp)
1260 {
1261 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
1262 if (responseFunction==NULL) continue;
1263 TString responseName=responseFunction->GetName();
1264
1265 if (!reg.MatchB(responseName)) continue;
1266
1267 TObjArray *arr=reg.MatchS(responseName); if (!arr) continue;
1268 TObject* tmp=NULL;
1269 tmp=arr->At(1); if (!tmp) continue;
1270 TString particleName=tmp->GetName();
1271 tmp=arr->At(3); if (!tmp) continue;
1272 TString gainScenarioName=tmp->GetName();
1273 delete arr;
1274 if (particleName.IsNull()) continue;
1275 if (!grAll && particleName=="ALL" && gainScenarioName==gainScenario) grAll=responseFunction;
1276 else
1277 {
1278 for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
1279 {
1280 TString particle=AliPID::ParticleName(ispec);
1281 particle.ToUpper();
1282 //std::cout<<responseName<<" "<<particle<<" "<<particleName<<" "<<gainScenario<<" "<<gainScenarioName<<std::endl;
1283 if ( particle == particleName && gainScenario == gainScenarioName )
1284 {
1285 fTPCResponse.SetResponseFunction( responseFunction,
1286 (AliPID::EParticleType)ispec,
1287 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1288 fTPCResponse.SetUseDatabase(kTRUE);
5a9dc560 1289 AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunction->GetName(),
1290 GetChecksum((TSpline3*)responseFunction).Data()));
644666df 1291 found=kTRUE;
644666df 1292 break;
1293 }
4ec8e76d 1294 }
1295 }
1296 }
bf26ce58 1297
1298 // Retrieve responsefunction for pions - will (if available) be used for muons if there are no dedicated muon splines.
1299 // For light nuclei, try to set the proton spline, if no dedicated splines are available.
1300 // In both cases: Use default splines, if no dedicated splines and no pion/proton splines are available.
1301 TObject* responseFunctionPion = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kPion,
1302 (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1303 TObject* responseFunctionProton = fTPCResponse.GetResponseFunction( (AliPID::EParticleType)AliPID::kProton,
1304 (AliTPCPIDResponse::ETPCgainScenario)igainScenario);
1305
1306 for (Int_t ispec=0; ispec<(AliTPCPIDResponse::fgkNumberOfParticleSpecies); ++ispec)
644666df 1307 {
bf26ce58 1308 if (!fTPCResponse.GetResponseFunction( (AliPID::EParticleType)ispec,
1309 (AliTPCPIDResponse::ETPCgainScenario)igainScenario))
644666df 1310 {
bf26ce58 1311 if (ispec == AliPID::kMuon) { // Muons
1312 if (responseFunctionPion) {
1313 fTPCResponse.SetResponseFunction( responseFunctionPion,
1314 (AliPID::EParticleType)ispec,
1315 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1316 fTPCResponse.SetUseDatabase(kTRUE);
5a9dc560 1317 AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionPion->GetName(),
1318 GetChecksum((TSpline3*)responseFunctionPion).Data()));
bf26ce58 1319 found=kTRUE;
1320 }
1321 else if (grAll) {
1322 fTPCResponse.SetResponseFunction( grAll,
1323 (AliPID::EParticleType)ispec,
1324 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1325 fTPCResponse.SetUseDatabase(kTRUE);
5a9dc560 1326 AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
1327 GetChecksum((TSpline3*)grAll).Data()));
bf26ce58 1328 found=kTRUE;
1329 }
1330 //else
1331 // AliError(Form("No splines found for muons (also no pion splines and no default splines) for gain scenario %d!", igainScenario));
1332 }
1333 else if (ispec >= AliPID::kSPECIES) { // Light nuclei
1334 if (responseFunctionProton) {
1335 fTPCResponse.SetResponseFunction( responseFunctionProton,
1336 (AliPID::EParticleType)ispec,
1337 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1338 fTPCResponse.SetUseDatabase(kTRUE);
5a9dc560 1339 AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,responseFunctionProton->GetName(),
1340 GetChecksum((TSpline3*)responseFunctionProton).Data()));
bf26ce58 1341 found=kTRUE;
1342 }
1343 else if (grAll) {
644666df 1344 fTPCResponse.SetResponseFunction( grAll,
1345 (AliPID::EParticleType)ispec,
1346 (AliTPCPIDResponse::ETPCgainScenario)igainScenario );
1347 fTPCResponse.SetUseDatabase(kTRUE);
5a9dc560 1348 AliInfo(Form("Adding graph: %d %d - %s (MD5(spline) = %s)",ispec,igainScenario,grAll->GetName(),
1349 GetChecksum((TSpline3*)grAll).Data()));
644666df 1350 found=kTRUE;
bf26ce58 1351 }
1352 //else
1353 // AliError(Form("No splines found for species %d (also no proton splines and no default splines) for gain scenario %d!",
1354 // ispec, igainScenario));
644666df 1355 }
4ec8e76d 1356 }
1357 }
1358 }
1359 }
644666df 1360 else AliInfo("no fArrPidResponseMaster");
4a527e08 1361
1362 if (!found){
f85a3764 1363 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
4a527e08 1364 }
644666df 1365
87da0205 1366
1367 //
22158469 1368 // Setup multiplicity correction (only used for non-pp collisions)
87da0205 1369 //
22158469
JW		 1370
1371 const Bool_t isPP = (fBeamType.CompareTo("PP") == 0);
1372
1373 // 2013 pPb data taking at low luminosity
1374 const Bool_t isPPb2013LowLuminosity = period.Contains("LHC13B") || period.Contains("LHC13C") || period.Contains("LHC13D");
1375 // PbPb 2010, period 10h.pass2
1376 //TODO Needs further development const Bool_t is10hpass2 = period.Contains("LHC10H") && recopass == 2;
1377
eafa51b5
JW		 1378
1379 // In case of MC without(!) tune on data activated for the TPC, don't use the multiplicity correction for the moment
1380 Bool_t isMCandNotTPCtuneOnData = fIsMC && !(fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC));
1381
22158469 1382 // If correction is available, but disabled (highly NOT recommended!), print warning
eafa51b5 1383 if (!fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
22158469
JW		 1384 //TODO: Needs further development if (is10hpass2 || isPPb2013LowLuminosity) {
1385 if (isPPb2013LowLuminosity) {
1386 AliWarning("Mulitplicity correction disabled, but correction parameters for this period exist. It is highly recommended to use enable the correction. Otherwise the splines might be off!");
1387 }
1388 }
1389
eafa51b5 1390 if (fUseTPCMultiplicityCorrection && !isPP && !isMCandNotTPCtuneOnData) {
87da0205 1391 AliInfo("Multiplicity correction enabled!");
1392
1393 //TODO After testing, load parameters from outside
22158469 1394 /*TODO no correction for MC
87da0205 1395 if (period.Contains("LHC11A10")) {//LHC11A10A
1396 AliInfo("Using multiplicity correction parameters for 11a10!");
1397 fTPCResponse.SetParameterMultiplicityCorrection(0, 6.90133e-06);
1398 fTPCResponse.SetParameterMultiplicityCorrection(1, -1.22123e-03);
1399 fTPCResponse.SetParameterMultiplicityCorrection(2, 1.80220e-02);
1400 fTPCResponse.SetParameterMultiplicityCorrection(3, 0.1);
1401 fTPCResponse.SetParameterMultiplicityCorrection(4, 6.45306e-03);
1402
1403 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -2.85505e-07);
1404 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, -1.31911e-06);
1405 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1406
1407 fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -4.29665e-05);
1408 fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 1.37023e-02);
1409 fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -6.36337e-01);
1410 fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.13479e-02);
1411 }
22158469
JW		 1412 else*/ if (isPPb2013LowLuminosity) {// 2013 pPb data taking at low luminosity
1413 AliInfo("Using multiplicity correction parameters for 13b.pass2 (at least also valid for 13{c,d} and pass 3)!");
ef7661fd 1414
1415 fTPCResponse.SetParameterMultiplicityCorrection(0, -5.906e-06);
1416 fTPCResponse.SetParameterMultiplicityCorrection(1, -5.064e-04);
1417 fTPCResponse.SetParameterMultiplicityCorrection(2, -3.521e-02);
1418 fTPCResponse.SetParameterMultiplicityCorrection(3, 2.469e-02);
1419 fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1420
1421 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.32e-06);
1422 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.177e-05);
1423 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1424
1425 fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 0.);
1426 fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 0.);
1427 fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 0.);
1428 fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 0.);
1429
1430 /* Not too bad, but far from perfect in the details
1431 fTPCResponse.SetParameterMultiplicityCorrection(0, -6.27187e-06);
1432 fTPCResponse.SetParameterMultiplicityCorrection(1, -4.60649e-04);
1433 fTPCResponse.SetParameterMultiplicityCorrection(2, -4.26450e-02);
1434 fTPCResponse.SetParameterMultiplicityCorrection(3, 2.40590e-02);
1435 fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1436
1437 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, -5.338e-06);
1438 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 1.220e-05);
1439 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1440
1441 fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, 7.89237e-05);
1442 fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, -1.30662e-02);
1443 fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, 8.91548e-01);
1444 fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.47931e-02);
1445 */
1446 }
22158469
JW		 1447 /*TODO: Needs further development
1448 else if (is10hpass2) {
87da0205 1449 AliInfo("Using multiplicity correction parameters for 10h.pass2!");
1450 fTPCResponse.SetParameterMultiplicityCorrection(0, 3.21636e-07);
1451 fTPCResponse.SetParameterMultiplicityCorrection(1, -6.65876e-04);
1452 fTPCResponse.SetParameterMultiplicityCorrection(2, 1.28786e-03);
1453 fTPCResponse.SetParameterMultiplicityCorrection(3, 1.47677e-02);
1454 fTPCResponse.SetParameterMultiplicityCorrection(4, 0);
1455
1456 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(0, 7.23591e-08);
1457 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(1, 2.7469e-06);
1458 fTPCResponse.SetParameterMultiplicityCorrectionTanTheta(2, -0.5);
1459
1460 fTPCResponse.SetParameterMultiplicitySigmaCorrection(0, -1.22590e-05);
1461 fTPCResponse.SetParameterMultiplicitySigmaCorrection(1, 6.88888e-03);
1462 fTPCResponse.SetParameterMultiplicitySigmaCorrection(2, -3.20788e-01);
1463 fTPCResponse.SetParameterMultiplicitySigmaCorrection(3, 1.07345e-02);
1464 }
22158469 1465 */
87da0205 1466 else {
1467 AliError(Form("Multiplicity correction is enabled, but no multiplicity correction parameters have been found for period %s.pass%d -> Mulitplicity correction DISABLED!", period.Data(), recopass));
1468 fUseTPCMultiplicityCorrection = kFALSE;
1469 fTPCResponse.ResetMultiplicityCorrectionFunctions();
1470 }
1471 }
1472 else {
1473 // Just set parameters such that overall correction factor is 1, i.e. no correction.
1474 // This is just a reasonable choice for the parameters for safety reasons. Disabling
1475 // the multiplicity correction will anyhow skip the calculation of the corresponding
1476 // correction factor inside THIS class. Nevertheless, experts can access the TPCPIDResponse
1477 // directly and use it for calculations - which should still give valid results, even if
1478 // the multiplicity correction is explicitely enabled in such expert calls.
1479
eafa51b5
JW		 1480 TString reasonForDisabling = "requested by user";
1481 if (fUseTPCMultiplicityCorrection) {
1482 if (isPP)
1483 reasonForDisabling = "pp collisions";
1484 else
1485 reasonForDisabling = "MC w/o tune on data";
1486 }
1487
87da0205 1488 AliInfo(Form("Multiplicity correction %sdisabled (%s)!", fUseTPCMultiplicityCorrection ? "automatically " : "",
eafa51b5 1489 reasonForDisabling.Data()));
87da0205 1490
1491 fUseTPCMultiplicityCorrection = kFALSE;
1492 fTPCResponse.ResetMultiplicityCorrectionFunctions();
1493 }
1494
22158469
JW		 1495 if (fUseTPCMultiplicityCorrection) {
1496 for (Int_t i = 0; i <= 4 + 1; i++) {
1497 AliInfo(Form("parMultCorr: %d, %e", i, fTPCResponse.GetMultiplicityCorrectionFunction()->GetParameter(i)));
1498 }
1499 for (Int_t j = 0; j <= 2 + 1; j++) {
1500 AliInfo(Form("parMultCorrTanTheta: %d, %e", j, fTPCResponse.GetMultiplicityCorrectionFunctionTanTheta()->GetParameter(j)));
1501 }
1502 for (Int_t j = 0; j <= 3 + 1; j++) {
1503 AliInfo(Form("parMultSigmaCorr: %d, %e", j, fTPCResponse.GetMultiplicitySigmaCorrectionFunction()->GetParameter(j)));
1504 }
87da0205 1505 }
1506
4ec8e76d 1507 //
87da0205 1508 // Setup old resolution parametrisation
4ec8e76d 1509 //
1510
1511 //default
1512 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
1513
3b3bf053 1514 if (fRun>=122195){ //LHC10d
4ec8e76d 1515 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
1516 }
3b3bf053 1517
1518 if (fRun>=170719){ // LHC12a
1519 fTPCResponse.SetSigma(2.95714e-03, 1.01953e+05);
1520 }
1521
1522 if (fRun>=177312){ // LHC12b
1523 fTPCResponse.SetSigma(3.74633e-03, 7.11829e+04 );
1524 }
1525
1526 if (fRun>=186346){ // LHC12e
723c4874 1527 fTPCResponse.SetSigma(8.62022e-04, 9.08156e+05);
1528 }
1529
23425eb2 1530 if (fArrPidResponseMaster)
f85a3764 1531 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
4ec8e76d 1532
5a9dc560 1533 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s (MD5(corr function) = %s)",
1534 fResolutionCorrection->GetName(), GetChecksum(fResolutionCorrection).Data()));
644666df 1535
1536 //read in the voltage map
12d3abbc 1537 TVectorF* gsm = 0x0;
1538 if (fOADBvoltageMaps) gsm=dynamic_cast<TVectorF*>(fOADBvoltageMaps->GetObject(fRun));
644666df 1539 if (gsm)
1540 {
1541 fTPCResponse.SetVoltageMap(*gsm);
1542 TString vals;
1543 AliInfo(Form("Reading the voltage map for run %d\n",fRun));
1544 vals="IROC A: "; for (Int_t i=0; i<18; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1545 AliInfo(vals.Data());
1546 vals="IROC C: "; for (Int_t i=18; i<36; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1547 AliInfo(vals.Data());
1548 vals="OROC A: "; for (Int_t i=36; i<54; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1549 AliInfo(vals.Data());
1550 vals="OROC C: "; for (Int_t i=54; i<72; i++){vals+=Form("%.2f ",(*gsm)[i]);}
1551 AliInfo(vals.Data());
1552 }
1553 else AliInfo("no voltage map, ideal default assumed");
4ec8e76d 1554}
1555
ea235c90 1556//______________________________________________________________________________
1557void AliPIDResponse::SetTRDPidResponseMaster()
1558{
1559 //
1560 // Load the TRD pid params and references from the OADB
1561 //
db0e2c5f 1562 if(fTRDPIDResponseObject) return;
53d016dc 1563 AliOADBContainer contParams("contParams");
1564
db0e2c5f 1565 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
1566 if(statusResponse){
1567 AliError("Failed initializing PID Response Object from OADB");
59a8e853 1568 } else {
db0e2c5f 1569 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
1570 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
1571 if(!fTRDPIDResponseObject){
1572 AliError(Form("TRD Response not found in run %d", fRun));
59a8e853 1573 }
1574 }
ea235c90 1575}
1576
1577//______________________________________________________________________________
1578void AliPIDResponse::InitializeTRDResponse(){
1579 //
1580 // Set PID Params and references to the TRD PID response
1581 //
db0e2c5f 1582 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
9c499471 1583 fTRDResponse.SetdEdxParams(fTRDdEdxParams);
f2762b1c 1584}
1585
bd58d4b9 1586//______________________________________________________________________________
1587void AliPIDResponse::SetTRDSlices(UInt_t TRDslicesForPID[2],AliTRDPIDResponse::ETRDPIDMethod method) const{
1588
72abc110 1589 if(fLHCperiod.Contains("LHC10D") || fLHCperiod.Contains("LHC10E")){
bd58d4b9 1590 // backward compatibility for setting with 8 slices
1591 TRDslicesForPID[0] = 0;
1592 TRDslicesForPID[1] = 7;
f2762b1c 1593 }
bd58d4b9 1594 else{
1595 if(method==AliTRDPIDResponse::kLQ1D){
1596 TRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
1597 TRDslicesForPID[1] = 0;
1598 }
1599 if(method==AliTRDPIDResponse::kLQ2D){
1600 TRDslicesForPID[0] = 1;
1601 TRDslicesForPID[1] = 7;
1602 }
db0e2c5f 1603 }
bd58d4b9 1604 AliDebug(1,Form("Slice Range set to %d - %d",TRDslicesForPID[0],TRDslicesForPID[1]));
ea235c90 1605}
9c499471 1606//______________________________________________________________________________
1607void AliPIDResponse::SetTRDdEdxParams()
1608{
1609 if(fTRDdEdxParams) return;
1610
1611 const TString containerName = "TRDdEdxParamsContainer";
1612 AliOADBContainer cont(containerName.Data());
1613
1614 const TString filePathNamePackage=Form("%s/COMMON/PID/data/TRDdEdxParams.root", fOADBPath.Data());
1615
1616 const Int_t statusCont = cont.InitFromFile(filePathNamePackage.Data(), cont.GetName());
1617 if (statusCont){
1618 AliFatal("Failed initializing settings from OADB");
1619 }
1620 else{
2ddf4abb 1621 AliInfo(Form("Loading %s from %s\n", cont.GetName(), filePathNamePackage.Data()));
9c499471 1622
1623 fTRDdEdxParams = (AliTRDdEdxParams*)(cont.GetObject(fRun, "default"));
1624 //fTRDdEdxParams->Print();
1625
1626 if(!fTRDdEdxParams){
1627 AliError(Form("TRD dEdx Params default not found"));
1628 }
1629 }
1630}
1631
b79db598 1632//______________________________________________________________________________
1633void AliPIDResponse::SetTOFPidResponseMaster()
1634{
1635 //
1636 // Load the TOF pid params from the OADB
1637 //
00a38d07 1638
1639 if (fTOFPIDParams) delete fTOFPIDParams;
644666df 1640 fTOFPIDParams=NULL;
00a38d07 1641
b79db598 1642 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
00a38d07 1643 if (oadbf && oadbf->IsOpen()) {
b79db598 1644 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
1645 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
00a38d07 1646 if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 1647 oadbf->Close();
1648 delete oadbc;
b79db598 1649 }
1650 delete oadbf;
1651
00a38d07 1652 if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
1653}
b79db598 1654
1655//______________________________________________________________________________
1656void AliPIDResponse::InitializeTOFResponse(){
1657 //
1658 // Set PID Params to the TOF PID response
00a38d07 1659 //
1660
1661 AliInfo("TOF PID Params loaded from OADB");
1662 AliInfo(Form(" TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
1663 AliInfo(Form(" StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
1664 AliInfo(Form(" TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
1665 fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
c53e310b 1666 AliInfo(Form(" Fraction of tracks within gaussian behaviour: %6.4f",fTOFPIDParams->GetTOFtail()));
1667 AliInfo(Form(" MC: Fraction of tracks (percentage) to cut to fit matching in data: %6.2f%%",fTOFPIDParams->GetTOFmatchingLossMC()));
1668 AliInfo(Form(" MC: Fraction of random hits (percentage) to add to fit mismatch in data: %6.2f%%",fTOFPIDParams->GetTOFadditionalMismForMC()));
1669 AliInfo(Form(" Start Time Offset %6.2f ps",fTOFPIDParams->GetTOFtimeOffset()));
1670
b79db598 1671 for (Int_t i=0;i<4;i++) {
1672 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
1673 }
1674 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
1675
78cbd205 1676 AliInfo("TZERO resolution loaded from ESDrun/AODheader");
1677 Float_t t0Spread[4];
1678 for (Int_t i=0;i<4;i++) t0Spread[i]=fCurrentEvent->GetT0spread(i);
1679 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]));
1680 Float_t a = t0Spread[1]*t0Spread[1]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1681 Float_t c = t0Spread[2]*t0Spread[2]-t0Spread[0]*t0Spread[0]+t0Spread[3]*t0Spread[3];
1682 if ( (t0Spread[0] > 50. && t0Spread[0] < 400.) && (a > 0.) && (c>0.)) {
1683 fResT0AC=t0Spread[3];
1684 fResT0A=TMath::Sqrt(a);
1685 fResT0C=TMath::Sqrt(c);
1686 } else {
1687 AliInfo(" TZERO spreads not present or inconsistent, loading default");
1688 fResT0A=75.;
1689 fResT0C=65.;
1690 fResT0AC=55.;
1691 }
1692 AliInfo(Form(" TZERO resolution set to: T0A: %f [ps] T0C: %f [ps] T0AC %f [ps]",fResT0A,fResT0C,fResT0AC));
1693
b79db598 1694}
1695
567624b5 1696//______________________________________________________________________________
1697void AliPIDResponse::SetHMPIDPidResponseMaster()
1698{
1699 //
1700 // Load the HMPID pid params from the OADB
1701 //
1702
1703 if (fHMPIDPIDParams) delete fHMPIDPIDParams;
1704 fHMPIDPIDParams=NULL;
1705
b2f22270 1706 TFile *oadbf;
1707 if(!fIsMC) oadbf = new TFile(Form("%s/COMMON/PID/data/HMPIDPIDParams.root",fOADBPath.Data()));
1708 else oadbf = new TFile(Form("%s/COMMON/PID/MC/HMPIDPIDParams.root",fOADBPath.Data()));
567624b5 1709 if (oadbf && oadbf->IsOpen()) {
1710 AliInfo(Form("Loading HMPID Params from %s/COMMON/PID/data/HMPIDPIDParams.root", fOADBPath.Data()));
1711 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("HMPoadb");
1712 if (oadbc) fHMPIDPIDParams = dynamic_cast<AliHMPIDPIDParams *>(oadbc->GetObject(fRun,"HMPparams"));
1713 oadbf->Close();
1714 delete oadbc;
1715 }
1716 delete oadbf;
1717
1718 if (!fHMPIDPIDParams) AliFatal("HMPIDPIDParams could not be retrieved");
1719}
1720
1721//______________________________________________________________________________
1722void AliPIDResponse::InitializeHMPIDResponse(){
1723 //
1724 // Set PID Params to the HMPID PID response
1725 //
1726
1727 fHMPIDResponse.SetRefIndexArray(fHMPIDPIDParams->GetHMPIDrefIndex());
1728}
b79db598 1729
1c9d11be 1730//______________________________________________________________________________
239fe91c 1731Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
1732 // old function for compatibility
1733 Int_t ntracklets=0;
1734 return IdentifiedAsElectronTRD(vtrack,ntracklets,efficiencyLevel,centrality,PIDmethod);
1735}
1736
1737//______________________________________________________________________________
1738Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Int_t &ntracklets,Double_t efficiencyLevel,Double_t centrality,AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const {
ea235c90 1739 //
1740 // Check whether track is identified as electron under a given electron efficiency hypothesis
bd58d4b9 1741 //
239fe91c 1742 // ntracklets is the number of tracklets that has been used to calculate the PID signal
bd58d4b9 1743
ea235c90 1744 Double_t probs[AliPID::kSPECIES];
ea235c90 1745
239fe91c 1746 ntracklets =CalculateTRDResponse(vtrack,probs,PIDmethod);
1747
99e9d5ec 1748 // Take mean of the TRD momenta in the given tracklets
1749 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
1750 Int_t nmomenta = 0;
ea235c90 1751 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
1752 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 1753 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 1754 }
1755 }
99e9d5ec 1756 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 1757
bd58d4b9 1758 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel,centrality,PIDmethod);
ea235c90 1759}
1760
b2138b40 1761//______________________________________________________________________________
1762void AliPIDResponse::SetEMCALPidResponseMaster()
1763{
1764 //
1765 // Load the EMCAL pid response functions from the OADB
1766 //
1767 TObjArray* fEMCALPIDParamsRun = NULL;
1768 TObjArray* fEMCALPIDParamsPass = NULL;
1769
1770 if(fEMCALPIDParams) return;
1771 AliOADBContainer contParams("contParams");
1772
1773 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
1774 if(statusPars){
1775 AliError("Failed initializing PID Params from OADB");
1776 }
1777 else {
1778 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
1779
1780 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
1781 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
1782 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1783
1784 if(!fEMCALPIDParams){
f8d39067 1785 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 1786 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 1787
1f631618 1788 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1789 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 1790 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1791
1792 if(!fEMCALPIDParams){
1f631618 1793 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 1794 }
1795 }
1796 }
1797}
1798
1799//______________________________________________________________________________
1800void AliPIDResponse::InitializeEMCALResponse(){
1801 //
1802 // Set PID Params to the EMCAL PID response
1803 //
1804 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1805
1806}
00a38d07 1807
1c9d11be 1808//______________________________________________________________________________
1809void AliPIDResponse::FillTrackDetectorPID(const AliVTrack *track, EDetector detector) const
00a38d07 1810{
1811 //
1812 // create detector PID information and setup the transient pointer in the track
1813 //
1c9d11be 1814
1815 // check if detector number is inside accepted range
1816 if (detector == kNdetectors) return;
1817
1818 // get detector pid
1819 AliDetectorPID *detPID=const_cast<AliDetectorPID*>(track->GetDetectorPID());
1820 if (!detPID) {
1821 detPID=new AliDetectorPID;
1822 (const_cast<AliVTrack*>(track))->SetDetectorPID(detPID);
1823 }
1824
1825 //check if values exist
355b831b 1826 if (detPID->HasRawProbability(detector) && detPID->HasNumberOfSigmas(detector)) return;
00a38d07 1827
1828 //TODO: which particles to include? See also the loops below...
1829 Double_t values[AliPID::kSPECIESC]={0};
1c9d11be 1830
355b831b 1831 //probabilities
1832 EDetPidStatus status=GetComputePIDProbability(detector,track,AliPID::kSPECIESC,values);
1833 detPID->SetRawProbability(detector, values, (Int_t)AliPID::kSPECIESC, status);
1834
1c9d11be 1835 //nsigmas
1836 for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
1837 values[ipart]=GetNumberOfSigmas(detector,track,(AliPID::EParticleType)ipart);
355b831b 1838 // the pid status is the same for probabilities and nSigmas, so it is
1839 // fine to use the one from the probabilities also here
1840 detPID->SetNumberOfSigmas(detector, values, (Int_t)AliPID::kSPECIESC, status);
1c9d11be 1841
1c9d11be 1842}
1843
1844//______________________________________________________________________________
1845void AliPIDResponse::FillTrackDetectorPID()
1846{
1847 //
1848 // create detector PID information and setup the transient pointer in the track
1849 //
1850
1851 if (!fCurrentEvent) return;
00a38d07 1852
1853 for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
1854 AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
1855 if (!track) continue;
1856
00a38d07 1857 for (Int_t idet=0; idet<kNdetectors; ++idet){
1c9d11be 1858 FillTrackDetectorPID(track, (EDetector)idet);
00a38d07 1859 }
00a38d07 1860 }
1861}
1862
1c9d11be 1863//______________________________________________________________________________
5f8db5fe 1864void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1865 //
1866 // Set TOF response function
1867 // Input option for event_time used
1868 //
c53e310b 1869
5f8db5fe 1870 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1871 if(t0spread < 10) t0spread = 80;
1872
c53e310b 1873 // T0-FILL and T0-TO offset (because of TOF misallignment
1874 Float_t starttimeoffset = 0;
1875 if(fTOFPIDParams && !(fIsMC)) starttimeoffset=fTOFPIDParams->GetTOFtimeOffset();
b3f687a1 1876 if(fTOFPIDParams){
1877 fTOFtail = fTOFPIDParams->GetTOFtail();
1878 GetTOFResponse().SetTOFtail(fTOFtail);
1879 }
5f8db5fe 1880
c53e310b 1881 // T0 from TOF algorithm
5f8db5fe 1882 Bool_t flagT0TOF=kFALSE;
1883 Bool_t flagT0T0=kFALSE;
1884 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1885 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1886 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1887
1888 // T0-TOF arrays
1889 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1890 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1891 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1892 estimatedT0event[i]=0.0;
1893 estimatedT0resolution[i]=0.0;
1894 startTimeMask[i] = 0;
1895 }
1896
78cbd205 1897 Float_t resT0A=fResT0A;
1898 Float_t resT0C=fResT0C;
1899 Float_t resT0AC=fResT0AC;
5f8db5fe 1900 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1901 flagT0T0=kTRUE;
1902 }
1903
1904
1905 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1906
1907 if (tofHeader) { // read global info and T0-TOF
6aff3c0b 1908 // fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution()); // read from OADB in the initialization
5f8db5fe 1909 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1910 if(t0spread < 10) t0spread = 80;
1911
1912 flagT0TOF=kTRUE;
1913 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1914 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1915 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1916 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
c53e310b 1917
1918 if(startTimeRes[i] > t0spread - 10 && TMath::Abs(startTime[i]) < 0.001) startTime[i] = -starttimeoffset; // apply offset for T0-fill
5f8db5fe 1919 }
1920
1921 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1922 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1923 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1924 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1925 Int_t icurrent = (Int_t)ibin->GetAt(j);
1926 startTime[icurrent]=t0Bin->GetAt(j);
1927 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1928 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
c53e310b 1929 if(startTimeRes[icurrent] > t0spread - 10 && TMath::Abs(startTime[icurrent]) < 0.001) startTime[icurrent] = -starttimeoffset; // apply offset for T0-fill
5f8db5fe 1930 }
1931 }
1932
1933 // for cut of 3 sigma on t0 spread
1934 Float_t t0cut = 3 * t0spread;
1935 if(t0cut < 500) t0cut = 500;
1936
1937 if(option == kFILL_T0){ // T0-FILL is used
1938 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
c53e310b 1939 estimatedT0event[i]=0.0-starttimeoffset;
5f8db5fe 1940 estimatedT0resolution[i]=t0spread;
1941 }
1942 fTOFResponse.SetT0event(estimatedT0event);
1943 fTOFResponse.SetT0resolution(estimatedT0resolution);
1944 }
1945
1946 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1947 if(flagT0TOF){
1948 fTOFResponse.SetT0event(startTime);
1949 fTOFResponse.SetT0resolution(startTimeRes);
1950 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1951 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1952 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1953 }
1954 }
1955 else{
1956 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
c53e310b 1957 estimatedT0event[i]=0.0-starttimeoffset;
5f8db5fe 1958 estimatedT0resolution[i]=t0spread;
1959 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1960 }
1961 fTOFResponse.SetT0event(estimatedT0event);
1962 fTOFResponse.SetT0resolution(estimatedT0resolution);
1963 }
1964 }
1965 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1966 Float_t t0AC=-10000;
1967 Float_t t0A=-10000;
1968 Float_t t0C=-10000;
1969 if(flagT0T0){
c53e310b 1970 t0A= vevent->GetT0TOF()[1] - starttimeoffset;
1971 t0C= vevent->GetT0TOF()[2] - starttimeoffset;
f84b18dd 1972 // t0AC= vevent->GetT0TOF()[0];
c53e310b 1973 t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
6aff3c0b 1974 resT0AC= 1./TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
1975 t0AC *= resT0AC*resT0AC;
5f8db5fe 1976 }
1977
1978 Float_t t0t0Best = 0;
1979 Float_t t0t0BestRes = 9999;
1980 Int_t t0used=0;
1981 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1982 t0t0Best = t0AC;
1983 t0t0BestRes = resT0AC;
1984 t0used=6;
1985 }
1986 else if(TMath::Abs(t0C) < t0cut){
1987 t0t0Best = t0C;
1988 t0t0BestRes = resT0C;
1989 t0used=4;
1990 }
1991 else if(TMath::Abs(t0A) < t0cut){
1992 t0t0Best = t0A;
1993 t0t0BestRes = resT0A;
1994 t0used=2;
1995 }
1996
1997 if(flagT0TOF){ // if T0-TOF info is available
1998 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1999 if(t0t0BestRes < 999){
2000 if(startTimeRes[i] < t0spread){
2001 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
2002 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
2003 estimatedT0event[i]=t0best / wtot;
2004 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
2005 startTimeMask[i] = t0used+1;
2006 }
2007 else {
2008 estimatedT0event[i]=t0t0Best;
2009 estimatedT0resolution[i]=t0t0BestRes;
2010 startTimeMask[i] = t0used;
2011 }
2012 }
2013 else{
2014 estimatedT0event[i]=startTime[i];
2015 estimatedT0resolution[i]=startTimeRes[i];
2016 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
2017 }
2018 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
2019 }
2020 fTOFResponse.SetT0event(estimatedT0event);
2021 fTOFResponse.SetT0resolution(estimatedT0resolution);
2022 }
2023 else{ // if no T0-TOF info is available
2024 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2025 fTOFResponse.SetT0binMask(i,t0used);
2026 if(t0t0BestRes < 999){
2027 estimatedT0event[i]=t0t0Best;
2028 estimatedT0resolution[i]=t0t0BestRes;
2029 }
2030 else{
c53e310b 2031 estimatedT0event[i]=0.0-starttimeoffset;
5f8db5fe 2032 estimatedT0resolution[i]=t0spread;
2033 }
2034 }
2035 fTOFResponse.SetT0event(estimatedT0event);
2036 fTOFResponse.SetT0resolution(estimatedT0resolution);
2037 }
2038 }
2039
2040 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
2041 Float_t t0AC=-10000;
2042 Float_t t0A=-10000;
2043 Float_t t0C=-10000;
2044 if(flagT0T0){
c53e310b 2045 t0A= vevent->GetT0TOF()[1] - starttimeoffset;
2046 t0C= vevent->GetT0TOF()[2] - starttimeoffset;
f84b18dd 2047 // t0AC= vevent->GetT0TOF()[0];
c53e310b 2048 t0AC= t0A/resT0A/resT0A + t0C/resT0C/resT0C;
6aff3c0b 2049 resT0AC= 1./TMath::Sqrt(1./resT0A/resT0A + 1./resT0C/resT0C);
2050 t0AC *= resT0AC*resT0AC;
5f8db5fe 2051 }
2052
2053 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
2054 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2055 estimatedT0event[i]=t0AC;
2056 estimatedT0resolution[i]=resT0AC;
2057 fTOFResponse.SetT0binMask(i,6);
2058 }
2059 }
2060 else if(TMath::Abs(t0C) < t0cut){
2061 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2062 estimatedT0event[i]=t0C;
2063 estimatedT0resolution[i]=resT0C;
2064 fTOFResponse.SetT0binMask(i,4);
2065 }
2066 }
2067 else if(TMath::Abs(t0A) < t0cut){
2068 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
2069 estimatedT0event[i]=t0A;
2070 estimatedT0resolution[i]=resT0A;
2071 fTOFResponse.SetT0binMask(i,2);
2072 }
2073 }
2074 else{
2075 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
c53e310b 2076 estimatedT0event[i]= 0.0 - starttimeoffset;
5f8db5fe 2077 estimatedT0resolution[i]=t0spread;
2078 fTOFResponse.SetT0binMask(i,0);
2079 }
2080 }
2081 fTOFResponse.SetT0event(estimatedT0event);
2082 fTOFResponse.SetT0resolution(estimatedT0resolution);
2083 }
c53e310b 2084
5f8db5fe 2085 delete [] startTime;
2086 delete [] startTimeRes;
2087 delete [] startTimeMask;
2088 delete [] estimatedT0event;
2089 delete [] estimatedT0resolution;
2090}
1c9d11be 2091
2092//______________________________________________________________________________
2093// private non cached versions of the PID calculation
2094//
2095
2096
2097//______________________________________________________________________________
355b831b 2098Float_t AliPIDResponse::GetNumberOfSigmas(EDetector detector, const AliVParticle *vtrack, AliPID::EParticleType type) const
1c9d11be 2099{
2100 //
2101 // NumberOfSigmas for 'detCode'
2102 //
355b831b 2103
2104 const AliVTrack *track=static_cast<const AliVTrack*>(vtrack);
1c9d11be 2105
355b831b 2106 switch (detector){
567624b5 2107 case kITS: return GetNumberOfSigmasITS(track, type); break;
2108 case kTPC: return GetNumberOfSigmasTPC(track, type); break;
5cd0300d 2109 case kTRD: return GetNumberOfSigmasTRD(track, type); break;
567624b5 2110 case kTOF: return GetNumberOfSigmasTOF(track, type); break;
2111 case kHMPID: return GetNumberOfSigmasHMPID(track, type); break;
1c9d11be 2112 case kEMCAL: return GetNumberOfSigmasEMCAL(track, type); break;
2113 default: return -999.;
2114 }
1c9d11be 2115
355b831b 2116 return -999.;
2117}
1c9d11be 2118
2119//______________________________________________________________________________
2120Float_t AliPIDResponse::GetNumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
2121{
2122 //
2123 // Calculate the number of sigmas in the ITS
2124 //
2125
2126 AliVTrack *track=(AliVTrack*)vtrack;
355b831b 2127
2128 const EDetPidStatus pidStatus=GetITSPIDStatus(track);
2129 if (pidStatus!=kDetPidOk) return -999.;
355b831b 2130
567624b5 2131 return fITSResponse.GetNumberOfSigmas(track,type);
1c9d11be 2132}
2133
2134//______________________________________________________________________________
2135Float_t AliPIDResponse::GetNumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
2136{
2137 //
2138 // Calculate the number of sigmas in the TPC
2139 //
2140
2141 AliVTrack *track=(AliVTrack*)vtrack;
355b831b 2142
2143 const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
a017c06a 2144 if (pidStatus==kDetNoSignal) return -999.;
1d59271b 2145
2146 // the following call is needed in order to fill the transient data member
2147 // fTPCsignalTuned which is used in the TPCPIDResponse to judge
2148 // if using tuned on data
87da0205 2149 if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
2150 this->GetTPCsignalTunedOnData(track);
1c9d11be 2151
87da0205 2152 return fTPCResponse.GetNumberOfSigmas(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
1c9d11be 2153}
2154
5cd0300d 2155//______________________________________________________________________________
2156Float_t AliPIDResponse::GetNumberOfSigmasTRD(const AliVParticle *vtrack, AliPID::EParticleType type) const
2157{
2158 //
2159 // Calculate the number of sigmas in the TRD
2160 //
2161
2162 AliVTrack *track=(AliVTrack*)vtrack;
2163
2164 const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
2165 if (pidStatus!=kDetPidOk) return -999.;
2166
2167 return fTRDResponse.GetNumberOfSigmas(track,type);
2168}
2169
1c9d11be 2170//______________________________________________________________________________
355b831b 2171Float_t AliPIDResponse::GetNumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
1c9d11be 2172{
2173 //
355b831b 2174 // Calculate the number of sigmas in the TOF
1c9d11be 2175 //
2176
2177 AliVTrack *track=(AliVTrack*)vtrack;
355b831b 2178
2179 const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
2180 if (pidStatus!=kDetPidOk) return -999.;
1c9d11be 2181
355b831b 2182 return GetNumberOfSigmasTOFold(vtrack, type);
2183}
567624b5 2184//______________________________________________________________________________
2185
2186Float_t AliPIDResponse::GetNumberOfSigmasHMPID(const AliVParticle *vtrack, AliPID::EParticleType type) const
2187{
2188 //
2189 // Calculate the number of sigmas in the HMPID
2190 //
2191 AliVTrack *track=(AliVTrack*)vtrack;
2192
2193 const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2194 if (pidStatus!=kDetPidOk) return -999.;
2195
2196 return fHMPIDResponse.GetNumberOfSigmas(track, type);
2197}
355b831b 2198
2199//______________________________________________________________________________
2200Float_t AliPIDResponse::GetNumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
2201{
2202 //
2203 // Calculate the number of sigmas in the EMCAL
2204 //
1c9d11be 2205
355b831b 2206 AliVTrack *track=(AliVTrack*)vtrack;
2207
2208 const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
2209 if (pidStatus!=kDetPidOk) return -999.;
2210
2211 const Int_t nMatchClus = track->GetEMCALcluster();
2212 AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1c9d11be 2213
355b831b 2214 const Double_t mom = track->P();
2215 const Double_t pt = track->Pt();
2216 const Int_t charge = track->Charge();
2217 const Double_t fClsE = matchedClus->E();
2218 const Double_t EovP = fClsE/mom;
1c9d11be 2219
355b831b 2220 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
1c9d11be 2221}
2222
567624b5 2223//______________________________________________________________________________
1d59271b 2224AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaITS(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
567624b5 2225{
2226 //
2227 // Signal minus expected Signal for ITS
2228 //
2229 AliVTrack *track=(AliVTrack*)vtrack;
1d59271b 2230 val=fITSResponse.GetSignalDelta(track,type,ratio);
567624b5 2231
2232 return GetITSPIDStatus(track);
2233}
2234
2235//______________________________________________________________________________
1d59271b 2236AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTPC(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
567624b5 2237{
2238 //
2239 // Signal minus expected Signal for TPC
2240 //
2241 AliVTrack *track=(AliVTrack*)vtrack;
1d59271b 2242
2243 // the following call is needed in order to fill the transient data member
2244 // fTPCsignalTuned which is used in the TPCPIDResponse to judge
2245 // if using tuned on data
87da0205 2246 if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC))
1d59271b 2247 this->GetTPCsignalTunedOnData(track);
2248
87da0205 2249 val=fTPCResponse.GetSignalDelta(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection, ratio);
567624b5 2250
2251 return GetTPCPIDStatus(track);
2252}
2253
5cd0300d 2254//______________________________________________________________________________
2255AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTRD(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
2256{
2257 //
2258 // Signal minus expected Signal for TRD
2259 //
2260 AliVTrack *track=(AliVTrack*)vtrack;
2261 val=fTRDResponse.GetSignalDelta(track,type,ratio);
2262
2263 return GetTRDPIDStatus(track);
2264}
2265
567624b5 2266//______________________________________________________________________________
1d59271b 2267AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaTOF(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
567624b5 2268{
2269 //
2270 // Signal minus expected Signal for TOF
2271 //
2272 AliVTrack *track=(AliVTrack*)vtrack;
1d59271b 2273 val=GetSignalDeltaTOFold(track, type, ratio);
87da0205 2274
567624b5 2275 return GetTOFPIDStatus(track);
2276}
2277
2278//______________________________________________________________________________
1d59271b 2279AliPIDResponse::EDetPidStatus AliPIDResponse::GetSignalDeltaHMPID(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &val, Bool_t ratio/*=kFALSE*/) const
567624b5 2280{
2281 //
2282 // Signal minus expected Signal for HMPID
2283 //
2284 AliVTrack *track=(AliVTrack*)vtrack;
1d59271b 2285 val=fHMPIDResponse.GetSignalDelta(track, type, ratio);
567624b5 2286
2287 return GetHMPIDPIDStatus(track);
2288}
1c9d11be 2289
2290//______________________________________________________________________________
2291AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePIDProbability (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2292{
2293 //
2294 // Compute PID response of 'detCode'
2295 //
2296
2297 switch (detCode){
2298 case kITS: return GetComputeITSProbability(track, nSpecies, p); break;
2299 case kTPC: return GetComputeTPCProbability(track, nSpecies, p); break;
2300 case kTRD: return GetComputeTRDProbability(track, nSpecies, p); break;
2301 case kTOF: return GetComputeTOFProbability(track, nSpecies, p); break;
2302 case kPHOS: return GetComputePHOSProbability(track, nSpecies, p); break;
2303 case kEMCAL: return GetComputeEMCALProbability(track, nSpecies, p); break;
2304 case kHMPID: return GetComputeHMPIDProbability(track, nSpecies, p); break;
2305 default: return kDetNoSignal;
2306 }
2307}
2308
2309//______________________________________________________________________________
2310AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2311{
2312 //
2313 // Compute PID response for the ITS
2314 //
2315
1c9d11be 2316 // set flat distribution (no decision)
2317 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2318
355b831b 2319 const EDetPidStatus pidStatus=GetITSPIDStatus(track);
2320 if (pidStatus!=kDetPidOk) return pidStatus;
2321
2322 if (track->GetDetectorPID()){
2323 return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
2324 }
1c9d11be 2325
2326 //check for ITS standalone tracks
2327 Bool_t isSA=kTRUE;
2328 if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
2329
2330 Double_t mom=track->P();
2331 Double_t dedx=track->GetITSsignal();
2332 Double_t momITS=mom;
2333 UChar_t clumap=track->GetITSClusterMap();
2334 Int_t nPointsForPid=0;
2335 for(Int_t i=2; i<6; i++){
2336 if(clumap&(1<<i)) ++nPointsForPid;
2337 }
2338
1c9d11be 2339 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
bf26ce58 2340 for (Int_t j=0; j<nSpecies; j++) {
1c9d11be 2341 const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
1c9d11be 2342 //TODO: in case of the electron, use the SA parametrisation,
2343 // this needs to be changed if ITS provides a parametrisation
2344 // for electrons also for ITS+TPC tracks
5f42450e 2345 Double_t bethe=fITSResponse.Bethe(momITS,(AliPID::EParticleType)j,isSA || (j==(Int_t)AliPID::kElectron))*chargeFactor;
1c9d11be 2346 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
04f7af3e 2347 Double_t nSigma=fITSResponse.GetNumberOfSigmas(track, (AliPID::EParticleType)j);
2348 if (TMath::Abs(nSigma) > fRange) {
1c9d11be 2349 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
2350 } else {
04f7af3e 2351 p[j]=TMath::Exp(-0.5*nSigma*nSigma)/sigma;
1c9d11be 2352 mismatch=kFALSE;
2353 }
1c9d11be 2354 }
2355
2356 if (mismatch){
bf26ce58 2357 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1c9d11be 2358 }
2359
1c9d11be 2360 return kDetPidOk;
2361}
2362//______________________________________________________________________________
2363AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2364{
2365 //
2366 // Compute PID response for the TPC
2367 //
2368
2369 // set flat distribution (no decision)
2370 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2371
355b831b 2372 const EDetPidStatus pidStatus=GetTPCPIDStatus(track);
a017c06a 2373 if (pidStatus==kDetNoSignal) return pidStatus;
1c9d11be 2374
1c9d11be 2375 Double_t dedx=track->GetTPCsignal();
2376 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
2377
87da0205 2378 if (fTuneMConData && ((fTuneMConDataMask & kDetTPC) == kDetTPC)) dedx = this->GetTPCsignalTunedOnData(track);
1c9d11be 2379
f84b18dd 2380 Double_t bethe = 0.;
2381 Double_t sigma = 0.;
2382
bf26ce58 2383 for (Int_t j=0; j<nSpecies; j++) {
1c9d11be 2384 AliPID::EParticleType type=AliPID::EParticleType(j);
f84b18dd 2385
87da0205 2386 bethe=fTPCResponse.GetExpectedSignal(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
2387 sigma=fTPCResponse.GetExpectedSigma(track, type, AliTPCPIDResponse::kdEdxDefault, fUseTPCEtaCorrection, fUseTPCMultiplicityCorrection);
f85a3764 2388
1c9d11be 2389 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
2390 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
2391 } else {
2392 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
2393 mismatch=kFALSE;
2394 }
2395 }
2396
2397 if (mismatch){
2398 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
1c9d11be 2399 }
2400
a017c06a 2401 return pidStatus;
1c9d11be 2402}
2403//______________________________________________________________________________
506c1482 2404AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],Bool_t kNoMism) const
1c9d11be 2405{
2406 //
57e985ed 2407 // Compute PID probabilities for TOF
1c9d11be 2408 //
4fc0f532 2409
42fcc729 2410 fgTOFmismatchProb = 1E-8;
2411
2412 // centrality --> fCurrCentrality
2413 // Beam type --> fBeamTypeNum
2414 // N TOF cluster --> TOF header --> to get the TOF header we need to add a virtual method in AliVTrack extended to ESD and AOD tracks
2415 // isMC --> fIsMC
506c1482 2416 Float_t pt = track->Pt();
2417 Float_t mismPropagationFactor[10] = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
6839b659 2418 if(! (kNoMism | fNoTOFmism)){ // this flag allows to disable mismatch for iterative procedure to get priors
506c1482 2419 mismPropagationFactor[3] = 1 + TMath::Exp(1 - 1.12*pt);// it has to be alligned with the one in AliPIDCombined
2420 mismPropagationFactor[4] = 1 + 1./(4.71114 - 5.72372*pt + 2.94715*pt*pt);// it has to be alligned with the one in AliPIDCombined
2421
42fcc729 2422 Int_t nTOFcluster = 0;
6aff3c0b 2423 if(track->GetTOFHeader() && track->GetTOFHeader()->GetTriggerMask() && track->GetTOFHeader()->GetNumberOfTOFclusters() > -1){ // N TOF clusters available
42fcc729 2424 nTOFcluster = track->GetTOFHeader()->GetNumberOfTOFclusters();
506c1482 2425 if(fIsMC) nTOFcluster = Int_t(nTOFcluster * 1.5); // +50% in MC
42fcc729 2426 }
2427 else{
2428 switch(fBeamTypeNum){
506c1482 2429 case kPP: // pp
2430 nTOFcluster = 80;
2431 break;
2432 case kPPB: // pPb 5.05 ATeV
2433 nTOFcluster = Int_t(308 - 2.12*fCurrCentrality + TMath::Exp(4.917 -0.1604*fCurrCentrality));
2434 break;
2435 case kPBPB: // PbPb 2.76 ATeV
2436 nTOFcluster = Int_t(TMath::Exp(9.4 - 0.022*fCurrCentrality));
2437 break;
2438 }
2439 }
2440
2441 switch(fBeamTypeNum){ // matching window factors for 3 cm and 10 cm (about (10/3)^2)
42fcc729 2442 case kPP: // pp 7 TeV
506c1482 2443 nTOFcluster *= 10;
42fcc729 2444 break;
2445 case kPPB: // pPb 5.05 ATeV
506c1482 2446 nTOFcluster *= 10;
42fcc729 2447 break;
506c1482 2448 case kPBPB: // pPb 5.05 ATeV
2449 // nTOFcluster *= 1;
42fcc729 2450 break;
2451 }
506c1482 2452
6aff3c0b 2453 if(nTOFcluster < 0) nTOFcluster = 10;
2454
506c1482 2455
2456 fgTOFmismatchProb=fTOFResponse.GetMismatchProbability(track->GetTOFsignal(),track->Eta()) * nTOFcluster *6E-6 * (1 + 2.90505e-01/pt/pt); // mism weight * tof occupancy (including matching window factor) * pt dependence
2457
42fcc729 2458 }
2459
1c9d11be 2460 // set flat distribution (no decision)
2461 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2462
355b831b 2463 const EDetPidStatus pidStatus=GetTOFPIDStatus(track);
2464 if (pidStatus!=kDetPidOk) return pidStatus;
2465
c5fb644a 2466 const Double_t meanCorrFactor = 0.07/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
1c9d11be 2467
bf26ce58 2468 for (Int_t j=0; j<nSpecies; j++) {
1c9d11be 2469 AliPID::EParticleType type=AliPID::EParticleType(j);
355b831b 2470 const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
1c9d11be 2471
355b831b 2472 const Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
2473 const Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
4fc0f532 2474
2475 if(nsigmas < fTOFtail)
2476 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
2477 else
2478 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
2479
506c1482 2480 p[j] += fgTOFmismatchProb*mismPropagationFactor[j];
1c9d11be 2481 }
2482
1c9d11be 2483 return kDetPidOk;
2484}
239fe91c 2485
2486Int_t AliPIDResponse::CalculateTRDResponse(const AliVTrack *track,Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
2487{
2488 // new function for backward compatibility
2489 // returns number of tracklets PID
2490
2491 UInt_t TRDslicesForPID[2];
2492 SetTRDSlices(TRDslicesForPID,PIDmethod);
2493
2494 Float_t mom[6]={0.};
2495 Double_t dedx[48]={0.}; // Allocate space for the maximum number of TRD slices
2496 Int_t nslices = TRDslicesForPID[1] - TRDslicesForPID[0] + 1;
2497 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", TRDslicesForPID[0], TRDslicesForPID[1], nslices));
2498 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
2499 mom[ilayer] = track->GetTRDmomentum(ilayer);
2500 for(UInt_t islice = TRDslicesForPID[0]; islice <= TRDslicesForPID[1]; islice++){
2501 dedx[ilayer*nslices+islice-TRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
2502 }
2503 }
2504
2505 return fTRDResponse.GetResponse(nslices, dedx, mom, p,PIDmethod);
2506
2507}
1c9d11be 2508//______________________________________________________________________________
239fe91c 2509AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[],AliTRDPIDResponse::ETRDPIDMethod PIDmethod) const
1c9d11be 2510{
2511 //
355b831b 2512 // Compute PID probabilities for the TRD
1c9d11be 2513 //
2514
1c9d11be 2515 // set flat distribution (no decision)
2516 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
355b831b 2517
2518 const EDetPidStatus pidStatus=GetTRDPIDStatus(track);
2519 if (pidStatus!=kDetPidOk) return pidStatus;
2520
239fe91c 2521 CalculateTRDResponse(track,p,PIDmethod);
2522
1c9d11be 2523 return kDetPidOk;
1c9d11be 2524}
355b831b 2525
1c9d11be 2526//______________________________________________________________________________
2527AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2528{
2529 //
2530 // Compute PID response for the EMCAL
2531 //
2532
2533 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
355b831b 2534
2535 const EDetPidStatus pidStatus=GetEMCALPIDStatus(track);
2536 if (pidStatus!=kDetPidOk) return pidStatus;
2537
2538 const Int_t nMatchClus = track->GetEMCALcluster();
2539 AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
1c9d11be 2540
355b831b 2541 const Double_t mom = track->P();
2542 const Double_t pt = track->Pt();
2543 const Int_t charge = track->Charge();
2544 const Double_t fClsE = matchedClus->E();
2545 const Double_t EovP = fClsE/mom;
1c9d11be 2546
355b831b 2547 // compute the probabilities
2548 fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p);
2549 return kDetPidOk;
1c9d11be 2550}
355b831b 2551
1c9d11be 2552//______________________________________________________________________________
2553AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
2554{
2555 //
2556 // Compute PID response for the PHOS
2557 //
2558
2559 // set flat distribution (no decision)
2560 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
2561 return kDetNoSignal;
2562}
355b831b 2563
1c9d11be 2564//______________________________________________________________________________
2565AliPIDResponse::EDetPidStatus AliPIDResponse::GetComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
2566{
2567 //
2568 // Compute PID response for the HMPID
2569 //
355b831b 2570
1c9d11be 2571 // set flat distribution (no decision)
2572 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
355b831b 2573
2574 const EDetPidStatus pidStatus=GetHMPIDPIDStatus(track);
2575 if (pidStatus!=kDetPidOk) return pidStatus;
1c9d11be 2576
567624b5 2577 fHMPIDResponse.GetProbability(track,nSpecies,p);
2578
1c9d11be 2579 return kDetPidOk;
2580}
355b831b 2581
2582//______________________________________________________________________________
2583AliPIDResponse::EDetPidStatus AliPIDResponse::GetITSPIDStatus(const AliVTrack *track) const
2584{
2585 // compute ITS pid status
2586
2587 // check status bits
2588 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
2589 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
2590
2591 const Float_t dEdx=track->GetITSsignal();
2592 if (dEdx<=0) return kDetNoSignal;
2593
2594 // requite at least 3 pid clusters
2595 const UChar_t clumap=track->GetITSClusterMap();
2596 Int_t nPointsForPid=0;
2597 for(Int_t i=2; i<6; i++){
2598 if(clumap&(1<<i)) ++nPointsForPid;
2599 }
2600
2601 if(nPointsForPid<3) {
2602 return kDetNoSignal;
2603 }
2604
2605 return kDetPidOk;
2606}
2607
2608//______________________________________________________________________________
2609AliPIDResponse::EDetPidStatus AliPIDResponse:: GetTPCPIDStatus(const AliVTrack *track) const
2610{
2611 // compute TPC pid status
2612
2613 // check quality of the track
2614 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
2615
2616 // check pid values
2617 const Double_t dedx=track->GetTPCsignal();
2618 const UShort_t signalN=track->GetTPCsignalN();
2619 if (signalN<10 || dedx<10) return kDetNoSignal;
2620
2621 if (!(fArrPidResponseMaster && fArrPidResponseMaster->At(AliPID::kPion))) return kDetNoParams;
2622
2623 return kDetPidOk;
2624}
2625
2626//______________________________________________________________________________
2627AliPIDResponse::EDetPidStatus AliPIDResponse::GetTRDPIDStatus(const AliVTrack *track) const
2628{
2629 // compute TRD pid status
2630
2631 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
2632 return kDetPidOk;
2633}
2634
2635//______________________________________________________________________________
2636AliPIDResponse::EDetPidStatus AliPIDResponse::GetTOFPIDStatus(const AliVTrack *track) const
2637{
2638 // compute TOF pid status
2639
2640 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
2641 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
2642
2643 return kDetPidOk;
2644}
2645
2646//______________________________________________________________________________
2647Float_t AliPIDResponse::GetTOFMismatchProbability(const AliVTrack *track) const
2648{
2649 // compute mismatch probability cross-checking at 5 sigmas with TPC
2650 // currently just implemented as a 5 sigma compatibility cut
2651
42fcc729 2652 if(!track) return fgTOFmismatchProb;
2653
355b831b 2654 // check pid status
2655 const EDetPidStatus tofStatus=GetTOFPIDStatus(track);
2656 if (tofStatus!=kDetPidOk) return 0.;
2657
2658 //mismatch
2659 const EDetPidStatus tpcStatus=GetTPCPIDStatus(track);
a017c06a 2660 if (tpcStatus==kDetNoSignal) return 0.;
355b831b 2661
2662 const Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
2663 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
2664 for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
2665 AliPID::EParticleType type=AliPID::EParticleType(j);
2666 const Double_t nsigmas=GetNumberOfSigmasTOFold(track,type) + meanCorrFactor;
2667
2668 if (TMath::Abs(nsigmas)<5.){
2669 const Double_t nsigmasTPC=GetNumberOfSigmasTPC(track,type);
2670 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
2671 }
2672 }
2673
2674 if (mismatch){
2675 return 1.;
2676 }
2677
2678 return 0.;
2679}
2680
355b831b 2681//______________________________________________________________________________
2682AliPIDResponse::EDetPidStatus AliPIDResponse:: GetHMPIDPIDStatus(const AliVTrack *track) const
2683{
2684 // compute HMPID pid status
567624b5 2685
2686 Int_t ch = track->GetHMPIDcluIdx()/1000000;
2687 Double_t HMPIDsignal = track->GetHMPIDsignal();
2688
2689 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0 || ch<0 || ch>6 || HMPIDsignal<0) return kDetNoSignal;
2690
355b831b 2691 return kDetPidOk;
2692}
2693
2694//______________________________________________________________________________
2695AliPIDResponse::EDetPidStatus AliPIDResponse:: GetPHOSPIDStatus(const AliVTrack */*track*/) const
2696{
2697 // compute PHOS pid status
2698 return kDetNoSignal;
2699}
2700
2701//______________________________________________________________________________
2702AliPIDResponse::EDetPidStatus AliPIDResponse:: GetEMCALPIDStatus(const AliVTrack *track) const
2703{
2704 // compute EMCAL pid status
2705
2706
2707 // Track matching
2708 const Int_t nMatchClus = track->GetEMCALcluster();
2709 if (nMatchClus<0) return kDetNoSignal;
2710
2711 AliVCluster *matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
2712
2713 if (!(matchedClus && matchedClus->IsEMCAL())) return kDetNoSignal;
2714
2715 const Int_t charge = track->Charge();
2716 if (TMath::Abs(charge)!=1) return kDetNoSignal;
2717
2718 if (!(fEMCALPIDParams && fEMCALPIDParams->At(AliPID::kElectron))) return kDetNoParams;
2719
2720 return kDetPidOk;
2721
2722}
2723
2724//______________________________________________________________________________
2725AliPIDResponse::EDetPidStatus AliPIDResponse::GetPIDStatus(EDetector detector, const AliVTrack *track) const
2726{
2727 //
2728 // check pid status for a track
2729 //
2730
2731 switch (detector){
2732 case kITS: return GetITSPIDStatus(track); break;
2733 case kTPC: return GetTPCPIDStatus(track); break;
2734 case kTRD: return GetTRDPIDStatus(track); break;
2735 case kTOF: return GetTOFPIDStatus(track); break;
2736 case kPHOS: return GetPHOSPIDStatus(track); break;
2737 case kEMCAL: return GetEMCALPIDStatus(track); break;
2738 case kHMPID: return GetHMPIDPIDStatus(track); break;
2739 default: return kDetNoSignal;
2740 }
2741 return kDetNoSignal;
2742
2743}
5a9dc560 2744
2745//______________________________________________________________________________
2746TString AliPIDResponse::GetChecksum(const TObject* obj) const
2747{
2748 // Return the checksum for an object obj (tested to work properly at least for histograms and TSplines).
2749
2750 TString fileName = Form("tempChecksum.C"); // File name must be fixed for data type "TSpline3", since the file name will end up in the file content!
2751
2752 // For parallel processing, a unique file pathname is required. Uniqueness can be guaranteed by using a unique directory name
2753 UInt_t index = 0;
2754 TString uniquePathName = Form("tempChecksum_%u", index);
2755
2756 // To get a unique path name, increase the index until no directory
2757 // of such a name exists.
2758 // NOTE: gSystem->AccessPathName(...) returns kTRUE, if the access FAILED!
2759 while (!gSystem->AccessPathName(uniquePathName.Data()))
2760 uniquePathName = Form("tempChecksum_%u", ++index);
2761
2762 if (gSystem->mkdir(uniquePathName.Data()) < 0) {
2763 AliError("Could not create temporary directory to store temp file for checksum determination!");
2764 return "ERROR";
2765 }
2766
2767 TString option = "";
2768
2769 // Save object as a macro, which will be deleted immediately after the checksum has been computed
2770 // (does not work for desired data types if saved as *.root for some reason) - one only wants to compare the content, not
2771 // the modification time etc. ...
2772 if (dynamic_cast<const TH1*>(obj))
2773 option = "colz"; // Histos need this option, since w/o this option, a counter is added to the filename
2774
2775
2776 // SaveAs must be called with the fixed fileName only, since the first argument goes into the file content
2777 // for some object types. Thus, change the directory, save the file and then go back
2778 TString oldDir = gSystem->pwd();
2779 gSystem->cd(uniquePathName.Data());
2780 obj->SaveAs(fileName.Data(), option.Data());
2781 gSystem->cd(oldDir.Data());
2782
2783 // Use the file to calculate the MD5 checksum
2784 TMD5* md5 = TMD5::FileChecksum(Form("%s/%s", uniquePathName.Data(), fileName.Data()));
2785 TString checksum = md5->AsString();
2786
2787 // Clean up
2788 delete md5;
2789 gSystem->Exec(Form("rm -rf %s", uniquePathName.Data()));
2790
2791 return checksum;
2792}