load properly the library in case of par files, update macro with possibility to...
[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";
4ec8e76d 690}
691
692//______________________________________________________________________________
693void AliPIDResponse::SetITSParametrisation()
694{
695 //
696 // Set the ITS parametrisation
697 //
698}
699
700//______________________________________________________________________________
701void AliPIDResponse::SetTPCPidResponseMaster()
702{
703 //
704 // Load the TPC pid response functions from the OADB
705 //
09b50a42 706 //don't load twice for the moment
707 if (fArrPidResponseMaster) return;
708
709
4ec8e76d 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
ea235c90 716 TFile *f=TFile::Open(fileName.Data());
717 if (f && f->IsOpen() && !f->IsZombie()){
718 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 719 }
ea235c90 720 delete f;
4ec8e76d 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//______________________________________________________________________________
730void 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 }
23425eb2 822 if (fArrPidResponseMaster)
4ec8e76d 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
ea235c90 828//______________________________________________________________________________
829void AliPIDResponse::SetTRDPidResponseMaster()
830{
831 //
832 // Load the TRD pid params and references from the OADB
833 //
834 if(fTRDPIDParams) return;
53d016dc 835 AliOADBContainer contParams("contParams");
836
59a8e853 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()));
b2138b40 842 fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
59a8e853 843 if(!fTRDPIDParams){
844 AliError(Form("TRD Params not found in run %d", fRun));
845 }
846 }
ea235c90 847
53d016dc 848 AliOADBContainer contRefs("contRefs");
59a8e853 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 }
ea235c90 859}
860
861//______________________________________________________________________________
862void 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
b79db598 874//______________________________________________________________________________
875void 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;
f57f1e18 885 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 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//______________________________________________________________________________
896void AliPIDResponse::InitializeTOFResponse(){
897 //
898 // Set PID Params to the TOF PID response
899 //
b79db598 900 for (Int_t i=0;i<4;i++) {
901 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
902 }
903 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
904
905}
906
907
ea235c90 908//_________________________________________________________________________
909Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
910 //
911 // Check whether track is identified as electron under a given electron efficiency hypothesis
912 //
913 Double_t probs[AliPID::kSPECIES];
914 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
915
99e9d5ec 916 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
917 // Take mean of the TRD momenta in the given tracklets
918 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
919 Int_t nmomenta = 0;
ea235c90 920 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
921 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 922 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 923 }
924 }
99e9d5ec 925 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 926
927 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
928}
929
b2138b40 930//______________________________________________________________________________
931void AliPIDResponse::SetEMCALPidResponseMaster()
932{
933 //
934 // Load the EMCAL pid response functions from the OADB
935 //
936 TObjArray* fEMCALPIDParamsRun = NULL;
937 TObjArray* fEMCALPIDParamsPass = NULL;
938
939 if(fEMCALPIDParams) return;
940 AliOADBContainer contParams("contParams");
941
942 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
943 if(statusPars){
944 AliError("Failed initializing PID Params from OADB");
945 }
946 else {
947 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
948
949 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
950 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
951 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
952
953 if(!fEMCALPIDParams){
f8d39067 954 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 955 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 956
1f631618 957 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
958 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 959 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
960
961 if(!fEMCALPIDParams){
1f631618 962 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 963 }
964 }
965 }
966}
967
968//______________________________________________________________________________
969void AliPIDResponse::InitializeEMCALResponse(){
970 //
971 // Set PID Params to the EMCAL PID response
972 //
973 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
974
975}
5f8db5fe 976//_________________________________________________________________________
977void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
978 //
979 // Set TOF response function
980 // Input option for event_time used
981 //
982
983 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
984 if(t0spread < 10) t0spread = 80;
985
986 // T0 from TOF algorithm
987
988 Bool_t flagT0TOF=kFALSE;
989 Bool_t flagT0T0=kFALSE;
990 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
991 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
992 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
993
994 // T0-TOF arrays
995 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
996 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
997 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
998 estimatedT0event[i]=0.0;
999 estimatedT0resolution[i]=0.0;
1000 startTimeMask[i] = 0;
1001 }
1002
1003 Float_t resT0A=75,resT0C=65,resT0AC=55;
1004 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1005 flagT0T0=kTRUE;
1006 }
1007
1008
1009 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1010
1011 if (tofHeader) { // read global info and T0-TOF
1012 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1013 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1014 if(t0spread < 10) t0spread = 80;
1015
1016 flagT0TOF=kTRUE;
1017 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1018 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1019 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1020 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1021 }
1022
1023 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1024 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1025 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1026 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1027 Int_t icurrent = (Int_t)ibin->GetAt(j);
1028 startTime[icurrent]=t0Bin->GetAt(j);
1029 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1030 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1031 }
1032 }
1033
1034 // for cut of 3 sigma on t0 spread
1035 Float_t t0cut = 3 * t0spread;
1036 if(t0cut < 500) t0cut = 500;
1037
1038 if(option == kFILL_T0){ // T0-FILL is used
1039 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1040 estimatedT0event[i]=0.0;
1041 estimatedT0resolution[i]=t0spread;
1042 }
1043 fTOFResponse.SetT0event(estimatedT0event);
1044 fTOFResponse.SetT0resolution(estimatedT0resolution);
1045 }
1046
1047 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1048 if(flagT0TOF){
1049 fTOFResponse.SetT0event(startTime);
1050 fTOFResponse.SetT0resolution(startTimeRes);
1051 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1052 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1053 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1054 }
1055 }
1056 else{
1057 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1058 estimatedT0event[i]=0.0;
1059 estimatedT0resolution[i]=t0spread;
1060 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1061 }
1062 fTOFResponse.SetT0event(estimatedT0event);
1063 fTOFResponse.SetT0resolution(estimatedT0resolution);
1064 }
1065 }
1066 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1067 Float_t t0AC=-10000;
1068 Float_t t0A=-10000;
1069 Float_t t0C=-10000;
1070 if(flagT0T0){
1071 t0AC= vevent->GetT0TOF()[0];
1072 t0A= vevent->GetT0TOF()[1];
1073 t0C= vevent->GetT0TOF()[2];
1074 }
1075
1076 Float_t t0t0Best = 0;
1077 Float_t t0t0BestRes = 9999;
1078 Int_t t0used=0;
1079 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1080 t0t0Best = t0AC;
1081 t0t0BestRes = resT0AC;
1082 t0used=6;
1083 }
1084 else if(TMath::Abs(t0C) < t0cut){
1085 t0t0Best = t0C;
1086 t0t0BestRes = resT0C;
1087 t0used=4;
1088 }
1089 else if(TMath::Abs(t0A) < t0cut){
1090 t0t0Best = t0A;
1091 t0t0BestRes = resT0A;
1092 t0used=2;
1093 }
1094
1095 if(flagT0TOF){ // if T0-TOF info is available
1096 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1097 if(t0t0BestRes < 999){
1098 if(startTimeRes[i] < t0spread){
1099 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1100 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1101 estimatedT0event[i]=t0best / wtot;
1102 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1103 startTimeMask[i] = t0used+1;
1104 }
1105 else {
1106 estimatedT0event[i]=t0t0Best;
1107 estimatedT0resolution[i]=t0t0BestRes;
1108 startTimeMask[i] = t0used;
1109 }
1110 }
1111 else{
1112 estimatedT0event[i]=startTime[i];
1113 estimatedT0resolution[i]=startTimeRes[i];
1114 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1115 }
1116 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1117 }
1118 fTOFResponse.SetT0event(estimatedT0event);
1119 fTOFResponse.SetT0resolution(estimatedT0resolution);
1120 }
1121 else{ // if no T0-TOF info is available
1122 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1123 fTOFResponse.SetT0binMask(i,t0used);
1124 if(t0t0BestRes < 999){
1125 estimatedT0event[i]=t0t0Best;
1126 estimatedT0resolution[i]=t0t0BestRes;
1127 }
1128 else{
1129 estimatedT0event[i]=0.0;
1130 estimatedT0resolution[i]=t0spread;
1131 }
1132 }
1133 fTOFResponse.SetT0event(estimatedT0event);
1134 fTOFResponse.SetT0resolution(estimatedT0resolution);
1135 }
1136 }
1137
1138 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1139 Float_t t0AC=-10000;
1140 Float_t t0A=-10000;
1141 Float_t t0C=-10000;
1142 if(flagT0T0){
1143 t0AC= vevent->GetT0TOF()[0];
1144 t0A= vevent->GetT0TOF()[1];
1145 t0C= vevent->GetT0TOF()[2];
1146 }
1147
1148 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1149 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1150 estimatedT0event[i]=t0AC;
1151 estimatedT0resolution[i]=resT0AC;
1152 fTOFResponse.SetT0binMask(i,6);
1153 }
1154 }
1155 else if(TMath::Abs(t0C) < t0cut){
1156 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1157 estimatedT0event[i]=t0C;
1158 estimatedT0resolution[i]=resT0C;
1159 fTOFResponse.SetT0binMask(i,4);
1160 }
1161 }
1162 else if(TMath::Abs(t0A) < t0cut){
1163 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1164 estimatedT0event[i]=t0A;
1165 estimatedT0resolution[i]=resT0A;
1166 fTOFResponse.SetT0binMask(i,2);
1167 }
1168 }
1169 else{
1170 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1171 estimatedT0event[i]=0.0;
1172 estimatedT0resolution[i]=t0spread;
1173 fTOFResponse.SetT0binMask(i,0);
1174 }
1175 }
1176 fTOFResponse.SetT0event(estimatedT0event);
1177 fTOFResponse.SetT0resolution(estimatedT0resolution);
1178 }
1179 delete [] startTime;
1180 delete [] startTimeRes;
1181 delete [] startTimeMask;
1182 delete [] estimatedT0event;
1183 delete [] estimatedT0resolution;
1184}