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

#include <G4HEKaonPlusInelastic.hh>

+ Inheritance diagram for G4HEKaonPlusInelastic:

Public Member Functions

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

Constructor & Destructor Documentation

◆ G4HEKaonPlusInelastic()

G4HEKaonPlusInelastic::G4HEKaonPlusInelastic ( )
inline

Definition at line 55 of file G4HEKaonPlusInelastic.hh.

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

◆ ~G4HEKaonPlusInelastic()

G4HEKaonPlusInelastic::~G4HEKaonPlusInelastic ( )
inline

Definition at line 66 of file G4HEKaonPlusInelastic.hh.

66{};

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 57 of file G4HEKaonPlusInelastic.cc.

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

◆ FirstIntInCasKaonPlus()

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

Definition at line 191 of file G4HEKaonPlusInelastic.cc.

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

Definition at line 75 of file G4HEKaonPlusInelastic.hh.

75{return vecLength;}

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 42 of file G4HEKaonPlusInelastic.cc.

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

Member Data Documentation

◆ vecLength

G4int G4HEKaonPlusInelastic::vecLength

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