From cmssw 12 6 0 pre2

DPGAnalysis/CaloNanoAOD/python/mergedSimClusters_cff.py

@@ -3,21 +3,22 @@
 from DPGAnalysis.CaloNanoAOD.simClusters_cff import *
 from DPGAnalysis.HGCalNanoAOD.hgcRecHits_cff import *
 
+
+closeBySimClusters = cms.EDProducer("CPtoSimClusters",
+                                    caloParticles = cms.InputTag("mix:MergedCaloTruth"),
+                                    simVertices = cms.InputTag("g4SimHits")
+                                    )
+
 hgcSimTruth = cms.EDProducer("SimClusterMerger",
 
-    useNLayers = cms.int32(2),
-    searchRadiusScale = cms.double(2.),
-    clusterRadiusScale = cms.double(1.),
-
-    mergeRadiusScale = cms.double(7.),#13 is about 10 layers in CE
-    energyContainment = cms.double(1.1),
+    ## only these three actually matter atm
 
-    smear = cms.double(-0.0),
-    highEfracThreshold = cms.double(0.85), 
-    connectThreshold = cms.double(.3),
-
-    relOverlapDistance = cms.double(.9),# dist/(merged radius + sensor radius)
+    useNLayers = cms.int32(10),
+    highEfracThreshold = cms.double(0.95), 
+    connectThreshold = cms.double(.1),
 
+    ###
+
     simClusters= cms.InputTag("mix:MergedCaloTruth"),
     simVertices= cms.InputTag("g4SimHits"),
     simTracks= cms.InputTag("g4SimHits"),
@@ -75,7 +76,8 @@
     docString = cms.string("SimCluster deposited reconstructed energy associated to SimCluster")
 )
 
-mergedSimClusterTables = cms.Sequence(hgcSimTruth+
+mergedSimClusterTables = cms.Sequence(closeBySimClusters+
+    hgcSimTruth+
 	mergedSimClusterTable+
 	mergedToUnmergedSCTable+
 	hgcRecHitsToMergedSimClusters+

DPGAnalysis/HGCalNanoAOD/plugins/HGCALHitPositionTableProducer.cc

@@ -18,8 +18,8 @@ class HGCalHitPositionTableProducer : public HitPositionTableProducer<edm::View<
 //class HGCalHitPositionTableProducer : public HitPositionTableProducer<edm::View<CaloRecHit>> {
 public:
   HGCalHitPositionTableProducer(edm::ParameterSet const& params)
-      : HitPositionTableProducer<edm::View<T>>(params), 
-        caloGeoToken_(edm::stream::EDProducer<>::esConsumes<edm::Transition::BeginRun>()) {}
+      : HitPositionTableProducer<edm::View<T>>(params),
+        caloGeoToken_(edm::stream::EDProducer<>::esConsumes<edm::Transition::BeginRun>()) { }
 
   ~HGCalHitPositionTableProducer() override {}
 
@@ -33,6 +33,10 @@ class HGCalHitPositionTableProducer : public HitPositionTableProducer<edm::View<
     return positionFromDetId(detId); 
   }
 
+  float radiusFromHit(const CaloRecHit& hit) { return radiusFromDetId(hit.detid()); }
+
+  float radiusFromHit(const PCaloHit& hit) { return radiusFromDetId(hit.id()); }
+
   void beginRun(const edm::Run&, const edm::EventSetup& iSetup) override {
     auto& geom = iSetup.getData(caloGeoToken_);
     rhtools_.setGeometry(geom);
@@ -47,6 +51,28 @@ class HGCalHitPositionTableProducer : public HitPositionTableProducer<edm::View<
     }
   }
 
