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