Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
g4n::NeutronSubstitutionProcess Class Reference

Transforms neutrons into their final states according to a library. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-28/g4nova/NeutronSubstitutionPhysics.h"

Inheritance diagram for g4n::NeutronSubstitutionProcess:

Public Member Functions

 NeutronSubstitutionProcess (double minKE, double maxKE, const std::string &fname)
 
 ~NeutronSubstitutionProcess ()
 
virtual double GetMeanFreePath (const G4Track &trk, G4double prevStep, G4ForceCondition *condition) override
 
virtual G4VParticleChange * PostStepDoIt (const G4Track &track, const G4Step &step) override
 

Protected Member Functions

bool InsideDetector (const G4ThreeVector &r) const
 

Protected Attributes

const geo::Geometrygeom
 cache service handle More...
 
double fMinKE
 
double fMaxKE
 In MeV. More...
 
std::vector< NeutronFatefFates
 
NeutronFate::ChildfChildArena
 

Detailed Description

Transforms neutrons into their final states according to a library.

Definition at line 38 of file NeutronSubstitutionPhysics.h.

Constructor & Destructor Documentation

g4n::NeutronSubstitutionProcess::NeutronSubstitutionProcess ( double  minKE,
double  maxKE,
const std::string fname 
)

Definition at line 60 of file NeutronSubstitutionPhysics.cxx.

References om::cout, allTimeWatchdog::endl, fChildArena, fFates, plot_validation_datamc::fname, cet::getenv(), g4n::LoadNeutronFates(), and string.

Referenced by g4n::NeutronSubstitutionPhysics::ConstructProcess().

62  : G4VDiscreteProcess("neutronSubstitution", fDecay),
64  fMinKE(minKE), fMaxKE(maxKE)
65  {
66  // Apparently geant insists the subtype is set to something. I haven't
67  // investigated what the different values mean.
68  SetProcessSubType(1);
69 
70  char* libPath = getenv("MCNP_LIBS_LIB_PATH");
71  if(!libPath){
72  std::cout << "Environment variable $MCNP_LIBS_LIB_PATH is not set" << std::endl;
73  abort();
74  }
75 
76  std::cout << "NeutronSubstitution range: " << minKE << " MeV < KE < " << maxKE << " MeV" << std::endl;
77  const std::string fatesFile = std::string(libPath)+"/"+fname;
78  std::cout << "NeutronSubstitution: Loading neutron fates from " << fatesFile << " ..." << std::endl;
79  LoadNeutronFates(fatesFile, fFates, fChildArena);
80  std::cout << "Loaded " << fFates.size() << " neutron fates" << std::endl;
81  }
const geo::Geometry * geom
cache service handle
void LoadNeutronFates(const std::string &name, std::vector< NeutronFate > &out, NeutronFate::Child *&childArena)
Fills out with list of neutrons from lowest to highest energy.
std::string getenv(std::string const &name)
OStream cout
Definition: OStream.cxx:6
enum BeamMode string
g4n::NeutronSubstitutionProcess::~NeutronSubstitutionProcess ( )

Definition at line 84 of file NeutronSubstitutionPhysics.cxx.

References fChildArena.

85  {
86  // Clean up the children of all the fates
87  delete[] fChildArena;
88  }

Member Function Documentation

double g4n::NeutronSubstitutionProcess::GetMeanFreePath ( const G4Track &  trk,
G4double  prevStep,
G4ForceCondition *  condition 
)
overridevirtual

Definition at line 100 of file NeutronSubstitutionPhysics.cxx.

References fMaxKE, fMinKE, InsideDetector(), and PostStepDoIt().

103  {
104  // Leave neutrons outside the configured energy window alone - ie geant
105  // will handle them
106  if((fMinKE > 0 && trk.GetKineticEnergy() < fMinKE) ||
107  (fMaxKE > 0 && trk.GetKineticEnergy() > fMaxKE)){
108  return std::numeric_limits<double>::infinity();
109  }
110 
111  // Also don't do anything to neutrons outside the detector
112  if(!InsideDetector(trk.GetPosition())){
113  return std::numeric_limits<double>::infinity();
114  }
115 
116  // Happens immediately
117  return 0;
118  }
Track finder for cosmic rays.
bool InsideDetector(const G4ThreeVector &r) const
bool g4n::NeutronSubstitutionProcess::InsideDetector ( const G4ThreeVector &  r) const
protected

