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

#include <G4HEPionMinusInelastic.hh>

+ Inheritance diagram for G4HEPionMinusInelastic:

Public Member Functions

 G4HEPionMinusInelastic ()
 
 ~G4HEPionMinusInelastic ()
 
virtual void ModelDescription (std::ostream &) const
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4int GetNumberOfSecondaries ()
 
void FirstIntInCasPionMinus (G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)
 
- 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 G4HEPionMinusInelastic.hh.

Constructor & Destructor Documentation

◆ G4HEPionMinusInelastic()

G4HEPionMinusInelastic::G4HEPionMinusInelastic ( )
inline

Definition at line 55 of file G4HEPionMinusInelastic.hh.

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

◆ ~G4HEPionMinusInelastic()

G4HEPionMinusInelastic::~G4HEPionMinusInelastic ( )
inline

Definition at line 66 of file G4HEPionMinusInelastic.hh.

66{ };

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 59 of file G4HEPionMinusInelastic.cc.

61{
62 G4HEVector* pv = new G4HEVector[MAXPART];
63 const G4HadProjectile* aParticle = &aTrack;
64 const G4double A = targetNucleus.GetA_asInt();
65 const G4double Z = targetNucleus.GetZ_asInt();
66 G4HEVector incidentParticle(aParticle);
67
68 G4double atomicNumber = Z;
69 G4double atomicWeight = A;
70
71 G4int incidentCode = incidentParticle.getCode();
72 G4double incidentMass = incidentParticle.getMass();
73 G4double incidentTotalEnergy = incidentParticle.getEnergy();
74
75 // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
76 // DHW 19 May 2011: variable set but not used
77
78 G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
79
80 if (incidentKineticEnergy < 1.)
81 G4cout << "GHEPionMinusInelastic: incident energy < 1 GeV" << G4endl;
82
83 if (verboseLevel > 1) {
84 G4cout << "G4HEPionMinusInelastic::ApplyYourself" << G4endl;
85 G4cout << "incident particle " << incidentParticle.getName()
86 << "mass " << incidentMass
87 << "kinetic energy " << incidentKineticEnergy
88 << G4endl;
89 G4cout << "target material with (A,Z) = ("
90 << atomicWeight << "," << atomicNumber << ")" << G4endl;
91 }
92
93 G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
94 atomicWeight, atomicNumber);
95 if (verboseLevel > 1)
96 G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
97
98 incidentKineticEnergy -= inelasticity;
99
100 G4double excitationEnergyGNP = 0.;
101 G4double excitationEnergyDTA = 0.;
102
103 G4double excitation = NuclearExcitation(incidentKineticEnergy,
104 atomicWeight, atomicNumber,
105 excitationEnergyGNP,
106 excitationEnergyDTA);
107 if (verboseLevel > 1)
108 G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
109 << excitationEnergyDTA << G4endl;
110
111 incidentKineticEnergy -= excitation;
112 incidentTotalEnergy = incidentKineticEnergy + incidentMass;
113 // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
114 // *(incidentTotalEnergy+incidentMass));
115 // DHW 19 May 2011: variable set but not used
116
117 G4HEVector targetParticle;
118
119 if (G4UniformRand() < atomicNumber/atomicWeight) {
120 targetParticle.setDefinition("Proton");
121 } else {
122 targetParticle.setDefinition("Neutron");
123 }
124
125 G4double targetMass = targetParticle.getMass();
126 G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
127 + targetMass*targetMass
128 + 2.0*targetMass*incidentTotalEnergy);
129 G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
130
131 // The value of the inElastic flag was originally defined in the Gheisha
132 // stand-alone code 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 FirstIntInCasPionMinus(inElastic, availableEnergy, pv, vecLength,
149 incidentParticle, targetParticle);
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 FirstIntInCasPionMinus(G4bool &inElastic, const G4double availableEnergy, G4HEVector pv[], G4int &vecLen, const G4HEVector &incidentParticle, const G4HEVector &targetParticle)
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

◆ FirstIntInCasPionMinus()

void G4HEPionMinusInelastic::FirstIntInCasPionMinus ( G4bool inElastic,
const G4double  availableEnergy,
G4HEVector  pv[],
G4int vecLen,
const G4HEVector incidentParticle,
const G4HEVector targetParticle 
)

Definition at line 201 of file G4HEPionMinusInelastic.cc.

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

Definition at line 75 of file G4HEPionMinusInelastic.hh.

75{return vecLength;}

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 44 of file G4HEPionMinusInelastic.cc.

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

Member Data Documentation

◆ vecLength

G4int G4HEPionMinusInelastic::vecLength

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