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 #include <TArrayF.h>
  33
  34 #include <AliVEvent.h>
  35 #include <AliVTrack.h>
  36 #include <AliLog.h>
  37 #include <AliPID.h>
  38 #include <AliOADBContainer.h>
  39 #include <AliTRDPIDResponseObject.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 fTRDPIDResponseObject(0x0),
  71 fTOFtail(1.1),
  72 fTOFPIDParams(0x0),
  73 fEMCALPIDParams(0x0),
  74 fCurrentEvent(0x0),
  75 fCurrCentrality(0.0)
  76 {
  77   //
  78   // default ctor
  79   //
  80   AliLog::SetClassDebugLevel("AliPIDResponse",0);
  81   AliLog::SetClassDebugLevel("AliESDpid",0);
  82   AliLog::SetClassDebugLevel("AliAODpidUtil",0);
  83
  84   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
  85 }
  86
  87 //______________________________________________________________________________
  88 AliPIDResponse::~AliPIDResponse()
  89 {
  90   //
  91   // dtor
  92   //
  93     delete fArrPidResponseMaster;
  94     delete fTRDPIDResponseObject;
  95   if (fTOFPIDParams) delete fTOFPIDParams;
  96 }
  97
  98 //______________________________________________________________________________
  99 AliPIDResponse::AliPIDResponse(const AliPIDResponse &other) :
 100 TNamed(other),
 101 fITSResponse(other.fITSResponse),
 102 fTPCResponse(other.fTPCResponse),
 103 fTRDResponse(other.fTRDResponse),
 104 fTOFResponse(other.fTOFResponse),
 105 fEMCALResponse(other.fEMCALResponse),
 106 fRange(other.fRange),
 107 fITSPIDmethod(other.fITSPIDmethod),
 108 fIsMC(other.fIsMC),
 109 fOADBPath(other.fOADBPath),
 110 fBeamType("PP"),
 111 fLHCperiod(),
 112 fMCperiodTPC(),
 113 fMCperiodUser(other.fMCperiodUser),
 114 fCurrentFile(),
 115 fRecoPass(0),
 116 fRecoPassUser(other.fRecoPassUser),
 117 fRun(0),
 118 fOldRun(0),
 119 fArrPidResponseMaster(0x0),
 120 fResolutionCorrection(0x0),
 121 fTRDPIDResponseObject(0x0),
 122 fTOFtail(1.1),
 123 fTOFPIDParams(0x0),
 124 fEMCALPIDParams(0x0),
 125 fCurrentEvent(0x0),
 126 fCurrCentrality(0.0)
 127 {
 128   //
 129   // copy ctor
 130   //
 131   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 132 }
 133
 134 //______________________________________________________________________________
 135 AliPIDResponse& 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;
 147     fEMCALResponse=other.fEMCALResponse;
 148     fRange=other.fRange;
 149     fITSPIDmethod=other.fITSPIDmethod;
 150     fOADBPath=other.fOADBPath;
 151     fIsMC=other.fIsMC;
 152     fBeamType="PP";
 153     fLHCperiod="";
 154     fMCperiodTPC="";
 155     fMCperiodUser=other.fMCperiodUser;
 156     fCurrentFile="";
 157     fRecoPass=0;
 158     fRecoPassUser=other.fRecoPassUser;
 159     fRun=0;
 160     fOldRun=0;
 161     fArrPidResponseMaster=0x0;
 162     fResolutionCorrection=0x0;
 163     fTRDPIDResponseObject=0x0;
 164     fEMCALPIDParams=0x0;
 165     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 166     fTOFtail=1.1;
 167     fTOFPIDParams=0x0;
 168     fCurrentEvent=other.fCurrentEvent;
 169   }
 170   return *this;
 171 }
 172
 173 //______________________________________________________________________________
 174 Float_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;
 186 //     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
 187 //     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
 188     default: return -999.;
 189   }
 190
 191 }
 192
 193 //______________________________________________________________________________
 194 Float_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){
 209
 210     mom    = track->P();
 211     pt     = track->Pt();
 212     charge = track->Charge();
 213
 214     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
 215
 216     if(matchedClus){
 217
 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       }
 228     }
 229   }
 230
 231   return -999;
 232
 233 }
 234
 235 //______________________________________________________________________________
 236 Float_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
 282 //______________________________________________________________________________
 283 AliPIDResponse::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 //______________________________________________________________________________
 302 AliPIDResponse::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 //______________________________________________________________________________
 357 AliPIDResponse::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 //______________________________________________________________________________
 400 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 401 {
 402   //
 403   // Compute PID response for the
 404   //
 405
 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
 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
 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   //  }
 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);
 425     Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
 426
 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));
 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
 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;
 447     */
 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 //______________________________________________________________________________
 464 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 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;
 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;
 486 }
 487 //______________________________________________________________________________
 488 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 489 {
 490   //
 491   // Compute PID response for the EMCAL
 492   //
 493
 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)
 536   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 537   return kDetNoSignal;
 538
 539 }
 540 //______________________________________________________________________________
 541 AliPIDResponse::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 //______________________________________________________________________________
 552 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 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;
 560   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
 561
 562   track->GetHMPIDpid(p);
 563
 564   return kDetPidOk;
 565 }
 566
 567 //______________________________________________________________________________
 568 void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
 569 {
 570   //
 571   // Apply settings for the current event
 572   //
 573   fRecoPass=pass;
 574
 575   fCurrentEvent=0x0;
 576   if (!event) return;
 577   fCurrentEvent=event;
 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
 592   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
 593
 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   }
 603 }
 604
 605 //______________________________________________________________________________
 606 void AliPIDResponse::ExecNewRun()
 607 {
 608   //
 609   // Things to Execute upon a new run
 610   //
 611   SetRecoInfo();
 612
 613   SetITSParametrisation();
 614
 615   SetTPCPidResponseMaster();
 616   SetTPCParametrisation();
 617
 618   SetTRDPidResponseMaster();
 619   InitializeTRDResponse();
 620
 621   SetEMCALPidResponseMaster();
 622   InitializeEMCALResponse();
 623
 624   SetTOFPidResponseMaster();
 625   InitializeTOFResponse();
 626 }
 627
 628 //_____________________________________________________
 629 Double_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 //______________________________________________________________________________
 647 void 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
 661
 662   TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
 663   TPRegexp reg12a17(".*(LHC12a17[a-z]+)/.*");
 664
 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"; }
 669   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 670   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 671   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 672   else if (fRun>=136851&&fRun<=139517) {
 673     fLHCperiod="LHC10H";
 674     fMCperiodTPC="LHC10H8";
 675     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 676     fBeamType="PBPB";
 677   }
 678   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 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"; }
 682   else if (fRun>=166529) {
 683     fLHCperiod="LHC11H";
 684     fMCperiodTPC="LHC11A10";
 685     fBeamType="PBPB";
 686     if (reg12a17.MatchB(fCurrentFile)) fMCperiodTPC="LHC12A17";
 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(fTRDPIDResponseObject) return;
 848   AliOADBContainer contParams("contParams");
 849
 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");
 853   } else {
 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));
 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 //______________________________________________________________________________
 876 void AliPIDResponse::InitializeTRDResponse(){
 877   //
 878   // Set PID Params and references to the TRD PID response
 879   //
 880     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
 881     SetTRDPIDmethod();
 882 }
 883
 884 void AliPIDResponse::SetTRDPIDmethod(AliTRDPIDResponse::ETRDPIDMethod method){
 885
 886   fTRDResponse.SetPIDmethod(method);
 887   if(fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
 888     // backward compatibility for setting with 8 slices
 889     fTRDslicesForPID[0] = 0;
 890     fTRDslicesForPID[1] = 7;
 891   }
 892   else{
 893     if(method==AliTRDPIDResponse::kLQ1D){
 894       fTRDslicesForPID[0] = 0; // first Slice contains normalized dEdx
 895       fTRDslicesForPID[1] = 0;
 896     }
 897     if(method==AliTRDPIDResponse::kLQ2D){
 898       fTRDslicesForPID[0] = 1;
 899       fTRDslicesForPID[1] = 7;
 900     }
 901   }
 902   AliDebug(1,Form("Slice Range set to %d - %d",fTRDslicesForPID[0],fTRDslicesForPID[1]));
 903 }
 904
 905 //______________________________________________________________________________
 906 void AliPIDResponse::SetTOFPidResponseMaster()
 907 {
 908   //
 909   // Load the TOF pid params from the OADB
 910   //
 911   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
 912   if (oadbf->IsOpen()) {
 913     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
 914     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
 915     if (fTOFPIDParams) delete fTOFPIDParams;
 916     fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
 917     oadbf->Close();
 918     delete oadbc;
 919   } else {
 920     AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
 921   }
 922   delete oadbf;
 923
 924   }
 925
 926 //______________________________________________________________________________
 927 void AliPIDResponse::InitializeTOFResponse(){
 928   //
 929   // Set PID Params to the TOF PID response
 930   //
 931   for (Int_t i=0;i<4;i++) {
 932     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
 933   }
 934   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
 935
 936 }
 937
 938
 939 //_________________________________________________________________________
 940 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
 941   //
 942   // Check whether track is identified as electron under a given electron efficiency hypothesis
 943   //
 944   Double_t probs[AliPID::kSPECIES];
 945   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
 946
 947   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
 948   // Take mean of the TRD momenta in the given tracklets
 949   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
 950   Int_t nmomenta = 0;
 951   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
 952     if(vtrack->GetTRDmomentum(iPl) > 0.){
 953       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
 954     }
 955   }
 956   p = TMath::Mean(nmomenta, trdmomenta);
 957
 958   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
 959 }
 960
 961 //______________________________________________________________________________
 962 void AliPIDResponse::SetEMCALPidResponseMaster()
 963 {
 964   //
 965   // Load the EMCAL pid response functions from the OADB
 966   //
 967   TObjArray* fEMCALPIDParamsRun      = NULL;
 968   TObjArray* fEMCALPIDParamsPass     = NULL;
 969
 970   if(fEMCALPIDParams) return;
 971   AliOADBContainer contParams("contParams");
 972
 973   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
 974   if(statusPars){
 975     AliError("Failed initializing PID Params from OADB");
 976   }
 977   else {
 978     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 979
 980     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
 981     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
 982     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 983
 984     if(!fEMCALPIDParams){
 985       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
 986       AliInfo("Will take the standard LHC11d instead ...");
 987
 988       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
 989       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
 990       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 991
 992       if(!fEMCALPIDParams){
 993         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 994       }
 995     }
 996   }
 997 }
 998
 999 //______________________________________________________________________________
