Initial implementation

DataFormats/MuonReco/src/classes.h

@@ -20,6 +20,7 @@
 #include "DataFormats/MuonReco/interface/MuonQuality.h"
 #include "DataFormats/MuonReco/interface/MuonCosmicCompatibility.h"
 #include "DataFormats/MuonReco/interface/MuonShower.h"
+#include "DataFormats/MuonReco/interface/MuonShowerRechitCluster.h"
 #include "DataFormats/MuonReco/interface/MuonToMuonMap.h"
 #include "DataFormats/TrackReco/interface/Track.h"
 #include "DataFormats/Common/interface/AssociationMap.h"

DataFormats/MuonReco/src/classes_def.xml

@@ -164,6 +164,13 @@ initial version number of a class which has never been stored before.
   </class>
   <class name="edm::Wrapper<edm::ValueMap<reco::MuonShower> >"/>
 
+  <class name="reco::MuonShowerRechitCluster" ClassVersion="4">
+    <version ClassVersion="3" checksum="1585606091"/>
+    <version ClassVersion="4" checksum="3217667733"/>
+  </class>
+  <class name="std::vector<reco::MuonShowerRechitCluster>"/>
+  <class name="edm::Wrapper<std::vector<reco::MuonShowerRechitCluster>>"/>
+
   <class name="std::vector<reco::MuonRef>"/>
   <class name="std::vector<reco::MuonRef>::const_iterator"/>
   <class name="edm::ValueMap<reco::MuonRef>"/>

RecoMuon/MuonShowerProducer/plugins/BuildFile.xml

@@ -0,0 +1,8 @@
+<use   name="FWCore/Framework"/>
+<use   name="FWCore/PluginManager"/>
+<use   name="FWCore/ParameterSet"/>
+<use   name="Geometry/CSCGeometry"/>
+<use   name="Geometry/Records"/>
+<use   name="DataFormats/MuonReco"/>
+<use   name="fastjet"/>
+<flags EDM_PLUGIN="1"/>

RecoMuon/MuonShowerProducer/plugins/CSCRechitClusterProducer.cc

