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

#include <G4HEXiZeroInelastic.hh>

+ Inheritance diagram for G4HEXiZeroInelastic:

Public Member Functions

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

Constructor & Destructor Documentation

◆ G4HEXiZeroInelastic()

G4HEXiZeroInelastic::G4HEXiZeroInelastic ( )
inline

Definition at line 55 of file G4HEXiZeroInelastic.hh.

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

◆ ~G4HEXiZeroInelastic()

G4HEXiZeroInelastic::~G4HEXiZeroInelastic ( )
inline

Definition at line 66 of file G4HEXiZeroInelastic.hh.

66{};

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 56 of file G4HEXiZeroInelastic.cc.

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

◆ FirstIntInCasXiZero()

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

Definition at line 189 of file G4HEXiZeroInelastic.cc.

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

Referenced by ApplyYourself().

◆ GetNumberOfSecondaries()

G4int G4HEXiZeroInelastic::GetNumberOfSecondaries ( )
inline

Definition at line 75 of file G4HEXiZeroInelastic.hh.

75{return vecLength;}

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 41 of file G4HEXiZeroInelastic.cc.

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

Member Data Documentation

◆ vecLength

G4int G4HEXiZeroInelastic::vecLength

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