1000 void AliPIDResponse::InitializeEMCALResponse(){
1001   //
1002   // Set PID Params to the EMCAL PID response
1003   //
1004   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
1005
1006 }
1007 //_________________________________________________________________________
1008 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
1009   //
1010   // Set TOF response function
1011   // Input option for event_time used
1012   //
1013
1014     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
1015     if(t0spread < 10) t0spread = 80;
1016
1017     // T0 from TOF algorithm
1018
1019     Bool_t flagT0TOF=kFALSE;
1020     Bool_t flagT0T0=kFALSE;
1021     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1022     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1023     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1024
1025     // T0-TOF arrays
1026     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1027     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1028     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1029       estimatedT0event[i]=0.0;
1030       estimatedT0resolution[i]=0.0;
1031       startTimeMask[i] = 0;
1032     }
1033
1034     Float_t resT0A=75,resT0C=65,resT0AC=55;
1035     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1036         flagT0T0=kTRUE;
1037     }
1038
1039
1040     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1041
1042     if (tofHeader) { // read global info and T0-TOF
1043       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1044       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1045       if(t0spread < 10) t0spread = 80;
1046
1047       flagT0TOF=kTRUE;
1048       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1049         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1050         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1051         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1052       }
1053
1054       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1055       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1056       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1057       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1058         Int_t icurrent = (Int_t)ibin->GetAt(j);
1059         startTime[icurrent]=t0Bin->GetAt(j);
1060         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1061         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1062       }
1063     }
1064
1065     // for cut of 3 sigma on t0 spread
1066     Float_t t0cut = 3 * t0spread;
1067     if(t0cut < 500) t0cut = 500;
1068
1069     if(option == kFILL_T0){ // T0-FILL is used
1070         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1071           estimatedT0event[i]=0.0;
1072           estimatedT0resolution[i]=t0spread;
1073         }
1074         fTOFResponse.SetT0event(estimatedT0event);
1075         fTOFResponse.SetT0resolution(estimatedT0resolution);
1076     }
1077
1078     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1079         if(flagT0TOF){
1080             fTOFResponse.SetT0event(startTime);
1081             fTOFResponse.SetT0resolution(startTimeRes);
1082             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1083               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1084               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1085             }
1086         }
1087         else{
1088             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1089               estimatedT0event[i]=0.0;
1090               estimatedT0resolution[i]=t0spread;
1091               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1092             }
1093             fTOFResponse.SetT0event(estimatedT0event);
1094             fTOFResponse.SetT0resolution(estimatedT0resolution);
1095         }
1096     }
1097     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1098         Float_t t0AC=-10000;
1099         Float_t t0A=-10000;
1100         Float_t t0C=-10000;
1101         if(flagT0T0){
1102             t0AC= vevent->GetT0TOF()[0];
1103             t0A= vevent->GetT0TOF()[1];
1104             t0C= vevent->GetT0TOF()[2];
1105         }
1106
1107         Float_t t0t0Best = 0;
1108         Float_t t0t0BestRes = 9999;
1109         Int_t t0used=0;
1110         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1111             t0t0Best = t0AC;
1112             t0t0BestRes = resT0AC;
1113             t0used=6;
1114         }
1115         else if(TMath::Abs(t0C) < t0cut){
1116             t0t0Best = t0C;
1117             t0t0BestRes = resT0C;
1118             t0used=4;
1119         }
1120         else if(TMath::Abs(t0A) < t0cut){
1121             t0t0Best = t0A;
1122             t0t0BestRes = resT0A;
1123             t0used=2;
1124         }
1125
1126         if(flagT0TOF){ // if T0-TOF info is available
1127             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1128                 if(t0t0BestRes < 999){
1129                   if(startTimeRes[i] < t0spread){
1130                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1131                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1132                     estimatedT0event[i]=t0best / wtot;
1133                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1134                     startTimeMask[i] = t0used+1;
1135                   }
1136                   else {
1137                     estimatedT0event[i]=t0t0Best;
1138                     estimatedT0resolution[i]=t0t0BestRes;
1139                     startTimeMask[i] = t0used;
1140                   }
1141                 }
1142                 else{
1143                   estimatedT0event[i]=startTime[i];
1144                   estimatedT0resolution[i]=startTimeRes[i];
1145                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1146                 }
1147                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1148             }
1149             fTOFResponse.SetT0event(estimatedT0event);
1150             fTOFResponse.SetT0resolution(estimatedT0resolution);
1151         }
1152         else{ // if no T0-TOF info is available
1153             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1154               fTOFResponse.SetT0binMask(i,t0used);
1155               if(t0t0BestRes < 999){
1156                 estimatedT0event[i]=t0t0Best;
1157                 estimatedT0resolution[i]=t0t0BestRes;
1158               }
1159               else{
1160                 estimatedT0event[i]=0.0;
1161                 estimatedT0resolution[i]=t0spread;
1162               }
1163             }
1164             fTOFResponse.SetT0event(estimatedT0event);
1165             fTOFResponse.SetT0resolution(estimatedT0resolution);
1166         }
1167     }
1168
1169     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1170         Float_t t0AC=-10000;
1171         Float_t t0A=-10000;
1172         Float_t t0C=-10000;
1173         if(flagT0T0){
1174             t0AC= vevent->GetT0TOF()[0];
1175             t0A= vevent->GetT0TOF()[1];
1176             t0C= vevent->GetT0TOF()[2];
1177         }
1178
1179         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1180             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1181               estimatedT0event[i]=t0AC;
1182               estimatedT0resolution[i]=resT0AC;
1183               fTOFResponse.SetT0binMask(i,6);
1184             }
1185         }
1186         else if(TMath::Abs(t0C) < t0cut){
1187             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1188               estimatedT0event[i]=t0C;
1189               estimatedT0resolution[i]=resT0C;
1190               fTOFResponse.SetT0binMask(i,4);
1191             }
1192         }
1193         else if(TMath::Abs(t0A) < t0cut){
1194             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1195               estimatedT0event[i]=t0A;
1196               estimatedT0resolution[i]=resT0A;
1197               fTOFResponse.SetT0binMask(i,2);
1198             }
1199         }
1200         else{
1201             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1202               estimatedT0event[i]=0.0;
1203               estimatedT0resolution[i]=t0spread;
1204               fTOFResponse.SetT0binMask(i,0);
1205             }
1206         }
1207         fTOFResponse.SetT0event(estimatedT0event);
1208         fTOFResponse.SetT0resolution(estimatedT0resolution);
1209     }
1210     delete [] startTime;
1211     delete [] startTimeRes;
1212     delete [] startTimeMask;
1213     delete [] estimatedT0event;
1214     delete [] estimatedT0resolution;
1215 }