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