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