STEER/STEERBase/AliPIDResponse.cxx

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