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",10);
  82   AliLog::SetClassDebugLevel("AliESDpid",10);
  83   AliLog::SetClassDebugLevel("AliAODpidUtil",10);
  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     if (TMath::Abs(nsigmas) > (fRange+2)) {
 433       if(nsigmas < fTOFtail)
 434         p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
 435       else
 436         p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
 437     } else{
 438       if(nsigmas < fTOFtail)
 439         p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
 440       else
 441         p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
 442     }
 443
 444     /* OLD Gaussian shape
 445     if (TMath::Abs(nsigmas) > (fRange+2)) {
 446       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
 447     } else
 448       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
 449     */
 450
 451     if (TMath::Abs(nsigmas)<5.){
 452       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
 453       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
 454     }
 455   }
 456
 457   if (mismatch){
 458     return kDetMismatch;
 459   }
 460
 461     // TODO: Light nuclei
 462
 463   return kDetPidOk;
 464 }
 465 //______________________________________________________________________________
 466 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 467 {
 468   //
 469   // Compute PID response for the
 470   //
 471
 472   // set flat distribution (no decision)
 473   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 474   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
 475
 476   Float_t mom[6];
 477   Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
 478   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
 479   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
 480   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
 481     mom[ilayer] = track->GetTRDmomentum(ilayer);
 482     for(UInt_t islice = fTRDslicesForPID[0]; islice <= fTRDslicesForPID[1]; islice++){
 483       dedx[ilayer*nslices+islice-fTRDslicesForPID[0]] = track->GetTRDslice(ilayer, islice);
 484     }
 485   }
 486   fTRDResponse.GetResponse(nslices, dedx, mom, p);
 487   return kDetPidOk;
 488 }
 489 //______________________________________________________________________________
 490 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 491 {
 492   //
 493   // Compute PID response for the EMCAL
 494   //
 495
 496   AliVCluster *matchedClus = NULL;
 497
 498   Double_t mom     = -1.;
 499   Double_t pt      = -1.;
 500   Double_t EovP    = -1.;
 501   Double_t fClsE   = -1.;
 502
 503   Int_t nMatchClus = -1;
 504   Int_t charge     = 0;
 505
 506   // Track matching
 507   nMatchClus = track->GetEMCALcluster();
 508
 509   if(nMatchClus > -1){
 510
 511     mom    = track->P();
 512     pt     = track->Pt();
 513     charge = track->Charge();
 514
 515     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 516
 517     if(matchedClus){
 518
 519     // matched cluster is EMCAL
 520     if(matchedClus->IsEMCAL()){
 521
 522       fClsE       = matchedClus->E();
 523       EovP        = fClsE/mom;
 524
 525
 526       // compute the probabilities
 527       if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){
 528
 529         // in case everything is OK
 530         return kDetPidOk;
 531
 532       }
 533     }
 534   }
 535   }
 536
 537   // in all other cases set flat distribution (no decision)
 538   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 539   return kDetNoSignal;
 540
 541 }
 542 //______________________________________________________________________________
 543 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliVTrack */*track*/, Int_t nSpecies, Double_t p[]) const
 544 {
 545   //
 546   // Compute PID response for the PHOS
 547   //
 548
 549   // set flat distribution (no decision)
 550   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 551   return kDetNoSignal;
 552 }
 553 //______________________________________________________________________________
 554 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 555 {
 556   //
 557   // Compute PID response for the HMPID
 558   //
 559
 560   // set flat distribution (no decision)
 561   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 562   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
 563
 564   track->GetHMPIDpid(p);
 565
 566   return kDetPidOk;
 567 }
 568
 569 //______________________________________________________________________________
 570 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
 571 {
 572   //
 573   // Apply settings for the current event
 574   //
 575   fRecoPass=pass;
 576
 577   fCurrentEvent=0x0;
 578   if (!event) return;
 579   fCurrentEvent=event;
 580   fRun=event->GetRunNumber();
 581
 582   if (fRun!=fOldRun){
 583     ExecNewRun();
 584     fOldRun=fRun;
 585   }
 586
 587   //TPC resolution parametrisation PbPb
 588   if ( fResolutionCorrection ){
 589     Double_t corrSigma=fResolutionCorrection->Eval(GetTPCMultiplicityBin(event));
 590     fTPCResponse.SetSigma(3.79301e-03*corrSigma, 2.21280e+04);
 591   }
 592
 593   //TOF resolution
 594   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 595
 596
 597   // Get and set centrality
 598   AliCentrality *centrality = event->GetCentrality();
 599   if(centrality){
 600     fCurrCentrality = centrality->GetCentralityPercentile("V0M");
 601   }
 602   else{
 603     fCurrCentrality = -1;
 604   }
 605 }
 606
 607 //______________________________________________________________________________
 608 void AliPIDResponse::ExecNewRun()
 609 {
 610   //
 611   // Things to Execute upon a new run
 612   //
 613   SetRecoInfo();
 614
 615   SetITSParametrisation();
 616
 617   SetTPCPidResponseMaster();
 618   SetTPCParametrisation();
 619
 620   SetTRDPidResponseMaster();
 621   InitializeTRDResponse();
 622
 623   SetEMCALPidResponseMaster();
 624   InitializeEMCALResponse();
 625
 626   SetTOFPidResponseMaster();
 627   InitializeTOFResponse();
 628 }
 629
 630 //_____________________________________________________
 631 Double_t AliPIDResponse::GetTPCMultiplicityBin(const AliVEvent * const event)
 632 {
 633   //
 634   // Get TPC multiplicity in bins of 150
 635   //
 636
 637   const AliVVertex* vertexTPC = event->GetPrimaryVertex();
 638   Double_t tpcMulti=0.;
 639   if(vertexTPC){
 640     Double_t vertexContribTPC=vertexTPC->GetNContributors();
 641     tpcMulti=vertexContribTPC/150.;
 642     if (tpcMulti>20.) tpcMulti=20.;
 643   }
 644
 645   return tpcMulti;
 646 }
 647
 648 //______________________________________________________________________________
 649 void AliPIDResponse::SetRecoInfo()
 650 {
 651   //
 652   // Set reconstruction information
 653   //
 654
 655   //reset information
 656   fLHCperiod="";
 657   fMCperiodTPC="";
 658
 659   fBeamType="";
 660
 661   fBeamType="PP";
 662
 663
 664   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
 665
 666   //find the period by run number (UGLY, but not stored in ESD and AOD... )
 667   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
 668   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
 669   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
 670   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 671   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 672   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 673   else if (fRun>=136851&&fRun<=139517) {
 674     fLHCperiod="LHC10H";
 675     fMCperiodTPC="LHC10H8";
 676     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 677     fBeamType="PBPB";
 678   }
 679   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 680   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
 681   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 682   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 683   else if (fRun>=166529) {
 684     fLHCperiod="LHC11H";
 685     fMCperiodTPC="LHC11A10";
 686     fBeamType="PBPB";
 687   }
 688
 689
 690   //exception new pp MC productions from 2011
 691   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 692   // exception for 11f1
 693   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 694 }
 695
 696 //______________________________________________________________________________
 697 void AliPIDResponse::SetITSParametrisation()
 698 {
 699   //
 700   // Set the ITS parametrisation
 701   //
 702 }
 703
 704 //______________________________________________________________________________
 705 void AliPIDResponse::SetTPCPidResponseMaster()
 706 {
 707   //
 708   // Load the TPC pid response functions from the OADB
 709   //
 710   //don't load twice for the moment
 711    if (fArrPidResponseMaster) return;
 712
 713
 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
 720   TFile *f=TFile::Open(fileName.Data());
 721   if (f && f->IsOpen() && !f->IsZombie()){
 722     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
 723   }
 724   delete f;
 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 //______________________________________________________________________________
 734 void 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   }
 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;
 770   //
 771   //set the new PID splines
 772   //
 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
 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()));
 803             found=kTRUE;
 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()));
 816           found=kTRUE;
 817         }
 818       }
 819     }
 820   }
 821
 822   if (!found){
 823     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 824   }
 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   }
 835   if (fArrPidResponseMaster)
 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
 841 //______________________________________________________________________________
 842 void AliPIDResponse::SetTRDPidResponseMaster()
 843 {
 844   //
 845   // Load the TRD pid params and references from the OADB
 846   //
 847   if(fTRDPIDParams) return;
 848   AliOADBContainer contParams("contParams");
 849
 850   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()), "AliTRDPIDParams");
 851   if(statusPars){
 852     AliError("Failed initializing PID Params from OADB");
 853   } else {
 854     AliInfo(Form("Loading TRD Params from %s/COMMON/PID/data/TRDPIDParams.root", fOADBPath.Data()));
 855     fTRDPIDParams = dynamic_cast<AliTRDPIDParams *>(contParams.GetObject(fRun));
 856     if(!fTRDPIDParams){
 857       AliError(Form("TRD Params not found in run %d", fRun));
 858     }
 859   }
 860
 861   AliOADBContainer contRefs("contRefs");
 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     }
 871   }
 872 }
 873
 874 //______________________________________________________________________________
 875 void AliPIDResponse::InitializeTRDResponse(){
 876   //
 877   // Set PID Params and references to the TRD PID response
 878   //
 879   fTRDResponse.SetPIDParams(fTRDPIDParams);
 880   fTRDResponse.Load(fTRDPIDReference);
 881   if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
 882     fTRDslicesForPID[0] = 0;
 883     fTRDslicesForPID[1] = 7;
 884   }
 885 }
 886
 887 //______________________________________________________________________________
 888 void 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;
 898     fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
 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 //______________________________________________________________________________
 909 void AliPIDResponse::InitializeTOFResponse(){
 910   //
 911   // Set PID Params to the TOF PID response
 912   //
 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
 921 //_________________________________________________________________________
 922 Bool_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
 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;
 933   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
 934     if(vtrack->GetTRDmomentum(iPl) > 0.){
 935       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
 936     }
 937   }
 938   p = TMath::Mean(nmomenta, trdmomenta);
 939
 940   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
 941 }
 942
 943 //______________________________________________________________________________
 944 void 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){
 967       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
 968       AliInfo("Will take the standard LHC11d instead ...");
 969
 970       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
 971       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
 972       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 973
 974       if(!fEMCALPIDParams){
 975         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 976       }
 977     }
 978   }
 979 }
 980
 981 //______________________________________________________________________________
 982 void AliPIDResponse::InitializeEMCALResponse(){
 983   //
 984   // Set PID Params to the EMCAL PID response
 985   //
 986   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
 987
 988 }
 989 //_________________________________________________________________________
 990 void 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 }