+  float radiusFromDetId(DetId id) {
+    DetId::Detector det = id.det();
+    if (det == DetId::HGCalEE || det == DetId::HGCalHSi || det == DetId::HGCalHSc) {
+      if (rhtools_.isSilicon(id)) {
+        return rhtools_.getRadiusToSide(id);
+      } else if (rhtools_.isScintillator(id)) {
+          //return r * sin(dphi) / 2.;
+       // auto detadphi = rhtools_.getScintDEtaDPhi(id);
+       // float dphi = detadphi.second;  //same in both
+        float dphi = 0.0211809f; //DEBUG FIXME; the geometry access segfaults for some hits probably close to side
+        auto pos = rhtools_.getPosition(id);
+        float r = pos.transverse();
+        return r * sin(dphi) / 2.;  //this is anyway approximate
+      } else {
+        return 0.;
+      }
+
+    } else {
+      return 0;  //no except here
+    }
+  }
+
 protected:
   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeoToken_;
   hgcal::RecHitTools rhtools_;

HGCSimTruth/HGCSimTruth/plugins/CPtoSimClusters.cc

@@ -0,0 +1,87 @@
+#include <memory>
+#include <vector>
+#include <cstdlib>
+#include <iostream>
+#include <stack>
+#include <unordered_map>
+#include <map>
+#include <sstream>
+#include <utility>
+#include <set>
+#include <cmath>
+#include <numeric>
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "DataFormats/Common/interface/Ref.h"
+
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticle.h"
+#include "SimDataFormats/CaloAnalysis/interface/CaloParticleFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimClusterFwd.h"
+#include "SimDataFormats/CaloAnalysis/interface/SimCluster.h"
+#include "SimDataFormats/Vertex/interface/SimVertex.h"
+
+
+class CPtoSimClusters : public edm::stream::EDProducer<> {
+    public:
+        explicit CPtoSimClusters(const edm::ParameterSet&);
+        ~CPtoSimClusters() {}
+
+
+    private:
+        void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
+
+        edm::EDGetTokenT<CaloParticleCollection> cpCollectionToken_;
+        edm::EDGetTokenT<std::vector<SimVertex> > svCollectionToken_;
+
+    };
+
+CPtoSimClusters::CPtoSimClusters(const edm::ParameterSet &pset) :
+              cpCollectionToken_(consumes<CaloParticleCollection>(pset.getParameter<edm::InputTag>("caloParticles"))),
+              svCollectionToken_(consumes<std::vector<SimVertex> >(pset.getParameter<edm::InputTag>("simVertices")))
+{
+    produces<SimClusterCollection>();
+}
+
+void CPtoSimClusters::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+
+
+    edm::Handle<CaloParticleCollection> cpCollection;
+    iEvent.getByToken(cpCollectionToken_, cpCollection);
+
+    edm::Handle<std::vector<SimVertex> > svCollection;
+    iEvent.getByToken(svCollectionToken_, svCollection);
+
+    auto output = std::make_unique<SimClusterCollection>();
+
+    for(const auto& cp: *cpCollection){
+        SimCluster cpsc;
+        bool first = true;
+        for(const auto& sc: cp.simClusters()){
+            if(first){
+                cpsc = *sc;
+                first = false;
+            }
+            else{
+                cpsc += *sc;
+            }
+        }
+
+        cpsc.setImpactMomentum(cp.p4());
+        auto vertex = svCollection->at(cp.g4Tracks().at(0).vertIndex());
+        cpsc.setImpactPoint(math::XYZTLorentzVectorF(vertex.position()));
+        cpsc.setPdgId(cp.pdgId());//done
+        //cp.
+        //cp.
+        output->push_back(cpsc);
+    }
+
+    iEvent.put(std::move(output));
+
+}
+
+DEFINE_FWK_MODULE(CPtoSimClusters);

HGCSimTruth/HGCSimTruth/plugins/SimClusterMerger.cc

@@ -75,14 +75,6 @@ class SimClusterMerger : public edm::stream::EDProducer<> {
 
         //merge config
         int cNLayers_;
-        float cEContainment_;
-        float cSearchRadius_;
-        float cClusterRadiusScale_;
-        float cMergeRadiusScale_;
-
-        float cMergeThreshold_;
-
-        float cSmear_;
         float cIsHighEfracThreshold_;
         float cConnectThreshold_;
     };
@@ -100,13 +92,6 @@ SimClusterMerger::SimClusterMerger(const edm::ParameterSet &pset) :
 
     //move to initialisers
     cNLayers_ = pset.getParameter<int32_t> ( "useNLayers" );
