Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HENeutronInelastic.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$
27 //
28 
29 // G4 Process: Gheisha High Energy Collision model.
30 // This includes the high energy cascading model, the two-body-resonance model
31 // and the low energy two-body model. Not included are the low energy stuff
32 // like nuclear reactions, nuclear fission without any cascading and all
33 // processes for particles at rest.
34 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.
35 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
36 
37 #include "G4HENeutronInelastic.hh"
38 #include "globals.hh"
39 #include "G4ios.hh"
40 #include "G4PhysicalConstants.hh"
41 
42 void G4HENeutronInelastic::ModelDescription(std::ostream& outFile) const
43 {
44  outFile << "G4HENeutronInelastic is one of the High Energy\n"
45  << "Parameterized (HEP) models used to implement inelastic\n"
46  << "neutron scattering from nuclei. It is a re-engineered\n"
47  << "version of the GHEISHA code of H. Fesefeldt. It divides the\n"
48  << "initial collision products into backward- and forward-going\n"
49  << "clusters which are then decayed into final state hadrons.\n"
50  << "The model does not conserve energy on an event-by-event\n"
51  << "basis. It may be applied to neutrons with initial energies\n"
52  << "above 20 GeV.\n";
53 }
54 
55 
56 G4HadFinalState*
57 G4HENeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
58  G4Nucleus& targetNucleus)
59 {
60  G4HEVector* pv = new G4HEVector[MAXPART];
61  const G4HadProjectile* aParticle = &aTrack;
62  const G4double A = targetNucleus.GetA_asInt();
63  const G4double Z = targetNucleus.GetZ_asInt();
64  G4HEVector incidentParticle(aParticle);
65 
66  G4double atomicNumber = Z;
67  G4double atomicWeight = A;
68 
69  G4int incidentCode = incidentParticle.getCode();
70  G4double incidentMass = incidentParticle.getMass();
71  G4double incidentTotalEnergy = incidentParticle.getEnergy();
72 
73  // G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
74  // DHW 19 May 2011: variable set but not used
75 
76  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
77 
78  if (incidentKineticEnergy < 1.)
79  G4cout << "GHENeutronInelastic: incident energy < 1 GeV" << G4endl;
80 
81  if (verboseLevel > 1) {
82  G4cout << "G4HENeutronInelastic::ApplyYourself" << G4endl;
83  G4cout << "incident particle " << incidentParticle.getName()
84  << "mass " << incidentMass
85  << "kinetic energy " << incidentKineticEnergy
86  << G4endl;
87  G4cout << "target material with (A,Z) = ("
88  << atomicWeight << "," << atomicNumber << ")" << G4endl;
89  }
90 
91  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy,
92  atomicWeight, atomicNumber);
93  if (verboseLevel > 1)
94  G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
95 
96  incidentKineticEnergy -= inelasticity;
97 
98  G4double excitationEnergyGNP = 0.;
99  G4double excitationEnergyDTA = 0.;
100 
101  G4double excitation = NuclearExcitation(incidentKineticEnergy,
102  atomicWeight, atomicNumber,
103  excitationEnergyGNP,
104  excitationEnergyDTA);
105  if (verboseLevel > 1)
106  G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP
107  << excitationEnergyDTA << G4endl;
108 
109  incidentKineticEnergy -= excitation;
110  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
111  // incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
112  // *(incidentTotalEnergy+incidentMass));
113  // DHW 19 May 2011: variables et but not used
114 
115  G4HEVector targetParticle;
116  if (G4UniformRand() < atomicNumber/atomicWeight) {
117  targetParticle.setDefinition("Proton");
118  } else {
119  targetParticle.setDefinition("Neutron");
120  }
121 
122  G4double targetMass = targetParticle.getMass();
123  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
124  + targetMass*targetMass
125  + 2.0*targetMass*incidentTotalEnergy);
126  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
127 
128  // In the original Gheisha code the inElastic flag was set as follows
129  // G4bool inElastic = InElasticCrossSectionInFirstInt
130  // (availableEnergy, incidentCode, incidentTotalMomentum);
131  // For unknown reasons, it has been replaced by the following code in Geant???
132  //
133  G4bool inElastic = true;
134  // if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;
135 
136  vecLength = 0;
137 
138  if (verboseLevel > 1)
139  G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
140  << incidentCode << G4endl;
141 
142  G4bool successful = false;
143 
144  FirstIntInCasNeutron(inElastic, availableEnergy, pv, vecLength,
145  incidentParticle, targetParticle, atomicWeight);
146 
147  if (verboseLevel > 1)
148  G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
149 
150  if ((vecLength > 0) && (availableEnergy > 1.))
151  StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
152  pv, vecLength,
153  incidentParticle, targetParticle);
154 
155  HighEnergyCascading(successful, pv, vecLength,
156  excitationEnergyGNP, excitationEnergyDTA,
157  incidentParticle, targetParticle,
158  atomicWeight, atomicNumber);
159  if (!successful)
160  HighEnergyClusterProduction(successful, pv, vecLength,
161  excitationEnergyGNP, excitationEnergyDTA,
162  incidentParticle, targetParticle,
163  atomicWeight, atomicNumber);
164  if (!successful)
165  MediumEnergyCascading(successful, pv, vecLength,
166  excitationEnergyGNP, excitationEnergyDTA,
167  incidentParticle, targetParticle,
168  atomicWeight, atomicNumber);
169 
170  if (!successful)
171  MediumEnergyClusterProduction(successful, pv, vecLength,
172  excitationEnergyGNP, excitationEnergyDTA,
173  incidentParticle, targetParticle,
174  atomicWeight, atomicNumber);
175  if (!successful)
176  QuasiElasticScattering(successful, pv, vecLength,
177  excitationEnergyGNP, excitationEnergyDTA,
178  incidentParticle, targetParticle,
179  atomicWeight, atomicNumber);
180  if (!successful)
181  ElasticScattering(successful, pv, vecLength,
182  incidentParticle,
183  atomicWeight, atomicNumber);
184 
185  if (!successful)
186  G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles"
187  << G4endl;
188 
189  FillParticleChange(pv, vecLength);
190  delete [] pv;
191  theParticleChange.SetStatusChange(stopAndKill);
192  return &theParticleChange;
193 }
194 
195 
196 void
197 G4HENeutronInelastic::FirstIntInCasNeutron(G4bool& inElastic,
198  const G4double availableEnergy,
199  G4HEVector pv[],
200  G4int& vecLen,
201  const G4HEVector& incidentParticle,
202  const G4HEVector& targetParticle,
203  const G4double atomicWeight)
204 
205 // Neutron undergoes interaction with nucleon within a nucleus. Check if it is
206 // energetically possible to produce pions/kaons. In not, assume nuclear excitation
207 // occurs and input particle is degraded in energy. No other particles are produced.
208 // If reaction is possible, find the correct number of pions/protons/neutrons
209 // produced using an interpolation to multiplicity data. Replace some pions or
210 // protons/neutrons by kaons or strange baryons according to the average
211 // multiplicity per inelastic reaction.
212 {
213  static const G4double expxu = 82.; // upper bound for arg. of exp
214  static const G4double expxl = -expxu; // lower bound for arg. of exp
215 
216  static const G4double protb = 0.35;
217  static const G4double neutb = 0.35;
218  static const G4double c = 1.25;
219 
220  static const G4int numMul = 1200;
221  static const G4int numSec = 60;
222 
223  G4int neutronCode = Neutron.getCode();
224  G4int protonCode = Proton.getCode();
225 
226  G4int targetCode = targetParticle.getCode();
227  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
228 
229  static G4bool first = true;
230  static G4double protmul[numMul], protnorm[numSec]; // proton constants
231  static G4double neutmul[numMul], neutnorm[numSec]; // neutron constants
232 
233  // misc. local variables
234  // npos = number of pi+, nneg = number of pi-, nzero = number of pi0
235 
236  G4int i, counter, nt, npos, nneg, nzero;
237 
238  if (first) {
239  // compute normalization constants, this will only be done once
240  first = false;
241  for (i = 0; i < numMul; i++) protmul[i] = 0.0;
242  for (i = 0; i < numSec; i++) protnorm[i] = 0.0;
243  counter = -1;
244  for (npos = 0; npos < (numSec/3); npos++) {
245  for (nneg = std::max(0,npos-1); nneg <= (npos+1); nneg++) {
246  for( nzero=0; nzero<numSec/3; nzero++ )
247  {
248  if( ++counter < numMul )
249  {
250  nt = npos+nneg+nzero;
251  if( (nt>0) && (nt<=numSec) )
252  {
253  protmul[counter] = pmltpc(npos,nneg,nzero,nt,protb,c)
254  /(Factorial(1-npos+nneg)*Factorial(1+npos-nneg)) ;
255  protnorm[nt-1] += protmul[counter];
256  }
257  }
258  }
259  }
260  }
261 
262  for( i=0; i<numMul; i++ )neutmul[i] = 0.0;
263  for( i=0; i<numSec; i++ )neutnorm[i] = 0.0;
264  counter = -1;
265  for( npos=0; npos<numSec/3; npos++ )
266  {
267  for( nneg=npos; nneg<=(npos+2); nneg++ )
268  {
269  for( nzero=0; nzero<numSec/3; nzero++ )
270  {
271  if( ++counter < numMul )
272  {
273  nt = npos+nneg+nzero;
274  if( (nt>0) && (nt<=numSec) )
275  {
276  neutmul[counter] = pmltpc(npos,nneg,nzero,nt,neutb,c)
277  /(Factorial(-npos+nneg)*Factorial(2+npos-nneg));
278  neutnorm[nt-1] += neutmul[counter];
279  }
280  }
281  }
282  }
283  }
284  for( i=0; i<numSec; i++ )
285  {
286  if( protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
287  if( neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
288  }
289  } // end of initialization
290 
291 
292  // initialize the first two places
293  // the same as beam and target
294  pv[0] = incidentParticle;
295  pv[1] = targetParticle;
296  vecLen = 2;
297 
298  if( !inElastic )
299  { // quasi-elastic scattering, no pions produced
300  if( targetCode == protonCode )
301  {
302  G4double cech[] = {0.50, 0.45, 0.40, 0.35, 0.30, 0.25, 0.06, 0.04, 0.005, 0.};
303  G4int iplab = G4int( std::min( 9.0, incidentTotalMomentum*2.5) );
304  if( G4UniformRand() < cech[iplab]/std::pow(atomicWeight,0.42) )
305  { // charge exchange n p -> p n
306  pv[0] = Proton;
307  pv[1] = Neutron;
308  }
309  }
310  return;
311  }
312  else if (availableEnergy <= PionPlus.getMass())
313  return;
314 
315  // inelastic scattering
316 
317  npos = 0, nneg = 0, nzero = 0;
318  G4double eab = availableEnergy;
319  G4int ieab = G4int( eab*5.0 );
320 
321  G4double supp[] = {0., 0.4, 0.55, 0.65, 0.75, 0.82, 0.86, 0.90, 0.94, 0.98};
322  if( (ieab <= 9) && (G4UniformRand() >= supp[ieab]) )
323  {
324  // suppress high multiplicity events at low momentum
325  // only one additional pion will be produced
326  G4double w0, wp, wm, wt, ran;
327  if( targetCode == neutronCode ) // target is a neutron
328  {
329  w0 = - sqr(1.+neutb)/(2.*c*c);
330  wm = w0 = std::exp(w0);
331  w0 = w0/2.;
332  if( G4UniformRand() < w0/(w0+wm) )
333  { npos = 0; nneg = 0; nzero = 1; }
334  else
335  { npos = 0; nneg = 1; nzero = 0; }
336  }
337  else
338  { // target is a proton
339  w0 = -sqr(1.+protb)/(2.*c*c);
340  w0 = std::exp(w0);
341  wp = w0/2.;
342  wm = -sqr(-1.+protb)/(2.*c*c);
343  wm = std::exp(wm)/2.;
344  wt = w0+wp+wm;
345  wp = w0+wp;
346  ran = G4UniformRand();
347  if( ran < w0/wt)
348  { npos = 0; nneg = 0; nzero = 1; }
349  else if( ran < wp/wt)
350  { npos = 1; nneg = 0; nzero = 0; }
351  else
352  { npos = 0; nneg = 1; nzero = 0; }
353  }
354  }
355  else
356  {
357  // number of total particles vs. centre of mass Energy - 2*proton mass
358 
359  G4double aleab = std::log(availableEnergy);
360  G4double n = 3.62567+aleab*(0.665843+aleab*(0.336514
361  + aleab*(0.117712+0.0136912*aleab))) - 2.0;
362 
363  // normalization constant for kno-distribution.
364  // calculate first the sum of all constants, check for numerical problems.
365  G4double test, dum, anpn = 0.0;
366 
367  for (nt=1; nt<=numSec; nt++) {
368  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
369  dum = pi*nt/(2.0*n*n);
370  if (std::fabs(dum) < 1.0) {
371  if( test >= 1.0e-10 )anpn += dum*test;
372  } else {
373  anpn += dum*test;
374  }
375  }
376 
377  G4double ran = G4UniformRand();
378  G4double excs = 0.0;
379  if( targetCode == protonCode )
380  {
381  counter = -1;
382  for (npos=0; npos<numSec/3; npos++) {
383  for (nneg=std::max(0,npos-1); nneg<=(npos+1); nneg++) {
384  for (nzero=0; nzero<numSec/3; nzero++) {
385  if (++counter < numMul) {
386  nt = npos+nneg+nzero;
387  if ( (nt>0) && (nt<=numSec) ) {
388  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
389  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
390  if (std::fabs(dum) < 1.0) {
391  if( test >= 1.0e-10 )excs += dum*test;
392  } else {
393  excs += dum*test;
394  }
395  if (ran < excs) goto outOfLoop; //----------------------->
396  }
397  }
398  }
399  }
400  }
401 
402  // 3 previous loops continued to the end
403  inElastic = false; // quasi-elastic scattering
404  return;
405  }
406  else
407  { // target must be a neutron
408  counter = -1;
409  for (npos=0; npos<numSec/3; npos++) {
410  for (nneg=npos; nneg<=(npos+2); nneg++) {
411  for (nzero=0; nzero<numSec/3; nzero++) {
412  if (++counter < numMul) {
413  nt = npos+nneg+nzero;
414  if ( (nt>=1) && (nt<=numSec) ) {
415  test = std::exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
416  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
417  if (std::fabs(dum) < 1.0) {
418  if( test >= 1.0e-10 )excs += dum*test;
419  } else {
420  excs += dum*test;
421  }
422  if (ran < excs) goto outOfLoop; // -------------------------->
423  }
424  }
425  }
426  }
427  }
428  // 3 previous loops continued to the end
429  inElastic = false; // quasi-elastic scattering.
430  return;
431  }
432  }
433  outOfLoop: // <----------------------------------------------
434 
435  if( targetCode == neutronCode)
436  {
437  if( npos == nneg)
438  {
439  }
440  else if (npos == (nneg-1))
441  {
442  if( G4UniformRand() < 0.5)
443  {
444  pv[1] = Proton;
445  }
446  else
447  {
448  pv[0] = Proton;
449  }
450  }
451  else
452  {
453  pv[0] = Proton;
454  pv[1] = Proton;
455  }
456  }
457  else
458  {
459  if( npos == nneg)
460  {
461  if( G4UniformRand() < 0.25)
462  {
463  pv[0] = Proton;
464  pv[1] = Neutron;
465  }
466  else
467  {
468  }
469  }
470  else if ( npos == (1+nneg))
471  {
472  pv[1] = Neutron;
473  }
474  else
475  {
476  pv[0] = Proton;
477  }
478  }
479 
480 
481  nt = npos + nneg + nzero;
482  while ( nt > 0)
483  {
484  G4double ran = G4UniformRand();
485  if ( ran < (G4double)npos/nt)
486  {
487  if( npos > 0 )
488  { pv[vecLen++] = PionPlus;
489  npos--;
490  }
491  }
492  else if ( ran < (G4double)(npos+nneg)/nt)
493  {
494  if( nneg > 0 )
495  {
496  pv[vecLen++] = PionMinus;
497  nneg--;
498  }
499  }
500  else
501  {
502  if( nzero > 0 )
503  {
504  pv[vecLen++] = PionZero;
505  nzero--;
506  }
507  }
508  nt = npos + nneg + nzero;
509  }
510  if (verboseLevel > 1)
511  {
512  G4cout << "Particles produced: " ;
513  G4cout << pv[0].getName() << " " ;
514  G4cout << pv[1].getName() << " " ;
515  for (i=2; i < vecLen; i++)
516  {
517  G4cout << pv[i].getName() << " " ;
518  }
519  G4cout << G4endl;
520  }
521  return;
522  }
523 
524 
525 
526 
527 
528 
529 
530 
531 
Generated on Sat May 25 2013 14:33:56 for Geant4 by   doxygen 1.8.4