PixelCluster/plugins/PixelClusterTagInfoProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoBTag/PixelCluster
0004 // Class:      PixelClusterTagInfoProducer
0005 //
0006 /**\class PixelClusterTagInfoProducer PixelCluster RecoBTag/PixelCluster/plugins/PixelClusterTagInfoProducer.cc
0007
0008  Description: Produces a collection of PixelClusterTagInfo objects,
0009   that contain the pixel cluster hit multiplicity in each pixel layer or disk
0010   in a narrow cone around the jet axis.
0011
0012  Implementation:
0013      If the event does not fulfill minimum conditions (at least one jet above threshold,
0014      and a valid primary vertex) and empty collection is filled. Otherwise, a loop over
0015      the pixel cluster collection fills a vector of reco::PixelClusterProperties that
0016      contains the geometrical position and the charge of the cluster above threshold.
0017      A second loop on jets performs the dR association, and fills the TagInfo collection.
0018 */
0019 //
0020 // Original Author:  Alberto Zucchetta, Manuel Sommerhalder (UniZ) [zucchett]
0021 //         Created:  Wed, 03 Jul 2019 12:37:30 GMT
0022 //
0023 //
0024
0025 // system include files
0026 #include <memory>
0027
0028 // user include files
0029 #include "FWCore/Framework/interface/Frameworkfwd.h"
0030 #include "FWCore/Framework/interface/global/EDProducer.h"
0031
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/Framework/interface/MakerMacros.h"
0034
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/Utilities/interface/StreamID.h"
0037
0038 #include "FWCore/Framework/interface/ESHandle.h"
0039 #include "FWCore/Framework/interface/EventSetup.h"
0040 #include "FWCore/Utilities/interface/ESGetToken.h"
0041
0042 // TagInfo
0043 #include "DataFormats/BTauReco/interface/PixelClusterTagInfo.h"
0044
0045 // For vertices
0046 #include "DataFormats/VertexReco/interface/Vertex.h"
0047 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0048
0049 // For jet
0050 #include "DataFormats/JetReco/interface/PFJet.h"
0051 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0052
0053 // For pixel clusters and topology
0054 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0055 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0056 #include "DataFormats/Common/interface/Ref.h"
0057 #include "DataFormats/DetId/interface/DetId.h"
0058 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0059 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0060
0061 // Geometry
0062 #include "DataFormats/Math/interface/deltaR.h"
0063 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0064 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0065 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0066 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0067
0068 class PixelClusterTagInfoProducer : public edm::global::EDProducer<> {
0069 public:
0070   explicit PixelClusterTagInfoProducer(const edm::ParameterSet&);
0071   ~PixelClusterTagInfoProducer() override;
0072
0073   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0074
0075 private:
0076   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0077
0078   const edm::ParameterSet iConfig;
0079
0080   const edm::EDGetTokenT<edm::View<reco::Jet> > m_jets;
0081   const edm::EDGetTokenT<reco::VertexCollection> m_vertices;
0082   const edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > m_pixelhit;
0083   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> m_geomToken;
0084   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> m_topoToken;
0085   const bool m_isPhase1;
0086   const bool m_addFPIX;
0087   const int m_minADC;
0088   const float m_minJetPt;
0089   const float m_maxJetEta;
0090   const float m_hadronMass;
0091   const unsigned int m_nLayers;
0092 };
0093
0094 PixelClusterTagInfoProducer::PixelClusterTagInfoProducer(const edm::ParameterSet& iConfig)
0095     : m_jets(consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jets"))),
0096       m_vertices(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0097       m_pixelhit(consumes<edmNew::DetSetVector<SiPixelCluster> >(iConfig.getParameter<edm::InputTag>("pixelhit"))),
0098       m_geomToken(esConsumes()),
0099       m_topoToken(esConsumes()),
0100       m_isPhase1(iConfig.getParameter<bool>("isPhase1")),
0101       m_addFPIX(iConfig.getParameter<bool>("addForward")),
0102       m_minADC(iConfig.getParameter<int>("minAdcCount")),
0103       m_minJetPt(iConfig.getParameter<double>("minJetPtCut")),
0104       m_maxJetEta(iConfig.getParameter<double>("maxJetEtaCut")),
0105       m_hadronMass(iConfig.getParameter<double>("hadronMass")),
0106       m_nLayers(m_isPhase1 ? 4 : 3) {
0107   produces<reco::PixelClusterTagInfoCollection>();
0108 }
0109
0110 PixelClusterTagInfoProducer::~PixelClusterTagInfoProducer() {}
0111
0112 // ------------ method called to produce the data  ------------
0113 void PixelClusterTagInfoProducer::produce(edm::StreamID iID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0114   // Declare produced collection
0115   auto collectionTagInfo = std::make_unique<reco::PixelClusterTagInfoCollection>();
0116
0117   // Open jet collection
0118   edm::Handle<edm::View<reco::Jet> > collectionJets;
0119   iEvent.getByToken(m_jets, collectionJets);
0120
0121   // Get primary vertex in the event
0122   edm::Handle<reco::VertexCollection> collectionPVs;
0123   iEvent.getByToken(m_vertices, collectionPVs);
0124
0125   // Count jets above threshold
0126   int nJets(0);
0127   for (auto jetIt = collectionJets->begin(); jetIt != collectionJets->end(); ++jetIt) {
0128     if (jetIt->pt() > m_minJetPt)
0129       nJets++;
0130   }
0131
0132   // If no suitable Jet and PV is available, skip the event without opening pixel collection
0133   if (collectionPVs->empty() || nJets <= 0) {
0134     iEvent.put(std::move(collectionTagInfo));
0135     return;
0136   }
0137
0138   // Get primary vertex 3D position
0139   reco::VertexCollection::const_iterator firstPV = collectionPVs->begin();
0140   GlobalPoint v3(firstPV->x(), firstPV->y(), firstPV->z());
0141
0142   // Open Pixel Cluster collection
0143   edm::Handle<edmNew::DetSetVector<SiPixelCluster> > collectionHandle;
0144   iEvent.getByToken(m_pixelhit, collectionHandle);
0145   const edmNew::DetSetVector<SiPixelCluster>& collectionClusters(*collectionHandle);
0146
0147   // Open Geometry
0148   const TrackerGeometry& theTracker = iSetup.getData(m_geomToken);
0149
0150   // Retrieve tracker topology from geometry
0151   const TrackerTopology* tTopo = &iSetup.getData(m_topoToken);
0152
0153   std::vector<reco::PixelClusterProperties> clusters;
0154
0155   // Get vector of detunit ids, and fill a vector of PixelClusterProperties in the loop
0156   for (auto const& detUnit : collectionClusters) {
0157     if (detUnit.empty())
0158       continue;
0159     DetId detId = DetId(detUnit.detId());  // Get the Detid object for pixel detector selection
0160     if (detId.det() != DetId::Tracker)
0161       continue;
0162     if (!(detId.subdetId() == PixelSubdetector::PixelBarrel ||
0163           (m_addFPIX && detId.subdetId() == PixelSubdetector::PixelEndcap)))
0164       continue;
0165     unsigned int layer = tTopo->layer(detId);  // The layer index is in range 1-4 or 1-3
0166     if (layer == 0 || layer > m_nLayers)
0167       continue;
0168
0169     // Get the geom-detector
0170     const auto* geomDet = theTracker.idToDet(detId);
0171     const auto* topol = &geomDet->topology();
0172
0173     for (auto const& clUnit : detUnit) {
0174       if (m_minADC > 0 and clUnit.charge() < m_minADC)
0175         continue;  // skip cluster if below threshold
0176       // Get global position of the cluster
0177       LocalPoint lp = topol->localPosition(MeasurementPoint(clUnit.x(), clUnit.y()));
0178       GlobalPoint gp = geomDet->surface().toGlobal(lp);
0179       // Fill PixelClusterProperties vector for matching
0180       reco::PixelClusterProperties cp = {gp.x(), gp.y(), gp.z(), clUnit.charge(), layer};
0181       clusters.push_back(cp);
0182     }
0183   }
0184
0185   // Loop over jets and perform geometrical matching with pixel clusters
0186   for (unsigned int j = 0, nj = collectionJets->size(); j < nj; j++) {
0187     if (collectionJets->at(j).pt() < m_minJetPt)
0188       continue;
0189
0190     float cR = m_hadronMass * 2. / (collectionJets->at(j).pt());  // 2 mX / pT
0191     edm::RefToBase<reco::Jet> jetRef = collectionJets->refAt(j);  // Get jet RefToBase
0192     reco::PixelClusterData data(m_nLayers);
0193
0194     for (auto const& cluster : clusters) {
0195       GlobalPoint c3(cluster.x, cluster.y, cluster.z);  // Get cluster 3D position
0196       float dR = reco::deltaR(c3 - v3, jetRef->momentum());
0197       // Match pixel clusters to jets and fill Data struct
0198       if (cluster.layer >= 1 && cluster.layer <= m_nLayers) {
0199         int idx(cluster.layer - 1);
0200         if (dR < 0.16) {
0201           data.r016[idx]++;
0202           if (dR < 0.10) {
0203             data.r010[idx]++;
0204             if (dR < 0.08) {
0205               data.r008[idx]++;
0206               if (dR < 0.06) {
0207                 data.r006[idx]++;
0208                 if (dR < 0.04)
0209                   data.r004[idx]++;
0210               }
0211             }
0212           }
0213         }
0214         if (dR < cR) {
0215           data.rvar[idx]++;
0216           data.rvwt[idx] += cluster.charge;
0217         }
0218       }
0219     }
0220     // Create tagInfo object and fill the collection
0221     reco::PixelClusterTagInfo tagInfo(data, jetRef);
0222     collectionTagInfo->push_back(tagInfo);
0223   }
0224
0225   // Put the TagInfo collection in the event
0226   iEvent.put(std::move(collectionTagInfo));
0227 }
0228
0229 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0230 void PixelClusterTagInfoProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0231   edm::ParameterSetDescription desc;
0232   desc.add<edm::InputTag>("jets", edm::InputTag("ak4PFJetsCHS"));
0233   desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices"));
0234   desc.add<edm::InputTag>("pixelhit", edm::InputTag("siPixelClusters"));
0235   desc.add<bool>("isPhase1", true);
0236   desc.add<bool>("addForward", true);
0237   desc.add<int>("minAdcCount", -1);
0238   desc.add<double>("minJetPtCut", 100.);
0239   desc.add<double>("maxJetEtaCut", 2.5);
0240   desc.add<double>("hadronMass", 12.);
0241   descriptions.add("pixelClusterTagInfos", desc);
0242 }
0243
0244 //define this as a plug-in
0245 DEFINE_FWK_MODULE(PixelClusterTagInfoProducer);