Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4RPGStrangeProduction Class Reference

#include <G4RPGStrangeProduction.hh>

Inheritance diagram for G4RPGStrangeProduction:
G4RPGReaction

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 49 of file G4RPGStrangeProduction.hh.

Constructor & Destructor Documentation

G4RPGStrangeProduction::G4RPGStrangeProduction ( )

Definition at line 37 of file G4RPGStrangeProduction.cc.

38  : G4RPGReaction() {}

Member Function Documentation

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 42 of file G4RPGStrangeProduction.cc.

References G4AntiLambda::AntiLambda(), G4AntiNeutron::AntiNeutron(), G4AntiProton::AntiProton(), G4AntiSigmaMinus::AntiSigmaMinus(), G4AntiSigmaPlus::AntiSigmaPlus(), G4AntiSigmaZero::AntiSigmaZero(), G4UniformRand, G4ReactionProduct::GetDefinition(), G4DynamicParticle::GetDefinition(), G4ReactionProduct::GetMass(), G4ParticleDefinition::GetPDGMass(), G4ReactionProduct::GetTotalEnergy(), python.hepunit::GeV, G4KaonMinus::KaonMinus(), G4KaonPlus::KaonPlus(), G4KaonZeroLong::KaonZeroLong(), G4KaonZeroShort::KaonZeroShort(), G4Lambda::Lambda(), G4INCL::Math::max(), G4Neutron::Neutron(), G4Proton::Proton(), G4ReactionProduct::SetDefinition(), G4ReactionProduct::SetDefinitionAndUpdateE(), G4FastVector< Type, N >::SetElement(), G4ReactionProduct::SetMayBeKilled(), G4ReactionProduct::SetSide(), G4SigmaMinus::SigmaMinus(), G4SigmaPlus::SigmaPlus(), and G4SigmaZero::SigmaZero().

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  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  {
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  {
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 }
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
void SetMayBeKilled(const G4bool f)
static G4KaonZeroLong * KaonZeroLong()
void SetSide(const G4int sid)
G4ParticleDefinition * GetDefinition() const
static G4AntiSigmaPlus * AntiSigmaPlus()
int G4int
Definition: G4Types.hh:78
void SetDefinitionAndUpdateE(G4ParticleDefinition *aParticleDefinition)
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:99
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
G4ParticleDefinition * GetDefinition() const
static G4AntiSigmaMinus * AntiSigmaMinus()
#define G4UniformRand()
Definition: Randomize.hh:87
static G4KaonZeroShort * KaonZeroShort()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
static G4SigmaMinus * SigmaMinus()
G4double GetTotalEnergy() const
G4double GetPDGMass() const
static G4AntiLambda * AntiLambda()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4AntiSigmaZero * AntiSigmaZero()
void SetDefinition(G4ParticleDefinition *aParticleDefinition)
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
double G4double
Definition: G4Types.hh:76
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
G4double GetMass() const
static G4AntiNeutron * AntiNeutron()

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