PFTracking/plugins/HGCalTrackCollectionProducer.cc

0001 // S. Zenz, 12 February 2015
0002 //
0003 // Splits a track collection into two, based on whether they propagate to the HGCal or not
0004 // Tracks with bad pt resolution (suspected fakes) are dropped and not in either collection
0005
0006 #include "FWCore/Framework/interface/MakerMacros.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h"
0011 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0014 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0015 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0016 #include "DataFormats/Common/interface/Handle.h"
0017 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0018 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020
0021 // for track propagation through HGC
0022 // N.B. we are only propogating to first layer, so check these later
0023 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0026 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "MagneticField/Engine/interface/MagneticField.h"
0029 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0030
0031 //geometry records
0032 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0033 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0034 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0035
0036 #include <memory>
0037 #include <unordered_map>
0038
0039 class HGCalTrackCollectionProducer : public edm::stream::EDProducer<> {
0040 public:
0041   HGCalTrackCollectionProducer(const edm::ParameterSet&);
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045   void produce(edm::Event&, const edm::EventSetup&) override;
0046   void beginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup&) override;
0047
0048   edm::EDGetTokenT<edm::View<reco::PFRecTrack>> src_;
0049
0050   // variables needed for copied goodPtResolution function
0051   // need to go back and figure out sensible values
0052
0053   const reco::TrackBase::TrackQuality trackQuality_;
0054   const std::vector<double> DPtovPtCut_;
0055   const std::vector<unsigned> NHitCut_;
0056   const bool useIterTracking_;
0057   //  const bool _useFirstLayerOnly; // always true now
0058
0059   // variables needed for copied extrapolation
0060   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0061   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkerGeomToken_;
0062   const MagneticField* bField_;
0063   const TrackerGeometry* tkGeom_;
0064   std::array<std::string, 1> hgc_names_;  // 3 --> 1; extrapolate to hgcee only
0065   std::array<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>, 1>
0066       hgcGeometryTokens_;                              // 3 --> 1; extrapolate to hgcee only
0067   std::array<const HGCalGeometry*, 1> hgcGeometries_;  // 3 --> 1; extrapolate to hgcee only
0068   std::array<std::vector<ReferenceCountingPointer<BoundDisk>>, 1> plusSurface_,
0069       minusSurface_;  // 3 --> 1; extrapolate to hgcee only
0070   std::unique_ptr<PropagatorWithMaterial> mat_prop_;
0071
0072   float diskOuterRadius_;
0073   float diskInnerRadius_;
0074 };
0075
0076 void HGCalTrackCollectionProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078   desc.add<edm::InputTag>("src", {"pfTrack"});
0079   desc.add<std::string>("trackQuality", "highPurity");
0080   desc.add<bool>("useIterativeTracking", true);
0081   desc.add<std::vector<double>>("DPtOverPtCuts_byTrackAlgo", {10.0, 10.0, 10.0, 10.0, 10.0, 5.0});
0082   desc.add<std::vector<uint32_t>>("NHitCuts_byTrackAlgo", {3, 3, 3, 3, 3, 3});
0083   edm::ParameterSetDescription pset;
0084   pset.add<std::string>("HGC_ECAL", "HGCalEESensitive");
0085   desc.add<edm::ParameterSetDescription>("hgcalGeometryNames", pset);
0086   descriptions.add("hgcalTrackCollection", desc);
0087 }
0088
0089 HGCalTrackCollectionProducer::HGCalTrackCollectionProducer(const edm::ParameterSet& iConfig)
0090     : src_(consumes<edm::View<reco::PFRecTrack>>(iConfig.getParameter<edm::InputTag>("src"))),
0091       trackQuality_(reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"))),
0092       DPtovPtCut_(iConfig.getParameter<std::vector<double>>("DPtOverPtCuts_byTrackAlgo")),
0093       NHitCut_(iConfig.getParameter<std::vector<unsigned>>("NHitCuts_byTrackAlgo")),
0094       useIterTracking_(iConfig.getParameter<bool>("useIterativeTracking")),
0095       magneticFieldToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()),
0096       tkerGeomToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()) {
0097   LogDebug("HGCalTrackCollectionProducer")
0098       << " HGCalTrackCollectionProducer::HGCalTrackCollectionProducer " << std::endl;
0099
0100   const edm::ParameterSet& geoconf = iConfig.getParameterSet("hgcalGeometryNames");
0101   hgc_names_[0] = geoconf.getParameter<std::string>("HGC_ECAL");
0102   for (unsigned i = 0; i < hgcGeometryTokens_.size(); i++) {
0103     hgcGeometryTokens_[i] = esConsumes<edm::Transition::BeginLuminosityBlock>(edm::ESInputTag("", hgc_names_[i]));
0104   }
0105
0106   produces<reco::PFRecTrackCollection>("TracksInHGCal");
0107   produces<reco::PFRecTrackCollection>("TracksNotInHGCal");
0108 }
0109
0110 // From https://github.com/cms-sw/cmssw/blob/CMSSW_6_2_X_SLHC/RecoParticleFlow/PFClusterProducer/src/HGCClusterizer.cc#L441-L447 and beyond
0111 // TODO: we only need the front of the calorimeter, so modify this
0112 void HGCalTrackCollectionProducer::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& es) {
0113   constexpr float m_pion = 0.1396;
0114   // get dependencies for setting up propagator
0115   bField_ = &es.getData(magneticFieldToken_);
0116   tkGeom_ = &es.getData(tkerGeomToken_);
0117   // get HGC geometries (assume that layers are ordered in Z!)
0118   for (unsigned i = 0; i < hgcGeometries_.size(); ++i) {
0119     hgcGeometries_[i] = &es.getData(hgcGeometryTokens_[i]);
0120   }
0121
0122   // make propagator
0123   mat_prop_ = std::make_unique<PropagatorWithMaterial>(alongMomentum, m_pion, bField_);
0124   // setup HGC layers for track propagation
0125   Surface::RotationType rot;  //unit rotation matrix
0126   for (unsigned i = 0; i < hgcGeometries_.size(); ++i) {
0127     minusSurface_[i].clear();
0128     plusSurface_[i].clear();
0129     const HGCalDDDConstants& dddCons = hgcGeometries_[i]->topology().dddConstants();
0130     std::map<float, float> zrhoCoord;
0131     std::map<float, float> innerRadiusCoord;
0132     auto theTrForms = dddCons.getTrForms();
0133     const auto& firstLayerIt = theTrForms.back();
0134     float Z(std::abs(firstLayerIt.h3v.z()));
0135     // use hardcoded radii for now (FIX ME)
0136     diskInnerRadius_ = 31.5;
0137     diskOuterRadius_ = 161.0f;
0138     LogDebug("HGCalTrackCollectionProducer") << "O HAI I'm making a bound disk with Outer R=" << diskOuterRadius_
0139                                              << " Inner R=" << diskInnerRadius_ << " and Z=" << Z << std::endl;
0140     minusSurface_[i].push_back(ReferenceCountingPointer<BoundDisk>(new BoundDisk(
0141         Surface::PositionType(0, 0, -Z), rot, new SimpleDiskBounds(diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
0142     plusSurface_[i].push_back(ReferenceCountingPointer<BoundDisk>(new BoundDisk(
0143         Surface::PositionType(0, 0, +Z), rot, new SimpleDiskBounds(diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
0144   }
0145 }
0146
0147 void HGCalTrackCollectionProducer::produce(edm::Event& evt, const edm::EventSetup& iSetup) {
0148   edm::Handle<edm::View<reco::PFRecTrack>> trackHandle;
0149   evt.getByToken(src_, trackHandle);
0150   const auto& tracks = *trackHandle;
0151
0152   auto outputInHGCal = std::make_unique<reco::PFRecTrackCollection>();
0153   auto outputNotInHGCal = std::make_unique<reco::PFRecTrackCollection>();
0154
0155   for (unsigned int i = 0; i < tracks.size(); i++) {
0156     const auto track = tracks.ptrAt(i);
0157     bool isGood =
0158         PFTrackAlgoTools::goodPtResolution(track->trackRef(), DPtovPtCut_, NHitCut_, useIterTracking_, trackQuality_);
0159     LogDebug("HGCalTrackCollectionProducer") << "HGCalTrackCollectionProducer Track number " << i
0160                                              << " has a goodPtResolution result of " << isGood << std::endl;
0161     if (!isGood)
0162       continue;
0163     bool found = false;
0164     const TrajectoryStateOnSurface myTSOS =
0165         trajectoryStateTransform::outerStateOnSurface(*(track->trackRef()), *tkGeom_, bField_);
0166     auto detbegin = myTSOS.globalPosition().z() > 0 ? plusSurface_.begin() : minusSurface_.begin();
0167     auto detend = myTSOS.globalPosition().z() > 0 ? plusSurface_.end() : minusSurface_.end();
0168     for (auto det = detbegin; det != detend; ++det) {
0169       LogDebug("HGCalTrackCollectionProducer") << "at HGC detector: " << std::distance(detbegin, det) << std::endl;
0170       unsigned layer_count = 1;
0171       for (const auto& layer : *det) {
0172         LogDebug("HGCalTrackCollectionProducer") << "  at DET layer: " << layer_count++ << std::endl;
0173         TrajectoryStateOnSurface piStateAtSurface = mat_prop_->propagate(myTSOS, *layer);
0174         if (piStateAtSurface.isValid()) {
0175           LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is valid!" << std::endl;
0176           GlobalPoint pt = piStateAtSurface.globalPosition();
0177           if (pt.perp() < diskOuterRadius_) {
0178             if (pt.perp() > diskInnerRadius_) {
0179               LogDebug("HGCalTrackCollectionProducer")
0180                   << "(x,y,z,r)=(" << pt.x() << ", " << pt.y() << ", " << pt.z() << ", "
0181                   << sqrt(pt.x() * pt.x() + pt.y() * pt.y()) << ")" << std::endl;
0182               if (std::abs(track->trackRef()->eta()) < 1.47)
0183                 LogDebug("HGCalTrackCollectionProducer") << " ETA IN BARREL REGION: " << track->trackRef()->eta()
0184                                                          << " (PT: " << track->trackRef()->pt() << ")" << std::endl;
0185               found = true;
0186             } else {
0187               LogDebug("HGCalTrackCollectionProducer")
0188                   << " but r=" << pt.perp() << " < diskInnerRadius=" << diskInnerRadius_ << " so skipping "
0189                   << std::endl;
0190             }
0191           } else {
0192             LogDebug("HGCalTrackCollectionProducer")
0193                 << " but r=" << pt.perp() << " > diskOuterRadius=" << diskOuterRadius_ << " so skipping " << std::endl;
0194           }
0195         } else {
0196           LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is NOT valid!" << std::endl;
0197         }
0198       }
0199     }
0200     if (found) {
0201       LogDebug("HGCalTrackCollectionProducer") << " Track going to outputInHGCal pt eta " << track->trackRef()->pt()
0202                                                << " " << track->trackRef()->eta() << std::endl;
0203       outputInHGCal->push_back(*track);
0204     } else {
0205       outputNotInHGCal->push_back(*track);
0206     }
0207   }  // track loop
0208
0209   evt.put(std::move(outputInHGCal), "TracksInHGCal");
0210   evt.put(std::move(outputNotInHGCal), "TracksNotInHGCal");
0211 }
0212
0213 DEFINE_FWK_MODULE(HGCalTrackCollectionProducer);