Definition at line 91 of file NeutronSubstitutionPhysics.cxx.

References std::abs(), geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), and geom.

Referenced by GetMeanFreePath(), and PostStepDoIt().

92  {
93  // Factor of 10 converts novasoft's cms to geant's mms
94  return (std::abs(r.X) < 10*geom->DetHalfWidth() &&
95  std::abs(r.Y) < 10*geom->DetHalfHeight() &&
96  r.Z > 0 && r.Z < 10*geom->DetLength());
97  }
double DetLength() const
float abs(float number)
Definition: d0nt_math.hpp:39
const geo::Geometry * geom
cache service handle
double DetHalfHeight() const
double DetHalfWidth() const
TRandom3 r(0)
G4VParticleChange * g4n::NeutronSubstitutionProcess::PostStepDoIt ( const G4Track &  track,
const G4Step &  step 
)
overridevirtual

Definition at line 122 of file NeutronSubstitutionPhysics.cxx.

References a, ana::assert(), b, g4n::NeutronFate::children, std::cos(), om::cout, CLHEP::Hep3Vector::cross(), g4n::NeutronFate::Child::dr, g4n::NeutronFate::Child::dt, allTimeWatchdog::endl, stan::math::fabs(), fFates, InsideDetector(), it, m, M_PI, CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::mag2(), g4n::NeutronFate::nchildren, g4n::NeutronFate::Child::p, g4n::NeutronFate::p, plot_validation_datamc::p1, plot_validation_datamc::p2, g4n::NeutronFate::Child::pdg, std::sin(), std::sqrt(), CLHEP::Hep3Vector::unit(), g4n::NeutronFate::XYZ::x, CLHEP::Hep3Vector::x(), g4n::NeutronFate::XYZ::y, and g4n::NeutronFate::XYZ::z.

Referenced by GetMeanFreePath().

