]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/STEERBase/AliPIDResponse.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
An additional label
[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_]*)/.*");
a78fd045 663 TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
1436d6bb 664
4ec8e76d 665 //find the period by run number (UGLY, but not stored in ESD and AOD... )
666 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
667 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
668 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
99e9d5ec 669 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 670 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
671 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
ea235c90 672 else if (fRun>=136851&&fRun<=139517) {
673 fLHCperiod="LHC10H";
674 fMCperiodTPC="LHC10H8";
675 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
676 fBeamType="PBPB";
677 }
3077a03d 678 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
2e97a211 679 //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
680 else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
681 else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
3077a03d 682 else if (fRun>=166529) {
683 fLHCperiod="LHC11H";
684 fMCperiodTPC="LHC11A10";
685 fBeamType="PBPB";
a78fd045 686 if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
3077a03d 687 }
688
80ab5635 689
690 //exception new pp MC productions from 2011
691 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
4a527e08 692 // exception for 11f1
693 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
4ec8e76d 694}
695
696//______________________________________________________________________________
697void AliPIDResponse::SetITSParametrisation()
698{
699 //
700 // Set the ITS parametrisation
701 //
702}
703
704//______________________________________________________________________________
705void AliPIDResponse::SetTPCPidResponseMaster()
706{
707 //
708 // Load the TPC pid response functions from the OADB
709 //
09b50a42 710 //don't load twice for the moment
711 if (fArrPidResponseMaster) return;
712
713
4ec8e76d 714 //reset the PID response functions
715 delete fArrPidResponseMaster;
716 fArrPidResponseMaster=0x0;
717
718 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
719
ea235c90 720 TFile *f=TFile::Open(fileName.Data());
721 if (f && f->IsOpen() && !f->IsZombie()){
722 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 723 }
ea235c90 724 delete f;
4ec8e76d 725
726 if (!fArrPidResponseMaster){
727 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
728 return;
729 }
730 fArrPidResponseMaster->SetOwner();
731}
732
733//______________________________________________________________________________
734void AliPIDResponse::SetTPCParametrisation()
735{
736 //
737 // Change BB parametrisation for current run
738 //
739
740 if (fLHCperiod.IsNull()) {
741 AliFatal("No period set, not changing parametrisation");
742 return;
743 }
744
745 //
746 // Set default parametrisations for data and MC
747 //
748
749 //data type
750 TString datatype="DATA";
751 //in case of mc fRecoPass is per default 1
752 if (fIsMC) {
753 datatype="MC";
754 fRecoPass=1;
755 }
756
757 //
758 //reset old splines
759 //
760 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
761 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
762 }
4a527e08 763
764 // period
765 TString period=fLHCperiod;
766 if (fIsMC) period=fMCperiodTPC;
767
768 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
769 Bool_t found=kFALSE;
4ec8e76d 770 //
771 //set the new PID splines
772 //
4ec8e76d 773 if (fArrPidResponseMaster){
774 TObject *grAll=0x0;
775 //for MC don't use period information
776// if (fIsMC) period="[A-Z0-9]*";
777 //for MC use MC period information
4ec8e76d 778//pattern for the default entry (valid for all particles)
779 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
780
781 //loop over entries and filter them
782 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
783 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
784 if (responseFunction==0x0) continue;
785 TString responseName=responseFunction->GetName();
786
787 if (!reg.MatchB(responseName)) continue;
788
789 TObjArray *arr=reg.MatchS(responseName);
790 TString particleName=arr->At(1)->GetName();
791 delete arr;
792 if (particleName.IsNull()) continue;
793 if (particleName=="ALL") grAll=responseFunction;
794 else {
795 //find particle id
796 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
797 TString particle=AliPID::ParticleName(ispec);
798 particle.ToUpper();
799 if ( particle == particleName ){
800 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
801 fTPCResponse.SetUseDatabase(kTRUE);
802 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
4a527e08 803 found=kTRUE;
4ec8e76d 804 break;
805 }
806 }
807 }
808 }
809
810 //set default response function to all particles which don't have a specific one
811 if (grAll){
812 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
813 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
814 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
815 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
4a527e08 816 found=kTRUE;
4ec8e76d 817 }
818 }
819 }
820 }
4a527e08 821
822 if (!found){
823 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
824 }
4ec8e76d 825 //
826 // Setup resolution parametrisation
827 //
828
829 //default
830 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
831
832 if (fRun>=122195){
833 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
834 }
23425eb2 835 if (fArrPidResponseMaster)
4ec8e76d 836 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
837
838 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
839}
840
ea235c90 841//______________________________________________________________________________
842void AliPIDResponse::SetTRDPidResponseMaster()
843{
844 //
845 // Load the TRD pid params and references from the OADB
846 //
db0e2c5f 847 if(fTRDPIDResponseObject) return;
53d016dc 848 AliOADBContainer contParams("contParams");
849
db0e2c5f 850 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
851 if(statusResponse){
852 AliError("Failed initializing PID Response Object from OADB");
59a8e853 853 } else {
db0e2c5f 854 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
855 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
856 if(!fTRDPIDResponseObject){
857 AliError(Form("TRD Response not found in run %d", fRun));
59a8e853 858 }
859 }
db0e2c5f 860 /*
53d016dc 861 AliOADBContainer contRefs("contRefs");
59a8e853 862 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
863 if(statusRefs){
864 AliInfo("Failed Loading References for TRD");
865 } else {
866 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
867 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
868 if(!fTRDPIDReference){
869 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
870 }
db0e2c5f 871 }
872 */
ea235c90 873}
874
875//______________________________________________________________________________
876void AliPIDResponse::InitializeTRDResponse(){
877 //
878 // Set PID Params and references to the TRD PID response
879 //
db0e2c5f 880 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
f2762b1c 881 SetTRDPIDmethod();
882}
883
884void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
885
886 fTRDResponse.SetPIDmethod(method);
887 if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
888 // backward compatibility for setting with 8 slices
889 fTRDslicesForPID[0] = 0;
890 fTRDslicesForPID[1] = 7;
891 }
892 else{
893 if(method==AliTRDPIDResponse::kLQ1D){
894 fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
895 fTRDslicesForPID[1] = 0;
896 }
897 if(method==AliTRDPIDResponse::kLQ2D){
898 fTRDslicesForPID[0] = 1;
899 fTRDslicesForPID[1] = 7;
db0e2c5f 900 }
f2762b1c 901 }
902 AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
ea235c90 903}
904
b79db598 905//______________________________________________________________________________
906void AliPIDResponse::SetTOFPidResponseMaster()
907{
908 //
909 // Load the TOF pid params from the OADB
910 //
911 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
912 if (oadbf->IsOpen()) {
913 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
914 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
915 if (fTOFPIDParams) delete fTOFPIDParams;
f57f1e18 916 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 917 oadbf->Close();
918 delete oadbc;
919 } else {
920 AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
921 }
922 delete oadbf;
923
924 }
925
926//______________________________________________________________________________
927void AliPIDResponse::InitializeTOFResponse(){
928 //
929 // Set PID Params to the TOF PID response
930 //
b79db598 931 for (Int_t i=0;i<4;i++) {
932 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
933 }
934 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
935
936}
937
938
ea235c90 939//_________________________________________________________________________
940Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
941 //
942 // Check whether track is identified as electron under a given electron efficiency hypothesis
943 //
944 Double_t probs[AliPID::kSPECIES];
945 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
946
99e9d5ec 947 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
948 // Take mean of the TRD momenta in the given tracklets
949 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
950 Int_t nmomenta = 0;
ea235c90 951 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
952 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 953 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 954 }
955 }
99e9d5ec 956 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 957
958 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
959}
960
b2138b40 961//______________________________________________________________________________
962void AliPIDResponse::SetEMCALPidResponseMaster()
963{
964 //
965 // Load the EMCAL pid response functions from the OADB
966 //
967 TObjArray* fEMCALPIDParamsRun = NULL;
968 TObjArray* fEMCALPIDParamsPass = NULL;
969
970 if(fEMCALPIDParams) return;
971 AliOADBContainer contParams("contParams");
972
973 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
974 if(statusPars){
975 AliError("Failed initializing PID Params from OADB");
976 }
977 else {
978 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
979
980 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
981 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
982 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
983
984 if(!fEMCALPIDParams){
f8d39067 985 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 986 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 987
1f631618 988 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
989 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 990 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
991
992 if(!fEMCALPIDParams){
1f631618 993 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 994 }
995 }
996 }
997}
998
999//______________________________________________________________________________
1000void AliPIDResponse::InitializeEMCALResponse(){
1001 //
1002 // Set PID Params to the EMCAL PID response
1003 //
1004 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1005
1006}
5f8db5fe 1007//_________________________________________________________________________
1008void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1009 //
1010 // Set TOF response function
1011 // Input option for event_time used
1012 //
1013
1014 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1015 if(t0spread < 10) t0spread = 80;
1016
1017 // T0 from TOF algorithm
1018
1019 Bool_t flagT0TOF=kFALSE;
1020 Bool_t flagT0T0=kFALSE;
1021 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1022 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1023 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1024
1025 // T0-TOF arrays
1026 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1027 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1028 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1029 estimatedT0event[i]=0.0;
1030 estimatedT0resolution[i]=0.0;
1031 startTimeMask[i] = 0;
1032 }
1033
1034 Float_t resT0A=75,resT0C=65,resT0AC=55;
1035 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1036 flagT0T0=kTRUE;
1037 }
1038
1039
1040 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1041
1042 if (tofHeader) { // read global info and T0-TOF
1043 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1044 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1045 if(t0spread < 10) t0spread = 80;
1046
1047 flagT0TOF=kTRUE;
1048 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1049 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1050 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1051 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1052 }
1053
1054 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1055 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1056 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1057 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1058 Int_t icurrent = (Int_t)ibin->GetAt(j);
1059 startTime[icurrent]=t0Bin->GetAt(j);
1060 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1061 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1062 }
1063 }
1064
1065 // for cut of 3 sigma on t0 spread
1066 Float_t t0cut = 3 * t0spread;
1067 if(t0cut < 500) t0cut = 500;
1068
1069 if(option == kFILL_T0){ // T0-FILL is used
1070 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1071 estimatedT0event[i]=0.0;
1072 estimatedT0resolution[i]=t0spread;
1073 }
1074 fTOFResponse.SetT0event(estimatedT0event);
1075 fTOFResponse.SetT0resolution(estimatedT0resolution);
1076 }
1077
1078 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1079 if(flagT0TOF){
1080 fTOFResponse.SetT0event(startTime);
1081 fTOFResponse.SetT0resolution(startTimeRes);
1082 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1083 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1084 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1085 }
1086 }
1087 else{
1088 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1089 estimatedT0event[i]=0.0;
1090 estimatedT0resolution[i]=t0spread;
1091 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1092 }
1093 fTOFResponse.SetT0event(estimatedT0event);
1094 fTOFResponse.SetT0resolution(estimatedT0resolution);
1095 }
1096 }
1097 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1098 Float_t t0AC=-10000;
1099 Float_t t0A=-10000;
1100 Float_t t0C=-10000;
1101 if(flagT0T0){
1102 t0AC= vevent->GetT0TOF()[0];
1103 t0A= vevent->GetT0TOF()[1];
1104 t0C= vevent->GetT0TOF()[2];
1105 }
1106
1107 Float_t t0t0Best = 0;
1108 Float_t t0t0BestRes = 9999;
1109 Int_t t0used=0;
1110 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1111 t0t0Best = t0AC;
1112 t0t0BestRes = resT0AC;
1113 t0used=6;
1114 }
1115 else if(TMath::Abs(t0C) < t0cut){
1116 t0t0Best = t0C;
1117 t0t0BestRes = resT0C;
1118 t0used=4;
1119 }
1120 else if(TMath::Abs(t0A) < t0cut){
1121 t0t0Best = t0A;
1122 t0t0BestRes = resT0A;
1123 t0used=2;
1124 }
1125
1126 if(flagT0TOF){ // if T0-TOF info is available
1127 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1128 if(t0t0BestRes < 999){
1129 if(startTimeRes[i] < t0spread){
1130 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1131 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1132 estimatedT0event[i]=t0best / wtot;
1133 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1134 startTimeMask[i] = t0used+1;
1135 }
1136 else {
1137 estimatedT0event[i]=t0t0Best;
1138 estimatedT0resolution[i]=t0t0BestRes;
1139 startTimeMask[i] = t0used;
1140 }
1141 }
1142 else{
1143 estimatedT0event[i]=startTime[i];
1144 estimatedT0resolution[i]=startTimeRes[i];
1145 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1146 }
1147 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1148 }
1149 fTOFResponse.SetT0event(estimatedT0event);
1150 fTOFResponse.SetT0resolution(estimatedT0resolution);
1151 }
1152 else{ // if no T0-TOF info is available
1153 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1154 fTOFResponse.SetT0binMask(i,t0used);
1155 if(t0t0BestRes < 999){
1156 estimatedT0event[i]=t0t0Best;
1157 estimatedT0resolution[i]=t0t0BestRes;
1158 }
1159 else{
1160 estimatedT0event[i]=0.0;
1161 estimatedT0resolution[i]=t0spread;
1162 }
1163 }
1164 fTOFResponse.SetT0event(estimatedT0event);
1165 fTOFResponse.SetT0resolution(estimatedT0resolution);
1166 }
1167 }
1168
1169 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1170 Float_t t0AC=-10000;
1171 Float_t t0A=-10000;
1172 Float_t t0C=-10000;
1173 if(flagT0T0){
1174 t0AC= vevent->GetT0TOF()[0];
1175 t0A= vevent->GetT0TOF()[1];
1176 t0C= vevent->GetT0TOF()[2];
1177 }
1178
1179 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1180 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1181 estimatedT0event[i]=t0AC;
1182 estimatedT0resolution[i]=resT0AC;
1183 fTOFResponse.SetT0binMask(i,6);
1184 }
1185 }
1186 else if(TMath::Abs(t0C) < t0cut){
1187 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1188 estimatedT0event[i]=t0C;
1189 estimatedT0resolution[i]=resT0C;
1190 fTOFResponse.SetT0binMask(i,4);
1191 }
1192 }
1193 else if(TMath::Abs(t0A) < t0cut){
1194 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1195 estimatedT0event[i]=t0A;
1196 estimatedT0resolution[i]=resT0A;
1197 fTOFResponse.SetT0binMask(i,2);
1198 }
1199 }
1200 else{
1201 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1202 estimatedT0event[i]=0.0;
1203 estimatedT0resolution[i]=t0spread;
1204 fTOFResponse.SetT0binMask(i,0);
1205 }
1206 }
1207 fTOFResponse.SetT0event(estimatedT0event);
1208 fTOFResponse.SetT0resolution(estimatedT0resolution);
1209 }
1210 delete [] startTime;
1211 delete [] startTimeRes;
1212 delete [] startTimeMask;
1213 delete [] estimatedT0event;
1214 delete [] estimatedT0resolution;
1215}