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