@@ -0,0 +1,221 @@
+// -*- C++ -*-
+//
+// Package:    RecoMuon/MuonShowerProducer
+// Class:      MuonShowerProducer
+// 
+/**\class MuonShowerProducer MuonShowerProducer.cc RecoMuon/MuonShowerProducer/plugins/MuonShowerProducer.cc
+
+ Description: [one line class summary]
+
+ Implementation:
+     [Notes on implementation]
+*/
+//
+// Original Author:  Ka Hei Martin Kwok
+//         Created:  Mon, 11 Oct 2021 19:52:09 GMT
+//
+//
+
+
+// system include files
+#include <memory>
+#include <cmath>
+
+// user include files
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/StreamID.h"
+
+#include "DataFormats/Common/interface/Ptr.h"
+#include "DataFormats/Math/interface/deltaPhi.h"
+#include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
+#include <DataFormats/CSCRecHit/interface/CSCRecHit2DCollection.h>
+#include "DataFormats/MuonReco/interface/MuonShowerRechitCluster.h"
+#include "Geometry/CSCGeometry/interface/CSCGeometry.h"
+#include "Geometry/Records/interface/MuonGeometryRecord.h"
+#include "DataFormats/MuonDetId/interface/CSCDetId.h"
+
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "fastjet/PseudoJet.hh"
+//
+// class declaration
+//
+
+class CSCRechitClusterProducer : public edm::stream::EDProducer<> {
+   public:
+      explicit CSCRechitClusterProducer(const edm::ParameterSet&);
+      ~CSCRechitClusterProducer();
+
+      static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
+
+   private:
+      virtual void produce(edm::Event&, const edm::EventSetup&) override;
+      edm::EDGetTokenT<CSCRecHit2DCollection> cscRechitInputToken_;
+      const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> m_cscGeometryToken;
+      typedef std::vector<reco::MuonShowerRechitCluster> MuonShowerRechitClusterCollection;
+
+      // ----------member data ---------------------------
+   protected:
+
+    std::vector<edm::Ptr<CSCRecHit2D> > inputs_;
+    std::vector<edm::Ptr<CSCRecHit2D> > getConstituents(const std::vector<fastjet::PseudoJet>& fjConstituents);
+};
+
+//
+// constants, enums and typedefs
+//
+
+//
+// static data member definitions
+//
+
+//
+// constructors and destructor
+//
+CSCRechitClusterProducer::CSCRechitClusterProducer(const edm::ParameterSet& iConfig)
+ : m_cscGeometryToken(esConsumes<CSCGeometry,MuonGeometryRecord>()) 
+{
+
+   cscRechitInputToken_ = consumes<CSCRecHit2DCollection>(edm::InputTag("csc2DRecHits")),
+
+   produces<MuonShowerRechitClusterCollection>();
+
+}
+
+
+CSCRechitClusterProducer::~CSCRechitClusterProducer(){
+}
+
+
+//
+// member functions
+//
+
+// ------------ method called to produce the data  ------------
+void
+CSCRechitClusterProducer::produce(edm::Event& ev, const edm::EventSetup& iSetup)
+{
+  edm::ESHandle<CSCGeometry> cscG;
+  cscG = iSetup.getHandle(m_cscGeometryToken);
+
+  edm::Handle<CSCRecHit2DCollection> cscRechits;
+
+  ev.getByToken(cscRechitInputToken_,cscRechits);
+
+  std::vector<fastjet::PseudoJet> fjInput;
+
+  fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, 0.4);
+
+  int nRecHits = cscRechits->size();
+  CSCRecHit2DCollection::const_iterator recIt;
+
+//  for (size_t i=0; i<  cscRechits->size(); ++i) {
+//    inputs_.push_back( (*cscRechits)[i] );
+//   }
+  for (recIt = cscRechits->begin(); recIt != cscRechits->end(); recIt++) {
+    auto cscRechit = (*recIt);
+    LocalPoint  cscRecHitLocalPosition       = cscRechit.localPosition();
+    CSCDetId cscdetid = cscRechit.cscDetId();
+    const CSCChamber* cscchamber = cscG->chamber(cscdetid);
+    if (cscchamber) {
+        GlobalPoint globalPosition = cscchamber->toGlobal(cscRecHitLocalPosition);
+        float x = globalPosition.x();
+        float y = globalPosition.y();
+        float z = globalPosition.z();
+        //TODO: Find a way to construct edm::Ptr
+        //edm::Ptr<CSCRecHit2D> ptr = edm::Ptr<CSCRecHit2D>(&cscRechit);
+        //inputs_.push_back(  ptr );
+        fjInput.push_back(fastjet::PseudoJet( x, y, z, 1.0));
+        fjInput.back().set_user_index(recIt-cscRechits->begin());
+    }
+  }
+  fastjet::ClusterSequence clus_seq ( fjInput , jet_def);
+
+  double ptmin = 0.0;
+  std::vector<fastjet::PseudoJet> fjJets = fastjet::sorted_by_pt(clus_seq.inclusive_jets(ptmin));
+
+  auto CSCclusters = std::make_unique<MuonShowerRechitClusterCollection>(fjJets.size());
+  if (!fjJets.empty()){
+    for (unsigned int ijet = 0; ijet < fjJets.size(); ++ijet) {
+
+      // get the fastjet jet
+      const fastjet::PseudoJet& fjJet = fjJets[ijet];
+      // get the constituents from fastjet
+      std::vector<fastjet::PseudoJet> const& fjConstituents = fastjet::sorted_by_pt(fjJet.constituents());
+
+      std::vector<edm::Ptr<CSCRecHit2D> > const& rechits = getConstituents(fjConstituents);
+
+      //Derive cluster properties
+      float time=0;float eta=0;float phi=0;
+      float x=0;float y=0;float z=0;
+      int nME11_12 =0;
+      for (auto & rechit : rechits){
+
+        LocalPoint  cscRecHitLocalPosition       = (*rechit).localPosition();
+        CSCDetId cscdetid = (*rechit).cscDetId();
+        const CSCChamber* cscchamber = cscG->chamber(cscdetid());
+        int endcap = CSCDetId::endcap(cscdetid) == 1 ? 1 : -1;
+        GlobalPoint globalPosition = cscchamber->toGlobal(cscRecHitLocalPosition);
+        x += globalPosition.x();
+        y += globalPosition.y();
+        z += globalPosition.z();
+
+        int  chamber =  endcap * (CSCDetId::station(cscdetid)*10 + CSCDetId::ring(cscdetid));
+        if( abs(chamber)==11 || abs(chamber)==12) nME11_12++;
+        time+= (*rechit).tpeak() + (*rechit).wireTime();
+      } 
+      x = x/rechits.size(); 
+      y = y/rechits.size(); 
+      z = z/rechits.size(); 
+      phi = std::atan(y/x);
+      if  (x < 0.0){
+        phi = M_PI + phi;
+      }
+      phi = deltaPhi(phi,0.0);
+      eta = asinhf( z / std::sqrt(x * x + y* y)); 
+      time = time/(2*rechits.size()); 
+
+      reco::MuonShowerRechitCluster cls(eta,phi,x,y,z,time,rechits.size(),nME11_12);
+      //for (auto & rechit : rechits){
+      //  cls.addDaughter(rechit);
+      //}
+
+      CSCclusters->push_back(cls);
+    }
+  }
+  ev.put(std::move(CSCclusters));
+
+}
+
+std::vector<edm::Ptr<CSCRecHit2D> > CSCRechitClusterProducer::getConstituents(const std::vector<fastjet::PseudoJet>& fjConstituents) {
+
+  std::vector<edm::Ptr<CSCRecHit2D>> result;
+  result.reserve(fjConstituents.size() / 2);
+  for (unsigned int i = 0; i < fjConstituents.size(); i++) {
+    auto index = fjConstituents[i].user_index();
+    if (index >= 0 && static_cast<unsigned int>(index) < inputs_.size()) {
+      result.emplace_back(inputs_[index]);
+    }
+  }
+  return result;
+}
+
+// ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
+void
+CSCRechitClusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  //The following says we do not know what parameters are allowed so do no validation
+  // Please change this to state exactly what you do use, even if it is no parameters
+  edm::ParameterSetDescription desc;
+  desc.setUnknown();
+  descriptions.addDefault(desc);
+}
+
+//define this as a plug-in
+DEFINE_FWK_MODULE(CSCRechitClusterProducer);

RecoMuon/MuonShowerProducer/python/cscRechitCluster_cfi.py

@@ -0,0 +1,3 @@
+import FWCore.ParameterSet.Config as cms
+
+cscRechitClusters = cms.EDProducer("CSCRechitClusterProducer")