124  {
125  // We're only supposed to be invoked for neutron tracks
126  assert(track.GetDefinition()->GetPDGEncoding() == 2112);
127 
128  G4VParticleChange* change = new G4VParticleChange;
129  change->Initialize(track);
130  // Stop the existing track
131  change->ProposeTrackStatus(fStopAndKill);
132 
133  const CLHEP::Hep3Vector neutMom = track.GetDynamicParticle()->Get4Momentum().vect();
134 
135  static auto prev_it = fFates.end();
136 
137  // Find the library entry closest to the input neutron's energy
138  auto it = std::lower_bound(fFates.begin(), fFates.end(), neutMom.mag(),
139  [](const NeutronFate& a, double b)
140  {
141  return a.p < b;
142  });
143  if(it == fFates.end()){
144  std::cout << "NeutronSubstitutionProcess: WARNING! no library match for neutron with p = " << neutMom.mag() << " MeV, killing." << std::endl;
145  prev_it = fFates.end();
146  return change;
147  }
148 
149  // If we're about to throw the same exact fate as we already used last
150  // time, go to the next lowest in energy to avoid an infinite loop.
151  if(it == prev_it) --it;
152  prev_it = it;
153 
154  const NeutronFate& fate = *it;
155 
156  if(fate.nchildren > 0){
157  // Setup the coordinate system the children will be produced in. Z is
158  // defined to be parallel to the initial momentum.
159  const CLHEP::Hep3Vector zUnit = neutMom.unit();
160  // Construct two perpendicular unit vectors
161  const CLHEP::Hep3Vector p1 = (fabs(zUnit.x()) < 0.5 ? CLHEP::Hep3Vector(1, 0, 0) : CLHEP::Hep3Vector(0, 1, 0)).cross(zUnit).unit();
162  const CLHEP::Hep3Vector p2 = p1.cross(zUnit);
163  // And mix them to form a random xy basis
164  const double ang = gRandom->Uniform(0, 2*M_PI);
165  const CLHEP::Hep3Vector xUnit = cos(ang)*p1 + sin(ang)*p2;
166  const CLHEP::Hep3Vector yUnit = -sin(ang)*p1 + cos(ang)*p2;
167 
168  for(int childIdx = 0; childIdx < fate.nchildren; ++childIdx){
169  const NeutronFate::Child& child = fate.children[childIdx];
170 
171  const G4ParticleDefinition* partdef;
172  if(child.pdg < 1000000000){
173  partdef = G4ParticleTable::GetParticleTable()->FindParticle(child.pdg);
174  }
175  else if(child.pdg < 2000000000){
176  partdef = G4IonTable::GetIonTable()->GetIon(child.pdg);
177  }
178  else{
179  // special GENIE particle? skip
180  //
181  // n+n | 2000000200
182  // n+p | 2000000201
183  // NucBindE | 2000000101
184  // HadrClus | 2000000300
185  // HadrSyst | 2000000001
186  // HadrBlob | 2000000002
187  continue;
188  }
189 
190  if(!partdef){
191  std::cout << "NeutronSubstitutionProcess: Warning: skipping unknown daughter particle with pdg = " << child.pdg << ", p = " << sqrt(child.p.x*child.p.x + child.p.y*child.p.y + child.p.z*child.p.z) << " MeV" << std::endl;
192  continue;
193  }
194 
195  // Transform the child particle position into the frame defined above
196  const CLHEP::Hep3Vector childr0 = step.GetPreStepPoint()->GetPosition() + child.dr.x*xUnit + child.dr.y*yUnit + child.dr.z*zUnit;
197 
198  // Don't create children outside the detector. This isn't perfect,
199  // since we still allow cases where the neutron would leave the
200  // detector, "scatter" outside, come back in and create a child, but
201  // it's the best we can reasonably do.
202  if(!InsideDetector(childr0)) continue;
203 
204  const double childt0 = step.GetPreStepPoint()->GetGlobalTime() + child.dt;
205 
206  // Transform the particle momentum into the same coordinate system and
207  // compute the energy.
208  const CLHEP::Hep3Vector childp = xUnit*child.p.x + yUnit*child.p.y + zUnit*child.p.z;
209 
210  const double m = partdef->GetPDGMass();
211  const double childE = sqrt(m*m + childp.mag2());
212  const CLHEP::HepLorentzVector childp4(childp, childE);
213 
214  // Create the new particle/track as a daughter of the neutron we killed
215  G4DynamicParticle* childpart = new G4DynamicParticle(partdef, childp4);
216  G4Track* childtrack = new G4Track(childpart, childt0, childr0);
217  change->AddSecondary(childtrack);
218  } // end for children
219  } // end if any children
220 
221  return change;
222  }
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
Definition: event.h:19
double z() const
#define M_PI
Definition: SbMath.h:34
const double a
bool InsideDetector(const G4ThreeVector &r) const
double mag2() const
Hep3Vector cross(const Hep3Vector &) const
Hep3Vector unit() const
OStream cout
Definition: OStream.cxx:6
double mag() const
T sin(T number)
Definition: d0nt_math.hpp:132
double x() const
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
double y() const

Member Data Documentation

NeutronFate::Child* g4n::NeutronSubstitutionProcess::fChildArena
protected
std::vector<NeutronFate> g4n::NeutronSubstitutionProcess::fFates
protected

Definition at line 58 of file NeutronSubstitutionPhysics.h.

Referenced by NeutronSubstitutionProcess(), and PostStepDoIt().

double g4n::NeutronSubstitutionProcess::fMaxKE
protected

In MeV.

Definition at line 56 of file NeutronSubstitutionPhysics.h.

Referenced by GetMeanFreePath().

double g4n::NeutronSubstitutionProcess::fMinKE
protected

Definition at line 56 of file NeutronSubstitutionPhysics.h.

Referenced by GetMeanFreePath().

const geo::Geometry* g4n::NeutronSubstitutionProcess::geom
protected

cache service handle

Definition at line 54 of file NeutronSubstitutionPhysics.h.

Referenced by InsideDetector().


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