Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-04-27 00:16:59 

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 
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   // variables needed for copied goodPtResolution function
0050   // need to go back and figure out sensible values
0051 
0052   const reco::TrackBase::TrackQuality trackQuality_;
0053   const std::vector<double> DPtovPtCut_;
0054   const std::vector<unsigned> NHitCut_;
0055   const bool useIterTracking_;
0056   //  const bool _useFirstLayerOnly; // always true now
0057 
0058   // variables needed for copied extrapolation
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_;  // 3 --> 1; extrapolate to hgcee only
0064   std::array<edm::ESGetToken<HGCalGeometry, IdealGeometryRecord>, 1>
0065       hgcGeometryTokens_;                              // 3 --> 1; extrapolate to hgcee only
0066   std::array<const HGCalGeometry*, 1> hgcGeometries_;  // 3 --> 1; extrapolate to hgcee only
0067   std::array<std::vector<ReferenceCountingPointer<BoundDisk> >, 1> plusSurface_,
0068       minusSurface_;  // 3 --> 1; extrapolate to hgcee only
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 // From https://github.com/cms-sw/cmssw/blob/CMSSW_6_2_X_SLHC/RecoParticleFlow/PFClusterProducer/src/HGCClusterizer.cc#L441-L447 and beyond
0099 // TODO: we only need the front of the calorimeter, so modify this
0100 void HGCalTrackCollectionProducer::beginLuminosityBlock(const edm::LuminosityBlock& lumi, const edm::EventSetup& es) {
0101   constexpr float m_pion = 0.1396;
0102   // get dependencies for setting up propagator
0103   bField_ = &es.getData(magneticFieldToken_);
0104   tkGeom_ = &es.getData(tkerGeomToken_);
0105   // get HGC geometries (assume that layers are ordered in Z!)
0106   for (unsigned i = 0; i < hgcGeometries_.size(); ++i) {
0107     hgcGeometries_[i] = &es.getData(hgcGeometryTokens_[i]);
0108   }
0109 
0110   // make propagator
0111   mat_prop_ = std::make_unique<PropagatorWithMaterial>(alongMomentum, m_pion, bField_);
0112   // setup HGC layers for track propagation
0113   Surface::RotationType rot;  //unit rotation matrix
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     // use hardcoded radii for now (FIX ME)
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   }  // track loop
0196 
0197   evt.put(std::move(outputInHGCal), "TracksInHGCal");
0198   evt.put(std::move(outputNotInHGCal), "TracksNotInHGCal");
0199 }
0200 
0201 DEFINE_FWK_MODULE(HGCalTrackCollectionProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!