File indexing completed on 2021-04-27 00:16:59
0001
0002
0003
0004
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
0022
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
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
0043 private:
0044 void produce(edm::Event&, const edm::EventSetup&) override;
0045 void beginLuminosityBlock(const edm::LuminosityBlock&, const edm::EventSetup&) override;
0046
0047 edm::EDGetTokenT<edm::View<reco::PFRecTrack> > src_;
0048
0049
0050
0051
0052 const reco::TrackBase::TrackQuality trackQuality_;
0053 const std::vector<double> DPtovPtCut_;
0054 const std::vector<unsigned> NHitCut_;
0055 const bool useIterTracking_;
0056
0057
0058
0059 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0060 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> tkerGeomToken_;
0061 const MagneticField* bField_;
0062 const TrackerGeometry* tkGeom_;
0063 std::array<std::string, 1> hgc_names_;
0064 std::array<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>, 1>
0065 hgcGeometryTokens_;
0066 std::array<const HGCalGeometry*, 1> hgcGeometries_;
0067 std::array<std::vector<ReferenceCountingPointer<BoundDisk> >, 1> plusSurface_,
0068 minusSurface_;
0069 std::unique_ptr<PropagatorWithMaterial> mat_prop_;
0070
0071 float diskOuterRadius_;
0072 float diskInnerRadius_;
0073 };
0074
0075 HGCalTrackCollectionProducer::HGCalTrackCollectionProducer(const edm::ParameterSet& iConfig)
0076 : src_(consumes<edm::View<reco::PFRecTrack> >(iConfig.getParameter<edm::InputTag>("src"))),
0077 trackQuality_((iConfig.existsAs<std::string>("trackQuality"))
0078 ? reco::TrackBase::qualityByName(iConfig.getParameter<std::string>("trackQuality"))
0079 : reco::TrackBase::highPurity),
0080 DPtovPtCut_(iConfig.getParameter<std::vector<double> >("DPtOverPtCuts_byTrackAlgo")),
0081 NHitCut_(iConfig.getParameter<std::vector<unsigned> >("NHitCuts_byTrackAlgo")),
0082 useIterTracking_(iConfig.getParameter<bool>("useIterativeTracking")),
0083 magneticFieldToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()),
0084 tkerGeomToken_(esConsumes<edm::Transition::BeginLuminosityBlock>()) {
0085 LogDebug("HGCalTrackCollectionProducer")
0086 << " HGCalTrackCollectionProducer::HGCalTrackCollectionProducer " << std::endl;
0087
0088 const edm::ParameterSet& geoconf = iConfig.getParameterSet("hgcalGeometryNames");
0089 hgc_names_[0] = geoconf.getParameter<std::string>("HGC_ECAL");
0090 for (unsigned i = 0; i < hgcGeometryTokens_.size(); i++) {
0091 hgcGeometryTokens_[i] = esConsumes<edm::Transition::BeginLuminosityBlock>(edm::ESInputTag("", hgc_names_[i]));
0092 }
0093
0094 produces<reco::PFRecTrackCollection>("TracksInHGCal");
0095 produces<reco::PFRecTrackCollection>("TracksNotInHGCal");
0096 }
0097
0098
0099
0100 void HGCalTrackCollectionProducer::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& es) {
0101 constexpr float m_pion = 0.1396;
0102
0103 bField_ = &es.getData(magneticFieldToken_);
0104 tkGeom_ = &es.getData(tkerGeomToken_);
0105
0106 for (unsigned i = 0; i < hgcGeometries_.size(); ++i) {
0107 hgcGeometries_[i] = &es.getData(hgcGeometryTokens_[i]);
0108 }
0109
0110
0111 mat_prop_ = std::make_unique<PropagatorWithMaterial>(alongMomentum, m_pion, bField_);
0112
0113 Surface::RotationType rot;
0114 for (unsigned i = 0; i < hgcGeometries_.size(); ++i) {
0115 minusSurface_[i].clear();
0116 plusSurface_[i].clear();
0117 const HGCalDDDConstants& dddCons = hgcGeometries_[i]->topology().dddConstants();
0118 std::map<float, float> zrhoCoord;
0119 std::map<float, float> innerRadiusCoord;
0120 auto theTrForms = dddCons.getTrForms();
0121 const auto& firstLayerIt = theTrForms.back();
0122 float Z(std::abs(firstLayerIt.h3v.z()));
0123
0124 diskInnerRadius_ = 31.5;
0125 diskOuterRadius_ = 161.0f;
0126 LogDebug("HGCalTrackCollectionProducer") << "O HAI I'm making a bound disk with Outer R=" << diskOuterRadius_
0127 << " Inner R=" << diskInnerRadius_ << " and Z=" << Z << std::endl;
0128 minusSurface_[i].push_back(ReferenceCountingPointer<BoundDisk>(new BoundDisk(
0129 Surface::PositionType(0, 0, -Z), rot, new SimpleDiskBounds(diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
0130 plusSurface_[i].push_back(ReferenceCountingPointer<BoundDisk>(new BoundDisk(
0131 Surface::PositionType(0, 0, +Z), rot, new SimpleDiskBounds(diskInnerRadius_, diskOuterRadius_, -0.001, 0.001))));
0132 }
0133 }
0134
0135 void HGCalTrackCollectionProducer::produce(edm::Event& evt, const edm::EventSetup& iSetup) {
0136 edm::Handle<edm::View<reco::PFRecTrack> > trackHandle;
0137 evt.getByToken(src_, trackHandle);
0138 const auto& tracks = *trackHandle;
0139
0140 auto outputInHGCal = std::make_unique<reco::PFRecTrackCollection>();
0141 auto outputNotInHGCal = std::make_unique<reco::PFRecTrackCollection>();
0142
0143 for (unsigned int i = 0; i < tracks.size(); i++) {
0144 const auto track = tracks.ptrAt(i);
0145 bool isGood =
0146 PFTrackAlgoTools::goodPtResolution(track->trackRef(), DPtovPtCut_, NHitCut_, useIterTracking_, trackQuality_);
0147 LogDebug("HGCalTrackCollectionProducer") << "HGCalTrackCollectionProducer Track number " << i
0148 << " has a goodPtResolution result of " << isGood << std::endl;
0149 if (!isGood)
0150 continue;
0151 bool found = false;
0152 const TrajectoryStateOnSurface myTSOS =
0153 trajectoryStateTransform::outerStateOnSurface(*(track->trackRef()), *tkGeom_, bField_);
0154 auto detbegin = myTSOS.globalPosition().z() > 0 ? plusSurface_.begin() : minusSurface_.begin();
0155 auto detend = myTSOS.globalPosition().z() > 0 ? plusSurface_.end() : minusSurface_.end();
0156 for (auto det = detbegin; det != detend; ++det) {
0157 LogDebug("HGCalTrackCollectionProducer") << "at HGC detector: " << std::distance(detbegin, det) << std::endl;
0158 unsigned layer_count = 1;
0159 for (const auto& layer : *det) {
0160 LogDebug("HGCalTrackCollectionProducer") << " at DET layer: " << layer_count++ << std::endl;
0161 TrajectoryStateOnSurface piStateAtSurface = mat_prop_->propagate(myTSOS, *layer);
0162 if (piStateAtSurface.isValid()) {
0163 LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is valid!" << std::endl;
0164 GlobalPoint pt = piStateAtSurface.globalPosition();
0165 if (pt.perp() < diskOuterRadius_) {
0166 if (pt.perp() > diskInnerRadius_) {
0167 LogDebug("HGCalTrackCollectionProducer")
0168 << "(x,y,z,r)=(" << pt.x() << ", " << pt.y() << ", " << pt.z() << ", "
0169 << sqrt(pt.x() * pt.x() + pt.y() * pt.y()) << ")" << std::endl;
0170 if (std::abs(track->trackRef()->eta()) < 1.47)
0171 LogDebug("HGCalTrackCollectionProducer") << " ETA IN BARREL REGION: " << track->trackRef()->eta()
0172 << " (PT: " << track->trackRef()->pt() << ")" << std::endl;
0173 found = true;
0174 } else {
0175 LogDebug("HGCalTrackCollectionProducer")
0176 << " but r=" << pt.perp() << " < diskInnerRadius=" << diskInnerRadius_ << " so skipping "
0177 << std::endl;
0178 }
0179 } else {
0180 LogDebug("HGCalTrackCollectionProducer")
0181 << " but r=" << pt.perp() << " > diskOuterRadius=" << diskOuterRadius_ << " so skipping " << std::endl;
0182 }
0183 } else {
0184 LogDebug("HGCalTrackCollectionProducer") << "Extrapolation is NOT valid!" << std::endl;
0185 }
0186 }
0187 }
0188 if (found) {
0189 LogDebug("HGCalTrackCollectionProducer") << " Track going to outputInHGCal pt eta " << track->trackRef()->pt()
0190 << " " << track->trackRef()->eta() << std::endl;
0191 outputInHGCal->push_back(*track);
0192 } else {
0193 outputNotInHGCal->push_back(*track);
0194 }
0195 }
0196
0197 evt.put(std::move(outputInHGCal), "TracksInHGCal");
0198 evt.put(std::move(outputNotInHGCal), "TracksNotInHGCal");
0199 }
0200
0201 DEFINE_FWK_MODULE(HGCalTrackCollectionProducer);