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