Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiNeutronInelastic.cc
Go to the documentation of this file.
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: G4RPGAntiNeutronInelastic.cc 94214 2015-11-09 08:18:05Z gcosmo $
27 //
28 
30 #include "G4Exp.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "Randomize.hh"
34 
37  G4Nucleus &targetNucleus )
38 {
39  const G4HadProjectile *originalIncident = &aTrack;
40 
41  // create the target particle
42 
43  G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
44 
45  if( verboseLevel > 1 )
46  {
47  const G4Material *targetMaterial = aTrack.GetMaterial();
48  G4cout << "G4RPGAntiNeutronInelastic::ApplyYourself called" << G4endl;
49  G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
50  G4cout << "target material = " << targetMaterial->GetName() << ", ";
51  G4cout << "target particle = " << originalTarget->GetDefinition()->GetParticleName()
52  << G4endl;
53  }
54  //
55  // Fermi motion and evaporation
56  // As of Geant3, the Fermi energy calculation had not been Done
57  //
58  G4double ek = originalIncident->GetKineticEnergy()/MeV;
59  G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
60  G4ReactionProduct modifiedOriginal;
61  modifiedOriginal = *originalIncident;
62 
63  G4double tkin = targetNucleus.Cinema( ek );
64  ek += tkin;
65  modifiedOriginal.SetKineticEnergy( ek*MeV );
66  G4double et = ek + amas;
67  G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
68  G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
69  if( pp > 0.0 )
70  {
71  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
72  modifiedOriginal.SetMomentum( momentum * (p/pp) );
73  }
74  //
75  // calculate black track energies
76  //
77  tkin = targetNucleus.EvaporationEffects( ek );
78  ek -= tkin;
79  modifiedOriginal.SetKineticEnergy( ek*MeV );
80  et = ek + amas;
81  p = std::sqrt( std::abs((et-amas)*(et+amas)) );
82  pp = modifiedOriginal.GetMomentum().mag()/MeV;
83  if( pp > 0.0 )
84  {
85  G4ThreeVector momentum = modifiedOriginal.GetMomentum();
86  modifiedOriginal.SetMomentum( momentum * (p/pp) );
87  }
88 
89  G4ReactionProduct currentParticle = modifiedOriginal;
90  G4ReactionProduct targetParticle;
91  targetParticle = *originalTarget;
92  currentParticle.SetSide( 1 ); // incident always goes in forward hemisphere
93  targetParticle.SetSide( -1 ); // target always goes in backward hemisphere
94  G4bool incidentHasChanged = false;
95  G4bool targetHasChanged = false;
96  G4bool quasiElastic = false;
97  G4FastVector<G4ReactionProduct,GHADLISTSIZE> vec; // vec will contain the secondary particles
98  G4int vecLen = 0;
99  vec.Initialize( 0 );
100 
101  const G4double cutOff = 0.1*MeV;
102  const G4double anni = std::min( 1.3*currentParticle.GetTotalMomentum()/GeV, 0.4 );
103 
104  if( (currentParticle.GetKineticEnergy()/MeV > cutOff) ||
105  (G4UniformRand() > anni) )
106  Cascade( vec, vecLen,
107  originalIncident, currentParticle, targetParticle,
108  incidentHasChanged, targetHasChanged, quasiElastic );
109  else
110  quasiElastic = true;
111 
112  CalculateMomenta( vec, vecLen,
113  originalIncident, originalTarget, modifiedOriginal,
114  targetNucleus, currentParticle, targetParticle,
115  incidentHasChanged, targetHasChanged, quasiElastic );
116 
117  SetUpChange( vec, vecLen,
118  currentParticle, targetParticle,
119  incidentHasChanged );
120 
121  delete originalTarget;
122  return &theParticleChange;
123 }
124 
125 
126 void G4RPGAntiNeutronInelastic::Cascade(
128  G4int& vecLen,
129  const G4HadProjectile *originalIncident,
130  G4ReactionProduct &currentParticle,
131  G4ReactionProduct &targetParticle,
132  G4bool &incidentHasChanged,
133  G4bool &targetHasChanged,
134  G4bool &quasiElastic )
135 {
136  // Derived from H. Fesefeldt's original FORTRAN code CASNB
137  // AntiNeutron undergoes interaction with nucleon within a nucleus. Check if it is
138  // energetically possible to produce pions/kaons. In not, assume nuclear excitation
139  // occurs and input particle is degraded in energy. No other particles are produced.
140  // If reaction is possible, find the correct number of pions/protons/neutrons
141  // produced using an interpolation to multiplicity data. Replace some pions or
142  // protons/neutrons by kaons or strange baryons according to the average
143  // multiplicity per Inelastic reaction.
144 
145  const G4double mOriginal = originalIncident->GetDefinition()->GetPDGMass()/MeV;
146  const G4double etOriginal = originalIncident->GetTotalEnergy()/MeV;
147  const G4double pOriginal = originalIncident->GetTotalMomentum()/MeV;
148  const G4double targetMass = targetParticle.GetMass()/MeV;
149  G4double centerofmassEnergy = std::sqrt( mOriginal*mOriginal +
150  targetMass*targetMass +
151  2.0*targetMass*etOriginal );
152  G4double availableEnergy = centerofmassEnergy-(targetMass+mOriginal);
153 
154  static G4ThreadLocal G4bool first = true;
155  const G4int numMul = 1200;
156  const G4int numMulA = 400;
157  const G4int numSec = 60;
158  static G4ThreadLocal G4double protmul[numMul], protnorm[numSec]; // proton constants
159  static G4ThreadLocal G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
160  static G4ThreadLocal G4double protmulA[numMulA], protnormA[numSec]; // proton constants
161  static G4ThreadLocal G4double neutmulA[numMulA], neutnormA[numSec]; // neutron constants
162 
163  // np = number of pi+, nneg = number of pi-, nz = number of pi0
164  G4int counter, nt=0, np=0, nneg=0, nz=0;
165  G4double test;
166  const G4double c = 1.25;
167  const G4double b[] = { 0.70, 0.70 };
168  if (first) { // Computation of normalization constants will only be done once
169  first = false;
170  G4int i;
171  for (i = 0; i < numMul; ++i) protmul[i] = 0.0;
172  for (i = 0; i < numSec; ++i) protnorm[i] = 0.0;
173  counter = -1;
174  for (np = 0; np < (numSec/3); ++np) {
175  for (nneg = std::max(0,np-2); nneg <= np; ++nneg) {
176  for (nz = 0; nz < numSec/3; ++nz) {
177  if (++counter < numMul) {
178  nt = np+nneg+nz;
179  if (nt > 0 && nt <= numSec) {
180  protmul[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
181  protnorm[nt-1] += protmul[counter];
182  }
183  }
184  }
185  }
186  }
187 
188  for( i=0; i<numMul; ++i )neutmul[i] = 0.0;
189  for( i=0; i<numSec; ++i )neutnorm[i] = 0.0;
190  counter = -1;
191  for (np = 0; np < numSec/3; ++np) {
192  for (nneg=std::max(0,np-1); nneg<=(np+1); ++nneg ) {
193  for (nz=0; nz<numSec/3; ++nz ) {
194  if (++counter < numMul ) {
195  nt = np+nneg+nz;
196  if ((nt>0) && (nt<=numSec) ) {
197  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
198  neutnorm[nt-1] += neutmul[counter];
199  }
200  }
201  }
202  }
203  }
204 
205  for (i=0; i<numSec; ++i ) {
206  if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
207  if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
208  }
209 
210  // do the same for annihilation channels
211 
212  for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
213  for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
214  counter = -1;
215 
216  for (np=1; np<(numSec/3); ++np ) {
217  nneg = np-1;
218  for (nz=0; nz<numSec/3; ++nz ) {
219  if (++counter < numMulA ) {
220  nt = np+nneg+nz;
221  if (nt>1 && nt<=numSec ) {
222  protmulA[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
223  protnormA[nt-1] += protmulA[counter];
224  }
225  }
226  }
227  }
228 
229  for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
230  for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
231  counter = -1;
232  for (np=0; np<numSec/3; ++np ) {
233  nneg = np;
234  for (nz=0; nz<numSec/3; ++nz) {
235  if (++counter < numMulA ) {
236  nt = np+nneg+nz;
237  if (nt>1 && nt<=numSec ) {
238  neutmulA[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
239  neutnormA[nt-1] += neutmulA[counter];
240  }
241  }
242  }
243  }
244 
245  for (i=0; i<numSec; ++i ) {
246  if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
247  if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
248  }
249  } // end of initialization
250 
251  const G4double expxu = 82.; // upper bound for arg. of exp
252  const G4double expxl = -expxu; // lower bound for arg. of exp
257 
258  // energetically possible to produce pion(s) --> inelastic scattering
259  // otherwise quasi-elastic scattering
260 
261  const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
262  0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
263  0.39,0.36,0.33,0.10,0.01};
264  G4int iplab = G4int( pOriginal/GeV*10.0 );
265  if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
266  if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
267  if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
268  if( iplab > 24 )iplab = 24;
269 
270  if (G4UniformRand() > anhl[iplab] ) {
271  if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
272  {
273  quasiElastic = true;
274  return;
275  }
276  G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
277  const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
278  G4double w0, wp, wt, wm;
279  if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
280  {
281  // suppress high multiplicity events at low momentum
282  // only one pion will be produced
283  //
284  np = nneg = nz = 0;
285  if( targetParticle.GetDefinition() == aProton )
286  {
287  test = G4Exp( std::min( expxu, std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
288  w0 = test;
289  wp = test;
290  if( G4UniformRand() < w0/(w0+wp) )
291  nz = 1;
292  else
293  np = 1;
294  }
295  else // target is a neutron
296  {
297  test = G4Exp( std::min( expxu, std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
298  w0 = test;
299  wp = test;
300  test = G4Exp( std::min( expxu, std::max( expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
301  wm = test;
302  wt = w0+wp+wm;
303  wp += w0;
304  G4double ran = G4UniformRand();
305  if( ran < w0/wt )
306  nz = 1;
307  else if( ran < wp/wt )
308  np = 1;
309  else
310  nneg = 1;
311  }
312  }
313  else // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
314  {
315  G4double n, anpn;
316  GetNormalizationConstant( availableEnergy, n, anpn );
317  G4double ran = G4UniformRand();
318  G4double dum, excs = 0.0;
319  if( targetParticle.GetDefinition() == aProton )
320  {
321  counter = -1;
322  for( np=0; np<numSec/3 && ran>=excs; ++np )
323  {
324  for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
325  {
326  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
327  {
328  if( ++counter < numMul )
329  {
330  nt = np+nneg+nz;
331  if( nt > 0 )
332  {
333  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
334  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
335  if( std::fabs(dum) < 1.0 )
336  {
337  if( test >= 1.0e-10 )excs += dum*test;
338  }
339  else
340  excs += dum*test;
341  }
342  }
343  }
344  }
345  }
346  if( ran >= excs ) // 3 previous loops continued to the end
347  {
348  quasiElastic = true;
349  return;
350  }
351  np--; nneg--; nz--;
352  }
353  else // target must be a neutron
354  {
355  counter = -1;
356  for( np=0; np<numSec/3 && ran>=excs; ++np )
357  {
358  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
359  {
360  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
361  {
362  if( ++counter < numMul )
363  {
364  nt = np+nneg+nz;
365  if( (nt>=1) && (nt<=numSec) )
366  {
367  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
368  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
369  if( std::fabs(dum) < 1.0 )
370  {
371  if( test >= 1.0e-10 )excs += dum*test;
372  }
373  else
374  excs += dum*test;
375  }
376  }
377  }
378  }
379  }
380  if( ran >= excs ) // 3 previous loops continued to the end
381  {
382  quasiElastic = true;
383  return;
384  }
385  np--; nneg--; nz--;
386  }
387  }
388  if( targetParticle.GetDefinition() == aProton )
389  {
390  switch( np-nneg )
391  {
392  case 1:
393  if( G4UniformRand() < 0.5 )
394  {
395  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
396  incidentHasChanged = true;
397  }
398  else
399  {
400  targetParticle.SetDefinitionAndUpdateE( aNeutron );
401  targetHasChanged = true;
402  }
403  break;
404  case 2:
405  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
406  targetParticle.SetDefinitionAndUpdateE( aNeutron );
407  incidentHasChanged = true;
408  targetHasChanged = true;
409  break;
410  default:
411  break;
412  }
413  }
414  else // target must be a neutron
415  {
416  switch( np-nneg )
417  {
418  case 0:
419  if( G4UniformRand() < 0.33 )
420  {
421  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
422  targetParticle.SetDefinitionAndUpdateE( aProton );
423  incidentHasChanged = true;
424  targetHasChanged = true;
425  }
426  break;
427  case 1:
428  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
429  incidentHasChanged = true;
430  break;
431  default:
432  targetParticle.SetDefinitionAndUpdateE( aProton );
433  targetHasChanged = true;
434  break;
435  }
436  }
437  } else { // random number <= anhl[iplab]
438  if (centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV ) {
439  quasiElastic = true;
440  return;
441  }
442 
443  // annihilation channels
444 
445  G4double n, anpn;
446  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
447  G4double ran = G4UniformRand();
448  G4double dum, excs = 0.0;
449 
450  if (targetParticle.GetDefinition() == aProton ) {
451  counter = -1;
452  for (np=1; (np<numSec/3) && (ran>=excs); ++np ) {
453  nneg = np-1;
454  for (nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
455  if (++counter < numMulA) {
456  nt = np+nneg+nz;
457  if (nt>1 && nt<=numSec ) {
458  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
459  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
460  if (std::fabs(dum) < 1.0 ) {
461  if (test >= 1.0e-10) excs += dum*test;
462  }
463  else
464  excs += dum*test;
465  }
466  }
467  }
468  }
469  } else { // target must be a neutron
470  counter = -1;
471  for (np=0; (np<numSec/3) && (ran>=excs); ++np ) {
472  nneg = np;
473  for (nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
474  if (++counter < numMulA ) {
475  nt = np+nneg+nz;
476  if ((nt>1) && (nt<=numSec) ) {
477  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
478  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
479  if (std::fabs(dum) < 1.0 ) {
480  if( test >= 1.0e-10 )excs += dum*test;
481  }
482  else
483  excs += dum*test;
484  }
485  }
486  }
487  }
488  }
489 
490  if (ran >= excs) { // 3 previous loops continued to the end
491  quasiElastic = true;
492  return;
493  }
494  np--; nz--;
495  currentParticle.SetMass( 0.0 );
496  targetParticle.SetMass( 0.0 );
497  }
498 
499  G4int loop = 0;
501  ed << " While count exceeded " << G4endl;
502  while(np+nneg+nz<3) { /* Loop checking, 01.09.2015, D.Wright */
503  nz++;
504  loop++;
505  if (loop > 1000) {
506  G4Exception("G4RPGAntiNeutronInelastic::Cascade()", "HAD_RPG_100", JustWarning, ed);
507  break;
508  }
509  }
510 
511  SetUpPions( np, nneg, nz, vec, vecLen );
512  return;
513 }
514 
515  /* end of file */
516 
G4double EvaporationEffects(G4double kineticEnergy)
Definition: G4Nucleus.cc:278
G4double GetTotalMomentum() const
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void SetUpChange(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
const char * p
Definition: xmltok.h:285
const G4String & GetName() const
Definition: G4Material.hh:178
void SetSide(const G4int sid)
void CalculateMomenta(G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen, const G4HadProjectile *originalIncident, const G4DynamicParticle *originalTarget, G4ReactionProduct &modifiedOriginal, G4Nucleus &targetNucleus, G4ReactionProduct &currentParticle, G4ReactionProduct &targetParticle, G4bool &incidentHasChanged, G4bool &targetHasChanged, G4bool quasiElastic)
G4ParticleDefinition * GetDefinition() const
#define G4ThreadLocal
Definition: tls.hh:89
void Initialize(G4int items)
Definition: G4FastVector.hh:63
int G4int
Definition: G4Types.hh:78
G4DynamicParticle * ReturnTargetParticle() const
Definition: G4Nucleus.cc:241
const G4String & GetParticleName() const
static unsigned wp
Definition: csz_inflate.cc:294
G4double Pmltpc(G4int np, G4int nm, G4int nz, G4int n, G4double b, G4double c)
tuple b
Definition: test.py:12
const G4ParticleDefinition * GetDefinition() const
void SetMass(const G4double mas)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
G4double ek
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4int n
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetKineticEnergy() const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double Cinema(G4double kineticEnergy)
Definition: G4Nucleus.cc:382
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
G4ThreeVector GetMomentum() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
const G4Material * GetMaterial() const
static constexpr double pi
Definition: G4SIunits.hh:75
void GetNormalizationConstant(const G4double availableEnergy, G4double &n, G4double &anpn)
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
void SetUpPions(const G4int np, const G4int nm, const G4int nz, G4FastVector< G4ReactionProduct, 256 > &vec, G4int &vecLen)
program test
Definition: Main_HIJING.f:1
double mag() const
G4double GetMass() const
G4double GetTotalMomentum() const
G4double GetTotalEnergy() const