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