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