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

#include <G4HEProtonInelastic.hh>

+ Inheritance diagram for G4HEProtonInelastic:

Public Member Functions

 G4HEProtonInelastic ()
 
 ~G4HEProtonInelastic ()
 
virtual void ModelDescription (std::ostream &) const
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4int GetNumberOfSecondaries ()
 
void FirstIntInCasProton (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 G4HEProtonInelastic.hh.

Constructor & Destructor Documentation

◆ G4HEProtonInelastic()

G4HEProtonInelastic::G4HEProtonInelastic ( )
inline

Definition at line 55 of file G4HEProtonInelastic.hh.

55 : G4HEInelastic("G4HEProtonInelastic")
56 {
57 vecLength = 0;
58 theMinEnergy = 45*CLHEP::GeV;
59 theMaxEnergy = 10*CLHEP::TeV;
60 MAXPART = 2048;
61 verboseLevel = 0;
62 G4cout << "WARNING: model G4HEProtonInelastic 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

◆ ~G4HEProtonInelastic()

G4HEProtonInelastic::~G4HEProtonInelastic ( )
inline

Definition at line 66 of file G4HEProtonInelastic.hh.

66{};

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 57 of file G4HEProtonInelastic.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 if (verboseLevel > 1)
70 G4cout << "Z , A = " << atomicNumber << " " << atomicWeight << G4endl;
71
72 G4int incidentCode = incidentParticle.getCode();
73 G4double incidentMass = incidentParticle.getMass();
74 G4double incidentTotalEnergy = incidentParticle.getEnergy();
75
76 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
77 // DHW 19 May 2011: variable set but not used
78
79 G4double incidentKineticEnergy = incidentParticle.getKineticEnergy();
80
81 if (incidentKineticEnergy < 1.)
82 G4cout << "GHEProtonInelastic: incident energy < 1 GeV" << G4endl;
83
84 if (verboseLevel > 1) {
85 G4cout << "G4HEProtonInelastic::ApplyYourself" << G4endl;
86 G4cout << "incident particle " << incidentParticle.getName() << " "
87 << "mass " << incidentMass << " "
88 << "kinetic energy " << incidentKineticEnergy
89 << G4endl;
90 G4cout << "target material with (A,Z) = ("
91 << atomicWeight << "," << atomicNumber << ")" << G4endl;
92 }
93
94 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
95 atomicWeight, atomicNumber);
96 if (verboseLevel > 1)
97 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
98
99 incidentKineticEnergy -= inelasticity;
100
101 G4double excitationEnergyGNP = 0.;
102 G4double excitationEnergyDTA = 0.;
103
104 G4double excitation = NuclearExcitation(incidentKineticEnergy,
105 atomicWeight, atomicNumber,
106 excitationEnergyGNP,
107 excitationEnergyDTA);
108
109 if (verboseLevel > 1)
110 G4cout << "nuclear excitation = " << excitation << " "
111 << excitationEnergyGNP << " " << excitationEnergyDTA << G4endl;
112
113 incidentKineticEnergy -= excitation;
114 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
115 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
116 // *(incidentTotalEnergy+incidentMass));
117 // DHW 19 May 2011: variable set but not used
118
119 G4HEVector targetParticle;
120 if (G4UniformRand() < atomicNumber/atomicWeight) {
121 targetParticle.setDefinition("Proton");
122 } else {
123 targetParticle.setDefinition("Neutron");
124 }
125
126 G4double targetMass = targetParticle.getMass();
127 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
128 + targetMass*targetMass
129 + 2.0*targetMass*incidentTotalEnergy);
130 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
131
132 // In the original Gheisha code, the inElastic flag was defined as follows:
133 // G4bool inElastic = InElasticCrossSectionInFirstInt
134 // (availableEnergy, incidentCode, incidentTotalMomentum);
135 // For unknown reasons, it was replaced by the following code in Geant???
136
137 G4bool inElastic = true;
138 // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
139
140 vecLength = 0;
141
142 if (verboseLevel > 1)
143 G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
144 << incidentCode << G4endl;
145
146 G4bool successful = false;
147
148 FirstIntInCasProton(inElastic, availableEnergy, pv, vecLength,
149 incidentParticle, targetParticle, atomicWeight);
150
151 if (verboseLevel > 1)
152 G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
153
154 if ((vecLength > 0) && (availableEnergy > 1.))
155 StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
156 pv, vecLength,
157 incidentParticle, targetParticle);
158
159 HighEnergyCascading(successful, pv, vecLength,
160 excitationEnergyGNP, excitationEnergyDTA,
161 incidentParticle, targetParticle,
162 atomicWeight, atomicNumber);
163 if (!successful)
165 excitationEnergyGNP, excitationEnergyDTA,
166 incidentParticle, targetParticle,
167 atomicWeight, atomicNumber);
168 if (!successful)
169 MediumEnergyCascading(successful, pv, vecLength,
170 excitationEnergyGNP, excitationEnergyDTA,
171 incidentParticle, targetParticle,
172 atomicWeight, atomicNumber);
173
174 if (!successful)
176 excitationEnergyGNP, excitationEnergyDTA,
177 incidentParticle, targetParticle,
178 atomicWeight, atomicNumber);
179 if (!successful)
180 QuasiElasticScattering(successful, pv, vecLength,
181 excitationEnergyGNP, excitationEnergyDTA,
182 incidentParticle, targetParticle,
183 atomicWeight, atomicNumber);
184 if (!successful)
185 ElasticScattering(successful, pv, vecLength,
186 incidentParticle,
187 atomicWeight, atomicNumber);
188
189 if (!successful)
190 G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
191 << G4endl;
192
194 delete [] pv;
196 return &theParticleChange;
197}
@ 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 FirstIntInCasProton(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

◆ FirstIntInCasProton()

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

Definition at line 201 of file G4HEProtonInelastic.cc.

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

Definition at line 75 of file G4HEProtonInelastic.hh.

75{return vecLength;}

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HEProtonInelastic.cc.

43{
44 outFile << "G4HEProtonInelastic is one of the High Energy\n"
45 << "Parameterized (HEP) models used to implement inelastic\n"
46 << "proton 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 protons with initial energies\n"
52 << "above 20 GeV.\n";
53}

Member Data Documentation

◆ vecLength

G4int G4HEProtonInelastic::vecLength

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