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