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

#include <G4HEAntiLambdaInelastic.hh>

+ Inheritance diagram for G4HEAntiLambdaInelastic:

Public Member Functions

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

Constructor & Destructor Documentation

◆ G4HEAntiLambdaInelastic()

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

Definition at line 43 of file G4HEAntiLambdaInelastic.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 G4HEAntiLambdaInelastic 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

◆ ~G4HEAntiLambdaInelastic()

G4HEAntiLambdaInelastic::~G4HEAntiLambdaInelastic ( )
inline

Definition at line 56 of file G4HEAntiLambdaInelastic.hh.

56{};

Member Function Documentation

◆ ApplyYourself()

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

Implements G4HadronicInteraction.

Definition at line 71 of file G4HEAntiLambdaInelastic.cc.

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

Referenced by G4HEAntiSigmaZeroInelastic::ApplyYourself().

◆ FirstIntInCasAntiLambda()

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

Definition at line 197 of file G4HEAntiLambdaInelastic.cc.

213{
214 static const G4double expxu = 82.; // upper bound for arg. of exp
215 static const G4double expxl = -expxu; // lower bound for arg. of exp
216
217 static const G4double protb = 0.7;
218 static const G4double neutb = 0.7;
219 static const G4double c = 1.25;
220
221 static const G4int numMul = 1200;
222 static const G4int numMulAn = 400;
223 static const G4int numSec = 60;
224
225 G4int protonCode = Proton.getCode();
226
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 protmulAn[numMulAn],protnormAn[numSec];
233 static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
234 static G4double neutmulAn[numMulAn],neutnormAn[numSec];
235
236 // misc. local variables
237 // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
238
239 G4int i, counter, nt, npos, nneg, nzero;
240
241 if( first )
242 { // compute normalization constants, this will only be done once
243 first = false;
244 for( i=0; i<numMul ; i++ ) protmul[i] = 0.0;
245 for( i=0; i<numSec ; i++ ) protnorm[i] = 0.0;
246 counter = -1;
247 for( npos=0; npos<(numSec/3); npos++ )
248 {
249 for( nneg=std::max(0,npos-2); nneg<=(npos+1); nneg++ )
250 {
251 for( nzero=0; nzero<numSec/3; nzero++ )
252 {
253 if( ++counter < numMul )
254 {
255 nt = npos+nneg+nzero;
256 if( (nt>0) && (nt<=numSec) )
257 {
258 protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
259 protnorm[nt-1] += protmul[counter];
260 }
261 }
262 }
263 }
264 }
265 for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
266 for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
267 counter = -1;
268 for( npos=0; npos<numSec/3; npos++ )
269 {
270 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
271 {
272 for( nzero=0; nzero<numSec/3; nzero++ )
273 {
274 if( ++counter < numMul )
275 {
276 nt = npos+nneg+nzero;
277 if( (nt>0) && (nt<=numSec) )
278 {
279 neutmul[counter] = pmltpc(npos,nneg,nzero,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 // annihilation
292 for( i=0; i<numMulAn ; i++ ) protmulAn[i] = 0.0;
293 for( i=0; i<numSec ; i++ ) protnormAn[i] = 0.0;
294 counter = -1;
295 for( npos=1; npos<(numSec/3); npos++ )
296 {
297 nneg = std::max(0,npos-1);
298 for( nzero=0; nzero<numSec/3; nzero++ )
299 {
300 if( ++counter < numMulAn )
301 {
302 nt = npos+nneg+nzero;
303 if( (nt>1) && (nt<=numSec) )
304 {
305 protmulAn[counter] = pmltpc(npos,nneg,nzero,nt,protb,c);
306 protnormAn[nt-1] += protmulAn[counter];
307 }
308 }
309 }
310 }
311 for( i=0; i<numMulAn; i++ ) neutmulAn[i] = 0.0;
312 for( i=0; i<numSec; i++ ) neutnormAn[i] = 0.0;
313 counter = -1;
314 for( npos=0; npos<numSec/3; npos++ )
315 {
316 nneg = npos;
317 for( nzero=0; nzero<numSec/3; nzero++ )
318 {
319 if( ++counter < numMulAn )
320 {
321 nt = npos+nneg+nzero;
322 if( (nt>1) && (nt<=numSec) )
323 {
324 neutmulAn[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c);
325 neutnormAn[nt-1] += neutmulAn[counter];
326 }
327 }
328 }
329 }
330 for( i=0; i<numSec; i++ )
331 {
332 if( protnormAn[i] > 0.0 )protnormAn[i] = 1.0/protnormAn[i];
333 if( neutnormAn[i] > 0.0 )neutnormAn[i] = 1.0/neutnormAn[i];
334 }
335 } // end of initialization
336
337
338 // initialize the first two places
339 // the same as beam and target
340 pv[0] = incidentParticle;
341 pv[1] = targetParticle;
342 vecLen = 2;
343
344 if( !inElastic )
345 { // some two-body reactions
346 G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
347
348 G4int iplab = std::min(9, G4int( incidentTotalMomentum*2.5));
349 if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
350 {
351 G4double ran = G4UniformRand();
352
353 if ( targetCode == protonCode)
354 {
355 if(ran < 0.2)
356 {
357 pv[0] = AntiSigmaZero;
358 }
359 else if (ran < 0.4)
360 {
361 pv[0] = AntiSigmaMinus;
362 pv[1] = Neutron;
363 }
364 else if (ran < 0.6)
365 {
366 pv[0] = Proton;
367 pv[1] = AntiLambda;
368 }
369 else if (ran < 0.8)
370 {
371 pv[0] = Proton;
372 pv[1] = AntiSigmaZero;
373 }
374 else
375 {
376 pv[0] = Neutron;
377 pv[1] = AntiSigmaMinus;
378 }
379 }
380 else
381 {
382 if (ran < 0.2)
383 {
384 pv[0] = AntiSigmaZero;
385 }
386 else if (ran < 0.4)
387 {
388 pv[0] = AntiSigmaPlus;
389 pv[1] = Proton;
390 }
391 else if (ran < 0.6)
392 {
393 pv[0] = Neutron;
394 pv[1] = AntiLambda;
395 }
396 else if (ran < 0.8)
397 {
398 pv[0] = Neutron;
399 pv[1] = AntiSigmaZero;
400 }
401 else
402 {
403 pv[0] = Proton;
404 pv[1] = AntiSigmaPlus;
405 }
406 }
407 }
408 return;
409 }
410 else if (availableEnergy <= PionPlus.getMass())
411 return;
412
413 // inelastic scattering
414
415 npos = 0; nneg = 0; nzero = 0;
416 G4double anhl[] = {1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.97, 0.88,
417 0.85, 0.81, 0.75, 0.64, 0.64, 0.55, 0.55, 0.45, 0.47, 0.40,
418 0.39, 0.36, 0.33, 0.10, 0.01};
419 G4int iplab = G4int( incidentTotalMomentum*10.);
420 if ( iplab > 9) iplab = 10 + G4int( (incidentTotalMomentum -1.)*5. );
421 if ( iplab > 14) iplab = 15 + G4int( incidentTotalMomentum -2. );
422 if ( iplab > 22) iplab = 23 + G4int( (incidentTotalMomentum -10.)/10.);
423 iplab = std::min(24, iplab);
424
425 if (G4UniformRand() > anhl[iplab]) { // non- annihilation channels
426
427 // number of total particles vs. centre of mass Energy - 2*proton mass
428 G4double aleab = std::log(availableEnergy);
429 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
430 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
431
432 // normalization constant for kno-distribution.
433 // calculate first the sum of all constants, check for numerical problems.
434 G4double test, dum, anpn = 0.0;
435
436 for (nt = 1; nt <= numSec; nt++) {
437 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
438 dum = pi*nt/(2.0*n*n);
439 if (std::fabs(dum) < 1.0) {
440 if (test >= 1.0e-10) anpn += dum*test;
441 } else {
442 anpn += dum*test;
443 }
444 }
445
446 G4double ran = G4UniformRand();
447 G4double excs = 0.0;
448 if (targetCode == protonCode) {
449 counter = -1;
450 for (npos = 0; npos < numSec/3; npos++) {
451 for (nneg = std::max(0,npos-2); nneg <= (npos+1); nneg++) {
452 for (nzero = 0; nzero < numSec/3; nzero++) {
453 if (++counter < numMul) {
454 nt = npos+nneg+nzero;
455 if ((nt > 0) && (nt <= numSec) ) {
456 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
457 dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
458
459 if (std::fabs(dum) < 1.0) {
460 if (test >= 1.0e-10) excs += dum*test;
461 } else {
462 excs += dum*test;
463 }
464
465 if (ran < excs) goto outOfLoop; //----------------------->
466 }
467 }
468 }
469 }
470 }
471 // 3 previous loops continued to the end
472 inElastic = false; // quasi-elastic scattering
473 return;
474 } else { // target must be a neutron
475 counter = -1;
476 for( npos=0; npos<numSec/3; npos++ )
477 {
478 for( nneg=std::max(0,npos-1); nneg<=(npos+2); nneg++ )
479 {
480 for( nzero=0; nzero<numSec/3; nzero++ )
481 {
482 if( ++counter < numMul )
483 {
484 nt = npos+nneg+nzero;
485 if( (nt>0) && (nt<=numSec) )
486 {
487 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
488 dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
489 if (std::fabs(dum) < 1.0) {
490 if( test >= 1.0e-10 )excs += dum*test;
491 } else {
492 excs += dum*test;
493 }
494
495 if (ran < excs) goto outOfLoop; // -------------------------->
496 }
497 }
498 }
499 }
500 }
501 // 3 previous loops continued to the end
502 inElastic = false; // quasi-elastic scattering.
503 return;
504 }
505
506 outOfLoop: // <------------------------------------------------------------------------
507
508 ran = G4UniformRand();
509
510 if( targetCode == protonCode)
511 {
512 if( npos == nneg)
513 {
514 if (ran < 0.40)
515 {
516 }
517 else if (ran < 0.8)
518 {
519 pv[0] = AntiSigmaZero;
520 }
521 else
522 {
523 pv[0] = AntiSigmaMinus;
524 pv[1] = Neutron;
525 }
526 }
527 else if (npos == (nneg+1))
528 {
529 if( ran < 0.25)
530 {
531 pv[1] = Neutron;
532 }
533 else if (ran < 0.5)
534 {
535 pv[0] = AntiSigmaZero;
536 pv[1] = Neutron;
537 }
538 else
539 {
540 pv[0] = AntiSigmaPlus;
541 }
542 }
543 else if (npos == (nneg-1))
544 {
545 pv[0] = AntiSigmaMinus;
546 }
547 else
548 {
549 pv[0] = AntiSigmaPlus;
550 pv[1] = Neutron;
551 }
552 }
553 else
554 {
555 if( npos == nneg)
556 {
557 if (ran < 0.4)
558 {
559 }
560 else if(ran < 0.8)
561 {
562 pv[0] = AntiSigmaZero;
563 }
564 else
565 {
566 pv[0] = AntiSigmaPlus;
567 pv[1] = Proton;
568 }
569 }
570 else if ( npos == (nneg-1))
571 {
572 if (ran < 0.5)
573 {
574 pv[0] = AntiSigmaMinus;
575 }
576 else if (ran < 0.75)
577 {
578 pv[1] = Proton;
579 }
580 else
581 {
582 pv[0] = AntiSigmaZero;
583 pv[1] = Proton;
584 }
585 }
586 else if (npos == (nneg+1))
587 {
588 pv[0] = AntiSigmaPlus;
589 }
590 else
591 {
592 pv[0] = AntiSigmaMinus;
593 pv[1] = Proton;
594 }
595 }
596
597 }
598 else // annihilation
599 {
600 if ( availableEnergy > 2. * PionPlus.getMass() )
601 {
602
603 G4double aleab = std::log(availableEnergy);
604 G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
605 + aleab*(0.117712+0.0136912*aleab))) - 2.0;
606
607 // normalization constant for kno-distribution.
608 // calculate first the sum of all constants, check for numerical problems.
609 G4double test, dum, anpn = 0.0;
610
611 for (nt=2; nt<=numSec; nt++) {
612 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
613 dum = pi*nt/(2.0*n*n);
614
615 if (std::fabs(dum) < 1.0) {
616 if( test >= 1.0e-10 )anpn += dum*test;
617 } else {
618 anpn += dum*test;
619 }
620 }
621
622 G4double ran = G4UniformRand();
623 G4double excs = 0.0;
624 if (targetCode == protonCode) {
625 counter = -1;
626 for (npos=1; npos<numSec/3; npos++) {
627 nneg = npos-1;
628 for( nzero=0; nzero<numSec/3; nzero++ )
629 {
630 if( ++counter < numMulAn )
631 {
632 nt = npos+nneg+nzero;
633 if( (nt>1) && (nt<=numSec) ) {
634 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
635 dum = (pi/anpn)*nt*protmulAn[counter]*protnormAn[nt-1]/(2.0*n*n);
636
637 if (std::fabs(dum) < 1.0) {
638 if( test >= 1.0e-10 )excs += dum*test;
639 } else {
640 excs += dum*test;
641 }
642
643 if (ran < excs) goto outOfLoopAn; //----------------------->
644 }
645 }
646 }
647 }
648 // 3 previous loops continued to the end
649 inElastic = false; // quasi-elastic scattering
650 return;
651
652 } else { // target must be a neutron
653 counter = -1;
654 for (npos=0; npos<numSec/3; npos++) {
655 nneg = npos;
656 for( nzero=0; nzero<numSec/3; nzero++ ) {
657 if (++counter < numMulAn) {
658 nt = npos+nneg+nzero;
659 if( (nt>1) && (nt<=numSec) ) {
660 test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
661 dum = (pi/anpn)*nt*neutmulAn[counter]*neutnormAn[nt-1]/(2.0*n*n);
662
663 if (std::fabs(dum) < 1.0) {
664 if( test >= 1.0e-10 )excs += dum*test;
665 } else {
666 excs += dum*test;
667 }
668
669 if (ran < excs) goto outOfLoopAn; // -------------------------->
670 }
671 }
672 }
673 }
674
675 inElastic = false; // quasi-elastic scattering.
676 return;
677 }
678 outOfLoopAn: // <---------------------------------------------------------
679 vecLen = 0;
680 }
681 }
682
683 nt = npos + nneg + nzero;
684 while ( nt > 0)
685 {
686 G4double ran = G4UniformRand();
687 if ( ran < (G4double)npos/nt)
688 {
689 if( npos > 0 )
690 { pv[vecLen++] = PionPlus;
691 npos--;
692 }
693 }
694 else if ( ran < (G4double)(npos+nneg)/nt)
695 {
696 if( nneg > 0 )
697 {
698 pv[vecLen++] = PionMinus;
699 nneg--;
700 }
701 }
702 else
703 {
704 if( nzero > 0 )
705 {
706 pv[vecLen++] = PionZero;
707 nzero--;
708 }
709 }
710 nt = npos + nneg + nzero;
711 }
712 if (verboseLevel > 1)
713 {
714 G4cout << "Particles produced: " ;
715 G4cout << pv[0].getName() << " " ;
716 G4cout << pv[1].getName() << " " ;
717 for (i=2; i < vecLen; i++)
718 {
719 G4cout << pv[i].getName() << " " ;
720 }
721 G4cout << G4endl;
722 }
723 return;
724 }
G4HEVector PionPlus
G4HEVector AntiSigmaZero
G4double pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
G4HEVector AntiSigmaPlus
G4HEVector Neutron
G4HEVector PionMinus
G4HEVector PionZero
G4HEVector AntiSigmaMinus
G4HEVector AntiLambda
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 G4HEAntiLambdaInelastic::GetNumberOfSecondaries ( )
inline

Definition at line 65 of file G4HEAntiLambdaInelastic.hh.

65{return vecLength;}

Referenced by G4HEAntiSigmaZeroInelastic::ApplyYourself().

◆ ModelDescription()

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

Reimplemented from G4HadronicInteraction.

Definition at line 56 of file G4HEAntiLambdaInelastic.cc.

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

Member Data Documentation

◆ vecLength

G4int G4HEAntiLambdaInelastic::vecLength

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