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