cms-pepr · kdlong · Sep 9, 2020 · Sep 29, 2020 · Sep 29, 2020 · Sep 30, 2020
diff --git a/DataFormats/HGCalReco/interface/TICLLayerTile.h b/DataFormats/HGCalReco/interface/TICLLayerTile.h
@@ -35,6 +35,8 @@ class TICLLayerTileT {
 
   int globalBin(double eta, double phi) const { return phiBin(phi) + etaBin(eta) * T::nPhiBins; }
 
+  int nBins() const { return T::nBins; }
+
   void clear() {
     auto nBins = T::nEtaBins * T::nPhiBins;
     for (int j = 0; j < nBins; ++j)

diff --git a/RecoHGCal/GraphReco/BuildFile.xml b/RecoHGCal/GraphReco/BuildFile.xml
@@ -13,4 +13,4 @@
   <use   name="TrackPropagation/RungeKutta"/>
 <export>
   <lib   name="1"/>
-</export>
+</export>
diff --git a/RecoHGCal/GraphReco/plugins/BuildFile.xml b/RecoHGCal/GraphReco/plugins/BuildFile.xml
@@ -2,9 +2,12 @@
 <use name="FWCore/PluginManager" />
 <use name="FWCore/ParameterSet" />
 <use name="FWCore/MessageLogger" />
+<use name="DataFormats/Common" />
 <use name="DataFormats/TrackReco" />
-<use name="RecoLocalCalo/HGCalRecAlgos" />
+<use name="DataFormats/ParticleFlowCandidate" />
+<use name="DataFormats/HGCalReco" />
 <use name="PhysicsTools/TensorFlow" />
+<use name="RecoLocalCalo/HGCalRecAlgos" />
 <use name="SimDataFormats/CaloAnalysis" />
 <use name="HGCSimTruth/HGCSimTruth" />
 <use name="root"/>

diff --git a/RecoHGCal/GraphReco/plugins/HGCClusterTrackLinker.cc b/RecoHGCal/GraphReco/plugins/HGCClusterTrackLinker.cc
@@ -0,0 +1,102 @@
+#include <memory>
+#include <iostream>
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/stream/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "DataFormats/TrackReco/interface/Track.h"
+#include "DataFormats/TrackReco/interface/TrackFwd.h"
+
+#include "DataFormats/Common/interface/View.h"
+
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
+#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h" 
+
+#include "DataFormats/HGCalReco/interface/TICLLayerTile.h"
+
+class HGCClusterTrackLinker : public edm::stream::EDProducer<> {
+public:
+    explicit HGCClusterTrackLinker(const edm::ParameterSet&);
+    ~HGCClusterTrackLinker() {};
+
+ private:
+    void produce(edm::Event&, const edm::EventSetup&) override;
+
+    edm::EDGetTokenT<edm::View<reco::Track>> tracksToken_;
+    edm::EDGetTokenT<reco::PFCandidateCollection> pfCandsToken_;
+
+};
+
+HGCClusterTrackLinker::HGCClusterTrackLinker(const edm::ParameterSet& config) :
+    tracksToken_(consumes<edm::View<reco::Track>>(config.getParameter<edm::InputTag>("tracks"))),
+    pfCandsToken_(consumes<reco::PFCandidateCollection>(config.getParameter<edm::InputTag>("pfCands"))) {
+  produces<reco::PFCandidateCollection>();
+}
+
+void HGCClusterTrackLinker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
+  edm::Handle<reco::PFCandidateCollection> pfCands;
+  iEvent.getByToken(pfCandsToken_, pfCands);
+
+  edm::Handle<edm::View<reco::Track>> tracks;
+  iEvent.getByToken(tracksToken_, tracks);
+
+  ticl::TICLLayerTile candTile;
+  ticl::TICLLayerTile tracksTile;
+
+  auto out = std::make_unique<reco::PFCandidateCollection>();
+
+  for (size_t i = 0; i < tracks->size(); i++) {
+    edm::RefToBase<reco::Track> track(tracks, i);
+    if (std::abs(track->eta()) > ticl::TileConstants::minEta && std::abs(track->eta()) < ticl::TileConstants::maxEta)
+      tracksTile.fill(track->eta(), track->phi(), i);
+  }
+
+  for (size_t ip = 0; ip < pfCands->size(); ip++) {
+    reco::PFCandidateRef pfCand(pfCands, ip);
+    if (pfCand->charge())
+        candTile.fill(pfCand->eta(), pfCand->phi(), ip);
+  }
+
+  for (size_t bin = 0; bin < static_cast<size_t>(candTile.nBins()); bin++) {
+    const std::vector<unsigned int>& candIndices = candTile[bin]; 
+    if (!candIndices.size())
+        continue;
+
+    std::vector<unsigned int> trackIndices = tracksTile[bin]; 
+    for (auto& tidx : trackIndices) {
+        edm::RefToBase<reco::Track> track(tracks, tidx);
+    }
+    for (auto& cidx : candIndices) {
+      reco::PFCandidateRef pfCand(pfCands, cidx);
+      out->push_back(*pfCand->clone());
+    }
+    // Sort now so we prioritizes matches for higher pt clusters
+    std::sort(out->begin(), out->end(), [](auto& a, auto& b) { return a.pt() > b.pt(); });
+
+    for (auto& cand : *out) {
+      if (trackIndices.size()) {
+        // Find best match. Current just closest in pt, can do 
+        // a smarter mix of closest in pt and eta/phi or something else
+        std::vector<float> ptdiffs;
+        ptdiffs.reserve(trackIndices.size());
+        for (size_t entry : trackIndices) {
+          edm::RefToBase<reco::Track> t(tracks, entry);
+          ptdiffs.push_back(std::abs(t->pt()-cand.pt()));
+        }
+
+        size_t tidx = std::distance(ptdiffs.begin(), std::min_element(ptdiffs.begin(), ptdiffs.end()));
+        edm::RefToBase<reco::Track> matchingTrack(tracks, trackIndices.at(tidx));
+        trackIndices.erase(trackIndices.begin()+tidx);
+
+        cand.setTrackRef(matchingTrack.castTo<reco::TrackRef>());
+      }
+    }
+  }
+
+  iEvent.put(std::move(out));
+}
+
+DEFINE_FWK_MODULE(HGCClusterTrackLinker);
diff --git a/SimCalorimetry/EcalElectronicsEmulation/src/EcalSimpleProducer.cc b/SimCalorimetry/EcalElectronicsEmulation/src/EcalSimpleProducer.cc
@@ -84,7 +84,7 @@ void EcalSimpleProducer::produce(edm::Event &evt, const edm::EventSetup &) {
         double em = simHitFormula_->Eval(iEta0, iPhi0, ievt - 1);
         double eh = 0.;
         double t = 0.;
-        const PCaloHit hit(EBDetId(iEta1, iPhi).rawId(), em, eh, t, 0);
+        const PCaloHit hit(EBDetId(iEta1, iPhi).rawId(), em, eh, t, 0, 0, 0);
         hits->push_back(hit);
       }
     }

