simcore::bertini::LDMXCascadeInterface Class Reference

G4CascadeInterface that captures cascade history via preprocessor hack. More...

#include <LDMXCascadeInterface.h>

Public Member Functions
	LDMXCascadeInterface (const G4String &name="LDMXBertiniCascade")
G4HadFinalState *	ApplyYourself (const G4HadProjectile &projectile, G4Nucleus &targetNucleus) override
void	setEnergyThreshold (double threshold)
double	getEnergyThreshold () const
const ldmx::CascadeHistory *	getLastCascadeHistory () const
	Returns nullptr if no history was captured.
ldmx::CascadeHistory	extractHistory ()
	Move the captured history out.
bool	hasHistory () const
void	setIncidentTrackId (int trackId)

Private Member Functions
void	ensureCascadeHistoryExists ()
	Ensure the G4CascadeHistory object exists in the cascader Geant4 only creates this if G4CASCADE_SHOW_HISTORY envvar is set, so we force-create it here to enable history capture.
void	captureHistory ()
	Extract history from the internal G4CascadeHistory Navigates: this->collider->theIntraNucleiCascader->theCascadeHistory.
void	captureDeexcitationProducts (G4HadFinalState *finalState)
	Capture de-excitation products from G4HadFinalState.

Private Attributes
double	energy_threshold_ {5000.0}
	Minimum photon energy threshold for recording [MeV].
int	incident_track_id_ {-1}
	Track ID of incident particle.
ldmx::CascadeHistory	last_history_
	Captured history from last cascade.

Detailed Description

G4CascadeInterface that captures cascade history via preprocessor hack.

After ApplyYourself(), call extractHistory() to get the cascade in LDMX format.

Definition at line 28 of file LDMXCascadeInterface.h.

Constructor & Destructor Documentation

LDMXCascadeInterface()

simcore::bertini::LDMXCascadeInterface::LDMXCascadeInterface

(

const G4String &

name = "LDMXBertiniCascade"

)

Definition at line 203 of file LDMXCascadeInterface.cxx.

    : G4CascadeInterface(name), incident_track_id_{-1} {}

Member Function Documentation

ApplyYourself()

G4HadFinalState * simcore::bertini::LDMXCascadeInterface::ApplyYourself ( const G4HadProjectile & projectile, G4Nucleus & targetNucleus )

override

Definition at line 208 of file LDMXCascadeInterface.cxx.

                                                                 {
  double photon_energy = projectile.GetTotalEnergy();
 
  ldmx_log(debug) << "LDMXCascadeInterface::ApplyYourself";
  ldmx_log(debug) << "  Track ID: " << incident_track_id_;
  ldmx_log(debug) << "  Photon energy: " << photon_energy << " MeV";
  ldmx_log(debug) << "  Energy threshold: " << energy_threshold_ << " MeV";
 
  last_history_.clear();
 
  bool should_record = (photon_energy >= energy_threshold_);
 
  if (should_record) {
    ldmx_log(debug) << "  Recording history (above threshold)";
    // Geant4 only creates the history object if G4CASCADE_SHOW_HISTORY is set,
    // so we force-create it here
    ensureCascadeHistoryExists();
  } else {
    ldmx_log(debug) << "  Skipping history (below threshold)";
  }
 
  G4HadFinalState* result =
      G4CascadeInterface::ApplyYourself(projectile, targetNucleus);
 
  if (should_record) {
    captureHistory();
    last_history_.setIncidentEnergy(photon_energy);
 
    ldmx_log(debug) << "  Captured " << last_history_.getSteps().size()
                    << " steps";
 
    // De-excitation (evaporation, gamma) happens after the cascade via
    // G4ExcitationHandler. These products appear in the final state but
    // are not in G4CascadeHistory, so we capture them separately.
    captureDeexcitationProducts(result);
 
    ldmx_log(debug) << "  Total steps: " << last_history_.getSteps().size();
 
    if (!last_history_.empty()) {
      CascadeHistoryStore::getInstance().addHistory(incident_track_id_,
                                                    last_history_);
    }
  }
 
  return result;
}

captureDeexcitationProducts()

void simcore::bertini::LDMXCascadeInterface::captureDeexcitationProducts ( G4HadFinalState * finalState )

