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