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