NeutronSubstitutionPhysics.cxx
Go to the documentation of this file.
1 // Replaces neutrons with their products as found in the mcnp_libs library
2 //
3 // Chris Backhouse - c.backhouse@ucl.ac.uk
4 
6 
7 #include "Geant4/G4IonTable.hh"
8 #include "Geant4/G4ProcessManager.hh"
9 
10 #include "Geant4/G4Neutron.hh"
11 #include "Geant4/G4Proton.hh"
12 
14 #include "Geometry/Geometry.h"
15 
16 #include "TRandom3.h"
17 #include "TVector3.h"
18 
19 namespace g4n
20 {
24 
25  //-------------------------------------------------------------
27  : G4VPhysicsConstructor(name), fProc(0)
28  {
29  }
30 
31  //-------------------------------------------------------------
33  {
34  delete fProc;
35  }
36 
37  //-------------------------------------------------------------
39  {
40  // Make sure neutrons exist (probably not necessary?)
41  G4Neutron::NeutronDefinition();
42  }
43 
44  //-------------------------------------------------------------
46  {
47  if(fProc) return; // already constructed
48 
50 
51  // These final three numbers are ordering parameters whose significance I
52  // don't fully understand, but because we're defining a "decay" process,
53  // only the last (post-step) should be set.
54  G4Neutron::Definition()->GetProcessManager()->AddProcess(fProc, -1, -1, 0);
55  }
56 
57 
58  //-------------------------------------------------------------
60  NeutronSubstitutionProcess(double minKE, double maxKE,
61  const std::string& fname)
62  : G4VDiscreteProcess("neutronSubstitution", fDecay),
63  geom(art::ServiceHandle<geo::Geometry>().get()),
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  }
82 
83  //-------------------------------------------------------------
85  {
86  // Clean up the children of all the fates
87  delete[] fChildArena;
88  }
89 
90  //-------------------------------------------------------------
91  bool NeutronSubstitutionProcess::InsideDetector(const G4ThreeVector& r) const
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  }
98 
99  //-------------------------------------------------------------
101  G4double prevStep,
102  G4ForceCondition* cond)
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  }
119 
120  //-------------------------------------------------------------
121  G4VParticleChange* NeutronSubstitutionProcess::
122  PostStepDoIt(const G4Track& track,
123  const G4Step& step)
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  }
223 }
224 
225 #include "nug4/G4Base/G4PhysicsProcessFactorySingleton.hh"
XYZ dr
Position of the child relative to parent start (mm)
Definition: NeutronFates.h:20
TruthSlim – remove generated objects that don&#39;t contribute.
const XML_Char * name
Definition: expat.h:151
NeutronSubstitutionPhysics(const G4String &name="NeutronSubstitution")
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int pdg
PDG code of child particle.
Definition: NeutronFates.h:23
T sqrt(T number)
Definition: d0nt_math.hpp:156
Definition: event.h:19
double DetLength() const
const Child * children
Definition: NeutronFates.h:30
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &step) override
virtual double GetMeanFreePath(const G4Track &trk, G4double prevStep, G4ForceCondition *condition) override
float p
Initial neutron momentum (MeV)
Definition: NeutronFates.h:25
#define M_PI
Definition: SbMath.h:34
float abs(float number)
Definition: d0nt_math.hpp:39
const geo::Geometry * geom
cache service handle
Track finder for cosmic rays.
void LoadNeutronFates(const std::string &name, std::vector< NeutronFate > &out, NeutronFate::Child *&childArena)
Fills out with list of neutrons from lowest to highest energy.
Registers NeutronSubstitutionProcess.
Transforms neutrons into their final states according to a library.
NeutronSubstitutionProcess(double minKE, double maxKE, const std::string &fname)
float dt
Time of the child relative to parent start (ns)
Definition: NeutronFates.h:21
const double a
std::string getenv(std::string const &name)
bool InsideDetector(const G4ThreeVector &r) const
double mag2() const
double DetHalfHeight() 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 DetHalfWidth() const
void geom(int which=0)
Definition: geom.C:163
XYZ p
three-momentum of child particle (MeV)
Definition: NeutronFates.h:22
double x() const
const hit & b
Definition: hits.cxx:21
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
Helper for AttenCurve.
Definition: Path.h:10
Encapsulate the geometry of one entire detector (near, far, ndos)
enum BeamMode string