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

#include <G4HESigmaPlusInelastic.hh>

+ Inheritance diagram for G4HESigmaPlusInelastic:

Public Member Functions

 G4HESigmaPlusInelastic (const G4String &name="G4HESigmaPlusInelastic")
 
 ~G4HESigmaPlusInelastic ()
 
virtual void ModelDescription (std::ostream &) const
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4int GetNumberOfSecondaries ()
 
void FirstIntInCasSigmaPlus (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 51 of file G4HESigmaPlusInelastic.hh.

Constructor & Destructor Documentation

◆ G4HESigmaPlusInelastic()

G4HESigmaPlusInelastic::G4HESigmaPlusInelastic ( const G4String name = "G4HESigmaPlusInelastic")

Definition at line 43 of file G4HESigmaPlusInelastic.cc.

44 : G4HEInelastic(name)
45{
46 vecLength = 0;
47 theMinEnergy = 20*GeV;
48 theMaxEnergy = 10*TeV;
49 MAXPART = 2048;
50 verboseLevel = 0;
51 G4cout << "WARNING: model G4HESigmaPlusInelastic is being deprecated and will\n"
52 << "disappear in Geant4 version 10.0" << G4endl;
53}
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout

◆ ~G4HESigmaPlusInelastic()

G4HESigmaPlusInelastic::~G4HESigmaPlusInelastic ( )
inline

Definition at line 56 of file G4HESigmaPlusInelastic.hh.

56{};

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 70 of file G4HESigmaPlusInelastic.cc.

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

◆ FirstIntInCasSigmaPlus()

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

Definition at line 203 of file G4HESigmaPlusInelastic.cc.

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

Referenced by ApplyYourself().

◆ GetNumberOfSecondaries()

G4int G4HESigmaPlusInelastic::GetNumberOfSecondaries ( )
inline

Definition at line 65 of file G4HESigmaPlusInelastic.hh.

65{return vecLength;}

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 55 of file G4HESigmaPlusInelastic.cc.

56{
57 outFile << "G4HESigmaPlusInelastic is one of the High Energy\n"
58 << "Parameterized (HEP) models used to implement inelastic\n"
59 << "Sigma+ scattering from nuclei. It is a re-engineered\n"
60 << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
61 << "initial collision products into backward- and forward-going\n"
62 << "clusters which are then decayed into final state hadrons.\n"
63 << "The model does not conserve energy on an event-by-event\n"
64 << "basis. It may be applied to Sigma+ with initial energies\n"
65 << "above 20 GeV.\n";
66}

Member Data Documentation

◆ vecLength

G4int G4HESigmaPlusInelastic::vecLength

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