PID updates
[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>
b79db598 40#include <AliTOFPIDParams.h>
29bf19f2 41
42#include "AliPIDResponse.h"
43
80f28562 44#include "AliCentrality.h"
45
29bf19f2 46ClassImp(AliPIDResponse);
47
4ec8e76d 48AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
49TNamed("PIDResponse","PIDResponse"),
50fITSResponse(isMC),
51fTPCResponse(),
52fTRDResponse(),
53fTOFResponse(),
e96b9916 54fEMCALResponse(),
fd21ec8d 55fRange(5.),
56fITSPIDmethod(kITSTruncMean),
4ec8e76d 57fIsMC(isMC),
58fOADBPath(),
59fBeamType("PP"),
60fLHCperiod(),
61fMCperiodTPC(),
fd21ec8d 62fMCperiodUser(),
ea235c90 63fCurrentFile(),
4ec8e76d 64fRecoPass(0),
fd21ec8d 65fRecoPassUser(-1),
4ec8e76d 66fRun(0),
67fOldRun(0),
68fArrPidResponseMaster(0x0),
69fResolutionCorrection(0x0),
ea235c90 70fTRDPIDParams(0x0),
71fTRDPIDReference(0x0),
0b39f221 72fTOFtail(1.1),
b79db598 73fTOFPIDParams(0x0),
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;
b79db598 97 if (fTOFPIDParams) delete fTOFPIDParams;
4ec8e76d 98}
99
100//______________________________________________________________________________
101AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
102TNamed(other),
103fITSResponse(other.fITSResponse),
104fTPCResponse(other.fTPCResponse),
105fTRDResponse(other.fTRDResponse),
106fTOFResponse(other.fTOFResponse),
e96b9916 107fEMCALResponse(other.fEMCALResponse),
fd21ec8d 108fRange(other.fRange),
109fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 110fIsMC(other.fIsMC),
111fOADBPath(other.fOADBPath),
112fBeamType("PP"),
113fLHCperiod(),
114fMCperiodTPC(),
fd21ec8d 115fMCperiodUser(other.fMCperiodUser),
ea235c90 116fCurrentFile(),
4ec8e76d 117fRecoPass(0),
fd21ec8d 118fRecoPassUser(other.fRecoPassUser),
4ec8e76d 119fRun(0),
120fOldRun(0),
121fArrPidResponseMaster(0x0),
122fResolutionCorrection(0x0),
ea235c90 123fTRDPIDParams(0x0),
124fTRDPIDReference(0x0),
0b39f221 125fTOFtail(1.1),
b79db598 126fTOFPIDParams(0x0),
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);
0b39f221 170 fTOFtail=1.1;
b79db598 171 fTOFPIDParams=0x0;
e96b9916 172 fCurrentEvent=other.fCurrentEvent;
4ec8e76d 173 }
174 return *this;
175}
176
177//______________________________________________________________________________
fd21ec8d 178Float_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;
e96b9916 190// case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
fd21ec8d 191// case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
192 default: return -999.;
193 }
194
195}
196
197//______________________________________________________________________________
e96b9916 198Float_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){
6d0064aa 213
e96b9916 214 mom = track->P();
215 pt = track->Pt();
216 charge = track->Charge();
217
218 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
219
220 if(matchedClus){
221
6d0064aa 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 }
e96b9916 232 }
233 }
6d0064aa 234
e96b9916 235 return -999;
236
237}
238
239//______________________________________________________________________________
6d0064aa 240Float_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//______________________________________________________________________________
fd21ec8d 287AliPIDResponse::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//______________________________________________________________________________
306AliPIDResponse::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//______________________________________________________________________________
361AliPIDResponse::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//______________________________________________________________________________
404AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
405{
406 //
407 // Compute PID response for the
408 //
409
0b39f221 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
fd21ec8d 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);
0b39f221 429 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
430
fd21ec8d 431 Double_t sig = sigma[j];
432 if (TMath::Abs(nsigmas) > (fRange+2)) {
0b39f221 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)) {
fd21ec8d 446 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
447 } else
448 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
0b39f221 449 */
fd21ec8d 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//______________________________________________________________________________
ea235c90 466AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 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;
ea235c90 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;
fd21ec8d 488}
489//______________________________________________________________________________
e96b9916 490AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 491{
492 //
493 // Compute PID response for the EMCAL
494 //
495
e96b9916 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)
fd21ec8d 538 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
539 return kDetNoSignal;
e96b9916 540
fd21ec8d 541}
542//______________________________________________________________________________
543AliPIDResponse::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//______________________________________________________________________________
ea235c90 554AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 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;
ea235c90 562 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
563
564 track->GetHMPIDpid(p);
565
566 return kDetPidOk;
fd21ec8d 567}
568
569//______________________________________________________________________________
4ec8e76d 570void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
571{
572 //
573 // Apply settings for the current event
574 //
575 fRecoPass=pass;
e96b9916 576
577 fCurrentEvent=0x0;
4ec8e76d 578 if (!event) return;
e96b9916 579 fCurrentEvent=event;
4ec8e76d 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
b79db598 594 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
595
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
b79db598 626 SetTOFPidResponseMaster();
627 InitializeTOFResponse();
4ec8e76d 628}
629
630//_____________________________________________________
631Double_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//______________________________________________________________________________
649void 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
53d016dc 663 TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
4ec8e76d 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"; }
99e9d5ec 668 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 669 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
670 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
ea235c90 671 else if (fRun>=136851&&fRun<=139517) {
672 fLHCperiod="LHC10H";
673 fMCperiodTPC="LHC10H8";
674 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
675 fBeamType="PBPB";
676 }
3077a03d 677 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
2e97a211 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"; }
3077a03d 681 else if (fRun>=166529) {
682 fLHCperiod="LHC11H";
683 fMCperiodTPC="LHC11A10";
684 fBeamType="PBPB";
685 }
686
80ab5635 687
688 //exception new pp MC productions from 2011
689 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
4a527e08 690 // exception for 11f1
691 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
4ec8e76d 692}
693
694//______________________________________________________________________________
695void AliPIDResponse::SetITSParametrisation()
696{
697 //
698 // Set the ITS parametrisation
699 //
700}
701
702//______________________________________________________________________________
703void AliPIDResponse::SetTPCPidResponseMaster()
704{
705 //
706 // Load the TPC pid response functions from the OADB
707 //
09b50a42 708 //don't load twice for the moment
709 if (fArrPidResponseMaster) return;
710
711
4ec8e76d 712 //reset the PID response functions
713 delete fArrPidResponseMaster;
714 fArrPidResponseMaster=0x0;
715
716 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
717
ea235c90 718 TFile *f=TFile::Open(fileName.Data());
719 if (f && f->IsOpen() && !f->IsZombie()){
720 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 721 }
ea235c90 722 delete f;
4ec8e76d 723
724 if (!fArrPidResponseMaster){
725 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
726 return;
727 }
728 fArrPidResponseMaster->SetOwner();
729}
730
731//______________________________________________________________________________
732void AliPIDResponse::SetTPCParametrisation()
733{
734 //
735 // Change BB parametrisation for current run
736 //
737
738 if (fLHCperiod.IsNull()) {
739 AliFatal("No period set, not changing parametrisation");
740 return;
741 }
742
743 //
744 // Set default parametrisations for data and MC
745 //
746
747 //data type
748 TString datatype="DATA";
749 //in case of mc fRecoPass is per default 1
750 if (fIsMC) {
751 datatype="MC";
752 fRecoPass=1;
753 }
754
755 //
756 //reset old splines
757 //
758 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
759 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
760 }
4a527e08 761
762 // period
763 TString period=fLHCperiod;
764 if (fIsMC) period=fMCperiodTPC;
765
766 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
767 Bool_t found=kFALSE;
4ec8e76d 768 //
769 //set the new PID splines
770 //
4ec8e76d 771 if (fArrPidResponseMaster){
772 TObject *grAll=0x0;
773 //for MC don't use period information
774// if (fIsMC) period="[A-Z0-9]*";
775 //for MC use MC period information
4ec8e76d 776//pattern for the default entry (valid for all particles)
777 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
778
779 //loop over entries and filter them
780 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
781 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
782 if (responseFunction==0x0) continue;
783 TString responseName=responseFunction->GetName();
784
785 if (!reg.MatchB(responseName)) continue;
786
787 TObjArray *arr=reg.MatchS(responseName);
788 TString particleName=arr->At(1)->GetName();
789 delete arr;
790 if (particleName.IsNull()) continue;
791 if (particleName=="ALL") grAll=responseFunction;
792 else {
793 //find particle id
794 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
795 TString particle=AliPID::ParticleName(ispec);
796 particle.ToUpper();
797 if ( particle == particleName ){
798 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
799 fTPCResponse.SetUseDatabase(kTRUE);
800 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
4a527e08 801 found=kTRUE;
4ec8e76d 802 break;
803 }
804 }
805 }
806 }
807
808 //set default response function to all particles which don't have a specific one
809 if (grAll){
810 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
811 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
812 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
813 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
4a527e08 814 found=kTRUE;
4ec8e76d 815 }
816 }
817 }
818 }
4a527e08 819
820 if (!found){
821 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
822 }
4ec8e76d 823 //
824 // Setup resolution parametrisation
825 //
826
827 //default
828 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
829
830 if (fRun>=122195){
831 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
832 }
23425eb2 833 if (fArrPidResponseMaster)
4ec8e76d 834 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
835
836 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
837}
838
ea235c90 839//______________________________________________________________________________
840void AliPIDResponse::SetTRDPidResponseMaster()
841{
842 //
843 // Load the TRD pid params and references from the OADB
844 //
845 if(fTRDPIDParams) return;
53d016dc 846 AliOADBContainer contParams("contParams");
847
59a8e853 848 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
849 if(statusPars){
850 AliError("Failed initializing PID Params from OADB");
851 } else {
852 AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
b2138b40 853 fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
59a8e853 854 if(!fTRDPIDParams){
855 AliError(Form("TRD Params not found in run %d", fRun));
856 }
857 }
ea235c90 858
53d016dc 859 AliOADBContainer contRefs("contRefs");
59a8e853 860 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
861 if(statusRefs){
862 AliInfo("Failed Loading References for TRD");
863 } else {
864 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
865 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
866 if(!fTRDPIDReference){
867 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
868 }
869 }
ea235c90 870}
871
872//______________________________________________________________________________
873void AliPIDResponse::InitializeTRDResponse(){
874 //
875 // Set PID Params and references to the TRD PID response
876 //
877 fTRDResponse.SetPIDParams(fTRDPIDParams);
878 fTRDResponse.Load(fTRDPIDReference);
879 if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
880 fTRDslicesForPID[0] = 0;
881 fTRDslicesForPID[1] = 7;
882 }
883}
884
b79db598 885//______________________________________________________________________________
886void AliPIDResponse::SetTOFPidResponseMaster()
887{
888 //
889 // Load the TOF pid params from the OADB
890 //
891 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
892 if (oadbf->IsOpen()) {
893 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
894 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
895 if (fTOFPIDParams) delete fTOFPIDParams;
f57f1e18 896 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 897 oadbf->Close();
898 delete oadbc;
899 } else {
900 AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
901 }
902 delete oadbf;
903
904 }
905
906//______________________________________________________________________________
907void AliPIDResponse::InitializeTOFResponse(){
908 //
909 // Set PID Params to the TOF PID response
910 //
b79db598 911 for (Int_t i=0;i<4;i++) {
912 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
913 }
914 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
915
916}
917
918
ea235c90 919//_________________________________________________________________________
920Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
921 //
922 // Check whether track is identified as electron under a given electron efficiency hypothesis
923 //
924 Double_t probs[AliPID::kSPECIES];
925 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
926
99e9d5ec 927 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
928 // Take mean of the TRD momenta in the given tracklets
929 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
930 Int_t nmomenta = 0;
ea235c90 931 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
932 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 933 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 934 }
935 }
99e9d5ec 936 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 937
938 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
939}
940
b2138b40 941//______________________________________________________________________________
942void AliPIDResponse::SetEMCALPidResponseMaster()
943{
944 //
945 // Load the EMCAL pid response functions from the OADB
946 //
947 TObjArray* fEMCALPIDParamsRun = NULL;
948 TObjArray* fEMCALPIDParamsPass = NULL;
949
950 if(fEMCALPIDParams) return;
951 AliOADBContainer contParams("contParams");
952
953 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
954 if(statusPars){
955 AliError("Failed initializing PID Params from OADB");
956 }
957 else {
958 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
959
960 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
961 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
962 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
963
964 if(!fEMCALPIDParams){
f8d39067 965 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 966 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 967
1f631618 968 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
969 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 970 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
971
972 if(!fEMCALPIDParams){
1f631618 973 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 974 }
975 }
976 }
977}
978
979//______________________________________________________________________________
980void AliPIDResponse::InitializeEMCALResponse(){
981 //
982 // Set PID Params to the EMCAL PID response
983 //
984 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
985
986}
5f8db5fe 987//_________________________________________________________________________
988void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
989 //
990 // Set TOF response function
991 // Input option for event_time used
992 //
993
994 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
995 if(t0spread < 10) t0spread = 80;
996
997 // T0 from TOF algorithm
998
999 Bool_t flagT0TOF=kFALSE;
1000 Bool_t flagT0T0=kFALSE;
1001 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1002 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1003 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1004
1005 // T0-TOF arrays
1006 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1007 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1008 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1009 estimatedT0event[i]=0.0;
1010 estimatedT0resolution[i]=0.0;
1011 startTimeMask[i] = 0;
1012 }
1013
1014 Float_t resT0A=75,resT0C=65,resT0AC=55;
1015 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1016 flagT0T0=kTRUE;
1017 }
1018
1019
1020 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1021
1022 if (tofHeader) { // read global info and T0-TOF
1023 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1024 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1025 if(t0spread < 10) t0spread = 80;
1026
1027 flagT0TOF=kTRUE;
1028 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1029 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1030 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1031 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1032 }
1033
1034 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1035 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1036 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1037 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1038 Int_t icurrent = (Int_t)ibin->GetAt(j);
1039 startTime[icurrent]=t0Bin->GetAt(j);
1040 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1041 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1042 }
1043 }
1044
1045 // for cut of 3 sigma on t0 spread
1046 Float_t t0cut = 3 * t0spread;
1047 if(t0cut < 500) t0cut = 500;
1048
1049 if(option == kFILL_T0){ // T0-FILL is used
1050 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1051 estimatedT0event[i]=0.0;
1052 estimatedT0resolution[i]=t0spread;
1053 }
1054 fTOFResponse.SetT0event(estimatedT0event);
1055 fTOFResponse.SetT0resolution(estimatedT0resolution);
1056 }
1057
1058 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1059 if(flagT0TOF){
1060 fTOFResponse.SetT0event(startTime);
1061 fTOFResponse.SetT0resolution(startTimeRes);
1062 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1063 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1064 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1065 }
1066 }
1067 else{
1068 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1069 estimatedT0event[i]=0.0;
1070 estimatedT0resolution[i]=t0spread;
1071 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1072 }
1073 fTOFResponse.SetT0event(estimatedT0event);
1074 fTOFResponse.SetT0resolution(estimatedT0resolution);
1075 }
1076 }
1077 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1078 Float_t t0AC=-10000;
1079 Float_t t0A=-10000;
1080 Float_t t0C=-10000;
1081 if(flagT0T0){
1082 t0AC= vevent->GetT0TOF()[0];
1083 t0A= vevent->GetT0TOF()[1];
1084 t0C= vevent->GetT0TOF()[2];
1085 }
1086
1087 Float_t t0t0Best = 0;
1088 Float_t t0t0BestRes = 9999;
1089 Int_t t0used=0;
1090 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1091 t0t0Best = t0AC;
1092 t0t0BestRes = resT0AC;
1093 t0used=6;
1094 }
1095 else if(TMath::Abs(t0C) < t0cut){
1096 t0t0Best = t0C;
1097 t0t0BestRes = resT0C;
1098 t0used=4;
1099 }
1100 else if(TMath::Abs(t0A) < t0cut){
1101 t0t0Best = t0A;
1102 t0t0BestRes = resT0A;
1103 t0used=2;
1104 }
1105
1106 if(flagT0TOF){ // if T0-TOF info is available
1107 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1108 if(t0t0BestRes < 999){
1109 if(startTimeRes[i] < t0spread){
1110 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1111 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1112 estimatedT0event[i]=t0best / wtot;
1113 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1114 startTimeMask[i] = t0used+1;
1115 }
1116 else {
1117 estimatedT0event[i]=t0t0Best;
1118 estimatedT0resolution[i]=t0t0BestRes;
1119 startTimeMask[i] = t0used;
1120 }
1121 }
1122 else{
1123 estimatedT0event[i]=startTime[i];
1124 estimatedT0resolution[i]=startTimeRes[i];
1125 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1126 }
1127 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1128 }
1129 fTOFResponse.SetT0event(estimatedT0event);
1130 fTOFResponse.SetT0resolution(estimatedT0resolution);
1131 }
1132 else{ // if no T0-TOF info is available
1133 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1134 fTOFResponse.SetT0binMask(i,t0used);
1135 if(t0t0BestRes < 999){
1136 estimatedT0event[i]=t0t0Best;
1137 estimatedT0resolution[i]=t0t0BestRes;
1138 }
1139 else{
1140 estimatedT0event[i]=0.0;
1141 estimatedT0resolution[i]=t0spread;
1142 }
1143 }
1144 fTOFResponse.SetT0event(estimatedT0event);
1145 fTOFResponse.SetT0resolution(estimatedT0resolution);
1146 }
1147 }
1148
1149 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1150 Float_t t0AC=-10000;
1151 Float_t t0A=-10000;
1152 Float_t t0C=-10000;
1153 if(flagT0T0){
1154 t0AC= vevent->GetT0TOF()[0];
1155 t0A= vevent->GetT0TOF()[1];
1156 t0C= vevent->GetT0TOF()[2];
1157 }
1158
1159 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1160 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1161 estimatedT0event[i]=t0AC;
1162 estimatedT0resolution[i]=resT0AC;
1163 fTOFResponse.SetT0binMask(i,6);
1164 }
1165 }
1166 else if(TMath::Abs(t0C) < t0cut){
1167 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1168 estimatedT0event[i]=t0C;
1169 estimatedT0resolution[i]=resT0C;
1170 fTOFResponse.SetT0binMask(i,4);
1171 }
1172 }
1173 else if(TMath::Abs(t0A) < t0cut){
1174 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1175 estimatedT0event[i]=t0A;
1176 estimatedT0resolution[i]=resT0A;
1177 fTOFResponse.SetT0binMask(i,2);
1178 }
1179 }
1180 else{
1181 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1182 estimatedT0event[i]=0.0;
1183 estimatedT0resolution[i]=t0spread;
1184 fTOFResponse.SetT0binMask(i,0);
1185 }
1186 }
1187 fTOFResponse.SetT0event(estimatedT0event);
1188 fTOFResponse.SetT0resolution(estimatedT0resolution);
1189 }
1190 delete [] startTime;
1191 delete [] startTimeRes;
1192 delete [] startTimeMask;
1193 delete [] estimatedT0event;
1194 delete [] estimatedT0resolution;
1195}