ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4AntiNuclElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4AntiNuclElastic.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 // Geant4 Header : G4AntiNuclElastic
28 //
29 //
30 
31 #include "G4AntiNuclElastic.hh"
32 
33 #include "G4PhysicalConstants.hh"
34 #include "G4SystemOfUnits.hh"
35 #include "G4ParticleTable.hh"
36 #include "G4ParticleDefinition.hh"
37 #include "G4IonTable.hh"
38 #include "Randomize.hh"
39 #include "G4AntiProton.hh"
40 #include "G4AntiNeutron.hh"
41 #include "G4AntiDeuteron.hh"
42 #include "G4AntiAlpha.hh"
43 #include "G4AntiTriton.hh"
44 #include "G4AntiHe3.hh"
45 #include "G4Proton.hh"
46 #include "G4Neutron.hh"
47 #include "G4Deuteron.hh"
48 #include "G4Alpha.hh"
49 #include "G4Pow.hh"
50 #include "G4Exp.hh"
51 #include "G4Log.hh"
52 
53 #include "G4NucleiProperties.hh"
54 #include "G4CrossSectionDataSetRegistry.hh"
55 
56 G4AntiNuclElastic::G4AntiNuclElastic()
57  : G4HadronElastic("AntiAElastic")
58 {
59  //V.Ivanchenko commented out
60  //SetMinEnergy( 0.1*GeV );
61  //SetMaxEnergy( 10.*TeV );
62 
63  theAProton = G4AntiProton::AntiProton();
64  theANeutron = G4AntiNeutron::AntiNeutron();
65  theADeuteron = G4AntiDeuteron::AntiDeuteron();
66  theATriton = G4AntiTriton::AntiTriton();
67  theAAlpha = G4AntiAlpha::AntiAlpha();
68  theAHe3 = G4AntiHe3::AntiHe3();
69 
70  theProton = G4Proton::Proton();
71  theNeutron = G4Neutron::Neutron();
72  theDeuteron = G4Deuteron::Deuteron();
73  theAlpha = G4Alpha::Alpha();
74 
75  G4CrossSectionDataSetRegistry* reg = G4CrossSectionDataSetRegistry::Instance();
76  cs = static_cast<G4ComponentAntiNuclNuclearXS*>(reg->GetComponentCrossSection("AntiAGlauber"));
77  if(!cs) { cs = new G4ComponentAntiNuclNuclearXS(); }
78 
79  fParticle = 0;
80  fWaveVector = 0.;
81  fBeta = 0.;
82  fZommerfeld = 0.;
83  fAm = 0.;
84  fTetaCMS = 0.;
85  fRa = 0.;
86  fRef = 0.;
87  fceff = 0.;
88  fptot = 0.;
89  fTmax = 0.;
90  fThetaLab = 0.;
91 }
92 
94 G4AntiNuclElastic::~G4AntiNuclElastic()
95 {}
96 
98 // sample momentum transfer in the CMS system
99 G4double G4AntiNuclElastic::SampleInvariantT(const G4ParticleDefinition* particle,
100  G4double Plab, G4int Z, G4int A)
101 {
102  G4double T;
103  G4double Mproj = particle->GetPDGMass();
104  G4LorentzVector Pproj(0.,0.,Plab,std::sqrt(Plab*Plab+Mproj*Mproj));
105  G4double ctet1 = GetcosTeta1(Plab, A);
106 
107  G4double energy=Pproj.e()-Mproj;
108 
109  const G4ParticleDefinition* theParticle = particle;
110 
111  G4ParticleDefinition * theDef = 0;
112 
113  if(Z == 1 && A == 1) theDef = theProton;
114  else if (Z == 1 && A == 2) theDef = theDeuteron;
115  else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
116  else if (Z == 2 && A == 3) theDef = G4He3::He3();
117  else if (Z == 2 && A == 4) theDef = theAlpha;
118 
119 
120  G4double TargMass =G4NucleiProperties::GetNuclearMass(A,Z);
121 
122  //transform to CMS
123 
124  G4LorentzVector lv(0.0,0.0,0.0,TargMass);
125  lv += Pproj;
126  G4double S = lv.mag2()/(GeV*GeV);
127 
128  G4ThreeVector bst = lv.boostVector();
129  Pproj.boost(-bst);
130 
131  G4ThreeVector p1 = Pproj.vect();
132  G4double ptot = p1.mag();
133 
134  fbst = bst;
135  fptot= ptot;
136  fTmax = 4.0*ptot*ptot;
137 
138  if(Plab < (std::abs(particle->GetBaryonNumber())*100)*MeV) // Uzhi 24 Nov. 2011
139  {return fTmax*G4UniformRand();} // Uzhi 24 Nov. 2011
140 
141  G4double Z1 = particle->GetPDGCharge();
142  G4double Z2 = Z;
143 
144  G4double beta = CalculateParticleBeta(particle, ptot);
145  G4double n = CalculateZommerfeld( beta, Z1, Z2 );
146  G4double Am = CalculateAm( ptot, n, Z2 );
147  fWaveVector = ptot; // /hbarc;
148 
149  G4LorentzVector Fproj(0.,0.,0.,0.);
150  G4double XsCoulomb = sqr(n/fWaveVector)*pi*(1+ctet1)/(1.+Am)/(1.+2.*Am-ctet1);
151  XsCoulomb=XsCoulomb*0.38938e+6;
152 
153  G4double XsElastHad =cs->GetElasticElementCrossSection(particle, energy, Z, (G4double)A);
154  G4double XstotalHad =cs->GetTotalElementCrossSection(particle, energy, Z, (G4double)A);
155 
156  XsElastHad/=millibarn; XstotalHad/=millibarn;
157 
158  G4double CoulombProb = XsCoulomb/(XsCoulomb+XsElastHad);
159 
160 // G4cout<<" XselastHadron " << XsElastHad << " XsCol "<< XsCoulomb <<G4endl;
161 // G4cout <<" XsTotal" << XstotalHad <<G4endl;
162 // G4cout<<"XsInel"<< XstotalHad-XsElastHad<<G4endl;
163 
164  if(G4UniformRand() < CoulombProb)
165  { // Simulation of Coulomb scattering
166 
167  G4double phi = twopi * G4UniformRand();
168  G4double Ksi = G4UniformRand();
169 
170  G4double par1 = 2.*(1.+Am)/(1.+ctet1);
171 
172 // ////sample ThetaCMS in Coulomb part
173 
174  G4double cosThetaCMS = (par1*ctet1- Ksi*(1.+2.*Am))/(par1-Ksi);
175 
176  G4double PtZ=ptot*cosThetaCMS;
177  Fproj.setPz(PtZ);
178  G4double PtProjCMS = ptot*std::sqrt(1.0 - cosThetaCMS*cosThetaCMS);
179  G4double PtX= PtProjCMS * std::cos(phi);
180  G4double PtY= PtProjCMS * std::sin(phi);
181  Fproj.setPx(PtX);
182  Fproj.setPy(PtY);
183  Fproj.setE(std::sqrt(PtX*PtX+PtY*PtY+PtZ*PtZ+Mproj*Mproj));
184  T = -(Pproj-Fproj).mag2();
185  } else
186 
187  {
189 
190 // G4double Qmax = 2.*ptot*197.33; // in fm^-1
191  G4double Qmax = 2.*3.0*197.33; // in fm^-1
192  G4double Amag = 70*70; // A1 in Magora funct:A1*exp(-q*A2)
193  G4double SlopeMag = 2.*3.0; // A2 in Magora funct:A1*exp(-q*A2)
194 
195  G4double sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy);
196 
197  fRa = 1.113*G4Pow::GetInstance()->Z13(A) -
198  0.227/G4Pow::GetInstance()->Z13(A);
199  if(A == 3) fRa=1.81;
200  if(A == 4) fRa=1.37;
201 
202  if((A>=12.) && (A<27) ) fRa=fRa*0.85;
203  if((A>=27.) && (A<48) ) fRa=fRa*0.90;
204  if((A>=48.) && (A<65) ) fRa=fRa*0.95;
205 
206  G4double Ref2 = 0;
207  G4double ceff2 =0;
208  G4double rho = 0;
209  if ((theParticle == theAProton) || (theParticle == theANeutron))
210  {
211  if(theDef == theProton)
212  {
213 // G4double Mp2=sqr(theDef->GetPDGMass()/GeV );
214 
215 // change 30 October
216 
217  if(Plab < 610.)
218  { rho = 1.3347-10.342*Plab/1000.+22.277*Plab/1000.*Plab/1000.-
219  13.634*Plab/1000.*Plab/1000.*Plab/1000. ;}
220  if((Plab < 5500.)&&(Plab >= 610.) )
221  { rho = 0.22; }
222  if((Plab >= 5500.)&&(Plab < 12300.) )
223  { rho = -0.32; }
224  if( Plab >= 12300.)
225  { rho = 0.135-2.26/(std::sqrt(S)) ;}
226 
227  Ref2 = 0.35 + 0.9/std::sqrt(std::sqrt(S-4.*0.88))+0.04*G4Log(S) ;
228  ceff2 = 0.375 - 2./S + 0.44/(sqr(S-4.)+1.5) ;
229 
230 /*
231  Ref2=0.8/std::sqrt(std::sqrt(S-4.*Mp2)) + 0.55;
232  if(S>1000.) Ref2=0.62+0.02*G4Log(S) ;
233  ceff2 = 0.035/(sqr(S-4.3)+0.4) + 0.085 * G4Log(S) ;
234  if(S>1000.) ceff2 = 0.005 * G4Log(S) + 0.29;
235 */
236 
237  Ref2=Ref2*Ref2;
238  ceff2 = ceff2*ceff2;
239 
240  SlopeMag = 0.5; // Uzhi
241  Amag= 1.; // Uzhi
242  }
243 
244  if(Z>2)
245  { Ref2 = fRa*fRa +2.48*0.01*sig_pbarp*fRa - 2.23e-6*sig_pbarp*sig_pbarp*fRa*fRa;
246  ceff2 = 0.16+3.3e-4*sig_pbarp+0.35*G4Exp(-0.03*sig_pbarp);
247  }
248  if( (Z==2)&&(A==4) )
249  { Ref2 = fRa*fRa -0.46 +0.03*sig_pbarp - 2.98e-6*sig_pbarp*sig_pbarp;
250  ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*G4Exp(-0.03*sig_pbarp);
251  }
252  if( (Z==1)&&(A==3) )
253  { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
254  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
255  }
256  if( (Z==2)&&(A==3) )
257  { Ref2 = fRa*fRa - 1.36 + 0.025 * sig_pbarp - 3.69e-7 * sig_pbarp*sig_pbarp;
258  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
259  }
260  if( (Z==1)&&(A==2) )
261  {
262  Ref2 = fRa*fRa - 0.28 + 0.019 * sig_pbarp + 2.06e-6 * sig_pbarp*sig_pbarp;
263  ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*G4Exp(-0.03*sig_pbarp);
264  }
265  }
266 
267  if (theParticle == theADeuteron)
268  {
269  sig_pbarp= cs->GetAntiHadronNucleonTotCrSc(particle,energy/2.);
270  Ref2 = XstotalHad/10./2./pi ;
271  if(Z>2)
272  {
273  ceff2 = 0.38 + 2.0e-4 *sig_pbarp + 0.5 * G4Exp(-0.03*sig_pbarp);
274  }
275  if(theDef == theProton)
276  {
277  ceff2 = 0.297 + 7.853e-04*sig_pbarp + 0.2899*G4Exp(-0.03*sig_pbarp);
278  }
279  if(theDef == theDeuteron)
280  {
281  ceff2 = 0.65 + 3.0e-4*sig_pbarp + 0.55 * G4Exp(-0.03*sig_pbarp);
282  }
283  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
284  {
285  ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * G4Exp(-0.02*sig_pbarp);
286  }
287  if(theDef == theAlpha)
288  {
289  ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * G4Exp(-0.02*sig_pbarp);
290  }
291  }
292 
293  if( (theParticle ==theAHe3) || (theParticle ==theATriton) )
294  {
295  sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
296  Ref2 = XstotalHad/10./2./pi ;
297  if(Z>2)
298  {
299  ceff2 = 0.26 + 2.2e-4*sig_pbarp + 0.33*G4Exp(-0.03*sig_pbarp);
300  }
301  if(theDef == theProton)
302  {
303  ceff2 = 0.149 + 7.091e-04*sig_pbarp + 0.3743*G4Exp(-0.03*sig_pbarp);
304  }
305  if(theDef == theDeuteron)
306  {
307  ceff2 = 0.57 + 2.5e-4*sig_pbarp + 0.65 * G4Exp(-0.02*sig_pbarp);
308  }
309  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
310  {
311  ceff2 = 0.39 + 2.7e-4*sig_pbarp + 0.7 * G4Exp(-0.02*sig_pbarp);
312  }
313  if(theDef == theAlpha)
314  {
315  ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * G4Exp(-0.03*sig_pbarp);
316  }
317  }
318 
319 
320  if (theParticle == theAAlpha)
321  {
322  sig_pbarp = cs->GetAntiHadronNucleonTotCrSc(particle,energy/3.);
323  Ref2 = XstotalHad/10./2./pi ;
324  if(Z>2)
325  {
326  ceff2 = 0.22 + 2.0e-4*sig_pbarp + 0.2 * G4Exp(-0.03*sig_pbarp);
327  }
328  if(theDef == theProton)
329  {
330  ceff2= 0.078 + 6.657e-4*sig_pbarp + 0.3359*G4Exp(-0.03*sig_pbarp);
331  }
332  if(theDef == theDeuteron)
333  {
334  ceff2 = 0.40 + 3.5e-4 *sig_pbarp + 0.45 * G4Exp(-0.02*sig_pbarp);
335  }
336  if( (theDef == G4Triton::Triton()) || (theDef == G4He3::He3() ) )
337  {
338  ceff2 = 0.24 + 3.5e-4*sig_pbarp + 0.75 * G4Exp(-0.03*sig_pbarp);
339  }
340  if(theDef == theAlpha)
341  {
342  ceff2 = 0.17 + 3.5e-4*sig_pbarp + 0.45 * G4Exp(-0.03*sig_pbarp);
343  }
344  }
345 
346  fRef=std::sqrt(Ref2);
347  fceff = std::sqrt(ceff2);
348 // G4cout<<" Ref "<<fRef<<" c_eff "<<fceff<< " rho "<< rho<<G4endl;
349 
350 
351  G4double Q = 0.0 ;
352  G4double BracFunct;
353  const G4int maxNumberOfLoops = 10000;
354  G4int loopCounter = 0;
355  do
356  {
357  Q = -G4Log(1.-(1.- G4Exp(-SlopeMag * Qmax))* G4UniformRand() )/SlopeMag;
358  G4double x = fRef * Q;
359  BracFunct = ( ( sqr(BesselOneByArg(x))+sqr(rho/2. * BesselJzero(x)) )
360 * sqr(DampFactor(pi*fceff*Q))) /(Amag*G4Exp(-SlopeMag*Q));
361 
362  BracFunct = BracFunct * Q * sqr(sqr(fRef));
363  }
364  while ( (G4UniformRand()>BracFunct) &&
365  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
366  if ( loopCounter >= maxNumberOfLoops ) {
367  fTetaCMS = 0.0;
368  return 0.0;
369  }
370 
371  T= sqr(Q);
372  T*=3.893913e+4; // fm -> MeV^2
373  }
374 
375  // VI: 29.04.2019 unnecessary computation of trigonometry
376  /*
377  G4double cosTet=1.0-T/(2.*ptot*ptot);
378  if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
379  if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
380  fTetaCMS=std::acos(cosTet);
381  */
382  return T;
383 }
384 
386 // Sample of Theta in CMS
387  G4double G4AntiNuclElastic::SampleThetaCMS(const G4ParticleDefinition* p, G4double plab,
388  G4int Z, G4int A)
389 {
390  G4double T;
391  T = SampleInvariantT( p, plab, Z, A);
392 
393  // NaN finder
394  if(!(T < 0.0 || T >= 0.0))
395  {
396  if (verboseLevel > 0)
397  {
398  G4cout << "G4DiffuseElastic:WARNING: A = " << A
399  << " mom(GeV)= " << plab/GeV
400  << " S-wave will be sampled"
401  << G4endl;
402  }
403  T = G4UniformRand()*fTmax;
404 
405  }
406 
407  if(fptot > 0.) // Uzhi 24 Nov. 2011
408  {
409  G4double cosTet=1.0-T/(2.*fptot*fptot);
410  if(cosTet > 1.0 ) cosTet= 1.; // Uzhi 30 Nov.
411  if(cosTet < -1.0 ) cosTet=-1.; // Uzhi 30 Nov.
412  fTetaCMS=std::acos(cosTet);
413  return fTetaCMS;
414  } else // Uzhi 24 Nov. 2011
415  { // Uzhi 24 Nov. 2011
416  return 2.*G4UniformRand()-1.; // Uzhi 24 Nov. 2011
417  } // Uzhi 24 Nov. 2011
418 }
419 
420 
422 // Sample of Theta in Lab System
423  G4double G4AntiNuclElastic::SampleThetaLab(const G4ParticleDefinition* p, G4double plab,
424  G4int Z, G4int A)
425 {
426  G4double T;
427  T = SampleInvariantT( p, plab, Z, A);
428 
429  // NaN finder
430  if(!(T < 0.0 || T >= 0.0))
431  {
432  if (verboseLevel > 0)
433  {
434  G4cout << "G4DiffuseElastic:WARNING: A = " << A
435  << " mom(GeV)= " << plab/GeV
436  << " S-wave will be sampled"
437  << G4endl;
438  }
439  T = G4UniformRand()*fTmax;
440  }
441 
442  G4double phi = G4UniformRand()*twopi;
443 
444  G4double cost(1.);
445  if(fTmax > 0.) {cost = 1. - 2.0*T/fTmax;} // Uzhi 24 Nov. 2011
446 
447  G4double sint;
448  if( cost >= 1.0 )
449  {
450  cost = 1.0;
451  sint = 0.0;
452  }
453  else if( cost <= -1.0)
454  {
455  cost = -1.0;
456  sint = 0.0;
457  }
458  else
459  {
460  sint = std::sqrt((1.0-cost)*(1.0+cost));
461  }
462 
463  G4double m1 = p->GetPDGMass();
464  G4ThreeVector v(sint*std::cos(phi),sint*std::sin(phi),cost);
465  v *= fptot;
466  G4LorentzVector nlv(v.x(),v.y(),v.z(),std::sqrt(fptot*fptot + m1*m1));
467 
468  nlv.boost(fbst);
469 
470  G4ThreeVector np = nlv.vect();
471  G4double theta = np.theta();
472  fThetaLab = theta;
473 
474  return theta;
475 }
476 
478 // Calculation of Damp factor
479  G4double G4AntiNuclElastic::DampFactor(G4double x)
480 {
481  G4double df;
482  G4double f3 = 6.; // first factorials
483 
484  if( std::fabs(x) < 0.01 )
485  {
486  df=1./(1.+x*x/f3);
487  }
488  else
489  {
490  df = x/std::sinh(x);
491  }
492  return df;
493 }
494 
495 
497 // Calculation of particle velocity Beta
498 
499  G4double G4AntiNuclElastic::CalculateParticleBeta( const G4ParticleDefinition* particle,
500  G4double momentum )
501 {
502  G4double mass = particle->GetPDGMass();
503  G4double a = momentum/mass;
504  fBeta = a/std::sqrt(1+a*a);
505 
506  return fBeta;
507 }
508 
509 
511 // Calculation of parameter Zommerfeld
512 
513  G4double G4AntiNuclElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
514 {
515  fZommerfeld = fine_structure_const*Z1*Z2/beta;
516 
517  return fZommerfeld;
518 }
519 
521 //
522 G4double G4AntiNuclElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
523 {
524  G4double k = momentum/hbarc;
525  G4double ch = 1.13 + 3.76*n*n;
526  G4double zn = 1.77*k/G4Pow::GetInstance()->A13(Z)*Bohr_radius;
527  G4double zn2 = zn*zn;
528  fAm = ch/zn2;
529 
530  return fAm;
531 }
532 
534 //
535 // Bessel J0 function based on rational approximation from
536 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
537 
538 G4double G4AntiNuclElastic::BesselJzero(G4double value)
539 {
540  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
541 
542  modvalue = std::fabs(value);
543 
544  if ( value < 8.0 && value > -8.0 )
545  {
546  value2 = value*value;
547 
548  fact1 = 57568490574.0 + value2*(-13362590354.0
549  + value2*( 651619640.7
550  + value2*(-11214424.18
551  + value2*( 77392.33017
552  + value2*(-184.9052456 ) ) ) ) );
553 
554  fact2 = 57568490411.0 + value2*( 1029532985.0
555  + value2*( 9494680.718
556  + value2*(59272.64853
557  + value2*(267.8532712
558  + value2*1.0 ) ) ) );
559 
560  bessel = fact1/fact2;
561  }
562  else
563  {
564  arg = 8.0/modvalue;
565 
566  value2 = arg*arg;
567 
568  shift = modvalue-0.785398164;
569 
570  fact1 = 1.0 + value2*(-0.1098628627e-2
571  + value2*(0.2734510407e-4
572  + value2*(-0.2073370639e-5
573  + value2*0.2093887211e-6 ) ) );
574  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
575  + value2*(-0.6911147651e-5
576  + value2*(0.7621095161e-6
577  - value2*0.934945152e-7 ) ) );
578 
579  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
580  }
581  return bessel;
582 }
583 
584 
586 // Bessel J1 function based on rational approximation from
587 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
588 
589  G4double G4AntiNuclElastic::BesselJone(G4double value)
590 {
591  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
592 
593  modvalue = std::fabs(value);
594 
595  if ( modvalue < 8.0 )
596  {
597  value2 = value*value;
598  fact1 = value*(72362614232.0 + value2*(-7895059235.0
599  + value2*( 242396853.1
600  + value2*(-2972611.439
601  + value2*( 15704.48260
602  + value2*(-30.16036606 ) ) ) ) ) );
603 
604  fact2 = 144725228442.0 + value2*(2300535178.0
605  + value2*(18583304.74
606  + value2*(99447.43394
607  + value2*(376.9991397
608  + value2*1.0 ) ) ) );
609  bessel = fact1/fact2;
610  }
611  else
612  {
613  arg = 8.0/modvalue;
614  value2 = arg*arg;
615 
616  shift = modvalue - 2.356194491;
617 
618  fact1 = 1.0 + value2*( 0.183105e-2
619  + value2*(-0.3516396496e-4
620  + value2*(0.2457520174e-5
621  + value2*(-0.240337019e-6 ) ) ) );
622 
623  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
624  + value2*( 0.8449199096e-5
625  + value2*(-0.88228987e-6
626  + value2*0.105787412e-6 ) ) );
627 
628  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
629  if (value < 0.0) bessel = -bessel;
630  }
631  return bessel;
632 }
633 
635 // return J1(x)/x with special case for small x
636 G4double G4AntiNuclElastic::BesselOneByArg(G4double x)
637 {
638  G4double x2, result;
639 
640  if( std::fabs(x) < 0.01 )
641  {
642  x *= 0.5;
643  x2 = x*x;
644  result = (2.- x2 + x2*x2/6.)/4.;
645  }
646  else
647  {
648  result = BesselJone(x)/x;
649  }
650  return result;
651 }
652 
654 // return angle from which Coulomb scattering is calculated
655 G4double G4AntiNuclElastic::GetcosTeta1(G4double plab, G4int A)
656 {
657 
658 // G4double p0 =G4LossTableManager::Instance()->FactorForAngleLimit()*CLHEP::hbarc/CLHEP::fermi;
659  G4double p0 = 1.*hbarc/fermi;
660 //G4double cteta1 = 1.0 - p0*p0/2.0 * pow(A,2./3.)/(plab*plab);
661  G4double cteta1 = 1.0 - p0*p0/2.0 * G4Pow::GetInstance()->Z23(A)/(plab*plab);
663  if(cteta1 < -1.) cteta1 = -1.0;
664  return cteta1;
665 }
666 
667 
668 
669 
670 
671 
672