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