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