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