ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGStrangeProduction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGStrangeProduction.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 
28 #include <iostream>
29 #include <signal.h>
30 
31 #include "G4RPGStrangeProduction.hh"
32 #include "G4Log.hh"
33 #include "Randomize.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4HadReentrentException.hh"
36 
37 G4RPGStrangeProduction::G4RPGStrangeProduction()
38  : G4RPGReaction() {}
39 
40 
41 G4bool G4RPGStrangeProduction::
42 ReactionStage(const G4HadProjectile* /*originalIncident*/,
43  G4ReactionProduct& modifiedOriginal,
44  G4bool& incidentHasChanged,
45  const G4DynamicParticle* originalTarget,
46  G4ReactionProduct& targetParticle,
47  G4bool& targetHasChanged,
48  const G4Nucleus& /*targetNucleus*/,
49  G4ReactionProduct& currentParticle,
50  G4FastVector<G4ReactionProduct,256>& vec,
51  G4int& vecLen,
52  G4bool /*leadFlag*/,
53  G4ReactionProduct& /*leadingStrangeParticle*/)
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 
109  G4int loop = 0;
110  G4ExceptionDescription ed;
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  {
178  G4ReactionProduct *p1 = new G4ReactionProduct;
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;
228  G4ExceptionDescription eda;
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  {
262  G4ReactionProduct *p1 = new G4ReactionProduct;
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;
314  G4ExceptionDescription edb;
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 }
467 
468 
469  /* end of file */