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