git://git.uio.no / u/mrichter/AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
In AliMUONResponseTriggerV1:
[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 //
a635821f 81 AliLog::SetClassDebugLevel("AliPIDResponse",0);
82 AliLog::SetClassDebugLevel("AliESDpid",0);
83 AliLog::SetClassDebugLevel("AliAODpidUtil",0);
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
67376d1d 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 // }
fd21ec8d 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
67376d1d 431 // Double_t sig = sigma[j];
432 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
433 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
fd21ec8d 434 if (TMath::Abs(nsigmas) > (fRange+2)) {
0b39f221 435 if(nsigmas < fTOFtail)
436 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
437 else
438 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
439 } else{
440 if(nsigmas < fTOFtail)
441 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
442 else
443 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
444 }
445
446 /* OLD Gaussian shape
447 if (TMath::Abs(nsigmas) > (fRange+2)) {
fd21ec8d 448 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
449 } else
450 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
0b39f221 451 */
fd21ec8d 452
453 if (TMath::Abs(nsigmas)<5.){
454 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
455 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
456 }
457 }
458
459 if (mismatch){
460 return kDetMismatch;
461 }
462
463 // TODO: Light nuclei
464
465 return kDetPidOk;
466}
467//______________________________________________________________________________
ea235c90 468AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 469{
470 //
471 // Compute PID response for the
472 //
473
474 // set flat distribution (no decision)
475 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 476 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
477
478 Float_t mom[6];
479 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
480 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
481 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
482 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
483 mom[ilayer] = track->GetTRDmomentum(ilayer);
484 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
485 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
486 }
487 }
488 fTRDResponse.GetResponse(nslices, dedx, mom, p);
489 return kDetPidOk;
fd21ec8d 490}
491//______________________________________________________________________________
e96b9916 492AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 493{
494 //
495 // Compute PID response for the EMCAL
496 //
497
e96b9916 498 AliVCluster *matchedClus = NULL;
499
500 Double_t mom = -1.;
501 Double_t pt = -1.;
502 Double_t EovP = -1.;
503 Double_t fClsE = -1.;
504
505 Int_t nMatchClus = -1;
506 Int_t charge = 0;
507
508 // Track matching
509 nMatchClus = track->GetEMCALcluster();
510
511 if(nMatchClus > -1){
512
513 mom = track->P();
514 pt = track->Pt();
515 charge = track->Charge();
516
517 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
518
519 if(matchedClus){
520
521 // matched cluster is EMCAL
522 if(matchedClus->IsEMCAL()){
523
524 fClsE = matchedClus->E();
525 EovP = fClsE/mom;
526
527
528 // compute the probabilities
529 if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
530
531 // in case everything is OK
532 return kDetPidOk;
533
534 }
535 }
536 }
537 }
538
539 // in all other cases set flat distribution (no decision)
fd21ec8d 540 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
541 return kDetNoSignal;
e96b9916 542
fd21ec8d 543}
544//______________________________________________________________________________
545AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
546{
547 //
548 // Compute PID response for the PHOS
549 //
550
551 // set flat distribution (no decision)
552 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
553 return kDetNoSignal;
554}
555//______________________________________________________________________________
ea235c90 556AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 557{
558 //
559 // Compute PID response for the HMPID
560 //
561
562 // set flat distribution (no decision)
563 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 564 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
565
566 track->GetHMPIDpid(p);
567
568 return kDetPidOk;
fd21ec8d 569}
570
571//______________________________________________________________________________
4ec8e76d 572void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
573{
574 //
575 // Apply settings for the current event
576 //
577 fRecoPass=pass;
e96b9916 578
579 fCurrentEvent=0x0;
4ec8e76d 580 if (!event) return;
e96b9916 581 fCurrentEvent=event;
4ec8e76d 582 fRun=event->GetRunNumber();
583
584 if (fRun!=fOldRun){
585 ExecNewRun();
586 fOldRun=fRun;
587 }
588
589 //TPC resolution parametrisation PbPb
590 if ( fResolutionCorrection ){
591 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
592 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
593 }
594
595 //TOF resolution
b79db598 596 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
597
80f28562 598
599 // Get and set centrality
600 AliCentrality *centrality = event->GetCentrality();
601 if(centrality){
602 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
603 }
604 else{
605 fCurrCentrality = -1;
606 }
4ec8e76d 607}
608
609//______________________________________________________________________________
610void AliPIDResponse::ExecNewRun()
611{
612 //
613 // Things to Execute upon a new run
614 //
615 SetRecoInfo();
616
617 SetITSParametrisation();
618
619 SetTPCPidResponseMaster();
620 SetTPCParametrisation();
53d016dc 621
622 SetTRDPidResponseMaster();
623 InitializeTRDResponse();
b2138b40 624
625 SetEMCALPidResponseMaster();
626 InitializeEMCALResponse();
4ec8e76d 627
b79db598 628 SetTOFPidResponseMaster();
629 InitializeTOFResponse();
4ec8e76d 630}
631
632//_____________________________________________________
633Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
634{
635 //
636 // Get TPC multiplicity in bins of 150
637 //
638
639 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
640 Double_t tpcMulti=0.;
641 if(vertexTPC){
642 Double_t vertexContribTPC=vertexTPC->GetNContributors();
643 tpcMulti=vertexContribTPC/150.;
644 if (tpcMulti>20.) tpcMulti=20.;
645 }
646
647 return tpcMulti;
648}
649
650//______________________________________________________________________________
651void AliPIDResponse::SetRecoInfo()
652{
653 //
654 // Set reconstruction information
655 //
656
657 //reset information
658 fLHCperiod="";
659 fMCperiodTPC="";
660
661 fBeamType="";
662
663 fBeamType="PP";
664
1436d6bb 665
666 TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
667
4ec8e76d 668 //find the period by run number (UGLY, but not stored in ESD and AOD... )
669 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
670 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
671 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
99e9d5ec 672 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 673 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
674 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
ea235c90 675 else if (fRun>=136851&&fRun<=139517) {
676 fLHCperiod="LHC10H";
677 fMCperiodTPC="LHC10H8";
678 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
679 fBeamType="PBPB";
680 }
3077a03d 681 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
2e97a211 682 //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
683 else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
684 else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
3077a03d 685 else if (fRun>=166529) {
686 fLHCperiod="LHC11H";
687 fMCperiodTPC="LHC11A10";
688 fBeamType="PBPB";
689 }
690
80ab5635 691
692 //exception new pp MC productions from 2011
693 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
4a527e08 694 // exception for 11f1
695 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
4ec8e76d 696}
697
698//______________________________________________________________________________
699void AliPIDResponse::SetITSParametrisation()
700{
701 //
702 // Set the ITS parametrisation
703 //
704}
705
706//______________________________________________________________________________
707void AliPIDResponse::SetTPCPidResponseMaster()
708{
709 //
710 // Load the TPC pid response functions from the OADB
711 //
09b50a42 712 //don't load twice for the moment
713 if (fArrPidResponseMaster) return;
714
715
4ec8e76d 716 //reset the PID response functions
717 delete fArrPidResponseMaster;
718 fArrPidResponseMaster=0x0;
719
720 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
721
ea235c90 722 TFile *f=TFile::Open(fileName.Data());
723 if (f && f->IsOpen() && !f->IsZombie()){
724 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 725 }
ea235c90 726 delete f;
4ec8e76d 727
728 if (!fArrPidResponseMaster){
729 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
730 return;
731 }
732 fArrPidResponseMaster->SetOwner();
733}
734
735//______________________________________________________________________________
736void AliPIDResponse::SetTPCParametrisation()
737{
738 //
739 // Change BB parametrisation for current run
740 //
741
742 if (fLHCperiod.IsNull()) {
743 AliFatal("No period set, not changing parametrisation");
744 return;
745 }
746
747 //
748 // Set default parametrisations for data and MC
749 //
750
751 //data type
752 TString datatype="DATA";
753 //in case of mc fRecoPass is per default 1
754 if (fIsMC) {
755 datatype="MC";
756 fRecoPass=1;
757 }
758
759 //
760 //reset old splines
761 //
762 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
763 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
764 }
4a527e08 765
766 // period
767 TString period=fLHCperiod;
768 if (fIsMC) period=fMCperiodTPC;
769
770 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
771 Bool_t found=kFALSE;
4ec8e76d 772 //
773 //set the new PID splines
774 //
4ec8e76d 775 if (fArrPidResponseMaster){
776 TObject *grAll=0x0;
777 //for MC don't use period information
778// if (fIsMC) period="[A-Z0-9]*";
779 //for MC use MC period information
4ec8e76d 780//pattern for the default entry (valid for all particles)
781 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
782
783 //loop over entries and filter them
784 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
785 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
786 if (responseFunction==0x0) continue;
787 TString responseName=responseFunction->GetName();
788
789 if (!reg.MatchB(responseName)) continue;
790
791 TObjArray *arr=reg.MatchS(responseName);
792 TString particleName=arr->At(1)->GetName();
793 delete arr;
794 if (particleName.IsNull()) continue;
795 if (particleName=="ALL") grAll=responseFunction;
796 else {
797 //find particle id
798 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
799 TString particle=AliPID::ParticleName(ispec);
800 particle.ToUpper();
801 if ( particle == particleName ){
802 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
803 fTPCResponse.SetUseDatabase(kTRUE);
804 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
4a527e08 805 found=kTRUE;
4ec8e76d 806 break;
807 }
808 }
809 }
810 }
811
812 //set default response function to all particles which don't have a specific one
813 if (grAll){
814 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
815 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
816 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
817 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
4a527e08 818 found=kTRUE;
4ec8e76d 819 }
820 }
821 }
822 }
4a527e08 823
824 if (!found){
825 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
826 }
4ec8e76d 827 //
828 // Setup resolution parametrisation
829 //
830
831 //default
832 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
833
834 if (fRun>=122195){
835 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
836 }
23425eb2 837 if (fArrPidResponseMaster)
4ec8e76d 838 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
839
840 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
841}
842
ea235c90 843//______________________________________________________________________________
844void AliPIDResponse::SetTRDPidResponseMaster()
845{
846 //
847 // Load the TRD pid params and references from the OADB
848 //
849 if(fTRDPIDParams) return;
53d016dc 850 AliOADBContainer contParams("contParams");
851
59a8e853 852 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
853 if(statusPars){
854 AliError("Failed initializing PID Params from OADB");
855 } else {
856 AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
b2138b40 857 fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
59a8e853 858 if(!fTRDPIDParams){
859 AliError(Form("TRD Params not found in run %d", fRun));
860 }
861 }
ea235c90 862
53d016dc 863 AliOADBContainer contRefs("contRefs");
59a8e853 864 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
865 if(statusRefs){
866 AliInfo("Failed Loading References for TRD");
867 } else {
868 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
869 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
870 if(!fTRDPIDReference){
871 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
872 }
873 }
ea235c90 874}
875
876//______________________________________________________________________________
877void AliPIDResponse::InitializeTRDResponse(){
878 //
879 // Set PID Params and references to the TRD PID response
880 //
881 fTRDResponse.SetPIDParams(fTRDPIDParams);
882 fTRDResponse.Load(fTRDPIDReference);
883 if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
884 fTRDslicesForPID[0] = 0;
885 fTRDslicesForPID[1] = 7;
886 }
887}
888
b79db598 889//______________________________________________________________________________
890void AliPIDResponse::SetTOFPidResponseMaster()
891{
892 //
893 // Load the TOF pid params from the OADB
894 //
895 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
896 if (oadbf->IsOpen()) {
897 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
898 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
899 if (fTOFPIDParams) delete fTOFPIDParams;
f57f1e18 900 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 901 oadbf->Close();
902 delete oadbc;
903 } else {
904 AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
905 }
906 delete oadbf;
907
908 }
909
910//______________________________________________________________________________
911void AliPIDResponse::InitializeTOFResponse(){
912 //
913 // Set PID Params to the TOF PID response
914 //
b79db598 915 for (Int_t i=0;i<4;i++) {
916 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
917 }
918 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
919
920}
921
922
ea235c90 923//_________________________________________________________________________
924Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
925 //
926 // Check whether track is identified as electron under a given electron efficiency hypothesis
927 //
928 Double_t probs[AliPID::kSPECIES];
929 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
930
99e9d5ec 931 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
932 // Take mean of the TRD momenta in the given tracklets
933 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
934 Int_t nmomenta = 0;
ea235c90 935 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
936 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 937 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 938 }
939 }
99e9d5ec 940 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 941
942 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
943}
944
b2138b40 945//______________________________________________________________________________
946void AliPIDResponse::SetEMCALPidResponseMaster()
947{
948 //
949 // Load the EMCAL pid response functions from the OADB
950 //
951 TObjArray* fEMCALPIDParamsRun = NULL;
952 TObjArray* fEMCALPIDParamsPass = NULL;
953
954 if(fEMCALPIDParams) return;
955 AliOADBContainer contParams("contParams");
956
957 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
958 if(statusPars){
959 AliError("Failed initializing PID Params from OADB");
960 }
961 else {
962 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
963
964 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
965 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
966 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
967
968 if(!fEMCALPIDParams){
f8d39067 969 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 970 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 971
1f631618 972 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
973 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 974 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
975
976 if(!fEMCALPIDParams){
1f631618 977 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 978 }
979 }
980 }
981}
982
983//______________________________________________________________________________
984void AliPIDResponse::InitializeEMCALResponse(){
985 //
986 // Set PID Params to the EMCAL PID response
987 //
988 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
989
990}
5f8db5fe 991//_________________________________________________________________________
992void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
993 //
994 // Set TOF response function
995 // Input option for event_time used
996 //
997
998 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
999 if(t0spread < 10) t0spread = 80;
1000
1001 // T0 from TOF algorithm
1002
1003 Bool_t flagT0TOF=kFALSE;
1004 Bool_t flagT0T0=kFALSE;
1005 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1006 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1007 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1008
1009 // T0-TOF arrays
1010 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1011 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1012 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1013 estimatedT0event[i]=0.0;
1014 estimatedT0resolution[i]=0.0;
1015 startTimeMask[i] = 0;
1016 }
1017
1018 Float_t resT0A=75,resT0C=65,resT0AC=55;
1019 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1020 flagT0T0=kTRUE;
1021 }
1022
1023
1024 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1025
1026 if (tofHeader) { // read global info and T0-TOF
1027 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1028 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1029 if(t0spread < 10) t0spread = 80;
1030
1031 flagT0TOF=kTRUE;
1032 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1033 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1034 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1035 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1036 }
1037
1038 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1039 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1040 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1041 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1042 Int_t icurrent = (Int_t)ibin->GetAt(j);
1043 startTime[icurrent]=t0Bin->GetAt(j);
1044 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1045 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1046 }
1047 }
1048
1049 // for cut of 3 sigma on t0 spread
1050 Float_t t0cut = 3 * t0spread;
1051 if(t0cut < 500) t0cut = 500;
1052
1053 if(option == kFILL_T0){ // T0-FILL is used
1054 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1055 estimatedT0event[i]=0.0;
1056 estimatedT0resolution[i]=t0spread;
1057 }
1058 fTOFResponse.SetT0event(estimatedT0event);
1059 fTOFResponse.SetT0resolution(estimatedT0resolution);
1060 }
1061
1062 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1063 if(flagT0TOF){
1064 fTOFResponse.SetT0event(startTime);
1065 fTOFResponse.SetT0resolution(startTimeRes);
1066 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1067 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1068 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1069 }
1070 }
1071 else{
1072 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1073 estimatedT0event[i]=0.0;
1074 estimatedT0resolution[i]=t0spread;
1075 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1076 }
1077 fTOFResponse.SetT0event(estimatedT0event);
1078 fTOFResponse.SetT0resolution(estimatedT0resolution);
1079 }
1080 }
1081 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1082 Float_t t0AC=-10000;
1083 Float_t t0A=-10000;
1084 Float_t t0C=-10000;
1085 if(flagT0T0){
1086 t0AC= vevent->GetT0TOF()[0];
1087 t0A= vevent->GetT0TOF()[1];
1088 t0C= vevent->GetT0TOF()[2];
1089 }
1090
1091 Float_t t0t0Best = 0;
1092 Float_t t0t0BestRes = 9999;
1093 Int_t t0used=0;
1094 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1095 t0t0Best = t0AC;
1096 t0t0BestRes = resT0AC;
1097 t0used=6;
1098 }
1099 else if(TMath::Abs(t0C) < t0cut){
1100 t0t0Best = t0C;
1101 t0t0BestRes = resT0C;
1102 t0used=4;
1103 }
1104 else if(TMath::Abs(t0A) < t0cut){
1105 t0t0Best = t0A;
1106 t0t0BestRes = resT0A;
1107 t0used=2;
1108 }
1109
1110 if(flagT0TOF){ // if T0-TOF info is available
1111 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1112 if(t0t0BestRes < 999){
1113 if(startTimeRes[i] < t0spread){
1114 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1115 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1116 estimatedT0event[i]=t0best / wtot;
1117 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1118 startTimeMask[i] = t0used+1;
1119 }
1120 else {
1121 estimatedT0event[i]=t0t0Best;
1122 estimatedT0resolution[i]=t0t0BestRes;
1123 startTimeMask[i] = t0used;
1124 }
1125 }
1126 else{
1127 estimatedT0event[i]=startTime[i];
1128 estimatedT0resolution[i]=startTimeRes[i];
1129 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1130 }
1131 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1132 }
1133 fTOFResponse.SetT0event(estimatedT0event);
1134 fTOFResponse.SetT0resolution(estimatedT0resolution);
1135 }
1136 else{ // if no T0-TOF info is available
1137 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1138 fTOFResponse.SetT0binMask(i,t0used);
1139 if(t0t0BestRes < 999){
1140 estimatedT0event[i]=t0t0Best;
1141 estimatedT0resolution[i]=t0t0BestRes;
1142 }
1143 else{
1144 estimatedT0event[i]=0.0;
1145 estimatedT0resolution[i]=t0spread;
1146 }
1147 }
1148 fTOFResponse.SetT0event(estimatedT0event);
1149 fTOFResponse.SetT0resolution(estimatedT0resolution);
1150 }
1151 }
1152
1153 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1154 Float_t t0AC=-10000;
1155 Float_t t0A=-10000;
1156 Float_t t0C=-10000;
1157 if(flagT0T0){
1158 t0AC= vevent->GetT0TOF()[0];
1159 t0A= vevent->GetT0TOF()[1];
1160 t0C= vevent->GetT0TOF()[2];
1161 }
1162
1163 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1164 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1165 estimatedT0event[i]=t0AC;
1166 estimatedT0resolution[i]=resT0AC;
1167 fTOFResponse.SetT0binMask(i,6);
1168 }
1169 }
1170 else if(TMath::Abs(t0C) < t0cut){
1171 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1172 estimatedT0event[i]=t0C;
1173 estimatedT0resolution[i]=resT0C;
1174 fTOFResponse.SetT0binMask(i,4);
1175 }
1176 }
1177 else if(TMath::Abs(t0A) < t0cut){
1178 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1179 estimatedT0event[i]=t0A;
1180 estimatedT0resolution[i]=resT0A;
1181 fTOFResponse.SetT0binMask(i,2);
1182 }
1183 }
1184 else{
1185 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1186 estimatedT0event[i]=0.0;
1187 estimatedT0resolution[i]=t0spread;
1188 fTOFResponse.SetT0binMask(i,0);
1189 }
1190 }
1191 fTOFResponse.SetT0event(estimatedT0event);
1192 fTOFResponse.SetT0resolution(estimatedT0resolution);
1193 }
1194 delete [] startTime;
1195 delete [] startTimeRes;
1196 delete [] startTimeMask;
1197 delete [] estimatedT0event;
1198 delete [] estimatedT0resolution;
1199}