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