Merge pull request #41181 from civanch/user_actions_and_mctruth

Geant4 user actions and CMS mctruth handling clean-up

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";
 #endif
       }
     }

SimG4Core/Application/interface/SteppingAction.h

@@ -14,7 +14,6 @@
 #include <string>
 #include <vector>
 
-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;

SimG4Core/Application/src/EventAction.cc

@@ -19,9 +19,7 @@ EventAction::EventAction(const edm::ParameterSet& p,
       m_SteppingVerbose(sv),
       m_stopFile(p.getParameter<std::string>("StopFile")),
       m_printRandom(p.getParameter<bool>("PrintRandomSeed")),
-      m_debug(p.getUntrackedParameter<bool>("debug", false)) {
-  m_trackManager->setCollapsePrimaryVertices(p.getParameter<bool>("CollapsePrimaryVertices"));
-}
+      m_debug(p.getUntrackedParameter<bool>("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(); }

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);

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<double>("CriticalEnergyForVacuum") * CLHEP::MeV);
   if (0.0 < theCriticalEnergyForVacuum) {
     killBeamPipe = true;
@@ -162,17 +162,10 @@ 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<math::XYZVectorD, math::XYZTLorentzVectorD> p(pos, mom);
-      trackManager_->addTkCaloStateInfo(id, p);
+      TrackInformation* trkinfo = static_cast<TrackInformation*>(theTrack->GetUserInformation());
+      if (!trkinfo->crossedBoundary()) {
+        trkinfo->setCrossedBoundary(theTrack);
+      }
     }
   } else {
     theTrack->SetTrackStatus(fStopAndKill);

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,72 +74,29 @@ 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<math::XYZVectorD, math::XYZTLorentzVectorD> 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);
-
-#ifdef EDM_ML_DEBUG
-    edm::LogVerbatim("SimTrackManager") << "TrackingAction addTrack " << id << "  "
-                                        << aTrack->GetDefinition()->GetParticleName() << " added= " << withAncestor
-                                        << " at " << aTrack->GetPosition();
-#endif
-
-  } else {
-    trackManager_->addTrack(currentTrack_, false, false);
-    delete currentTrack_;
+  trackManager_->addTrack(currentTrack_, aTrack, isInHistory, withAncestor);
 
 #ifdef EDM_ML_DEBUG
-    edm::LogVerbatim("SimTrackManager") << "TrackingAction addTrack " << id << " added with false flags and deleted";
+  edm::LogVerbatim("TrackingAction") << "TrackingAction end track=" << id << "  "
+                                     << aTrack->GetDefinition()->GetParticleName() << " ansestor= " << withAncestor
+                                     << " saved= " << currentTrack_->saved() << " end point " << aTrack->GetPosition();
 #endif
-  }
 
   EndOfTrack et(aTrack);
   m_endOfTrackSignal(&et);

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<float> param_;
   std::vector<G4SimTrack*> g4tracks_;
   std::vector<G4SimVertex*> g4vertices_;

SimG4Core/Notification/interface/G4SimTrack.h

@@ -3,38 +3,16 @@
 
 #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 <cmath>
 
 class G4SimTrack {
 public:
-  G4SimTrack() {}
-
   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),
-        parentMomentum_(math::XYZVectorD(0., 0., 0.)),
-        tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)),
-        tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)),
-        crossedBoundary_(false) {}
+      : 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),
-        tkSurfacePosition_(math::XYZVectorD(0., 0., 0.)),
-        tkSurfaceMomentum_(math::XYZTLorentzVectorD(0., 0., 0., 0.)),
-        crossedBoundary_(false) {}
+      : id_(iid), ipart_(ipart), ip_(ip), ie_(ie), ivert_(iv), igenpart_(ig), parentMomentum_(ipmom) {}
 
   G4SimTrack(int iid,
              int ipart,
@@ -53,8 +31,7 @@ class G4SimTrack {
         igenpart_(ig),
         parentMomentum_(ipmom),
         tkSurfacePosition_(tkpos),
-        tkSurfaceMomentum_(tkmom),
-        crossedBoundary_(false) {}
+        tkSurfaceMomentum_(tkmom) {}
 
   ~G4SimTrack() = default;
 
@@ -94,13 +71,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

SimG4Core/Notification/interface/GenParticleInfo.h

@@ -6,6 +6,7 @@
 class GenParticleInfo : public G4VUserPrimaryParticleInformation {
 public:
   explicit GenParticleInfo(int id) : id_(id) {}
+  ~GenParticleInfo() override = default;
   int id() const { return id_; }
   void Print() const override {}
 

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: