]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - STEER/STEERBase/AliPIDResponse.cxx
Changes for #95543: request to commit to trunk: TOF code supporting PID for light...
[u/mrichter/AliRoot.git] / STEER / STEERBase / AliPIDResponse.cxx
index d5ac77cc420453ccd076aadcd5d75fbf03460213..511475dccb91a9926f94bdaceba7524f7eda6346 100644 (file)
 #include <AliOADBContainer.h>
 #include <AliTRDPIDParams.h>
 #include <AliTRDPIDReference.h>
+#include <AliTOFPIDParams.h>
 
 #include "AliPIDResponse.h"
 
+#include "AliCentrality.h"
+
 ClassImp(AliPIDResponse);
 
 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
@@ -66,17 +69,18 @@ fArrPidResponseMaster(0x0),
 fResolutionCorrection(0x0),
 fTRDPIDParams(0x0),
 fTRDPIDReference(0x0),
-fTOFTimeZeroType(kBest_T0),
-fTOFres(100.),
+fTOFtail(1.1),
+fTOFPIDParams(0x0),
 fEMCALPIDParams(0x0),
-fCurrentEvent(0x0)
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
 {
   //
   // default ctor
   //
-  AliLog::SetClassDebugLevel("AliPIDResponse",10);
-  AliLog::SetClassDebugLevel("AliESDpid",10);
-  AliLog::SetClassDebugLevel("AliAODpidUtil",10);
+  AliLog::SetClassDebugLevel("AliPIDResponse",0);
+  AliLog::SetClassDebugLevel("AliESDpid",0);
+  AliLog::SetClassDebugLevel("AliAODpidUtil",0);
 
   memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
 }
@@ -90,6 +94,7 @@ AliPIDResponse::~AliPIDResponse()
   delete fArrPidResponseMaster;
   delete fTRDPIDParams;
   delete fTRDPIDReference;
+  if (fTOFPIDParams) delete fTOFPIDParams;
 }
 
 //______________________________________________________________________________
@@ -117,10 +122,11 @@ fArrPidResponseMaster(0x0),
 fResolutionCorrection(0x0),
 fTRDPIDParams(0x0),
 fTRDPIDReference(0x0),
-fTOFTimeZeroType(AliPIDResponse::kBest_T0),
-fTOFres(100.),
+fTOFtail(1.1),
+fTOFPIDParams(0x0),
 fEMCALPIDParams(0x0),
-fCurrentEvent(0x0)
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
 {
   //
   // copy ctor
@@ -161,8 +167,8 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fTRDPIDReference=0x0;
     fEMCALPIDParams=0x0;
     memset(fTRDslicesForPID,0,sizeof(UInt_t)*2);
-    fTOFTimeZeroType=AliPIDResponse::kBest_T0;
-    fTOFres=100.;
+    fTOFtail=1.1;
+    fTOFPIDParams=0x0;
     fCurrentEvent=other.fCurrentEvent;
   }
   return *this;
@@ -204,7 +210,7 @@ Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EPar
   // Track matching
   nMatchClus = track->GetEMCALcluster();
   if(nMatchClus > -1){
-
+    
     mom    = track->P();
     pt     = track->Pt();
     charge = track->Charge();
@@ -213,23 +219,70 @@ Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EPar
     
     if(matchedClus){
       
-    // matched cluster is EMCAL
-    if(matchedClus->IsEMCAL()){
-      
-      fClsE       = matchedClus->E();
-      EovP        = fClsE/mom;
-      
-      
-      // NSigma value really meaningful only for electrons!
-      return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
+      // matched cluster is EMCAL
+      if(matchedClus->IsEMCAL()){
+       
+       fClsE       = matchedClus->E();
+       EovP        = fClsE/mom;
+       
+       
+       // NSigma value really meaningful only for electrons!
+       return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
+      }
     }
   }
-  }
-
+  
   return -999;
   
 }
 
