Timing hgcal realistic sim pf

RecoParticleFlow/PFClusterProducer/plugins/Cluster3DPCACalculator.cc

@@ -17,8 +17,9 @@
 class Cluster3DPCACalculator : public PFCPositionCalculatorBase {
  public:
   Cluster3DPCACalculator(const edm::ParameterSet& conf) :
-    PFCPositionCalculatorBase(conf),    
-    pca_(new TPrincipal(3,"D")){ 
+    PFCPositionCalculatorBase(conf),
+    updateTiming_(conf.getParameter<bool>("updateTiming")),
+    pca_(new TPrincipal(3,"D")){
   }
   Cluster3DPCACalculator(const Cluster3DPCACalculator&) = delete;
   Cluster3DPCACalculator& operator=(const Cluster3DPCACalculator&) = delete;
@@ -27,9 +28,10 @@ class Cluster3DPCACalculator : public PFCPositionCalculatorBase {
   void calculateAndSetPositions(reco::PFClusterCollection&) override;
 
  private:
+  const bool updateTiming_;
   std::unique_ptr<TPrincipal> pca_;
 
-  void showerParameters(const reco::PFCluster&, math::XYZPoint&, 
+  void showerParameters(const reco::PFCluster&, math::XYZPoint&,
 			math::XYZVector& );
 
   void calculateAndSetPositionActual(reco::PFCluster&);
@@ -54,18 +56,18 @@ calculateAndSetPositions(reco::PFClusterCollection& clusters) {
 }
 
 void Cluster3DPCACalculator::
-calculateAndSetPositionActual(reco::PFCluster& cluster) {  
+calculateAndSetPositionActual(reco::PFCluster& cluster) {
   if( !cluster.seed() ) {
     throw cms::Exception("ClusterWithNoSeed")
       << " Found a cluster with no seed: " << cluster;
-  }  				
-  double cl_energy = 0;  
+  }
+  double cl_energy = 0;
   double max_e = 0.0;
   double avg_time = 0.0;
   double time_norm = 0.0;
   PFLayer::Layer max_e_layer = PFLayer::NONE;
-  reco::PFRecHitRef refseed;  
-  double pcavars[3];  
+  reco::PFRecHitRef refseed;
+  double pcavars[3];
 
   for( const reco::PFRecHitFraction& rhf : cluster.recHitFractions() ) {
     const reco::PFRecHitRef& refhit = rhf.recHitRef();
@@ -76,33 +78,33 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) {
       // all times are offset by one nanosecond in digitizer
       // remove that here so all times of flight
       // are with respect to (0,0,0)
-      avg_time += (rh_time - 1.0); 
+      avg_time += (rh_time - 1.0);
       time_norm += 1.0;
     }
     if( rh_energy > max_e ) {
       max_e = rh_energy;
       max_e_layer = rhf.recHitRef()->layer();
-    }  
+    }
     if( refhit->detId() == cluster.seed() ) refseed = refhit;
     const double rh_fraction = rhf.fraction();
     rh_energy = refhit->energy()*rh_fraction;
     if( edm::isNotFinite(rh_energy) ) {
 //temporarily changed exception to warning
 //      throw cms::Exception("PFClusterAlgo")
       edm::LogWarning("PFClusterAlgo")
-      <<"rechit " << refhit->detId() << " has a NaN energy... " 
+      <<"rechit " << refhit->detId() << " has a NaN energy... "
       << "The input of the particle flow clustering seems to be corrupted.";
       continue;
-    }    
+    }
     pcavars[0] = refhit->position().x();
     pcavars[1] = refhit->position().y();
-    pcavars[2] = refhit->position().z();     
+    pcavars[2] = refhit->position().z();
     int nhit = int( rh_energy*100 ); // put rec_hit energy in units of 10 MeV
 
     for( int i = 0; i < nhit; ++i ) {
       pca_->AddRow(pcavars);
     }
-      
+
   }
   cluster.setEnergy(cl_energy);
   cluster.setLayer(max_e_layer);
@@ -111,7 +113,7 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) {
   pca_->MakePrincipals();
   const TVectorD& means = *(pca_->GetMeanValues());
   const TMatrixD& eigens = *(pca_->GetEigenVectors());
-  
+
   math::XYZPoint  barycenter(means[0],means[1],means[2]);
   math::XYZVector axis(eigens(0,0),eigens(1,0),eigens(2,0));
 
@@ -124,8 +126,10 @@ calculateAndSetPositionActual(reco::PFCluster& cluster) {
   if( axis.z()*barycenter.z() < 0.0 ) {
     axis = math::XYZVector(-eigens(0,0),-eigens(1,0),-eigens(2,0));
   }
-
-  cluster.setTime(avg_time);
+
+  if (updateTiming_) {
+    cluster.setTime(avg_time);
+  }
   cluster.setPosition(barycenter);
   cluster.calculatePositionREP();
 

RecoParticleFlow/PFClusterProducer/plugins/SimMappers/ComputeClusterTime.h

@@ -0,0 +1,128 @@
+#ifndef RecoParticleFlow_PFClusterProducer_ComputeClusterTime_h
+#define RecoParticleFlow_PFClusterProducer_ComputeClusterTime_h
+
+// user include files
+#include <vector>
+
+
+// functions to select the hits to compute the time of a given cluster
+// start with the only hits with timing information
+// average among the hits contained in the chosen time interval
+
+// N.B. time is corrected wrt vtx-calorimeter distance 
+// with straight line and light speed hypothesis
+// for charged tracks or heavy particles (longer track length or beta < 1) 
+// need to correct the offset at analysis level
+
+
+
+namespace hgcalsimclustertime {
+
+  //time-interval based on that ~210ps wide and with the highest number of hits
+  float fixSizeHighestDensity(std::vector<float>& t, float deltaT=0.210 /*time window in ns*/, float timeWidthBy=0.5){
+
+    float tolerance = 0.05f;
+    std::sort(t.begin(), t.end());
+
+    int max_elements = 0;
+    int start_el = 0;
+    int end_el = 0;
+    float timeW = 0.f;
+
+    for(auto start = t.begin(); start != t.end(); ++start) {
+      const auto startRef = *start;
+      int c = count_if(start, t.end(), [&](float el) {
+	  return el - startRef <= deltaT + tolerance;
+	});
+      if (c > max_elements) {
+	max_elements = c;
+	auto last_el = find_if_not(start, t.end(), [&](float el) {
+	    return el - startRef <= deltaT + tolerance;
+	  });
+	auto val = *(--last_el);
+	if (std::abs(deltaT - (val - startRef)) < tolerance) {
+	  tolerance = std::abs(deltaT - (val - startRef));
+	}
+	start_el = distance(t.begin(), start);
+	end_el = distance(t.begin(), last_el);
+	timeW = val - startRef;
+      }
+    }
+
+    // further adjust time width around the chosen one based on the hits density
+    // proved to improve the resolution: get as many hits as possible provided they are close in time
+
+    float HalfTimeDiff = timeW * timeWidthBy;
+    float sum = 0.;
+    int num = 0;
+    int totSize = t.size();
+
+    for(int ij=0; ij<=start_el; ++ij){
+      if(t[ij] > (t[start_el] - HalfTimeDiff) ){
+	for(int kl=ij; kl<totSize; ++kl){
+	  if(t[kl] < (t[end_el] + HalfTimeDiff) ){
+	    sum += t[kl];
+	    ++num;
+	  }
+	  else  break;
+	}
+	break;
+      }
+    }
+
+    if(num == 0) return -99.;
+    return sum/num;
+  }
+
+
+  //useful for future developments - baseline for 0PU
+  /*
+  //time-interval based on the smallest one containing a minimum fraction of hits
+  // vector with time values of the hit; fraction between 0 and 1; how much furher enlarge the selected time window
+  float highestDensityFraction(std::vector<float>& hitTimes, float fractionToKeep=0.68, float timeWidthBy=0.5){
+
+    std::sort(hitTimes.begin(), hitTimes.end());
+    int totSize = hitTimes.size();
+    int num = 0.;
+    float sum = 0.;
+    float minTimeDiff = 999.;
+    int startTimeBin = 0;
+
+    int totToKeep = int(totSize*fractionToKeep);
+    int maxStart = totSize - totToKeep;
+
+    for(int ij=0; ij<maxStart; ++ij){
+      float localDiff = fabs(hitTimes[ij] - hitTimes[int(ij+totToKeep)]);
+      if(localDiff < minTimeDiff){
+	minTimeDiff = localDiff;
+	startTimeBin = ij;
+      }
+    }
+
+    // further adjust time width around the chosen one based on the hits density
+    // proved to improve the resolution: get as many hits as possible provided they are close in time
+
+    int startBin = startTimeBin;
+    int endBin = startBin+totToKeep;
+    float HalfTimeDiff = std::abs(hitTimes[startBin] - hitTimes[endBin]) * timeWidthBy;
+
+    for(int ij=0; ij<startBin; ++ij){
+      if(hitTimes[ij] > (hitTimes[startBin] - HalfTimeDiff) ){
+	for(int kl=ij; kl<totSize; ++kl){
+	  if(hitTimes[kl] < (hitTimes[endBin] + HalfTimeDiff) ){
+	    sum += hitTimes[kl];
+	    ++num;
+	  }
+	  else  break;
+	}
+	break;
+      }
+    }
+
+    if(num == 0) return -99.;
+    return sum/num;
+  }
+  */
+}
+
+#endif

RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.cc

@@ -10,7 +10,9 @@
 #include "FWCore/Framework/interface/Event.h"
 #include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
-#include "RealisticHitToClusterAssociator.h"
+#include "RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticHitToClusterAssociator.h"
+#include "RecoParticleFlow/PFClusterProducer/plugins/SimMappers/ComputeClusterTime.h"
+#include "RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticCluster.h"
 
 #ifdef PFLOW_DEBUG
 #define LOGVERB(x) edm::LogVerbatim(x)
@@ -109,13 +111,17 @@ void RealisticSimClusterMapper::buildClusters(const edm::Handle<reco::PFRecHitCo
 
     const auto& realisticClusters = realisticAssociator.realisticClusters();
     unsigned int nClusters = realisticClusters.size();
+    float bin_norm = 1. / (calibMaxEta_ - calibMinEta_);
+    float egamma_bin_norm = egammaCalib_.size() * bin_norm;
+    float hadron_bin_norm = hadronCalib_.size() * bin_norm;
     for (unsigned ic = 0; ic < nClusters; ++ic)
     {
         float highest_energy = 0.0f;
         output.emplace_back();
         reco::PFCluster& back = output.back();
         edm::Ref < std::vector<reco::PFRecHit> > seed;
         float energyCorrection = 1.f;
+        float timeRealisticSC = -99.;
         if (realisticClusters[ic].isVisible())
         {
             int pdgId = simClusters[ic].pdgId();
@@ -125,20 +131,18 @@ void RealisticSimClusterMapper::buildClusters(const edm::Handle<reco::PFRecHitCo
                 if ((isEGamma(pdgId) or isPi0(pdgId)) and !egammaCalib_.empty())
                 {
                     unsigned int etabin = std::floor(
-                            ((abseta - calibMinEta_) * egammaCalib_.size())
-                                    / (calibMaxEta_ - calibMinEta_));
-
+                            ((abseta - calibMinEta_) * egamma_bin_norm));
                     energyCorrection = egammaCalib_[etabin];
                 }
                 else if (isHadron(pdgId) and !(isPi0(pdgId)) and !hadronCalib_.empty()) // this function is expensive.. should we treat as hadron everything which is not egamma?
                 {
                     unsigned int etabin = std::floor(
-                            ((abseta - calibMinEta_) * hadronCalib_.size())
-                                    / (calibMaxEta_ - calibMinEta_));
+                            ((abseta - calibMinEta_) * hadron_bin_norm));
                     energyCorrection = hadronCalib_[etabin];
                 }
 
             }
+	    std::vector<float> timeHits;
             const auto& hitsIdsAndFractions = realisticClusters[ic].hitsIdsAndFractions();
             for (const auto& idAndF : hitsIdsAndFractions)
             {
@@ -153,14 +157,26 @@ void RealisticSimClusterMapper::buildClusters(const edm::Handle<reco::PFRecHitCo
                         highest_energy = hit_energy;
                         seed = ref;
                     }
+                    //select hits good for timing
+                    if(ref->time() > -1. ){
+		      int rhLayer = rhtools_.getLayerWithOffset(ref->detId());
+		      std::array<float,3> scPosition = realisticClusters[ic].getCenterOfGravity(rhLayer);
+		      float distanceSquared = std::pow((ref->position().x() - scPosition[0]),2) + std::pow((ref->position().y() - scPosition[1]),2);
+		      if(distanceSquared < maxDforTimingSquared_){
+			timeHits.push_back(ref->time() - timeOffset_);
+		      }
+		    }
                 }
             }
-        }
+	    //assign time if minimum number of hits
+	    if(timeHits.size() >= minNHitsforTiming_) timeRealisticSC = hgcalsimclustertime::fixSizeHighestDensity(timeHits);
+	}
         if (!back.hitsAndFractions().empty())
         {
             back.setSeed(seed->detId());
             back.setEnergy(realisticClusters[ic].getEnergy());
             back.setCorrectedEnergy(energyCorrection*realisticClusters[ic].getEnergy()); //applying energy correction
+	    back.setTime(timeRealisticSC);
         }
         else
         {

RecoParticleFlow/PFClusterProducer/plugins/SimMappers/RealisticSimClusterMapper.h

@@ -20,6 +20,9 @@ class RealisticSimClusterMapper : public InitialClusteringStepBase {
     exclusiveFraction_(conf.getParameter<double>("exclusiveFraction")),
     maxDistanceFilter_(conf.getParameter<bool>("maxDistanceFilter")),
     maxDistance_(conf.getParameter<double>("maxDistance")),
+    maxDforTimingSquared_(conf.getParameter<double>("maxDforTimingSquared")),
+    timeOffset_(conf.getParameter<double>("timeOffset")),
+    minNHitsforTiming_(conf.getParameter<unsigned int>("minNHitsforTiming")),
     useMCFractionsForExclEnergy_(conf.getParameter<bool>("useMCFractionsForExclEnergy")),
     calibMinEta_(conf.getParameter<double>("calibMinEta")),
     calibMaxEta_(conf.getParameter<double>("calibMaxEta"))
@@ -47,6 +50,9 @@ class RealisticSimClusterMapper : public InitialClusteringStepBase {
   const float exclusiveFraction_ = 0.7f;
   const bool maxDistanceFilter_ = false;
   const float maxDistance_ = 10.f;
+  const float maxDforTimingSquared_ = 4.0f;
+  const float timeOffset_;
+  const unsigned int minNHitsforTiming_ = 3;
   const bool useMCFractionsForExclEnergy_ = false;
   const float calibMinEta_ = 1.4;
   const float calibMaxEta_ = 3.0;

RecoParticleFlow/PFClusterProducer/python/particleFlowClusterHGC_cfi.py

@@ -1,5 +1,6 @@
 import FWCore.ParameterSet.Config as cms
 from RecoParticleFlow.PFClusterProducer.particleFlowRealisticSimClusterHGCCalibrations_cfi import *
+from SimCalorimetry.HGCalSimProducers.hgcalDigitizer_cfi import *
 #### PF CLUSTER HGCal ####
 
 #cleaning (none for now)
@@ -20,21 +21,25 @@
     maxDistanceFilter = cms.bool(True),
     #filtering out hits outside a cylinder of 10cm radius, built around the center of gravity per each layer
     maxDistance =  cms.double(10.0),
-    useMCFractionsForExclEnergy = cms.bool(False),    
+    maxDforTimingSquared = cms.double(4.0),
+    timeOffset = hgceeDigitizer.tofDelay,
+    minNHitsforTiming = cms.uint32(3),
+    useMCFractionsForExclEnergy = cms.bool(False),
     thresholdsByDetector = cms.VPSet(
     ),
     hadronCalib = hadronCorrections,
-    egammaCalib = egammaCorrections, 
-    calibMinEta = minEtaCorrection, 
-    calibMaxEta = maxEtaCorrection,                                   
+    egammaCalib = egammaCorrections,
+    calibMinEta = minEtaCorrection,
+    calibMaxEta = maxEtaCorrection,
     simClusterSrc = cms.InputTag("mix:MergedCaloTruth")
 )
 
 
 #position calculations
 _positionCalcPCA_HGCal = cms.PSet(
         algoName = cms.string("Cluster3DPCACalculator"),
-        minFractionInCalc = cms.double(1e-9)
+        minFractionInCalc = cms.double(1e-9),
+        updateTiming = cms.bool(False)
 )
 
 _hgcalMultiClusterMapper_HGCal = cms.PSet(