Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4HENeutronInelastic Class Reference

#include <G4HENeutronInelastic.hh>

+ Inheritance diagram for G4HENeutronInelastic:

Public Member Functions

 G4HENeutronInelastic ()
 
 ~G4HENeutronInelastic ()
 
virtual void ModelDescription (std::ostream &) const
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4int GetNumberOfSecondaries ()
 
void FirstIntInCasNeutron (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
 
- Public Member Functions inherited from G4HEInelastic
 G4HEInelastic (const G4String &modelName="HEInelastic")
 
 ~G4HEInelastic ()
 
void SetMaxNumberOfSecondaries (const G4int maxnumber)
 
void SetVerboseLevel (const G4int level)
 
void ForceEnergyConservation (G4bool energyConservation)
 
G4bool EnergyConservation (void)
 
virtual const std::pair< G4double, G4doubleGetFatalEnergyCheckLevels () const
 
G4double Amin (G4double a, G4double b)
 
G4double Amax (G4double a, G4double b)
 
G4int Imin (G4int a, G4int b)
 
G4int Imax (G4int a, G4int b)
 
void FillParticleChange (G4HEVector pv[], G4int aVecLength)
 
G4double pmltpc (G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
 
G4int Factorial (G4int n)
 
G4double NuclearInelasticity (G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber)
 
G4double NuclearExcitation (G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber, G4double &excitationEnergyCascade, G4double &excitationEnergyEvaporation)
 
void HighEnergyCascading (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void HighEnergyClusterProduction (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void TuningOfHighEnergyCascading (G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void MediumEnergyCascading (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void MediumEnergyClusterProduction (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void QuasiElasticScattering (G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
 
void ElasticScattering (G4bool &successful, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, G4double atomicWeight, G4double atomicNumber)
 
G4int rtmi (G4double *x, G4double xli, G4double xri, G4double eps, G4int iend, G4double aa, G4double bb, G4double cc, G4double dd, G4double rr)
 
G4double fctcos (G4double t, G4double aa, G4double bb, G4double cc, G4double dd, G4double rr)
 
void StrangeParticlePairProduction (const G4double availableEnergy, const G4double centerOfMassEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)
 
G4double NBodyPhaseSpace (const G4double totalEnergy, const G4bool constantCrossSection, G4HEVector pv[], G4int &vecLen)
 
G4double NBodyPhaseSpace (G4int npart, G4HEVector pv[], G4double wmax, G4double wfcn, G4int maxtrial, G4int ntrial)
 
G4double gpdk (G4double a, G4double b, G4double c)
 
void QuickSort (G4double arr[], const G4int lidx, const G4int ridx)
 
G4double Alam (G4double a, G4double b, G4double c)
 
G4double CalculatePhaseSpaceWeight (G4int npart)
 
G4double normal (void)
 
G4double GammaRand (G4double avalue)
 
G4double Erlang (G4int mvalue)
 
G4int Poisson (G4double x)
 
void SetParticles (void)
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
const G4HadronicInteractionGetMyPointer () const
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair< G4double, G4doubleGetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4doubleGetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 

Public Attributes

G4int vecLength
 
- Public Attributes inherited from G4HEInelastic
G4int verboseLevel
 
G4int MAXPART
 
G4bool conserveEnergy
 
G4HEVector PionPlus
 
G4HEVector PionZero
 
G4HEVector PionMinus
 
G4HEVector KaonPlus
 
G4HEVector KaonZero
 
G4HEVector AntiKaonZero
 
G4HEVector KaonMinus
 
G4HEVector KaonZeroShort
 
G4HEVector KaonZeroLong
 
G4HEVector Proton
 
G4HEVector AntiProton
 
G4HEVector Neutron
 
G4HEVector AntiNeutron
 
G4HEVector Lambda
 
G4HEVector AntiLambda
 
G4HEVector SigmaPlus
 
G4HEVector SigmaZero
 
G4HEVector SigmaMinus
 
G4HEVector AntiSigmaPlus
 
G4HEVector AntiSigmaZero
 
G4HEVector AntiSigmaMinus
 
G4HEVector XiZero
 
G4HEVector XiMinus
 
G4HEVector AntiXiZero
 
G4HEVector AntiXiMinus
 
G4HEVector OmegaMinus
 
G4HEVector AntiOmegaMinus
 
G4HEVector Deuteron
 
G4HEVector Triton
 
G4HEVector Alpha
 
G4HEVector Gamma
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 52 of file G4HENeutronInelastic.hh.

Constructor & Destructor Documentation

◆ G4HENeutronInelastic()

G4HENeutronInelastic::G4HENeutronInelastic ( )
inline

Definition at line 55 of file G4HENeutronInelastic.hh.

55 : G4HEInelastic("G4HENeutronInelastic")
56 {
57 vecLength = 0;
58 theMinEnergy = 45*CLHEP::GeV;
59 theMaxEnergy = 10*CLHEP::TeV;
60 MAXPART = 2048;
61 verboseLevel = 0;
62 G4cout << "WARNING: model G4HENeutronInelastic is being deprecated and will\n"
63 << "disappear in Geant4 version 10.0" << G4endl;
64 }
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout

◆ ~G4HENeutronInelastic()

G4HENeutronInelastic::~G4HENeutronInelastic ( )
inline

Definition at line 66 of file G4HENeutronInelastic.hh.

66{};

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4HENeutronInelastic::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
)
virtual

Implements G4HadronicInteraction.

Definition at line 57 of file G4HENeutronInelastic.cc.

59{
60 G4HEVector* pv = new G4HEVector[MAXPART];
61 const G4HadProjectile* aParticle = &aTrack;
62 const G4double A = targetNucleus.GetA_asInt();
63 const G4double Z = targetNucleus.GetZ_asInt();
64 G4HEVector incidentParticle(aParticle);
65
66 G4double atomicNumber = Z;
67 G4double atomicWeight = A;
68
69 G4int incidentCode = incidentParticle.getCode();
70 G4double incidentMass = incidentParticle.getMass();
71 G4double incidentTotalEnergy = incidentParticle.getEnergy();
72
73 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
74 // DHW 19 May 2011: variable set but not used
75
76 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
77
78 if (incidentKineticEnergy < 1.)
79 G4cout << "GHENeutronInelastic: incident energy < 1 GeV" << G4endl;
80
81 if (verboseLevel > 1) {
82 G4cout << "G4HENeutronInelastic::ApplyYourself" << G4endl;
83 G4cout << "incident particle " << incidentParticle.getName()
84 << "mass " << incidentMass
85 << "kinetic energy " << incidentKineticEnergy
86 << G4endl;
87 G4cout << "target material with (A,Z) = ("
88 << atomicWeight << "," << atomicNumber << ")" << G4endl;
89 }
90
91 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
92 atomicWeight, atomicNumber);
93 if (verboseLevel > 1)
94 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
95
96 incidentKineticEnergy -= inelasticity;
97
98 G4double excitationEnergyGNP = 0.;
99 G4double excitationEnergyDTA = 0.;
100
101 G4double excitation = NuclearExcitation(incidentKineticEnergy,
102 atomicWeight, atomicNumber,
103 excitationEnergyGNP,
104 excitationEnergyDTA);
105 if (verboseLevel > 1)
106 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
107 << excitationEnergyDTA << G4endl;
108
109 incidentKineticEnergy -= excitation;
110 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
111 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
112 // *(incidentTotalEnergy+incidentMass));
113 // DHW 19 May 2011: variables et but not used
114
115 G4HEVector targetParticle;
116 if (G4UniformRand() < atomicNumber/atomicWeight) {
117 targetParticle.setDefinition("Proton");
118 } else {
119 targetParticle.setDefinition("Neutron");
120 }
121
122 G4double targetMass = targetParticle.getMass();
123 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
124 + targetMass*targetMass
125 + 2.0*targetMass*incidentTotalEnergy);
126 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
127
128 // In the original Gheisha code the inElastic flag was set as follows
129 // G4bool inElastic = InElasticCrossSectionInFirstInt
130 // (availableEnergy, incidentCode, incidentTotalMomentum);
131 // For unknown reasons, it has been replaced by the following code in Geant???
132 //
133 G4bool inElastic = true;
134 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
135
136 vecLength = 0;
137
138 if (verboseLevel > 1)
139 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
140 << incidentCode << G4endl;
141
142 G4bool successful = false;
143
144 FirstIntInCasNeutron(inElastic, availableEnergy, pv, vecLength,
145 incidentParticle, targetParticle, atomicWeight);
146
147 if (verboseLevel > 1)
148 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
149
150 if ((vecLength > 0) && (availableEnergy > 1.))
151 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
152 pv, vecLength,
153 incidentParticle, targetParticle);
154
155 HighEnergyCascading(successful, pv, vecLength,
156 excitationEnergyGNP, excitationEnergyDTA,
157 incidentParticle, targetParticle,
158 atomicWeight, atomicNumber);
159 if (!successful)
161 excitationEnergyGNP, excitationEnergyDTA,
162 incidentParticle, targetParticle,
163 atomicWeight, atomicNumber);
164 if (!successful)
165 MediumEnergyCascading(successful, pv, vecLength,
166 excitationEnergyGNP, excitationEnergyDTA,
167 incidentParticle, targetParticle,
168 atomicWeight, atomicNumber);
169
170 if (!successful)
172 excitationEnergyGNP, excitationEnergyDTA,
173 incidentParticle, targetParticle,
174 atomicWeight, atomicNumber);
175 if (!successful)
176 QuasiElasticScattering(successful, pv, vecLength,
177 excitationEnergyGNP, excitationEnergyDTA,
178 incidentParticle, targetParticle,
179 atomicWeight, atomicNumber);
180 if (!successful)
181 ElasticScattering(successful, pv, vecLength,
182 incidentParticle,
183 atomicWeight, atomicNumber);
184
185 if (!successful)
186 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
187 << G4endl;
188
190 delete [] pv;
192 return &theParticleChange;
193}
@ stopAndKill
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
bool G4bool
Definition: G4Types.hh:67
#define G4UniformRand()
Definition: Randomize.hh:53
void MediumEnergyClusterProduction(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
void ElasticScattering(G4bool &successful, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, G4double atomicWeight, G4double atomicNumber)
void QuasiElasticScattering(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
void FillParticleChange(G4HEVector pv[], G4int aVecLength)
void HighEnergyClusterProduction(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4double NuclearExcitation(G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber, G4double &excitationEnergyCascade, G4double &excitationEnergyEvaporation)
void MediumEnergyCascading(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
G4double NuclearInelasticity(G4double incidentKineticEnergy, G4double atomicWeight, G4double atomicNumber)
void StrangeParticlePairProduction(const G4double availableEnergy, const G4double centerOfMassEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)
void HighEnergyCascading(G4bool &successful, G4HEVector pv[], G4int &vecLen, G4double &excitationEnergyGNP, G4double &excitationEnergyDTA, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, G4double atomicWeight, G4double atomicNumber)
void FirstIntInCasNeutron(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle, const G4double atomicWeight)
G4double getMass() const
Definition: G4HEVector.cc:361
void setDefinition(G4String name)
Definition: G4HEVector.cc:812
void SetStatusChange(G4HadFinalStateStatus aS)
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115

◆ FirstIntInCasNeutron()

void G4HENeutronInelastic::FirstIntInCasNeutron ( G4bool inElastic,
const G4double  availableEnergy,
G4HEVector  pv[],
G4int vecLen,
const G4HEVector incidentParticle,
const G4HEVector targetParticle,
const G4double  atomicWeight 
)

Definition at line 197 of file G4HENeutronInelastic.cc.

212{
213 static const G4double expxu = 82.; // upper bound for arg. of exp
214 static const G4double expxl = -expxu; // lower bound for arg. of exp
215
216 static const G4double protb = 0.35;
217 static const G4double neutb = 0.35;
218 static const G4double c = 1.25;
219
220 static const G4int numMul = 1200;
221 static const G4int numSec = 60;
222
224 G4int protonCode = Proton.getCode();
225
226 G4int targetCode = targetParticle.getCode();
227 G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
228
229 static G4bool first = true;
230 static G4double protmul[numMul], protnorm[numSec]; // proton constants
231 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
232
233 // misc. local variables
234 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
235
236 G4int i, counter, nt, npos, nneg, nzero;
237
238 if (first) {
239 // compute normalization constants, this will only be done once
240 first = false;
241 for (i = 0; i < numMul; i++) protmul[i] = 0.0;
242 for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
243 counter = -1;
244 for (npos = 0; npos < (numSec/3); npos++) {
245 for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
246 for( nzero=0; nzero<numSec/3; nzero++ )
247 {
248 if( ++counter < numMul )
249 {
250 nt = npos+nneg+nzero;
251 if( (nt>0) && (nt<=numSec) )
252 {
253 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c)
254 /(Factorial(1-npos+nneg)*Factorial(1+npos-nneg)) ;
255 protnorm[nt-1] += protmul[counter];
256 }
257 }
258 }
259 }
260 }
261
262 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
263 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
264 counter = -1;
265 for( npos=0; npos<numSec/3; npos++ )
266 {
267 for( nneg=npos; nneg<=(npos+2); nneg++ )
268 {
269 for( nzero=0; nzero<numSec/3; nzero++ )
270 {
271 if( ++counter < numMul )
272 {
273 nt = npos+nneg+nzero;
274 if( (nt>0) && (nt<=numSec) )
275 {
276 neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c)
277 /(Factorial(-npos+nneg)*Factorial(2+npos-nneg));
278 neutnorm[nt-1] += neutmul[counter];
279 }
280 }
281 }
282 }
283 }
284 for( i=0; i<numSec; i++ )
285 {
286 if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
287 if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
288 }
289 } // end of initialization
290
291
292 // initialize the first two places
293 // the same as beam and target
294 pv[0] = incidentParticle;
295 pv[1] = targetParticle;
296 vecLen = 2;
297
298 if( !inElastic )
299 { // quasi-elastic scattering, no pions produced
300 if( targetCode == protonCode )
301 {
302 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
303 G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5) );
304 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
305 { // charge exchange n p -> p n
306 pv[0] = Proton;
307 pv[1] = Neutron;
308 }
309 }
310 return;
311 }
312 else if (availableEnergy <= PionPlus.getMass())
313 return;
314
315 // inelastic scattering
316
317 npos = 0, nneg = 0, nzero = 0;
318 G4double eab = availableEnergy;
319 G4int ieab = G4int( eab*5.0 );
320
321 G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
322 if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
323 {
324 // suppress high multiplicity events at low momentum
325 // only one additional pion will be produced
326 G4double w0, wp, wm, wt, ran;
327 if( targetCode == neutronCode ) // target is a neutron
328 {
329 w0 = - sqr(1.+neutb)/(2.*c*c);
330 wm = w0 = std::exp(w0);
331 w0 = w0/2.;
332 if( G4UniformRand() < w0/(w0+wm) )
333 { npos = 0; nneg = 0; nzero = 1; }
334 else
335 { npos = 0; nneg = 1; nzero = 0; }
336 }
337 else
338 { // target is a proton
339 w0 = -sqr(1.+protb)/(2.*c*c);
340 w0 = std::exp(w0);
341 wp = w0/2.;
342 wm = -sqr(-1.+protb)/(2.*c*c);
343 wm = std::exp(wm)/2.;
344 wt = w0+wp+wm;
345 wp = w0+wp;
346 ran = G4UniformRand();
347 if( ran < w0/wt)
348 { npos = 0; nneg = 0; nzero = 1; }
349 else if( ran < wp/wt)
350 { npos = 1; nneg = 0; nzero = 0; }
351 else
352 { npos = 0; nneg = 1; nzero = 0; }
353 }
354 }
355 else
356 {
357 // number of total particles vs. centre of mass Energy - 2*proton mass
358
359 G4double aleab = std::log(availableEnergy);
360 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
361 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
362
363 // normalization constant for kno-distribution.
364 // calculate first the sum of all constants, check for numerical problems.
365 G4double test, dum, anpn = 0.0;
366
367 for (nt=1; nt<=numSec; nt++) {
368 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
369 dum = pi*nt/(2.0*n*n);
370 if (std::fabs(dum) < 1.0) {
371 if( test >= 1.0e-10 )anpn += dum*test;
372 } else {
373 anpn += dum*test;
374 }
375 }
376
377 G4double ran = G4UniformRand();
378 G4double excs = 0.0;
379 if( targetCode == protonCode )
380 {
381 counter = -1;
382 for (npos=0; npos<numSec/3; npos++) {
383 for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
384 for (nzero=0; nzero<numSec/3; nzero++) {
385 if (++counter < numMul) {
386 nt = npos+nneg+nzero;
387 if ( (nt>0) && (nt<=numSec) ) {
388 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
389 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
390 if (std::fabs(dum) < 1.0) {
391 if( test >= 1.0e-10 )excs += dum*test;
392 } else {
393 excs += dum*test;
394 }
395 if (ran < excs) goto outOfLoop; //----------------------->
396 }
397 }
398 }
399 }
400 }
401
402 // 3 previous loops continued to the end
403 inElastic = false; // quasi-elastic scattering
404 return;
405 }
406 else
407 { // target must be a neutron
408 counter = -1;
409 for (npos=0; npos<numSec/3; npos++) {
410 for (nneg=npos; nneg<=(npos+2); nneg++) {
411 for (nzero=0; nzero<numSec/3; nzero++) {
412 if (++counter < numMul) {
413 nt = npos+nneg+nzero;
414 if ( (nt>=1) && (nt<=numSec) ) {
415 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
416 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
417 if (std::fabs(dum) < 1.0) {
418 if( test >= 1.0e-10 )excs += dum*test;
419 } else {
420 excs += dum*test;
421 }
422 if (ran < excs) goto outOfLoop; // -------------------------->
423 }
424 }
425 }
426 }
427 }
428 // 3 previous loops continued to the end
429 inElastic = false; // quasi-elastic scattering.
430 return;
431 }
432 }
433 outOfLoop: // <----------------------------------------------
434
435 if( targetCode == neutronCode)
436 {
437 if( npos == nneg)
438 {
439 }
440 else if (npos == (nneg-1))
441 {
442 if( G4UniformRand() < 0.5)
443 {
444 pv[1] = Proton;
445 }
446 else
447 {
448 pv[0] = Proton;
449 }
450 }
451 else
452 {
453 pv[0] = Proton;
454 pv[1] = Proton;
455 }
456 }
457 else
458 {
459 if( npos == nneg)
460 {
461 if( G4UniformRand() < 0.25)
462 {
463 pv[0] = Proton;
464 pv[1] = Neutron;
465 }
466 else
467 {
468 }
469 }
470 else if ( npos == (1+nneg))
471 {
472 pv[1] = Neutron;
473 }
474 else
475 {
476 pv[0] = Proton;
477 }
478 }
479
480
481 nt = npos + nneg + nzero;
482 while ( nt > 0)
483 {
484 G4double ran = G4UniformRand();
485 if ( ran < (G4double)npos/nt)
486 {
487 if( npos > 0 )
488 { pv[vecLen++] = PionPlus;
489 npos--;
490 }
491 }
492 else if ( ran < (G4double)(npos+nneg)/nt)
493 {
494 if( nneg > 0 )
495 {
496 pv[vecLen++] = PionMinus;
497 nneg--;
498 }
499 }
500 else
501 {
502 if( nzero > 0 )
503 {
504 pv[vecLen++] = PionZero;
505 nzero--;
506 }
507 }
508 nt = npos + nneg + nzero;
509 }
510 if (verboseLevel > 1)
511 {
512 G4cout << "Particles produced: " ;
513 G4cout << pv[0].getName() << " " ;
514 G4cout << pv[1].getName() << " " ;
515 for (i=2; i < vecLen; i++)
516 {
517 G4cout << pv[i].getName() << " " ;
518 }
519 G4cout << G4endl;
520 }
521 return;
522 }
@ neutronCode
G4HEVector PionPlus
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4int Factorial(G4int n)
G4HEVector Neutron
G4HEVector PionMinus
G4HEVector PionZero
G4HEVector Proton
G4int getCode() const
Definition: G4HEVector.cc:426
G4double getTotalMomentum() const
Definition: G4HEVector.cc:166
G4String getName() const
Definition: G4HEVector.cc:431
const G4double pi
T sqr(const T &x)
Definition: templates.hh:145

Referenced by ApplyYourself().

◆ GetNumberOfSecondaries()

G4int G4HENeutronInelastic::GetNumberOfSecondaries ( )
inline

Definition at line 75 of file G4HENeutronInelastic.hh.

75{return vecLength;}

◆ ModelDescription()

void G4HENeutronInelastic::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HENeutronInelastic.cc.

43{
44 outFile << "G4HENeutronInelastic is one of the High Energy\n"
45 << "Parameterized (HEP) models used to implement inelastic\n"
46 << "neutron scattering from nuclei. It is a re-engineered\n"
47 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
48 << "initial collision products into backward- and forward-going\n"
49 << "clusters which are then decayed into final state hadrons.\n"
50 << "The model does not conserve energy on an event-by-event\n"
51 << "basis. It may be applied to neutrons with initial energies\n"
52 << "above 20 GeV.\n";
53}

Member Data Documentation

◆ vecLength

G4int G4HENeutronInelastic::vecLength

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