private

Capture de-excitation products from G4HadFinalState.

De-excitation (evaporation, gamma, fission) is handled by G4ExcitationHandler after the cascade and isn't recorded in G4CascadeHistory. Identifies them by comparing final state secondaries against cascade escapees.

Definition at line 439 of file LDMXCascadeInterface.cxx.

                                 {
  if (!finalState) return;
 
  int n_secondaries = finalState->GetNumberOfSecondaries();
  if (n_secondaries == 0) return;
 
  // Match final state particles against cascade escapees to identify
  // de-excitation products. Energy matching has 1 MeV tolerance for recoil.
  struct EscapedParticle {
    int pdg_;
    double energy_;
    bool matched_;
  };
  std::vector<EscapedParticle> cascade_escaped;
 
  for (const auto& step : last_history_.getSteps()) {
    if (step.didEscape()) {
      cascade_escaped.push_back({step.getPdgId(), step.getEnergy(), false});
    }
  }
 
  int next_history_id = 0;
  for (const auto& step : last_history_.getSteps()) {
    if (step.getHistoryId() >= next_history_id) {
      next_history_id = step.getHistoryId() + 1;
    }
  }
 
  for (int i = 0; i < n_secondaries; ++i) {
    G4HadSecondary* secondary = finalState->GetSecondary(i);
    if (!secondary) continue;
 
    const G4DynamicParticle* dyn_particle = secondary->GetParticle();
    if (!dyn_particle) continue;
 
    int pdg = dyn_particle->GetPDGcode();
    double energy = dyn_particle->GetTotalEnergy();  // MeV
 
    bool is_deexcitation = true;
    const double energy_tolerance = 1.0;
 
    for (auto& escaped : cascade_escaped) {
      if (!escaped.matched_ && escaped.pdg_ == pdg &&
          std::abs(escaped.energy_ - energy) < energy_tolerance) {
        escaped.matched_ = true;
        is_deexcitation = false;
        break;
      }
    }
 
    if (is_deexcitation) {
      ldmx::CascadeStep step;
      step.setHistoryId(next_history_id++);
      step.setParentId(-2);  // marker for de-excitation origin
      step.setPdgId(pdg);
 
      G4ThreeVector mom = dyn_particle->GetMomentum();
      step.setMomentum(mom.x(), mom.y(), mom.z(), energy);
      step.setPosition(0, 0, 0);
      step.setGeneration(-1);  // marker for de-excitation
      step.setZone(0);
      step.setInteracted(false);
      step.setEscaped(true);
      step.setStage(ldmx::CascadeStage::DEEXCITATION);
 
      last_history_.addStep(std::move(step));
    }
  }
}

captureHistory()

void simcore::bertini::LDMXCascadeInterface::captureHistory ( )

private

Extract history from the internal G4CascadeHistory Navigates: this->collider->theIntraNucleiCascader->theCascadeHistory.