diff --git a/SimDataFormats/CaloHit/interface/PCaloHit.h b/SimDataFormats/CaloHit/interface/PCaloHit.h
@@ -14,6 +14,7 @@ class PCaloHit {
       : myEnergy(e), myEMFraction(emFraction), myTime(t), myItra(i), detId(id), myDepth(d) {}
   PCaloHit(float eEM, float eHad, float t, int i = 0, uint16_t d = 0);
   PCaloHit(unsigned int id, float eEM, float eHad, float t, int i = 0, uint16_t d = 0);
+  PCaloHit(unsigned int id, float eEM, float eHad, float t, int i = 0, int iFine = 0, uint16_t d = 0);
 
   //Names
   static const char *name() { return "Hit"; }
@@ -31,6 +32,7 @@ class PCaloHit {
 
   //Geant track number
   int geantTrackId() const { return myItra; }
+  int geantFineTrackId() const { return myFineItra; }
 
   //DetId where the Hit is recorded
   void setID(unsigned int id) { detId = id; }
@@ -67,6 +69,7 @@ class PCaloHit {
   float myEMFraction;
   float myTime;
   int myItra;
+  int myFineItra;
   unsigned int detId;
   uint16_t myDepth;
   EncodedEventId theEventId;

diff --git a/SimDataFormats/CaloHit/src/PCaloHit.cc b/SimDataFormats/CaloHit/src/PCaloHit.cc
@@ -12,6 +12,12 @@ PCaloHit::PCaloHit(unsigned int id, float eEM, float eHad, float t, int i, uint1
   myEMFraction = (myEnergy <= 0. ? 1. : eEM / myEnergy);
 }
 
+PCaloHit::PCaloHit(unsigned int id, float eEM, float eHad, float t, int i, int iFine, uint16_t d)
+    : myTime(t), myItra(i), myFineItra(iFine), detId(id), myDepth(d) {
+  myEnergy = eEM + eHad;
+  myEMFraction = (myEnergy <= 0. ? 1. : eEM / myEnergy);
+}
+
 std::ostream& operator<<(std::ostream& o, const PCaloHit& hit) {
   o << "0x" << std::hex << hit.id() << std::dec << ": Energy (EM) " << hit.energyEM() << " GeV "
     << ": Energy (Had) " << hit.energyHad() << " GeV "

diff --git a/SimDataFormats/CaloHit/src/classes_def.xml b/SimDataFormats/CaloHit/src/classes_def.xml
@@ -1,5 +1,6 @@
 <lcgdict>
-  <class name="PCaloHit" ClassVersion="10">
+  <class name="PCaloHit" ClassVersion="11">
+   <version ClassVersion="11" checksum="2195704462"/>
    <version ClassVersion="10" checksum="682759659"/>
   </class>
   <class name="std::vector<PCaloHit>"/>

diff --git a/SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h b/SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h
@@ -24,6 +24,7 @@ class CaloSlaveSD {
   virtual void Initialize();
   std::string name() const { return name_; }
   virtual bool processHits(uint32_t, double, double, double, int, uint16_t depth = 0);
+  virtual bool processHits(uint32_t, double, double, double, int, int, uint16_t depth = 0);
   virtual bool format();
   Collection &hits() { return hits_; }
   std::string type() { return "calo"; }

diff --git a/SimDataFormats/SimHitMaker/src/CaloSlaveSD.cc b/SimDataFormats/SimHitMaker/src/CaloSlaveSD.cc
@@ -33,6 +33,13 @@ bool CaloSlaveSD::processHits(uint32_t unitID, double eDepEM, double eDepHad, do
   return true;
 }
 
+bool CaloSlaveSD::processHits(uint32_t unitID, double eDepEM, double eDepHad, double tSlice, int tkID, int fineTkId, uint16_t depth) {
+  PCaloHit aCal = PCaloHit(unitID, eDepEM, eDepHad, tSlice, tkID, fineTkId, depth);
+  LogDebug("HitBuildInfo") << " Sent Hit " << aCal << " to ROU " << name_;
+  hits_.push_back(aCal);
+  return true;
+}
+
 void CaloSlaveSD::Clean() {
   LogDebug("HitBuildIndo") << "CaloSlaveSD " << name_ << " cleaning the collection";
   Collection().swap(hits_);

diff --git a/SimDataFormats/Track/interface/SimTrack.h b/SimDataFormats/Track/interface/SimTrack.h
@@ -2,6 +2,9 @@
 #define SimTrack_H
 
 #include "SimDataFormats/Track/interface/CoreSimTrack.h"
+#include "DataFormats/Math/interface/Vector3D.h"
+#include "DataFormats/Math/interface/LorentzVector.h"
+#include "FWCore/Utilities/interface/Exception.h"
 
 class SimTrack : public CoreSimTrack {
 public:
@@ -44,12 +47,44 @@ class SimTrack : public CoreSimTrack {
 
   inline void setVertexIndex(const int v) { ivert = v; }
 
+  // Boundary crossing variables
+  void setCrossedBoundaryPosMom(int id, const math::XYZVectorD position, const math::XYZTLorentzVectorD momentum){
+    crossedBoundary_ = true;
+    idAtBoundary_ = id;
+    positionAtBoundary_ = position;
+    momentumAtBoundary_ = momentum;
+    }
+  bool crossedBoundary() const { return crossedBoundary_; }
+  math::XYZVectorD getPositionAtBoundary() const {
+    assertCrossedBoundary();
+    return positionAtBoundary_;
+    }
+  math::XYZTLorentzVectorD getMomentumAtBoundary() const {
+    assertCrossedBoundary();
+    return momentumAtBoundary_;
+    }
+  int getIDAtBoundary() const {
+    assertCrossedBoundary();
+    return idAtBoundary_;
+    }
+
 private:
   int ivert;
   int igenpart;
 
   math::XYZVectorD tkposition;
   math::XYZTLorentzVectorD tkmomentum;
+
+  bool crossedBoundary_;
+  int idAtBoundary_;
+  math::XYZVectorD positionAtBoundary_;
+  math::XYZTLorentzVectorD momentumAtBoundary_;
+  void assertCrossedBoundary() const {
+    if (!crossedBoundary_){
+      throw cms::Exception("Unknown", "SimTrack")
+        << "Assert crossed boundary failed for track " << trackId();
+      }
+    }
 };
 
 #include <iosfwd>

diff --git a/SimDataFormats/Track/src/SimTrack.cc b/SimDataFormats/Track/src/SimTrack.cc
@@ -1,21 +1,21 @@
 #include "SimDataFormats/Track/interface/SimTrack.h"
 
-SimTrack::SimTrack() : ivert(-1), igenpart(-1) {}
+SimTrack::SimTrack() : ivert(-1), igenpart(-1), crossedBoundary_(false) {}
 
-SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p) : Core(ipart, p), ivert(-1), igenpart(-1) {}
+SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p) : Core(ipart, p), ivert(-1), igenpart(-1), crossedBoundary_(false) {}
 
 SimTrack::SimTrack(int ipart, const math::XYZTLorentzVectorD& p, int iv, int ig)
-    : Core(ipart, p), ivert(iv), igenpart(ig) {}
+    : Core(ipart, p), ivert(iv), igenpart(ig), crossedBoundary_(false) {}
 
 SimTrack::SimTrack(int ipart,
                    const math::XYZTLorentzVectorD& p,
                    int iv,
                    int ig,
                    const math::XYZVectorD& tkp,
                    const math::XYZTLorentzVectorD& tkm)
-    : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm) {}
+    : Core(ipart, p), ivert(iv), igenpart(ig), tkposition(tkp), tkmomentum(tkm), crossedBoundary_(false) {}
 
-SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig) {}
+SimTrack::SimTrack(const CoreSimTrack& t, int iv, int ig) : Core(t), ivert(iv), igenpart(ig), crossedBoundary_(false) {}
 
 std::ostream& operator<<(std::ostream& o, const SimTrack& t) {
   return o << (SimTrack::Core)(t) << ", " << t.vertIndex() << ", " << t.genpartIndex();

diff --git a/SimDataFormats/Track/src/classes_def.xml b/SimDataFormats/Track/src/classes_def.xml
@@ -4,7 +4,8 @@
   <class name="CoreSimTrack" ClassVersion="10">
    <version ClassVersion="10" checksum="3936841839"/>
   </class>
-  <class name="SimTrack" ClassVersion="11">
+  <class name="SimTrack" ClassVersion="12">
+   <version ClassVersion="12" checksum="392307072"/>
    <version ClassVersion="11" checksum="1785575744"/>
    <version ClassVersion="10" checksum="1430205451"/>
   </class>

diff --git a/SimG4CMS/Calo/interface/CaloG4Hit.h b/SimG4CMS/Calo/interface/CaloG4Hit.h
@@ -62,6 +62,7 @@ class CaloG4Hit : public G4VHit {
   void setIncidentEnergy(double e) { theIncidentEnergy = e; }
 
   int getTrackID() const { return hitID.trackID(); }
+
   uint32_t getUnitID() const { return hitID.unitID(); }
   double getTimeSlice() const { return hitID.timeSlice(); }
   int getTimeSliceID() const { return hitID.timeSliceID(); }
@@ -97,6 +98,8 @@ class CaloG4HitLess {
       return (a->getUnitID() < b->getUnitID());
     } else if (a->getDepth() != b->getDepth()) {
       return (a->getDepth() < b->getDepth());
+    } else if (a->getID().getFineTrackID() != b->getID().getFineTrackID()) {
+      return (a->getID().getFineTrackID() < b->getID().getFineTrackID());
     } else {
       return (a->getTimeSliceID() < b->getTimeSliceID());
     }
@@ -107,7 +110,7 @@ class CaloG4HitEqual {
 public:
   bool operator()(const CaloG4Hit* a, const CaloG4Hit* b) {
     return (a->getTrackID() == b->getTrackID() && a->getUnitID() == b->getUnitID() && a->getDepth() == b->getDepth() &&
-            a->getTimeSliceID() == b->getTimeSliceID());
+            a->getTimeSliceID() == b->getTimeSliceID() && a->getID().getFineTrackID() == b->getID().getFineTrackID());
   }
 };
 

diff --git a/SimG4CMS/Calo/interface/CaloHitID.h b/SimG4CMS/Calo/interface/CaloHitID.h
@@ -25,6 +25,13 @@ class CaloHitID {
   void setID(uint32_t unitID, double timeSlice, int trackID, uint16_t depth = 0);
   void reset();
 
+  bool hasFineTrackID() const { return hasFineTrackID_; }
+  void setFineTrackID(int fineTrackID){
+    theFineTrackID = fineTrackID;
+    hasFineTrackID_ = true;
+    }
+  int getFineTrackID() const { return hasFineTrackID_ ? theFineTrackID : theTrackID; }
+
   bool operator==(const CaloHitID&) const;
   bool operator<(const CaloHitID&) const;
   bool operator>(const CaloHitID&) const;
@@ -37,6 +44,8 @@ class CaloHitID {
   uint16_t theDepth;
   float timeSliceUnit;
   bool ignoreTrackID;
+  bool hasFineTrackID_;
+  int theFineTrackID;
 };
 
 std::ostream& operator<<(std::ostream&, const CaloHitID&);

diff --git a/SimG4CMS/Calo/interface/CaloSD.h b/SimG4CMS/Calo/interface/CaloSD.h
@@ -162,7 +162,8 @@ class CaloSD : public SensitiveCaloDetector,
   float timeSlice;
   double eminHitD;
   double correctT;
-  bool useFineCaloID_;
+  bool doFineCalo_;
+  double eMinFine_;
 
   std::map<CaloHitID, CaloG4Hit*> hitMap;
   std::map<int, TrackWithHistory*> tkMap;

diff --git a/SimG4CMS/Calo/interface/CaloTrkProcessing.h b/SimG4CMS/Calo/interface/CaloTrkProcessing.h
@@ -57,7 +57,7 @@ class CaloTrkProcessing : public SensitiveCaloDetector,
   G4LogicalVolume* detLV(const G4VTouchable*, int) const;
   void detectorLevel(const G4VTouchable*, int&, int*, G4String*) const;
 
-  bool testBeam_, putHistory_, doFineCalo_, storeHGCBoundaryCross_;
+  bool testBeam_, putHistory_, doFineCalo_, storeAllTracksCalo_;
   double eMin_, eMinFine_, eMinFinePhoton_;
   int lastTrackID_;
   std::vector<Detector> detectors_, fineDetectors_;

diff --git a/SimG4CMS/Calo/src/CaloG4Hit.cc b/SimG4CMS/Calo/src/CaloG4Hit.cc
@@ -38,7 +38,6 @@ const CaloG4Hit& CaloG4Hit::operator=(const CaloG4Hit& right) {
   hadr = right.hadr;
   theIncidentEnergy = right.theIncidentEnergy;
   hitID = right.hitID;
-
   return *this;
 }