Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
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, G4int > GetFinalStateNucleons (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 49 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
81 G4ParticleDefinition *aProton = G4Proton::Proton();
82 G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
83 G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
84 G4ParticleDefinition *anAntiNeutron = G4AntiNeutron::AntiNeutron();
85 G4ParticleDefinition *aSigmaMinus = G4SigmaMinus::SigmaMinus();
86 G4ParticleDefinition *aSigmaPlus = G4SigmaPlus::SigmaPlus();
87 G4ParticleDefinition *aSigmaZero = G4SigmaZero::SigmaZero();
88 G4ParticleDefinition *anAntiSigmaMinus = G4AntiSigmaMinus::AntiSigmaMinus();
89 G4ParticleDefinition *anAntiSigmaPlus = G4AntiSigmaPlus::AntiSigmaPlus();
90 G4ParticleDefinition *anAntiSigmaZero = G4AntiSigmaZero::AntiSigmaZero();
91 G4ParticleDefinition *aKaonMinus = G4KaonMinus::KaonMinus();
92 G4ParticleDefinition *aKaonPlus = G4KaonPlus::KaonPlus();
93 G4ParticleDefinition *aKaonZL = G4KaonZeroLong::KaonZeroLong();
94 G4ParticleDefinition *aKaonZS = G4KaonZeroShort::KaonZeroShort();
95 G4ParticleDefinition *aLambda = G4Lambda::Lambda();
96 G4ParticleDefinition *anAntiLambda = G4AntiLambda::AntiLambda();
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 while( (i<12) && (centerofmassEnergy>avrs[i]) )++i;
109 if( i == 12 )
110 ibin = 11;
111 else
112 ibin = i;
113 //
114 // the fortran code chooses a random replacement of produced kaons
115 // but does not take into account charge conservation
116 //
117 if (vecLen == 1) { // we know that vecLen > 0
118 i3 = 0;
119 i4 = 1; // note that we will be adding a new secondary particle in this case only
120 } else { // otherwise 0 <= i3,i4 < vecLen
121 G4double rnd = G4UniformRand();
122 while (rnd == 1.0) rnd = G4UniformRand();
123 i4 = i3 = G4int(vecLen*rnd);
124 while(i3 == i4) {
125 rnd = G4UniformRand();
126 while( rnd == 1.0 ) rnd = G4UniformRand();
127 i4 = G4int(vecLen*rnd);
128 }
129 }
130
131 // use linear interpolation or extrapolation by y=centerofmassEnergy*x+b
132 //
133 const G4double avkkb[] = { 0.0015, 0.005, 0.012, 0.0285, 0.0525, 0.075,
134 0.0975, 0.123, 0.28, 0.398, 0.495, 0.573 };
135 const G4double avky[] = { 0.005, 0.03, 0.064, 0.095, 0.115, 0.13,
136 0.145, 0.155, 0.20, 0.205, 0.210, 0.212 };
137 const G4double avnnb[] = { 0.00001, 0.0001, 0.0006, 0.0025, 0.01, 0.02,
138 0.04, 0.05, 0.12, 0.15, 0.18, 0.20 };
139
140 avk = (std::log(avkkb[ibin])-std::log(avkkb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
141 /(avrs[ibin]-avrs[ibin-1]) + std::log(avkkb[ibin-1]);
142 avk = std::exp(avk);
143
144 avy = (std::log(avky[ibin])-std::log(avky[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
145 /(avrs[ibin]-avrs[ibin-1]) + std::log(avky[ibin-1]);
146 avy = std::exp(avy);
147
148 avn = (std::log(avnnb[ibin])-std::log(avnnb[ibin-1]))*(centerofmassEnergy-avrs[ibin-1])
149 /(avrs[ibin]-avrs[ibin-1]) + std::log(avnnb[ibin-1]);
150 avn = std::exp(avn);
151
152 if( avk+avy+avn <= 0.0 )return true;
153
154 if( currentMass < protonMass )avy /= 2.0;
155 if( targetMass < protonMass )avy = 0.0;
156 avy += avk+avn;
157 avk += avn;
158 ran = G4UniformRand();
159 if( ran < avn )
160 {
161 if( availableEnergy < 2.0 )return true;
162 if( vecLen == 1 ) // add a new secondary
163 {
164 G4ReactionProduct *p1 = new G4ReactionProduct;
165 if( G4UniformRand() < 0.5 )
166 {
167 vec[0]->SetDefinition( aNeutron );
168 p1->SetDefinition( anAntiNeutron );
169 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
170 vec[0]->SetMayBeKilled(false);
171 p1->SetMayBeKilled(false);
172 }
173 else
174 {
175 vec[0]->SetDefinition( aProton );
176 p1->SetDefinition( anAntiProton );
177 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
178 vec[0]->SetMayBeKilled(false);
179 p1->SetMayBeKilled(false);
180 }
181 vec.SetElement( vecLen++, p1 );
182 // DEBUGGING --> DumpFrames::DumpFrame(vec, vecLen);
183 }
184 else
185 { // replace two secondaries
186 if( G4UniformRand() < 0.5 )
187 {
188 vec[i3]->SetDefinition( aNeutron );
189 vec[i4]->SetDefinition( anAntiNeutron );
190 vec[i3]->SetMayBeKilled(false);
191 vec[i4]->SetMayBeKilled(false);
192 }
193 else
194 {
195 vec[i3]->SetDefinition( aProton );
196 vec[i4]->SetDefinition( anAntiProton );
197 vec[i3]->SetMayBeKilled(false);
198 vec[i4]->SetMayBeKilled(false);
199 }
200 }
201 }
202 else if( ran < avk )
203 {
204 if( availableEnergy < 1.0 )return true;
205
206 const G4double kkb[] = { 0.2500, 0.3750, 0.5000, 0.5625, 0.6250,
207 0.6875, 0.7500, 0.8750, 1.000 };
208 const G4int ipakkb1[] = { 10, 10, 10, 11, 11, 12, 12, 11, 12 };
209 const G4int ipakkb2[] = { 13, 11, 12, 11, 12, 11, 12, 13, 13 };
210 ran = G4UniformRand();
211 i = 0;
212 while( (i<9) && (ran>=kkb[i]) )++i;
213 if( i == 9 )return true;
214 //
215 // ipakkb[] = { 10,13, 10,11, 10,12, 11,11, 11,12, 12,11, 12,12, 11,13, 12,13 };
216 // charge + - + 0 + 0 0 0 0 0 0 0 0 0 0 - 0 -
217 //
218 switch( ipakkb1[i] )
219 {
220 case 10:
221 vec[i3]->SetDefinition( aKaonPlus );
222 vec[i3]->SetMayBeKilled(false);
223 break;
224 case 11:
225 vec[i3]->SetDefinition( aKaonZS );
226 vec[i3]->SetMayBeKilled(false);
227 break;
228 case 12:
229 vec[i3]->SetDefinition( aKaonZL );
230 vec[i3]->SetMayBeKilled(false);
231 break;
232 }
233
234 if( vecLen == 1 ) // add a secondary
235 {
236 G4ReactionProduct *p1 = new G4ReactionProduct;
237 switch( ipakkb2[i] )
238 {
239 case 11:
240 p1->SetDefinition( aKaonZS );
241 p1->SetMayBeKilled(false);
242 break;
243 case 12:
244 p1->SetDefinition( aKaonZL );
245 p1->SetMayBeKilled(false);
246 break;
247 case 13:
248 p1->SetDefinition( aKaonMinus );
249 p1->SetMayBeKilled(false);
250 break;
251 }
252 (G4UniformRand() < 0.5) ? p1->SetSide( -1 ) : p1->SetSide( 1 );
253 vec.SetElement( vecLen++, p1 );
254
255 }
256 else // replace
257 {
258 switch( ipakkb2[i] )
259 {
260 case 11:
261 vec[i4]->SetDefinition( aKaonZS );
262 vec[i4]->SetMayBeKilled(false);
263 break;
264 case 12:
265 vec[i4]->SetDefinition( aKaonZL );
266 vec[i4]->SetMayBeKilled(false);
267 break;
268 case 13:
269 vec[i4]->SetDefinition( aKaonMinus );
270 vec[i4]->SetMayBeKilled(false);
271 break;
272 }
273 }
274 }
275 else if( ran < avy )
276 {
277 if( availableEnergy < 1.6 )return true;
278
279 const G4double ky[] = { 0.200, 0.300, 0.400, 0.550, 0.625, 0.700,
280 0.800, 0.850, 0.900, 0.950, 0.975, 1.000 };
281 const G4int ipaky1[] = { 18, 18, 18, 20, 20, 20, 21, 21, 21, 22, 22, 22 };
282 const G4int ipaky2[] = { 10, 11, 12, 10, 11, 12, 10, 11, 12, 10, 11, 12 };
283 const G4int ipakyb1[] = { 19, 19, 19, 23, 23, 23, 24, 24, 24, 25, 25, 25 };
284 const G4int ipakyb2[] = { 13, 12, 11, 13, 12, 11, 13, 12, 11, 13, 12, 11 };
285 ran = G4UniformRand();
286 i = 0;
287 while( (i<12) && (ran>ky[i]) )++i;
288 if( i == 12 )return true;
289 if( (currentMass<protonMass) || (G4UniformRand()<0.5) )
290 {
291 // ipaky[] = { 18,10, 18,11, 18,12, 20,10, 20,11, 20,12,
292 // 0 + 0 0 0 0 + + + 0 + 0
293 //
294 // 21,10, 21,11, 21,12, 22,10, 22,11, 22,12 }
295 // 0 + 0 0 0 0 - + - 0 - 0
296 switch( ipaky1[i] )
297 {
298 case 18:
299 targetParticle.SetDefinition( aLambda );
300 break;
301 case 20:
302 targetParticle.SetDefinition( aSigmaPlus );
303 break;
304 case 21:
305 targetParticle.SetDefinition( aSigmaZero );
306 break;
307 case 22:
308 targetParticle.SetDefinition( aSigmaMinus );
309 break;
310 }
311 targetHasChanged = true;
312 switch( ipaky2[i] )
313 {
314 case 10:
315 vec[i3]->SetDefinition( aKaonPlus );
316 vec[i3]->SetMayBeKilled(false);
317 break;
318 case 11:
319 vec[i3]->SetDefinition( aKaonZS );
320 vec[i3]->SetMayBeKilled(false);
321 break;
322 case 12:
323 vec[i3]->SetDefinition( aKaonZL );
324 vec[i3]->SetMayBeKilled(false);
325 break;
326 }
327 }
328 else // (currentMass >= protonMass) && (G4UniformRand() >= 0.5)
329 {
330 // ipakyb[] = { 19,13, 19,12, 19,11, 23,13, 23,12, 23,11,
331 // 24,13, 24,12, 24,11, 25,13, 25,12, 25,11 };
332 if( (currentParticle.GetDefinition() == anAntiProton) ||
333 (currentParticle.GetDefinition() == anAntiNeutron) ||
334 (currentParticle.GetDefinition() == anAntiLambda) ||
335 (currentMass > sigmaMinusMass) )
336 {
337 switch( ipakyb1[i] )
338 {
339 case 19:
340 currentParticle.SetDefinitionAndUpdateE( anAntiLambda );
341 break;
342 case 23:
343 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaPlus );
344 break;
345 case 24:
346 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaZero );
347 break;
348 case 25:
349 currentParticle.SetDefinitionAndUpdateE( anAntiSigmaMinus );
350 break;
351 }
352 incidentHasChanged = true;
353 switch( ipakyb2[i] )
354 {
355 case 11:
356 vec[i3]->SetDefinition( aKaonZS );
357 vec[i3]->SetMayBeKilled(false);
358 break;
359 case 12:
360 vec[i3]->SetDefinition( aKaonZL );
361 vec[i3]->SetMayBeKilled(false);
362 break;
363 case 13:
364 vec[i3]->SetDefinition( aKaonMinus );
365 vec[i3]->SetMayBeKilled(false);
366 break;
367 }
368 }
369 else
370 {
371 switch( ipaky1[i] )
372 {
373 case 18:
374 currentParticle.SetDefinitionAndUpdateE( aLambda );
375 break;
376 case 20:
377 currentParticle.SetDefinitionAndUpdateE( aSigmaPlus );
378 break;
379 case 21:
380 currentParticle.SetDefinitionAndUpdateE( aSigmaZero );
381 break;
382 case 22:
383 currentParticle.SetDefinitionAndUpdateE( aSigmaMinus );
384 break;
385 }
386 incidentHasChanged = true;
387 switch( ipaky2[i] )
388 {
389 case 10:
390 vec[i3]->SetDefinition( aKaonPlus );
391 vec[i3]->SetMayBeKilled(false);
392 break;
393 case 11:
394 vec[i3]->SetDefinition( aKaonZS );
395 vec[i3]->SetMayBeKilled(false);
396 break;
397 case 12:
398 vec[i3]->SetDefinition( aKaonZL );
399 vec[i3]->SetMayBeKilled(false);
400 break;
401 }
402 }
403 }
404 }
405 else return true;
406
407 //
408 // check the available energy
409 // if there is not enough energy for kkb/ky pair production
410 // then reduce the number of secondary particles
411 // NOTE:
412 // the number of secondaries may have been changed
413 // the incident and/or target particles may have changed
414 // charge conservation is ignored (as well as strangness conservation)
415 //
416 currentMass = currentParticle.GetMass()/GeV;
417 targetMass = targetParticle.GetMass()/GeV;
418
419 G4double energyCheck = centerofmassEnergy-(currentMass+targetMass);
420 for( i=0; i<vecLen; ++i )
421 {
422 energyCheck -= vec[i]->GetMass()/GeV;
423 if( energyCheck < 0.0 ) // chop off the secondary List
424 {
425 vecLen = std::max( 0, --i ); // looks like a memory leak @@@@@@@@@@@@
426 G4int j;
427 for(j=i; j<vecLen; j++) delete vec[j];
428 break;
429 }
430 }
431
432 return true;
433}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4AntiLambda::AntiLambda
static G4AntiLambda * AntiLambda()
Definition: G4AntiLambda.cc:108
G4AntiNeutron::AntiNeutron
static G4AntiNeutron * AntiNeutron()
Definition: G4AntiNeutron.cc:103
G4AntiProton::AntiProton
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4AntiSigmaMinus::AntiSigmaMinus
static G4AntiSigmaMinus * AntiSigmaMinus()
Definition: G4AntiSigmaMinus.cc:106
G4AntiSigmaPlus::AntiSigmaPlus
static G4AntiSigmaPlus * AntiSigmaPlus()
Definition: G4AntiSigmaPlus.cc:108
G4AntiSigmaZero::AntiSigmaZero
static G4AntiSigmaZero * AntiSigmaZero()
Definition: G4AntiSigmaZero.cc:99
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4FastVector::SetElement
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
G4KaonMinus::KaonMinus
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4KaonPlus::KaonPlus
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4KaonZeroLong::KaonZeroLong
static G4KaonZeroLong * KaonZeroLong()
Definition: G4KaonZeroLong.cc:115
G4KaonZeroShort::KaonZeroShort
static G4KaonZeroShort * KaonZeroShort()
Definition: G4KaonZeroShort.cc:104
G4Lambda::Lambda
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4Proton::Proton
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::SetMayBeKilled
void SetMayBeKilled(const G4bool f)
Definition: G4ReactionProduct.hh:168
G4ReactionProduct::SetSide
void SetSide(const G4int sid)
Definition: G4ReactionProduct.hh:156
G4ReactionProduct::SetDefinitionAndUpdateE
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:142
G4ReactionProduct::GetDefinition
G4ParticleDefinition * GetDefinition() const
Definition: G4ReactionProduct.hh:104
G4ReactionProduct::SetDefinition
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
Definition: G4ReactionProduct.cc:155
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:147
G4SigmaMinus::SigmaMinus
static G4SigmaMinus * SigmaMinus()
Definition: G4SigmaMinus.cc:106
G4SigmaPlus::SigmaPlus
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
G4SigmaZero::SigmaZero
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
The documentation for this class was generated from the following files: