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