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