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

#include <G4RPGStrangeProduction.hh>

+ Inheritance diagram for G4RPGStrangeProduction:

Public Member Functions

 G4RPGStrangeProduction ()
 
G4bool ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
 
- Public Member Functions inherited from G4RPGReaction
 G4RPGReaction ()
 
virtual ~G4RPGReaction ()
 
G4bool ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &)
 
void AddBlackTrackParticles (const G4double, const G4int, const G4double, const G4int, const G4ReactionProduct &, G4int, G4int, const G4Nucleus &, G4FastVector< G4ReactionProduct, 256 > &, G4int &)
 
G4double GenerateNBodyEvent (const G4double totalEnergy, const G4bool constantCrossSection, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double GenerateNBodyEventT (const G4double totalEnergy, const G4bool constantCrossSection, std::vector< G4ReactionProduct * > &list)
 
void NuclearReaction (G4FastVector< G4ReactionProduct, 4 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4Nucleus &aNucleus, const G4double theAtomicMass, const G4double *massVec)
 

Additional Inherited Members

- Protected Member Functions inherited from G4RPGReaction
void Rotate (const G4double numberofFinalStateNucleons, const G4ThreeVector &temp, const G4ReactionProduct &modifiedOriginal, const G4HadProjectile *originalIncident, const G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
void Defs1 (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
std::pair< G4int, G4intGetFinalStateNucleons (const G4DynamicParticle *originalTarget, const G4FastVector< G4ReactionProduct, 256 > &vec, const G4int &vecLen)
 
void MomentumCheck (const G4ReactionProduct &modifiedOriginal, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
 
G4double normal ()
 
G4ThreeVector Isotropic (const G4double &)
 

Detailed Description

Definition at line 48 of file G4RPGStrangeProduction.hh.

Constructor & Destructor Documentation

◆ G4RPGStrangeProduction()

G4RPGStrangeProduction::G4RPGStrangeProduction ( )

Definition at line 37 of file G4RPGStrangeProduction.cc.

38 : G4RPGReaction() {}

Member Function Documentation

◆ ReactionStage()

G4bool G4RPGStrangeProduction::ReactionStage ( const G4HadProjectile ,
G4ReactionProduct modifiedOriginal,
G4bool incidentHasChanged,
const G4DynamicParticle originalTarget,
G4ReactionProduct targetParticle,
G4bool targetHasChanged,
const G4Nucleus ,
G4ReactionProduct currentParticle,
G4FastVector< G4ReactionProduct, 256 > &  vec,
G4int vecLen,
G4bool  ,
G4ReactionProduct  
)

Definition at line 41 of file G4RPGStrangeProduction.cc.

54{
55 // Derived from H. Fesefeldt's original FORTRAN code STPAIR
56 //
57 // Choose charge combinations K+ K-, K+ K0B, K0 K0B, K0 K-,
58 // K+ Y0, K0 Y+, K0 Y-
59 // For antibaryon induced reactions half of the cross sections KB YB
60 // pairs are produced. Charge is not conserved, no experimental data available
61 // for exclusive reactions, therefore some average behaviour assumed.
62 // The ratio L/SIGMA is taken as 3:1 (from experimental low energy)
63 //
64
65 if( vecLen == 0 )return true;
66 //
67 // the following protects against annihilation processes
68 //
69 if( currentParticle.GetMass() == 0.0 || targetParticle.GetMass() == 0.0 )return true;
70
71 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV;
72 const G4double mOriginal = modifiedOriginal.GetDefinition()->GetPDGMass()/GeV;
73 G4double targetMass = originalTarget->GetDefinition()->GetPDGMass()/GeV;
74 G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
75 targetMass*targetMass +
76 2.0*targetMass*etOriginal ); // GeV
77 G4double currentMass = currentParticle.GetMass()/GeV;
78 G4double availableEnergy = centerofmassEnergy-(targetMass+currentMass);
79 if( availableEnergy <= 1.0 )return true;
80
97
98 const G4double protonMass = aProton->GetPDGMass()/GeV;
99 const G4double sigmaMinusMass = aSigmaMinus->GetPDGMass()/GeV;
100 //
101 // determine the center of mass energy bin
102 //
103 const G4double avrs[] = {3.,4.,5.,6.,7.,8.,9.,10.,20.,30.,40.,50.};
104
105 G4int ibin, i3, i4;
106 G4double avk, avy, avn, ran;
107 G4int i = 1;
108
109 G4int loop = 0;
111 ed << " While count exceeded " << G4endl;
112 while ((i<12) && (centerofmassEnergy>avrs[i]) ) { /* Loop checking, 01.09.2015, D.Wright */
113 i++;
114 loop++;
115 if (loop > 1000) {
116 G4Exception("G4RPGStrangeProduction::ReactionStage()", "HAD_RPG_100", JustWarning, ed);
117 break;
118 }
119 }
120
121
122 if( i == 12 )
123 ibin = 11;
124 else
125 ibin = i;
126 //
127 // the fortran code chooses a random replacement of produced kaons
128 // but does not take into account charge conservation
129 //
130 if (vecLen == 1) { // we know that vecLen > 0
131 i3 = 0;
132 i4 = 1; // note that we will be adding a new secondary particle in this case only
133 } else { // otherwise 0 <= i3,i4 < vecLen
134 G4double rnd = G4UniformRand();
135 while (rnd == 1.0) rnd = G4UniformRand(); /* Loop checking, 01.09.2015, D.Wright */
136 i4 = i3 = G4int(vecLen*rnd);
137
138 while(i3 == i4) { /* Loop checking, 01.09.2015, D.Wright */
139 rnd = G4UniformRand();
140 while( rnd == 1.0 ) rnd = G4UniformRand(); /* Loop checking, 01.09.2015, D.Wright */
141 i4 = G4int(vecLen*rnd);
142 }
143 }
144
145 // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
146 //
147 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
148 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
149 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
150 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
151 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
152 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
153
154 avk = (G4Log(avkkb[ibin])-G4Log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
155 /(avrs[ibin]-avrs[ibin-1]) + G4Log(avkkb[ibin-1]);
156 avk = G4Exp(avk);
157
158 avy = (G4Log(avky[ibin])-G4Log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
159 /(avrs[ibin]-avrs[ibin-1]) + G4Log(avky[ibin-1]);
160 avy = G4Exp(avy);
161
162 avn = (G4Log(avnnb[ibin])-G4Log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
163 /(avrs[ibin]-avrs[ibin-1]) + G4Log(avnnb[ibin-1]);
164 avn = G4Exp(avn);
165
166 if( avk+avy+avn <= 0.0 )return true;
167
168 if( currentMass < protonMass )avy /= 2.0;
169 if( targetMass < protonMass )avy = 0.0;
170 avy += avk+avn;
171 avk += avn;
172 ran = G4UniformRand();
173 if( ran < avn )
174 {
175 if( availableEnergy < 2.0 )return true;
176 if( vecLen == 1 ) // add a new secondary
177 {
179 if( G4UniformRand() < 0.5 )
180 {
181 vec[0]->SetDefinition( aNeutron );
182 p1->SetDefinition( anAntiNeutron );
183 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
184 vec[0]->SetMayBeKilled(false);
185 p1->SetMayBeKilled(false);
186 }
187 else
188 {
189 vec[0]->SetDefinition( aProton );
190 p1->SetDefinition( anAntiProton );
191 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
192 vec[0]->SetMayBeKilled(false);
193 p1->SetMayBeKilled(false);
194 }
195 vec.SetElement( vecLen++, p1 );
196 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
197 }
198 else
199 { // replace two secondaries
200 if( G4UniformRand() < 0.5 )
201 {
202 vec[i3]->SetDefinition( aNeutron );
203 vec[i4]->SetDefinition( anAntiNeutron );
204 vec[i3]->SetMayBeKilled(false);
205 vec[i4]->SetMayBeKilled(false);
206 }
207 else
208 {
209 vec[i3]->SetDefinition( aProton );
210 vec[i4]->SetDefinition( anAntiProton );
211 vec[i3]->SetMayBeKilled(false);
212 vec[i4]->SetMayBeKilled(false);
213 }
214 }
215 }
216 else if( ran < avk )
217 {
218 if( availableEnergy < 1.0 )return true;
219
220 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
221 0.6875, 0.7500, 0.8750, 1.000 };
222 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
223 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
224 ran = G4UniformRand();
225 i = 0;
226
227 loop = 0;
229 eda << " While count exceeded " << G4endl;
230 while( (i<9) && (ran>=kkb[i]) ) { /* Loop checking, 01.09.2015, D.Wright */
231 ++i;
232 loop++;
233 if (loop > 1000) {
234 G4Exception("G4RPGStrangeProduction::ReactionStage()", "HAD_RPG_100", JustWarning, eda);
235 break;
236 }
237 }
238
239 if( i == 9 )return true;
240 //
241 // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
242 // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
243 //
244 switch( ipakkb1[i] )
245 {
246 case 10:
247 vec[i3]->SetDefinition( aKaonPlus );
248 vec[i3]->SetMayBeKilled(false);
249 break;
250 case 11:
251 vec[i3]->SetDefinition( aKaonZS );
252 vec[i3]->SetMayBeKilled(false);
253 break;
254 case 12:
255 vec[i3]->SetDefinition( aKaonZL );
256 vec[i3]->SetMayBeKilled(false);
257 break;
258 }
259
260 if( vecLen == 1 ) // add a secondary
261 {
263 switch( ipakkb2[i] )
264 {
265 case 11:
266 p1->SetDefinition( aKaonZS );
267 p1->SetMayBeKilled(false);
268 break;
269 case 12:
270 p1->SetDefinition( aKaonZL );
271 p1->SetMayBeKilled(false);
272 break;
273 case 13:
274 p1->SetDefinition( aKaonMinus );
275 p1->SetMayBeKilled(false);
276 break;
277 }
278 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
279 vec.SetElement( vecLen++, p1 );
280
281 }
282 else // replace
283 {
284 switch( ipakkb2[i] )
285 {
286 case 11:
287 vec[i4]->SetDefinition( aKaonZS );
288 vec[i4]->SetMayBeKilled(false);
289 break;
290 case 12:
291 vec[i4]->SetDefinition( aKaonZL );
292 vec[i4]->SetMayBeKilled(false);
293 break;
294 case 13:
295 vec[i4]->SetDefinition( aKaonMinus );
296 vec[i4]->SetMayBeKilled(false);
297 break;
298 }
299 }
300
301 } else if( ran < avy ) {
302 if( availableEnergy < 1.6 )return true;
303
304 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
305 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
306 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
307 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
308 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
309 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
310 ran = G4UniformRand();
311 i = 0;
312
313 loop = 0;
315 edb << " While count exceeded " << G4endl;
316 while( (i<12) && (ran>ky[i]) ) { /* Loop checking, 01.09.2015, D.Wright */
317 ++i;
318 loop++;
319 if (loop > 1000) {
320 G4Exception("G4RPGStrangeProduction::ReactionStage()", "HAD_RPG_100", JustWarning, edb);
321 break;
322 }
323 }
324
325 if( i == 12 )return true;
326 if ( (currentMass<protonMass) || (G4UniformRand()<0.5) ) {
327 // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
328 // 0 + 0 0 0 0 + + + 0 + 0
329 //
330 // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
331 // 0 + 0 0 0 0 - + - 0 - 0
332 switch( ipaky1[i] )
333 {
334 case 18:
335 targetParticle.SetDefinition( aLambda );
336 break;
337 case 20:
338 targetParticle.SetDefinition( aSigmaPlus );
339 break;
340 case 21:
341 targetParticle.SetDefinition( aSigmaZero );
342 break;
343 case 22:
344 targetParticle.SetDefinition( aSigmaMinus );
345 break;
346 }
347 targetHasChanged = true;
348 switch( ipaky2[i] )
349 {
350 case 10:
351 vec[i3]->SetDefinition( aKaonPlus );
352 vec[i3]->SetMayBeKilled(false);
353 break;
354 case 11:
355 vec[i3]->SetDefinition( aKaonZS );
356 vec[i3]->SetMayBeKilled(false);
357 break;
358 case 12:
359 vec[i3]->SetDefinition( aKaonZL );
360 vec[i3]->SetMayBeKilled(false);
361 break;
362 }
363
364 } else { // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
365 // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
366 // 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
367 if ( (currentParticle.GetDefinition() == anAntiProton) ||
368 (currentParticle.GetDefinition() == anAntiNeutron) ||
369 (currentParticle.GetDefinition() == anAntiLambda) ||
370 (currentMass > sigmaMinusMass) ) {
371 switch( ipakyb1[i] )
372 {
373 case 19:
374 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
375 break;
376 case 23:
377 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
378 break;
379 case 24:
380 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
381 break;
382 case 25:
383 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
384 break;
385 }
386 incidentHasChanged = true;
387 switch( ipakyb2[i] )
388 {
389 case 11:
390 vec[i3]->SetDefinition( aKaonZS );
391 vec[i3]->SetMayBeKilled(false);
392 break;
393 case 12:
394 vec[i3]->SetDefinition( aKaonZL );
395 vec[i3]->SetMayBeKilled(false);
396 break;
397 case 13:
398 vec[i3]->SetDefinition( aKaonMinus );
399 vec[i3]->SetMayBeKilled(false);
400 break;
401 }
402
403 } else {
404 switch( ipaky1[i] )
405 {
406 case 18:
407 currentParticle.SetDefinitionAndUpdateE( aLambda );
408 break;
409 case 20:
410 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
411 break;
412 case 21:
413 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
414 break;
415 case 22:
416 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
417 break;
418 }
419 incidentHasChanged = true;
420 switch( ipaky2[i] )
421 {
422 case 10:
423 vec[i3]->SetDefinition( aKaonPlus );
424 vec[i3]->SetMayBeKilled(false);
425 break;
426 case 11:
427 vec[i3]->SetDefinition( aKaonZS );
428 vec[i3]->SetMayBeKilled(false);
429 break;
430 case 12:
431 vec[i3]->SetDefinition( aKaonZL );
432 vec[i3]->SetMayBeKilled(false);
433 break;
434 }
435 }
436 }
437 }
438 else return true;
439
440 //
441 // check the available energy
442 // if there is not enough energy for kkb/ky pair production
443 // then reduce the number of secondary particles
444 // NOTE:
445 // the number of secondaries may have been changed
446 // the incident and/or target particles may have changed
447 // charge conservation is ignored (as well as strangness conservation)
448 //
449 currentMass = currentParticle.GetMass()/GeV;
450 targetMass = targetParticle.GetMass()/GeV;
451
452 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
453 for( i=0; i<vecLen; ++i )
454 {
455 energyCheck -= vec[i]->GetMass()/GeV;
456 if( energyCheck < 0.0 ) // chop off the secondary List
457 {
458 vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
459 G4int j;
460 for(j=i; j<vecLen; j++) delete vec[j];
461 break;
462 }
463 }
464
465 return true;
466}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
G4double G4Log(G4double x)
Definition: G4Log.hh:226
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
#define G4UniformRand()
Definition: Randomize.hh:52
static G4AntiLambda * AntiLambda()
static G4AntiNeutron * AntiNeutron()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:92
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4AntiSigmaPlus * AntiSigmaPlus()
static G4AntiSigmaZero * AntiSigmaZero()
G4ParticleDefinition * GetDefinition() const
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:72
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:112
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:112
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4Lambda * Lambda()
Definition: G4Lambda.cc:107
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
static G4Proton * Proton()
Definition: G4Proton.cc:92
const G4ParticleDefinition * GetDefinition() const
G4double GetTotalEnergy() const
void SetMayBeKilled(const G4bool f)
void SetSide(const G4int sid)
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double GetMass() const
static G4SigmaMinus * SigmaMinus()
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:107
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:101

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