+//______________________________________________________________________________
+Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
+
+  AliVCluster *matchedClus = NULL;
+
+  Double_t mom     = -1.; 
+  Double_t pt      = -1.; 
+  Double_t EovP    = -1.;
+  Double_t fClsE   = -1.;
+  
+  Int_t nMatchClus = -1;
+  Int_t charge     = 0;
+  
+  // Track matching
+  nMatchClus = track->GetEMCALcluster();
+  if(nMatchClus > -1){
+
+    mom    = track->P();
+    pt     = track->Pt();
+    charge = track->Charge();
+    
+    matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
+    
+    if(matchedClus){
+      
+      // matched cluster is EMCAL
+      if(matchedClus->IsEMCAL()){
+       
+       fClsE       = matchedClus->E();
+       EovP        = fClsE/mom;
+       
+       // fill used EMCAL variables here
+       eop            = EovP; // E/p
+       showershape[0] = matchedClus->GetNCells(); // number of cells in cluster
+       showershape[1] = matchedClus->GetM02(); // long axis
+       showershape[2] = matchedClus->GetM20(); // short axis
+       showershape[3] = matchedClus->GetDispersion(); // dispersion
+       
+       // NSigma value really meaningful only for electrons!
+       return fEMCALResponse.GetNumberOfSigmas(pt,EovP,type,charge); 
+      }
+    }
+  }
+  return -999;
+}
+
+
 //______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 {
@@ -354,6 +407,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliV
   // Compute PID response for the
   //
 
+  Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
+
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
   
@@ -363,21 +418,37 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliV
   Double_t time[AliPID::kSPECIESN];
   track->GetIntegratedTimes(time);
   
-  Double_t sigma[AliPID::kSPECIES];
-  for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
-    sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
-  }
+  //  Double_t sigma[AliPID::kSPECIES];
+  // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
+  //    sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
+  //  }
   
   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
     AliPID::EParticleType type=AliPID::EParticleType(j);
-    Double_t nsigmas=NumberOfSigmasTOF(track,type);
-    
-    Double_t sig = sigma[j];
+    Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
+
+    //    Double_t sig = sigma[j];
+    Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
+    Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
+    if (TMath::Abs(nsigmas) > (fRange+2)) {
+      if(nsigmas < fTOFtail)
+       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
+      else
+       p[j] = TMath::Exp(-(fRange+2 - fTOFtail*0.5)*fTOFtail)/sig;
+    } else{
+      if(nsigmas < fTOFtail)
+       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+      else
+       p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
+    }
+
+    /* OLD Gaussian shape
     if (TMath::Abs(nsigmas) > (fRange+2)) {
       p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
     } else
       p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
+    */
 
     if (TMath::Abs(nsigmas)<5.){
       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
@@ -522,8 +593,17 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
   }
   
   //TOF resolution
-  SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
-  
+  SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFPIDParams->GetStartTimeMethod());
+
+
+  // Get and set centrality
+  AliCentrality *centrality = event->GetCentrality();
+  if(centrality){
+    fCurrCentrality = centrality->GetCentralityPercentile("V0M");
+  }
+  else{
+    fCurrCentrality = -1;
+  }
 }
 
 //______________________________________________________________________________
@@ -545,7 +625,8 @@ void AliPIDResponse::ExecNewRun()
   SetEMCALPidResponseMaster(); 
   InitializeEMCALResponse();
   
-  fTOFResponse.SetTimeResolution(fTOFres);
+  SetTOFPidResponseMaster();
+  InitializeTOFResponse();
 }
 
 //_____________________________________________________
@@ -581,7 +662,9 @@ void AliPIDResponse::SetRecoInfo()
     
   fBeamType="PP";
   
-  TPRegexp reg(".*(LHC11[a-z]+[0-9]+[a-z_]*)/.*");
+
+  TPRegexp reg(".*(LHC1[1-2][a-z]+[0-9]+[a-z_]*)/.*");
+
   //find the period by run number (UGLY, but not stored in ESD and AOD... )
   if (fRun>=114737&&fRun<=117223)      { fLHCperiod="LHC10B"; fMCperiodTPC="LHC10D1";  }
   else if (fRun>=118503&&fRun<=121040) { fLHCperiod="LHC10C"; fMCperiodTPC="LHC10D1";  }
@@ -596,6 +679,9 @@ void AliPIDResponse::SetRecoInfo()
     fBeamType="PBPB";
   }
   else if (fRun>=139699&&fRun<=146860) { fLHCperiod="LHC11A"; fMCperiodTPC="LHC10F6A"; }
+  //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
+  else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+  else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
   else if (fRun>=166529) {
     fLHCperiod="LHC11H";
     fMCperiodTPC="LHC11A10";
@@ -605,6 +691,8 @@ void AliPIDResponse::SetRecoInfo()
 
   //exception new pp MC productions from 2011
   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
+  // exception for 11f1
+  if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";                                                                                                                
 }
 
 //______________________________________________________________________________
@@ -674,17 +762,21 @@ void AliPIDResponse::SetTPCParametrisation()
   for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
   }
-  
+
+  // period
+  TString period=fLHCperiod;
+  if (fIsMC) period=fMCperiodTPC;
+
+  AliInfo(Form("Searching splines for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+  Bool_t found=kFALSE;
   //
   //set the new PID splines
   //
-  TString period=fLHCperiod;
   if (fArrPidResponseMaster){
     TObject *grAll=0x0;
     //for MC don't use period information
 //     if (fIsMC) period="[A-Z0-9]*";
     //for MC use MC period information
-    if (fIsMC) period=fMCperiodTPC;
 //pattern for the default entry (valid for all particles)
     TPRegexp reg(Form("TSPLINE3_%s_([A-Z]*)_%s_PASS%d_%s_MEAN",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
     
@@ -710,6 +802,7 @@ void AliPIDResponse::SetTPCParametrisation()
             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
             fTPCResponse.SetUseDatabase(kTRUE);
             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
+            found=kTRUE;
             break;
           }
         }
@@ -722,11 +815,15 @@ void AliPIDResponse::SetTPCParametrisation()
         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
+          found=kTRUE;
         }
       }
     }
   }
-  
+
+  if (!found){
+    AliError(Form("No splines found for: %s %s PASS%d %s",datatype.Data(),period.Data(),fRecoPass,fBeamType.Data()));
+  }
   //
   // Setup resolution parametrisation
   //
@@ -789,6 +886,40 @@ void AliPIDResponse::InitializeTRDResponse(){
   }
 }
 
