From 5b667d332cfd73ac7a086855ccf45f08bce6341e Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 24 Mar 2023 17:52:51 +0100 Subject: [PATCH 1/5] Updated design of Geant4 user actions and MC truth handling --- SimG4CMS/Calo/src/CaloSD.cc | 23 +- .../Application/interface/SteppingAction.h | 4 +- SimG4Core/Application/src/EventAction.cc | 7 +- .../Application/src/RunManagerMTWorker.cc | 3 +- SimG4Core/Application/src/SteppingAction.cc | 19 +- SimG4Core/Application/src/TrackingAction.cc | 73 ++---- SimG4Core/Notification/interface/G4SimEvent.h | 8 +- SimG4Core/Notification/interface/G4SimTrack.h | 30 +-- .../Notification/interface/GenParticleInfo.h | 1 + .../Notification/interface/NewTrackAction.h | 5 +- .../Notification/interface/SimTrackManager.h | 98 +++----- .../Notification/interface/TrackInformation.h | 21 +- .../Notification/interface/TrackWithHistory.h | 26 +-- SimG4Core/Notification/src/G4SimEvent.cc | 50 ++--- SimG4Core/Notification/src/NewTrackAction.cc | 6 - SimG4Core/Notification/src/SimTrackManager.cc | 209 +++++++++--------- .../Notification/src/TrackWithHistory.cc | 102 +++------ 17 files changed, 256 insertions(+), 429 deletions(-) diff --git a/SimG4CMS/Calo/src/CaloSD.cc b/SimG4CMS/Calo/src/CaloSD.cc index 271f1adea7a4b..ecaca60900821 100644 --- a/SimG4CMS/Calo/src/CaloSD.cc +++ b/SimG4CMS/Calo/src/CaloSD.cc @@ -544,7 +544,7 @@ unsigned int CaloSD::findBoundaryCrossingParent(const G4Track* track, bool markA TrackWithHistory* parentTrack = m_trackManager->getTrackByID(parentID, true); if (parentTrack->crossedBoundary()) { if (markAsSaveable) - parentTrack->save(); + parentTrack->setToBeSaved(); decayChain.push_back(parentID); // Record this boundary crossing parent for all traversed ancestors for (auto ancestorID : decayChain) @@ -623,7 +623,7 @@ CaloG4Hit* CaloSD::createNewHit(const G4Step* aStep, const G4Track* theTrack) { if (trkh != nullptr) { etrack = sqrt(trkh->momentum().Mag2()); if (etrack >= energyCut) { - trkh->save(); + trkh->setToBeSaved(); #ifdef EDM_ML_DEBUG edm::LogVerbatim("CaloSim") << "CaloSD: set save the track " << currentID.trackID() << " with Hit"; #endif @@ -701,19 +701,14 @@ void CaloSD::update(const EndOfTrack* trk) { if (trkI) lastTrackID = trkI->getIDonCaloSurface(); if (id == lastTrackID) { - const TrackContainer* trksForThisEvent = m_trackManager->trackContainer(); - if (trksForThisEvent != nullptr) { - int it = (int)(trksForThisEvent->size()) - 1; - if (it >= 0) { - TrackWithHistory* trkH = (*trksForThisEvent)[it]; - if (trkH->trackID() == (unsigned int)(id)) - tkMap[id] = trkH; + auto trksForThisEvent = m_trackManager->trackContainer(); + if (!trksForThisEvent->empty()) { + TrackWithHistory* trkH = trksForThisEvent->back(); + if (trkH->trackID() == (unsigned int)(id)) { + tkMap[id] = trkH; #ifdef EDM_ML_DEBUG - edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size " - << trksForThisEvent->size() << " with ID " << trkH->trackID(); - } else { - edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size " - << trksForThisEvent->size() << " with no ID"; + edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size " + << trksForThisEvent->size() << " with ID " << trkH->trackID(); #endif } } diff --git a/SimG4Core/Application/interface/SteppingAction.h b/SimG4Core/Application/interface/SteppingAction.h index 8bd8d4ccf317b..28f4a5c5e7b94 100644 --- a/SimG4Core/Application/interface/SteppingAction.h +++ b/SimG4Core/Application/interface/SteppingAction.h @@ -14,7 +14,6 @@ #include #include -class SimTrackManager; class CMSSteppingVerbose; enum TrackStatus { @@ -31,7 +30,7 @@ enum TrackStatus { class SteppingAction : public G4UserSteppingAction { public: - explicit SteppingAction(SimTrackManager*, const CMSSteppingVerbose*, const edm::ParameterSet&, bool hasW); + explicit SteppingAction(const CMSSteppingVerbose*, const edm::ParameterSet&, bool hasW); ~SteppingAction() override = default; void UserSteppingAction(const G4Step* aStep) final; @@ -47,7 +46,6 @@ class SteppingAction : public G4UserSteppingAction { bool isLowEnergy(const G4LogicalVolume*, const G4Track*) const; void PrintKilledTrack(const G4Track*, const TrackStatus&) const; - SimTrackManager* trackManager_; const G4VPhysicalVolume* tracker{nullptr}; const G4VPhysicalVolume* calo{nullptr}; const CMSSteppingVerbose* steppingVerbose; diff --git a/SimG4Core/Application/src/EventAction.cc b/SimG4Core/Application/src/EventAction.cc index 6d65fd4640e32..1f9f25cb01116 100644 --- a/SimG4Core/Application/src/EventAction.cc +++ b/SimG4Core/Application/src/EventAction.cc @@ -19,9 +19,7 @@ EventAction::EventAction(const edm::ParameterSet& p, m_SteppingVerbose(sv), m_stopFile(p.getParameter("StopFile")), m_printRandom(p.getParameter("PrintRandomSeed")), - m_debug(p.getUntrackedParameter("debug", false)) { - m_trackManager->setCollapsePrimaryVertices(p.getParameter("CollapsePrimaryVertices")); -} + m_debug(p.getUntrackedParameter("debug", false)) {} void EventAction::BeginOfEventAction(const G4Event* anEvent) { m_trackManager->reset(); @@ -63,8 +61,7 @@ void EventAction::EndOfEventAction(const G4Event* anEvent) { m_endOfEventSignal(&e); // delete transient objects - m_trackManager->deleteTracks(); - m_trackManager->cleanTkCaloStateInfoMap(); + m_trackManager->reset(); } void EventAction::abortEvent() { m_runInterface->abortEvent(); } diff --git a/SimG4Core/Application/src/RunManagerMTWorker.cc b/SimG4Core/Application/src/RunManagerMTWorker.cc index 719f7d7ce2467..d56f5fd172eec 100644 --- a/SimG4Core/Application/src/RunManagerMTWorker.cc +++ b/SimG4Core/Application/src/RunManagerMTWorker.cc @@ -415,8 +415,7 @@ void RunManagerMTWorker::initializeUserActions() { m_evtManager->SetUserAction(userTrackingAction); } - auto userSteppingAction = - new SteppingAction(m_tls->trackManager.get(), m_sVerbose.get(), m_pSteppingAction, m_hasWatchers); + auto userSteppingAction = new SteppingAction(m_sVerbose.get(), m_pSteppingAction, m_hasWatchers); Connect(userSteppingAction); if (m_UseG4EventManager) { eventManager->SetUserAction(userSteppingAction); diff --git a/SimG4Core/Application/src/SteppingAction.cc b/SimG4Core/Application/src/SteppingAction.cc index 51a63c5028eae..0f8fd7c40af8f 100644 --- a/SimG4Core/Application/src/SteppingAction.cc +++ b/SimG4Core/Application/src/SteppingAction.cc @@ -1,6 +1,6 @@ #include "SimG4Core/Application/interface/SteppingAction.h" -#include "SimG4Core/Notification/interface/SimTrackManager.h" +#include "SimG4Core/Notification/interface/TrackInformation.h" #include "SimG4Core/Notification/interface/CMSSteppingVerbose.h" #include "G4LogicalVolumeStore.hh" @@ -15,8 +15,8 @@ //#define DebugLog -SteppingAction::SteppingAction(SimTrackManager* stm, const CMSSteppingVerbose* sv, const edm::ParameterSet& p, bool hasW) - : trackManager_(stm), steppingVerbose(sv), hasWatcher(hasW) { +SteppingAction::SteppingAction(const CMSSteppingVerbose* sv, const edm::ParameterSet& p, bool hasW) + : steppingVerbose(sv), hasWatcher(hasW) { theCriticalEnergyForVacuum = (p.getParameter("CriticalEnergyForVacuum") * CLHEP::MeV); if (0.0 < theCriticalEnergyForVacuum) { killBeamPipe = true; @@ -162,17 +162,8 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) { bool isKilled = false; if (sAlive == tstat || sVeryForward == tstat) { if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == calo) { - math::XYZVectorD pos((postStep->GetPosition()).x(), (postStep->GetPosition()).y(), (postStep->GetPosition()).z()); - - math::XYZTLorentzVectorD mom((postStep->GetMomentum()).x(), - (postStep->GetMomentum()).y(), - (postStep->GetMomentum()).z(), - postStep->GetTotalEnergy()); - - // record intersection - uint32_t id = theTrack->GetTrackID(); - std::pair p(pos, mom); - trackManager_->addTkCaloStateInfo(id, p); + TrackInformation* trkinfo = static_cast(theTrack->GetUserInformation()); + if(!trkinfo->crossedBoundary()) { trkinfo->setCrossedBoundary(theTrack); } } } else { theTrack->SetTrackStatus(fStopAndKill); diff --git a/SimG4Core/Application/src/TrackingAction.cc b/SimG4Core/Application/src/TrackingAction.cc index 1effced693bd0..3894e3753d982 100644 --- a/SimG4Core/Application/src/TrackingAction.cc +++ b/SimG4Core/Application/src/TrackingAction.cc @@ -46,7 +46,7 @@ void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) { // default. The primary is the earliest ancestor, and it must be saved. if (trkInfo_->isPrimary()) { trackManager_->cleanTracksWithHistory(); - currentTrack_->save(); + currentTrack_->setToBeSaved(); } if (nullptr != steppingVerbose_) { steppingVerbose_->trackStarted(aTrack, false); @@ -74,73 +74,30 @@ void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) { } void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { + // 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. int id = aTrack->GetTrackID(); - const auto& ppos = aTrack->GetStep()->GetPostStepPoint()->GetPosition(); - math::XYZVectorD pos(ppos.x(), ppos.y(), ppos.z()); - math::XYZTLorentzVectorD mom; - std::pair p(pos, mom); - -#ifdef EDM_ML_DEBUG - edm::LogVerbatim("DoFineCalo") << "PostUserTrackingAction:" - << " aTrack->GetTrackID()=" << id - << " currentTrack_->saved()=" << currentTrack_->saved(); -#endif - - bool boundary = false; - if (doFineCalo_) { - // 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. - trackManager_->addTkCaloStateInfo(id, p); - boundary = true; - } else if (trkInfo_->storeTrack() || currentTrack_->saved() || - (saveCaloBoundaryInformation_ && trkInfo_->crossedBoundary())) { - currentTrack_->save(); - trackManager_->addTkCaloStateInfo(id, p); - boundary = true; - } - if (boundary && trkInfo_->crossedBoundary()) { + bool ok = (trkInfo_->storeTrack() || currentTrack_->saved()); + if (trkInfo_->crossedBoundary()) { currentTrack_->setCrossedBoundaryPosMom(id, trkInfo_->getPositionAtBoundary(), trkInfo_->getMomentumAtBoundary()); - currentTrack_->save(); - -#ifdef EDM_ML_DEBUG - edm::LogVerbatim("DoFineCalo") << "PostUserTrackingAction:" - << " Track " << id << " crossed boundary; pos=(" - << trkInfo_->getPositionAtBoundary().x() << "," - << trkInfo_->getPositionAtBoundary().y() << "," - << trkInfo_->getPositionAtBoundary().z() << ")" - << " mom[GeV]=(" << trkInfo_->getMomentumAtBoundary().x() << "," - << trkInfo_->getMomentumAtBoundary().y() << "," - << trkInfo_->getMomentumAtBoundary().z() << "," - << trkInfo_->getMomentumAtBoundary().e() << ")"; -#endif + ok = (ok || saveCaloBoundaryInformation_ || doFineCalo_); } + if(ok) { currentTrack_->setToBeSaved(); } bool withAncestor = (trkInfo_->getIDonCaloSurface() == id || trkInfo_->isAncestor()); + bool isInHistory = trkInfo_->isInHistory(); - if (trkInfo_->isInHistory()) { - // check with end-of-track information - if (checkTrack_) { - currentTrack_->checkAtEnd(aTrack); - } - - trackManager_->addTrack(currentTrack_, true, withAncestor); + trackManager_->addTrack(currentTrack_, aTrack, isInHistory, withAncestor); #ifdef EDM_ML_DEBUG - edm::LogVerbatim("SimTrackManager") << "TrackingAction addTrack " << id << " " - << aTrack->GetDefinition()->GetParticleName() << " added= " << withAncestor - << " at " << aTrack->GetPosition(); + edm::LogVerbatim("TrackingAction") << "TrackingAction end track=" << id + << " " << aTrack->GetDefinition()->GetParticleName() + << " ansestor= " << withAncestor + << " saved= " << currentTrack_->saved() + << " end point " << aTrack->GetPosition(); #endif - } else { - trackManager_->addTrack(currentTrack_, false, false); - delete currentTrack_; - -#ifdef EDM_ML_DEBUG - edm::LogVerbatim("SimTrackManager") << "TrackingAction addTrack " << id << " added with false flags and deleted"; -#endif - } - EndOfTrack et(aTrack); m_endOfTrackSignal(&et); } diff --git a/SimG4Core/Notification/interface/G4SimEvent.h b/SimG4Core/Notification/interface/G4SimEvent.h index 5219967e16b2d..18371d8439f9d 100644 --- a/SimG4Core/Notification/interface/G4SimEvent.h +++ b/SimG4Core/Notification/interface/G4SimEvent.h @@ -35,10 +35,10 @@ class G4SimEvent { void clear(); private: - const HepMC::GenEvent* hepMCEvent_; - float weight_; - math::XYZTLorentzVectorD collisionPoint_; - int nparam_; + const HepMC::GenEvent* hepMCEvent_{nullptr}; + float weight_{0.f}; + math::XYZTLorentzVectorD collisionPoint_{math::XYZTLorentzVectorD(0.,0.,0.,0.)}; + int nparam_{0}; std::vector param_; std::vector g4tracks_; std::vector g4vertices_; diff --git a/SimG4Core/Notification/interface/G4SimTrack.h b/SimG4Core/Notification/interface/G4SimTrack.h index 2fb4a41b02041..dd19c2c2d3b14 100644 --- a/SimG4Core/Notification/interface/G4SimTrack.h +++ b/SimG4Core/Notification/interface/G4SimTrack.h @@ -3,13 +3,11 @@ #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Math/interface/LorentzVector.h" -#include "FWCore/Utilities/interface/Exception.h" #include "SimG4Core/Notification/interface/TrackWithHistory.h" #include class G4SimTrack { public: - G4SimTrack() {} G4SimTrack(int iid, int ipart, const math::XYZVectorD& ip, double ie) : id_(iid), @@ -18,11 +16,7 @@ class G4SimTrack { ie_(ie), ivert_(-1), igenpart_(-1), - parentID_(-1), - parentMomentum_(math::XYZVectorD(0., 0., 0.)), - tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)), - tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)), - crossedBoundary_(false) {} + parentID_(-1) {} G4SimTrack(int iid, int ipart, const math::XYZVectorD& ip, double ie, int iv, int ig, const math::XYZVectorD& ipmom) : id_(iid), @@ -31,10 +25,7 @@ class G4SimTrack { ie_(ie), ivert_(iv), igenpart_(ig), - parentMomentum_(ipmom), - tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)), - tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)), - crossedBoundary_(false) {} + parentMomentum_(ipmom) {} G4SimTrack(int iid, int ipart, @@ -53,8 +44,7 @@ class G4SimTrack { igenpart_(ig), parentMomentum_(ipmom), tkSurfacePosition_(tkpos), - tkSurfaceMomentum_(tkmom), - crossedBoundary_(false) {} + tkSurfaceMomentum_(tkmom) {} ~G4SimTrack() = default; @@ -94,13 +84,13 @@ class G4SimTrack { int ivert_; int igenpart_; int parentID_; - math::XYZVectorD parentMomentum_; - math::XYZVectorD tkSurfacePosition_; - math::XYZTLorentzVectorD tkSurfaceMomentum_; - bool crossedBoundary_; - int idAtBoundary_; - math::XYZTLorentzVectorF positionAtBoundary_; - math::XYZTLorentzVectorF momentumAtBoundary_; + math::XYZVectorD parentMomentum_{math::XYZVectorD(0.,0.,0.)}; + math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0.,0.,0.)}; + math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0.,0.,0.,0.)}; + bool crossedBoundary_{false}; + int idAtBoundary_{-1}; + math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; + math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; }; #endif diff --git a/SimG4Core/Notification/interface/GenParticleInfo.h b/SimG4Core/Notification/interface/GenParticleInfo.h index 8da448f26a90c..b7cc2ff17ab38 100644 --- a/SimG4Core/Notification/interface/GenParticleInfo.h +++ b/SimG4Core/Notification/interface/GenParticleInfo.h @@ -6,6 +6,7 @@ class GenParticleInfo : public G4VUserPrimaryParticleInformation { public: explicit GenParticleInfo(int id) : id_(id) {} + ~GenParticleInfo() = default; int id() const { return id_; } void Print() const override {} diff --git a/SimG4Core/Notification/interface/NewTrackAction.h b/SimG4Core/Notification/interface/NewTrackAction.h index b1fe4823f374a..6b50459f26ab9 100644 --- a/SimG4Core/Notification/interface/NewTrackAction.h +++ b/SimG4Core/Notification/interface/NewTrackAction.h @@ -14,9 +14,8 @@ class TrackInformation; class NewTrackAction { public: NewTrackAction(); - void primary(const G4Track* aSecondary) const; - void primary(G4Track* aSecondary) const; - void secondary(const G4Track* aSecondary, const G4Track& mother, int) const; + + void primary(G4Track* aPrimary) const; void secondary(G4Track* aSecondary, const G4Track& mother, int) const; private: diff --git a/SimG4Core/Notification/interface/SimTrackManager.h b/SimG4Core/Notification/interface/SimTrackManager.h index e0dd47a5c5d07..aeba0351f0981 100644 --- a/SimG4Core/Notification/interface/SimTrackManager.h +++ b/SimG4Core/Notification/interface/SimTrackManager.h @@ -8,15 +8,9 @@ /**\class SimTrackManager SimTrackManager.h SimG4Core/Notification/interface/SimTrackManager.h Description: Holds tracking information used by the sensitive detectors - - Usage: - + Created: Fri Nov 25 17:36:41 EST 2005 */ -// -// Original Author: -// Created: Fri Nov 25 17:36:41 EST 2005 -// // system include files #include @@ -24,13 +18,12 @@ // 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 class G4SimEvent; +class G4Track; class SimTrackManager { public: @@ -38,33 +31,32 @@ class SimTrackManager { public: bool operator()(TrackWithHistory*& p, const unsigned int& i) const { return p->trackID() < i; } }; - // enum SpecialNumbers {InvalidID = 65535}; /// this map contains association between vertex number and position - typedef std::pair MapVertexPosition; - typedef std::vector > MapVertexPositionVector; - typedef std::map MotherParticleToVertexMap; - typedef MotherParticleToVertexMap VertexMap; + typedef std::pair VertexPosition; + typedef std::vector > VertexPositionVector; + typedef std::map VertexMap; - SimTrackManager(bool iCollapsePrimaryVertices = false); + SimTrackManager(); virtual ~SimTrackManager(); // ---------- const member functions --------------------- - const TrackContainer* trackContainer() const { return m_trksForThisEvent; } + const std::vector* trackContainer() const { return &m_trackContainer; } // ---------- member functions --------------------------- void storeTracks(G4SimEvent* simEvent); void reset(); void deleteTracks(); - void cleanTkCaloStateInfoMap(); - void cleanTracksWithHistory(); - void addTrack(TrackWithHistory* iTrack, bool inHistory, bool withAncestor) { + void addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor) { std::pair thePair(iTrack->trackID(), iTrack->parentID()); idsave.push_back(thePair); if (inHistory) { - m_trksForThisEvent->push_back(iTrack); + m_trackContainer.push_back(iTrack); + const auto& pos = track->GetStep()->GetPostStepPoint()->GetPosition(); + std::pair p(iTrack->trackID(), math::XYZVectorD(pos.x(), pos.y(), pos.z())); + m_endPoints.push_back(p); } if (withAncestor) { std::pair thisPair(iTrack->trackID(), 0); @@ -72,29 +64,21 @@ class SimTrackManager { } } - void addTkCaloStateInfo(uint32_t t, const std::pair& p) { - std::map >::const_iterator it = - mapTkCaloStateInfo.find(t); - - if (it == mapTkCaloStateInfo.end()) { - mapTkCaloStateInfo.insert(std::pair >(t, p)); - } - } - void setCollapsePrimaryVertices(bool iSet) { m_collapsePrimaryVertices = iSet; } int giveMotherNeeded(int i) const { int theResult = 0; - for (unsigned int itr = 0; itr < idsave.size(); itr++) { - if ((idsave[itr]).first == i) { - theResult = (idsave[itr]).second; + for (auto & p : idsave) { + if (p.first == i) { + theResult = p.second; break; } } return theResult; } + bool trackExists(unsigned int i) const { bool flag = false; - for (unsigned int itr = 0; itr < (*m_trksForThisEvent).size(); ++itr) { - if ((*m_trksForThisEvent)[itr]->trackID() == i) { + for (auto & ptr : m_trackContainer) { + if (ptr->trackID() == i) { flag = true; break; } @@ -102,29 +86,17 @@ class SimTrackManager { return flag; } TrackWithHistory* getTrackByID(unsigned int trackID, bool strict = false) 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; + TrackWithHistory* track = nullptr; + for (auto & ptr : m_trackContainer) { + if (ptr->trackID() == trackID) { + track = ptr; break; } } - if (!trackFound) { - if (strict) { - 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 nullptr; - } + if (nullptr == track && strict) { ReportException(trackID); } return track; } + void setLHCTransportLink(const edm::LHCTransportLinkContainer* thisLHCTlink) { theLHCTlink = thisLHCTlink; } // stop default @@ -132,31 +104,31 @@ class SimTrackManager { const SimTrackManager& operator=(const SimTrackManager&) = delete; private: + void saveTrackAndItsBranch(TrackWithHistory*); int getOrCreateVertex(TrackWithHistory*, int, G4SimEvent* simEvent); void cleanVertexMap(); void reallyStoreTracks(G4SimEvent* simEvent); void fillMotherList(); int idSavedTrack(int) const; + void ReportException(unsigned int id) const; // to restore the pre-LHC Transport GenParticle id link to a SimTrack void resetGenID(); // ---------- member data -------------------------------- - TrackContainer* m_trksForThisEvent; - bool m_SaveSimTracks; - MotherParticleToVertexMap m_vertexMap; - int m_nVertices; - bool m_collapsePrimaryVertices; - std::map > mapTkCaloStateInfo; - std::vector > idsave; - std::vector > ancestorList; + int m_nVertices{0}; + unsigned int lastTrack{0}; + unsigned int lastHist{0}; - unsigned int lastTrack; - unsigned int lastHist; + const edm::LHCTransportLinkContainer* theLHCTlink{nullptr}; - const edm::LHCTransportLinkContainer* theLHCTlink; + VertexMap m_vertexMap; + std::vector > idsave; + std::vector > ancestorList; + std::vector > m_endPoints; + std::vector m_trackContainer; }; class trkIDLess { diff --git a/SimG4Core/Notification/interface/TrackInformation.h b/SimG4Core/Notification/interface/TrackInformation.h index 2f6211c88ef69..f54b3c4a363d5 100644 --- a/SimG4Core/Notification/interface/TrackInformation.h +++ b/SimG4Core/Notification/interface/TrackInformation.h @@ -1,17 +1,18 @@ #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" -#include "FWCore/MessageLogger/interface/MessageLogger.h" + +constexpr double invcm = 1.0/CLHEP::cm; +constexpr double invgev = 1.0/CLHEP::GeV; class TrackInformation : public G4VUserTrackInformation { public: - ~TrackInformation() override{}; + ~TrackInformation() override = default; inline void *operator new(size_t); inline void operator delete(void *TrackInformation); @@ -58,14 +59,14 @@ class TrackInformation : public G4VUserTrackInformation { // Boundary crossing variables void setCrossedBoundary(const G4Track *track) { crossedBoundary_ = true; - positionAtBoundary_ = math::XYZTLorentzVectorF(track->GetPosition().x() / CLHEP::cm, - track->GetPosition().y() / CLHEP::cm, - track->GetPosition().z() / CLHEP::cm, + positionAtBoundary_ = math::XYZTLorentzVectorF(track->GetPosition().x()*invcm, + track->GetPosition().y()*invcm, + track->GetPosition().z()*invcm, track->GetGlobalTime()); - momentumAtBoundary_ = math::XYZTLorentzVectorF(track->GetMomentum().x() / CLHEP::GeV, - track->GetMomentum().y() / CLHEP::GeV, - track->GetMomentum().z() / CLHEP::GeV, - track->GetTotalEnergy() / CLHEP::GeV); + momentumAtBoundary_ = math::XYZTLorentzVectorF(track->GetMomentum().x()*invgev, + track->GetMomentum().y()*invgev, + track->GetMomentum().z()*invgev, + track->GetTotalEnergy()*invgev); } bool crossedBoundary() const { return crossedBoundary_; } const math::XYZTLorentzVectorF &getPositionAtBoundary() const { return positionAtBoundary_; } diff --git a/SimG4Core/Notification/interface/TrackWithHistory.h b/SimG4Core/Notification/interface/TrackWithHistory.h index 3070abe99622f..17299b8621dd0 100644 --- a/SimG4Core/Notification/interface/TrackWithHistory.h +++ b/SimG4Core/Notification/interface/TrackWithHistory.h @@ -19,12 +19,12 @@ class TrackWithHistory { * when some of the information is not available yet. */ TrackWithHistory(const G4Track *g4track); - ~TrackWithHistory() {} + ~TrackWithHistory() = default; inline void *operator new(size_t); inline void operator delete(void *TrackWithHistory); - void save() { saved_ = true; } + void setToBeSaved() { saved_ = true; } unsigned int trackID() const { return trackID_; } int particleID() const { return particleID_; } int parentID() const { return parentID_; } @@ -56,18 +56,13 @@ class TrackWithHistory { const math::XYZTLorentzVectorF &getPositionAtBoundary() const { return positionAtBoundary_; } const math::XYZTLorentzVectorF &getMomentumAtBoundary() const { return momentumAtBoundary_; } int getIDAtBoundary() const { return idAtBoundary_; } - /** Internal consistency check (optional). - * Method called at PostUserTrackingAction time, to check - * if the information is consistent with that provided - * to the constructor. - */ - void checkAtEnd(const G4Track *); private: + unsigned int trackID_; int particleID_; int parentID_; - int genParticleID_; + int genParticleID_{-1}; math::XYZVectorD momentum_; double totalEnergy_; math::XYZVectorD vertexPosition_; @@ -77,15 +72,12 @@ class TrackWithHistory { const G4VProcess *creatorProcess_; double weight_; bool storeTrack_; - bool saved_; - - bool isPrimary_; - bool crossedBoundary_; - int idAtBoundary_; - math::XYZTLorentzVectorF positionAtBoundary_; - math::XYZTLorentzVectorF momentumAtBoundary_; + bool saved_{false}; - int extractGenID(const G4Track *gt) const; + bool crossedBoundary_{false}; + int idAtBoundary_{-1}; + math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; + math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; }; extern G4ThreadLocal G4Allocator *fpTrackWithHistoryAllocator; diff --git a/SimG4Core/Notification/src/G4SimEvent.cc b/SimG4Core/Notification/src/G4SimEvent.cc index db2faff234d22..4b50d6813eba0 100644 --- a/SimG4Core/Notification/src/G4SimEvent.cc +++ b/SimG4Core/Notification/src/G4SimEvent.cc @@ -8,21 +8,16 @@ class IdSort { bool operator()(const SimTrack& a, const SimTrack& b) { return a.trackId() < b.trackId(); } }; -G4SimEvent::G4SimEvent() - : hepMCEvent_(nullptr), - weight_(0), - collisionPoint_(math::XYZTLorentzVectorD(0., 0., 0., 0.)), - nparam_(0), - param_(0) { +G4SimEvent::G4SimEvent() { g4vertices_.reserve(2000); g4tracks_.reserve(4000); } -G4SimEvent::~G4SimEvent() { clear(); } +G4SimEvent::~G4SimEvent() { + clear(); +} -void G4SimEvent::clear() { - // per suggestion by Chris Jones, it's faster - // that delete back() and pop_back() +void G4SimEvent::clear() { for (auto& ptr : g4tracks_) { delete ptr; } @@ -34,46 +29,39 @@ void G4SimEvent::clear() { } void G4SimEvent::load(edm::SimTrackContainer& c) const { + const double invgev = 1.0/CLHEP::GeV; for (auto& trk : g4tracks_) { int ip = trk->part(); - math::XYZTLorentzVectorD p( - trk->momentum().x() / GeV, trk->momentum().y() / GeV, trk->momentum().z() / GeV, trk->energy() / GeV); + const math::XYZVectorD& mom = trk->momentum(); + math::XYZTLorentzVectorD p(mom.x()*invgev, mom.y()*invgev, mom.z()*invgev, trk->energy()*invgev); int iv = trk->ivert(); int ig = trk->igenpart(); int id = trk->id(); - math::XYZVectorD tkpos(trk->trackerSurfacePosition().x() / cm, - trk->trackerSurfacePosition().y() / cm, - trk->trackerSurfacePosition().z() / cm); - math::XYZTLorentzVectorD tkmom(trk->trackerSurfaceMomentum().x() / GeV, - trk->trackerSurfaceMomentum().y() / GeV, - trk->trackerSurfaceMomentum().z() / GeV, - trk->trackerSurfaceMomentum().e() / GeV); // ip = particle ID as PDG - // pp = 4-momentum + // pp = 4-momentum in GeV // iv = corresponding G4SimVertex index // ig = corresponding GenParticle index - SimTrack t = SimTrack(ip, p, iv, ig, tkpos, tkmom); + SimTrack t = SimTrack(ip, p, iv, ig, trk->trackerSurfacePosition(), trk->trackerSurfaceMomentum()); t.setTrackId(id); t.setEventId(EncodedEventId(0)); - if (trk->crossedBoundary()) - t.setCrossedBoundaryVars( - trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); + t.setCrossedBoundaryVars(trk->crossedBoundary(), trk->getIDAtBoundary(), + trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); c.push_back(t); } std::stable_sort(c.begin(), c.end(), IdSort()); } void G4SimEvent::load(edm::SimVertexContainer& c) const { + const double invcm = 1.0/CLHEP::cm; + // index of the vertex is needed to make SimVertex object for (unsigned int i = 0; i < g4vertices_.size(); ++i) { G4SimVertex* vtx = g4vertices_[i]; - // - // starting 1_1_0_pre3, SimVertex stores in cm !!! - // - math::XYZVectorD v3(vtx->vertexPosition().x() / cm, vtx->vertexPosition().y() / cm, vtx->vertexPosition().z() / cm); - float t = vtx->vertexGlobalTime() / second; + auto pos = vtx->vertexPosition(); + math::XYZVectorD v3(pos.x()*invcm, pos.y()*invcm, pos.z()*invcm); + float t = vtx->vertexGlobalTime() / CLHEP::second; int iv = vtx->parentIndex(); - // vv = position - // t = global time + // v3 = position in cm + // t = global time in second // iv = index of the parent in the SimEvent SimTrack container (-1 if no parent) SimVertex v = SimVertex(v3, t, iv, i); v.setProcessType(vtx->processType()); diff --git a/SimG4Core/Notification/src/NewTrackAction.cc b/SimG4Core/Notification/src/NewTrackAction.cc index 16c9f21e1744f..5fa71e5261632 100644 --- a/SimG4Core/Notification/src/NewTrackAction.cc +++ b/SimG4Core/Notification/src/NewTrackAction.cc @@ -8,14 +8,8 @@ NewTrackAction::NewTrackAction() {} -void NewTrackAction::primary(const G4Track *aTrack) const { primary(const_cast(aTrack)); } - void NewTrackAction::primary(G4Track *aTrack) const { addUserInfoToPrimary(aTrack); } -void NewTrackAction::secondary(const G4Track *aSecondary, const G4Track &mother, int flag) const { - secondary(const_cast(aSecondary), mother, flag); -} - void NewTrackAction::secondary(G4Track *aSecondary, const G4Track &mother, int flag) const { const TrackInformation *motherInfo = static_cast(mother.GetUserInformation()); addUserInfoToSecondary(aSecondary, *motherInfo, flag); diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index 140c84a7dcd96..2f3678dd730ad 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -18,70 +18,65 @@ #include "SimG4Core/Notification/interface/G4SimTrack.h" #include "SimG4Core/Notification/interface/G4SimVertex.h" #include "SimG4Core/Notification/interface/G4SimEvent.h" +#include "SimG4Core/Notification/interface/TrackInformation.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/Exception.h" #include "G4VProcess.hh" +#include "G4Track.hh" +#include "G4ThreeVector.hh" //#define DebugLog -SimTrackManager::SimTrackManager(bool iCollapsePrimaryVertices) - : m_trksForThisEvent(nullptr), - m_nVertices(0), - m_collapsePrimaryVertices(iCollapsePrimaryVertices), - lastTrack(0), - lastHist(0), - theLHCTlink(nullptr) {} +SimTrackManager::SimTrackManager() { + idsave.reserve(1000); + ancestorList.reserve(1000); + m_trackContainer.reserve(1000); + m_endPoints.reserve(1000); +} SimTrackManager::~SimTrackManager() { - if (m_trksForThisEvent != nullptr) - deleteTracks(); + reset(); } -// -// member functions -// void SimTrackManager::reset() { - if (m_trksForThisEvent == nullptr) { - m_trksForThisEvent = new TrackContainer(); - } else { - for (unsigned int i = 0; i < m_trksForThisEvent->size(); i++) { - delete (*m_trksForThisEvent)[i]; - } - delete m_trksForThisEvent; - m_trksForThisEvent = new TrackContainer(); - } + deleteTracks(); cleanVertexMap(); - cleanTkCaloStateInfoMap(); - std::vector >().swap(idsave); + idsave.clear(); ancestorList.clear(); lastTrack = 0; lastHist = 0; } void SimTrackManager::deleteTracks() { - for (unsigned int i = 0; i < m_trksForThisEvent->size(); ++i) { - delete (*m_trksForThisEvent)[i]; + if(!m_trackContainer.empty()) { + for (auto & ptr : m_trackContainer) { delete ptr; } + m_trackContainer.clear(); + m_endPoints.clear(); } - delete m_trksForThisEvent; - m_trksForThisEvent = nullptr; +} + +void SimTrackManager::cleanVertexMap() { + m_vertexMap.clear(); + m_vertexMap.swap(m_vertexMap); + m_nVertices = 0; } /// this saves a track and all its parents looping over the non ordered vector void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) { - using namespace std; TrackWithHistory* trkH = trkWHist; if (trkH == nullptr) { edm::LogError("SimTrackManager") << " SimTrackManager::saveTrackAndItsBranch got 0 pointer "; throw cms::Exception("SimTrackManager::saveTrackAndItsBranch") << " cannot handle hits for tracking"; } - trkH->save(); + trkH->setToBeSaved(); unsigned int parent = trkH->parentID(); - TrackContainer::const_iterator tk_itr = std::lower_bound( - (*m_trksForThisEvent).begin(), (*m_trksForThisEvent).end(), parent, SimTrackManager::StrictWeakOrdering()); + auto tk_itr = std::lower_bound( + m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering()); - if (tk_itr != m_trksForThisEvent->end() && (*tk_itr)->trackID() == parent) { + if (tk_itr != m_trackContainer.end() && (*tk_itr)->trackID() == parent) { saveTrackAndItsBranch(*tk_itr); } } @@ -91,12 +86,12 @@ void SimTrackManager::storeTracks(G4SimEvent* simEvent) { // fill the map with the final mother-daughter relationship idsave.swap(ancestorList); - stable_sort(idsave.begin(), idsave.end()); + std::stable_sort(idsave.begin(), idsave.end()); std::vector >().swap(ancestorList); // to get a backward compatible order - stable_sort(m_trksForThisEvent->begin(), m_trksForThisEvent->end(), trkIDLess()); + std::stable_sort(m_trackContainer.begin(), m_trackContainer.end(), trkIDLess()); // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value resetGenID(); @@ -107,41 +102,55 @@ void SimTrackManager::storeTracks(G4SimEvent* simEvent) { void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) { // loop over the (now ordered) vector and really save the tracks #ifdef DebugLog - LogDebug("SimTrackManager") << "Inside the reallyStoreTracks method object to be stored = " - << m_trksForThisEvent->size(); + edm::LogVerbatim("SimTrackManager") + << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size(); #endif - for (auto& trkH : *m_trksForThisEvent) { + int nn = m_endPoints.size(); + for (auto& trkH : m_trackContainer) { // at this stage there is one vertex per track, // so the vertex id of track N is also N - int ivertex = -1; - int ig; - - math::XYZVectorD pm(0., 0., 0.); unsigned int iParentID = trkH->parentID(); - for (auto& trk : *m_trksForThisEvent) { - if (trk->trackID() == iParentID) { - pm = trk->momentum(); - break; + int ig = trkH->genParticleID(); + int ivertex = getOrCreateVertex(trkH, iParentID, simEvent); + + auto ptr = trkH; + if(0 < iParentID) { + for (auto& trk : m_trackContainer) { + if (trk->trackID() == iParentID) { + ptr = trk; + break; + } } } - ig = trkH->genParticleID(); - ivertex = getOrCreateVertex(trkH, iParentID, simEvent); - std::map >::const_iterator cit = - mapTkCaloStateInfo.find(trkH->trackID()); - std::pair tcinfo; - if (cit != mapTkCaloStateInfo.end()) { - tcinfo = cit->second; + // Track at surface is the track at intersection point between tracker and calo + // envelops if exist. If not exist the momentum is zero, position is the end of + // the track + const math::XYZVectorD& pm = ptr->momentum(); + math::XYZVectorD spos(0.,0.,0.); + math::XYZTLorentzVectorD smom(0.,0.,0.,0.); + int id = trkH->trackID(); + if(trkH->crossedBoundary()) { + spos = trkH->getPositionAtBoundary(); + smom = trkH->getMomentumAtBoundary(); + } else { + for(int i = 0; itrackID(), + + G4SimTrack* g4simtrack = new G4SimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, - tcinfo.first, - tcinfo.second); + spos, + smom); g4simtrack->copyCrossedBoundaryVars(trkH); simEvent->add(g4simtrack); } @@ -149,7 +158,7 @@ void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) { int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID, G4SimEvent* simEvent) { int parent = -1; - for (auto& trk : *m_trksForThisEvent) { + for (auto& trk : m_trackContainer) { int id = trk->trackID(); if (id == iParentID) { parent = id; @@ -173,22 +182,11 @@ int SimTrackManager::getOrCreateVertex(TrackWithHistory* trkH, int iParentID, G4 ptype = pr->GetProcessSubType(); } simEvent->add(new G4SimVertex(trkH->vertexPosition(), trkH->globalTime(), parent, ptype)); - m_vertexMap[parent].push_back(MapVertexPosition(m_nVertices, trkH->vertexPosition())); + m_vertexMap[parent].push_back(VertexPosition(m_nVertices, trkH->vertexPosition())); ++m_nVertices; return (m_nVertices - 1); } -void SimTrackManager::cleanVertexMap() { - m_vertexMap.clear(); - MotherParticleToVertexMap().swap(m_vertexMap); - m_nVertices = 0; -} - -void SimTrackManager::cleanTkCaloStateInfoMap() { - mapTkCaloStateInfo.clear(); - std::map >().swap(mapTkCaloStateInfo); -} - int SimTrackManager::idSavedTrack(int id) const { int idMother = id; if (id > 0) { @@ -253,22 +251,18 @@ void SimTrackManager::fillMotherList() { lastHist = ancestorList.size(); edm::LogError("SimTrackManager") << " SimTrackManager::fillMotherList track index corrupted"; } - /* - edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: " - << idsave.size() << " saved; ancestor: " << lastHist - << " " << ancestorList.size() << std::endl; +#ifdef DebugLog + edm::LogVerbatim("SimTrackManager") + << "### SimTrackManager::fillMotherList: " + << idsave.size() << " saved; ancestor: " << lastHist << " " << ancestorList.size(); for (unsigned int i = 0; i< idsave.size(); ++i) { - edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first - << " Mother ID = " << (idsave[i]).second << std::endl; + edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first + << " Mother ID = " << (idsave[i]).second; } - */ +#endif for (unsigned int n = lastHist; n < ancestorList.size(); ++n) { int theMotherId = idSavedTrack((ancestorList[n]).first); ancestorList[n].second = theMotherId; - /* - std::cout << " ANC: Track ID = " << (ancestorList[n]).first - << " Mother ID = " << (ancestorList[n]).second << std::endl; - */ #ifdef DebugLog LogDebug("SimTrackManager") << "Track ID = " << (ancestorList[n]).first << " Mother ID = " << (ancestorList[n]).second; @@ -276,12 +270,11 @@ void SimTrackManager::fillMotherList() { } lastHist = ancestorList.size(); - idsave.clear(); } void SimTrackManager::cleanTracksWithHistory() { - if ((*m_trksForThisEvent).empty() && idsave.empty()) { + if (m_trackContainer.empty() && idsave.empty()) { return; } @@ -290,38 +283,36 @@ void SimTrackManager::cleanTracksWithHistory() { << " mother-daughter relationships stored with lastTrack = " << lastTrack; #endif - if (lastTrack > 0 && lastTrack >= (*m_trksForThisEvent).size()) { + if (lastTrack > 0 && lastTrack >= m_trackContainer.size()) { lastTrack = 0; edm::LogError("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory track index corrupted"; } - stable_sort(m_trksForThisEvent->begin() + lastTrack, m_trksForThisEvent->end(), trkIDLess()); - - stable_sort(idsave.begin(), idsave.end()); + std::stable_sort(m_trackContainer.begin() + lastTrack, m_trackContainer.end(), trkIDLess()); + std::stable_sort(idsave.begin(), idsave.end()); #ifdef DebugLog LogDebug("SimTrackManager") << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size() << " tracks with history before branching"; - for (unsigned int it = 0; it < (*m_trksForThisEvent).size(); it++) { + for (unsigned int it = 0; it < m_trackContainer.size(); it++) { LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number " - << (*m_trksForThisEvent)[it]->trackID() << " mother " - << (*m_trksForThisEvent)[it]->parentID() << " status " - << (*m_trksForThisEvent)[it]->saved(); + << m_trackContainer[it]->trackID() << " mother " + << m_trackContainer[it]->parentID() << " status " + << m_trackContainer[it]->saved(); } #endif - for (auto& t : *m_trksForThisEvent) { + for (auto & t : m_trackContainer) { if (t->saved()) { saveTrackAndItsBranch(t); } } unsigned int num = lastTrack; - for (unsigned int it = lastTrack; it < m_trksForThisEvent->size(); ++it) { - TrackWithHistory* t = (*m_trksForThisEvent)[it]; - int g4ID = t->trackID(); - if (t->saved() == true) { - if (it > num) - (*m_trksForThisEvent)[num] = t; + for (unsigned int it = lastTrack; it < m_trackContainer.size(); ++it) { + auto t = m_trackContainer[it]; + int g4ID = m_trackContainer[it]->trackID(); + if (t->saved()) { + if (it > num) { m_trackContainer[num] = t; } ++num; for (auto& xx : idsave) { if (xx.first == g4ID) { @@ -334,31 +325,30 @@ void SimTrackManager::cleanTracksWithHistory() { } } - m_trksForThisEvent->resize(num); + m_trackContainer.resize(num); #ifdef DebugLog - LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << (*m_trksForThisEvent).size() + LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << m_trackContainer.size() << " tracks to be saved persistently"; - for (unsigned int it = 0; it < (*m_trksForThisEvent).size(); ++it) { + for (unsigned int it = 0; it < m_trackContainer.size(); ++it) { LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " - << (*m_trksForThisEvent)[it]->trackID() << " mother " - << (*m_trksForThisEvent)[it]->parentID() << " Status " - << (*m_trksForThisEvent)[it]->saved() << " id " - << (*m_trksForThisEvent)[it]->particleID() - << " E(MeV)= " << (*m_trksForThisEvent)[it]->totalEnergy(); + << m_trackContainer[it]->trackID() << " mother " + << m_trackContainer[it]->parentID() << " Status " + << m_trackContainer[it]->saved() << " id " + << m_trackContainer[it]->particleID() + << " E(MeV)= " << m_trackContainer[it]->totalEnergy(); } #endif fillMotherList(); - - lastTrack = (*m_trksForThisEvent).size(); + lastTrack = m_trackContainer.size(); } void SimTrackManager::resetGenID() { if (theLHCTlink == nullptr) return; - for (auto& trkH : *m_trksForThisEvent) { + for (auto& trkH : m_trackContainer) { int genParticleID = trkH->genParticleID(); if (genParticleID == -1) { continue; @@ -374,3 +364,8 @@ void SimTrackManager::resetGenID() { theLHCTlink = nullptr; } + +void SimTrackManager::ReportException(unsigned int id) const { + throw cms::Exception("Unknown", "SimTrackManager::getTrackByID") + << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size(); +} diff --git a/SimG4Core/Notification/src/TrackWithHistory.cc b/SimG4Core/Notification/src/TrackWithHistory.cc index 24fa2c438e486..0140609849b56 100644 --- a/SimG4Core/Notification/src/TrackWithHistory.cc +++ b/SimG4Core/Notification/src/TrackWithHistory.cc @@ -1,11 +1,14 @@ #include "SimG4Core/Notification/interface/TrackWithHistory.h" #include "SimG4Core/Notification/interface/G4TrackToParticleID.h" -#include "SimG4Core/Notification/interface/TrackInformationExtractor.h" -#include "SimG4Core/Notification/interface/GenParticleInfoExtractor.h" +#include "SimG4Core/Notification/interface/TrackInformation.h" +#include "SimG4Core/Notification/interface/GenParticleInfo.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "G4VProcess.hh" +#include "G4DynamicParticle.hh" +#include "G4PrimaryParticle.hh" +#include "G4ThreeVector.hh" #include @@ -13,76 +16,31 @@ G4ThreadLocal G4Allocator* fpTrackWithHistoryAllocator = nullp //#define DEBUG -TrackWithHistory::TrackWithHistory(const G4Track* g4trk) - : trackID_(0), - particleID_(0), - parentID_(0), - momentum_(math::XYZVectorD(0., 0., 0.)), - totalEnergy_(0), - vertexPosition_(math::XYZVectorD(0., 0., 0.)), - globalTime_(0), - localTime_(0), - properTime_(0), - creatorProcess_(nullptr), - weight_(0), - storeTrack_(false), - saved_(false) { - if (g4trk != nullptr) { - TrackInformationExtractor extractor; - trackID_ = g4trk->GetTrackID(); - particleID_ = G4TrackToParticleID::particleID(g4trk); - parentID_ = g4trk->GetParentID(); - momentum_ = math::XYZVectorD(g4trk->GetMomentum().x(), g4trk->GetMomentum().y(), g4trk->GetMomentum().z()); - totalEnergy_ = g4trk->GetTotalEnergy(); - vertexPosition_ = math::XYZVectorD(g4trk->GetPosition().x(), g4trk->GetPosition().y(), g4trk->GetPosition().z()); - globalTime_ = g4trk->GetGlobalTime(); - localTime_ = g4trk->GetLocalTime(); - properTime_ = g4trk->GetProperTime(); - 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() - weight_ = 10000 * genParticleID_; +TrackWithHistory::TrackWithHistory(const G4Track* g4trk) { + trackID_ = g4trk->GetTrackID(); + particleID_ = G4TrackToParticleID::particleID(g4trk); + parentID_ = g4trk->GetParentID(); + auto mom = g4trk->GetMomentum(); + momentum_ = math::XYZVectorD(mom.x(), mom.y(), mom.z()); + totalEnergy_ = g4trk->GetTotalEnergy(); + auto pos = g4trk->GetPosition(); + vertexPosition_ = math::XYZVectorD(pos.x(), pos.y(), pos.z()); + globalTime_ = g4trk->GetGlobalTime(); + localTime_ = g4trk->GetLocalTime(); + properTime_ = g4trk->GetProperTime(); + creatorProcess_ = g4trk->GetCreatorProcess(); + TrackInformation* trkinfo = static_cast(g4trk->GetUserInformation()); + storeTrack_ = trkinfo->storeTrack(); + auto vgprimary = g4trk->GetDynamicParticle()->GetPrimaryParticle(); + if (vgprimary != nullptr) { + auto priminfo = static_cast(vgprimary->GetUserInformation()); + if(nullptr != priminfo) { genParticleID_ = priminfo->id(); } + } + // V.I. weight is computed in the same way as before + // without usage of G4Track::GetWeight() + weight_ = 10000 * genParticleID_; #ifdef DEBUG - LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother " - << parentID_; + LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother " + << parentID_; #endif - } -} - -void TrackWithHistory::checkAtEnd(const G4Track* gt) { - math::XYZVectorD vposdir(gt->GetVertexPosition().x(), gt->GetVertexPosition().y(), gt->GetVertexPosition().z()); - math::XYZVectorD vmomdir( - gt->GetVertexMomentumDirection().x(), gt->GetVertexMomentumDirection().y(), gt->GetVertexMomentumDirection().z()); - bool ok = true; - double epsilon = 1.e-6; - double eps2 = epsilon * epsilon; - if ((vertexPosition_ - vposdir).Mag2() > eps2) { - edm::LogWarning("TrackInformation") << "TrackWithHistory vertex position check failed" - << "\nAt construction: " << vertexPosition_ << "\nAt end: " << vposdir; - ok = false; - } - math::XYZVectorD dirDiff = momentum_.Unit() - vmomdir; - if (dirDiff.Mag2() > eps2 && momentum_.Unit().R() > eps2) { - edm::LogWarning("TrackInformation") << "TrackWithHistory momentum direction check failed" - << "\nAt construction: " << momentum_.Unit() - << "\nAt end: " << vmomdir; - ok = false; - } - - if (!ok) - G4Exception("TrackWithHistory::checkAtEnd()", "mc001", FatalException, "check at track end failed"); -} - -int TrackWithHistory::extractGenID(const G4Track* gt) const { - void* vgprimary = gt->GetDynamicParticle()->GetPrimaryParticle(); - if (vgprimary == nullptr) - return -1; - // replace old-style cast with appropriate new-style cast... - G4PrimaryParticle* gprimary = (G4PrimaryParticle*)vgprimary; - GenParticleInfoExtractor ext; - return ext(gprimary).id(); } From 2ad4ac7e669c702d618d442b00d9db820012ca42 Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 24 Mar 2023 17:56:41 +0100 Subject: [PATCH 2/5] fode format --- SimG4CMS/Calo/src/CaloSD.cc | 6 +- SimG4Core/Application/src/SteppingAction.cc | 4 +- SimG4Core/Application/src/TrackingAction.cc | 12 +-- SimG4Core/Notification/interface/G4SimEvent.h | 2 +- SimG4Core/Notification/interface/G4SimTrack.h | 27 ++---- .../Notification/interface/GenParticleInfo.h | 2 +- .../Notification/interface/SimTrackManager.h | 11 +-- .../Notification/interface/TrackInformation.h | 18 ++-- .../Notification/interface/TrackWithHistory.h | 5 +- SimG4Core/Notification/src/G4SimEvent.cc | 18 ++-- SimG4Core/Notification/src/SimTrackManager.cc | 82 ++++++++----------- .../Notification/src/TrackWithHistory.cc | 6 +- 12 files changed, 86 insertions(+), 107 deletions(-) diff --git a/SimG4CMS/Calo/src/CaloSD.cc b/SimG4CMS/Calo/src/CaloSD.cc index ecaca60900821..2bb6ac2130326 100644 --- a/SimG4CMS/Calo/src/CaloSD.cc +++ b/SimG4CMS/Calo/src/CaloSD.cc @@ -705,10 +705,10 @@ void CaloSD::update(const EndOfTrack* trk) { if (!trksForThisEvent->empty()) { TrackWithHistory* trkH = trksForThisEvent->back(); if (trkH->trackID() == (unsigned int)(id)) { - tkMap[id] = trkH; + tkMap[id] = trkH; #ifdef EDM_ML_DEBUG - edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size " - << trksForThisEvent->size() << " with ID " << trkH->trackID(); + edm::LogVerbatim("CaloSim") << "CaloSD: get track " << it << " from Container of size " + << trksForThisEvent->size() << " with ID " << trkH->trackID(); #endif } } diff --git a/SimG4Core/Application/src/SteppingAction.cc b/SimG4Core/Application/src/SteppingAction.cc index 0f8fd7c40af8f..10b39543770e2 100644 --- a/SimG4Core/Application/src/SteppingAction.cc +++ b/SimG4Core/Application/src/SteppingAction.cc @@ -163,7 +163,9 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) { if (sAlive == tstat || sVeryForward == tstat) { if (preStep->GetPhysicalVolume() == tracker && postStep->GetPhysicalVolume() == calo) { TrackInformation* trkinfo = static_cast(theTrack->GetUserInformation()); - if(!trkinfo->crossedBoundary()) { trkinfo->setCrossedBoundary(theTrack); } + if (!trkinfo->crossedBoundary()) { + trkinfo->setCrossedBoundary(theTrack); + } } } else { theTrack->SetTrackStatus(fStopAndKill); diff --git a/SimG4Core/Application/src/TrackingAction.cc b/SimG4Core/Application/src/TrackingAction.cc index 3894e3753d982..cc3bd47e5c3ab 100644 --- a/SimG4Core/Application/src/TrackingAction.cc +++ b/SimG4Core/Application/src/TrackingAction.cc @@ -83,7 +83,9 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { currentTrack_->setCrossedBoundaryPosMom(id, trkInfo_->getPositionAtBoundary(), trkInfo_->getMomentumAtBoundary()); ok = (ok || saveCaloBoundaryInformation_ || doFineCalo_); } - if(ok) { currentTrack_->setToBeSaved(); } + if (ok) { + currentTrack_->setToBeSaved(); + } bool withAncestor = (trkInfo_->getIDonCaloSurface() == id || trkInfo_->isAncestor()); bool isInHistory = trkInfo_->isInHistory(); @@ -91,11 +93,9 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { trackManager_->addTrack(currentTrack_, aTrack, isInHistory, withAncestor); #ifdef EDM_ML_DEBUG - edm::LogVerbatim("TrackingAction") << "TrackingAction end track=" << id - << " " << aTrack->GetDefinition()->GetParticleName() - << " ansestor= " << withAncestor - << " saved= " << currentTrack_->saved() - << " end point " << aTrack->GetPosition(); + edm::LogVerbatim("TrackingAction") << "TrackingAction end track=" << id << " " + << aTrack->GetDefinition()->GetParticleName() << " ansestor= " << withAncestor + << " saved= " << currentTrack_->saved() << " end point " << aTrack->GetPosition(); #endif EndOfTrack et(aTrack); diff --git a/SimG4Core/Notification/interface/G4SimEvent.h b/SimG4Core/Notification/interface/G4SimEvent.h index 18371d8439f9d..fdffd790f0230 100644 --- a/SimG4Core/Notification/interface/G4SimEvent.h +++ b/SimG4Core/Notification/interface/G4SimEvent.h @@ -37,7 +37,7 @@ class G4SimEvent { private: const HepMC::GenEvent* hepMCEvent_{nullptr}; float weight_{0.f}; - math::XYZTLorentzVectorD collisionPoint_{math::XYZTLorentzVectorD(0.,0.,0.,0.)}; + math::XYZTLorentzVectorD collisionPoint_{math::XYZTLorentzVectorD(0., 0., 0., 0.)}; int nparam_{0}; std::vector param_; std::vector g4tracks_; diff --git a/SimG4Core/Notification/interface/G4SimTrack.h b/SimG4Core/Notification/interface/G4SimTrack.h index dd19c2c2d3b14..b97354f00b706 100644 --- a/SimG4Core/Notification/interface/G4SimTrack.h +++ b/SimG4Core/Notification/interface/G4SimTrack.h @@ -8,24 +8,11 @@ class G4SimTrack { public: - G4SimTrack(int iid, int ipart, const math::XYZVectorD& ip, double ie) - : id_(iid), - ipart_(ipart), - ip_(ip), - ie_(ie), - ivert_(-1), - igenpart_(-1), - parentID_(-1) {} + : id_(iid), ipart_(ipart), ip_(ip), ie_(ie), ivert_(-1), igenpart_(-1), parentID_(-1) {} G4SimTrack(int iid, int ipart, const math::XYZVectorD& ip, double ie, int iv, int ig, const math::XYZVectorD& ipmom) - : id_(iid), - ipart_(ipart), - ip_(ip), - ie_(ie), - ivert_(iv), - igenpart_(ig), - parentMomentum_(ipmom) {} + : id_(iid), ipart_(ipart), ip_(ip), ie_(ie), ivert_(iv), igenpart_(ig), parentMomentum_(ipmom) {} G4SimTrack(int iid, int ipart, @@ -84,13 +71,13 @@ class G4SimTrack { int ivert_; int igenpart_; int parentID_; - math::XYZVectorD parentMomentum_{math::XYZVectorD(0.,0.,0.)}; - math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0.,0.,0.)}; - math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0.,0.,0.,0.)}; + math::XYZVectorD parentMomentum_{math::XYZVectorD(0., 0., 0.)}; + math::XYZVectorD tkSurfacePosition_{math::XYZVectorD(0., 0., 0.)}; + math::XYZTLorentzVectorD tkSurfaceMomentum_{math::XYZTLorentzVectorD(0., 0., 0., 0.)}; bool crossedBoundary_{false}; int idAtBoundary_{-1}; - math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; - math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; + math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; + math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; }; #endif diff --git a/SimG4Core/Notification/interface/GenParticleInfo.h b/SimG4Core/Notification/interface/GenParticleInfo.h index b7cc2ff17ab38..eb9ceb86aadfb 100644 --- a/SimG4Core/Notification/interface/GenParticleInfo.h +++ b/SimG4Core/Notification/interface/GenParticleInfo.h @@ -6,7 +6,7 @@ class GenParticleInfo : public G4VUserPrimaryParticleInformation { public: explicit GenParticleInfo(int id) : id_(id) {} - ~GenParticleInfo() = default; + ~GenParticleInfo() override = default; int id() const { return id_; } void Print() const override {} diff --git a/SimG4Core/Notification/interface/SimTrackManager.h b/SimG4Core/Notification/interface/SimTrackManager.h index aeba0351f0981..0adf91a13477f 100644 --- a/SimG4Core/Notification/interface/SimTrackManager.h +++ b/SimG4Core/Notification/interface/SimTrackManager.h @@ -66,7 +66,7 @@ class SimTrackManager { int giveMotherNeeded(int i) const { int theResult = 0; - for (auto & p : idsave) { + for (auto& p : idsave) { if (p.first == i) { theResult = p.second; break; @@ -77,7 +77,7 @@ class SimTrackManager { bool trackExists(unsigned int i) const { bool flag = false; - for (auto & ptr : m_trackContainer) { + for (auto& ptr : m_trackContainer) { if (ptr->trackID() == i) { flag = true; break; @@ -87,13 +87,15 @@ class SimTrackManager { } TrackWithHistory* getTrackByID(unsigned int trackID, bool strict = false) const { TrackWithHistory* track = nullptr; - for (auto & ptr : m_trackContainer) { + for (auto& ptr : m_trackContainer) { if (ptr->trackID() == trackID) { track = ptr; break; } } - if (nullptr == track && strict) { ReportException(trackID); } + if (nullptr == track && strict) { + ReportException(trackID); + } return track; } @@ -104,7 +106,6 @@ class SimTrackManager { const SimTrackManager& operator=(const SimTrackManager&) = delete; private: - void saveTrackAndItsBranch(TrackWithHistory*); int getOrCreateVertex(TrackWithHistory*, int, G4SimEvent* simEvent); void cleanVertexMap(); diff --git a/SimG4Core/Notification/interface/TrackInformation.h b/SimG4Core/Notification/interface/TrackInformation.h index f54b3c4a363d5..009c5d4ab6e03 100644 --- a/SimG4Core/Notification/interface/TrackInformation.h +++ b/SimG4Core/Notification/interface/TrackInformation.h @@ -7,8 +7,8 @@ #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Math/interface/LorentzVector.h" -constexpr double invcm = 1.0/CLHEP::cm; -constexpr double invgev = 1.0/CLHEP::GeV; +constexpr double invcm = 1.0 / CLHEP::cm; +constexpr double invgev = 1.0 / CLHEP::GeV; class TrackInformation : public G4VUserTrackInformation { public: @@ -59,14 +59,14 @@ class TrackInformation : public G4VUserTrackInformation { // Boundary crossing variables void setCrossedBoundary(const G4Track *track) { crossedBoundary_ = true; - positionAtBoundary_ = math::XYZTLorentzVectorF(track->GetPosition().x()*invcm, - track->GetPosition().y()*invcm, - track->GetPosition().z()*invcm, + positionAtBoundary_ = math::XYZTLorentzVectorF(track->GetPosition().x() * invcm, + track->GetPosition().y() * invcm, + track->GetPosition().z() * invcm, track->GetGlobalTime()); - momentumAtBoundary_ = math::XYZTLorentzVectorF(track->GetMomentum().x()*invgev, - track->GetMomentum().y()*invgev, - track->GetMomentum().z()*invgev, - track->GetTotalEnergy()*invgev); + momentumAtBoundary_ = math::XYZTLorentzVectorF(track->GetMomentum().x() * invgev, + track->GetMomentum().y() * invgev, + track->GetMomentum().z() * invgev, + track->GetTotalEnergy() * invgev); } bool crossedBoundary() const { return crossedBoundary_; } const math::XYZTLorentzVectorF &getPositionAtBoundary() const { return positionAtBoundary_; } diff --git a/SimG4Core/Notification/interface/TrackWithHistory.h b/SimG4Core/Notification/interface/TrackWithHistory.h index 17299b8621dd0..f9b64c7f9a189 100644 --- a/SimG4Core/Notification/interface/TrackWithHistory.h +++ b/SimG4Core/Notification/interface/TrackWithHistory.h @@ -58,7 +58,6 @@ class TrackWithHistory { int getIDAtBoundary() const { return idAtBoundary_; } private: - unsigned int trackID_; int particleID_; int parentID_; @@ -76,8 +75,8 @@ class TrackWithHistory { bool crossedBoundary_{false}; int idAtBoundary_{-1}; - math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; - math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f,0.f,0.f,0.f)}; + math::XYZTLorentzVectorF positionAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; + math::XYZTLorentzVectorF momentumAtBoundary_{math::XYZTLorentzVectorF(0.f, 0.f, 0.f, 0.f)}; }; extern G4ThreadLocal G4Allocator *fpTrackWithHistoryAllocator; diff --git a/SimG4Core/Notification/src/G4SimEvent.cc b/SimG4Core/Notification/src/G4SimEvent.cc index 4b50d6813eba0..1a7680da71acd 100644 --- a/SimG4Core/Notification/src/G4SimEvent.cc +++ b/SimG4Core/Notification/src/G4SimEvent.cc @@ -13,11 +13,9 @@ G4SimEvent::G4SimEvent() { g4tracks_.reserve(4000); } -G4SimEvent::~G4SimEvent() { - clear(); -} +G4SimEvent::~G4SimEvent() { clear(); } -void G4SimEvent::clear() { +void G4SimEvent::clear() { for (auto& ptr : g4tracks_) { delete ptr; } @@ -29,11 +27,11 @@ void G4SimEvent::clear() { } void G4SimEvent::load(edm::SimTrackContainer& c) const { - const double invgev = 1.0/CLHEP::GeV; + const double invgev = 1.0 / CLHEP::GeV; for (auto& trk : g4tracks_) { int ip = trk->part(); const math::XYZVectorD& mom = trk->momentum(); - math::XYZTLorentzVectorD p(mom.x()*invgev, mom.y()*invgev, mom.z()*invgev, trk->energy()*invgev); + math::XYZTLorentzVectorD p(mom.x() * invgev, mom.y() * invgev, mom.z() * invgev, trk->energy() * invgev); int iv = trk->ivert(); int ig = trk->igenpart(); int id = trk->id(); @@ -44,20 +42,20 @@ void G4SimEvent::load(edm::SimTrackContainer& c) const { SimTrack t = SimTrack(ip, p, iv, ig, trk->trackerSurfacePosition(), trk->trackerSurfaceMomentum()); t.setTrackId(id); t.setEventId(EncodedEventId(0)); - t.setCrossedBoundaryVars(trk->crossedBoundary(), trk->getIDAtBoundary(), - trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); + t.setCrossedBoundaryVars( + trk->crossedBoundary(), trk->getIDAtBoundary(), trk->getPositionAtBoundary(), trk->getMomentumAtBoundary()); c.push_back(t); } std::stable_sort(c.begin(), c.end(), IdSort()); } void G4SimEvent::load(edm::SimVertexContainer& c) const { - const double invcm = 1.0/CLHEP::cm; + const double invcm = 1.0 / CLHEP::cm; // index of the vertex is needed to make SimVertex object for (unsigned int i = 0; i < g4vertices_.size(); ++i) { G4SimVertex* vtx = g4vertices_[i]; auto pos = vtx->vertexPosition(); - math::XYZVectorD v3(pos.x()*invcm, pos.y()*invcm, pos.z()*invcm); + math::XYZVectorD v3(pos.x() * invcm, pos.y() * invcm, pos.z() * invcm); float t = vtx->vertexGlobalTime() / CLHEP::second; int iv = vtx->parentIndex(); // v3 = position in cm diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index 2f3678dd730ad..104fd5c30055f 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -36,9 +36,7 @@ SimTrackManager::SimTrackManager() { m_endPoints.reserve(1000); } -SimTrackManager::~SimTrackManager() { - reset(); -} +SimTrackManager::~SimTrackManager() { reset(); } void SimTrackManager::reset() { deleteTracks(); @@ -50,8 +48,10 @@ void SimTrackManager::reset() { } void SimTrackManager::deleteTracks() { - if(!m_trackContainer.empty()) { - for (auto & ptr : m_trackContainer) { delete ptr; } + if (!m_trackContainer.empty()) { + for (auto& ptr : m_trackContainer) { + delete ptr; + } m_trackContainer.clear(); m_endPoints.clear(); } @@ -73,8 +73,8 @@ void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) { trkH->setToBeSaved(); unsigned int parent = trkH->parentID(); - auto tk_itr = std::lower_bound( - m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering()); + auto tk_itr = + std::lower_bound(m_trackContainer.begin(), m_trackContainer.end(), parent, SimTrackManager::StrictWeakOrdering()); if (tk_itr != m_trackContainer.end() && (*tk_itr)->trackID() == parent) { saveTrackAndItsBranch(*tk_itr); @@ -102,8 +102,7 @@ void SimTrackManager::storeTracks(G4SimEvent* simEvent) { void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) { // loop over the (now ordered) vector and really save the tracks #ifdef DebugLog - edm::LogVerbatim("SimTrackManager") - << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size(); + edm::LogVerbatim("SimTrackManager") << "reallyStoreTracks() NtracksWithHistory= " << m_trackContainer->size(); #endif int nn = m_endPoints.size(); @@ -115,42 +114,35 @@ void SimTrackManager::reallyStoreTracks(G4SimEvent* simEvent) { int ivertex = getOrCreateVertex(trkH, iParentID, simEvent); auto ptr = trkH; - if(0 < iParentID) { + if (0 < iParentID) { for (auto& trk : m_trackContainer) { - if (trk->trackID() == iParentID) { - ptr = trk; - break; - } + if (trk->trackID() == iParentID) { + ptr = trk; + break; + } } } // Track at surface is the track at intersection point between tracker and calo // envelops if exist. If not exist the momentum is zero, position is the end of // the track const math::XYZVectorD& pm = ptr->momentum(); - math::XYZVectorD spos(0.,0.,0.); - math::XYZTLorentzVectorD smom(0.,0.,0.,0.); + math::XYZVectorD spos(0., 0., 0.); + math::XYZTLorentzVectorD smom(0., 0., 0., 0.); int id = trkH->trackID(); - if(trkH->crossedBoundary()) { + if (trkH->crossedBoundary()) { spos = trkH->getPositionAtBoundary(); smom = trkH->getMomentumAtBoundary(); } else { - for(int i = 0; iparticleID(), - trkH->momentum(), - trkH->totalEnergy(), - ivertex, - ig, - pm, - spos, - smom); + G4SimTrack* g4simtrack = + new G4SimTrack(id, trkH->particleID(), trkH->momentum(), trkH->totalEnergy(), ivertex, ig, pm, spos, smom); g4simtrack->copyCrossedBoundaryVars(trkH); simEvent->add(g4simtrack); } @@ -252,11 +244,10 @@ void SimTrackManager::fillMotherList() { edm::LogError("SimTrackManager") << " SimTrackManager::fillMotherList track index corrupted"; } #ifdef DebugLog - edm::LogVerbatim("SimTrackManager") - << "### SimTrackManager::fillMotherList: " - << idsave.size() << " saved; ancestor: " << lastHist << " " << ancestorList.size(); - for (unsigned int i = 0; i< idsave.size(); ++i) { - edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first + edm::LogVerbatim("SimTrackManager") << "### SimTrackManager::fillMotherList: " << idsave.size() + << " saved; ancestor: " << lastHist << " " << ancestorList.size(); + for (unsigned int i = 0; i < idsave.size(); ++i) { + edm::LogVerbatim("SimTrackManager") << " ISV: Track ID = " << (idsave[i]).first << " Mother ID = " << (idsave[i]).second; } #endif @@ -296,13 +287,12 @@ void SimTrackManager::cleanTracksWithHistory() { << " tracks with history before branching"; for (unsigned int it = 0; it < m_trackContainer.size(); it++) { LogDebug("SimTrackManager") << " 1 - Track in position " << it << " G4 track number " - << m_trackContainer[it]->trackID() << " mother " - << m_trackContainer[it]->parentID() << " status " - << m_trackContainer[it]->saved(); + << m_trackContainer[it]->trackID() << " mother " << m_trackContainer[it]->parentID() + << " status " << m_trackContainer[it]->saved(); } #endif - for (auto & t : m_trackContainer) { + for (auto& t : m_trackContainer) { if (t->saved()) { saveTrackAndItsBranch(t); } @@ -312,7 +302,9 @@ void SimTrackManager::cleanTracksWithHistory() { auto t = m_trackContainer[it]; int g4ID = m_trackContainer[it]->trackID(); if (t->saved()) { - if (it > num) { m_trackContainer[num] = t; } + if (it > num) { + m_trackContainer[num] = t; + } ++num; for (auto& xx : idsave) { if (xx.first == g4ID) { @@ -331,11 +323,9 @@ void SimTrackManager::cleanTracksWithHistory() { LogDebug("SimTrackManager") << " AFTER CLEANING, I GET " << m_trackContainer.size() << " tracks to be saved persistently"; for (unsigned int it = 0; it < m_trackContainer.size(); ++it) { - LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " - << m_trackContainer[it]->trackID() << " mother " - << m_trackContainer[it]->parentID() << " Status " - << m_trackContainer[it]->saved() << " id " - << m_trackContainer[it]->particleID() + LogDebug("SimTrackManager") << " Track in position " << it << " G4 track number " << m_trackContainer[it]->trackID() + << " mother " << m_trackContainer[it]->parentID() << " Status " + << m_trackContainer[it]->saved() << " id " << m_trackContainer[it]->particleID() << " E(MeV)= " << m_trackContainer[it]->totalEnergy(); } #endif @@ -367,5 +357,5 @@ void SimTrackManager::resetGenID() { void SimTrackManager::ReportException(unsigned int id) const { throw cms::Exception("Unknown", "SimTrackManager::getTrackByID") - << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size(); + << "Fail to get track " << id << " from SimTrackManager, container size= " << m_trackContainer.size(); } diff --git a/SimG4Core/Notification/src/TrackWithHistory.cc b/SimG4Core/Notification/src/TrackWithHistory.cc index 0140609849b56..5f17cba753087 100644 --- a/SimG4Core/Notification/src/TrackWithHistory.cc +++ b/SimG4Core/Notification/src/TrackWithHistory.cc @@ -34,13 +34,15 @@ TrackWithHistory::TrackWithHistory(const G4Track* g4trk) { auto vgprimary = g4trk->GetDynamicParticle()->GetPrimaryParticle(); if (vgprimary != nullptr) { auto priminfo = static_cast(vgprimary->GetUserInformation()); - if(nullptr != priminfo) { genParticleID_ = priminfo->id(); } + if (nullptr != priminfo) { + genParticleID_ = priminfo->id(); + } } // V.I. weight is computed in the same way as before // without usage of G4Track::GetWeight() weight_ = 10000 * genParticleID_; #ifdef DEBUG LogDebug("TrackInformation") << " TrackWithHistory : created history for " << trackID_ << " with mother " - << parentID_; + << parentID_; #endif } From 4f913727cb4684bd957373fbd122bbfb65d664bb Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 24 Mar 2023 18:25:45 +0100 Subject: [PATCH 3/5] fixed units --- .../Notification/interface/SimTrackManager.h | 21 ++++--------------- SimG4Core/Notification/src/SimTrackManager.cc | 18 ++++++++++++++++ 2 files changed, 22 insertions(+), 17 deletions(-) diff --git a/SimG4Core/Notification/interface/SimTrackManager.h b/SimG4Core/Notification/interface/SimTrackManager.h index 0adf91a13477f..551963b844977 100644 --- a/SimG4Core/Notification/interface/SimTrackManager.h +++ b/SimG4Core/Notification/interface/SimTrackManager.h @@ -31,7 +31,7 @@ class SimTrackManager { public: bool operator()(TrackWithHistory*& p, const unsigned int& i) const { return p->trackID() < i; } }; - /// this map contains association between vertex number and position + typedef std::pair VertexPosition; typedef std::vector > VertexPositionVector; typedef std::map VertexMap; @@ -39,30 +39,15 @@ class SimTrackManager { SimTrackManager(); virtual ~SimTrackManager(); - // ---------- const member functions --------------------- const std::vector* trackContainer() const { return &m_trackContainer; } - // ---------- member functions --------------------------- void storeTracks(G4SimEvent* simEvent); void reset(); void deleteTracks(); void cleanTracksWithHistory(); - void addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor) { - std::pair thePair(iTrack->trackID(), iTrack->parentID()); - idsave.push_back(thePair); - if (inHistory) { - m_trackContainer.push_back(iTrack); - const auto& pos = track->GetStep()->GetPostStepPoint()->GetPosition(); - std::pair p(iTrack->trackID(), math::XYZVectorD(pos.x(), pos.y(), pos.z())); - m_endPoints.push_back(p); - } - if (withAncestor) { - std::pair thisPair(iTrack->trackID(), 0); - ancestorList.push_back(thisPair); - } - } + void addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor); int giveMotherNeeded(int i) const { int theResult = 0; @@ -85,6 +70,7 @@ class SimTrackManager { } return flag; } + TrackWithHistory* getTrackByID(unsigned int trackID, bool strict = false) const { TrackWithHistory* track = nullptr; for (auto& ptr : m_trackContainer) { @@ -106,6 +92,7 @@ class SimTrackManager { const SimTrackManager& operator=(const SimTrackManager&) = delete; private: + void saveTrackAndItsBranch(TrackWithHistory*); int getOrCreateVertex(TrackWithHistory*, int, G4SimEvent* simEvent); void cleanVertexMap(); diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index 104fd5c30055f..cbcf3101bca6c 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -26,6 +26,7 @@ #include "G4VProcess.hh" #include "G4Track.hh" #include "G4ThreeVector.hh" +#include "G4SystemOfUnits.hh" //#define DebugLog @@ -63,6 +64,23 @@ void SimTrackManager::cleanVertexMap() { m_nVertices = 0; } +void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, + bool inHistory, bool withAncestor) { + std::pair thePair(iTrack->trackID(), iTrack->parentID()); + idsave.push_back(thePair); + if (inHistory) { + m_trackContainer.push_back(iTrack); + const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition(); + const double invcm = 1.0/CLHEP::cm; + std::pair p(iTrack->trackID(), math::XYZVectorD(v.x()*invcm, v.y()*invcm, v.z()*invcm)); + m_endPoints.push_back(p); + } + if (withAncestor) { + std::pair thisPair(iTrack->trackID(), 0); + ancestorList.push_back(thisPair); + } +} + /// this saves a track and all its parents looping over the non ordered vector void SimTrackManager::saveTrackAndItsBranch(TrackWithHistory* trkWHist) { TrackWithHistory* trkH = trkWHist; From da69d4fbfdee2e67df1955e34b779d65a6b50c41 Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 24 Mar 2023 18:30:41 +0100 Subject: [PATCH 4/5] code format --- SimG4Core/Notification/interface/SimTrackManager.h | 1 - SimG4Core/Notification/src/SimTrackManager.cc | 8 ++++---- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/SimG4Core/Notification/interface/SimTrackManager.h b/SimG4Core/Notification/interface/SimTrackManager.h index 551963b844977..83833442a817e 100644 --- a/SimG4Core/Notification/interface/SimTrackManager.h +++ b/SimG4Core/Notification/interface/SimTrackManager.h @@ -92,7 +92,6 @@ class SimTrackManager { const SimTrackManager& operator=(const SimTrackManager&) = delete; private: - void saveTrackAndItsBranch(TrackWithHistory*); int getOrCreateVertex(TrackWithHistory*, int, G4SimEvent* simEvent); void cleanVertexMap(); diff --git a/SimG4Core/Notification/src/SimTrackManager.cc b/SimG4Core/Notification/src/SimTrackManager.cc index cbcf3101bca6c..2e494d1bfca68 100644 --- a/SimG4Core/Notification/src/SimTrackManager.cc +++ b/SimG4Core/Notification/src/SimTrackManager.cc @@ -64,15 +64,15 @@ void SimTrackManager::cleanVertexMap() { m_nVertices = 0; } -void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, - bool inHistory, bool withAncestor) { +void SimTrackManager::addTrack(TrackWithHistory* iTrack, const G4Track* track, bool inHistory, bool withAncestor) { std::pair thePair(iTrack->trackID(), iTrack->parentID()); idsave.push_back(thePair); if (inHistory) { m_trackContainer.push_back(iTrack); const auto& v = track->GetStep()->GetPostStepPoint()->GetPosition(); - const double invcm = 1.0/CLHEP::cm; - std::pair p(iTrack->trackID(), math::XYZVectorD(v.x()*invcm, v.y()*invcm, v.z()*invcm)); + const double invcm = 1.0 / CLHEP::cm; + std::pair p(iTrack->trackID(), + math::XYZVectorD(v.x() * invcm, v.y() * invcm, v.z() * invcm)); m_endPoints.push_back(p); } if (withAncestor) { From 0ce793a36637c35b8d353dafb10c5a1bdabeb482 Mon Sep 17 00:00:00 2001 From: Vladimir Date: Fri, 24 Mar 2023 19:23:09 +0100 Subject: [PATCH 5/5] fixed constant shadowing --- .../Notification/interface/TrackInformation.h | 15 +-------------- SimG4Core/Notification/src/TrackInformation.cc | 18 +++++++++++++++++- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/SimG4Core/Notification/interface/TrackInformation.h b/SimG4Core/Notification/interface/TrackInformation.h index 009c5d4ab6e03..62fa4a5680499 100644 --- a/SimG4Core/Notification/interface/TrackInformation.h +++ b/SimG4Core/Notification/interface/TrackInformation.h @@ -7,9 +7,6 @@ #include "DataFormats/Math/interface/Vector3D.h" #include "DataFormats/Math/interface/LorentzVector.h" -constexpr double invcm = 1.0 / CLHEP::cm; -constexpr double invgev = 1.0 / CLHEP::GeV; - class TrackInformation : public G4VUserTrackInformation { public: ~TrackInformation() override = default; @@ -57,17 +54,7 @@ class TrackInformation : public G4VUserTrackInformation { void setCaloSurfaceParticleP(double p) { caloSurfaceParticleP_ = p; } // Boundary crossing variables - void setCrossedBoundary(const G4Track *track) { - crossedBoundary_ = true; - positionAtBoundary_ = math::XYZTLorentzVectorF(track->GetPosition().x() * invcm, - track->GetPosition().y() * invcm, - track->GetPosition().z() * invcm, - track->GetGlobalTime()); - momentumAtBoundary_ = math::XYZTLorentzVectorF(track->GetMomentum().x() * invgev, - track->GetMomentum().y() * invgev, - track->GetMomentum().z() * invgev, - track->GetTotalEnergy() * invgev); - } + void setCrossedBoundary(const G4Track *track); bool crossedBoundary() const { return crossedBoundary_; } const math::XYZTLorentzVectorF &getPositionAtBoundary() const { return positionAtBoundary_; } const math::XYZTLorentzVectorF &getMomentumAtBoundary() const { return momentumAtBoundary_; } diff --git a/SimG4Core/Notification/src/TrackInformation.cc b/SimG4Core/Notification/src/TrackInformation.cc index f17e947ad5442..318074d820c1a 100644 --- a/SimG4Core/Notification/src/TrackInformation.cc +++ b/SimG4Core/Notification/src/TrackInformation.cc @@ -1,9 +1,25 @@ #include "SimG4Core/Notification/interface/TrackInformation.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "G4ThreeVector.hh" +#include "G4SystemOfUnits.hh" + #include -G4ThreadLocal G4Allocator *fpTrackInformationAllocator = nullptr; +G4ThreadLocal G4Allocator* fpTrackInformationAllocator = nullptr; + +void TrackInformation::setCrossedBoundary(const G4Track* track) { + const double invcm = 1.0 / CLHEP::cm; + const double invgev = 1.0 / CLHEP::GeV; + const double invsec = 1.0 / CLHEP::second; + crossedBoundary_ = true; + const G4ThreeVector& v = track->GetPosition(); + positionAtBoundary_ = + math::XYZTLorentzVectorF(v.x() * invcm, v.y() * invcm, v.z() * invcm, track->GetGlobalTime() * invsec); + const G4ThreeVector& p = track->GetMomentum(); + momentumAtBoundary_ = + math::XYZTLorentzVectorF(p.x() * invgev, p.y() * invgev, p.z() * invgev, track->GetTotalEnergy() * invgev); +} void TrackInformation::Print() const { LogDebug("TrackInformation") << " TrackInformation : storeTrack = " << storeTrack_ << "\n"