ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclNuclDiffuseElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclNuclDiffuseElastic.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 // Physics model class G4NuclNuclDiffuseElastic
29 //
30 //
31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
32 //
33 // 24-May-07 V. Grichine
34 //
35 
37 #include "G4ParticleTable.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4IonTable.hh"
40 #include "G4NucleiProperties.hh"
41 
42 #include "Randomize.hh"
43 #include "G4Integrator.hh"
44 #include "globals.hh"
45 #include "G4PhysicalConstants.hh"
46 #include "G4SystemOfUnits.hh"
47 
48 #include "G4Proton.hh"
49 #include "G4Neutron.hh"
50 #include "G4Deuteron.hh"
51 #include "G4Alpha.hh"
52 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh"
54 
55 #include "G4Element.hh"
56 #include "G4ElementTable.hh"
57 #include "G4NistManager.hh"
58 #include "G4PhysicsTable.hh"
59 #include "G4PhysicsLogVector.hh"
60 #include "G4PhysicsFreeVector.hh"
61 #include "G4HadronicParameters.hh"
62 
64 //
65 // Test Constructor. Just to check xsc
66 
67 
69  : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
70 {
71  SetMinEnergy( 50*MeV );
73  verboseLevel = 0;
74  lowEnergyRecoilLimit = 100.*keV;
75  lowEnergyLimitQ = 0.0*GeV;
76  lowEnergyLimitHE = 0.0*GeV;
77  lowestEnergyLimit= 0.0*keV;
78  plabLowLimit = 20.0*MeV;
79 
86 
87  fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV
88  fAngleBin = 200;
89 
91  fAngleTable = 0;
92 
93  fParticle = 0;
94  fWaveVector = 0.;
95  fAtomicWeight = 0.;
96  fAtomicNumber = 0.;
97  fNuclearRadius = 0.;
98  fBeta = 0.;
99  fZommerfeld = 0.;
100  fAm = 0.;
101  fAddCoulomb = false;
102  // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
103 
104  // Empirical parameters
105 
106  fCofAlphaMax = 1.5;
107  fCofAlphaCoulomb = 0.5;
108 
109  fProfileDelta = 1.;
110  fProfileAlpha = 0.5;
111 
112  fCofLambda = 1.0;
113  fCofDelta = 0.04;
114  fCofAlpha = 0.095;
115 
119  = fSumSigma = fEtaRatio = fReZ = 0.0;
120  fMaxL = 0;
121 
122  fNuclearRadiusCof = 1.0;
123  fCoulombMuC = 0.0;
124 }
125 
127 //
128 // Destructor
129 
131 {
132  if ( fEnergyVector ) {
133  delete fEnergyVector;
134  fEnergyVector = 0;
135  }
136 
137  for ( std::vector<G4PhysicsTable*>::iterator it = fAngleBank.begin();
138  it != fAngleBank.end(); ++it ) {
139  if ( (*it) ) (*it)->clearAndDestroy();
140  delete *it;
141  *it = 0;
142  }
143  fAngleTable = 0;
144 }
145 
147 //
148 // Initialisation for given particle using element table of application
149 
151 {
152 
153  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
154 
155  const G4ElementTable* theElementTable = G4Element::GetElementTable();
156  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
157 
158  // projectile radius
159 
161  G4double R1 = CalculateNuclearRad(A1);
162 
163  for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
164  {
165  fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
167 
169  fNuclearRadius += R1;
170 
171  if(verboseLevel > 0)
172  {
173  G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
174  <<(*theElementTable)[jEl]->GetName()<<G4endl;
175  }
177  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
178 
179  BuildAngleTable();
180  fAngleBank.push_back(fAngleTable);
181  }
182 }
183 
184 
186 //
187 // return differential elastic cross section d(sigma)/d(omega)
188 
189 G4double
191  G4double theta,
193  G4double A )
194 {
196  fWaveVector = momentum/hbarc;
197  fAtomicWeight = A;
198  fAddCoulomb = false;
200 
202 
203  return sigma;
204 }
205 
207 //
208 // return invariant differential elastic cross section d(sigma)/d(tMand)
209 
210 G4double
212  G4double tMand,
213  G4double plab,
214  G4double A, G4double Z )
215 {
216  G4double m1 = particle->GetPDGMass();
217  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
218 
219  G4int iZ = static_cast<G4int>(Z+0.5);
220  G4int iA = static_cast<G4int>(A+0.5);
221  G4ParticleDefinition * theDef = 0;
222 
223  if (iZ == 1 && iA == 1) theDef = theProton;
224  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
225  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
226  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
227  else if (iZ == 2 && iA == 4) theDef = theAlpha;
228  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
229 
230  G4double tmass = theDef->GetPDGMass();
231 
232  G4LorentzVector lv(0.0,0.0,0.0,tmass);
233  lv += lv1;
234 
235  G4ThreeVector bst = lv.boostVector();
236  lv1.boost(-bst);
237 
238  G4ThreeVector p1 = lv1.vect();
239  G4double ptot = p1.mag();
240  G4double ptot2 = ptot*ptot;
241  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
242 
243  if( cost >= 1.0 ) cost = 1.0;
244  else if( cost <= -1.0) cost = -1.0;
245 
246  G4double thetaCMS = std::acos(cost);
247 
248  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
249 
250  sigma *= pi/ptot2;
251 
252  return sigma;
253 }
254 
256 //
257 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
258 // correction
259 
260 G4double
262  G4double theta,
264  G4double A, G4double Z )
265 {
267  fWaveVector = momentum/hbarc;
268  fAtomicWeight = A;
269  fAtomicNumber = Z;
271  fAddCoulomb = false;
272 
273  G4double z = particle->GetPDGCharge();
274 
276  G4double kRtC = 1.9;
277 
278  if( z && (kRt > kRtC) )
279  {
280  fAddCoulomb = true;
281  fBeta = CalculateParticleBeta( particle, momentum);
283  fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
284  }
286 
287  return sigma;
288 }
289 
291 //
292 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
293 // correction
294 
295 G4double
297  G4double tMand,
298  G4double plab,
299  G4double A, G4double Z )
300 {
301  G4double m1 = particle->GetPDGMass();
302 
303  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
304 
305  G4int iZ = static_cast<G4int>(Z+0.5);
306  G4int iA = static_cast<G4int>(A+0.5);
307 
308  G4ParticleDefinition* theDef = 0;
309 
310  if (iZ == 1 && iA == 1) theDef = theProton;
311  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
312  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
313  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
314  else if (iZ == 2 && iA == 4) theDef = theAlpha;
315  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
316 
317  G4double tmass = theDef->GetPDGMass();
318 
319  G4LorentzVector lv(0.0,0.0,0.0,tmass);
320  lv += lv1;
321 
322  G4ThreeVector bst = lv.boostVector();
323  lv1.boost(-bst);
324 
325  G4ThreeVector p1 = lv1.vect();
326  G4double ptot = p1.mag();
327  G4double ptot2 = ptot*ptot;
328  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
329 
330  if( cost >= 1.0 ) cost = 1.0;
331  else if( cost <= -1.0) cost = -1.0;
332 
333  G4double thetaCMS = std::acos(cost);
334 
335  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
336 
337  sigma *= pi/ptot2;
338 
339  return sigma;
340 }
341 
343 //
344 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
345 // correction
346 
347 G4double
349  G4double tMand,
350  G4double plab,
351  G4double A, G4double Z )
352 {
353  G4double m1 = particle->GetPDGMass();
354  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
355 
356  G4int iZ = static_cast<G4int>(Z+0.5);
357  G4int iA = static_cast<G4int>(A+0.5);
358  G4ParticleDefinition * theDef = 0;
359 
360  if (iZ == 1 && iA == 1) theDef = theProton;
361  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
362  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
363  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
364  else if (iZ == 2 && iA == 4) theDef = theAlpha;
365  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
366 
367  G4double tmass = theDef->GetPDGMass();
368 
369  G4LorentzVector lv(0.0,0.0,0.0,tmass);
370  lv += lv1;
371 
372  G4ThreeVector bst = lv.boostVector();
373  lv1.boost(-bst);
374 
375  G4ThreeVector p1 = lv1.vect();
376  G4double ptot = p1.mag();
377  G4double ptot2 = ptot*ptot;
378  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
379 
380  if( cost >= 1.0 ) cost = 1.0;
381  else if( cost <= -1.0) cost = -1.0;
382 
383  G4double thetaCMS = std::acos(cost);
384 
385  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
386 
387  sigma *= pi/ptot2;
388 
389  return sigma;
390 }
391 
393 //
394 // return differential elastic probability d(probability)/d(omega)
395 
396 G4double
397 G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
398  G4double theta
399  // G4double momentum,
400  // G4double A
401  )
402 {
403  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
404  G4double delta, diffuse, gamma;
405  G4double e1, e2, bone, bone2;
406 
407  // G4double wavek = momentum/hbarc; // wave vector
408  // G4double r0 = 1.08*fermi;
409  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
410 
411  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
412  G4double kr2 = kr*kr;
413  G4double krt = kr*theta;
414 
415  bzero = BesselJzero(krt);
416  bzero2 = bzero*bzero;
417  bone = BesselJone(krt);
418  bone2 = bone*bone;
419  bonebyarg = BesselOneByArg(krt);
420  bonebyarg2 = bonebyarg*bonebyarg;
421 
422  // VI - Coverity complains
423  /*
424  if (fParticle == theProton)
425  {
426  diffuse = 0.63*fermi;
427  gamma = 0.3*fermi;
428  delta = 0.1*fermi*fermi;
429  e1 = 0.3*fermi;
430  e2 = 0.35*fermi;
431  }
432  else // as proton, if were not defined
433  {
434  */
435  diffuse = 0.63*fermi;
436  gamma = 0.3*fermi;
437  delta = 0.1*fermi*fermi;
438  e1 = 0.3*fermi;
439  e2 = 0.35*fermi;
440  //}
441  G4double lambda = 15.; // 15 ok
442 
443  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
444 
445  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
446  G4double kgamma2 = kgamma*kgamma;
447 
448  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
449  // G4double dk2t2 = dk2t*dk2t;
450  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
451 
452  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
453 
454  damp = DampFactor(pikdt);
455  damp2 = damp*damp;
456 
457  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
458  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
459 
460 
461  sigma = kgamma2;
462  // sigma += dk2t2;
463  sigma *= bzero2;
464  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465  sigma += kr2*bonebyarg2;
466  sigma *= damp2; // *rad*rad;
467 
468  return sigma;
469 }
470 
472 //
473 // return differential elastic probability d(probability)/d(omega) with
474 // Coulomb correction
475 
476 G4double
477 G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
478  G4double theta
479  // G4double momentum,
480  // G4double A
481  )
482 {
483  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
484  G4double delta, diffuse, gamma;
485  G4double e1, e2, bone, bone2;
486 
487  // G4double wavek = momentum/hbarc; // wave vector
488  // G4double r0 = 1.08*fermi;
489  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
490 
491  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
492  G4double kr2 = kr*kr;
493  G4double krt = kr*theta;
494 
495  bzero = BesselJzero(krt);
496  bzero2 = bzero*bzero;
497  bone = BesselJone(krt);
498  bone2 = bone*bone;
499  bonebyarg = BesselOneByArg(krt);
500  bonebyarg2 = bonebyarg*bonebyarg;
501 
502  if (fParticle == theProton)
503  {
504  diffuse = 0.63*fermi;
505  // diffuse = 0.6*fermi;
506  gamma = 0.3*fermi;
507  delta = 0.1*fermi*fermi;
508  e1 = 0.3*fermi;
509  e2 = 0.35*fermi;
510  }
511  else // as proton, if were not defined
512  {
513  diffuse = 0.63*fermi;
514  gamma = 0.3*fermi;
515  delta = 0.1*fermi*fermi;
516  e1 = 0.3*fermi;
517  e2 = 0.35*fermi;
518  }
519  G4double lambda = 15.; // 15 ok
520  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
521  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
522 
523  // G4cout<<"kgamma = "<<kgamma<<G4endl;
524 
525  if(fAddCoulomb) // add Coulomb correction
526  {
527  G4double sinHalfTheta = std::sin(0.5*theta);
528  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
529 
530  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
531  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
532  }
533 
534  G4double kgamma2 = kgamma*kgamma;
535 
536  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
537  // G4cout<<"dk2t = "<<dk2t<<G4endl;
538  // G4double dk2t2 = dk2t*dk2t;
539  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
540 
541  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
542 
543  // G4cout<<"pikdt = "<<pikdt<<G4endl;
544 
545  damp = DampFactor(pikdt);
546  damp2 = damp*damp;
547 
548  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
549  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
550 
551  sigma = kgamma2;
552  // sigma += dk2t2;
553  sigma *= bzero2;
554  sigma += mode2k2*bone2;
555  sigma += e2dk3t*bzero*bone;
556 
557  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
558  sigma += kr2*bonebyarg2; // correction at J1()/()
559 
560  sigma *= damp2; // *rad*rad;
561 
562  return sigma;
563 }
564 
565 
567 //
568 // return differential elastic probability d(probability)/d(t) with
569 // Coulomb correction
570 
571 G4double
573 {
574  G4double theta;
575 
576  theta = std::sqrt(alpha);
577 
578  // theta = std::acos( 1 - alpha/2. );
579 
580  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
581  G4double delta, diffuse, gamma;
582  G4double e1, e2, bone, bone2;
583 
584  // G4double wavek = momentum/hbarc; // wave vector
585  // G4double r0 = 1.08*fermi;
586  // G4double rad = r0*G4Pow::GetInstance()->A13(A);
587 
588  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
589  G4double kr2 = kr*kr;
590  G4double krt = kr*theta;
591 
592  bzero = BesselJzero(krt);
593  bzero2 = bzero*bzero;
594  bone = BesselJone(krt);
595  bone2 = bone*bone;
596  bonebyarg = BesselOneByArg(krt);
597  bonebyarg2 = bonebyarg*bonebyarg;
598 
599  if (fParticle == theProton)
600  {
601  diffuse = 0.63*fermi;
602  // diffuse = 0.6*fermi;
603  gamma = 0.3*fermi;
604  delta = 0.1*fermi*fermi;
605  e1 = 0.3*fermi;
606  e2 = 0.35*fermi;
607  }
608  else // as proton, if were not defined
609  {
610  diffuse = 0.63*fermi;
611  gamma = 0.3*fermi;
612  delta = 0.1*fermi*fermi;
613  e1 = 0.3*fermi;
614  e2 = 0.35*fermi;
615  }
616  G4double lambda = 15.; // 15 ok
617  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
618  G4double kgamma = lambda*(1.-G4Exp(-fWaveVector*gamma/lambda)); // wavek*delta;
619 
620  // G4cout<<"kgamma = "<<kgamma<<G4endl;
621 
622  if(fAddCoulomb) // add Coulomb correction
623  {
624  G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
625  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
626 
627  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
628  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
629  }
630 
631  G4double kgamma2 = kgamma*kgamma;
632 
633  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
634  // G4cout<<"dk2t = "<<dk2t<<G4endl;
635  // G4double dk2t2 = dk2t*dk2t;
636  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
637 
638  G4double pikdt = lambda*(1.-G4Exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
639 
640  // G4cout<<"pikdt = "<<pikdt<<G4endl;
641 
642  damp = DampFactor(pikdt);
643  damp2 = damp*damp;
644 
645  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
646  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
647 
648  sigma = kgamma2;
649  // sigma += dk2t2;
650  sigma *= bzero2;
651  sigma += mode2k2*bone2;
652  sigma += e2dk3t*bzero*bone;
653 
654  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
655  sigma += kr2*bonebyarg2; // correction at J1()/()
656 
657  sigma *= damp2; // *rad*rad;
658 
659  return sigma;
660 }
661 
662 
664 //
665 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
666 
667 G4double
669 {
670  G4double result;
671 
672  result = GetDiffElasticSumProbA(alpha);
673 
674  // result *= 2*pi*std::sin(theta);
675 
676  return result;
677 }
678 
680 //
681 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
682 
683 G4double
685  G4double theta,
687  G4double A )
688 {
689  G4double result;
691  fWaveVector = momentum/hbarc;
692  fAtomicWeight = A;
693 
695 
696 
698 
699  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
700  result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
701 
702  return result;
703 }
704 
706 //
707 // Return inv momentum transfer -t > 0
708 
710  G4double p, G4double A)
711 {
712  G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
713  G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
714  return t;
715 }
716 
718 //
719 // Return scattering angle sampled in cms
720 
721 
722 G4double
725 {
726  G4int i, iMax = 100;
727  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
728 
730  fWaveVector = momentum/hbarc;
731  fAtomicWeight = A;
732 
734 
735  thetaMax = 10.174/fWaveVector/fNuclearRadius;
736 
737  if (thetaMax > pi) thetaMax = pi;
738 
740 
741  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
742  norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
743 
744  norm *= G4UniformRand();
745 
746  for(i = 1; i <= iMax; i++)
747  {
748  theta1 = (i-1)*thetaMax/iMax;
749  theta2 = i*thetaMax/iMax;
750  sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
751 
752  if ( sum >= norm )
753  {
754  result = 0.5*(theta1 + theta2);
755  break;
756  }
757  }
758  if (i > iMax ) result = 0.5*(theta1 + theta2);
759 
760  G4double sigma = pi*thetaMax/iMax;
761 
762  result += G4RandGauss::shoot(0.,sigma);
763 
764  if(result < 0.) result = 0.;
765  if(result > thetaMax) result = thetaMax;
766 
767  return result;
768 }
769 
772 
774 //
775 // Interface function. Return inv momentum transfer -t > 0 from initialisation table
776 
778  G4int Z, G4int A)
779 {
780  fParticle = aParticle;
781  fAtomicNumber = G4double(Z);
782  fAtomicWeight = G4double(A);
783 
784  G4double t(0.), m1 = fParticle->GetPDGMass();
785  G4double totElab = std::sqrt(m1*m1+p*p);
787  G4LorentzVector lv1(p,0.0,0.0,totElab);
788  G4LorentzVector lv(0.0,0.0,0.0,mass2);
789  lv += lv1;
790 
791  G4ThreeVector bst = lv.boostVector();
792  lv1.boost(-bst);
793 
794  G4ThreeVector p1 = lv1.vect();
795  G4double momentumCMS = p1.mag();
796 
797  // t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
798 
799  t = SampleCoulombMuCMS( aParticle, momentumCMS); // sample theta2 in cms
800 
801 
802  return t;
803 }
804 
806 //
807 // Return inv momentum transfer -t > 0 as Coulomb scattering <= thetaC
808 
810 {
811  G4double t(0.), rand(0.), mu(0.);
812 
813  G4double A1 = G4double( aParticle->GetBaryonNumber() );
814  G4double R1 = CalculateNuclearRad(A1);
815 
817  fNuclearRadius += R1;
818 
820 
822 
823  rand = G4UniformRand();
824 
825  // sample (1-cosTheta) in 0 < Theta < ThetaC as Coulomb scattering
826 
827  mu = fCoulombMuC*rand*fAm;
828  mu /= fAm + fCoulombMuC*( 1. - rand );
829 
830  t = 4.*p*p*mu;
831 
832  return t;
833 }
834 
835 
837 //
838 // Return inv momentum transfer -t > 0 from initialisation table
839 
841  G4double Z, G4double A)
842 {
843  G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
844  // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
845  G4double t = p*p*alpha; // -t !!!
846  return t;
847 }
848 
850 //
851 // Return scattering angle2 sampled in cms according to precalculated table.
852 
853 
854 G4double
857 {
858  size_t iElement;
859  G4int iMomentum, iAngle;
860  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
861  G4double m1 = particle->GetPDGMass();
862 
863  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
864  {
865  if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
866  }
867  if ( iElement == fElementNumberVector.size() )
868  {
869  InitialiseOnFly(Z,A); // table preparation, if needed
870 
871  // iElement--;
872 
873  // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
874  // << " is not found, return zero angle" << G4endl;
875  // return 0.; // no table for this element
876  }
877  // G4cout<<"iElement = "<<iElement<<G4endl;
878 
879  fAngleTable = fAngleBank[iElement];
880 
881  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882 
883  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
884  {
885  // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
886 
887  if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
888  }
889  // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
890 
891 
892  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
893  if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
894 
895 
896  if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
897  {
898  position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
899 
900  // G4cout<<"position = "<<position<<G4endl;
901 
902  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
903  {
904  if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
905  }
906  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
907 
908  // G4cout<<"iAngle = "<<iAngle<<G4endl;
909 
910  randAngle = GetScatteringAngle(iMomentum, iAngle, position);
911 
912  // G4cout<<"randAngle = "<<randAngle<<G4endl;
913  }
914  else // kinE inside between energy table edges
915  {
916  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
917  position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
918 
919  // G4cout<<"position = "<<position<<G4endl;
920 
921  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
922  {
923  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
924  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
925  }
926  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
927 
928  // G4cout<<"iAngle = "<<iAngle<<G4endl;
929 
930  theta2 = GetScatteringAngle(iMomentum, iAngle, position);
931 
932  // G4cout<<"theta2 = "<<theta2<<G4endl;
933 
934  E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
935 
936  // G4cout<<"E2 = "<<E2<<G4endl;
937 
938  iMomentum--;
939 
940  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
941 
942  // G4cout<<"position = "<<position<<G4endl;
943 
944  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
945  {
946  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
947  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
948  }
949  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
950 
951  theta1 = GetScatteringAngle(iMomentum, iAngle, position);
952 
953  // G4cout<<"theta1 = "<<theta1<<G4endl;
954 
955  E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
956 
957  // G4cout<<"E1 = "<<E1<<G4endl;
958 
959  W = 1.0/(E2 - E1);
960  W1 = (E2 - kinE)*W;
961  W2 = (kinE - E1)*W;
962 
963  randAngle = W1*theta1 + W2*theta2;
964 
965  // randAngle = theta2;
966  }
967  // G4double angle = randAngle;
968  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
969  // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
970 
971  return randAngle;
972 }
973 
975 //
976 // Initialisation for given particle on fly using new element number
977 
979 {
980  fAtomicNumber = Z; // atomic number
981  fAtomicWeight = G4NistManager::Instance()->GetAtomicMassAmu( static_cast< G4int >( Z ) );
982 
984  G4double R1 = CalculateNuclearRad(A1);
985 
987 
988  if( verboseLevel > 0 )
989  {
990  G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
991  <<Z<<"; and A = "<<A<<G4endl;
992  }
994 
995  BuildAngleTable();
996 
997  fAngleBank.push_back(fAngleTable);
998 
999  return;
1000 }
1001 
1003 //
1004 // Build for given particle and element table of momentum, angle probability.
1005 // For the moment in lab system.
1006 
1008 {
1009  G4int i, j;
1010  G4double partMom, kinE, m1 = fParticle->GetPDGMass();
1011  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1012 
1013  // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
1014 
1016 
1018 
1019  for( i = 0; i < fEnergyBin; i++)
1020  {
1021  kinE = fEnergyVector->GetLowEdgeEnergy(i);
1022 
1023  // G4cout<<G4endl;
1024  // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
1025 
1026  partMom = std::sqrt( kinE*(kinE + 2*m1) );
1027 
1028  InitDynParameters(fParticle, partMom);
1029 
1030  alphaMax = fRutherfordTheta*fCofAlphaMax;
1031 
1032  if(alphaMax > pi) alphaMax = pi;
1033 
1034  // VI: Coverity complain
1035  //alphaMax = pi2;
1036 
1037  alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
1038 
1039  // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1040 
1041  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1042 
1043  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1044 
1045  G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1046 
1047  sum = 0.;
1048 
1049  // fAddCoulomb = false;
1050  fAddCoulomb = true;
1051 
1052  // for(j = 1; j < fAngleBin; j++)
1053  for(j = fAngleBin-1; j >= 1; j--)
1054  {
1055  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1056  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1057 
1058  // alpha1 = alphaMax*(j-1)/fAngleBin;
1059  // alpha2 = alphaMax*( j )/fAngleBin;
1060 
1061  alpha1 = alphaCoulomb + delth*(j-1);
1062  // if(alpha1 < kRlim2) alpha1 = kRlim2;
1063  alpha2 = alpha1 + delth;
1064 
1065  delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelIntegrandXsc, alpha1, alpha2);
1066  // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1067 
1068  sum += delta;
1069 
1070  angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1071  // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1072  }
1073  fAngleTable->insertAt(i,angleVector);
1074 
1075  // delete[] angleVector;
1076  // delete[] angleBins;
1077  }
1078  return;
1079 }
1080 
1082 //
1083 //
1084 
1085 G4double
1087 {
1088  G4double x1, x2, y1, y2, randAngle;
1089 
1090  if( iAngle == 0 )
1091  {
1092  randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1093  // iAngle++;
1094  }
1095  else
1096  {
1097  if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1098  {
1099  iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1100  }
1101  y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1102  y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1103 
1104  x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1105  x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1106 
1107  if ( x1 == x2 ) randAngle = x2;
1108  else
1109  {
1110  if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1111  else
1112  {
1113  randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1114  }
1115  }
1116  }
1117  return randAngle;
1118 }
1119 
1120 
1121 
1123 //
1124 // Return scattering angle sampled in lab system (target at rest)
1125 
1126 
1127 
1128 G4double
1130  G4double tmass, G4double A)
1131 {
1132  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1133  G4double m1 = theParticle->GetPDGMass();
1134  G4double plab = aParticle->GetTotalMomentum();
1135  G4LorentzVector lv1 = aParticle->Get4Momentum();
1136  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1137  lv += lv1;
1138 
1139  G4ThreeVector bst = lv.boostVector();
1140  lv1.boost(-bst);
1141 
1142  G4ThreeVector p1 = lv1.vect();
1143  G4double ptot = p1.mag();
1144  G4double tmax = 4.0*ptot*ptot;
1145  G4double t = 0.0;
1146 
1147 
1148  //
1149  // Sample t
1150  //
1151 
1152  t = SampleT( theParticle, ptot, A);
1153 
1154  // NaN finder
1155  if(!(t < 0.0 || t >= 0.0))
1156  {
1157  if (verboseLevel > 0)
1158  {
1159  G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1160  << " mom(GeV)= " << plab/GeV
1161  << " S-wave will be sampled"
1162  << G4endl;
1163  }
1164  t = G4UniformRand()*tmax;
1165  }
1166  if(verboseLevel>1)
1167  {
1168  G4cout <<" t= " << t << " tmax= " << tmax
1169  << " ptot= " << ptot << G4endl;
1170  }
1171  // Sampling of angles in CM system
1172 
1174  G4double cost = 1. - 2.0*t/tmax;
1175  G4double sint;
1176 
1177  if( cost >= 1.0 )
1178  {
1179  cost = 1.0;
1180  sint = 0.0;
1181  }
1182  else if( cost <= -1.0)
1183  {
1184  cost = -1.0;
1185  sint = 0.0;
1186  }
1187  else
1188  {
1189  sint = std::sqrt((1.0-cost)*(1.0+cost));
1190  }
1191  if (verboseLevel>1)
1192  {
1193  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1194  }
1195  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1196  v1 *= ptot;
1197  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1198 
1199  nlv1.boost(bst);
1200 
1201  G4ThreeVector np1 = nlv1.vect();
1202 
1203  // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1204 
1205  G4double theta = np1.theta();
1206 
1207  return theta;
1208 }
1209 
1211 //
1212 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1213 
1214 
1215 
1216 G4double
1218  G4double tmass, G4double thetaCMS)
1219 {
1220  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1221  G4double m1 = theParticle->GetPDGMass();
1222  // G4double plab = aParticle->GetTotalMomentum();
1223  G4LorentzVector lv1 = aParticle->Get4Momentum();
1224  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1225 
1226  lv += lv1;
1227 
1228  G4ThreeVector bst = lv.boostVector();
1229 
1230  lv1.boost(-bst);
1231 
1232  G4ThreeVector p1 = lv1.vect();
1233  G4double ptot = p1.mag();
1234 
1236  G4double cost = std::cos(thetaCMS);
1237  G4double sint;
1238 
1239  if( cost >= 1.0 )
1240  {
1241  cost = 1.0;
1242  sint = 0.0;
1243  }
1244  else if( cost <= -1.0)
1245  {
1246  cost = -1.0;
1247  sint = 0.0;
1248  }
1249  else
1250  {
1251  sint = std::sqrt((1.0-cost)*(1.0+cost));
1252  }
1253  if (verboseLevel>1)
1254  {
1255  G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1256  }
1257  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1258  v1 *= ptot;
1259  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1260 
1261  nlv1.boost(bst);
1262 
1263  G4ThreeVector np1 = nlv1.vect();
1264 
1265 
1266  G4double thetaLab = np1.theta();
1267 
1268  return thetaLab;
1269 }
1270 
1272 //
1273 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1274 
1275 
1276 
1277 G4double
1279  G4double tmass, G4double thetaLab)
1280 {
1281  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1282  G4double m1 = theParticle->GetPDGMass();
1283  G4double plab = aParticle->GetTotalMomentum();
1284  G4LorentzVector lv1 = aParticle->Get4Momentum();
1285  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1286 
1287  lv += lv1;
1288 
1289  G4ThreeVector bst = lv.boostVector();
1290 
1291  // lv1.boost(-bst);
1292 
1293  // G4ThreeVector p1 = lv1.vect();
1294  // G4double ptot = p1.mag();
1295 
1297  G4double cost = std::cos(thetaLab);
1298  G4double sint;
1299 
1300  if( cost >= 1.0 )
1301  {
1302  cost = 1.0;
1303  sint = 0.0;
1304  }
1305  else if( cost <= -1.0)
1306  {
1307  cost = -1.0;
1308  sint = 0.0;
1309  }
1310  else
1311  {
1312  sint = std::sqrt((1.0-cost)*(1.0+cost));
1313  }
1314  if (verboseLevel>1)
1315  {
1316  G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1317  }
1318  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1319  v1 *= plab;
1320  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1321 
1322  nlv1.boost(-bst);
1323 
1324  G4ThreeVector np1 = nlv1.vect();
1325 
1326 
1327  G4double thetaCMS = np1.theta();
1328 
1329  return thetaCMS;
1330 }
1331 
1333 //
1334 // Test for given particle and element table of momentum, angle probability.
1335 // For the moment in lab system.
1336 
1338  G4double Z, G4double A)
1339 {
1340  fAtomicNumber = Z; // atomic number
1341  fAtomicWeight = A; // number of nucleons
1343 
1344 
1345 
1346  G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1347  <<Z<<"; and A = "<<A<<G4endl;
1348 
1350 
1351 
1352 
1353 
1354  G4int i=0, j;
1355  G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1356  G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1357  G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1358  G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1359  G4double epsilon = 0.001;
1360 
1362 
1364 
1365  fWaveVector = partMom/hbarc;
1366 
1368  G4double kR2 = kR*kR;
1369  G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1370  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1371 
1372  alphaMax = kRmax*kRmax/kR2;
1373 
1374  if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1375 
1376  alphaCoulomb = kRcoul*kRcoul/kR2;
1377 
1378  if( z )
1379  {
1380  a = partMom/m1; // beta*gamma for m1
1381  fBeta = a/std::sqrt(1+a*a);
1383  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1384  }
1385  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1386 
1387  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1388 
1389 
1390  fAddCoulomb = false;
1391 
1392  for(j = 1; j < fAngleBin; j++)
1393  {
1394  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1395  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1396 
1397  alpha1 = alphaMax*(j-1)/fAngleBin;
1398  alpha2 = alphaMax*( j )/fAngleBin;
1399 
1400  if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1401 
1402  deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1403  deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1405  alpha1, alpha2,epsilon);
1406 
1407  // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1408  // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1409 
1410  sumL10 += deltaL10;
1411  sumL96 += deltaL96;
1412  sumAG += deltaAG;
1413 
1414  G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1415  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1416 
1417  angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1418  }
1419  fAngleTable->insertAt(i,angleVector);
1420  fAngleBank.push_back(fAngleTable);
1421 
1422  /*
1423  // Integral over all angle range - Bad accuracy !!!
1424 
1425  sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1426  sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1427  sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1428  0., alpha2,epsilon);
1429  G4cout<<G4endl;
1430  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1431  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1432  */
1433  return;
1434 }
1435 
1437 //
1438 //
1439 
1441 {
1442  G4double legPol, epsilon = 1.e-6;
1443  G4double x = std::cos(theta);
1444 
1445  if ( n < 0 ) legPol = 0.;
1446  else if( n == 0 ) legPol = 1.;
1447  else if( n == 1 ) legPol = x;
1448  else if( n == 2 ) legPol = (3.*x*x-1.)/2.;
1449  else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/2.;
1450  else if( n == 4 ) legPol = (35.*x*x*x*x-30.*x*x+3.)/8.;
1451  else if( n == 5 ) legPol = (63.*x*x*x*x*x-70.*x*x*x+15.*x)/8.;
1452  else if( n == 6 ) legPol = (231.*x*x*x*x*x*x-315.*x*x*x*x+105.*x*x-5.)/16.;
1453  else
1454  {
1455  // legPol = ( (2*n-1)*x*GetLegendrePol(n-1,x) - (n-1)*GetLegendrePol(n-2,x) )/n;
1456 
1457  legPol = std::sqrt( 2./(n*CLHEP::pi*std::sin(theta+epsilon)) )*std::sin( (n+0.5)*theta+0.25*CLHEP::pi );
1458  }
1459  return legPol;
1460 }
1461 
1463 //
1464 //
1465 
1467 {
1468  G4int n;
1469  G4double n2, cofn, shny, chny, fn, gn;
1470 
1471  G4double x = z.real();
1472  G4double y = z.imag();
1473 
1474  G4double outRe = 0., outIm = 0.;
1475 
1476  G4double twox = 2.*x;
1477  G4double twoxy = twox*y;
1478  G4double twox2 = twox*twox;
1479 
1480  G4double cof1 = G4Exp(-x*x)/CLHEP::pi;
1481 
1482  G4double cos2xy = std::cos(twoxy);
1483  G4double sin2xy = std::sin(twoxy);
1484 
1485  G4double twoxcos2xy = twox*cos2xy;
1486  G4double twoxsin2xy = twox*sin2xy;
1487 
1488  for( n = 1; n <= nMax; n++)
1489  {
1490  n2 = n*n;
1491 
1492  cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n2+0.5*twox2);
1493 
1494  chny = std::cosh(n*y);
1495  shny = std::sinh(n*y);
1496 
1497  fn = twox - twoxcos2xy*chny + n*sin2xy*shny;
1498  gn = twoxsin2xy*chny + n*cos2xy*shny;
1499 
1500  fn *= cofn;
1501  gn *= cofn;
1502 
1503  outRe += fn;
1504  outIm += gn;
1505  }
1506  outRe *= 2*cof1;
1507  outIm *= 2*cof1;
1508 
1509  if(std::abs(x) < 0.0001)
1510  {
1511  outRe += GetErf(x);
1512  outIm += cof1*y;
1513  }
1514  else
1515  {
1516  outRe += GetErf(x) + cof1*(1-cos2xy)/twox;
1517  outIm += cof1*sin2xy/twox;
1518  }
1519  return G4complex(outRe, outIm);
1520 }
1521 
1522 
1524 //
1525 //
1526 
1528 {
1529  G4double outRe, outIm;
1530 
1531  G4double x = z.real();
1532  G4double y = z.imag();
1533  fReZ = x;
1534 
1536 
1537  outRe = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpSin, 0., y );
1538  outIm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetExpCos, 0., y );
1539 
1540  outRe *= 2./std::sqrt(CLHEP::pi);
1541  outIm *= 2./std::sqrt(CLHEP::pi);
1542 
1543  outRe += GetErf(x);
1544 
1545  return G4complex(outRe, outIm);
1546 }
1547 
1549 //
1550 //
1551 
1552 
1554 {
1555  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1556  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1557 
1558  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1559  G4double kappa = u/std::sqrt(CLHEP::pi);
1560  G4double dTheta = theta - fRutherfordTheta;
1561  u *= dTheta;
1562  G4double u2 = u*u;
1563  G4double u2m2p3 = u2*2./3.;
1564 
1565  G4complex im = G4complex(0.,1.);
1566  G4complex order = G4complex(u,u);
1567  order /= std::sqrt(2.);
1568 
1569  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(-order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1570  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1571  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1572  G4complex out = gamma*(1. - a1*dTheta) - a0;
1573 
1574  return out;
1575 }
1576 
1578 //
1579 //
1580 
1582 {
1583  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1584  G4double cosHalfThetaR2 = 1./(1. + fHalfRutThetaTg2);
1585 
1586  G4double u = std::sqrt(0.5*fProfileLambda/sinThetaR);
1587  G4double kappa = u/std::sqrt(CLHEP::pi);
1588  G4double dTheta = theta - fRutherfordTheta;
1589  u *= dTheta;
1590  G4double u2 = u*u;
1591  G4double u2m2p3 = u2*2./3.;
1592 
1593  G4complex im = G4complex(0.,1.);
1594  G4complex order = G4complex(u,u);
1595  order /= std::sqrt(2.);
1596  G4complex gamma = CLHEP::pi*kappa*GetErfcInt(order)*std::exp(im*(u*u+0.25*CLHEP::pi));
1597  G4complex a0 = 0.5*(1. + 4.*(1.+im*u2)*cosHalfThetaR2/3.)/sinThetaR;
1598  G4complex a1 = 0.5*(1. + 2.*(1.+im*u2m2p3)*cosHalfThetaR2)/sinThetaR;
1599  G4complex out = -gamma*(1. - a1*dTheta) - a0;
1600 
1601  return out;
1602 }
1603 
1605 //
1606 //
1607 
1609 {
1610  G4double kappa = std::sqrt(0.5*fProfileLambda/std::sin(theta)/CLHEP::pi);
1611  G4complex out = G4complex(kappa/fWaveVector,0.);
1612 
1613  out *= PhaseNear(theta);
1614 
1615  if( theta <= fRutherfordTheta )
1616  {
1617  out *= GammaLess(theta) + ProfileNear(theta);
1618  // out *= GammaMore(theta) + ProfileNear(theta);
1619  out += CoulombAmplitude(theta);
1620  }
1621  else
1622  {
1623  out *= GammaMore(theta) + ProfileNear(theta);
1624  // out *= GammaLess(theta) + ProfileNear(theta);
1625  }
1626  return out;
1627 }
1628 
1630 //
1631 //
1632 
1634 {
1635  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1636  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1637  G4double sindTheta = std::sin(dTheta);
1638  G4double persqrt2 = std::sqrt(0.5);
1639 
1640  G4complex order = G4complex(persqrt2,persqrt2);
1641  order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*sindTheta;
1642  // order *= std::sqrt(0.5*fProfileLambda/sinThetaR)*2.*dTheta;
1643 
1644  G4complex out;
1645 
1646  if ( theta <= fRutherfordTheta )
1647  {
1648  out = 1. - 0.5*GetErfcInt(-order)*ProfileNear(theta);
1649  }
1650  else
1651  {
1652  out = 0.5*GetErfcInt(order)*ProfileNear(theta);
1653  }
1654 
1655  out *= CoulombAmplitude(theta);
1656 
1657  return out;
1658 }
1659 
1661 //
1662 //
1663 
1665 {
1666  G4int n;
1667  G4double T12b, b, b2; // cosTheta = std::cos(theta);
1668  G4complex out = G4complex(0.,0.), shiftC, shiftN;
1669  G4complex im = G4complex(0.,1.);
1670 
1671  for( n = 0; n < fMaxL; n++)
1672  {
1673  shiftC = std::exp( im*2.*CalculateCoulombPhase(n) );
1674  // b = ( fZommerfeld + std::sqrt( fZommerfeld*fZommerfeld + n*(n+1) ) )/fWaveVector;
1675  b = ( std::sqrt( G4double(n*(n+1)) ) )/fWaveVector;
1676  b2 = b*b;
1678  shiftN = std::exp( -0.5*(1.-im*fEtaRatio)*T12b ) - 1.;
1679  out += (2.*n+1.)*shiftC*shiftN*GetLegendrePol(n, theta);
1680  }
1681  out /= 2.*im*fWaveVector;
1682  out += CoulombAmplitude(theta);
1683  return out;
1684 }
1685 
1686 
1688 //
1689 //
1690 
1692 {
1693  G4int n;
1694  G4double T12b, a, aTemp, b2, sinThetaH = std::sin(0.5*theta);
1695  G4double sinThetaH2 = sinThetaH*sinThetaH;
1696  G4complex out = G4complex(0.,0.);
1697  G4complex im = G4complex(0.,1.);
1698 
1701 
1702  aTemp = a;
1703 
1704  for( n = 1; n < fMaxL; n++)
1705  {
1706  T12b = aTemp*G4Exp(-b2/n)/n;
1707  aTemp *= a;
1708  out += T12b;
1709  G4cout<<"out = "<<out<<G4endl;
1710  }
1711  out *= -4.*im*fWaveVector/CLHEP::pi;
1712  out += CoulombAmplitude(theta);
1713  return out;
1714 }
1715 
1716 
1718 //
1719 // Test for given particle and element table of momentum, angle probability.
1720 // For the partMom in CMS.
1721 
1723  G4double partMom, G4double Z, G4double A)
1724 {
1725  fAtomicNumber = Z; // atomic number
1726  fAtomicWeight = A; // number of nucleons
1727 
1729  G4double A1 = G4double( theParticle->GetBaryonNumber() );
1731  // fNuclearRadius = std::sqrt(fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2);
1732  fNuclearRadius = fNuclearRadius1 + fNuclearRadius2;
1733 
1734  G4double a = 0.;
1735  G4double z = theParticle->GetPDGCharge();
1736  G4double m1 = theParticle->GetPDGMass();
1737 
1738  fWaveVector = partMom/CLHEP::hbarc;
1739 
1741  G4cout<<"kR = "<<lambda<<G4endl;
1742 
1743  if( z )
1744  {
1745  a = partMom/m1; // beta*gamma for m1
1746  fBeta = a/std::sqrt(1+a*a);
1749  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1750  }
1751  G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
1752  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1753  G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
1756 
1759 
1760  return;
1761 }
1763 //
1764 // Test for given particle and element table of momentum, angle probability.
1765 // For the partMom in CMS.
1766 
1768  G4double partMom)
1769 {
1770  G4double a = 0.;
1771  G4double z = theParticle->GetPDGCharge();
1772  G4double m1 = theParticle->GetPDGMass();
1773 
1774  fWaveVector = partMom/CLHEP::hbarc;
1775 
1777 
1778  if( z )
1779  {
1780  a = partMom/m1; // beta*gamma for m1
1781  fBeta = a/std::sqrt(1+a*a);
1784  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1785  }
1786  fProfileLambda = lambda; // *std::sqrt(1.-2*fZommerfeld/lambda);
1789 
1792 
1793  return;
1794 }
1795 
1796 
1798 //
1799 // Test for given particle and element table of momentum, angle probability.
1800 // For the partMom in CMS.
1801 
1803  G4double partMom, G4double Z, G4double A)
1804 {
1805  fAtomicNumber = Z; // target atomic number
1806  fAtomicWeight = A; // target number of nucleons
1807 
1808  fNuclearRadius2 = CalculateNuclearRad(fAtomicWeight); // target nucleus radius
1809  G4double A1 = G4double( aParticle->GetDefinition()->GetBaryonNumber() );
1810  fNuclearRadius1 = CalculateNuclearRad(A1); // projectile nucleus radius
1811  fNuclearRadiusSquare = fNuclearRadius1*fNuclearRadius1+fNuclearRadius2*fNuclearRadius2;
1812 
1813 
1814  G4double a = 0., kR12;
1815  G4double z = aParticle->GetDefinition()->GetPDGCharge();
1816  G4double m1 = aParticle->GetDefinition()->GetPDGMass();
1817 
1818  fWaveVector = partMom/CLHEP::hbarc;
1819 
1820  G4double pN = A1 - z;
1821  if( pN < 0. ) pN = 0.;
1822 
1823  G4double tN = A - Z;
1824  if( tN < 0. ) tN = 0.;
1825 
1826  G4double pTkin = aParticle->GetKineticEnergy();
1827  pTkin /= A1;
1828 
1829 
1830  fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
1831  (z*tN+pN*Z)*GetHadronNucleonXscNS(theProton, pTkin, theNeutron);
1832 
1833  G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::millibarn<<" mb"<<G4endl;
1835  kR12 = fWaveVector*std::sqrt(fNuclearRadiusSquare);
1836  G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl;
1837  fMaxL = (G4int(kR12)+1)*4;
1838  G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl;
1839 
1840  if( z )
1841  {
1842  a = partMom/m1; // beta*gamma for m1
1843  fBeta = a/std::sqrt(1+a*a);
1845  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1846  }
1847 
1849 
1850 
1851  return;
1852 }
1853 
1854 
1856 //
1857 // Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
1858 // data from mainly http://wwwppds.ihep.su:8001/c5-6A.html database
1859 // projectile nucleon is pParticle with pTkin shooting target nucleon tParticle
1860 
1861 G4double
1863  G4double pTkin,
1864  G4ParticleDefinition* tParticle)
1865 {
1866  G4double xsection(0), /*Delta,*/ A0, B0;
1867  G4double hpXsc(0);
1868  G4double hnXsc(0);
1869 
1870 
1871  G4double targ_mass = tParticle->GetPDGMass();
1872  G4double proj_mass = pParticle->GetPDGMass();
1873 
1874  G4double proj_energy = proj_mass + pTkin;
1875  G4double proj_momentum = std::sqrt(pTkin*(pTkin+2*proj_mass));
1876 
1877  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
1878 
1879  sMand /= CLHEP::GeV*CLHEP::GeV; // in GeV for parametrisation
1880  proj_momentum /= CLHEP::GeV;
1881  proj_energy /= CLHEP::GeV;
1882  proj_mass /= CLHEP::GeV;
1883  G4double logS = G4Log(sMand);
1884 
1885  // General PDG fit constants
1886 
1887 
1888  // fEtaRatio=Re[f(0)]/Im[f(0)]
1889 
1890  if( proj_momentum >= 1.2 )
1891  {
1892  fEtaRatio = 0.13*(logS - 5.8579332)*G4Pow::GetInstance()->powA(sMand,-0.18);
1893  }
1894  else if( proj_momentum >= 0.6 )
1895  {
1896  fEtaRatio = -75.5*(G4Pow::GetInstance()->powA(proj_momentum,0.25)-0.95)/
1897  (G4Pow::GetInstance()->powA(3*proj_momentum,2.2)+1);
1898  }
1899  else
1900  {
1901  fEtaRatio = 15.5*proj_momentum/(27*proj_momentum*proj_momentum*proj_momentum+2);
1902  }
1903  G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl;
1904 
1905  // xsc
1906 
1907  if( proj_momentum >= 10. ) // high energy: pp = nn = np
1908  // if( proj_momentum >= 2.)
1909  {
1910  //Delta = 1.;
1911 
1912  //if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
1913 
1914  //AR-12Aug2016 if( proj_momentum >= 10.)
1915  {
1916  B0 = 7.5;
1917  A0 = 100. - B0*G4Log(3.0e7);
1918 
1919  xsection = A0 + B0*G4Log(proj_energy) - 11
1920  + 103*G4Pow::GetInstance()->powA(2*0.93827*proj_energy + proj_mass*proj_mass+
1921  0.93827*0.93827,-0.165); // mb
1922  }
1923  }
1924  else // low energy pp = nn != np
1925  {
1926  if(pParticle == tParticle) // pp or nn // nn to be pp
1927  {
1928  if( proj_momentum < 0.73 )
1929  {
1930  hnXsc = 23 + 50*( G4Pow::GetInstance()->powA( G4Log(0.73/proj_momentum), 3.5 ) );
1931  }
1932  else if( proj_momentum < 1.05 )
1933  {
1934  hnXsc = 23 + 40*(G4Log(proj_momentum/0.73))*
1935  (G4Log(proj_momentum/0.73));
1936  }
1937  else // if( proj_momentum < 10. )
1938  {
1939  hnXsc = 39.0 +
1940  75*(proj_momentum - 1.2)/(G4Pow::GetInstance()->powA(proj_momentum,3.0) + 0.15);
1941  }
1942  xsection = hnXsc;
1943  }
1944  else // pn to be np
1945  {
1946  if( proj_momentum < 0.8 )
1947  {
1948  hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/1.3),4.0);
1949  }
1950  else if( proj_momentum < 1.4 )
1951  {
1952  hpXsc = 33+30*G4Pow::GetInstance()->powA(G4Log(proj_momentum/0.95),2.0);
1953  }
1954  else // if( proj_momentum < 10. )
1955  {
1956  hpXsc = 33.3+
1957  20.8*(G4Pow::GetInstance()->powA(proj_momentum,2.0)-1.35)/
1958  (G4Pow::GetInstance()->powA(proj_momentum,2.50)+0.95);
1959  }
1960  xsection = hpXsc;
1961  }
1962  }
1963  xsection *= CLHEP::millibarn; // parametrised in mb
1964  G4cout<<"xsection = "<<xsection/CLHEP::millibarn<<" mb"<<G4endl;
1965  return xsection;
1966 }
1967 
1969 //
1970 // The ratio el/ruth for Fresnel smooth nucleus profile
1971 
1973 {
1974  G4double sinThetaR = 2.*fHalfRutThetaTg/(1. + fHalfRutThetaTg2);
1975  G4double dTheta = 0.5*(theta - fRutherfordTheta);
1976  G4double sindTheta = std::sin(dTheta);
1977 
1978  G4double prof = Profile(theta);
1979  G4double prof2 = prof*prof;
1980  // G4double profmod = std::abs(prof);
1981  G4double order = std::sqrt(fProfileLambda/sinThetaR/CLHEP::pi)*2.*sindTheta;
1982 
1983  order = std::abs(order); // since sin changes sign!
1984  // G4cout<<"order = "<<order<<G4endl;
1985 
1986  G4double cosFresnel = GetCint(order);
1987  G4double sinFresnel = GetSint(order);
1988 
1989  G4double out;
1990 
1991  if ( theta <= fRutherfordTheta )
1992  {
1993  out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1994  out += ( cosFresnel + sinFresnel - 1. )*prof;
1995  }
1996  else
1997  {
1998  out = 0.5*( (0.5-cosFresnel)*(0.5-cosFresnel)+(0.5-sinFresnel)*(0.5-sinFresnel) )*prof2;
1999  }
2000 
2001  return out;
2002 }
2003 
2005 //
2006 // For the calculation of arg Gamma(z) one needs complex extension
2007 // of ln(Gamma(z))
2008 
2010 {
2011  const G4double cof[6] = { 76.18009172947146, -86.50532032941677,
2012  24.01409824083091, -1.231739572450155,
2013  0.1208650973866179e-2, -0.5395239384953e-5 } ;
2014  G4int j;
2015  G4complex z = zz - 1.0;
2016  G4complex tmp = z + 5.5;
2017  tmp -= (z + 0.5) * std::log(tmp);
2018  G4complex ser = G4complex(1.000000000190015,0.);
2019 
2020  for ( j = 0; j <= 5; j++ )
2021  {
2022  z += 1.0;
2023  ser += cof[j]/z;
2024  }
2025  return -tmp + std::log(2.5066282746310005*ser);
2026 }
2027 
2029 //
2030 // Bessel J0 function based on rational approximation from
2031 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2032 
2034 {
2035  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2036 
2037  modvalue = std::fabs(value);
2038 
2039  if ( value < 8.0 && value > -8.0 )
2040  {
2041  value2 = value*value;
2042 
2043  fact1 = 57568490574.0 + value2*(-13362590354.0
2044  + value2*( 651619640.7
2045  + value2*(-11214424.18
2046  + value2*( 77392.33017
2047  + value2*(-184.9052456 ) ) ) ) );
2048 
2049  fact2 = 57568490411.0 + value2*( 1029532985.0
2050  + value2*( 9494680.718
2051  + value2*(59272.64853
2052  + value2*(267.8532712
2053  + value2*1.0 ) ) ) );
2054 
2055  bessel = fact1/fact2;
2056  }
2057  else
2058  {
2059  arg = 8.0/modvalue;
2060 
2061  value2 = arg*arg;
2062 
2063  shift = modvalue-0.785398164;
2064 
2065  fact1 = 1.0 + value2*(-0.1098628627e-2
2066  + value2*(0.2734510407e-4
2067  + value2*(-0.2073370639e-5
2068  + value2*0.2093887211e-6 ) ) );
2069 
2070  fact2 = -0.1562499995e-1 + value2*(0.1430488765e-3
2071  + value2*(-0.6911147651e-5
2072  + value2*(0.7621095161e-6
2073  - value2*0.934945152e-7 ) ) );
2074 
2075  bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
2076  }
2077  return bessel;
2078 }
2079 
2081 //
2082 // Bessel J1 function based on rational approximation from
2083 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141
2084 
2086 {
2087  G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
2088 
2089  modvalue = std::fabs(value);
2090 
2091  if ( modvalue < 8.0 )
2092  {
2093  value2 = value*value;
2094 
2095  fact1 = value*(72362614232.0 + value2*(-7895059235.0
2096  + value2*( 242396853.1
2097  + value2*(-2972611.439
2098  + value2*( 15704.48260
2099  + value2*(-30.16036606 ) ) ) ) ) );
2100 
2101  fact2 = 144725228442.0 + value2*(2300535178.0
2102  + value2*(18583304.74
2103  + value2*(99447.43394
2104  + value2*(376.9991397
2105  + value2*1.0 ) ) ) );
2106  bessel = fact1/fact2;
2107  }
2108  else
2109  {
2110  arg = 8.0/modvalue;
2111 
2112  value2 = arg*arg;
2113 
2114  shift = modvalue - 2.356194491;
2115 
2116  fact1 = 1.0 + value2*( 0.183105e-2
2117  + value2*(-0.3516396496e-4
2118  + value2*(0.2457520174e-5
2119  + value2*(-0.240337019e-6 ) ) ) );
2120 
2121  fact2 = 0.04687499995 + value2*(-0.2002690873e-3
2122  + value2*( 0.8449199096e-5
2123  + value2*(-0.88228987e-6
2124  + value2*0.105787412e-6 ) ) );
2125 
2126  bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
2127 
2128  if (value < 0.0) bessel = -bessel;
2129  }
2130  return bessel;
2131 }
2132 
2133 //
2134 //