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