]> 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
In AddTaskPHOSPi0Flow.C set Cent. Bin past event buffers/lists,
[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),
539a5a59 75fCurrCentrality(0.0),
76fTuneMConData(kFALSE)
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 //
db0e2c5f 94 delete fArrPidResponseMaster;
95 delete fTRDPIDResponseObject;
b79db598 96 if (fTOFPIDParams) delete fTOFPIDParams;
4ec8e76d 97}
98
99//______________________________________________________________________________
100AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
101TNamed(other),
102fITSResponse(other.fITSResponse),
103fTPCResponse(other.fTPCResponse),
104fTRDResponse(other.fTRDResponse),
105fTOFResponse(other.fTOFResponse),
e96b9916 106fEMCALResponse(other.fEMCALResponse),
fd21ec8d 107fRange(other.fRange),
108fITSPIDmethod(other.fITSPIDmethod),
4ec8e76d 109fIsMC(other.fIsMC),
110fOADBPath(other.fOADBPath),
111fBeamType("PP"),
112fLHCperiod(),
113fMCperiodTPC(),
fd21ec8d 114fMCperiodUser(other.fMCperiodUser),
ea235c90 115fCurrentFile(),
4ec8e76d 116fRecoPass(0),
fd21ec8d 117fRecoPassUser(other.fRecoPassUser),
4ec8e76d 118fRun(0),
119fOldRun(0),
120fArrPidResponseMaster(0x0),
121fResolutionCorrection(0x0),
db0e2c5f 122fTRDPIDResponseObject(0x0),
0b39f221 123fTOFtail(1.1),
b79db598 124fTOFPIDParams(0x0),
b2138b40 125fEMCALPIDParams(0x0),
80f28562 126fCurrentEvent(0x0),
539a5a59 127fCurrCentrality(0.0),
128fTuneMConData(kFALSE)
4ec8e76d 129{
130 //
131 // copy ctor
132 //
ea235c90 133 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
4ec8e76d 134}
135
136//______________________________________________________________________________
137AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
138{
139 //
140 // copy ctor
141 //
142 if(this!=&other) {
143 delete fArrPidResponseMaster;
144 TNamed::operator=(other);
145 fITSResponse=other.fITSResponse;
146 fTPCResponse=other.fTPCResponse;
147 fTRDResponse=other.fTRDResponse;
148 fTOFResponse=other.fTOFResponse;
e96b9916 149 fEMCALResponse=other.fEMCALResponse;
fd21ec8d 150 fRange=other.fRange;
151 fITSPIDmethod=other.fITSPIDmethod;
4ec8e76d 152 fOADBPath=other.fOADBPath;
153 fIsMC=other.fIsMC;
154 fBeamType="PP";
155 fLHCperiod="";
156 fMCperiodTPC="";
fd21ec8d 157 fMCperiodUser=other.fMCperiodUser;
ea235c90 158 fCurrentFile="";
4ec8e76d 159 fRecoPass=0;
fd21ec8d 160 fRecoPassUser=other.fRecoPassUser;
4ec8e76d 161 fRun=0;
162 fOldRun=0;
163 fArrPidResponseMaster=0x0;
164 fResolutionCorrection=0x0;
db0e2c5f 165 fTRDPIDResponseObject=0x0;
b2138b40 166 fEMCALPIDParams=0x0;
ea235c90 167 memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
0b39f221 168 fTOFtail=1.1;
b79db598 169 fTOFPIDParams=0x0;
e96b9916 170 fCurrentEvent=other.fCurrentEvent;
4ec8e76d 171 }
172 return *this;
173}
174
fd21ec8d 175//______________________________________________________________________________
176Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const
177{
178 //
179 // NumberOfSigmas for 'detCode'
180 //
181
182 switch (detCode){
183 case kDetITS: return NumberOfSigmasITS(track, type); break;
184 case kDetTPC: return NumberOfSigmasTPC(track, type); break;
185 case kDetTOF: return NumberOfSigmasTOF(track, type); break;
186// case kDetTRD: return ComputeTRDProbability(track, type); break;
187// case kDetPHOS: return ComputePHOSProbability(track, type); break;
e96b9916 188// case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
fd21ec8d 189// case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
190 default: return -999.;
191 }
192
193}
194
e96b9916 195//______________________________________________________________________________
196Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
197
198 AliVCluster *matchedClus = NULL;
199
200 Double_t mom = -1.;
201 Double_t pt = -1.;
202 Double_t EovP = -1.;
203 Double_t fClsE = -1.;
204
205 Int_t nMatchClus = -1;
206 Int_t charge = 0;
207
208 // Track matching
209 nMatchClus = track->GetEMCALcluster();
210 if(nMatchClus > -1){
6d0064aa 211
e96b9916 212 mom = track->P();
213 pt = track->Pt();
214 charge = track->Charge();
215
216 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
217
218 if(matchedClus){
219
6d0064aa 220 // matched cluster is EMCAL
221 if(matchedClus->IsEMCAL()){
222
223 fClsE = matchedClus->E();
224 EovP = fClsE/mom;
225
226
227 // NSigma value really meaningful only for electrons!
228 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
229 }
e96b9916 230 }
231 }
6d0064aa 232
e96b9916 233 return -999;
234
235}
236
6d0064aa 237//______________________________________________________________________________
238Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
239
240 AliVCluster *matchedClus = NULL;
241
242 Double_t mom = -1.;
243 Double_t pt = -1.;
244 Double_t EovP = -1.;
245 Double_t fClsE = -1.;
32fa24d6 246
247 // initialize eop and shower shape parameters
248 eop = -1.;
249 for(Int_t i = 0; i < 4; i++){
250 showershape[i] = -1.;
251 }
6d0064aa 252
253 Int_t nMatchClus = -1;
254 Int_t charge = 0;
255
256 // Track matching
257 nMatchClus = track->GetEMCALcluster();
258 if(nMatchClus > -1){
259
260 mom = track->P();
261 pt = track->Pt();
262 charge = track->Charge();
263
264 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
265
266 if(matchedClus){
267
268 // matched cluster is EMCAL
269 if(matchedClus->IsEMCAL()){
270
271 fClsE = matchedClus->E();
272 EovP = fClsE/mom;
273
274 // fill used EMCAL variables here
275 eop = EovP; // E/p
276 showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
277 showershape[1] = matchedClus->GetM02(); // long axis
278 showershape[2] = matchedClus->GetM20(); // short axis
279 showershape[3] = matchedClus->GetDispersion(); // dispersion
280
281 // NSigma value really meaningful only for electrons!
282 return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge);
283 }
284 }
285 }
286 return -999;
287}
288
289
fd21ec8d 290//______________________________________________________________________________
291AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
292{
293 //
294 // Compute PID response of 'detCode'
295 //
296
297 switch (detCode){
298 case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
299 case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
300 case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
301 case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
302 case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
303 case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
304 case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
305 default: return kDetNoSignal;
306 }
307}
308
309//______________________________________________________________________________
310AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
311{
312 //
313 // Compute PID response for the ITS
314 //
315
316 // set flat distribution (no decision)
317 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
318
319 if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
320 (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
321
322 Double_t mom=track->P();
323 Double_t dedx=track->GetITSsignal();
324 Bool_t isSA=kTRUE;
325 Double_t momITS=mom;
326 ULong_t trStatus=track->GetStatus();
327 if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
328 UChar_t clumap=track->GetITSClusterMap();
329 Int_t nPointsForPid=0;
330 for(Int_t i=2; i<6; i++){
331 if(clumap&(1<<i)) ++nPointsForPid;
332 }
333
334 if(nPointsForPid<3) { // track not to be used for combined PID purposes
335 // track->ResetStatus(AliVTrack::kITSpid);
336 return kDetNoSignal;
337 }
338
339 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
340 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
341 Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
342 Double_t bethe=fITSResponse.Bethe(momITS,mass);
343 Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
344 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
345 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
346 } else {
347 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
348 mismatch=kFALSE;
349 }
350
351 // Check for particles heavier than (AliPID::kSPECIES - 1)
352 // if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
353
354 }
355
356 if (mismatch){
357 for (Int_t j=0; j<AliPID::kSPECIES; j++) p[j]=1./AliPID::kSPECIES;
358 return kDetNoSignal;
359 }
360
361
362 return kDetPidOk;
363}
364//______________________________________________________________________________
365AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
366{
367 //
368 // Compute PID response for the TPC
369 //
370
371 // set flat distribution (no decision)
372 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
373
374 // check quality of the track
375 if ( (track->GetStatus()&AliVTrack::kTPCin )==0 && (track->GetStatus()&AliVTrack::kTPCout)==0 ) return kDetNoSignal;
376
377 Double_t mom = track->GetTPCmomentum();
378
379 Double_t dedx=track->GetTPCsignal();
380 Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
381
539a5a59 382 if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
383
fd21ec8d 384 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
385 AliPID::EParticleType type=AliPID::EParticleType(j);
386 Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
387 Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
388 if (TMath::Abs(dedx-bethe) > fRange*sigma) {
389 p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
390 } else {
391 p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
392 mismatch=kFALSE;
393 }
394
395 // TODO: Light nuclei, also in TPC pid response
396
397 // Check for particles heavier than (AliPID::kSPECIES - 1)
398// if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
399
400 }
401
402 if (mismatch){
403 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
404 return kDetNoSignal;
405 }
406
407 return kDetPidOk;
408}
409//______________________________________________________________________________
410AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
411{
412 //
413 // Compute PID response for the
414 //
415
0b39f221 416 Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
417
fd21ec8d 418 // set flat distribution (no decision)
419 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
420
421 if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
422 if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
423
424 Double_t time[AliPID::kSPECIESN];
425 track->GetIntegratedTimes(time);
426
67376d1d 427 // Double_t sigma[AliPID::kSPECIES];
428 // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
429 // sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
430 // }
fd21ec8d 431
432 Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
433 for (Int_t j=0; j<AliPID::kSPECIES; j++) {
434 AliPID::EParticleType type=AliPID::EParticleType(j);
0b39f221 435 Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
436
67376d1d 437 // Double_t sig = sigma[j];
438 Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
439 Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
0b39f221 440 if (TMath::Abs(nsigmas) > (fRange+2)) {
441 if(nsigmas < fTOFtail)
442 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
443 else
444 p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
445 } else{
446 if(nsigmas < fTOFtail)
447 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
448 else
449 p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
450 }
451
452 /* OLD Gaussian shape
fd21ec8d 453 if (TMath::Abs(nsigmas) > (fRange+2)) {
454 p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
455 } else
456 p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
0b39f221 457 */
fd21ec8d 458
459 if (TMath::Abs(nsigmas)<5.){
460 Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
461 if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
462 }
463 }
464
465 if (mismatch){
466 return kDetMismatch;
467 }
468
469 // TODO: Light nuclei
470
471 return kDetPidOk;
472}
473//______________________________________________________________________________
ea235c90 474AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 475{
476 //
477 // Compute PID response for the
478 //
479
480 // set flat distribution (no decision)
481 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 482 if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
483
484 Float_t mom[6];
485 Double_t dedx[48]; // Allocate space for the maximum number of TRD slices
486 Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
487 AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d", fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
488 for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
489 mom[ilayer] = track->GetTRDmomentum(ilayer);
490 for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
491 dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
492 }
493 }
494 fTRDResponse.GetResponse(nslices, dedx, mom, p);
495 return kDetPidOk;
fd21ec8d 496}
497//______________________________________________________________________________
e96b9916 498AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 499{
500 //
501 // Compute PID response for the EMCAL
502 //
503
e96b9916 504 AliVCluster *matchedClus = NULL;
505
506 Double_t mom = -1.;
507 Double_t pt = -1.;
508 Double_t EovP = -1.;
509 Double_t fClsE = -1.;
510
511 Int_t nMatchClus = -1;
512 Int_t charge = 0;
513
514 // Track matching
515 nMatchClus = track->GetEMCALcluster();
516
517 if(nMatchClus > -1){
518
519 mom = track->P();
520 pt = track->Pt();
521 charge = track->Charge();
522
523 matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
524
525 if(matchedClus){
526
527 // matched cluster is EMCAL
528 if(matchedClus->IsEMCAL()){
529
530 fClsE = matchedClus->E();
531 EovP = fClsE/mom;
532
533
534 // compute the probabilities
535 if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
536
537 // in case everything is OK
538 return kDetPidOk;
539
540 }
541 }
542 }
543 }
544
545 // in all other cases set flat distribution (no decision)
fd21ec8d 546 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
547 return kDetNoSignal;
e96b9916 548
fd21ec8d 549}
550//______________________________________________________________________________
551AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
552{
553 //
554 // Compute PID response for the PHOS
555 //
556
557 // set flat distribution (no decision)
558 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
559 return kDetNoSignal;
560}
561//______________________________________________________________________________
ea235c90 562AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
fd21ec8d 563{
564 //
565 // Compute PID response for the HMPID
566 //
567
568 // set flat distribution (no decision)
569 for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
ea235c90 570 if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
571
572 track->GetHMPIDpid(p);
573
574 return kDetPidOk;
fd21ec8d 575}
576
4ec8e76d 577//______________________________________________________________________________
578void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
579{
580 //
581 // Apply settings for the current event
582 //
583 fRecoPass=pass;
e96b9916 584
585 fCurrentEvent=0x0;
4ec8e76d 586 if (!event) return;
e96b9916 587 fCurrentEvent=event;
4ec8e76d 588 fRun=event->GetRunNumber();
589
590 if (fRun!=fOldRun){
591 ExecNewRun();
592 fOldRun=fRun;
593 }
594
595 //TPC resolution parametrisation PbPb
596 if ( fResolutionCorrection ){
597 Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
598 fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
599 }
600
601 //TOF resolution
b79db598 602 SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
603
80f28562 604
605 // Get and set centrality
606 AliCentrality *centrality = event->GetCentrality();
607 if(centrality){
608 fCurrCentrality = centrality->GetCentralityPercentile("V0M");
609 }
610 else{
611 fCurrCentrality = -1;
612 }
4ec8e76d 613}
614
615//______________________________________________________________________________
616void AliPIDResponse::ExecNewRun()
617{
618 //
619 // Things to Execute upon a new run
620 //
621 SetRecoInfo();
622
623 SetITSParametrisation();
624
625 SetTPCPidResponseMaster();
626 SetTPCParametrisation();
53d016dc 627
628 SetTRDPidResponseMaster();
629 InitializeTRDResponse();
b2138b40 630
631 SetEMCALPidResponseMaster();
632 InitializeEMCALResponse();
4ec8e76d 633
b79db598 634 SetTOFPidResponseMaster();
635 InitializeTOFResponse();
4ec8e76d 636}
637
638//_____________________________________________________
639Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
640{
641 //
642 // Get TPC multiplicity in bins of 150
643 //
644
645 const AliVVertex* vertexTPC = event->GetPrimaryVertex();
646 Double_t tpcMulti=0.;
647 if(vertexTPC){
648 Double_t vertexContribTPC=vertexTPC->GetNContributors();
649 tpcMulti=vertexContribTPC/150.;
650 if (tpcMulti>20.) tpcMulti=20.;
651 }
652
653 return tpcMulti;
654}
655
656//______________________________________________________________________________
657void AliPIDResponse::SetRecoInfo()
658{
659 //
660 // Set reconstruction information
661 //
662
663 //reset information
664 fLHCperiod="";
665 fMCperiodTPC="";
666
667 fBeamType="";
668
669 fBeamType="PP";
670
1436d6bb 671
672 TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
a78fd045 673 TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
1436d6bb 674
4ec8e76d 675 //find the period by run number (UGLY, but not stored in ESD and AOD... )
676 if (fRun>=114737&&fRun<=117223) { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1"; }
677 else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1"; }
678 else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
99e9d5ec 679 else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
4ec8e76d 680 else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
681 else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
ea235c90 682 else if (fRun>=136851&&fRun<=139517) {
683 fLHCperiod="LHC10H";
684 fMCperiodTPC="LHC10H8";
685 if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
686 fBeamType="PBPB";
687 }
3077a03d 688 else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
2e97a211 689 //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
690 else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
691 else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
3077a03d 692 else if (fRun>=166529) {
693 fLHCperiod="LHC11H";
694 fMCperiodTPC="LHC11A10";
695 fBeamType="PBPB";
a78fd045 696 if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
3077a03d 697 }
698
80ab5635 699
700 //exception new pp MC productions from 2011
701 if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
4a527e08 702 // exception for 11f1
703 if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
4ec8e76d 704}
705
706//______________________________________________________________________________
707void AliPIDResponse::SetITSParametrisation()
708{
709 //
710 // Set the ITS parametrisation
711 //
712}
713
714//______________________________________________________________________________
715void AliPIDResponse::SetTPCPidResponseMaster()
716{
717 //
718 // Load the TPC pid response functions from the OADB
719 //
09b50a42 720 //don't load twice for the moment
721 if (fArrPidResponseMaster) return;
722
723
4ec8e76d 724 //reset the PID response functions
725 delete fArrPidResponseMaster;
726 fArrPidResponseMaster=0x0;
727
728 TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
729
ea235c90 730 TFile *f=TFile::Open(fileName.Data());
731 if (f && f->IsOpen() && !f->IsZombie()){
732 fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
4ec8e76d 733 }
ea235c90 734 delete f;
4ec8e76d 735
736 if (!fArrPidResponseMaster){
737 AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
738 return;
739 }
740 fArrPidResponseMaster->SetOwner();
741}
742
743//______________________________________________________________________________
744void AliPIDResponse::SetTPCParametrisation()
745{
746 //
747 // Change BB parametrisation for current run
748 //
749
750 if (fLHCperiod.IsNull()) {
751 AliFatal("No period set, not changing parametrisation");
752 return;
753 }
754
755 //
756 // Set default parametrisations for data and MC
757 //
758
759 //data type
760 TString datatype="DATA";
761 //in case of mc fRecoPass is per default 1
762 if (fIsMC) {
539a5a59 763 if(!fTuneMConData) datatype="MC";
764 fRecoPass=1;
4ec8e76d 765 }
766
767 //
768 //reset old splines
769 //
770 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
771 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
772 }
4a527e08 773
774 // period
775 TString period=fLHCperiod;
539a5a59 776 if (fIsMC && !fTuneMConData) period=fMCperiodTPC;
4a527e08 777
778 AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
779 Bool_t found=kFALSE;
4ec8e76d 780 //
781 //set the new PID splines
782 //
4ec8e76d 783 if (fArrPidResponseMaster){
539a5a59 784 Int_t recopass = fRecoPass;
785 if(fTuneMConData) recopass = fRecoPassUser;
4ec8e76d 786 TObject *grAll=0x0;
787 //for MC don't use period information
788// if (fIsMC) period="[A-Z0-9]*";
789 //for MC use MC period information
4ec8e76d 790//pattern for the default entry (valid for all particles)
539a5a59 791 TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),recopass,fBeamType.Data()));
792
4ec8e76d 793 //loop over entries and filter them
794 for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
795 TObject *responseFunction=fArrPidResponseMaster->At(iresp);
796 if (responseFunction==0x0) continue;
797 TString responseName=responseFunction->GetName();
798
799 if (!reg.MatchB(responseName)) continue;
800
801 TObjArray *arr=reg.MatchS(responseName);
802 TString particleName=arr->At(1)->GetName();
803 delete arr;
804 if (particleName.IsNull()) continue;
805 if (particleName=="ALL") grAll=responseFunction;
806 else {
807 //find particle id
808 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
809 TString particle=AliPID::ParticleName(ispec);
810 particle.ToUpper();
811 if ( particle == particleName ){
812 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
813 fTPCResponse.SetUseDatabase(kTRUE);
814 AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
4a527e08 815 found=kTRUE;
4ec8e76d 816 break;
817 }
818 }
819 }
820 }
821
822 //set default response function to all particles which don't have a specific one
823 if (grAll){
824 for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
825 if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
826 fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
827 AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
4a527e08 828 found=kTRUE;
4ec8e76d 829 }
830 }
831 }
832 }
4a527e08 833
834 if (!found){
835 AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
836 }
4ec8e76d 837 //
838 // Setup resolution parametrisation
839 //
840
841 //default
842 fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
843
844 if (fRun>=122195){
845 fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
846 }
23425eb2 847 if (fArrPidResponseMaster)
4ec8e76d 848 fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
849
850 if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
851}
852
ea235c90 853//______________________________________________________________________________
854void AliPIDResponse::SetTRDPidResponseMaster()
855{
856 //
857 // Load the TRD pid params and references from the OADB
858 //
db0e2c5f 859 if(fTRDPIDResponseObject) return;
53d016dc 860 AliOADBContainer contParams("contParams");
861
db0e2c5f 862 Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
863 if(statusResponse){
864 AliError("Failed initializing PID Response Object from OADB");
59a8e853 865 } else {
db0e2c5f 866 AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
867 fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
868 if(!fTRDPIDResponseObject){
869 AliError(Form("TRD Response not found in run %d", fRun));
59a8e853 870 }
871 }
db0e2c5f 872 /*
53d016dc 873 AliOADBContainer contRefs("contRefs");
59a8e853 874 Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
875 if(statusRefs){
876 AliInfo("Failed Loading References for TRD");
877 } else {
878 AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
879 fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
880 if(!fTRDPIDReference){
881 AliError(Form("TRD References not found in OADB Container for run %d", fRun));
882 }
db0e2c5f 883 }
884 */
ea235c90 885}
886
887//______________________________________________________________________________
888void AliPIDResponse::InitializeTRDResponse(){
889 //
890 // Set PID Params and references to the TRD PID response
891 //
db0e2c5f 892 fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
f2762b1c 893 SetTRDPIDmethod();
894}
895
896void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
897
898 fTRDResponse.SetPIDmethod(method);
899 if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
900 // backward compatibility for setting with 8 slices
901 fTRDslicesForPID[0] = 0;
902 fTRDslicesForPID[1] = 7;
903 }
904 else{
905 if(method==AliTRDPIDResponse::kLQ1D){
906 fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
907 fTRDslicesForPID[1] = 0;
908 }
909 if(method==AliTRDPIDResponse::kLQ2D){
910 fTRDslicesForPID[0] = 1;
911 fTRDslicesForPID[1] = 7;
db0e2c5f 912 }
f2762b1c 913 }
914 AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
ea235c90 915}
916
b79db598 917//______________________________________________________________________________
918void AliPIDResponse::SetTOFPidResponseMaster()
919{
920 //
921 // Load the TOF pid params from the OADB
922 //
923 TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
924 if (oadbf->IsOpen()) {
925 AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
926 AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
927 if (fTOFPIDParams) delete fTOFPIDParams;
f57f1e18 928 fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
b79db598 929 oadbf->Close();
930 delete oadbc;
931 } else {
932 AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
933 }
934 delete oadbf;
935
936 }
937
938//______________________________________________________________________________
939void AliPIDResponse::InitializeTOFResponse(){
940 //
941 // Set PID Params to the TOF PID response
942 //
b79db598 943 for (Int_t i=0;i<4;i++) {
944 fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
945 }
946 fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
947
948}
949
950
ea235c90 951//_________________________________________________________________________
952Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
953 //
954 // Check whether track is identified as electron under a given electron efficiency hypothesis
955 //
956 Double_t probs[AliPID::kSPECIES];
957 ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
958
99e9d5ec 959 Int_t ntracklets = vtrack->GetTRDntrackletsPID();
960 // Take mean of the TRD momenta in the given tracklets
961 Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
962 Int_t nmomenta = 0;
ea235c90 963 for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
964 if(vtrack->GetTRDmomentum(iPl) > 0.){
99e9d5ec 965 trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
ea235c90 966 }
967 }
99e9d5ec 968 p = TMath::Mean(nmomenta, trdmomenta);
ea235c90 969
970 return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
971}
972
b2138b40 973//______________________________________________________________________________
974void AliPIDResponse::SetEMCALPidResponseMaster()
975{
976 //
977 // Load the EMCAL pid response functions from the OADB
978 //
979 TObjArray* fEMCALPIDParamsRun = NULL;
980 TObjArray* fEMCALPIDParamsPass = NULL;
981
982 if(fEMCALPIDParams) return;
983 AliOADBContainer contParams("contParams");
984
985 Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
986 if(statusPars){
987 AliError("Failed initializing PID Params from OADB");
988 }
989 else {
990 AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
991
992 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
993 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
994 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
995
996 if(!fEMCALPIDParams){
f8d39067 997 AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
1f631618 998 AliInfo("Will take the standard LHC11d instead ...");
b2138b40 999
1f631618 1000 fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
1001 if(fEMCALPIDParamsRun) fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
b2138b40 1002 if(fEMCALPIDParamsPass) fEMCALPIDParams = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
1003
1004 if(!fEMCALPIDParams){
1f631618 1005 AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
b2138b40 1006 }
1007 }
1008 }
1009}
1010
1011//______________________________________________________________________________
1012void AliPIDResponse::InitializeEMCALResponse(){
1013 //
1014 // Set PID Params to the EMCAL PID response
1015 //
1016 fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1017
1018}
5f8db5fe 1019//_________________________________________________________________________
1020void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1021 //
1022 // Set TOF response function
1023 // Input option for event_time used
1024 //
1025
1026 Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1027 if(t0spread < 10) t0spread = 80;
1028
1029 // T0 from TOF algorithm
1030
1031 Bool_t flagT0TOF=kFALSE;
1032 Bool_t flagT0T0=kFALSE;
1033 Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1034 Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1035 Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1036
1037 // T0-TOF arrays
1038 Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1039 Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1040 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1041 estimatedT0event[i]=0.0;
1042 estimatedT0resolution[i]=0.0;
1043 startTimeMask[i] = 0;
1044 }
1045
1046 Float_t resT0A=75,resT0C=65,resT0AC=55;
1047 if(vevent->GetT0TOF()){ // check if T0 detector information is available
1048 flagT0T0=kTRUE;
1049 }
1050
1051
1052 AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1053
1054 if (tofHeader) { // read global info and T0-TOF
1055 fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1056 t0spread = tofHeader->GetT0spread(); // read t0 sprad
1057 if(t0spread < 10) t0spread = 80;
1058
1059 flagT0TOF=kTRUE;
1060 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1061 startTime[i]=tofHeader->GetDefaultEventTimeVal();
1062 startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1063 if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1064 }
1065
1066 TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1067 TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1068 TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1069 for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1070 Int_t icurrent = (Int_t)ibin->GetAt(j);
1071 startTime[icurrent]=t0Bin->GetAt(j);
1072 startTimeRes[icurrent]=t0ResBin->GetAt(j);
1073 if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1074 }
1075 }
1076
1077 // for cut of 3 sigma on t0 spread
1078 Float_t t0cut = 3 * t0spread;
1079 if(t0cut < 500) t0cut = 500;
1080
1081 if(option == kFILL_T0){ // T0-FILL is used
1082 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1083 estimatedT0event[i]=0.0;
1084 estimatedT0resolution[i]=t0spread;
1085 }
1086 fTOFResponse.SetT0event(estimatedT0event);
1087 fTOFResponse.SetT0resolution(estimatedT0resolution);
1088 }
1089
1090 if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1091 if(flagT0TOF){
1092 fTOFResponse.SetT0event(startTime);
1093 fTOFResponse.SetT0resolution(startTimeRes);
1094 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1095 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1096 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1097 }
1098 }
1099 else{
1100 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1101 estimatedT0event[i]=0.0;
1102 estimatedT0resolution[i]=t0spread;
1103 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1104 }
1105 fTOFResponse.SetT0event(estimatedT0event);
1106 fTOFResponse.SetT0resolution(estimatedT0resolution);
1107 }
1108 }
1109 else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1110 Float_t t0AC=-10000;
1111 Float_t t0A=-10000;
1112 Float_t t0C=-10000;
1113 if(flagT0T0){
1114 t0AC= vevent->GetT0TOF()[0];
1115 t0A= vevent->GetT0TOF()[1];
1116 t0C= vevent->GetT0TOF()[2];
1117 }
1118
1119 Float_t t0t0Best = 0;
1120 Float_t t0t0BestRes = 9999;
1121 Int_t t0used=0;
1122 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1123 t0t0Best = t0AC;
1124 t0t0BestRes = resT0AC;
1125 t0used=6;
1126 }
1127 else if(TMath::Abs(t0C) < t0cut){
1128 t0t0Best = t0C;
1129 t0t0BestRes = resT0C;
1130 t0used=4;
1131 }
1132 else if(TMath::Abs(t0A) < t0cut){
1133 t0t0Best = t0A;
1134 t0t0BestRes = resT0A;
1135 t0used=2;
1136 }
1137
1138 if(flagT0TOF){ // if T0-TOF info is available
1139 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1140 if(t0t0BestRes < 999){
1141 if(startTimeRes[i] < t0spread){
1142 Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1143 Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1144 estimatedT0event[i]=t0best / wtot;
1145 estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1146 startTimeMask[i] = t0used+1;
1147 }
1148 else {
1149 estimatedT0event[i]=t0t0Best;
1150 estimatedT0resolution[i]=t0t0BestRes;
1151 startTimeMask[i] = t0used;
1152 }
1153 }
1154 else{
1155 estimatedT0event[i]=startTime[i];
1156 estimatedT0resolution[i]=startTimeRes[i];
1157 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1158 }
1159 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1160 }
1161 fTOFResponse.SetT0event(estimatedT0event);
1162 fTOFResponse.SetT0resolution(estimatedT0resolution);
1163 }
1164 else{ // if no T0-TOF info is available
1165 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1166 fTOFResponse.SetT0binMask(i,t0used);
1167 if(t0t0BestRes < 999){
1168 estimatedT0event[i]=t0t0Best;
1169 estimatedT0resolution[i]=t0t0BestRes;
1170 }
1171 else{
1172 estimatedT0event[i]=0.0;
1173 estimatedT0resolution[i]=t0spread;
1174 }
1175 }
1176 fTOFResponse.SetT0event(estimatedT0event);
1177 fTOFResponse.SetT0resolution(estimatedT0resolution);
1178 }
1179 }
1180
1181 else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1182 Float_t t0AC=-10000;
1183 Float_t t0A=-10000;
1184 Float_t t0C=-10000;
1185 if(flagT0T0){
1186 t0AC= vevent->GetT0TOF()[0];
1187 t0A= vevent->GetT0TOF()[1];
1188 t0C= vevent->GetT0TOF()[2];
1189 }
1190
1191 if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1192 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1193 estimatedT0event[i]=t0AC;
1194 estimatedT0resolution[i]=resT0AC;
1195 fTOFResponse.SetT0binMask(i,6);
1196 }
1197 }
1198 else if(TMath::Abs(t0C) < t0cut){
1199 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1200 estimatedT0event[i]=t0C;
1201 estimatedT0resolution[i]=resT0C;
1202 fTOFResponse.SetT0binMask(i,4);
1203 }
1204 }
1205 else if(TMath::Abs(t0A) < t0cut){
1206 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1207 estimatedT0event[i]=t0A;
1208 estimatedT0resolution[i]=resT0A;
1209 fTOFResponse.SetT0binMask(i,2);
1210 }
1211 }
1212 else{
1213 for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1214 estimatedT0event[i]=0.0;
1215 estimatedT0resolution[i]=t0spread;
1216 fTOFResponse.SetT0binMask(i,0);
1217 }
1218 }
1219 fTOFResponse.SetT0event(estimatedT0event);
1220 fTOFResponse.SetT0resolution(estimatedT0resolution);
1221 }
1222 delete [] startTime;
1223 delete [] startTimeRes;
1224 delete [] startTimeMask;
1225 delete [] estimatedT0event;
1226 delete [] estimatedT0resolution;
1227}