diff --git a/DataFormats/HGCalReco/interface/TICLLayerTile.h b/DataFormats/HGCalReco/interface/TICLLayerTile.h index 01f0b2e3ca7b0..50fcc679d9c4f 100644 --- 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 index fc30404ff2366..d5b82a4d2211e 100644 --- a/RecoHGCal/GraphReco/BuildFile.xml +++ b/RecoHGCal/GraphReco/BuildFile.xml @@ -13,4 +13,4 @@ - \ No newline at end of file + diff --git a/RecoHGCal/GraphReco/plugins/BuildFile.xml b/RecoHGCal/GraphReco/plugins/BuildFile.xml index 932fce8cc0b24..c29cc0fe18a70 100644 --- a/RecoHGCal/GraphReco/plugins/BuildFile.xml +++ b/RecoHGCal/GraphReco/plugins/BuildFile.xml @@ -2,9 +2,12 @@ + - + + + diff --git a/RecoHGCal/GraphReco/plugins/HGCClusterTrackLinker.cc b/RecoHGCal/GraphReco/plugins/HGCClusterTrackLinker.cc new file mode 100644 index 0000000000000..525ca57965cd9 --- /dev/null +++ b/RecoHGCal/GraphReco/plugins/HGCClusterTrackLinker.cc @@ -0,0 +1,102 @@ +#include +#include + +#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> tracksToken_; + edm::EDGetTokenT pfCandsToken_; + +}; + +HGCClusterTrackLinker::HGCClusterTrackLinker(const edm::ParameterSet& config) : + tracksToken_(consumes>(config.getParameter("tracks"))), + pfCandsToken_(consumes(config.getParameter("pfCands"))) { + produces(); +} + +void HGCClusterTrackLinker::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + edm::Handle pfCands; + iEvent.getByToken(pfCandsToken_, pfCands); + + edm::Handle> tracks; + iEvent.getByToken(tracksToken_, tracks); + + ticl::TICLLayerTile candTile; + ticl::TICLLayerTile tracksTile; + + auto out = std::make_unique(); + + for (size_t i = 0; i < tracks->size(); i++) { + edm::RefToBase 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(candTile.nBins()); bin++) { + const std::vector& candIndices = candTile[bin]; + if (!candIndices.size()) + continue; + + std::vector trackIndices = tracksTile[bin]; + for (auto& tidx : trackIndices) { + edm::RefToBase 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 ptdiffs; + ptdiffs.reserve(trackIndices.size()); + for (size_t entry : trackIndices) { + edm::RefToBase 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 matchingTrack(tracks, trackIndices.at(tidx)); + trackIndices.erase(trackIndices.begin()+tidx); + + cand.setTrackRef(matchingTrack.castTo()); + } + } + } + + iEvent.put(std::move(out)); +} + +DEFINE_FWK_MODULE(HGCClusterTrackLinker); diff --git a/SimCalorimetry/EcalElectronicsEmulation/src/EcalSimpleProducer.cc b/SimCalorimetry/EcalElectronicsEmulation/src/EcalSimpleProducer.cc index f619ff7563d5a..36f36c56d875b 100644 --- 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 index e1184ce78aef6..78f51c75e4f8a 100644 --- 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 index bf2c76c9055c0..7f512794645e7 100644 --- 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 index 6e7cbbd897498..b63aa2b1e7868 100644 --- a/SimDataFormats/CaloHit/src/classes_def.xml +++ b/SimDataFormats/CaloHit/src/classes_def.xml @@ -1,5 +1,6 @@ - + + diff --git a/SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h b/SimDataFormats/SimHitMaker/interface/CaloSlaveSD.h index e89a3ceee7742..06ceee2727e70 100644 --- 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 index e8bb371bbc75f..d3a97c10bd694 100644 --- 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 index 3cedee725b3ee..ec5adc8d86e1a 100644 --- 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 diff --git a/SimDataFormats/Track/src/SimTrack.cc b/SimDataFormats/Track/src/SimTrack.cc index 0fe3d889c10cd..2fac53edfd7de 100644 --- a/SimDataFormats/Track/src/SimTrack.cc +++ b/SimDataFormats/Track/src/SimTrack.cc @@ -1,11 +1,11 @@ #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, @@ -13,9 +13,9 @@ SimTrack::SimTrack(int ipart, 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 index a2a52b4a404be..8c668e7b57815 100644 --- a/SimDataFormats/Track/src/classes_def.xml +++ b/SimDataFormats/Track/src/classes_def.xml @@ -4,7 +4,8 @@ - + + diff --git a/SimG4CMS/Calo/interface/CaloG4Hit.h b/SimG4CMS/Calo/interface/CaloG4Hit.h index 55e22de003429..bca9f65083b3e 100644 --- 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 index 5f644f93d6421..4abffa1fb3b7b 100644 --- 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 index 7b02d993f2904..c6d016ef77484 100644 --- 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 hitMap; std::map tkMap; diff --git a/SimG4CMS/Calo/interface/CaloTrkProcessing.h b/SimG4CMS/Calo/interface/CaloTrkProcessing.h index c7d182c99657f..74ea8f06dc32e 100644 --- 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 detectors_, fineDetectors_; diff --git a/SimG4CMS/Calo/src/CaloG4Hit.cc b/SimG4CMS/Calo/src/CaloG4Hit.cc index 370921a223060..92eacba34df98 100644 --- 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; } diff --git a/SimG4CMS/Calo/src/CaloHitID.cc b/SimG4CMS/Calo/src/CaloHitID.cc index fa64c6f07e136..b05f8dc91bee7 100644 --- a/SimG4CMS/Calo/src/CaloHitID.cc +++ b/SimG4CMS/Calo/src/CaloHitID.cc @@ -21,6 +21,8 @@ CaloHitID::CaloHitID(const CaloHitID& id) { theDepth = id.theDepth; timeSliceUnit = id.timeSliceUnit; ignoreTrackID = id.ignoreTrackID; + hasFineTrackID_ = id.hasFineTrackID_; + theFineTrackID = id.theFineTrackID; } const CaloHitID& CaloHitID::operator=(const CaloHitID& id) { @@ -31,7 +33,8 @@ const CaloHitID& CaloHitID::operator=(const CaloHitID& id) { theDepth = id.theDepth; timeSliceUnit = id.timeSliceUnit; ignoreTrackID = id.ignoreTrackID; - + hasFineTrackID_ = id.hasFineTrackID_; + theFineTrackID = id.theFineTrackID; return *this; } @@ -51,11 +54,13 @@ void CaloHitID::reset() { theTrackID = -2; theTimeSliceID = (int)(theTimeSlice / timeSliceUnit); theDepth = 0; + theFineTrackID = -3; + hasFineTrackID_ = false; } bool CaloHitID::operator==(const CaloHitID& id) const { return ((theUnitID == id.unitID()) && (theTrackID == id.trackID() || ignoreTrackID) && - (theTimeSliceID == id.timeSliceID()) && (theDepth == id.depth())) + (theTimeSliceID == id.timeSliceID()) && (theDepth == id.depth()) && (getFineTrackID() == id.getFineTrackID())) ? true : false; } @@ -63,6 +68,8 @@ bool CaloHitID::operator==(const CaloHitID& id) const { bool CaloHitID::operator<(const CaloHitID& id) const { if (theTrackID != id.trackID()) { return (theTrackID > id.trackID()); + } else if (getFineTrackID() != id.getFineTrackID()) { + return (getFineTrackID() > id.getFineTrackID()); } else if (theUnitID != id.unitID()) { return (theUnitID > id.unitID()); } else if (theDepth != id.depth()) { @@ -75,6 +82,8 @@ bool CaloHitID::operator<(const CaloHitID& id) const { bool CaloHitID::operator>(const CaloHitID& id) const { if (theTrackID != id.trackID()) { return (theTrackID < id.trackID()); + } else if (getFineTrackID() != id.getFineTrackID()) { + return (getFineTrackID() < id.getFineTrackID()); } else if (theUnitID != id.unitID()) { return (theUnitID < id.unitID()); } else if (theDepth != id.depth()) { @@ -87,5 +96,6 @@ bool CaloHitID::operator>(const CaloHitID& id) const { std::ostream& operator<<(std::ostream& os, const CaloHitID& id) { os << "UnitID 0x" << std::hex << id.unitID() << std::dec << " Depth " << std::setw(6) << id.depth() << " Time " << std::setw(6) << id.timeSlice() << " TrackID " << std::setw(8) << id.trackID(); + if (id.hasFineTrackID()) os << " fineTrackID " << id.getFineTrackID(); return os; } diff --git a/SimG4CMS/Calo/src/CaloSD.cc b/SimG4CMS/Calo/src/CaloSD.cc index dc3b6a2a138cf..325c443740764 100644 --- a/SimG4CMS/Calo/src/CaloSD.cc +++ b/SimG4CMS/Calo/src/CaloSD.cc @@ -8,6 +8,7 @@ #include "SimG4Core/Notification/interface/TrackInformation.h" #include "SimG4Core/Notification/interface/G4TrackToParticleID.h" #include "SimG4Core/Notification/interface/SimTrackManager.h" +#include "FWCore/Utilities/interface/Exception.h" #include "G4EventManager.hh" #include "G4SDManager.hh" @@ -21,8 +22,9 @@ #include "G4PhysicalConstants.hh" #include +#include -//#define EDM_ML_DEBUG +#define EDM_ML_DEBUG CaloSD::CaloSD(const std::string& name, const edm::EventSetup& es, @@ -60,7 +62,8 @@ CaloSD::CaloSD(const std::string& name, corrTOFBeam = m_CaloSD.getParameter("CorrectTOFBeam"); double beamZ = m_CaloSD.getParameter("BeamPosition") * CLHEP::cm; correctT = beamZ / CLHEP::c_light / CLHEP::nanosecond; - useFineCaloID_ = m_CaloSD.getParameter("UseFineCaloID"); + doFineCalo_ = m_CaloSD.getParameter("DoFineCalo"); + eMinFine_ = m_CaloSD.getParameter("EMinFine"); SetVerboseLevel(verbn); meanResponse.reset(nullptr); @@ -99,7 +102,7 @@ CaloSD::CaloSD(const std::string& name, << " ns and if energy is above " << eminHit / CLHEP::MeV << " MeV (for depth 0) or " << eminHitD / CLHEP::MeV << " MeV (for nonzero depths);\n Time Slice Unit " << timeSlice << "\nIgnore TrackID Flag " << ignoreTrackID << " UseFineCaloID flag " - << useFineCaloID_; + << doFineCalo_; } CaloSD::~CaloSD() {} @@ -139,6 +142,7 @@ G4bool CaloSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { int primaryID = getTrackID(theTrack); if (unitID > 0) { currentID.setID(unitID, time, primaryID, depth); + if (doFineCalo_) currentID.setFineTrackID(theTrack->GetTrackID()); } else { if (aStep->GetTotalEnergyDeposit() > 0.0) { const G4TouchableHistory* touch = static_cast(theTrack->GetTouchable()); @@ -209,6 +213,7 @@ bool CaloSD::ProcessHits(G4GFlashSpot* aSpot, G4TouchableHistory*) { int primaryID = getTrackID(track); uint16_t depth = getDepth(&fFakeStep); currentID.setID(unitID, time, primaryID, depth); + if (doFineCalo_) currentID.setFineTrackID(track->GetTrackID()); #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD:: GetSpotInfo for Unit 0x" << std::hex << currentID.unitID() << std::dec << " Edeposit = " << edepositEM << " " << edepositHAD; @@ -373,36 +378,163 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) { updateHit(aHit); storeHit(aHit); - double etrack = 0; - if (currentID.trackID() == primIDSaved) { // The track is saved; nothing to be done - } else if (currentID.trackID() == theTrack->GetTrackID()) { - etrack = theTrack->GetKineticEnergy(); + + TrackInformation* trkInfo = cmsTrackInformation(theTrack); + #ifdef EDM_ML_DEBUG - edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " etrack " << etrack - << " eCut " << energyCut << " force: " << forceSave - << " save: " << (etrack >= energyCut || forceSave); -#endif - if (etrack >= energyCut || forceSave) { - TrackInformation* trkInfo = cmsTrackInformation(theTrack); + edm::LogInfo("DoFineCalo") + << "Creating new hit " << aHit->getUnitID() + << "; currentID.trackID=" << currentID.trackID() + << " currentID.getFineTrackID=" << currentID.getFineTrackID() + << " theTrack=" << theTrack->GetTrackID() + << " (parentTrackId=" << theTrack->GetParentID() + << " getIDfineCalo=" << trkInfo->getIDfineCalo() + << " getIDonCaloSurface=" << trkInfo->getIDonCaloSurface() + << ")" + << " primIDSaved=" << primIDSaved + ; +#endif + + if (doFineCalo_){ + // theTrack itself won't be available in m_trackManager.m_trksForThisEvent, + // so check the track itself first without getting TrackWithHistory from the track manager + // if (theTrack->GetKineticEnergy() > eMinFine_){ + if (trkInfo->crossedBoundary()){ +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "theTrack " << theTrack->GetTrackID() + << " crossed the boundary itself; recording it for hit " << aHit->getUnitID(); + // << " itself passes eMinFine_ (" + // << theTrack->GetKineticEnergy() << ">" << eMinFine_ << ")"; +#endif + currentID.setFineTrackID(theTrack->GetTrackID()); + aHit->setID(currentID); // Actually overwrite the ID for the hit trkInfo->storeTrack(true); - trkInfo->putInHistory(); + } + // theTrack itself does not pass thresholds / does not cross boundary; go through its history to find a track that does + else{ + // double recordTrackEnergy; + bool foundViableTrack = false; + TrackWithHistory* recordTrackWithHistory; + // Keep track of decay chain of this track for debugging purposes + std::vector decayChain; + decayChain.push_back(theTrack->GetTrackID()); + // Lambda to convert the vector to a nice printable string for debugging + auto decayChainToStr = [](std::vector decayChain) { + // std::string str = ""; + std::stringstream ss; + for(long unsigned int i=0; i < decayChain.size(); i++){ + if (i>0) ss << " <- "; + ss << decayChain[i]; + } + return ss.str(); + }; + // Find the first parent of this track that passes the energy threshold + // Start from first parent + unsigned int recordTrackID = theTrack->GetParentID(); +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Trying to find the first parent of hit " << aHit->getUnitID() + // << " that passes energy > " << eMinFine_ << " cut" + << " that crosses the boundary" + << "; starting with first parent track " << recordTrackID; +#endif + while(true){ + decayChain.push_back(recordTrackID); + recordTrackWithHistory = m_trackManager->getTrackByID(recordTrackID); + // recordTrackEnergy = sqrt(recordTrackWithHistory->momentum().Mag2()); + // if (recordTrackEnergy > eMinFine_){ + if (recordTrackWithHistory->crossedBoundary() && recordTrackWithHistory->getIDAtBoundary() == (int)recordTrackID){ + foundViableTrack = true; +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Recording track " << recordTrackID + // << " with E=" << recordTrackEnergy + << " as source of hit " << aHit->getUnitID() + << "; crossed boundary at pos=(" + << recordTrackWithHistory->getPositionAtBoundary().x() << "," + << recordTrackWithHistory->getPositionAtBoundary().y() << "," + << recordTrackWithHistory->getPositionAtBoundary().z() << ")" + << " mom=(" + << recordTrackWithHistory->getMomentumAtBoundary().x() << "," + << recordTrackWithHistory->getMomentumAtBoundary().y() << "," + << recordTrackWithHistory->getMomentumAtBoundary().z() << "," + << recordTrackWithHistory->getMomentumAtBoundary().e() << ")" + << "id@boundary=" << recordTrackWithHistory->getIDAtBoundary() + << "; decayChain: " << decayChainToStr(decayChain); +#endif + break; + } + else{ + // Go to the next parent +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Track " << recordTrackID + // << " does not pass energy threshold, E=" << recordTrackEnergy << " < " << eMinFine_; + << " did not cross the boundary"; +#endif + recordTrackID = recordTrackWithHistory->parentID(); + if (recordTrackID <= 0){ + throw cms::Exception("Unknown", "CaloSD") + << "Hit " << aHit->getUnitID() + << " does not have any parent track that passes the energy threshold!" + << " decayChain so far: " << decayChainToStr(decayChain) + ; + } + } + } + // Store the first viable track, or raise if it couldn't be found + if (!foundViableTrack){ + throw cms::Exception("Unknown", "CaloSD") + << "Something went wrong determining the track for hit " << aHit->getUnitID() + << "; no viable track was found." + << " decayChain so far: " << decayChainToStr(decayChain); + } + recordTrackWithHistory->save(); + currentID.setFineTrackID(recordTrackID); + aHit->setID(currentID); // Actually overwrite the ID for the hit +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "currentID.getFineTrackID()=" << currentID.getFineTrackID() + << " recordTrackWithHistory->trackID()=" << recordTrackWithHistory->trackID() + << " recordTrackWithHistory->saved()=" << recordTrackWithHistory->saved(); +#endif + } } - } else { - TrackWithHistory* trkh = tkMap[currentID.trackID()]; + // Non-fine history bookkeeping + else{ + double etrack = 0; + if (currentID.trackID() == primIDSaved) { // The track is saved; nothing to be done + } else if (currentID.trackID() == theTrack->GetTrackID()) { + etrack = theTrack->GetKineticEnergy(); +#ifdef EDM_ML_DEBUG + edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " etrack " << etrack + << " eCut " << energyCut << " force: " << forceSave + << " save: " << (etrack >= energyCut || forceSave); +#endif + if (etrack >= energyCut || forceSave) { + TrackInformation* trkInfo = cmsTrackInformation(theTrack); + trkInfo->storeTrack(true); + trkInfo->putInHistory(); + } + } else { + TrackWithHistory* trkh = tkMap[currentID.trackID()]; #ifdef EDM_ML_DEBUG - edm::LogVerbatim("CaloSim") << "CaloSD : TrackwithHistory pointer for " << currentID.trackID() << " is " << trkh; + edm::LogVerbatim("CaloSim") << "CaloSD : TrackwithHistory pointer for " << currentID.trackID() << " is " << trkh; #endif - if (trkh != nullptr) { - etrack = sqrt(trkh->momentum().Mag2()); - if (etrack >= energyCut) { - trkh->save(); + if (trkh != nullptr) { + etrack = sqrt(trkh->momentum().Mag2()); + if (etrack >= energyCut) { + trkh->save(); #ifdef EDM_ML_DEBUG - edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " with Hit"; + edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " with Hit"; #endif + } } } + primIDSaved = currentID.trackID(); } - primIDSaved = currentID.trackID(); + if (useMap) ++totalHits; return aHit; @@ -572,10 +704,9 @@ void CaloSD::endEvent() {} int CaloSD::getTrackID(const G4Track* aTrack) { int primaryID = 0; - forceSave = false; TrackInformation* trkInfo = cmsTrackInformation(aTrack); if (trkInfo) { - primaryID = (useFineCaloID_) ? trkInfo->getIDfineCalo() : trkInfo->getIDonCaloSurface(); + primaryID = (doFineCalo_) ? trkInfo->getIDfineCalo() : trkInfo->getIDonCaloSurface(); #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "Track ID: " << trkInfo->getIDfineCalo() << ":" << trkInfo->getIDonCaloSurface() << ":" << aTrack->GetTrackID() << ":" << primaryID; @@ -592,7 +723,7 @@ int CaloSD::getTrackID(const G4Track* aTrack) { int CaloSD::setTrackID(const G4Step* aStep) { auto const theTrack = aStep->GetTrack(); TrackInformation* trkInfo = cmsTrackInformation(theTrack); - int primaryID = (useFineCaloID_) ? trkInfo->getIDfineCalo() : trkInfo->getIDonCaloSurface(); + int primaryID = (doFineCalo_) ? trkInfo->getIDfineCalo() : trkInfo->getIDonCaloSurface(); if (primaryID <= 0) { primaryID = theTrack->GetTrackID(); } @@ -647,31 +778,79 @@ void CaloSD::storeHit(CaloG4Hit* hit) { bool CaloSD::saveHit(CaloG4Hit* aHit) { int tkID; + int fineTrackID; bool ok = true; - if (m_trackManager) { - tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID()); - if (tkID == 0) { - if (m_trackManager->trackExists(aHit->getTrackID())) - tkID = (aHit->getTrackID()); - else { + + double time = aHit->getTimeSlice(); + if (corrTOFBeam) + time += correctT; + + // Do track bookkeeping a little differently for fine tracking + if (doFineCalo_){ + tkID = aHit->getTrackID(); + fineTrackID = aHit->getID().getFineTrackID(); + // Check if the track is actually in the trackManager + if (m_trackManager){ + if (!m_trackManager->trackExists(tkID)){ ok = false; + throw cms::Exception("Unknown", "CaloSD") + << "aHit " << aHit->getUnitID() + << " has primary trackID " << tkID + << ", which is NOT IN THE TRACK MANAGER" + ; + } + if (!m_trackManager->trackExists(fineTrackID)){ + ok = false; + throw cms::Exception("Unknown", "CaloSD") + << "aHit " << aHit->getUnitID() + << " has fineTrackID " << fineTrackID + << ", which is NOT IN THE TRACK MANAGER" + ; + } } + else { + ok = false; + throw cms::Exception("Unknown", "CaloSD") << "m_trackManager not set, saveHit ok=false!"; + } +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Saving hit " << aHit->getUnitID() + << " with trackID=" << tkID << " fineTrackID=" << fineTrackID; +#endif + slave.get()->processHits( + aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, fineTrackID, aHit->getDepth()); } - } else { - tkID = aHit->getTrackID(); - ok = false; + // Regular, not-fine way: + else { + if (m_trackManager) { + tkID = m_trackManager->giveMotherNeeded(aHit->getTrackID()); + if (tkID == 0) { + if (m_trackManager->trackExists(aHit->getTrackID())) + tkID = (aHit->getTrackID()); + else { + ok = false; + } + } + } else { + tkID = aHit->getTrackID(); + ok = false; + } +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Saving hit " << aHit->getUnitID() + << " with trackID=" << tkID << " (no fineTrackID)"; +#endif + slave.get()->processHits( + aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, aHit->getDepth()); } + #ifdef EDM_ML_DEBUG if (!ok) edm::LogWarning("CaloSim") << "CaloSD:Cannot find track ID for " << aHit->getTrackID(); edm::LogVerbatim("CaloSim") << "CalosD: Track ID " << aHit->getTrackID() << " changed to " << tkID << " by SimTrackManager Status " << ok; #endif - double time = aHit->getTimeSlice(); - if (corrTOFBeam) - time += correctT; - slave.get()->processHits( - aHit->getUnitID(), aHit->getEM() / CLHEP::GeV, aHit->getHadr() / CLHEP::GeV, time, tkID, aHit->getDepth()); + #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD: Store Hit at " << std::hex << aHit->getUnitID() << std::dec << " " << aHit->getDepth() << " due to " << tkID << " in time " << time << " of energy " @@ -722,8 +901,7 @@ void CaloSD::cleanHitCollection() { #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: sort hits in buffer starting from " << "element = " << cleanIndex; - for (unsigned int i = 0; i < hitvec.size(); ++i) - edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i]; + for (unsigned int i = 0; i < hitvec.size(); ++i) edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i]; #endif CaloG4HitEqual equal; for (unsigned int i = cleanIndex; i < hitvec.size(); ++i) { @@ -741,8 +919,6 @@ void CaloSD::cleanHitCollection() { } #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD: cleanHitCollection merge the hits in buffer "; - for (unsigned int i = 0; i < hitvec.size(); ++i) - edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i]; #endif //move all nullptr to end of list and then remove them hitvec.erase( @@ -751,8 +927,8 @@ void CaloSD::cleanHitCollection() { #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD::cleanHitCollection: remove the merged hits in buffer," << " new size = " << hitvec.size(); - for (unsigned int i = 0; i < hitvec.size(); ++i) - edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i]; + // for (unsigned int i = 0; i < hitvec.size(); ++i) + // edm::LogVerbatim("CaloSim") << i << " " << *hitvec[i]; #endif hitvec.swap(*theCollection); totalHits = theHC->entries(); diff --git a/SimG4CMS/Calo/src/CaloTrkProcessing.cc b/SimG4CMS/Calo/src/CaloTrkProcessing.cc index a327b4cfde1a8..6132b067a814f 100644 --- a/SimG4CMS/Calo/src/CaloTrkProcessing.cc +++ b/SimG4CMS/Calo/src/CaloTrkProcessing.cc @@ -20,7 +20,7 @@ #include "G4SystemOfUnits.hh" -//#define EDM_ML_DEBUG +// #define EDM_ML_DEBUG CaloTrkProcessing::CaloTrkProcessing(const std::string& name, const edm::EventSetup& es, @@ -34,7 +34,7 @@ CaloTrkProcessing::CaloTrkProcessing(const std::string& name, eMin_ = m_p.getParameter("EminTrack") * CLHEP::MeV; putHistory_ = m_p.getParameter("PutHistory"); doFineCalo_ = m_p.getParameter("DoFineCalo"); - storeHGCBoundaryCross_ = m_p.getParameter("StoreHGCBoundaryCross"); + storeAllTracksCalo_ = m_p.getParameter("StoreAllTracksCalo"); eMinFine_ = m_p.getParameter("EminFineTrack") * CLHEP::MeV; eMinFinePhoton_ = m_p.getParameter("EminFinePhoton") * CLHEP::MeV; @@ -187,6 +187,68 @@ void CaloTrkProcessing::update(const G4Step* aStep) { throw cms::Exception("Unknown", "CaloTrkProcessing") << "cannot get trkInfo for Track " << id << "\n"; } + if (doFineCalo_ || storeAllTracksCalo_) { + // Boundary-crossing logic + // const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable(); + int prestepLV = isItCalo(aStep->GetPreStepPoint()->GetTouchable(), fineDetectors_); + int poststepLV = isItCalo(aStep->GetPostStepPoint()->GetTouchable(), fineDetectors_); + if ( + prestepLV < 0 && poststepLV >= 0 + // Require abs(pre z position) < abs(current z position) to prevent back scattering tracks from being counted + && std::abs(theTrack->GetStep()->GetPreStepPoint()->GetPosition().z()) < std::abs(theTrack->GetPosition().z()) + ) { + edm::LogInfo("DoFineCalo") + << "Crossed boundary:" + << " Track " << id + << " prestepLV=" << prestepLV + << " poststepLV=" << poststepLV + << " GetKineticEnergy=" << theTrack->GetKineticEnergy() + << " GetVertexKineticEnergy=" << theTrack->GetVertexKineticEnergy() + << " prestepPosition=(" + << theTrack->GetStep()->GetPreStepPoint()->GetPosition().x() << "," + << theTrack->GetStep()->GetPreStepPoint()->GetPosition().y() << "," + << theTrack->GetStep()->GetPreStepPoint()->GetPosition().z() << ")" + << " poststepPosition=(" + << theTrack->GetStep()->GetPostStepPoint()->GetPosition().x() << "," + << theTrack->GetStep()->GetPostStepPoint()->GetPosition().y() << "," + << theTrack->GetStep()->GetPostStepPoint()->GetPosition().z() << ")" + << " position=(" + << theTrack->GetPosition().x() << "," + << theTrack->GetPosition().y() << "," + << theTrack->GetPosition().z() << ")" + << " vertex_position=(" + << theTrack->GetVertexPosition().x() << "," + << theTrack->GetVertexPosition().y() << "," + << theTrack->GetVertexPosition().z() << ")" + ; + trkInfo->setCrossedBoundary(theTrack); + } + // Decide whether to store in history + if (!trkInfo->isInHistory()){ + // For fine calo, put every single track in history + trkInfo->putInHistory(); + if (storeAllTracksCalo_) trkInfo->storeTrack(true); +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "Putting in history:" + << " Track " << id + << " vertex=(" + << theTrack->GetVertexPosition().x() << "," + << theTrack->GetVertexPosition().y() << "," + << theTrack->GetVertexPosition().z() << ")" + << " position=(" + << theTrack->GetPosition().x() << "," + << theTrack->GetPosition().y() << "," + << theTrack->GetPosition().z() << ")" + << " energy=" << theTrack->GetKineticEnergy() + << " getIDfineCalo=" << trkInfo->getIDfineCalo() + << " getIDonCaloSurface=" << trkInfo->getIDonCaloSurface() + << " parentID=" << theTrack->GetParentID() + ; +#endif + } + } + if (testBeam_) { if (trkInfo->getIDonCaloSurface() == 0) { #ifdef EDM_ML_DEBUG @@ -243,19 +305,12 @@ void CaloTrkProcessing::update(const G4Step* aStep) { } } } - if ((doFineCalo_ || storeHGCBoundaryCross_) && (!trkInfo->isInHistory())) { + if (doFineCalo_ && (!trkInfo->isInHistory())) { const G4VTouchable* pre_touch = aStep->GetPreStepPoint()->GetTouchable(); - - bool crossedHGCBoundary = false; - if (storeHGCBoundaryCross_) { - const G4VTouchable* post_touch = aStep->GetPostStepPoint()->GetTouchable(); - crossedHGCBoundary = isItCalo(pre_touch, fineDetectors_) < 0 && isItCalo(post_touch, fineDetectors_) >= 0; - } - - if (isItCalo(pre_touch, fineDetectors_) >= 0 || crossedHGCBoundary) { + if (isItCalo(pre_touch, fineDetectors_) >= 0) { int pdg = aStep->GetTrack()->GetDefinition()->GetPDGEncoding(); double cut = (pdg == 22) ? eMinFinePhoton_ : eMinFine_; - if (crossedHGCBoundary || aStep->GetTrack()->GetKineticEnergy() / CLHEP::MeV > cut) { + if (aStep->GetTrack()->GetKineticEnergy() / CLHEP::MeV > cut) { trkInfo->putInHistory(); trkInfo->setIDfineCalo(id); } diff --git a/SimG4Core/Application/interface/TrackingAction.h b/SimG4Core/Application/interface/TrackingAction.h index 059b3c07a8b16..3a823306c6a95 100644 --- a/SimG4Core/Application/interface/TrackingAction.h +++ b/SimG4Core/Application/interface/TrackingAction.h @@ -35,6 +35,7 @@ class TrackingAction : public G4UserTrackingAction { CMSSteppingVerbose* steppingVerbose_; const G4Track* g4Track_; bool checkTrack_; + bool doFineCalo_; }; #endif diff --git a/SimG4Core/Application/python/g4SimHits_cfi.py b/SimG4Core/Application/python/g4SimHits_cfi.py index 562902502a093..3085995ac17ae 100644 --- a/SimG4Core/Application/python/g4SimHits_cfi.py +++ b/SimG4Core/Application/python/g4SimHits_cfi.py @@ -254,7 +254,8 @@ ), TrackingAction = cms.PSet( DetailedTiming = cms.untracked.bool(False), - CheckTrack = cms.untracked.bool(False) + CheckTrack = cms.untracked.bool(False), + DoFineCalo = cms.untracked.bool(False) ), SteppingAction = cms.PSet( common_maximum_time, @@ -287,7 +288,8 @@ UseResponseTables = cms.vint32(0,0,0,0,0), BeamPosition = cms.double(0.0), CorrectTOFBeam = cms.bool(False), - UseFineCaloID = cms.bool(False), + DoFineCalo = cms.bool(False), + EMinFine = cms.double(10.), DetailedTiming = cms.untracked.bool(False), UseMap = cms.untracked.bool(False), Verbosity = cms.untracked.int32(0), @@ -354,6 +356,7 @@ EminTrack = cms.double(0.01), PutHistory = cms.bool(False), DoFineCalo = cms.bool(False), + StoreAllTracksCalo = cms.bool(False), EminFineTrack = cms.double(10000.0), EminFinePhoton = cms.double(5000.0), StoreHGCBoundaryCross = cms.bool(True), diff --git a/SimG4Core/Application/src/TrackingAction.cc b/SimG4Core/Application/src/TrackingAction.cc index ea899fb9ecee4..65eaa67f7fedc 100644 --- a/SimG4Core/Application/src/TrackingAction.cc +++ b/SimG4Core/Application/src/TrackingAction.cc @@ -20,7 +20,8 @@ TrackingAction::TrackingAction(EventAction* e, const edm::ParameterSet& p, CMSSt currentTrack_(nullptr), steppingVerbose_(sv), g4Track_(nullptr), - checkTrack_(p.getUntrackedParameter("CheckTrack", false)) {} + checkTrack_(p.getUntrackedParameter("CheckTrack", false)), + doFineCalo_(p.getUntrackedParameter("DoFineCalo", false)) {} TrackingAction::~TrackingAction() {} @@ -42,17 +43,49 @@ void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) { void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { if (eventAction_->trackContainer() != nullptr) { - if (extractor_(aTrack).storeTrack()) { - currentTrack_->save(); + uint32_t id = aTrack->GetTrackID(); + math::XYZVectorD pos( + aTrack->GetStep()->GetPostStepPoint()->GetPosition().x(), + aTrack->GetStep()->GetPostStepPoint()->GetPosition().y(), + aTrack->GetStep()->GetPostStepPoint()->GetPosition().z() + ); + math::XYZTLorentzVectorD mom; + std::pair p(pos, mom); + +#ifdef EDM_ML_DEBUG + edm::LogInfo("DoFineCalo") + << "PostUserTrackingAction:" + << " aTrack->GetTrackID()=" << aTrack->GetTrackID() + << " currentTrack_->saved()=" << currentTrack_->saved() + << " PostStepPosition=(" + << pos.x() << "," << pos.y() << "," << pos.z() << ")" + ; +#endif - math::XYZVectorD pos((aTrack->GetStep()->GetPostStepPoint()->GetPosition()).x(), - (aTrack->GetStep()->GetPostStepPoint()->GetPosition()).y(), - (aTrack->GetStep()->GetPostStepPoint()->GetPosition()).z()); - math::XYZTLorentzVectorD mom; + if (doFineCalo_){ + TrackInformation* trkInfo = (TrackInformation*)aTrack->GetUserInformation(); + if (trkInfo->crossedBoundary()){ + currentTrack_->save(); + currentTrack_->setCrossedBoundaryPosMom(id, trkInfo->getPositionAtBoundary(), trkInfo->getMomentumAtBoundary()); + edm::LogInfo("DoFineCalo") + << "PostUserTrackingAction:" + << " Track " << aTrack->GetTrackID() + << " crossed boundary; pos=(" + << trkInfo->getPositionAtBoundary().x() << "," + << trkInfo->getPositionAtBoundary().y() << "," + << trkInfo->getPositionAtBoundary().z() << ")" + << " mom=(" + << trkInfo->getMomentumAtBoundary().x() << "," + << trkInfo->getMomentumAtBoundary().y() << "," + << trkInfo->getMomentumAtBoundary().z() << "," + << trkInfo->getMomentumAtBoundary().e() << ")" + ; + } + } - uint32_t id = aTrack->GetTrackID(); + if (extractor_(aTrack).storeTrack() || currentTrack_->saved()) { + currentTrack_->save(); - std::pair p(pos, mom); eventAction_->addTkCaloStateInfo(id, p); #ifdef EDM_ML_DEBUG edm::LogVerbatim("SimTrackManager") @@ -65,11 +98,18 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { (extractor_(aTrack).getIDfineCalo() == aTrack->GetTrackID()) || (extractor_(aTrack).isAncestor())); if (extractor_(aTrack).isInHistory()) { + // check with end-of-track information if (checkTrack_) { currentTrack_->checkAtEnd(aTrack); } + if (doFineCalo_) + // SHOULD ONLY BE DONE FOR FINE CALO: Add the post-step position for _every_ track + // in history to the TrackManager. Tracks in history _may_ be upgraded to stored + // tracks, at which point the post-step position is needed again. + eventAction_->addTkCaloStateInfo(id, p); + eventAction_->addTrack(currentTrack_, true, withAncestor); #ifdef EDM_ML_DEBUG @@ -77,10 +117,6 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { << "TrackingAction addTrack " << currentTrack_->trackID() << " E(GeV)= " << aTrack->GetKineticEnergy() << " " << aTrack->GetDefinition()->GetParticleName() << " added= " << withAncestor << " at " << aTrack->GetPosition(); - - math::XYZVectorD pos((aTrack->GetStep()->GetPostStepPoint()->GetPosition()).x(), - (aTrack->GetStep()->GetPostStepPoint()->GetPosition()).y(), - (aTrack->GetStep()->GetPostStepPoint()->GetPosition()).z()); edm::LogVerbatim("SimTrackManager") << "TrackingAction addTrack " << currentTrack_->trackID() << " added with " << true << " and " << withAncestor << " at " << pos; #endif diff --git a/SimG4Core/Notification/interface/G4SimTrack.h b/SimG4Core/Notification/interface/G4SimTrack.h index 7de9040681a70..3d4949583ca6f 100644 --- a/SimG4Core/Notification/interface/G4SimTrack.h +++ b/SimG4Core/Notification/interface/G4SimTrack.h @@ -3,6 +3,7 @@ #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Math/interface/LorentzVector.h" +#include "FWCore/Utilities/interface/Exception.h" #include #include @@ -20,7 +21,8 @@ class G4SimTrack { parentID_(-1), parentMomentum_(math::XYZVectorD(0., 0., 0.)), tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)), - tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)) {} + tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)), + crossedBoundary_(false) {} G4SimTrack(int iid, int ipart, const math::XYZVectorD& ip, double ie, int iv, int ig, const math::XYZVectorD& ipmom) : id_(iid), @@ -31,7 +33,8 @@ class G4SimTrack { igenpart_(ig), parentMomentum_(ipmom), tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)), - tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)) {} + tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)), + crossedBoundary_(false) {} G4SimTrack(int iid, int ipart, @@ -50,7 +53,8 @@ class G4SimTrack { igenpart_(ig), parentMomentum_(ipmom), tkSurfacePosition_(tkpos), - tkSurfaceMomentum_(tkmom) {} + tkSurfaceMomentum_(tkmom), + crossedBoundary_(false) {} ~G4SimTrack() {} @@ -69,6 +73,27 @@ class G4SimTrack { // is stored, else = -1) int parentID() const { return parentID_; } + // 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 id_; int ipart_; @@ -80,6 +105,16 @@ class G4SimTrack { math::XYZVectorD parentMomentum_; math::XYZVectorD tkSurfacePosition_; math::XYZTLorentzVectorD tkSurfaceMomentum_; + bool crossedBoundary_; + int idAtBoundary_; + math::XYZVectorD positionAtBoundary_; + math::XYZTLorentzVectorD momentumAtBoundary_; + void assertCrossedBoundary() const { + if (!crossedBoundary_){ + throw cms::Exception("Unknown", "G4SimTrack") + << "Assert crossed boundary failed for track " << id_; + } + } }; #endif diff --git a/SimG4Core/Notification/interface/SimTrackManager.h b/SimG4Core/Notification/interface/SimTrackManager.h index 7df75d4b73cb7..efb10756c0442 100644 --- a/SimG4Core/Notification/interface/SimTrackManager.h +++ b/SimG4Core/Notification/interface/SimTrackManager.h @@ -25,8 +25,8 @@ // user include files #include "SimG4Core/Notification/interface/TrackWithHistory.h" #include "SimG4Core/Notification/interface/TrackContainer.h" - #include "SimDataFormats/Forward/interface/LHCTransportLinkContainer.h" +#include "FWCore/Utilities/interface/Exception.h" // forward declarations @@ -101,6 +101,29 @@ class SimTrackManager { } return flag; } + TrackWithHistory* getTrackByID(unsigned int trackID) const { + bool trackFound = false; + TrackWithHistory* track; + if (m_trksForThisEvent == nullptr){ + throw cms::Exception("Unknown", "SimTrackManager") + << "m_trksForThisEvent is a nullptr, cannot get any track!"; + } + for (unsigned int itr = 0; itr < (*m_trksForThisEvent).size(); ++itr) { + if ((*m_trksForThisEvent)[itr]->trackID() == trackID) { + track = (*m_trksForThisEvent)[itr]; + trackFound = true; + break; + } + } + if (!trackFound){ + throw cms::Exception("Unknown", "SimTrackManager") + << "Attempted to get track " << trackID + << " from SimTrackManager, but it's not in m_trksForThisEvent (" + << (*m_trksForThisEvent).size() << " tracks in m_trksForThisEvent)" + << "\n"; + } + return track; + } void setLHCTransportLink(const edm::LHCTransportLinkContainer* thisLHCTlink) { theLHCTlink = thisLHCTlink; } private: diff --git a/SimG4Core/Notification/interface/TrackInformation.h b/SimG4Core/Notification/interface/TrackInformation.h index be840de56ebf5..1ecdf1729b318 100644 --- a/SimG4Core/Notification/interface/TrackInformation.h +++ b/SimG4Core/Notification/interface/TrackInformation.h @@ -1,9 +1,12 @@ #ifndef SimG4Core_TrackInformation_H #define SimG4Core_TrackInformation_H +#include "FWCore/Utilities/interface/Exception.h" #include "G4VUserTrackInformation.hh" - #include "G4Allocator.hh" +#include "G4Track.hh" +#include "DataFormats/Math/interface/Vector3D.h" +#include "DataFormats/Math/interface/LorentzVector.h" class TrackInformation : public G4VUserTrackInformation { public: @@ -55,6 +58,34 @@ class TrackInformation : public G4VUserTrackInformation { int getIDfineCalo() const { return ((idFineCalo_ > 0) ? idFineCalo_ : idOnCaloSurface_); } void setIDfineCalo(int id) { idFineCalo_ = id; } + // Boundary crossing variables + void setCrossedBoundary(const G4Track* track){ + crossedBoundary_ = true; + // Double-check units! Any conversions necessary? Is x,y,z,E fine for XYZTLorentzVectorD? + idAtBoundary_ = track->GetTrackID(); + positionAtBoundary_ = math::XYZVectorD( + track->GetPosition().x(), + track->GetPosition().y(), + track->GetPosition().z() + ); + momentumAtBoundary_ = math::XYZTLorentzVectorD( + track->GetMomentum().x(), track->GetMomentum().y(), track->GetMomentum().z(), track->GetKineticEnergy() + ); + } + 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_; + } + // Generator information int genParticlePID() const { return genParticlePID_; } void setGenParticlePID(int id) { genParticlePID_ = id; } @@ -84,12 +115,26 @@ class TrackInformation : public G4VUserTrackInformation { int idLastVolume_; bool caloIDChecked_; int idFineCalo_; + bool crossedBoundary_; + bool idAtBoundary_; + math::XYZVectorD positionAtBoundary_; + math::XYZTLorentzVectorD momentumAtBoundary_; + int genParticlePID_, caloSurfaceParticlePID_; double genParticleP_, caloSurfaceParticleP_; bool hasCastorHit_; int castorHitPID_; + void assertCrossedBoundary() const { + if (!crossedBoundary_){ + throw cms::Exception("Unknown", "TrackInformation") + << "Assert crossed boundary failed for track " + << getIDonCaloSurface() << " (fine: " << getIDfineCalo() << ")" + ; + } + } + // Restrict construction to friends TrackInformation() : G4VUserTrackInformation(), @@ -104,6 +149,7 @@ class TrackInformation : public G4VUserTrackInformation { idLastVolume_(-1), caloIDChecked_(false), idFineCalo_(-1), + crossedBoundary_(false), genParticlePID_(-1), caloSurfaceParticlePID_(0), genParticleP_(0), diff --git a/SimG4Core/Notification/interface/TrackWithHistory.h b/SimG4Core/Notification/interface/TrackWithHistory.h index 70709cc9f3c30..381799b03a22a 100644 --- a/SimG4Core/Notification/interface/TrackWithHistory.h +++ b/SimG4Core/Notification/interface/TrackWithHistory.h @@ -2,6 +2,7 @@ #define SimG4Core_TrackWithHistory_H #include "G4Track.hh" +#include "FWCore/Utilities/interface/Exception.h" #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Math/interface/LorentzVector.h" @@ -42,6 +43,28 @@ class TrackWithHistory { void setGenParticleID(int i) { genParticleID_ = i; } bool storeTrack() const { return storeTrack_; } bool saved() const { return saved_; } + + // 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_; + } + /** Internal consistency check (optional). * Method called at PostUserTrackingAction time, to check * if the information is consistent with that provided @@ -64,7 +87,20 @@ class TrackWithHistory { double weight_; bool storeTrack_; bool saved_; + + bool crossedBoundary_; + int idAtBoundary_; + math::XYZVectorD positionAtBoundary_; + math::XYZTLorentzVectorD momentumAtBoundary_; + int extractGenID(const G4Track *gt) const; + + void assertCrossedBoundary() const { + if (!crossedBoundary_){ + throw cms::Exception("Unknown", "TrackWithHistory") + << "Assert crossed boundary failed for track " << trackID_; + } + } }; extern G4ThreadLocal G4Allocator *fpTrackWithHistoryAllocator; diff --git a/SimG4Core/Notification/src/G4SimEvent.cc b/SimG4Core/Notification/src/G4SimEvent.cc index b8106bb6a8270..58a7578de926b 100644 --- a/SimG4Core/Notification/src/G4SimEvent.cc +++ b/SimG4Core/Notification/src/G4SimEvent.cc @@ -57,6 +57,8 @@ void G4SimEvent::load(edm::SimTrackContainer& c) const { SimTrack t = SimTrack(ip, p, iv, ig, tkpos, tkmom); t.setTrackId(id); t.setEventId(EncodedEventId(0)); + if (trk->crossedBoundary()) + t.setCrossedBoundaryPosMom(trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); c.push_back(t); } std::stable_sort(c.begin(), c.end(), IdSort()); diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index d23ca80cd98d9..d6d6989eb74d2 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -136,15 +136,19 @@ void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) { if (cit != mapTkCaloStateInfo.end()) { tcinfo = cit->second; } - simEvent->add(new G4SimTrack(trkH->trackID(), - trkH->particleID(), - trkH->momentum(), - trkH->totalEnergy(), - ivertex, - ig, - pm, - tcinfo.first, - tcinfo.second)); + G4SimTrack * g4simtrack = new G4SimTrack( + trkH->trackID(), + trkH->particleID(), + trkH->momentum(), + trkH->totalEnergy(), + ivertex, + ig, + pm, + tcinfo.first, + tcinfo.second); + if (trkH->crossedBoundary()) + g4simtrack->setCrossedBoundaryPosMom(trkH->getIDAtBoundary(), trkH->getPositionAtBoundary(), trkH->getMomentumAtBoundary()); + simEvent->add(g4simtrack); } } diff --git a/SimG4Core/Notification/src/TrackWithHistory.cc b/SimG4Core/Notification/src/TrackWithHistory.cc index d1cafdf02a477..cf46fac3f8934 100644 --- a/SimG4Core/Notification/src/TrackWithHistory.cc +++ b/SimG4Core/Notification/src/TrackWithHistory.cc @@ -42,6 +42,7 @@ TrackWithHistory::TrackWithHistory(const G4Track* g4trk) creatorProcess_ = g4trk->GetCreatorProcess(); storeTrack_ = extractor(g4trk).storeTrack(); saved_ = false; + crossedBoundary_ = false; genParticleID_ = extractGenID(g4trk); // V.I. weight is computed in the same way as before // without usage of G4Track::GetWeight()