MuonRechitClusterProducer/plugins/RechitClusterProducer.cc

0001 #include <map>
0002 #include <memory>
0003 #include <string>
0004 #include <utility>
0005 #include <vector>
0006 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0007 #include "DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h"
0008 #include "DataFormats/Common/interface/RangeMap.h"
0009 #include "DataFormats/Common/interface/Ref.h"
0010 #include "DataFormats/Common/interface/RefVector.h"
0011 #include "DataFormats/Common/interface/RefVectorIterator.h"
0012 #include "DataFormats/DTRecHit/interface/DTRecHit1DPair.h"
0013 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0017 #include "DataFormats/Math/interface/PtEtaPhiMass.h"
0018 #include "DataFormats/Math/interface/Vector3D.h"
0019 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0020 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0021 #include "DataFormats/MuonReco/interface/MuonRecHitCluster.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/Framework/interface/Event.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Framework/interface/global/EDProducer.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSetDescriptionFiller.h"
0029 #include "FWCore/Utilities/interface/EDGetToken.h"
0030 #include "FWCore/Utilities/interface/ESGetToken.h"
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032 #include "FWCore/Utilities/interface/StreamID.h"
0033 #include "FWCore/Utilities/interface/typelookup.h"
0034 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0035 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0036 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0037 #include "fastjet/ClusterSequence.hh"
0038 #include "fastjet/JetDefinition.hh"
0039 #include "fastjet/PseudoJet.hh"
0040
0041 typedef std::vector<reco::MuonRecHitCluster> RecHitClusterCollection;
0042
0043 template <typename Trait>
0044 class RechitClusterProducerT : public edm::global::EDProducer<> {
0045   typedef typename Trait::RecHitRef RechitRef;
0046   typedef typename Trait::RecHitRefVector RecHitRefVector;
0047
0048 public:
0049   explicit RechitClusterProducerT(const edm::ParameterSet&);
0050   ~RechitClusterProducerT() override = default;
0051
0052   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0053
0054 private:
0055   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0056   const edm::ESGetToken<typename Trait::GeometryType, MuonGeometryRecord> geometryToken_;
0057   edm::EDGetTokenT<typename Trait::InputType> inputToken_;
0058   typedef std::vector<reco::MuonRecHitCluster> RecHitClusterCollection;
0059
0060   const double rParam_;      // distance paramter
0061   const int nRechitMin_;     // min number of rechits
0062   const int nStationThres_;  // min number of rechits to count towards nStation
0063 };
0064
0065 template <typename Trait>
0066 RechitClusterProducerT<Trait>::RechitClusterProducerT(const edm::ParameterSet& iConfig)
0067     : geometryToken_(esConsumes<typename Trait::GeometryType, MuonGeometryRecord>()),
0068       inputToken_(consumes<typename Trait::InputType>(iConfig.getParameter<edm::InputTag>("recHitLabel"))),
0069       rParam_(iConfig.getParameter<double>("rParam")),
0070       nRechitMin_(iConfig.getParameter<int>("nRechitMin")),
0071       nStationThres_(iConfig.getParameter<int>("nStationThres")) {
0072   produces<RecHitClusterCollection>();
0073 }
0074
0075 //
0076 // member functions
0077 //
0078
0079 // ------------ method called to produce the data  ------------
0080 template <typename Trait>
0081 void RechitClusterProducerT<Trait>::produce(edm::StreamID, edm::Event& ev, const edm::EventSetup& iSetup) const {
0082   auto const& geo = iSetup.getData(geometryToken_);
0083   auto const& rechits = ev.get(inputToken_);
0084
0085   std::vector<fastjet::PseudoJet> fjInput;
0086   RecHitRefVector inputs;
0087
0088   fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, rParam_);
0089
0090   Trait recHitTrait;
0091
0092   int recIt = 0;
0093   for (auto const& rechit : rechits) {
0094     LocalPoint recHitLocalPosition = rechit.localPosition();
0095     auto detid = recHitTrait.detid(rechit);
0096     auto thischamber = geo.chamber(detid);
0097     if (thischamber) {
0098       GlobalPoint globalPosition = thischamber->toGlobal(recHitLocalPosition);
0099       float x = globalPosition.x();
0100       float y = globalPosition.y();
0101       float z = globalPosition.z();
0102       RechitRef ref = RechitRef(&rechits, recIt);
0103       inputs.push_back(ref);
0104       fjInput.push_back(fastjet::PseudoJet(x, y, z, globalPosition.mag()));
0105       fjInput.back().set_user_index(recIt);
0106     }
0107     recIt++;
0108   }
0109   fastjet::ClusterSequence clus_seq(fjInput, jet_def);
0110
0111   //keep all the clusters
0112   double ptmin = 0.0;
0113   std::vector<fastjet::PseudoJet> fjJets = clus_seq.inclusive_jets(ptmin);
0114
0115   auto clusters = std::make_unique<RecHitClusterCollection>();
0116   for (auto const& fjJet : fjJets) {
0117     // skip if the cluster has too few rechits
0118     if (int(fjJet.constituents().size()) < nRechitMin_)
0119       continue;
0120     // get the constituents from fastjet
0121     RecHitRefVector rechits;
0122     for (auto const& constituent : fjJet.constituents()) {
0123       auto index = constituent.user_index();
0124       if (index >= 0 && static_cast<unsigned int>(index) < inputs.size()) {
0125         rechits.push_back(inputs[index]);
0126       }
0127     }
0128
0129     //Derive cluster properties
0130     int nStation = 0;
0131     int totStation = 0;
0132     float avgStation = 0.0;
0133     std::map<int, int> station_count_map;
0134     for (auto const& rechit : rechits) {
0135       station_count_map[recHitTrait.station(*rechit)]++;
0136     }
0137     //station statistics
0138     std::map<int, int>::iterator it;
0139     for (auto const& [station, count] : station_count_map) {
0140       if (count >= nStationThres_) {
0141         nStation++;
0142         avgStation += station * count;
0143         totStation += count;
0144       }
0145     }
0146     if (totStation != 0) {
0147       avgStation = avgStation / totStation;
0148     }
0149     float invN = 1.f / rechits.size();
0150     // cluster position is the average position of the constituent rechits
0151     float jetX = fjJet.px() * invN;
0152     float jetY = fjJet.py() * invN;
0153     float jetZ = fjJet.pz() * invN;
0154
0155     math::RhoEtaPhiVectorF position(
0156         std::sqrt(jetX * jetX + jetY * jetY), etaFromXYZ(jetX, jetY, jetZ), std::atan2(jetY, jetX));
0157     Trait::emplace_back(clusters.get(), position, nStation, avgStation, rechits);
0158   }
0159   ev.put(std::move(clusters));
0160 }
0161
0162 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0163 template <typename Trait>
0164 void RechitClusterProducerT<Trait>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0165   edm::ParameterSetDescription desc;
0166   desc.add<int>("nRechitMin", 50);
0167   desc.add<double>("rParam", 0.4);
0168   desc.add<int>("nStationThres", 10);
0169   desc.add<edm::InputTag>("recHitLabel", edm::InputTag(Trait::recHitLabel()));
0170   descriptions.add(Trait::producerName(), desc);
0171 }
0172 struct DTRecHitTrait {
0173   using GeometryType = DTGeometry;
0174   using InputType = DTRecHitCollection;
0175   using RecHitRef = edm::Ref<DTRecHitCollection>;
0176   using RecHitRefVector = edm::RefVector<DTRecHitCollection>;
0177   static std::string recHitLabel() { return "dt1DRecHits"; }
0178   static std::string producerName() { return "dtRechitClusterProducer"; }
0179
0180   static int station(const DTRecHit1DPair& dtRechit) { return detid(dtRechit).station(); }
0181   static DTChamberId detid(const DTRecHit1DPair& dtRechit) {
0182     DetId geoid = dtRechit.geographicalId();
0183     DTChamberId dtdetid = DTChamberId(geoid);
0184     return dtdetid;
0185   }
0186   static void emplace_back(RecHitClusterCollection* clusters,
0187                            math::RhoEtaPhiVectorF const& position,
0188                            int nStation,
0189                            float avgStation,
0190                            RecHitRefVector const& rechits) {
0191     // compute nMB1, nMB2
0192     int nMB1 = 0;
0193     int nMB2 = 0;
0194     for (auto const& rechit : rechits) {
0195       DetId geoid = rechit->geographicalId();
0196       DTChamberId dtdetid = DTChamberId(geoid);
0197       if (dtdetid.station() == 1)
0198         nMB1++;
0199       if (dtdetid.station() == 2)
0200         nMB2++;
0201     }
0202     //set time, timespread, nME11,nME12 to 0
0203     reco::MuonRecHitCluster cls(position, rechits.size(), nStation, avgStation, 0.0, 0.0, 0, 0, 0, 0, nMB1, nMB2);
0204     clusters->emplace_back(cls);
0205   }
0206 };
0207
0208 struct CSCRecHitTrait {
0209   using GeometryType = CSCGeometry;
0210   using InputType = CSCRecHit2DCollection;
0211   using RecHitRef = edm::Ref<CSCRecHit2DCollection>;
0212   using RecHitRefVector = edm::RefVector<CSCRecHit2DCollection>;
0213   static std::string recHitLabel() { return "csc2DRecHits"; }
0214   static std::string producerName() { return "cscRechitClusterProducer"; }
0215
0216   static int station(const CSCRecHit2D& cscRechit) { return CSCDetId::station(detid(cscRechit)); }
0217   static CSCDetId detid(const CSCRecHit2D& cscRechit) { return cscRechit.cscDetId(); }
0218   static void emplace_back(RecHitClusterCollection* clusters,
0219                            math::RhoEtaPhiVectorF const& position,
0220                            int nStation,
0221                            float avgStation,
0222                            RecHitRefVector const& rechits) {
0223     int nME11 = 0;
0224     int nME12 = 0;
0225     int nME41 = 0;
0226     int nME42 = 0;
0227     float timeSpread = 0.0;
0228     float time = 0.0;
0229     float time_strip = 0.0;  // for timeSpread calculation
0230     for (auto const& rechit : rechits) {
0231       CSCDetId cscdetid = rechit->cscDetId();
0232       int stationRing = (CSCDetId::station(cscdetid) * 10 + CSCDetId::ring(cscdetid));
0233       if (CSCDetId::ring(cscdetid) == 4)
0234         stationRing = (CSCDetId::station(cscdetid) * 10 + 1);  // ME1/a has ring==4
0235       if (stationRing == 11)
0236         nME11++;
0237       if (stationRing == 12)
0238         nME12++;
0239       if (stationRing == 41)
0240         nME41++;
0241       if (stationRing == 42)
0242         nME42++;
0243       time += (rechit->tpeak() + rechit->wireTime());
0244       time_strip += rechit->tpeak();
0245     }
0246     float invN = 1.f / rechits.size();
0247     time = (time / 2.f) * invN;
0248     time_strip = time_strip * invN;
0249
0250     //derive cluster statistics
0251     for (auto& rechit : rechits) {
0252       timeSpread += (rechit->tpeak() - time_strip) * (rechit->tpeak() - time_strip);
0253     }
0254     timeSpread = std::sqrt(timeSpread * invN);
0255
0256     //set nMB1,nMB2 to 0
0257     reco::MuonRecHitCluster cls(
0258         position, rechits.size(), nStation, avgStation, time, timeSpread, nME11, nME12, nME41, nME42, 0, 0);
0259     clusters->emplace_back(cls);
0260   }
0261 };
0262
0263 using DTrechitClusterProducer = RechitClusterProducerT<DTRecHitTrait>;
0264 using CSCrechitClusterProducer = RechitClusterProducerT<CSCRecHitTrait>;
0265 DEFINE_FWK_MODULE(DTrechitClusterProducer);
0266 DEFINE_FWK_MODULE(CSCrechitClusterProducer);