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
 664   //find the period by run number (UGLY, but not stored in ESD and AOD... )
 665   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
 666   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
 667   else if (fRun>=122195&&fRun<=126437) { fLHCperiod="LHC10D"; fMCperiodTPC="LHC10F6A"; }
 668   else if (fRun>=127710&&fRun<=130850) { fLHCperiod="LHC10E"; fMCperiodTPC="LHC10F6A"; }
 669   else if (fRun>=133004&&fRun<=135029) { fLHCperiod="LHC10F"; fMCperiodTPC="LHC10F6A"; }
 670   else if (fRun>=135654&&fRun<=136377) { fLHCperiod="LHC10G"; fMCperiodTPC="LHC10F6A"; }
 671   else if (fRun>=136851&&fRun<=139517) {
 672     fLHCperiod="LHC10H";
 673     fMCperiodTPC="LHC10H8";
 674     if (reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11A10";
 675     fBeamType="PBPB";
 676   }
 677   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
 678   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
 679   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 680   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
 681   else if (fRun>=166529) {
 682     fLHCperiod="LHC11H";
 683     fMCperiodTPC="LHC11A10";
 684     fBeamType="PBPB";
 685   }
 686
 687
 688   //exception new pp MC productions from 2011
 689   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
 690   // exception for 11f1
 691   if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 692 }
 693
 694 //______________________________________________________________________________
 695 void AliPIDResponse::SetITSParametrisation()
 696 {
 697   //
 698   // Set the ITS parametrisation
 699   //
 700 }
 701
 702 //______________________________________________________________________________
 703 void AliPIDResponse::SetTPCPidResponseMaster()
 704 {
 705   //
 706   // Load the TPC pid response functions from the OADB
 707   //
 708   //don't load twice for the moment
 709    if (fArrPidResponseMaster) return;
 710
 711
 712   //reset the PID response functions
 713   delete fArrPidResponseMaster;
 714   fArrPidResponseMaster=0x0;
 715
 716   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
 717
 718   TFile *f=TFile::Open(fileName.Data());
 719   if (f && f->IsOpen() && !f->IsZombie()){
 720     fArrPidResponseMaster=dynamic_cast<TObjArray*>(f->Get("TPCPIDResponse"));
 721   }
 722   delete f;
 723
 724   if (!fArrPidResponseMaster){
 725     AliFatal(Form("Could not retrieve the TPC pid response from: %s",fileName.Data()));
 726     return;
 727   }
 728   fArrPidResponseMaster->SetOwner();
 729 }
 730
 731 //______________________________________________________________________________
 732 void AliPIDResponse::SetTPCParametrisation()
 733 {
 734   //
 735   // Change BB parametrisation for current run
 736   //
 737
 738   if (fLHCperiod.IsNull()) {
 739     AliFatal("No period set, not changing parametrisation");
 740     return;
 741   }
 742
 743   //
 744   // Set default parametrisations for data and MC
 745   //
 746
 747   //data type
 748   TString datatype="DATA";
 749   //in case of mc fRecoPass is per default 1
 750   if (fIsMC) {
 751     datatype="MC";
 752     fRecoPass=1;
 753   }
 754
 755   //
 756   //reset old splines
 757   //
 758   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 759     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
 760   }
 761
 762   // period
 763   TString period=fLHCperiod;
 764   if (fIsMC) period=fMCperiodTPC;
 765
 766   AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 767   Bool_t found=kFALSE;
 768   //
 769   //set the new PID splines
 770   //
 771   if (fArrPidResponseMaster){
 772     TObject *grAll=0x0;
 773     //for MC don't use period information
 774 //     if (fIsMC) period="[A-Z0-9]*";
 775     //for MC use MC period information
 776 //pattern for the default entry (valid for all particles)
 777     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 778
 779     //loop over entries and filter them
 780     for (Int_t iresp=0; iresp<fArrPidResponseMaster->GetEntriesFast();++iresp){
 781       TObject *responseFunction=fArrPidResponseMaster->At(iresp);
 782       if (responseFunction==0x0) continue;
 783       TString responseName=responseFunction->GetName();
 784
 785       if (!reg.MatchB(responseName)) continue;
 786
 787       TObjArray *arr=reg.MatchS(responseName);
 788       TString particleName=arr->At(1)->GetName();
 789       delete arr;
 790       if (particleName.IsNull()) continue;
 791       if (particleName=="ALL") grAll=responseFunction;
 792       else {
 793         //find particle id
 794         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 795           TString particle=AliPID::ParticleName(ispec);
 796           particle.ToUpper();
 797           if ( particle == particleName ){
 798             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
 799             fTPCResponse.SetUseDatabase(kTRUE);
 800             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
 801             found=kTRUE;
 802             break;
 803           }
 804         }
 805       }
 806     }
 807
 808     //set default response function to all particles which don't have a specific one
 809     if (grAll){
 810       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
 811         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
 812           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
 813           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
 814           found=kTRUE;
 815         }
 816       }
 817     }
 818   }
 819
 820   if (!found){
 821     AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 822   }
 823   //
 824   // Setup resolution parametrisation
 825   //
 826
 827   //default
 828   fTPCResponse.SetSigma(3.79301e-03, 2.21280e+04);
 829
 830   if (fRun>=122195){
 831     fTPCResponse.SetSigma(2.30176e-02, 5.60422e+02);
 832   }
 833   if (fArrPidResponseMaster)
 834   fResolutionCorrection=(TF1*)fArrPidResponseMaster->FindObject(Form("TF1_%s_ALL_%s_PASS%d_%s_SIGMA",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
 835
 836   if (fResolutionCorrection) AliInfo(Form("Setting multiplicity correction function: %s",fResolutionCorrection->GetName()));
 837 }
 838
 839 //______________________________________________________________________________
 840 void AliPIDResponse::SetTRDPidResponseMaster()
 841 {
 842   //
 843   // Load the TRD pid params and references from the OADB
 844   //
 845   if(fTRDPIDResponseObject) return;
 846   AliOADBContainer contParams("contParams");
 847
 848   Int_t statusResponse = contParams.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()), "AliTRDPIDResponseObject");
 849   if(statusResponse){
 850     AliError("Failed initializing PID Response Object from OADB");
 851   } else {
 852     AliInfo(Form("Loading TRD Response from %s/COMMON/PID/data/TRDPIDResponse.root", fOADBPath.Data()));
 853     fTRDPIDResponseObject = dynamic_cast<AliTRDPIDResponseObject *>(contParams.GetObject(fRun));
 854     if(!fTRDPIDResponseObject){
 855       AliError(Form("TRD Response not found in run %d", fRun));
 856     }
 857   }
 858   /*
 859   AliOADBContainer contRefs("contRefs");
 860   Int_t statusRefs = contRefs.InitFromFile(Form("%s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()), "AliTRDPIDReference");
 861   if(statusRefs){
 862     AliInfo("Failed Loading References for TRD");
 863   } else {
 864     AliInfo(Form("Loading TRD References from %s/COMMON/PID/data/TRDPIDReferenceLQ1D.root", fOADBPath.Data()));
 865     fTRDPIDReference = dynamic_cast<AliTRDPIDReference *>(contRefs.GetObject(fRun));
 866     if(!fTRDPIDReference){
 867       AliError(Form("TRD References not found in OADB Container for run %d", fRun));
 868     }
 869     }
 870     */
 871 }
 872
 873 //______________________________________________________________________________
 874 void AliPIDResponse::InitializeTRDResponse(){
 875   //
 876   // Set PID Params and references to the TRD PID response
 877   //
 878     fTRDResponse.SetPIDResponseObject(fTRDPIDResponseObject);
 879     if(fLHCperiod == "LHC10b" || fLHCperiod == "LHC10c" || fLHCperiod == "LHC10d" || fLHCperiod == "LHC10e"){
 880         fTRDslicesForPID[0] = 0;
 881         fTRDslicesForPID[1] = 7;
 882     }
 883 }
 884
 885 //______________________________________________________________________________
 886 void AliPIDResponse::SetTOFPidResponseMaster()
 887 {
 888   //
 889   // Load the TOF pid params from the OADB
 890   //
 891   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
 892   if (oadbf->IsOpen()) {
 893     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
 894     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
 895     if (fTOFPIDParams) delete fTOFPIDParams;
 896     fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
 897     oadbf->Close();
 898     delete oadbc;
 899   } else {
 900     AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
 901   }
 902   delete oadbf;
 903
 904   }
 905
 906 //______________________________________________________________________________
 907 void AliPIDResponse::InitializeTOFResponse(){
 908   //
 909   // Set PID Params to the TOF PID response
 910   //
 911   for (Int_t i=0;i<4;i++) {
 912     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
 913   }
 914   fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
 915
 916 }
 917
 918
 919 //_________________________________________________________________________
 920 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
 921   //
 922   // Check whether track is identified as electron under a given electron efficiency hypothesis
 923   //
 924   Double_t probs[AliPID::kSPECIES];
 925   ComputeTRDProbability(vtrack, AliPID::kSPECIES, probs);
 926
 927   Int_t ntracklets = vtrack->GetTRDntrackletsPID();
 928   // Take mean of the TRD momenta in the given tracklets
 929   Float_t p = 0, trdmomenta[AliVTrack::kTRDnPlanes];
 930   Int_t nmomenta = 0;
 931   for(Int_t iPl=0;iPl<AliVTrack::kTRDnPlanes;iPl++){
 932     if(vtrack->GetTRDmomentum(iPl) > 0.){
 933       trdmomenta[nmomenta++] = vtrack->GetTRDmomentum(iPl);
 934     }
 935   }
 936   p = TMath::Mean(nmomenta, trdmomenta);
 937
 938   return fTRDResponse.IdentifiedAsElectron(ntracklets, probs, p, efficiencyLevel);
 939 }
 940
 941 //______________________________________________________________________________
 942 void AliPIDResponse::SetEMCALPidResponseMaster()
 943 {
 944   //
 945   // Load the EMCAL pid response functions from the OADB
 946   //
 947   TObjArray* fEMCALPIDParamsRun      = NULL;
 948   TObjArray* fEMCALPIDParamsPass     = NULL;
 949
 950   if(fEMCALPIDParams) return;
 951   AliOADBContainer contParams("contParams");
 952
 953   Int_t statusPars = contParams.InitFromFile(Form("%s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()), "AliEMCALPIDParams");
 954   if(statusPars){
 955     AliError("Failed initializing PID Params from OADB");
 956   }
 957   else {
 958     AliInfo(Form("Loading EMCAL Params from %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 959
 960     fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(fRun));
 961     if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",fRecoPass)));
 962     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 963
 964     if(!fEMCALPIDParams){
 965       AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
 966       AliInfo("Will take the standard LHC11d instead ...");
 967
 968       fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
 969       if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
 970       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 971
 972       if(!fEMCALPIDParams){
 973         AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));
 974       }
 975     }
 976   }
 977 }
 978
 979 //______________________________________________________________________________
 980 void AliPIDResponse::InitializeEMCALResponse(){
 981   //
 982   // Set PID Params to the EMCAL PID response
 983   //
 984   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
 985
 986 }
 987 //_________________________________________________________________________
 988 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
 989   //
 990   // Set TOF response function
 991   // Input option for event_time used
 992   //
 993
 994     Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
 995     if(t0spread < 10) t0spread = 80;
 996
 997     // T0 from TOF algorithm
 998
 999     Bool_t flagT0TOF=kFALSE;
1000     Bool_t flagT0T0=kFALSE;
1001     Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
1002     Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
1003     Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
1004
1005     // T0-TOF arrays
1006     Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
1007     Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
1008     for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1009       estimatedT0event[i]=0.0;
1010       estimatedT0resolution[i]=0.0;
1011       startTimeMask[i] = 0;
1012     }
1013
1014     Float_t resT0A=75,resT0C=65,resT0AC=55;
1015     if(vevent->GetT0TOF()){ // check if T0 detector information is available
1016         flagT0T0=kTRUE;
1017     }
1018
1019
1020     AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
1021
1022     if (tofHeader) { // read global info and T0-TOF
1023       fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
1024       t0spread = tofHeader->GetT0spread(); // read t0 sprad
1025       if(t0spread < 10) t0spread = 80;
1026
1027       flagT0TOF=kTRUE;
1028       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
1029         startTime[i]=tofHeader->GetDefaultEventTimeVal();
1030         startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
1031         if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
1032       }
1033
1034       TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
1035       TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
1036       TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
1037       for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
1038         Int_t icurrent = (Int_t)ibin->GetAt(j);
1039         startTime[icurrent]=t0Bin->GetAt(j);
1040         startTimeRes[icurrent]=t0ResBin->GetAt(j);
1041         if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
1042       }
1043     }
1044
1045     // for cut of 3 sigma on t0 spread
1046     Float_t t0cut = 3 * t0spread;
1047     if(t0cut < 500) t0cut = 500;
1048
1049     if(option == kFILL_T0){ // T0-FILL is used
1050         for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1051           estimatedT0event[i]=0.0;
1052           estimatedT0resolution[i]=t0spread;
1053         }
1054         fTOFResponse.SetT0event(estimatedT0event);
1055         fTOFResponse.SetT0resolution(estimatedT0resolution);
1056     }
1057
1058     if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
1059         if(flagT0TOF){
1060             fTOFResponse.SetT0event(startTime);
1061             fTOFResponse.SetT0resolution(startTimeRes);
1062             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1063               if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1064               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1065             }
1066         }
1067         else{
1068             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1069               estimatedT0event[i]=0.0;
1070               estimatedT0resolution[i]=t0spread;
1071               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1072             }
1073             fTOFResponse.SetT0event(estimatedT0event);
1074             fTOFResponse.SetT0resolution(estimatedT0resolution);
1075         }
1076     }
1077     else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
1078         Float_t t0AC=-10000;
1079         Float_t t0A=-10000;
1080         Float_t t0C=-10000;
1081         if(flagT0T0){
1082             t0AC= vevent->GetT0TOF()[0];
1083             t0A= vevent->GetT0TOF()[1];
1084             t0C= vevent->GetT0TOF()[2];
1085         }
1086
1087         Float_t t0t0Best = 0;
1088         Float_t t0t0BestRes = 9999;
1089         Int_t t0used=0;
1090         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1091             t0t0Best = t0AC;
1092             t0t0BestRes = resT0AC;
1093             t0used=6;
1094         }
1095         else if(TMath::Abs(t0C) < t0cut){
1096             t0t0Best = t0C;
1097             t0t0BestRes = resT0C;
1098             t0used=4;
1099         }
1100         else if(TMath::Abs(t0A) < t0cut){
1101             t0t0Best = t0A;
1102             t0t0BestRes = resT0A;
1103             t0used=2;
1104         }
1105
1106         if(flagT0TOF){ // if T0-TOF info is available
1107             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1108                 if(t0t0BestRes < 999){
1109                   if(startTimeRes[i] < t0spread){
1110                     Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
1111                     Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
1112                     estimatedT0event[i]=t0best / wtot;
1113                     estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
1114                     startTimeMask[i] = t0used+1;
1115                   }
1116                   else {
1117                     estimatedT0event[i]=t0t0Best;
1118                     estimatedT0resolution[i]=t0t0BestRes;
1119                     startTimeMask[i] = t0used;
1120                   }
1121                 }
1122                 else{
1123                   estimatedT0event[i]=startTime[i];
1124                   estimatedT0resolution[i]=startTimeRes[i];
1125                   if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
1126                 }
1127                 fTOFResponse.SetT0binMask(i,startTimeMask[i]);
1128             }
1129             fTOFResponse.SetT0event(estimatedT0event);
1130             fTOFResponse.SetT0resolution(estimatedT0resolution);
1131         }
1132         else{ // if no T0-TOF info is available
1133             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1134               fTOFResponse.SetT0binMask(i,t0used);
1135               if(t0t0BestRes < 999){
1136                 estimatedT0event[i]=t0t0Best;
1137                 estimatedT0resolution[i]=t0t0BestRes;
1138               }
1139               else{
1140                 estimatedT0event[i]=0.0;
1141                 estimatedT0resolution[i]=t0spread;
1142               }
1143             }
1144             fTOFResponse.SetT0event(estimatedT0event);
1145             fTOFResponse.SetT0resolution(estimatedT0resolution);
1146         }
1147     }
1148
1149     else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
1150         Float_t t0AC=-10000;
1151         Float_t t0A=-10000;
1152         Float_t t0C=-10000;
1153         if(flagT0T0){
1154             t0AC= vevent->GetT0TOF()[0];
1155             t0A= vevent->GetT0TOF()[1];
1156             t0C= vevent->GetT0TOF()[2];
1157         }
1158
1159         if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
1160             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1161               estimatedT0event[i]=t0AC;
1162               estimatedT0resolution[i]=resT0AC;
1163               fTOFResponse.SetT0binMask(i,6);
1164             }
1165         }
1166         else if(TMath::Abs(t0C) < t0cut){
1167             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1168               estimatedT0event[i]=t0C;
1169               estimatedT0resolution[i]=resT0C;
1170               fTOFResponse.SetT0binMask(i,4);
1171             }
1172         }
1173         else if(TMath::Abs(t0A) < t0cut){
1174             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1175               estimatedT0event[i]=t0A;
1176               estimatedT0resolution[i]=resT0A;
1177               fTOFResponse.SetT0binMask(i,2);
1178             }
1179         }
1180         else{
1181             for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
1182               estimatedT0event[i]=0.0;
1183               estimatedT0resolution[i]=t0spread;
1184               fTOFResponse.SetT0binMask(i,0);
1185             }
1186         }
1187         fTOFResponse.SetT0event(estimatedT0event);
1188         fTOFResponse.SetT0resolution(estimatedT0resolution);
1189     }
1190     delete [] startTime;
1191     delete [] startTimeRes;
1192     delete [] startTimeMask;
1193     delete [] estimatedT0event;
1194     delete [] estimatedT0resolution;
1195 }