ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RPGAntiNeutronInelastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RPGAntiNeutronInelastic.cc
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 //
27 
28 #include "G4RPGAntiNeutronInelastic.hh"
29 #include "G4Exp.hh"
30 #include "G4PhysicalConstants.hh"
31 #include "G4SystemOfUnits.hh"
32 #include "Randomize.hh"
33 
34 G4HadFinalState*
35 G4RPGAntiNeutronInelastic::ApplyYourself( const G4HadProjectile &aTrack,
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(
126  G4FastVector<G4ReactionProduct,GHADLISTSIZE> &vec,
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  for (nz=0; nz<numSec/3; ++nz ) {
193  if (++counter < numMul ) {
194  nt = np+nneg+nz;
195  if ((nt>0) && (nt<=numSec) ) {
196  neutmul[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
197  neutnorm[nt-1] += neutmul[counter];
198  }
199  }
200  }
201  }
202  }
203 
204  for (i=0; i<numSec; ++i ) {
205  if (protnorm[i] > 0.0 )protnorm[i] = 1.0/protnorm[i];
206  if (neutnorm[i] > 0.0 )neutnorm[i] = 1.0/neutnorm[i];
207  }
208 
209  // do the same for annihilation channels
210 
211  for( i=0; i<numMulA; ++i )protmulA[i] = 0.0;
212  for( i=0; i<numSec; ++i )protnormA[i] = 0.0;
213  counter = -1;
214 
215  for (np=1; np<(numSec/3); ++np ) {
216  nneg = np-1;
217  for (nz=0; nz<numSec/3; ++nz ) {
218  if (++counter < numMulA ) {
219  nt = np+nneg+nz;
220  if (nt>1 && nt<=numSec ) {
221  protmulA[counter] = Pmltpc(np,nneg,nz,nt,b[0],c);
222  protnormA[nt-1] += protmulA[counter];
223  }
224  }
225  }
226  }
227 
228  for( i=0; i<numMulA; ++i )neutmulA[i] = 0.0;
229  for( i=0; i<numSec; ++i )neutnormA[i] = 0.0;
230  counter = -1;
231  for (np=0; np<numSec/3; ++np ) {
232  nneg = np;
233  for (nz=0; nz<numSec/3; ++nz) {
234  if (++counter < numMulA ) {
235  nt = np+nneg+nz;
236  if (nt>1 && nt<=numSec ) {
237  neutmulA[counter] = Pmltpc(np,nneg,nz,nt,b[1],c);
238  neutnormA[nt-1] += neutmulA[counter];
239  }
240  }
241  }
242  }
243 
244  for (i=0; i<numSec; ++i ) {
245  if( protnormA[i] > 0.0 )protnormA[i] = 1.0/protnormA[i];
246  if( neutnormA[i] > 0.0 )neutnormA[i] = 1.0/neutnormA[i];
247  }
248  } // end of initialization
249 
250  const G4double expxu = 82.; // upper bound for arg. of exp
251  const G4double expxl = -expxu; // lower bound for arg. of exp
252  G4ParticleDefinition *aNeutron = G4Neutron::Neutron();
253  G4ParticleDefinition *aProton = G4Proton::Proton();
254  G4ParticleDefinition *anAntiProton = G4AntiProton::AntiProton();
255  G4ParticleDefinition *aPiPlus = G4PionPlus::PionPlus();
256 
257  // energetically possible to produce pion(s) --> inelastic scattering
258  // otherwise quasi-elastic scattering
259 
260  const G4double anhl[] = {1.00,1.00,1.00,1.00,1.00,1.00,1.00,1.00,0.97,0.88,
261  0.85,0.81,0.75,0.64,0.64,0.55,0.55,0.45,0.47,0.40,
262  0.39,0.36,0.33,0.10,0.01};
263  G4int iplab = G4int( pOriginal/GeV*10.0 );
264  if( iplab > 9 )iplab = G4int( (pOriginal/GeV- 1.0)*5.0 ) + 10;
265  if( iplab > 14 )iplab = G4int( pOriginal/GeV- 2.0 ) + 15;
266  if( iplab > 22 )iplab = G4int( (pOriginal/GeV-10.0)/10.0 ) + 23;
267  if( iplab > 24 )iplab = 24;
268 
269  if (G4UniformRand() > anhl[iplab] ) {
270  if( availableEnergy <= aPiPlus->GetPDGMass()/MeV )
271  {
272  quasiElastic = true;
273  return;
274  }
275  G4int ieab = static_cast<G4int>(availableEnergy*5.0/GeV);
276  const G4double supp[] = {0.,0.4,0.55,0.65,0.75,0.82,0.86,0.90,0.94,0.98};
277  G4double w0, wp, wt, wm;
278  if( (availableEnergy < 2.0*GeV) && (G4UniformRand() >= supp[ieab]) )
279  {
280  // suppress high multiplicity events at low momentum
281  // only one pion will be produced
282  //
283  np = nneg = nz = 0;
284  if( targetParticle.GetDefinition() == aProton )
285  {
286  test = G4Exp( std::min( expxu, std::max( expxl, -(1.0+b[0])*(1.0+b[0])/(2.0*c*c) ) ) );
287  w0 = test;
288  wp = test;
289  if( G4UniformRand() < w0/(w0+wp) )
290  nz = 1;
291  else
292  np = 1;
293  }
294  else // target is a neutron
295  {
296  test = G4Exp( std::min( expxu, std::max( expxl, -(1.0+b[1])*(1.0+b[1])/(2.0*c*c) ) ) );
297  w0 = test;
298  wp = test;
299  test = G4Exp( std::min( expxu, std::max( expxl, -(-1.0+b[1])*(-1.0+b[1])/(2.0*c*c) ) ) );
300  wm = test;
301  wt = w0+wp+wm;
302  wp += w0;
303  G4double ran = G4UniformRand();
304  if( ran < w0/wt )
305  nz = 1;
306  else if( ran < wp/wt )
307  np = 1;
308  else
309  nneg = 1;
310  }
311  }
312  else // (availableEnergy >= 2.0*GeV) || (random number < supp[ieab])
313  {
314  G4double n, anpn;
315  GetNormalizationConstant( availableEnergy, n, anpn );
316  G4double ran = G4UniformRand();
317  G4double dum, excs = 0.0;
318  if( targetParticle.GetDefinition() == aProton )
319  {
320  counter = -1;
321  for( np=0; np<numSec/3 && ran>=excs; ++np )
322  {
323  for( nneg=std::max(0,np-2); nneg<=np && ran>=excs; ++nneg )
324  {
325  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
326  {
327  if( ++counter < numMul )
328  {
329  nt = np+nneg+nz;
330  if( nt > 0 )
331  {
332  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
333  dum = (pi/anpn)*nt*protmul[counter]*protnorm[nt-1]/(2.0*n*n);
334  if( std::fabs(dum) < 1.0 )
335  {
336  if( test >= 1.0e-10 )excs += dum*test;
337  }
338  else
339  excs += dum*test;
340  }
341  }
342  }
343  }
344  }
345  if( ran >= excs ) // 3 previous loops continued to the end
346  {
347  quasiElastic = true;
348  return;
349  }
350  np--; nneg--; nz--;
351  }
352  else // target must be a neutron
353  {
354  counter = -1;
355  for( np=0; np<numSec/3 && ran>=excs; ++np )
356  {
357  for( nneg=std::max(0,np-1); nneg<=(np+1) && ran>=excs; ++nneg )
358  {
359  for( nz=0; nz<numSec/3 && ran>=excs; ++nz )
360  {
361  if( ++counter < numMul )
362  {
363  nt = np+nneg+nz;
364  if( (nt>=1) && (nt<=numSec) )
365  {
366  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
367  dum = (pi/anpn)*nt*neutmul[counter]*neutnorm[nt-1]/(2.0*n*n);
368  if( std::fabs(dum) < 1.0 )
369  {
370  if( test >= 1.0e-10 )excs += dum*test;
371  }
372  else
373  excs += dum*test;
374  }
375  }
376  }
377  }
378  }
379  if( ran >= excs ) // 3 previous loops continued to the end
380  {
381  quasiElastic = true;
382  return;
383  }
384  np--; nneg--; nz--;
385  }
386  }
387  if( targetParticle.GetDefinition() == aProton )
388  {
389  switch( np-nneg )
390  {
391  case 1:
392  if( G4UniformRand() < 0.5 )
393  {
394  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
395  incidentHasChanged = true;
396  }
397  else
398  {
399  targetParticle.SetDefinitionAndUpdateE( aNeutron );
400  targetHasChanged = true;
401  }
402  break;
403  case 2:
404  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
405  targetParticle.SetDefinitionAndUpdateE( aNeutron );
406  incidentHasChanged = true;
407  targetHasChanged = true;
408  break;
409  default:
410  break;
411  }
412  }
413  else // target must be a neutron
414  {
415  switch( np-nneg )
416  {
417  case 0:
418  if( G4UniformRand() < 0.33 )
419  {
420  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
421  targetParticle.SetDefinitionAndUpdateE( aProton );
422  incidentHasChanged = true;
423  targetHasChanged = true;
424  }
425  break;
426  case 1:
427  currentParticle.SetDefinitionAndUpdateE( anAntiProton );
428  incidentHasChanged = true;
429  break;
430  default:
431  targetParticle.SetDefinitionAndUpdateE( aProton );
432  targetHasChanged = true;
433  break;
434  }
435  }
436  } else { // random number <= anhl[iplab]
437  if (centerofmassEnergy <= 2*aPiPlus->GetPDGMass()/MeV ) {
438  quasiElastic = true;
439  return;
440  }
441 
442  // annihilation channels
443 
444  G4double n, anpn;
445  GetNormalizationConstant( -centerofmassEnergy, n, anpn );
446  G4double ran = G4UniformRand();
447  G4double dum, excs = 0.0;
448 
449  if (targetParticle.GetDefinition() == aProton ) {
450  counter = -1;
451  for (np=1; (np<numSec/3) && (ran>=excs); ++np ) {
452  nneg = np-1;
453  for (nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
454  if (++counter < numMulA) {
455  nt = np+nneg+nz;
456  if (nt>1 && nt<=numSec ) {
457  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
458  dum = (pi/anpn)*nt*protmulA[counter]*protnormA[nt-1]/(2.0*n*n);
459  if (std::fabs(dum) < 1.0 ) {
460  if (test >= 1.0e-10) excs += dum*test;
461  }
462  else
463  excs += dum*test;
464  }
465  }
466  }
467  }
468  } else { // target must be a neutron
469  counter = -1;
470  for (np=0; (np<numSec/3) && (ran>=excs); ++np ) {
471  nneg = np;
472  for (nz=0; (nz<numSec/3) && (ran>=excs); ++nz ) {
473  if (++counter < numMulA ) {
474  nt = np+nneg+nz;
475  if ((nt>1) && (nt<=numSec) ) {
476  test = G4Exp( std::min( expxu, std::max( expxl, -(pi/4.0)*(nt*nt)/(n*n) ) ) );
477  dum = (pi/anpn)*nt*neutmulA[counter]*neutnormA[nt-1]/(2.0*n*n);
478  if (std::fabs(dum) < 1.0 ) {
479  if( test >= 1.0e-10 )excs += dum*test;
480  }
481  else
482  excs += dum*test;
483  }
484  }
485  }
486  }
487  }
488 
489  if (ran >= excs) { // 3 previous loops continued to the end
490  quasiElastic = true;
491  return;
492  }
493  np--; nz--;
494  currentParticle.SetMass( 0.0 );
495  targetParticle.SetMass( 0.0 );
496  }
497 
498  G4int loop = 0;
499  G4ExceptionDescription ed;
500  ed << " While count exceeded " << G4endl;
501  while(np+nneg+nz<3) { /* Loop checking, 01.09.2015, D.Wright */
502  nz++;
503  loop++;
504  if (loop > 1000) {
505  G4Exception("G4RPGAntiNeutronInelastic::Cascade()", "HAD_RPG_100", JustWarning, ed);
506  break;
507  }
508  }
509 
510  SetUpPions( np, nneg, nz, vec, vecLen );
511  return;
512 }
513 
514  /* end of file */
515