+//______________________________________________________________________________
+void AliPIDResponse::SetTOFPidResponseMaster()
+{
+  //
+  // Load the TOF pid params from the OADB
+  //
+  TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
+  if (oadbf->IsOpen()) {
+    AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
+    AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
+    if (fTOFPIDParams) delete fTOFPIDParams;
+    fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
+    oadbf->Close();
+    delete oadbc;
+  } else {
+    AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
+  }
+  delete oadbf;
+
+  }
+
+//______________________________________________________________________________
+void AliPIDResponse::InitializeTOFResponse(){
+  //
+  // Set PID Params to the TOF PID response
+  // 
+  for (Int_t i=0;i<4;i++) {
+    fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
+  }
+  fTOFResponse.SetTimeResolution(fTOFPIDParams->GetTOFresolution());
+
+}
+
+
 //_________________________________________________________________________
 Bool_t AliPIDResponse::IdentifiedAsElectronTRD(const AliVTrack *vtrack, Double_t efficiencyLevel) const {
   //
@@ -835,15 +966,15 @@ void AliPIDResponse::SetEMCALPidResponseMaster()
     if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 
     if(!fEMCALPIDParams){
-      AliError(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
-      AliInfo("Will take the standard LHC11a pass2 instead ...");
+      AliInfo(Form("EMCAL Params not found in run %d pass %d", fRun, fRecoPass));
+      AliInfo("Will take the standard LHC11d instead ...");
 
-      fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(144871));
-      if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",2)));
+      fEMCALPIDParamsRun = dynamic_cast<TObjArray *>(contParams.GetObject(156477));
+      if(fEMCALPIDParamsRun)  fEMCALPIDParamsPass = dynamic_cast<TObjArray *>(fEMCALPIDParamsRun->FindObject(Form("pass%d",1)));
       if(fEMCALPIDParamsPass) fEMCALPIDParams     = dynamic_cast<TObjArray *>(fEMCALPIDParamsPass->FindObject(Form("EMCALPIDParams_Particles")));
 
       if(!fEMCALPIDParams){
-       AliError(Form("DEFAULT EMCAL Params (LHC11a pass2) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));       
+       AliError(Form("DEFAULT EMCAL Params (LHC11d) not found in file %s/COMMON/PID/data/EMCALPIDParams.root", fOADBPath.Data()));     
       }
     }
   }
@@ -857,3 +988,212 @@ void AliPIDResponse::InitializeEMCALResponse(){
   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
 
 }
+//_________________________________________________________________________
+void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
+  //
+  // Set TOF response function
+  // Input option for event_time used
+  //
+  
+    Float_t t0spread = 0.; //vevent->GetEventTimeSpread();
+    if(t0spread < 10) t0spread = 80;
+
+    // T0 from TOF algorithm
+
+    Bool_t flagT0TOF=kFALSE;
+    Bool_t flagT0T0=kFALSE;
+    Float_t *startTime = new Float_t[fTOFResponse.GetNmomBins()];
+    Float_t *startTimeRes = new Float_t[fTOFResponse.GetNmomBins()];
+    Int_t *startTimeMask = new Int_t[fTOFResponse.GetNmomBins()];
+
+    // T0-TOF arrays
+    Float_t *estimatedT0event = new Float_t[fTOFResponse.GetNmomBins()];
+    Float_t *estimatedT0resolution = new Float_t[fTOFResponse.GetNmomBins()];
+    for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+      estimatedT0event[i]=0.0;
+      estimatedT0resolution[i]=0.0;
+      startTimeMask[i] = 0;
+    }
+
+    Float_t resT0A=75,resT0C=65,resT0AC=55;
+    if(vevent->GetT0TOF()){ // check if T0 detector information is available
+       flagT0T0=kTRUE;
+    }
+
+
+    AliTOFHeader *tofHeader = (AliTOFHeader*)vevent->GetTOFHeader();
+
+    if (tofHeader) { // read global info and T0-TOF
+      fTOFResponse.SetTimeResolution(tofHeader->GetTOFResolution());
+      t0spread = tofHeader->GetT0spread(); // read t0 sprad
+      if(t0spread < 10) t0spread = 80;
+
+      flagT0TOF=kTRUE;
+      for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){ // read T0-TOF default value
+       startTime[i]=tofHeader->GetDefaultEventTimeVal();
+       startTimeRes[i]=tofHeader->GetDefaultEventTimeRes();
+       if(startTimeRes[i] < 1.e-5) startTimeRes[i] = t0spread;
+      }
+
+      TArrayI *ibin=(TArrayI*)tofHeader->GetNvalues();
+      TArrayF *t0Bin=(TArrayF*)tofHeader->GetEventTimeValues();
+      TArrayF *t0ResBin=(TArrayF*)tofHeader->GetEventTimeRes();
+      for(Int_t j=0;j < tofHeader->GetNbins();j++){ // fill T0-TOF in p-bins
+       Int_t icurrent = (Int_t)ibin->GetAt(j);
+       startTime[icurrent]=t0Bin->GetAt(j);
+       startTimeRes[icurrent]=t0ResBin->GetAt(j);
+       if(startTimeRes[icurrent] < 1.e-5) startTimeRes[icurrent] = t0spread;
+      }
+    }
+
+    // for cut of 3 sigma on t0 spread
+    Float_t t0cut = 3 * t0spread;
+    if(t0cut < 500) t0cut = 500;
+
+    if(option == kFILL_T0){ // T0-FILL is used
+       for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+         estimatedT0event[i]=0.0;
+         estimatedT0resolution[i]=t0spread;
+       }
+       fTOFResponse.SetT0event(estimatedT0event);
+       fTOFResponse.SetT0resolution(estimatedT0resolution);
+    }
+
+    if(option == kTOF_T0){ // T0-TOF is used when available (T0-FILL otherwise) from ESD
+       if(flagT0TOF){
+           fTOFResponse.SetT0event(startTime);
+           fTOFResponse.SetT0resolution(startTimeRes);
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
+             fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+           }
+       }
+       else{
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             estimatedT0event[i]=0.0;
+             estimatedT0resolution[i]=t0spread;
+             fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+           }
+           fTOFResponse.SetT0event(estimatedT0event);
+           fTOFResponse.SetT0resolution(estimatedT0resolution);
+       }
+    }
+    else if(option == kBest_T0){ // T0-T0 or T0-TOF are used when available (T0-FILL otherwise) from ESD
+       Float_t t0AC=-10000;
+       Float_t t0A=-10000;
+       Float_t t0C=-10000;
+       if(flagT0T0){
+           t0AC= vevent->GetT0TOF()[0];
+           t0A= vevent->GetT0TOF()[1];
+           t0C= vevent->GetT0TOF()[2];
+       }
+
+       Float_t t0t0Best = 0;
+       Float_t t0t0BestRes = 9999;
+       Int_t t0used=0;
+       if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
+           t0t0Best = t0AC;
+           t0t0BestRes = resT0AC;
+           t0used=6;
+       }
+       else if(TMath::Abs(t0C) < t0cut){
+           t0t0Best = t0C;
+           t0t0BestRes = resT0C;
+           t0used=4;
+       }
+       else if(TMath::Abs(t0A) < t0cut){
+           t0t0Best = t0A;
+           t0t0BestRes = resT0A;
+           t0used=2;
+       }
+
+       if(flagT0TOF){ // if T0-TOF info is available
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+               if(t0t0BestRes < 999){
+                 if(startTimeRes[i] < t0spread){
+                   Double_t wtot = 1./startTimeRes[i]/startTimeRes[i] + 1./t0t0BestRes/t0t0BestRes;
+                   Double_t t0best = startTime[i]/startTimeRes[i]/startTimeRes[i] + t0t0Best/t0t0BestRes/t0t0BestRes;
+                   estimatedT0event[i]=t0best / wtot;
+                   estimatedT0resolution[i]=1./TMath::Sqrt(wtot);
+                   startTimeMask[i] = t0used+1;
+                 }
+                 else {
+                   estimatedT0event[i]=t0t0Best;
+                   estimatedT0resolution[i]=t0t0BestRes;
+                   startTimeMask[i] = t0used;
+                 }
+               }
+               else{
+                 estimatedT0event[i]=startTime[i];
+                 estimatedT0resolution[i]=startTimeRes[i];
+                 if(startTimeRes[i]<t0spread) startTimeMask[i]=1;
+               }
+               fTOFResponse.SetT0binMask(i,startTimeMask[i]);
+           }
+           fTOFResponse.SetT0event(estimatedT0event);
+           fTOFResponse.SetT0resolution(estimatedT0resolution);
+       }
+       else{ // if no T0-TOF info is available
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             fTOFResponse.SetT0binMask(i,t0used);
+             if(t0t0BestRes < 999){
+               estimatedT0event[i]=t0t0Best;
+               estimatedT0resolution[i]=t0t0BestRes;
+             }
+             else{
+               estimatedT0event[i]=0.0;
+               estimatedT0resolution[i]=t0spread;
+             }
+           }
+           fTOFResponse.SetT0event(estimatedT0event);
+           fTOFResponse.SetT0resolution(estimatedT0resolution);
+       }
+    }
+
+    else if(option == kT0_T0){ // T0-T0 is used when available (T0-FILL otherwise)
+       Float_t t0AC=-10000;
+       Float_t t0A=-10000;
+       Float_t t0C=-10000;
+       if(flagT0T0){
+           t0AC= vevent->GetT0TOF()[0];
+           t0A= vevent->GetT0TOF()[1];
+           t0C= vevent->GetT0TOF()[2];
+       }
+
+       if(TMath::Abs(t0A) < t0cut && TMath::Abs(t0C) < t0cut && TMath::Abs(t0C-t0A) < 500){
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             estimatedT0event[i]=t0AC;
+             estimatedT0resolution[i]=resT0AC;
+             fTOFResponse.SetT0binMask(i,6);
+           }
+       }
+       else if(TMath::Abs(t0C) < t0cut){
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             estimatedT0event[i]=t0C;
+             estimatedT0resolution[i]=resT0C;
+             fTOFResponse.SetT0binMask(i,4);
+           }
+       }
+       else if(TMath::Abs(t0A) < t0cut){
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             estimatedT0event[i]=t0A;
+             estimatedT0resolution[i]=resT0A;
+             fTOFResponse.SetT0binMask(i,2);
+           }
+       }
+       else{
+           for(Int_t i=0;i<fTOFResponse.GetNmomBins();i++){
+             estimatedT0event[i]=0.0;
+             estimatedT0resolution[i]=t0spread;
+             fTOFResponse.SetT0binMask(i,0);
+           }
+       }
+       fTOFResponse.SetT0event(estimatedT0event);
+       fTOFResponse.SetT0resolution(estimatedT0resolution);
+    }
+    delete [] startTime;
+    delete [] startTimeRes;
+    delete [] startTimeMask;
+    delete [] estimatedT0event;
+    delete [] estimatedT0resolution;
+}