PFSimProducer/plugins/EcalBarrelClusterFastTimer.cc

0001 // This producer eats standard PF ECAL clusters
0002 // finds the corresponding fast-timing det IDs and attempts to
0003 // assign a reasonable time guess
0004
0005 #include "CLHEP/Units/SystemOfUnits.h"
0006 #include "DataFormats/Common/interface/ValueMap.h"
0007 #include "DataFormats/Common/interface/View.h"
0008 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0009 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0010 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0011 #include "DataFormats/VertexReco/interface/Vertex.h"
0012 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/global/EDProducer.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0020 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0023 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0024 #include "SimTracker/TrackAssociation/interface/ResolutionModel.h"
0025
0026 #include <random>
0027 #include <memory>
0028
0029 class EcalBarrelClusterFastTimer : public edm::global::EDProducer<> {
0030 public:
0031   EcalBarrelClusterFastTimer(const edm::ParameterSet&);
0032
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0036
0037 private:
0038   // inputs
0039   edm::EDGetTokenT<EcalRecHitCollection> ebTimeHitsToken_;
0040   edm::EDGetTokenT<std::vector<reco::PFCluster>> ebClustersToken_;
0041   // options
0042   std::vector<std::unique_ptr<const ResolutionModel>> _resolutions;
0043   const float minFraction_, minEnergy_;
0044   const unsigned ecalDepth_;
0045   // functions
0046   std::pair<float, DetId> getTimeForECALPFCluster(const EcalRecHitCollection&, const reco::PFCluster&) const;
0047   float correctTimeToVertex(const float intime,
0048                             const DetId& timeDet,
0049                             const reco::Vertex& vtx,
0050                             const CaloSubdetectorGeometry* ecalGeom) const;
0051 };
0052
0053 #include "FWCore/Framework/interface/MakerMacros.h"
0054 DEFINE_FWK_MODULE(EcalBarrelClusterFastTimer);
0055
0056 void EcalBarrelClusterFastTimer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0057   // ecalBarrelClusterFastTimer
0058   edm::ParameterSetDescription desc;
0059   desc.add<double>("ecalDepth", 7.0);
0060   {
0061     edm::ParameterSetDescription vpsd1;
0062     vpsd1.add<std::string>("modelName", "PerfectResolutionModel");
0063     std::vector<edm::ParameterSet> temp1;
0064     temp1.reserve(1);
0065     {
0066       edm::ParameterSet temp2;
0067       temp2.addParameter<std::string>("modelName", "PerfectResolutionModel");
0068       temp1.push_back(temp2);
0069     }
0070     desc.addVPSet("resolutionModels", vpsd1, temp1);
0071   }
0072   desc.add<edm::InputTag>("timedVertices", edm::InputTag("offlinePrimaryVertices4D"));
0073   desc.add<double>("minFractionToConsider", 0.1);
0074   desc.add<edm::InputTag>("ebClusters", edm::InputTag("particleFlowClusterECALUncorrected"));
0075   desc.add<double>("minEnergyToConsider", 0.0);
0076   desc.add<edm::InputTag>("ebTimeHits", edm::InputTag("ecalDetailedTimeRecHit", "EcalRecHitsEB"));
0077   descriptions.add("ecalBarrelClusterFastTimer", desc);
0078 }
0079
0080 namespace {
0081   const std::string resolution("Resolution");
0082
0083   constexpr float cm_per_ns = 29.9792458;
0084
0085   template <typename T>
0086   void writeValueMap(edm::Event& iEvent,
0087                      const edm::Handle<std::vector<reco::PFCluster>>& handle,
0088                      const std::vector<T>& values,
0089                      const std::string& label) {
0090     auto valMap = std::make_unique<edm::ValueMap<T>>();
0091     typename edm::ValueMap<T>::Filler filler(*valMap);
0092     filler.insert(handle, values.begin(), values.end());
0093     filler.fill();
0094     iEvent.put(std::move(valMap), label);
0095   }
0096 }  // namespace
0097
0098 EcalBarrelClusterFastTimer::EcalBarrelClusterFastTimer(const edm::ParameterSet& conf)
0099     : ebTimeHitsToken_(consumes<EcalRecHitCollection>(conf.getParameter<edm::InputTag>("ebTimeHits"))),
0100       ebClustersToken_(consumes<std::vector<reco::PFCluster>>(conf.getParameter<edm::InputTag>("ebClusters"))),
0101       minFraction_(conf.getParameter<double>("minFractionToConsider")),
0102       minEnergy_(conf.getParameter<double>("minEnergyToConsider")),
0103       ecalDepth_(conf.getParameter<double>("ecalDepth")) {
0104   // setup resolution models
0105   const std::vector<edm::ParameterSet>& resos = conf.getParameterSetVector("resolutionModels");
0106   for (const auto& reso : resos) {
0107     const std::string& name = reso.getParameter<std::string>("modelName");
0108     _resolutions.emplace_back(ResolutionModelFactory::get()->create(name, reso));
0109
0110     // times and time resolutions for general tracks
0111     produces<edm::ValueMap<float>>(name);
0112     produces<edm::ValueMap<float>>(name + resolution);
0113   }
0114 }
0115
0116 void EcalBarrelClusterFastTimer::produce(edm::StreamID sid, edm::Event& evt, const edm::EventSetup& es) const {
0117   edm::Handle<std::vector<reco::PFCluster>> clustersH;
0118   edm::Handle<EcalRecHitCollection> timehitsH;
0119
0120   evt.getByToken(ebClustersToken_, clustersH);
0121   evt.getByToken(ebTimeHitsToken_, timehitsH);
0122
0123   const auto& clusters = *clustersH;
0124   const auto& timehits = *timehitsH;
0125
0126   // get event-based seed for RNG
0127   unsigned int runNum_uint = static_cast<unsigned int>(evt.id().run());
0128   unsigned int lumiNum_uint = static_cast<unsigned int>(evt.id().luminosityBlock());
0129   unsigned int evNum_uint = static_cast<unsigned int>(evt.id().event());
0130   std::uint32_t seed = (lumiNum_uint << 10) + (runNum_uint << 20) + evNum_uint;
0131   std::mt19937 rng(seed);
0132
0133   std::vector<std::pair<float, DetId>> times;  // perfect times keyed to cluster index
0134   times.reserve(clusters.size());
0135
0136   for (const auto& cluster : clusters) {
0137     times.push_back(getTimeForECALPFCluster(timehits, cluster));
0138   }
0139
0140   for (const auto& reso : _resolutions) {
0141     const std::string& name = reso->name();
0142     std::vector<float> resolutions;
0143     std::vector<float> smeared_times;
0144     resolutions.reserve(clusters.size());
0145     smeared_times.reserve(clusters.size());
0146
0147     // smear once then correct to multiple vertices
0148     for (unsigned i = 0; i < clusters.size(); ++i) {
0149       const float theresolution = reso->getTimeResolution(clusters[i]);
0150       std::normal_distribution<float> gausTime(times[i].first, theresolution);
0151
0152       smeared_times.emplace_back(gausTime(rng));
0153       resolutions.push_back(theresolution);
0154     }
0155
0156     writeValueMap(evt, clustersH, smeared_times, name);
0157     writeValueMap(evt, clustersH, resolutions, name + resolution);
0158   }
0159 }
0160
0161 std::pair<float, DetId> EcalBarrelClusterFastTimer::getTimeForECALPFCluster(const EcalRecHitCollection& timehits,
0162                                                                             const reco::PFCluster& cluster) const {
0163   const auto& rhfs = cluster.recHitFractions();
0164   unsigned best_hit = 0;
0165   float best_energy = -1.f;
0166   for (const auto& rhf : rhfs) {
0167     const auto& hitref = rhf.recHitRef();
0168     const unsigned detid = hitref->detId();
0169     const auto fraction = rhf.fraction();
0170     const auto energy = hitref->energy();
0171     if (fraction < minFraction_ || energy < minEnergy_)
0172       continue;
0173     auto timehit = timehits.find(detid);
0174     float e_hit = rhf.fraction() * hitref->energy();
0175     if (e_hit > best_energy && timehit->isTimeValid()) {
0176       best_energy = e_hit;
0177       best_hit = detid;
0178     }
0179   }
0180
0181   float best_time_guess = std::numeric_limits<float>::max();
0182   if (best_energy > 0.) {
0183     best_time_guess = timehits.find(best_hit)->time();
0184   }
0185
0186   //std::cout << "EcalBarrelFastTimer: " << best_time_guess << ' ' << best_energy << ' ' << best_hit << std::endl;
0187
0188   return std::make_pair(best_time_guess, DetId(best_hit));
0189 }
0190
0191 float EcalBarrelClusterFastTimer::correctTimeToVertex(const float intime,
0192                                                       const DetId& timeDet,
0193                                                       const reco::Vertex& vtx,
0194                                                       const CaloSubdetectorGeometry* ecalGeom) const {
0195   if (timeDet.rawId() == 0)
0196     return -999.;
0197   // correct the cluster time from 0,0,0 to the primary vertex given
0198   auto cellGeometry = ecalGeom->getGeometry(timeDet);
0199   if (nullptr == cellGeometry) {
0200     throw cms::Exception("BadECALBarrelCell")
0201         << std::hex << timeDet.rawId() << std::dec << " is not a valid ECAL Barrel DetId!";
0202   }
0203   //depth in mm in the middle of the layer position;
0204   GlobalPoint layerPos = cellGeometry->getPosition(ecalDepth_ + 0.5);
0205   const math::XYZPoint layerPos_cm(layerPos.x(), layerPos.y(), layerPos.z());
0206   const math::XYZVector to_center = layerPos_cm - math::XYZPoint(0., 0., 0.);
0207   const math::XYZVector to_vtx = layerPos_cm - vtx.position();
0208
0209   /*
0210   std::cout << intime << ' ' << to_center.r()/cm_per_ns << ' ' << to_vtx.r()/cm_per_ns
0211             << ' ' << intime + to_center.r()/cm_per_ns - to_vtx.r()/cm_per_ns << std::endl;
0212   */
0213
0214   return intime + to_center.r() / cm_per_ns - to_vtx.r() / cm_per_ns;
0215 }