Definition at line 280 of file LDMXCascadeInterface.cxx.

                                          {
  // Navigate: this->collider->theIntraNucleiCascader->theCascadeHistory
  if (!collider) return;
 
  G4IntraNucleiCascader* cascader = collider->theIntraNucleiCascader;
  if (!cascader) return;
 
  G4CascadeHistory* g4_history = cascader->theCascadeHistory;
  if (!g4_history) return;
 
  last_history_.setIncidentTrackId(incident_track_id_);
 
  if (cascader->tnuclei) {
    last_history_.setTargetNucleus(cascader->tnuclei->getA(),
                                   cascader->tnuclei->getZ());
  }
 
  const std::vector<G4CascadeHistory::HistoryEntry>& entries =
      g4_history->theHistory;
 
  if (entries.empty()) {
    return;
  }
 
  // Build parent ID map from daughter lists
  std::vector<int> parent_ids(entries.size(), -1);
 
  for (size_t i = 0; i < entries.size(); ++i) {
    const auto& entry = entries[i];
    for (int d = 0; d < entry.n && d < 10; ++d) {
      int daughter_id = entry.dId[d];
      if (daughter_id >= 0 &&
          static_cast<size_t>(daughter_id) < entries.size()) {
        parent_ids[daughter_id] = static_cast<int>(i);
      }
    }
  }
 
  // First pass: convert G4 entries to LDMX steps
  std::vector<ldmx::CascadeStep> steps;
  steps.reserve(entries.size());
 
  for (size_t i = 0; i < entries.size(); ++i) {
    const auto& entry = entries[i];
    int parent_id = parent_ids[i];
 
    ldmx::CascadeStep step;
 
    const G4CascadParticle& cpart = entry.cpart;
    const G4InuclElementaryParticle& particle = cpart.getParticle();
 
    step.setHistoryId(cpart.getHistoryId());
    step.setParentId(parent_id);
    step.setPdgId(getPdgCode(particle.type()));
 
    // Geant4 Bertini uses GeV for momentum, fm for position
    G4LorentzVector mom = cpart.getMomentum();
    step.setMomentum(mom.px() * 1000.0, mom.py() * 1000.0, mom.pz() * 1000.0,
                     mom.e() * 1000.0);
 
    const G4ThreeVector& pos = cpart.getPosition();
    step.setPosition(pos.x(), pos.y(), pos.z());
 
    step.setGeneration(cpart.getGeneration());
    step.setZone(cpart.getCurrentZone());
 
    std::vector<int> daughter_ids;
    for (int d = 0; d < entry.n && d < 10; ++d) {
      daughter_ids.push_back(entry.dId[d]);
    }
    step.setDaughterIds(daughter_ids);
 
    bool interacted = (entry.n > 0);
    step.setInteracted(interacted);
 
    // Escaped if didn't interact (simplification - true escape needs more
    // tracking)
    bool escaped = !interacted && cpart.getGeneration() >= 0;
    step.setEscaped(escaped);
    step.setTargetPdgId(0);  // inferred in second pass
 
    ldmx::CascadeStage stage = ldmx::CascadeStage::UNKNOWN;
    int generation = cpart.getGeneration();
 
    if (generation == 0) {
      stage = ldmx::CascadeStage::INCIDENT;
    } else if (generation == 1) {
      stage = ldmx::CascadeStage::PRIMARY;
    } else if (!interacted && !escaped) {
      stage = ldmx::CascadeStage::ABSORBED;
    } else if (escaped && generation <= 2) {
      stage = ldmx::CascadeStage::PREEQUILIBRIUM;
    } else {
      stage = ldmx::CascadeStage::CASCADE;
    }
    step.setStage(stage);
 
    steps.push_back(std::move(step));
  }
 
  // Second pass: infer target nucleon from charge/baryon conservation
  for (auto& step : steps) {
    if (!step.didInteract()) continue;
 
    int bullet_charge = getCharge(step.getPdgId());
    int bullet_baryon = getBaryonNumber(step.getPdgId());
    int daughter_charge = 0;
    int daughter_baryon = 0;
 
    for (int daughter_id : step.getDaughterIds()) {
      for (const auto& s : steps) {
        if (s.getHistoryId() == daughter_id) {
          daughter_charge += getCharge(s.getPdgId());
          daughter_baryon += getBaryonNumber(s.getPdgId());
          break;
        }
      }
    }
 
    int delta_charge = daughter_charge - bullet_charge;
    int delta_baryon = daughter_baryon - bullet_baryon;
    step.setTargetPdgId(inferTargetPdg(delta_charge, delta_baryon));
  }
 
  last_history_.reserve(steps.size());
  for (auto& step : steps) {
    last_history_.addStep(std::move(step));
  }
 
  // Excitation energy approximation: E_incident - sum(KE of escaped)
  double total_escaped_energy = 0.0;
  int escaped_protons = 0;
  int escaped_neutrons = 0;
 
  for (const auto& step : last_history_.getSteps()) {
    if (step.didEscape()) {
      total_escaped_energy += step.getKineticEnergy();
      int pdg = step.getPdgId();
      if (pdg == 2212)
        escaped_protons++;
      else if (pdg == 2112)
        escaped_neutrons++;
    }
  }
 
  double excitation_energy =
      last_history_.getIncidentEnergy() - total_escaped_energy;
  if (excitation_energy < 0) excitation_energy = 0;
  last_history_.setExcitationEnergy(excitation_energy);
 
  int target_a = last_history_.getTargetA();
  int target_z = last_history_.getTargetZ();
  int residual_a = target_a - escaped_protons - escaped_neutrons;
  int residual_z = target_z - escaped_protons;
  if (residual_a < 0) residual_a = 0;
  if (residual_z < 0) residual_z = 0;
  last_history_.setResidualNucleus(residual_a, residual_z);
}

