Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4NuclNuclDiffuseElastic.cc 106722 2017-10-20 09:48:19Z gcosmo $
27 //
28 //
29 // Physics model class G4NuclNuclDiffuseElastic
30 //
31 //
32 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
33 //
34 // 24-May-07 V. Grichine
35 //
36 
38 #include "G4ParticleTable.hh"
39 #include "G4ParticleDefinition.hh"
40 #include "G4IonTable.hh"
41 #include "G4NucleiProperties.hh"
42 
43 #include "Randomize.hh"
44 #include "G4Integrator.hh"
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 
49 #include "G4Proton.hh"
50 #include "G4Neutron.hh"
51 #include "G4Deuteron.hh"
52 #include "G4Alpha.hh"
53 #include "G4PionPlus.hh"
54 #include "G4PionMinus.hh"
55 
56 #include "G4Element.hh"
57 #include "G4ElementTable.hh"
58 #include "G4NistManager.hh"
59 #include "G4PhysicsTable.hh"
60 #include "G4PhysicsLogVector.hh"
61 #include "G4PhysicsFreeVector.hh"
62 
64 //
65 // Test Constructor. Just to check xsc
66 
67 
69  : G4HadronElastic("NNDiffuseElastic"), fParticle(0)
70 {
71  SetMinEnergy( 50*MeV );
72  SetMaxEnergy( 1.*TeV );
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 = 200;
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,
192  G4double momentum,
193  G4double A )
194 {
195  fParticle = particle;
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,
263  G4double momentum,
264  G4double A, G4double Z )
265 {
266  fParticle = particle;
267  fWaveVector = momentum/hbarc;
268  fAtomicWeight = A;
269  fAtomicNumber = Z;
271  fAddCoulomb = false;
272 
273  G4double z = particle->GetPDGCharge();
274 
275  G4double kRt = fWaveVector*fNuclearRadius*theta;
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 {
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,
686  G4double momentum,
687  G4double A )
688 {
690  fParticle = particle;
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
724  G4double momentum, G4double A)
725 {
726  G4int i, iMax = 100;
727  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
728 
729  fParticle = particle;
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
856  G4double momentum, G4double Z, G4double A)
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 
1173  G4double phi = G4UniformRand()*twopi;
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 
1235  G4double phi = G4UniformRand()*twopi;
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 
1296  G4double phi = G4UniformRand()*twopi;
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 //
Float_t x
Definition: compare.C:6
const G4ParticleDefinition * thePionPlus
G4double CalculateNuclearRad(G4double A)
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
static G4He3 * He3()
Definition: G4He3.cc:94
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4complex GetErfcInt(G4complex z)
G4double GetDiffElasticProb(G4double theta)
G4complex AmplitudeGG(G4double theta)
void SetMinEnergy(G4double anEnergy)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void SetMaxEnergy(const G4double anEnergy)
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
Definition: G4SIunits.hh:214
G4double GetLowEdgeEnergy(size_t binNumber) const
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
G4double GetDiffElasticSumProbA(G4double alpha)
Float_t y1[n_points_granero]
Definition: compare.C:5
G4complex GetErfComp(G4complex z, G4int nMax)
G4complex GammaMore(G4double theta)
void InitialiseOnFly(G4double Z, G4double A)
const G4ParticleDefinition * fParticle
static constexpr double keV
Definition: G4SIunits.hh:216
Float_t x1[n_points_granero]
Definition: compare.C:5
G4complex CoulombAmplitude(G4double theta)
#define G4endl
Definition: G4ios.hh:61
Float_t y
Definition: compare.C:6
const char * p
Definition: xmltok.h:285
static constexpr double hbarc
Double_t z
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4double GetDiffElasticSumProb(G4double theta)
G4IonTable * GetIonTable() const
G4double GetPDGCharge() const
G4double GetTotalMomentum() const
double z() const
Double_t zz
void InitDynParameters(const G4ParticleDefinition *theParticle, G4double partMom)
Float_t tmp
G4complex GammaLogarithm(G4complex xx)
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
static G4Proton * Proton()
Definition: G4Proton.cc:93
void insertAt(size_t, G4PhysicsVector *)
G4ParticleDefinition * theDeuteron
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
std::vector< G4PhysicsTable * > fAngleBank
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
G4double GetAtomicMassAmu(const G4String &symb) const
G4double GetPDGMass() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double theta() const
std::vector< G4String > fElementNameVector
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4complex GammaLess(G4double theta)
Float_t y2[n_points_geant4]
Definition: compare.C:26
G4double GetRatioGen(G4double theta)
Float_t Z
G4ParticleDefinition * theNeutron
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
static constexpr double TeV
Definition: G4SIunits.hh:218
G4double SampleCoulombMuCMS(const G4ParticleDefinition *aParticle, G4double p)
G4double GetHadronNucleonXscNS(G4ParticleDefinition *pParticle, G4double pTkin, G4ParticleDefinition *tParticle)
std::complex< G4double > G4complex
Definition: G4Types.hh:81
double G4double
Definition: G4Types.hh:76
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
static constexpr double fermi
Definition: G4SIunits.hh:103
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
G4ParticleDefinition * GetDefinition() const
#define position
Definition: xmlparse.cc:622
const XML_Char int const XML_Char * value
Definition: expat.h:331
void InitParametersGla(const G4DynamicParticle *aParticle, G4double partMom, G4double Z, G4double A)
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
G4double GetFresnelIntegrandXsc(G4double alpha)
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
std::vector< G4Element * > G4ElementTable
static const G4double alpha
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
G4ParticleDefinition * theProton
G4LorentzVector Get4Momentum() const
double A(double temperature)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double CalcMandelstamS(const G4double mp, const G4double mt, const G4double Plab)
ThreeVector shoot(const G4int Ap, const G4int Af)
const G4LorentzVector & Get4Momentum() const
G4double G4ParticleHPJENDLHEData::G4double result
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
double epsilon(double density, double temperature)
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4complex AmplitudeNear(G4double theta)
G4complex PhaseNear(G4double theta)
const G4ParticleDefinition * GetDefinition() const
const G4double a0
double mag() const
int G4int
Definition: G4Types.hh:78
void InitParameters(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
Float_t norm
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double ProfileNear(G4double theta)
static constexpr double degree
Definition: G4SIunits.hh:144
G4double GetKineticEnergy() const
G4GLOB_DLL std::ostream G4cout
double x() const
Hep3Vector boostVector() const
std::vector< G4double > fElementNumberVector
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
Hep3Vector vect() const
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
Char_t n[5]
static size_t GetNumberOfElements()
Definition: G4Element.cc:405
void PutValue(size_t index, G4double energy, G4double dataValue)
G4double GetTotalMomentum() const
static constexpr double pi
Definition: G4SIunits.hh:75
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:504
double y() const
G4double Profile(G4double theta)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
Float_t x2[n_points_geant4]
Definition: compare.C:26
const G4ParticleDefinition * thePionMinus
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
static constexpr double GeV
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
static constexpr double GeV
Definition: G4SIunits.hh:217
G4double GetIntegrandFunction(G4double theta)
static constexpr double pi
Definition: SystemOfUnits.h:54
static constexpr double millibarn
Definition: SystemOfUnits.h:86
G4double GetLegendrePol(G4int n, G4double x)
HepLorentzVector & boost(double, double, double)
G4complex AmplitudeGla(G4double theta)
static G4NistManager * Instance()
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:398
G4complex AmplitudeSim(G4double theta)