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