-    cSearchRadius_ = pset.getParameter<double> ( "searchRadiusScale" );
-    cClusterRadiusScale_ = pset.getParameter<double> ( "clusterRadiusScale" );
-    cMergeRadiusScale_ = pset.getParameter<double> ( "mergeRadiusScale" );
-    cEContainment_ = pset.getParameter<double> ( "energyContainment" );
-    cMergeThreshold_ = pset.getParameter<double> ( "relOverlapDistance" );
-
-    cSmear_ = pset.getParameter<double> ( "smear" );
     cIsHighEfracThreshold_ = pset.getParameter<double> ( "highEfracThreshold" );
     cConnectThreshold_ = pset.getParameter<double> ( "connectThreshold" );
 
@@ -140,20 +125,25 @@ void SimClusterMerger::produce(edm::Event& iEvent, const edm::EventSetup& iSetup
     HGCalSimClusterMerger merger(*caloRecHitCollection.product(),
             &hgcalRecHitToolInstance_, &histtool);
 
-    merger.setCClusterRadiusScale(cClusterRadiusScale_);
-    merger.setCEContainment(cEContainment_);
-    merger.setCMergeRadiusScale(cMergeRadiusScale_);
+
     merger.setCNLayers(cNLayers_);
-    merger.setCSearchRadius(cSearchRadius_);
+    merger.setHighEfracThreshold(cIsHighEfracThreshold_);
+    merger.setConnectThreshold(cConnectThreshold_);
 
     std::vector<const SimCluster*> tobemerged;
     for(size_t i = 0; i < scCollection->size(); i++)
         tobemerged.push_back(&scCollection->at(i));
 
+    //reduce muon impact momenta
+    auto sc_tobemerged = merger.downScaleMuons(tobemerged);
+    tobemerged.clear();
+    for(size_t i = 0; i < scCollection->size(); i++)
+        tobemerged.push_back(&sc_tobemerged.at(i));
+
     //cMergeThreshold_
     //FIXME remove hard coded
     std::vector<std::vector<size_t> > idxs; //// <--- @Shah Rukh: that's the one to check if something is merged
-    auto mergedSC = merger.merge(tobemerged, cIsHighEfracThreshold_, cConnectThreshold_, cSmear_, idxs,true); //now idxs is filled
+    auto mergedSC = merger.merge(tobemerged, idxs); //now idxs is filled
 
     std::cout << "SCmerger initial: " << scCollection->size() << " merged: " << idxs.size() << std::endl;//DEBUG
     size_t largestgroup=0;

IOMC/ParticleGuns/interface/FlatEtaRangeNoTrackerGunProducer.h

@@ -13,7 +13,7 @@
 #include "HepMC/GenEvent.h"
 #include "HepPDT/ParticleDataTable.hh"
 
-#include "RecoHGCal/GraphReco/interface/HGCalParticlePropagator.h"
+#include "RecoHGCal/GraphReco/interface/HGCalTrackPropagator.h"
 
 
 namespace edm {
@@ -60,14 +60,8 @@ namespace edm {
     // debug flag
     bool debug_;
 
-    // pointer to the current event
-    HepMC::GenEvent* genEvent_;
 
-    // pdg table
-    ESHandle<HepPDT::ParticleDataTable> pdgTable_;
 
-    //propagator
-    HGCalParticlePropagator prop_;
 
     //smear vertex
     double timeSmear_;
@@ -76,6 +70,14 @@ namespace edm {
     double momSmear_;
     double minDistDR_;
 
+    // pdg table
+    const ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> pdgTableToken_;
+    ESHandle<HepPDT::ParticleDataTable> pdgTable_;
+    // pointer to the current event
+    HepMC::GenEvent* genEvent_;
+    //propagator
+    HGCalObjectPropagator<int> prop_; //just dummy template argument - not used
+
   };
 
 }  // namespace edm

IOMC/ParticleGuns/src/FlatEtaRangeNoTrackerGunProducer.cc

@@ -61,10 +61,12 @@ edm::FlatEtaRangeNoTrackerGunProducer::FlatEtaRangeNoTrackerGunProducer(const ed
       phiMin_(params.getParameter<double>("phiMin")),
       phiMax_(params.getParameter<double>("phiMax")),
       debug_(params.getUntrackedParameter<bool>("debug")),
-      genEvent_(nullptr),
       timeSmear_(params.getParameter<double>("timeSmearInPs")*1e-12),
       momSmear_(params.getParameter<double>("momSmear")),
-      minDistDR_(params.getParameter<double>("minDistDR")){
+      minDistDR_(params.getParameter<double>("minDistDR")),
+      pdgTableToken_(esConsumes<Transition::BeginRun>()),
+      genEvent_(nullptr),
+      prop_(consumesCollector()){
   produces<edm::HepMCProduct>("unsmeared");
   produces<GenEventInfoProduct>();
   produces<GenRunInfoProduct, edm::Transition::EndRun>();
@@ -73,8 +75,8 @@ edm::FlatEtaRangeNoTrackerGunProducer::FlatEtaRangeNoTrackerGunProducer(const ed
 edm::FlatEtaRangeNoTrackerGunProducer::~FlatEtaRangeNoTrackerGunProducer() {}
 
 void edm::FlatEtaRangeNoTrackerGunProducer::beginRun(const edm::Run& run, const edm::EventSetup& setup) {
-  setup.getData(pdgTable_);
-  prop_.setEventSetup(setup);
+  pdgTable_ = setup.getHandle(pdgTableToken_);
+  prop_.setupRun(setup);
 }
 
 void edm::FlatEtaRangeNoTrackerGunProducer::endRun(const edm::Run& run, const edm::EventSetup& setup) {}
@@ -143,13 +145,11 @@ void edm::FlatEtaRangeNoTrackerGunProducer::produce(edm::Event& event, const edm
         if(debug_)
             LogDebug("FlatEtaRangeNoTrackerGunProducer") << " : Initial vertex position " << vertexPos << std::endl;
         //propagate
-        prop_.propagate(vertexPos,momentum,(pData->charge() > 0) ? 1 : ((pData->charge() < 0) ? -1 : 0));
+        prop_.propagateVectors(vertexPos,momentum,(pData->charge() > 0) ? 1 : ((pData->charge() < 0) ? -1 : 0));
+
+        if(debug_)
+            std::cout << "FlatEtaRangeNoTrackerGunProducer" << " : final vertex position " << vertexPos << std::endl;
 
-        //just move a tad (1mm) towards beamspot
-        if(vertexPos.z()>0)
-            vertexPos.SetXYZT(vertexPos.x(),vertexPos.y(),vertexPos.z()-0.1,vertexPos.t());
-        else
-            vertexPos.SetXYZT(vertexPos.x(),vertexPos.y(),vertexPos.z()+0.1,vertexPos.t());
 
         if(minDistDR_>0){
             bool next=false;

PhysicsTools/NanoAOD/interface/HitPositionTableProducer.h

@@ -28,6 +28,10 @@ class HitPositionTableProducer : public edm::stream::EDProducer<> {
     throw cms::Exception("HitPositionTableProducer") << "Virtual function positionFromHist is not implemented!"; 
   }
 
+  float virtual radiusFromHit(const typename T::value_type& hit) {
+    throw cms::Exception("HitPositionTableProducer") << "Virtual function radiusFromHit is not implemented!";
+  }
+
   void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override {
     edm::Handle<T> objs;
     iEvent.getByToken(src_, objs);
@@ -42,7 +46,7 @@ class HitPositionTableProducer : public edm::stream::EDProducer<> {
         xvals.emplace_back(position.x());
         yvals.emplace_back(position.y());
         zvals.emplace_back(position.z());
-        //hitrvals.emplace_back(radiusFromHit(obj));
+        hitrvals.emplace_back(radiusFromHit(obj));
       }
     }
 
@@ -51,7 +55,7 @@ class HitPositionTableProducer : public edm::stream::EDProducer<> {
     tab->addColumn<float>("y", yvals, "y position");
     tab->addColumn<float>("z", zvals, "z position");
     //addExtraColumns(tab);
-    //tab->addColumn<float>("hitr", hitrvals, "radius");
+    tab->addColumn<float>("hitr", hitrvals, "radius");
 
     iEvent.put(std::move(tab));
   }