git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates in order to enable the '2D' PID for the TRD developed by Daniel Lohner.
[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>
db0e2c5f 32#include <TArrayF.h>
4ec8e76d 33
34#include <AliVEvent.h>
fd21ec8d 35#include <AliVTrack.h>
4ec8e76d 36#include <AliLog.h>
37#include <AliPID.h>
ea235c90 38#include <AliOADBContainer.h>
db0e2c5f 39#include <AliTRDPIDResponseObject.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),
db0e2c5f 70fTRDPIDResponseObject(0x0),
0b39f221 71fTOFtail(1.1),
b79db598 72fTOFPIDParams(0x0),
b2138b40 73fEMCALPIDParams(0x0),
80f28562 74fCurrentEvent(0x0),
75fCurrCentrality(0.0)
4ec8e76d 76{
77 //
78 // default ctor
79 //
a635821f 80 AliLog::SetClassDebugLevel("AliPIDResponse",0);
81 AliLog::SetClassDebugLevel("AliESDpid",0);
82 AliLog::SetClassDebugLevel("AliAODpidUtil",0);
ea235c90 83
84 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 85}
86
87//______________________________________________________________________________
88AliPIDResponse::~AliPIDResponse()
89{
90 //
91 // dtor
92 //
db0e2c5f 93 delete fArrPidResponseMaster;
94 delete fTRDPIDResponseObject;
b79db598 95 if (fTOFPIDParams) delete fTOFPIDParams;
4ec8e76d 96}
97
98//______________________________________________________________________________
99AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
100TNamed(other),
101fITSResponse(other.fITSResponse),
102fTPCResponse(other.fTPCResponse),
103fTRDResponse(other.fTRDResponse),
104fTOFResponse(other.fTOFResponse),
e96b9916 105fEMCALResponse(other.fEMCALResponse),
fd21ec8d 106fRange(other.fRange),
107fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 108fIsMC(other.fIsMC),
109fOADBPath(other.fOADBPath),
110fBeamType("PP"),
111fLHCperiod(),
112fMCperiodTPC(),
fd21ec8d 113fMCperiodUser(other.fMCperiodUser),
ea235c90 114fCurrentFile(),
4ec8e76d 115fRecoPass(0),
fd21ec8d 116fRecoPassUser(other.fRecoPassUser),
4ec8e76d 117fRun(0),
118fOldRun(0),
119fArrPidResponseMaster(0x0),
120fResolutionCorrection(0x0),
db0e2c5f 121fTRDPIDResponseObject(0x0),
0b39f221 122fTOFtail(1.1),
b79db598 123fTOFPIDParams(0x0),
b2138b40 124fEMCALPIDParams(0x0),
80f28562 125fCurrentEvent(0x0),
126fCurrCentrality(0.0)
4ec8e76d 127{
128 //
129 // copy ctor
130 //
ea235c90 131 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 132}
133
134//______________________________________________________________________________
135AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
136{
137 //
138 // copy ctor
139 //
140 if(this!=&other) {
141 delete fArrPidResponseMaster;
142 TNamed::operator=(other);
143 fITSResponse=other.fITSResponse;
144 fTPCResponse=other.fTPCResponse;
145 fTRDResponse=other.fTRDResponse;
146 fTOFResponse=other.fTOFResponse;
e96b9916 147 fEMCALResponse=other.fEMCALResponse;
fd21ec8d 148 fRange=other.fRange;
149 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 150 fOADBPath=other.fOADBPath;
151 fIsMC=other.fIsMC;
152 fBeamType="PP";
153 fLHCperiod="";
154 fMCperiodTPC="";
fd21ec8d 155 fMCperiodUser=other.fMCperiodUser;
ea235c90 156 fCurrentFile="";
4ec8e76d 157 fRecoPass=0;
fd21ec8d 158 fRecoPassUser=other.fRecoPassUser;
4ec8e76d 159 fRun=0;
160 fOldRun=0;
161 fArrPidResponseMaster=0x0;
162 fResolutionCorrection=0x0;
db0e2c5f 163 fTRDPIDResponseObject=0x0;
b2138b40 164 fEMCALPIDParams=0x0;
ea235c90 165 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
0b39f221 166 fTOFtail=1.1;
b79db598 167 fTOFPIDParams=0x0;
e96b9916 168 fCurrentEvent=other.fCurrentEvent;
4ec8e76d 169 }
170 return *this;
171}
172
fd21ec8d 173//______________________________________________________________________________
174Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
175{
176 //
177 // NumberOfSigmas for 'detCode'
178 //
179
180 switch (detCode){
181 case kDetITS: return NumberOfSigmasITS(track, type); break;
182 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
183 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
184// case kDetTRD: return ComputeTRDProbability(track, type); break;
185// case kDetPHOS: return ComputePHOSProbability(track, type); break;
e96b9916 186// case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
fd21ec8d 187// case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
188 default: return -999.;
189 }
190
191}
192
e96b9916 193//______________________________________________________________________________
194Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
195
196 AliVCluster *matchedClus = NULL;
197
198 Double_t mom = -1.;
199 Double_t pt = -1.;
200 Double_t EovP = -1.;
201 Double_t fClsE = -1.;
202
203 Int_t nMatchClus = -1;
204 Int_t charge = 0;
205
206 // Track matching
207 nMatchClus = track->GetEMCALcluster();
208 if(nMatchClus > -1){
6d0064aa 209
e96b9916 210 mom = track->P();
211 pt = track->Pt();
212 charge = track->Charge();
213
214 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
215
216 if(matchedClus){
217
6d0064aa 218 // matched cluster is EMCAL
219 if(matchedClus->IsEMCAL()){
220
221 fClsE = matchedClus->E();
222 EovP = fClsE/mom;
223
224
225 // NSigma value really meaningful only for electrons!
226 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
227 }
e96b9916 228 }
229 }
6d0064aa 230
e96b9916 231 return -999;
232
233}
234
6d0064aa 235//______________________________________________________________________________
236Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
237
238 AliVCluster *matchedClus = NULL;
239
240 Double_t mom = -1.;
241 Double_t pt = -1.;
242 Double_t EovP = -1.;
243 Double_t fClsE = -1.;
244
245 Int_t nMatchClus = -1;
246 Int_t charge = 0;
247
248 // Track matching
249 nMatchClus = track->GetEMCALcluster();
250 if(nMatchClus > -1){
251
252 mom = track->P();
253 pt = track->Pt();
254 charge = track->Charge();
255
256 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
257
258 if(matchedClus){
259
260 // matched cluster is EMCAL
261 if(matchedClus->IsEMCAL()){
262
263 fClsE = matchedClus->E();
264 EovP = fClsE/mom;
265
266 // fill used EMCAL variables here
267 eop = EovP; // E/p
268 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
269 showershape[1] = matchedClus->GetM02(); // long axis
270 showershape[2] = matchedClus->GetM20(); // short axis
271 showershape[3] = matchedClus->GetDispersion(); // dispersion
272
273 // NSigma value really meaningful only for electrons!
274 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
275 }
276 }
277 }
278 return -999;
279}
280
281
fd21ec8d 282//______________________________________________________________________________
283AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
284{
285 //
286 // Compute PID response of 'detCode'
287 //
288
289 switch (detCode){
290 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
291 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
292 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
293 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
294 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
295 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
296 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
297 default: return kDetNoSignal;
298 }
299}
300
301//______________________________________________________________________________
302AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
303{
304 //
305 // Compute PID response for the ITS
306 //
307
308 // set flat distribution (no decision)
309 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
310
311 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
312 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
313
314 Double_t mom=track->P();
315 Double_t dedx=track->GetITSsignal();
316 Bool_t isSA=kTRUE;
317 Double_t momITS=mom;
318 ULong_t trStatus=track->GetStatus();
319 if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
320 UChar_t clumap=track->GetITSClusterMap();
321 Int_t nPointsForPid=0;
322 for(Int_t i=2; i<6; i++){
323 if(clumap&(1<<i)) ++nPointsForPid;
324 }
325
326 if(nPointsForPid<3) { // track not to be used for combined PID purposes
327 // track->ResetStatus(AliVTrack::kITSpid);
328 return kDetNoSignal;
329 }
330
331 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
332 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
333 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
334 Double_t bethe=fITSResponse.Bethe(momITS,mass);
335 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
336 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
337 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
338 } else {
339 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
340 mismatch=kFALSE;
341 }
342
343 // Check for particles heavier than (AliPID::kSPECIES - 1)
344 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
345
346 }
347
348 if (mismatch){
349 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
350 return kDetNoSignal;
351 }
352
353
354 return kDetPidOk;
355}
356//______________________________________________________________________________
357AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
358{
359 //
360 // Compute PID response for the TPC
361 //
362
363 // set flat distribution (no decision)
364 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
365
366 // check quality of the track
367 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
368
369 Double_t mom = track->GetTPCmomentum();
370
371 Double_t dedx=track->GetTPCsignal();
372 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
373
374 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
375 AliPID::EParticleType type=AliPID::EParticleType(j);
376 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
377 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
378 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
379 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
380 } else {
381 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
382 mismatch=kFALSE;
383 }
384
385 // TODO: Light nuclei, also in TPC pid response
386
387 // Check for particles heavier than (AliPID::kSPECIES - 1)
388// if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
389
390 }
391
392 if (mismatch){
393 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
394 return kDetNoSignal;
395 }
396
397 return kDetPidOk;
398}
399//______________________________________________________________________________
400AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
401{
402 //
403 // Compute PID response for the
404 //
405
0b39f221 406 Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
407
fd21ec8d 408 // set flat distribution (no decision)
409 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
410
411 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
412 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
413
414 Double_t time[AliPID::kSPECIESN];
415 track->GetIntegratedTimes(time);
416
67376d1d 417 // Double_t sigma[AliPID::kSPECIES];
418 // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
419 // sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
420 // }
fd21ec8d 421
422 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
423 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
424 AliPID::EParticleType type=AliPID::EParticleType(j);
0b39f221 425 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
426
67376d1d 427 // Double_t sig = sigma[j];
428 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
429 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
0b39f221 430 if (TMath::Abs(nsigmas) > (fRange+2)) {
431 if(nsigmas < fTOFtail)
432 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
433 else
434 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
435 } else{
436 if(nsigmas < fTOFtail)
437 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
438 else
439 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
440 }
441
442 /* OLD Gaussian shape
fd21ec8d 443 if (TMath::Abs(nsigmas) > (fRange+2)) {
444 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
445 } else
446 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
0b39f221 447 */
fd21ec8d 448
449 if (TMath::Abs(nsigmas)<5.){
450 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
451 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
452 }
453 }
454
455 if (mismatch){
456 return kDetMismatch;
457 }
458
459 // TODO: Light nuclei
460
461 return kDetPidOk;
462}
463//______________________________________________________________________________
ea235c90 464AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 465{
466 //
467 // Compute PID response for the
468 //
469
470 // set flat distribution (no decision)
471 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 472 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
473
474 Float_t mom[6];
475 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
476 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
477 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
478 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
479 mom[ilayer] = track->GetTRDmomentum(ilayer);
480 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
481 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
482 }
483 }
484 fTRDResponse.GetResponse(nslices, dedx, mom, p);
485 return kDetPidOk;
fd21ec8d 486}
487//______________________________________________________________________________
e96b9916 488AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 489{
490 //
491 // Compute PID response for the EMCAL
492 //
493
e96b9916 494 AliVCluster *matchedClus = NULL;
495
496 Double_t mom = -1.;
497 Double_t pt = -1.;
498 Double_t EovP = -1.;
499 Double_t fClsE = -1.;
500
501 Int_t nMatchClus = -1;
502 Int_t charge = 0;
503
504 // Track matching
505 nMatchClus = track->GetEMCALcluster();
506
507 if(nMatchClus > -1){
508
509 mom = track->P();
510 pt = track->Pt();
511 charge = track->Charge();
512
513 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
514
515 if(matchedClus){
516
517 // matched cluster is EMCAL
518 if(matchedClus->IsEMCAL()){
519
520 fClsE = matchedClus->E();
521 EovP = fClsE/mom;
522
523
524 // compute the probabilities
525 if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
526
527 // in case everything is OK
528 return kDetPidOk;
529
530 }
531 }
532 }
533 }
534
535 // in all other cases set flat distribution (no decision)
fd21ec8d 536 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
537 return kDetNoSignal;
e96b9916 538
fd21ec8d 539}
540//______________________________________________________________________________
541AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
542{
543 //
544 // Compute PID response for the PHOS
545 //
546
547 // set flat distribution (no decision)
548 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
549 return kDetNoSignal;
550}
551//______________________________________________________________________________
ea235c90 552AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 553{
554 //
555 // Compute PID response for the HMPID
556 //
557
558 // set flat distribution (no decision)
559 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 560 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
561
562 track->GetHMPIDpid(p);
563
564 return kDetPidOk;
fd21ec8d 565}
566
4ec8e76d 567//______________________________________________________________________________
568void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
569{
570 //
571 // Apply settings for the current event
572 //
573 fRecoPass=pass;
e96b9916 574
575 fCurrentEvent=0x0;
4ec8e76d 576 if (!event) return;
e96b9916 577 fCurrentEvent=event;
4ec8e76d 578 fRun=event->GetRunNumber();
579
580 if (fRun!=fOldRun){
581 ExecNewRun();
582 fOldRun=fRun;
583 }
584
585 //TPC resolution parametrisation PbPb
586 if ( fResolutionCorrection ){
587 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
588 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
589 }
590
591 //TOF resolution
b79db598 592 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
593
80f28562 594
595 // Get and set centrality
596 AliCentrality *centrality = event->GetCentrality();
597 if(centrality){
598 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
599 }
600 else{
601 fCurrCentrality = -1;
602 }
4ec8e76d 603}
604
605//______________________________________________________________________________
606void AliPIDResponse::ExecNewRun()
607{
608 //
609 // Things to Execute upon a new run
610 //
611 SetRecoInfo();
612
613 SetITSParametrisation();
614
615 SetTPCPidResponseMaster();
616 SetTPCParametrisation();
53d016dc 617
618 SetTRDPidResponseMaster();
619 InitializeTRDResponse();
b2138b40 620
621 SetEMCALPidResponseMaster();
622 InitializeEMCALResponse();
4ec8e76d 623
b79db598 624 SetTOFPidResponseMaster();
625 InitializeTOFResponse();
4ec8e76d 626}
627
628//_____________________________________________________
629Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
630{
631 //
632 // Get TPC multiplicity in bins of 150
633 //
634
635 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
636 Double_t tpcMulti=0.;
637 if(vertexTPC){
638 Double_t vertexContribTPC=vertexTPC->GetNContributors();
639 tpcMulti=vertexContribTPC/150.;
640 if (tpcMulti>20.) tpcMulti=20.;
641 }
642
643 return tpcMulti;
644}
645
646//______________________________________________________________________________
647void AliPIDResponse::SetRecoInfo()
648{
649 //
650 // Set reconstruction information
651 //
652
653 //reset information
654 fLHCperiod="";
655 fMCperiodTPC="";
656
657 fBeamType="";
658
659 fBeamType="PP";
660
1436d6bb 661
662 TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
663
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 //
db0e2c5f 845 if(fTRDPIDResponseObject) return;
53d016dc 846 AliOADBContainer contParams("contParams");
847
db0e2c5f 848 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
849 if(statusResponse){
850 AliError("Failed initializing PID Response Object from OADB");
59a8e853 851 } else {
db0e2c5f 852 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
853 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
854 if(!fTRDPIDResponseObject){
855 AliError(Form("TRD Response not found in run %d", fRun));
59a8e853 856 }
857 }
db0e2c5f 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 }
db0e2c5f 869 }
870 */
ea235c90 871}
872
873//______________________________________________________________________________
874void AliPIDResponse::InitializeTRDResponse(){
875 //
876 // Set PID Params and references to the TRD PID response
877 //
db0e2c5f 878 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
879 if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
880 fTRDslicesForPID[0] = 0;
881 fTRDslicesForPID[1] = 7;
882 }
ea235c90 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}