ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DiffuseElastic.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 //
27 //
28 // Physics model class G4DiffuseElastic
29 //
30 //
31 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
32 //
33 // 24-May-07 V. Grichine
34 //
35 // 21.10.15 V. Grichine
36 // Bug fixed in BuildAngleTable, improving accuracy for
37 // angle bins at high energies > 50 GeV for pions.
38 //
39 
40 #include "G4DiffuseElastic.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh"
44 #include "G4NucleiProperties.hh"
45 
46 #include "Randomize.hh"
47 #include "G4Integrator.hh"
48 #include "globals.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 
52 #include "G4Proton.hh"
53 #include "G4Neutron.hh"
54 #include "G4Deuteron.hh"
55 #include "G4Alpha.hh"
56 #include "G4PionPlus.hh"
57 #include "G4PionMinus.hh"
58 
59 #include "G4Element.hh"
60 #include "G4ElementTable.hh"
61 #include "G4NistManager.hh"
62 #include "G4PhysicsTable.hh"
63 #include "G4PhysicsLogVector.hh"
64 #include "G4PhysicsFreeVector.hh"
65 
66 #include "G4Exp.hh"
67 
68 #include "G4HadronicParameters.hh"
69 
71 //
72 // Test Constructor. Just to check xsc
73 
74 
76  : G4HadronElastic("DiffuseElastic"), fParticle(0)
77 {
78  SetMinEnergy( 0.01*MeV );
80 
81  verboseLevel = 0;
82  lowEnergyRecoilLimit = 100.*keV;
83  lowEnergyLimitQ = 0.0*GeV;
84  lowEnergyLimitHE = 0.0*GeV;
85  lowestEnergyLimit = 0.0*keV;
86  plabLowLimit = 20.0*MeV;
87 
94 
95  fEnergyBin = 300; // Increased from the original 200 to have no wider log-energy-bins up to 10 PeV
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,
176  G4double A )
177 {
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,
247  G4double A, G4double Z )
248 {
250  fWaveVector = momentum/hbarc;
251  fAtomicWeight = A;
252  fAtomicNumber = Z;
254  fAddCoulomb = false;
255 
256  G4double z = particle->GetPDGCharge();
257 
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 {
682  G4double result;
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,
699  G4double A )
700 {
701  G4double result;
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
736 {
737  G4int i, iMax = 100;
738  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
739 
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
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 
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 
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 
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 //