TargetBremFilter.cxx

/*~~~~~~~~~~~~~*/
/*   Biasing   */
/*~~~~~~~~~~~~~*/
#include "Biasing/TargetBremFilter.h"
 
/*~~~~~~~~~~~~*/
/*   Geant4   */
/*~~~~~~~~~~~~*/
#include "G4Electron.hh"
#include "G4EventManager.hh"
#include "G4RunManager.hh"
 
/*~~~~~~~~~~~~~*/
/*   SimCore   */
/*~~~~~~~~~~~~~*/
#include "SimCore/G4User/PtrRetrieval.h"
#include "SimCore/G4User/UserEventInformation.h"
#include "SimCore/G4User/UserTrackInformation.h"
 
/*~~~~~~~~~~~~~*/
/*  C++ StdLib */
/*~~~~~~~~~~~~~*/
#include <cmath>
 
namespace biasing {
 
TargetBremFilter::TargetBremFilter(const std::string& name,
                                   framework::config::Parameters& parameters)
    : simcore::UserAction(name, parameters) {
  recoil_max_p_threshold_ = parameters.get<double>("recoil_max_p_threshold");
  brem_energy_threshold_ = parameters.get<double>("brem_min_energy_threshold");
  brem_theta_min_ = parameters.get<double>("brem_theta_min");
  brem_theta_max_ = parameters.get<double>("brem_theta_max");
  dral_min_ = parameters.get<double>("dral_min");
  dral_max_ = parameters.get<double>("dral_max");
  kill_recoil_ = parameters.get<bool>("kill_recoil_track");
 
  ldmx_log(info) << "TargetBremFilter configured:"
                 << " recoil_max_p_threshold=" << recoil_max_p_threshold_
                 << " brem_min_energy_threshold=" << brem_energy_threshold_
                 << " brem_theta=[" << brem_theta_min_ << ", "
                 << brem_theta_max_ << "]"
                 << " dral=[" << dral_min_ << ", " << dral_max_ << "]";
}
 
TargetBremFilter::~TargetBremFilter() {}
 
G4ClassificationOfNewTrack TargetBremFilter::classifyNewTrack(
    const G4Track* track, const G4ClassificationOfNewTrack& currentTrackClass) {
  // get the PDGID of the track.
  G4int pdg_id = track->GetParticleDefinition()->GetPDGEncoding();
 
  // Get the particle type.
  G4String particle_name = track->GetParticleDefinition()->GetParticleName();
 
  // Use current classification by default so values from other plugins are not
  // overridden.
  G4ClassificationOfNewTrack classification = currentTrackClass;
 
  if (track->GetTrackID() == 1 && pdg_id == 11) {
    return fWaiting;
  }
 
  return classification;
}
 
void TargetBremFilter::stepping(const G4Step* step) {
  double last_dral = -1;
  double last_theta = -1;
  // Get the track associated with this step.
  auto track{step->GetTrack()};
 
  // Only process the primary electron track
  if (track->GetParentID() != 0) return;
 
  // Get the PDG ID of the track and make sure it's an electron. If
  // another particle type is found, thrown an exception.
  if (auto pdg_id{track->GetParticleDefinition()->GetPDGEncoding()};
      pdg_id != 11)
    return;
 
  // Get the region the particle is currently in.  Continue processing
  // the particle only if it's in the target region.
  static auto target_region =
      simcore::g4user::ptrretrieval::getRegion("target");
  if (!target_region) {
    ldmx_log(warn) << "Region 'target' not found in Geant4 region store";
  }
  auto phy_vol{track->GetVolume()};
  auto log_vol{phy_vol ? phy_vol->GetLogicalVolume() : nullptr};
  auto track_region{log_vol ? log_vol->GetRegion() : nullptr};
  if (track_region != target_region) return;
 
  /*
  std::cout << "[TargetBremFilter] : Stepping primary electron in 'target'
  region into "
    << track->GetNextVolume()->GetName()
    << std::endl;
   */
 
  auto recoil_physical_volume =
      simcore::g4user::ptrretrieval::getPhysicalVolume("recoil_PV");
  auto world_physical_volume =
      simcore::g4user::ptrretrieval::getPhysicalVolume("World_PV");
  if (!recoil_physical_volume) {
    ldmx_log(warn) << "Volume 'recoil_PV' not found in Geant4 volume store";
  }
  if (!world_physical_volume) {
    ldmx_log(warn) << "Volume 'World_PV' not found in Geant4 volume store";
  }
  auto track_volume = track->GetNextVolume();
  if (track_volume == recoil_physical_volume or
      track_volume == world_physical_volume) {
    // If the recoil electron momentum is too high, abort
    if (track->GetMomentum().mag() >= recoil_max_p_threshold_) {
      ldmx_log(trace) << "Abort: recoil p=" << track->GetMomentum().mag()
                      << " MeV >= threshold=" << recoil_max_p_threshold_;
      track->SetTrackStatus(fKillTrackAndSecondaries);
      G4RunManager::GetRunManager()->AbortEvent();
      return;
    }
 
    // Get the electron secondaries
    bool has_brem_candidate = false;
    if (auto secondaries = step->GetSecondary(); secondaries->size() == 0) {
      ldmx_log(trace) << "Abort: no secondaries produced";
      track->SetTrackStatus(fKillTrackAndSecondaries);
      G4RunManager::GetRunManager()->AbortEvent();
      return;
    } else {
      ldmx_log(trace) << "Exiting target: recoil p ="
                      << track->GetMomentum().mag() << " MeV, "
                      << secondaries->size() << " secondaries";
 
      for (auto& secondary_track : *secondaries) {
        auto electron = G4Electron::Definition();
        auto ebrem_process =
            simcore::g4user::ptrretrieval::getProcess(electron, "eBrem");
        if (!ebrem_process) {
          ldmx_log(warn) << "Process 'eBrem' not found in Geant4 process store";
        }
 
        if (ebrem_process &&
            secondary_track->GetKineticEnergy() > brem_energy_threshold_) {
          // Check if secondary is photon
          auto secondary_pdg_id =
              secondary_track->GetParticleDefinition()->GetPDGEncoding();
          if (secondary_pdg_id != 22) continue;
 
          // Brem angle
          auto momentum = secondary_track->GetMomentum();
          double theta = std::atan2(std::sqrt(momentum.x() * momentum.x() +
                                              momentum.y() * momentum.y()),
                                    momentum.z());
          bool pass_brem_theta =
              theta >= brem_theta_min_ && theta <= brem_theta_max_;
 
          // Maximum and Minimum angle between outgoing Brem photons and
          // electrons
          auto gamma_mom = secondary_track->GetMomentum();
          auto electron_mom = track->GetMomentum();
          double gamma_eta = gamma_mom.eta();
          double gamma_phi = gamma_mom.phi();
          double electron_eta = electron_mom.eta();
          double electron_phi = electron_mom.phi();
 
          double dphi = std::atan2(std::sin(electron_phi - gamma_phi),
                                   std::cos(electron_phi - gamma_phi));
          double dral = std::sqrt((electron_eta - gamma_eta) *
                                      (electron_eta - gamma_eta) +
                                  dphi * dphi);
          bool pass_dral = dral >= dral_min_ && dral <= dral_max_;
 
          last_theta = theta;
          last_dral = dral;
 
          ldmx_log(trace) << "  photon E="
                          << secondary_track->GetKineticEnergy()
                          << " MeV, theta=" << theta
                          << " (pass=" << pass_brem_theta << ")"
                          << ", dR=" << dral << " (pass=" << pass_dral << ")";
 
          if (pass_brem_theta && pass_dral) {
            auto track_info{
                simcore::UserTrackInformation::get(secondary_track)};
            track_info->tagBremCandidate();
 
            getEventInfo()->incBremCandidateCount();
 
            has_brem_candidate = true;
          }
        }
      }
    }
 
    if (!has_brem_candidate) {
      ldmx_log(trace) << "Abort: no brem candidate passed cuts";
      track->SetTrackStatus(fKillTrackAndSecondaries);
      G4RunManager::GetRunManager()->AbortEvent();
      return;
    }
 
    ldmx_log(info) << "ACCEPTED: brem candidate with theta=" << last_theta
                   << ", dR=" << last_dral;
 
    // Check if the recoil electron should be killed.  If not, postpone
    // its processing until the brem gamma has been processed.
    if (kill_recoil_)
      track->SetTrackStatus(fStopAndKill);
    else
      track->SetTrackStatus(fSuspend);
 
  } else if (step->GetPostStepPoint()->GetKineticEnergy() == 0) {
    track->SetTrackStatus(fKillTrackAndSecondaries);
    G4RunManager::GetRunManager()->AbortEvent();
    return;
  }
}
 
void TargetBremFilter::endOfEventAction(const G4Event*) {}
}  // namespace biasing
 
DECLARE_ACTION(biasing::TargetBremFilter)