ensureCascadeHistoryExists()

void simcore::bertini::LDMXCascadeInterface::ensureCascadeHistoryExists ( )

private

Ensure the G4CascadeHistory object exists in the cascader Geant4 only creates this if G4CASCADE_SHOW_HISTORY envvar is set, so we force-create it here to enable history capture.

Definition at line 256 of file LDMXCascadeInterface.cxx.

                                                      {
  // Force-create the cascade history object if it doesn't exist
  // Geant4 normally only creates this if G4CASCADE_SHOW_HISTORY envvar is set
 
  if (!collider) {
    ldmx_log(debug) << "  ensureCascadeHistoryExists: no collider yet";
    return;
  }
 
  G4IntraNucleiCascader* cascader = collider->theIntraNucleiCascader;
  if (!cascader) {
    ldmx_log(debug) << "  ensureCascadeHistoryExists: no cascader yet";
    return;
  }
 
  if (!cascader->theCascadeHistory) {
    ldmx_log(debug)
        << "  ensureCascadeHistoryExists: creating G4CascadeHistory";
    cascader->theCascadeHistory = new G4CascadeHistory;
  } else {
    ldmx_log(debug) << "  ensureCascadeHistoryExists: history already exists";
  }
}

extractHistory()

ldmx::CascadeHistory simcore::bertini::LDMXCascadeInterface::extractHistory ( )

inline

Move the captured history out.

Definition at line 45 of file LDMXCascadeInterface.h.

{ return std::move(last_history_); }

References last_history_.

getEnergyThreshold()

double simcore::bertini::LDMXCascadeInterface::getEnergyThreshold ( ) const

inline

Definition at line 37 of file LDMXCascadeInterface.h.

{ return energy_threshold_; }

getLastCascadeHistory()

const ldmx::CascadeHistory * simcore::bertini::LDMXCascadeInterface::getLastCascadeHistory ( ) const

inline

Returns nullptr if no history was captured.

Definition at line 40 of file LDMXCascadeInterface.h.

                                                          {
    return last_history_.empty() ? nullptr : &last_history_;
  }

References last_history_.

hasHistory()

bool simcore::bertini::LDMXCascadeInterface::hasHistory ( ) const

inline

Definition at line 47 of file LDMXCascadeInterface.h.

{ return !last_history_.empty(); }

setEnergyThreshold()

void simcore::bertini::LDMXCascadeInterface::setEnergyThreshold ( double threshold )

inline

Definition at line 36 of file LDMXCascadeInterface.h.

{ energy_threshold_ = threshold; }

setIncidentTrackId()

void simcore::bertini::LDMXCascadeInterface::setIncidentTrackId ( int trackId )

inline

Definition at line 49 of file LDMXCascadeInterface.h.

{ incident_track_id_ = trackId; }

Member Data Documentation

energy_threshold_

double simcore::bertini::LDMXCascadeInterface::energy_threshold_ {5000.0}

private

Minimum photon energy threshold for recording [MeV].

Definition at line 76 of file LDMXCascadeInterface.h.

{5000.0};

incident_track_id_

int simcore::bertini::LDMXCascadeInterface::incident_track_id_ {-1}

private

Track ID of incident particle.

Definition at line 79 of file LDMXCascadeInterface.h.

{-1};

last_history_

ldmx::CascadeHistory simcore::bertini::LDMXCascadeInterface::last_history_

private

Captured history from last cascade.

Definition at line 82 of file LDMXCascadeInterface.h.

Referenced by extractHistory(), and getLastCascadeHistory().

---

The documentation for this class was generated from the following files:

• SimCore/include/SimCore/Bertini/LDMXCascadeInterface.h
• SimCore/src/SimCore/Bertini/LDMXCascadeInterface.cxx