Geant4_10
G4DiffuseElastic.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: G4DiffuseElastic.cc 70682 2013-06-04 07:57:01Z gcosmo $
27 //
28 //
29 // Physics model class G4DiffuseElastic
30 //
31 //
32 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
33 //
34 // 24-May-07 V. Grichine
35 //
36 
37 #include "G4DiffuseElastic.hh"
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 "G4PhysicsTable.hh"
59 #include "G4PhysicsLogVector.hh"
60 #include "G4PhysicsFreeVector.hh"
61 
63 //
64 // Test Constructor. Just to check xsc
65 
66 
68  : G4HadronElastic("DiffuseElastic"), fParticle(0)
69 {
70  SetMinEnergy( 0.01*GeV );
71  SetMaxEnergy( 1.*TeV );
72  verboseLevel = 0;
73  lowEnergyRecoilLimit = 100.*keV;
74  lowEnergyLimitQ = 0.0*GeV;
75  lowEnergyLimitHE = 0.0*GeV;
76  lowestEnergyLimit= 0.0*keV;
77  plabLowLimit = 20.0*MeV;
78 
79  theProton = G4Proton::Proton();
80  theNeutron = G4Neutron::Neutron();
81  theDeuteron = G4Deuteron::Deuteron();
82  theAlpha = G4Alpha::Alpha();
83  thePionPlus = G4PionPlus::PionPlus();
84  thePionMinus= G4PionMinus::PionMinus();
85 
86  fEnergyBin = 200;
87  fAngleBin = 200;
88 
89  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
90  fAngleTable = 0;
91 
92  fParticle = 0;
93  fWaveVector = 0.;
94  fAtomicWeight = 0.;
95  fAtomicNumber = 0.;
96  fNuclearRadius = 0.;
97  fBeta = 0.;
98  fZommerfeld = 0.;
99  fAm = 0.;
100  fAddCoulomb = false;
101 }
102 
104 //
105 // Destructor
106 
108 {
109  if(fEnergyVector) delete fEnergyVector;
110 
111  if( fAngleTable )
112  {
113  fAngleTable->clearAndDestroy();
114  delete fAngleTable ;
115  }
116 }
117 
119 //
120 // Initialisation for given particle using element table of application
121 
123 {
124 
125  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
126 
127  const G4ElementTable* theElementTable = G4Element::GetElementTable();
128 
129  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
130 
131  for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
132  {
133  fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
134  fAtomicWeight = (*theElementTable)[jEl]->GetN(); // number of nucleons
135  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
136 
137  if(verboseLevel > 0)
138  {
139  G4cout<<"G4DiffuseElastic::Initialise() the element: "
140  <<(*theElementTable)[jEl]->GetName()<<G4endl;
141  }
142  fElementNumberVector.push_back(fAtomicNumber);
143  fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
144 
145  BuildAngleTable();
146  fAngleBank.push_back(fAngleTable);
147  }
148  return;
149 }
150 
152 //
153 // return differential elastic cross section d(sigma)/d(omega)
154 
155 G4double
157  G4double theta,
158  G4double momentum,
159  G4double A )
160 {
161  fParticle = particle;
162  fWaveVector = momentum/hbarc;
163  fAtomicWeight = A;
164  fAddCoulomb = false;
165  fNuclearRadius = CalculateNuclearRad(A);
166 
167  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
168 
169  return sigma;
170 }
171 
173 //
174 // return invariant differential elastic cross section d(sigma)/d(tMand)
175 
176 G4double
178  G4double tMand,
179  G4double plab,
180  G4double A, G4double Z )
181 {
182  G4double m1 = particle->GetPDGMass();
183  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
184 
185  G4int iZ = static_cast<G4int>(Z+0.5);
186  G4int iA = static_cast<G4int>(A+0.5);
187  G4ParticleDefinition * theDef = 0;
188 
189  if (iZ == 1 && iA == 1) theDef = theProton;
190  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
191  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
192  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
193  else if (iZ == 2 && iA == 4) theDef = theAlpha;
194  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
195 
196  G4double tmass = theDef->GetPDGMass();
197 
198  G4LorentzVector lv(0.0,0.0,0.0,tmass);
199  lv += lv1;
200 
201  G4ThreeVector bst = lv.boostVector();
202  lv1.boost(-bst);
203 
204  G4ThreeVector p1 = lv1.vect();
205  G4double ptot = p1.mag();
206  G4double ptot2 = ptot*ptot;
207  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
208 
209  if( cost >= 1.0 ) cost = 1.0;
210  else if( cost <= -1.0) cost = -1.0;
211 
212  G4double thetaCMS = std::acos(cost);
213 
214  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
215 
216  sigma *= pi/ptot2;
217 
218  return sigma;
219 }
220 
222 //
223 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
224 // correction
225 
226 G4double
228  G4double theta,
229  G4double momentum,
230  G4double A, G4double Z )
231 {
232  fParticle = particle;
233  fWaveVector = momentum/hbarc;
234  fAtomicWeight = A;
235  fAtomicNumber = Z;
236  fNuclearRadius = CalculateNuclearRad(A);
237  fAddCoulomb = false;
238 
239  G4double z = particle->GetPDGCharge();
240 
241  G4double kRt = fWaveVector*fNuclearRadius*theta;
242  G4double kRtC = 1.9;
243 
244  if( z && (kRt > kRtC) )
245  {
246  fAddCoulomb = true;
247  fBeta = CalculateParticleBeta( particle, momentum);
248  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
249  fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
250  }
251  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
252 
253  return sigma;
254 }
255 
257 //
258 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
259 // correction
260 
261 G4double
263  G4double tMand,
264  G4double plab,
265  G4double A, G4double Z )
266 {
267  G4double m1 = particle->GetPDGMass();
268 
269  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
270 
271  G4int iZ = static_cast<G4int>(Z+0.5);
272  G4int iA = static_cast<G4int>(A+0.5);
273 
274  G4ParticleDefinition* theDef = 0;
275 
276  if (iZ == 1 && iA == 1) theDef = theProton;
277  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
278  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
279  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
280  else if (iZ == 2 && iA == 4) theDef = theAlpha;
281  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
282 
283  G4double tmass = theDef->GetPDGMass();
284 
285  G4LorentzVector lv(0.0,0.0,0.0,tmass);
286  lv += lv1;
287 
288  G4ThreeVector bst = lv.boostVector();
289  lv1.boost(-bst);
290 
291  G4ThreeVector p1 = lv1.vect();
292  G4double ptot = p1.mag();
293  G4double ptot2 = ptot*ptot;
294  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
295 
296  if( cost >= 1.0 ) cost = 1.0;
297  else if( cost <= -1.0) cost = -1.0;
298 
299  G4double thetaCMS = std::acos(cost);
300 
301  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
302 
303  sigma *= pi/ptot2;
304 
305  return sigma;
306 }
307 
309 //
310 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
311 // correction
312 
313 G4double
315  G4double tMand,
316  G4double plab,
317  G4double A, G4double Z )
318 {
319  G4double m1 = particle->GetPDGMass();
320  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
321 
322  G4int iZ = static_cast<G4int>(Z+0.5);
323  G4int iA = static_cast<G4int>(A+0.5);
324  G4ParticleDefinition * theDef = 0;
325 
326  if (iZ == 1 && iA == 1) theDef = theProton;
327  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
328  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
329  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
330  else if (iZ == 2 && iA == 4) theDef = theAlpha;
331  else theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(iZ,iA,0);
332 
333  G4double tmass = theDef->GetPDGMass();
334 
335  G4LorentzVector lv(0.0,0.0,0.0,tmass);
336  lv += lv1;
337 
338  G4ThreeVector bst = lv.boostVector();
339  lv1.boost(-bst);
340 
341  G4ThreeVector p1 = lv1.vect();
342  G4double ptot = p1.mag();
343  G4double ptot2 = ptot*ptot;
344  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
345 
346  if( cost >= 1.0 ) cost = 1.0;
347  else if( cost <= -1.0) cost = -1.0;
348 
349  G4double thetaCMS = std::acos(cost);
350 
351  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
352 
353  sigma *= pi/ptot2;
354 
355  return sigma;
356 }
357 
359 //
360 // return differential elastic probability d(probability)/d(omega)
361 
362 G4double
363 G4DiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
364  G4double theta
365  // G4double momentum,
366  // G4double A
367  )
368 {
369  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
370  G4double delta, diffuse, gamma;
371  G4double e1, e2, bone, bone2;
372 
373  // G4double wavek = momentum/hbarc; // wave vector
374  // G4double r0 = 1.08*fermi;
375  // G4double rad = r0*std::pow(A, 1./3.);
376 
377  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
378  G4double kr2 = kr*kr;
379  G4double krt = kr*theta;
380 
381  bzero = BesselJzero(krt);
382  bzero2 = bzero*bzero;
383  bone = BesselJone(krt);
384  bone2 = bone*bone;
385  bonebyarg = BesselOneByArg(krt);
386  bonebyarg2 = bonebyarg*bonebyarg;
387 
388  if (fParticle == theProton)
389  {
390  diffuse = 0.63*fermi;
391  gamma = 0.3*fermi;
392  delta = 0.1*fermi*fermi;
393  e1 = 0.3*fermi;
394  e2 = 0.35*fermi;
395  }
396  else // as proton, if were not defined
397  {
398  diffuse = 0.63*fermi;
399  gamma = 0.3*fermi;
400  delta = 0.1*fermi*fermi;
401  e1 = 0.3*fermi;
402  e2 = 0.35*fermi;
403  }
404  G4double lambda = 15.; // 15 ok
405 
406  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
407 
408  G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
409  G4double kgamma2 = kgamma*kgamma;
410 
411  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
412  // G4double dk2t2 = dk2t*dk2t;
413  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
414 
415  G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
416 
417  damp = DampFactor(pikdt);
418  damp2 = damp*damp;
419 
420  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
421  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
422 
423 
424  sigma = kgamma2;
425  // sigma += dk2t2;
426  sigma *= bzero2;
427  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
428  sigma += kr2*bonebyarg2;
429  sigma *= damp2; // *rad*rad;
430 
431  return sigma;
432 }
433 
435 //
436 // return differential elastic probability d(probability)/d(omega) with
437 // Coulomb correction
438 
439 G4double
440 G4DiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
441  G4double theta
442  // G4double momentum,
443  // G4double A
444  )
445 {
446  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
447  G4double delta, diffuse, gamma;
448  G4double e1, e2, bone, bone2;
449 
450  // G4double wavek = momentum/hbarc; // wave vector
451  // G4double r0 = 1.08*fermi;
452  // G4double rad = r0*std::pow(A, 1./3.);
453 
454  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
455  G4double kr2 = kr*kr;
456  G4double krt = kr*theta;
457 
458  bzero = BesselJzero(krt);
459  bzero2 = bzero*bzero;
460  bone = BesselJone(krt);
461  bone2 = bone*bone;
462  bonebyarg = BesselOneByArg(krt);
463  bonebyarg2 = bonebyarg*bonebyarg;
464 
465  if (fParticle == theProton)
466  {
467  diffuse = 0.63*fermi;
468  // diffuse = 0.6*fermi;
469  gamma = 0.3*fermi;
470  delta = 0.1*fermi*fermi;
471  e1 = 0.3*fermi;
472  e2 = 0.35*fermi;
473  }
474  else // as proton, if were not defined
475  {
476  diffuse = 0.63*fermi;
477  gamma = 0.3*fermi;
478  delta = 0.1*fermi*fermi;
479  e1 = 0.3*fermi;
480  e2 = 0.35*fermi;
481  }
482  G4double lambda = 15.; // 15 ok
483  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
484  G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
485 
486  // G4cout<<"kgamma = "<<kgamma<<G4endl;
487 
488  if(fAddCoulomb) // add Coulomb correction
489  {
490  G4double sinHalfTheta = std::sin(0.5*theta);
491  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
492 
493  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
494  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
495  }
496 
497  G4double kgamma2 = kgamma*kgamma;
498 
499  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
500  // G4cout<<"dk2t = "<<dk2t<<G4endl;
501  // G4double dk2t2 = dk2t*dk2t;
502  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
503 
504  G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
505 
506  // G4cout<<"pikdt = "<<pikdt<<G4endl;
507 
508  damp = DampFactor(pikdt);
509  damp2 = damp*damp;
510 
511  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
512  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
513 
514  sigma = kgamma2;
515  // sigma += dk2t2;
516  sigma *= bzero2;
517  sigma += mode2k2*bone2;
518  sigma += e2dk3t*bzero*bone;
519 
520  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
521  sigma += kr2*bonebyarg2; // correction at J1()/()
522 
523  sigma *= damp2; // *rad*rad;
524 
525  return sigma;
526 }
527 
528 
530 //
531 // return differential elastic probability d(probability)/d(t) with
532 // Coulomb correction
533 
534 G4double
536 {
537  G4double theta;
538 
539  theta = std::sqrt(alpha);
540 
541  // theta = std::acos( 1 - alpha/2. );
542 
543  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
544  G4double delta, diffuse, gamma;
545  G4double e1, e2, bone, bone2;
546 
547  // G4double wavek = momentum/hbarc; // wave vector
548  // G4double r0 = 1.08*fermi;
549  // G4double rad = r0*std::pow(A, 1./3.);
550 
551  G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
552  G4double kr2 = kr*kr;
553  G4double krt = kr*theta;
554 
555  bzero = BesselJzero(krt);
556  bzero2 = bzero*bzero;
557  bone = BesselJone(krt);
558  bone2 = bone*bone;
559  bonebyarg = BesselOneByArg(krt);
560  bonebyarg2 = bonebyarg*bonebyarg;
561 
562  if (fParticle == theProton)
563  {
564  diffuse = 0.63*fermi;
565  // diffuse = 0.6*fermi;
566  gamma = 0.3*fermi;
567  delta = 0.1*fermi*fermi;
568  e1 = 0.3*fermi;
569  e2 = 0.35*fermi;
570  }
571  else // as proton, if were not defined
572  {
573  diffuse = 0.63*fermi;
574  gamma = 0.3*fermi;
575  delta = 0.1*fermi*fermi;
576  e1 = 0.3*fermi;
577  e2 = 0.35*fermi;
578  }
579  G4double lambda = 15.; // 15 ok
580  // G4double kgamma = fWaveVector*gamma; // wavek*delta;
581  G4double kgamma = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
582 
583  // G4cout<<"kgamma = "<<kgamma<<G4endl;
584 
585  if(fAddCoulomb) // add Coulomb correction
586  {
587  G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
588  G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
589 
590  kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
591  // kgamma += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
592  }
593 
594  G4double kgamma2 = kgamma*kgamma;
595 
596  // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
597  // G4cout<<"dk2t = "<<dk2t<<G4endl;
598  // G4double dk2t2 = dk2t*dk2t;
599  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
600 
601  G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
602 
603  // G4cout<<"pikdt = "<<pikdt<<G4endl;
604 
605  damp = DampFactor(pikdt);
606  damp2 = damp*damp;
607 
608  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
609  G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
610 
611  sigma = kgamma2;
612  // sigma += dk2t2;
613  sigma *= bzero2;
614  sigma += mode2k2*bone2;
615  sigma += e2dk3t*bzero*bone;
616 
617  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
618  sigma += kr2*bonebyarg2; // correction at J1()/()
619 
620  sigma *= damp2; // *rad*rad;
621 
622  return sigma;
623 }
624 
625 
627 //
628 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
629 
630 G4double
632 {
634 
635  result = GetDiffElasticSumProbA(alpha);
636 
637  // result *= 2*pi*std::sin(theta);
638 
639  return result;
640 }
641 
643 //
644 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
645 
646 G4double
648  G4double theta,
649  G4double momentum,
650  G4double A )
651 {
653  fParticle = particle;
654  fWaveVector = momentum/hbarc;
655  fAtomicWeight = A;
656 
657  fNuclearRadius = CalculateNuclearRad(A);
658 
659 
661 
662  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
663  result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
664 
665  return result;
666 }
667 
669 //
670 // Return inv momentum transfer -t > 0
671 
673 {
674  G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
675  G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
676  return t;
677 }
678 
680 //
681 // Return scattering angle sampled in cms
682 
683 
684 G4double
686  G4double momentum, G4double A)
687 {
688  G4int i, iMax = 100;
689  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
690 
691  fParticle = particle;
692  fWaveVector = momentum/hbarc;
693  fAtomicWeight = A;
694 
695  fNuclearRadius = CalculateNuclearRad(A);
696 
697  thetaMax = 10.174/fWaveVector/fNuclearRadius;
698 
699  if (thetaMax > pi) thetaMax = pi;
700 
702 
703  // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
704  norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
705 
706  norm *= G4UniformRand();
707 
708  for(i = 1; i <= iMax; i++)
709  {
710  theta1 = (i-1)*thetaMax/iMax;
711  theta2 = i*thetaMax/iMax;
712  sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
713 
714  if ( sum >= norm )
715  {
716  result = 0.5*(theta1 + theta2);
717  break;
718  }
719  }
720  if (i > iMax ) result = 0.5*(theta1 + theta2);
721 
722  G4double sigma = pi*thetaMax/iMax;
723 
724  result += G4RandGauss::shoot(0.,sigma);
725 
726  if(result < 0.) result = 0.;
727  if(result > thetaMax) result = thetaMax;
728 
729  return result;
730 }
731 
735 //
736 // Return inv momentum transfer -t > 0 from initialisation table
737 
739  G4int Z, G4int A)
740 {
741  fParticle = aParticle;
742  G4double m1 = fParticle->GetPDGMass();
743  G4double totElab = std::sqrt(m1*m1+p*p);
745  G4LorentzVector lv1(p,0.0,0.0,totElab);
746  G4LorentzVector lv(0.0,0.0,0.0,mass2);
747  lv += lv1;
748 
749  G4ThreeVector bst = lv.boostVector();
750  lv1.boost(-bst);
751 
752  G4ThreeVector p1 = lv1.vect();
753  G4double momentumCMS = p1.mag();
754 
755  G4double t = SampleTableT( aParticle, momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
756  return t;
757 }
758 
760 //
761 // Return inv momentum transfer -t > 0 from initialisation table
762 
764  G4double Z, G4double A)
765 {
766  G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
767  // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
768  G4double t = p*p*alpha; // -t !!!
769  return t;
770 }
771 
773 //
774 // Return scattering angle2 sampled in cms according to precalculated table.
775 
776 
777 G4double
779  G4double momentum, G4double Z, G4double A)
780 {
781  size_t iElement;
782  G4int iMomentum, iAngle;
783  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
784  G4double m1 = particle->GetPDGMass();
785 
786  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
787  {
788  if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
789  }
790  if ( iElement == fElementNumberVector.size() )
791  {
792  InitialiseOnFly(Z,A); // table preparation, if needed
793 
794  // iElement--;
795 
796  // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
797  // << " is not found, return zero angle" << G4endl;
798  // return 0.; // no table for this element
799  }
800  // G4cout<<"iElement = "<<iElement<<G4endl;
801 
802  fAngleTable = fAngleBank[iElement];
803 
804  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
805 
806  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
807  {
808  if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
809  }
810  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
811  if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
812 
813  // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
814 
815  if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
816  {
817  position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
818 
819  // G4cout<<"position = "<<position<<G4endl;
820 
821  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
822  {
823  if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
824  }
825  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
826 
827  // G4cout<<"iAngle = "<<iAngle<<G4endl;
828 
829  randAngle = GetScatteringAngle(iMomentum, iAngle, position);
830 
831  // G4cout<<"randAngle = "<<randAngle<<G4endl;
832  }
833  else // kinE inside between energy table edges
834  {
835  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
836  position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
837 
838  // G4cout<<"position = "<<position<<G4endl;
839 
840  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
841  {
842  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
843  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
844  }
845  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
846 
847  // G4cout<<"iAngle = "<<iAngle<<G4endl;
848 
849  theta2 = GetScatteringAngle(iMomentum, iAngle, position);
850 
851  // G4cout<<"theta2 = "<<theta2<<G4endl;
852  E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
853 
854  // G4cout<<"E2 = "<<E2<<G4endl;
855 
856  iMomentum--;
857 
858  // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
859 
860  // G4cout<<"position = "<<position<<G4endl;
861 
862  for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
863  {
864  // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
865  if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
866  }
867  if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
868 
869  theta1 = GetScatteringAngle(iMomentum, iAngle, position);
870 
871  // G4cout<<"theta1 = "<<theta1<<G4endl;
872 
873  E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
874 
875  // G4cout<<"E1 = "<<E1<<G4endl;
876 
877  W = 1.0/(E2 - E1);
878  W1 = (E2 - kinE)*W;
879  W2 = (kinE - E1)*W;
880 
881  randAngle = W1*theta1 + W2*theta2;
882 
883  // randAngle = theta2;
884  // G4cout<<"randAngle = "<<randAngle<<G4endl;
885  }
886  // G4double angle = randAngle;
887  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
888 
889  return randAngle;
890 }
891 
893 //
894 // Initialisation for given particle on fly using new element number
895 
897 {
898  fAtomicNumber = Z; // atomic number
899  fAtomicWeight = A; // number of nucleons
900 
901  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
902 
903  if( verboseLevel > 0 )
904  {
905  G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "
906  <<Z<<"; and A = "<<A<<G4endl;
907  }
908  fElementNumberVector.push_back(fAtomicNumber);
909 
910  BuildAngleTable();
911 
912  fAngleBank.push_back(fAngleTable);
913 
914  return;
915 }
916 
918 //
919 // Build for given particle and element table of momentum, angle probability.
920 // For the moment in lab system.
921 
923 {
924  G4int i, j;
925  G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
926  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
927 
929 
930  fAngleTable = new G4PhysicsTable(fEnergyBin);
931 
932  for( i = 0; i < fEnergyBin; i++)
933  {
934  kinE = fEnergyVector->GetLowEdgeEnergy(i);
935  partMom = std::sqrt( kinE*(kinE + 2*m1) );
936 
937  fWaveVector = partMom/hbarc;
938 
939  G4double kR = fWaveVector*fNuclearRadius;
940  G4double kR2 = kR*kR;
941  G4double kRmax = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
942  G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
943  // G4double kRlim = 1.2;
944  // G4double kRlim2 = kRlim*kRlim/kR2;
945 
946  alphaMax = kRmax*kRmax/kR2;
947 
948  if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
949 
950  alphaCoulomb = kRcoul*kRcoul/kR2;
951 
952  if( z )
953  {
954  a = partMom/m1; // beta*gamma for m1
955  fBeta = a/std::sqrt(1+a*a);
956  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
957  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
958  }
959  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
960 
961  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
962 
963  G4double delth = alphaMax/fAngleBin;
964 
965  sum = 0.;
966 
967  // fAddCoulomb = false;
968  fAddCoulomb = true;
969 
970  // for(j = 1; j < fAngleBin; j++)
971  for(j = fAngleBin-1; j >= 1; j--)
972  {
973  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
974  // alpha2 = angleBins->GetLowEdgeEnergy(j);
975 
976  // alpha1 = alphaMax*(j-1)/fAngleBin;
977  // alpha2 = alphaMax*( j )/fAngleBin;
978 
979  alpha1 = delth*(j-1);
980  // if(alpha1 < kRlim2) alpha1 = kRlim2;
981  alpha2 = alpha1 + delth;
982 
983  // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
984  if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
985 
986  delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
987  // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
988 
989  sum += delta;
990 
991  angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
992  // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl;
993  }
994  fAngleTable->insertAt(i,angleVector);
995 
996  // delete[] angleVector;
997  // delete[] angleBins;
998  }
999  return;
1000 }
1001 
1003 //
1004 //
1005 
1006 G4double
1008 {
1009  G4double x1, x2, y1, y2, randAngle;
1010 
1011  if( iAngle == 0 )
1012  {
1013  randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1014  // iAngle++;
1015  }
1016  else
1017  {
1018  if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1019  {
1020  iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1021  }
1022  y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1023  y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1024 
1025  x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1026  x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1027 
1028  if ( x1 == x2 ) randAngle = x2;
1029  else
1030  {
1031  if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1032  else
1033  {
1034  randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1035  }
1036  }
1037  }
1038  return randAngle;
1039 }
1040 
1041 
1042 
1044 //
1045 // Return scattering angle sampled in lab system (target at rest)
1046 
1047 
1048 
1049 G4double
1051  G4double tmass, G4double A)
1052 {
1053  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1054  G4double m1 = theParticle->GetPDGMass();
1055  G4double plab = aParticle->GetTotalMomentum();
1056  G4LorentzVector lv1 = aParticle->Get4Momentum();
1057  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1058  lv += lv1;
1059 
1060  G4ThreeVector bst = lv.boostVector();
1061  lv1.boost(-bst);
1062 
1063  G4ThreeVector p1 = lv1.vect();
1064  G4double ptot = p1.mag();
1065  G4double tmax = 4.0*ptot*ptot;
1066  G4double t = 0.0;
1067 
1068 
1069  //
1070  // Sample t
1071  //
1072 
1073  t = SampleT( theParticle, ptot, A);
1074 
1075  // NaN finder
1076  if(!(t < 0.0 || t >= 0.0))
1077  {
1078  if (verboseLevel > 0)
1079  {
1080  G4cout << "G4DiffuseElastic:WARNING: A = " << A
1081  << " mom(GeV)= " << plab/GeV
1082  << " S-wave will be sampled"
1083  << G4endl;
1084  }
1085  t = G4UniformRand()*tmax;
1086  }
1087  if(verboseLevel>1)
1088  {
1089  G4cout <<" t= " << t << " tmax= " << tmax
1090  << " ptot= " << ptot << G4endl;
1091  }
1092  // Sampling of angles in CM system
1093 
1094  G4double phi = G4UniformRand()*twopi;
1095  G4double cost = 1. - 2.0*t/tmax;
1096  G4double sint;
1097 
1098  if( cost >= 1.0 )
1099  {
1100  cost = 1.0;
1101  sint = 0.0;
1102  }
1103  else if( cost <= -1.0)
1104  {
1105  cost = -1.0;
1106  sint = 0.0;
1107  }
1108  else
1109  {
1110  sint = std::sqrt((1.0-cost)*(1.0+cost));
1111  }
1112  if (verboseLevel>1)
1113  {
1114  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1115  }
1116  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1117  v1 *= ptot;
1118  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1119 
1120  nlv1.boost(bst);
1121 
1122  G4ThreeVector np1 = nlv1.vect();
1123 
1124  // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1125 
1126  G4double theta = np1.theta();
1127 
1128  return theta;
1129 }
1130 
1132 //
1133 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1134 
1135 
1136 
1137 G4double
1139  G4double tmass, G4double thetaCMS)
1140 {
1141  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1142  G4double m1 = theParticle->GetPDGMass();
1143  // G4double plab = aParticle->GetTotalMomentum();
1144  G4LorentzVector lv1 = aParticle->Get4Momentum();
1145  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1146 
1147  lv += lv1;
1148 
1149  G4ThreeVector bst = lv.boostVector();
1150 
1151  lv1.boost(-bst);
1152 
1153  G4ThreeVector p1 = lv1.vect();
1154  G4double ptot = p1.mag();
1155 
1156  G4double phi = G4UniformRand()*twopi;
1157  G4double cost = std::cos(thetaCMS);
1158  G4double sint;
1159 
1160  if( cost >= 1.0 )
1161  {
1162  cost = 1.0;
1163  sint = 0.0;
1164  }
1165  else if( cost <= -1.0)
1166  {
1167  cost = -1.0;
1168  sint = 0.0;
1169  }
1170  else
1171  {
1172  sint = std::sqrt((1.0-cost)*(1.0+cost));
1173  }
1174  if (verboseLevel>1)
1175  {
1176  G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1177  }
1178  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1179  v1 *= ptot;
1180  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1181 
1182  nlv1.boost(bst);
1183 
1184  G4ThreeVector np1 = nlv1.vect();
1185 
1186 
1187  G4double thetaLab = np1.theta();
1188 
1189  return thetaLab;
1190 }
1192 //
1193 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1194 
1195 
1196 
1197 G4double
1199  G4double tmass, G4double thetaLab)
1200 {
1201  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1202  G4double m1 = theParticle->GetPDGMass();
1203  G4double plab = aParticle->GetTotalMomentum();
1204  G4LorentzVector lv1 = aParticle->Get4Momentum();
1205  G4LorentzVector lv(0.0,0.0,0.0,tmass);
1206 
1207  lv += lv1;
1208 
1209  G4ThreeVector bst = lv.boostVector();
1210 
1211  // lv1.boost(-bst);
1212 
1213  // G4ThreeVector p1 = lv1.vect();
1214  // G4double ptot = p1.mag();
1215 
1216  G4double phi = G4UniformRand()*twopi;
1217  G4double cost = std::cos(thetaLab);
1218  G4double sint;
1219 
1220  if( cost >= 1.0 )
1221  {
1222  cost = 1.0;
1223  sint = 0.0;
1224  }
1225  else if( cost <= -1.0)
1226  {
1227  cost = -1.0;
1228  sint = 0.0;
1229  }
1230  else
1231  {
1232  sint = std::sqrt((1.0-cost)*(1.0+cost));
1233  }
1234  if (verboseLevel>1)
1235  {
1236  G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1237  }
1238  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1239  v1 *= plab;
1240  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1241 
1242  nlv1.boost(-bst);
1243 
1244  G4ThreeVector np1 = nlv1.vect();
1245 
1246 
1247  G4double thetaCMS = np1.theta();
1248 
1249  return thetaCMS;
1250 }
1251 
1253 //
1254 // Test for given particle and element table of momentum, angle probability.
1255 // For the moment in lab system.
1256 
1258  G4double Z, G4double A)
1259 {
1260  fAtomicNumber = Z; // atomic number
1261  fAtomicWeight = A; // number of nucleons
1262  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1263 
1264 
1265 
1266  G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1267  <<Z<<"; and A = "<<A<<G4endl;
1268 
1269  fElementNumberVector.push_back(fAtomicNumber);
1270 
1271 
1272 
1273 
1274  G4int i=0, j;
1275  G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1276  G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1277  G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1278  G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1279  G4double epsilon = 0.001;
1280 
1282 
1283  fAngleTable = new G4PhysicsTable(fEnergyBin);
1284 
1285  fWaveVector = partMom/hbarc;
1286 
1287  G4double kR = fWaveVector*fNuclearRadius;
1288  G4double kR2 = kR*kR;
1289  G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1290  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1291 
1292  alphaMax = kRmax*kRmax/kR2;
1293 
1294  if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1295 
1296  alphaCoulomb = kRcoul*kRcoul/kR2;
1297 
1298  if( z )
1299  {
1300  a = partMom/m1; // beta*gamma for m1
1301  fBeta = a/std::sqrt(1+a*a);
1302  fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1303  fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1304  }
1305  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1306 
1307  // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1308 
1309 
1310  fAddCoulomb = false;
1311 
1312  for(j = 1; j < fAngleBin; j++)
1313  {
1314  // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1315  // alpha2 = angleBins->GetLowEdgeEnergy(j);
1316 
1317  alpha1 = alphaMax*(j-1)/fAngleBin;
1318  alpha2 = alphaMax*( j )/fAngleBin;
1319 
1320  if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1321 
1322  deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1323  deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1324  deltaAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1325  alpha1, alpha2,epsilon);
1326 
1327  // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1328  // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1329 
1330  sumL10 += deltaL10;
1331  sumL96 += deltaL96;
1332  sumAG += deltaAG;
1333 
1334  G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1335  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1336 
1337  angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1338  }
1339  fAngleTable->insertAt(i,angleVector);
1340  fAngleBank.push_back(fAngleTable);
1341 
1342  /*
1343  // Integral over all angle range - Bad accuracy !!!
1344 
1345  sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1346  sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1347  sumAG = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction,
1348  0., alpha2,epsilon);
1349  G4cout<<G4endl;
1350  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1351  <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1352  */
1353  return;
1354 }
1355 
1356 //
1357 //
G4double ThetaCMStoThetaLab(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaCMS)
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
Double_t y2[nxs]
Definition: Style.C:21
ThreeVector shoot(const G4int Ap, const G4int Af)
Hep3Vector boostVector() const
G4double IntegralElasticProb(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
G4double CalculateNuclearRad(G4double A)
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
Double_t y1[nxs]
Definition: Style.C:20
G4double CalculateZommerfeld(G4double beta, G4double Z1, G4double Z2)
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
Double_t x2[nxs]
Definition: Style.C:19
tuple a
Definition: test.py:11
G4double BesselJzero(G4double z)
double x() const
void InitialiseOnFly(G4double Z, G4double A)
G4double GetInvElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
const char * p
Definition: xmltok.h:285
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:449
Float_t norm
G4double G4NeutronHPJENDLHEData::G4double result
G4ParticleDefinition * GetDefinition() const
G4double GetDiffElasticSumProb(G4double theta)
G4double SampleThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double A)
G4double GetLowEdgeEnergy(size_t binNumber) const
G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
G4double BesselJone(G4double z)
G4double SampleT(const G4ParticleDefinition *aParticle, G4double p, G4double A)
int G4int
Definition: G4Types.hh:78
G4double SampleTableT(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
G4double SampleThetaLab(const G4HadProjectile *aParticle, G4double tmass, G4double A)
G4double GetCoulombElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double Z)
double z() const
G4double AdaptiveGauss(T &typeT, F f, G4double a, G4double b, G4double e)
G4double GetTotalMomentum() const
G4double CalculateAm(G4double momentum, G4double n, G4double Z)
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
static size_t GetNumberOfElements()
Definition: G4Element.cc:402
const G4ParticleDefinition * GetDefinition() const
virtual ~G4DiffuseElastic()
G4double GetDiffElasticProb(G4double theta)
tuple degree
Definition: hepunit.py:69
Float_t Z
Definition: plot.C:39
HepLorentzVector & boost(double, double, double)
Double_t x1[nxs]
Definition: Style.C:18
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
G4double GetDiffuseElasticXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A)
void TestAngleTable(const G4ParticleDefinition *theParticle, G4double partMom, G4double Z, G4double A)
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
const G4LorentzVector & Get4Momentum() const
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
G4LorentzVector Get4Momentum() const
double theta() const
G4double DampFactor(G4double z)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
int position
Definition: filter.cc:7
G4double GetInvElasticSumXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
double y() const
G4double GetDiffElasticSumProbA(G4double alpha)
tuple z
Definition: test.py:28
void insertAt(size_t, G4PhysicsVector *)
void SetMaxEnergy(const G4double anEnergy)
G4double GetDiffuseElasticSumXsc(const G4ParticleDefinition *particle, G4double theta, G4double momentum, G4double A, G4double Z)
G4double CalculateParticleBeta(const G4ParticleDefinition *particle, G4double momentum)
#define G4endl
Definition: G4ios.hh:61
G4double BesselOneByArg(G4double z)
G4double ThetaLabToThetaCMS(const G4DynamicParticle *aParticle, G4double tmass, G4double thetaLab)
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4double SampleTableThetaCMS(const G4ParticleDefinition *aParticle, G4double p, G4double Z, G4double A)
virtual G4double SampleInvariantT(const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
std::vector< G4Element * > G4ElementTable
double G4double
Definition: G4Types.hh:76
static G4ElementTable * GetElementTable()
Definition: G4Element.cc:395
G4double GetPDGCharge() const
static G4He3 * He3()
Definition: G4He3.cc:94
double mag() const
G4double GetInvCoulombElasticXsc(const G4ParticleDefinition *particle, G4double tMand, G4double momentum, G4double A, G4double Z)
void clearAndDestroy()
G4double GetIntegrandFunction(G4double theta)
G4double GetTotalMomentum() const