Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Generator2BN Class Reference

#include <G4Generator2BN.hh>

Inheritance diagram for G4Generator2BN:
Collaboration diagram for G4Generator2BN:

Public Member Functions

 G4Generator2BN (const G4String &name="")
 
virtual ~G4Generator2BN ()
 
virtual G4ThreeVectorSampleDirection (const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
 
void PrintGeneratorInformation () const
 
void SetInterpolationThetaIncrement (G4double increment)
 
G4double GetInterpolationThetaIncrement ()
 
void SetGammaCutValue (G4double cutValue)
 
G4double GetGammaCutValue ()
 
void ConstructMajorantSurface ()
 
- Public Member Functions inherited from G4VEmAngularDistribution
 G4VEmAngularDistribution (const G4String &name)
 
virtual ~G4VEmAngularDistribution ()
 
virtual G4ThreeVectorSampleDirectionForShell (const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
 
const G4StringGetName () const
 

Protected Member Functions

G4double CalculateFkt (G4double k, G4double theta, G4double A, G4double c) const
 
G4double Calculatedsdkdt (G4double kout, G4double theta, G4double Eel) const
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmAngularDistribution
G4ThreeVector fLocalDirection
 

Detailed Description

Definition at line 62 of file G4Generator2BN.hh.

Constructor & Destructor Documentation

G4Generator2BN::G4Generator2BN ( const G4String name = "")

Definition at line 157 of file G4Generator2BN.cc.

158  : G4VEmAngularDistribution("AngularGen2BN")
159 {
160  b = 1.2;
161  index_min = -300;
162  index_max = 319;
163 
164  // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
165  kmin = 100*eV;
166  Ekmin = 250*eV;
167  kcut = 100*eV;
168 
169  // Increment Theta value for surface interpolation
170  dtheta = 0.1*rad;
171 
172  nwarn = 0;
173 
174  // Construct Majorant Surface to 2BN cross-section
175  // (we dont need this. Pre-calculated values are used instead due to performance issues
176  // ConstructMajorantSurface();
177 }
static constexpr double rad
Definition: G4SIunits.hh:149
tuple b
Definition: test.py:12
static constexpr double eV
Definition: G4SIunits.hh:215
G4VEmAngularDistribution(const G4String &name)
G4Generator2BN::~G4Generator2BN ( )
virtual

Definition at line 181 of file G4Generator2BN.cc.

182 {}

Member Function Documentation

G4double G4Generator2BN::Calculatedsdkdt ( G4double  kout,
G4double  theta,
G4double  Eel 
) const
protected

Definition at line 270 of file G4Generator2BN.cc.

271 {
272 
273  G4double dsdkdt_value = 0.;
274  G4double Z = 1;
275  // classic radius (in cm)
276  G4double r0 = 2.82E-13;
277  //squared classic radius (in barn)
278  G4double r02 = r0*r0*1.0E+24;
279 
280  // Photon energy cannot be greater than electron kinetic energy
281  if(kout > (Eel-electron_mass_c2)){
282  dsdkdt_value = 0;
283  return dsdkdt_value;
284  }
285 
286  G4double E0 = Eel/electron_mass_c2;
287  G4double k = kout/electron_mass_c2;
288  G4double E = E0 - k;
289 
290  if(E <= 1*MeV ){ // Kinematic limit at 1 MeV !!
291  dsdkdt_value = 0;
292  return dsdkdt_value;
293  }
294 
295 
296  G4double p0 = std::sqrt(E0*E0-1);
297  G4double p = std::sqrt(E*E-1);
298  G4double LL = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
299  G4double delta0 = E0 - p0*std::cos(theta);
300  G4double epsilon = std::log((E+p)/(E-p));
301  G4double Z2 = Z*Z;
302  G4double sintheta2 = std::sin(theta)*std::sin(theta);
303  G4double E02 = E0*E0;
304  G4double E2 = E*E;
305  G4double p02 = E0*E0-1;
306  G4double k2 = k*k;
307  G4double delta02 = delta0*delta0;
308  G4double delta04 = delta02* delta02;
309  G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
310  G4double Q2 = Q*Q;
311  G4double epsilonQ = std::log((Q+p)/(Q-p));
312 
313 
314  dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
315  ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
316  ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
317  ((2*(p02-k2))/((Q2*delta02))) +
318  ((4*E)/(p02*delta0)) +
319  (LL/(p*p0))*(
320  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
321  ((4*E02*(E02+E2))/(p02*delta02)) +
322  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
323  (2*k*(E02+E*E0-1))/((p02*delta0))
324  ) -
325  ((4*epsilon)/(p*delta0)) +
326  ((epsilonQ)/(p*Q))*
327  (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
328  );
329 
330 
331 
332  dsdkdt_value = dsdkdt_value*std::sin(theta);
333  return dsdkdt_value;
334 }
const char * p
Definition: xmltok.h:285
static double Q[]
float electron_mass_c2
Definition: hepunit.py:274
static const G4int LL[nN]
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
double epsilon(double density, double temperature)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4Generator2BN::CalculateFkt ( G4double  k,
G4double  theta,
G4double  A,
G4double  c 
) const
protected

Definition at line 263 of file G4Generator2BN.cc.

264 {
265  G4double Fkt_value = 0;
266  Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
267  return Fkt_value;
268 }
tuple b
Definition: test.py:12
double A(double temperature)
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the caller graph for this function:

void G4Generator2BN::ConstructMajorantSurface ( )

Definition at line 336 of file G4Generator2BN.cc.

337 {
338 
339  G4double Eel;
340  G4int vmax;
341  G4double Ek;
342  G4double k, theta, k0, theta0, thetamax;
343  G4double ds, df, dsmax, dk, dt;
344  G4double ratmin;
345  G4double ratio = 0;
346  G4double Vds, Vdf;
347  G4double A, c;
348 
349  G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
350 
351  if(kcut > kmin) kmin = kcut;
352 
353  G4int i = 0;
354  // Loop on energy index
355  for(G4int index = index_min; index < index_max; index++){
356 
357  G4double fraction = index/100.;
358  Ek = std::pow(10.,fraction);
359  Eel = Ek + electron_mass_c2;
360 
361  // find x-section maximum at k=kmin
362  dsmax = 0.;
363  thetamax = 0.;
364 
365  for(theta = 0.; theta < pi; theta = theta + dtheta){
366 
367  ds = Calculatedsdkdt(kmin, theta, Eel);
368 
369  if(ds > dsmax){
370  dsmax = ds;
371  thetamax = theta;
372  }
373  }
374 
375  // Compute surface paremeters at kmin
376  if(Ek < kmin || thetamax == 0){
377  c = 0;
378  A = 0;
379  }else{
380  c = 1/(thetamax*thetamax);
381  A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
382  }
383 
384  // look for correction factor to normalization at kmin
385  ratmin = 1.;
386 
387  // Volume under surfaces
388  Vds = 0.;
389  Vdf = 0.;
390  k0 = 0.;
391  theta0 = 0.;
392 
393  vmax = G4int(100.*std::log10(Ek/kmin));
394 
395  for(G4int v = 0; v < vmax; v++){
396  G4double fractionLocal = (v/100.);
397  k = std::pow(10.,fractionLocal)*kmin;
398 
399  for(theta = 0.; theta < pi; theta = theta + dtheta){
400  dk = k - k0;
401  dt = theta - theta0;
402  ds = Calculatedsdkdt(k,theta, Eel);
403  Vds = Vds + ds*dk*dt;
404  df = CalculateFkt(k,theta, A, c);
405  Vdf = Vdf + df*dk*dt;
406 
407  if(ds != 0.){
408  if(df != 0.) ratio = df/ds;
409  }
410 
411  if(ratio < ratmin && ratio != 0.){
412  ratmin = ratio;
413  }
414  }
415  }
416 
417 
418  // sampling efficiency
419  Atab[i] = A/ratmin * 1.04;
420  ctab[i] = c;
421 
422 // G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
423  i++;
424  }
425 
426 }
int G4int
Definition: G4Types.hh:78
G4double Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
tuple b
Definition: test.py:12
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
float electron_mass_c2
Definition: hepunit.py:274
tuple v
Definition: test.py:18
G4double CalculateFkt(G4double k, G4double theta, G4double A, G4double c) const
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the call graph for this function:

G4double G4Generator2BN::GetGammaCutValue ( )
inline

Definition at line 84 of file G4Generator2BN.hh.

84 {return kcut;};
G4double G4Generator2BN::GetInterpolationThetaIncrement ( )
inline

Definition at line 81 of file G4Generator2BN.hh.

81 {return dtheta;};
void G4Generator2BN::PrintGeneratorInformation ( ) const

Definition at line 428 of file G4Generator2BN.cc.

429 {
430  G4cout << "\n" << G4endl;
431  G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
432  G4cout << "\n" << G4endl;
433 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
G4ThreeVector & G4Generator2BN::SampleDirection ( const G4DynamicParticle dp,
G4double  out_energy,
G4int  Z,
const G4Material mat = 0 
)
virtual

Implements G4VEmAngularDistribution.

Definition at line 184 of file G4Generator2BN.cc.

188 {
189  G4double Ek = dp->GetKineticEnergy();
190  G4double Eel = dp->GetTotalEnergy();
191  if(Eel > 2*MeV) {
192  return fGenerator2BS.SampleDirection(dp, out_energy, 0);
193  }
194 
195  G4double k = Eel - out_energy;
196 
197  G4double t;
198  G4double cte2;
199  G4double y, u;
200  G4double ds;
201  G4double dmax;
202 
203  G4int trials = 0;
204 
205  // find table index
206  G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
207  if(index > index_max) { index = index_max; }
208  else if(index < 0) { index = 0; }
209 
210  //kmax = Ek;
211  //kmin2 = kcut;
212 
213  G4double c = ctab[index];
214  G4double A = Atab[index];
215  if(index < index_max) { A = std::max(A, Atab[index+1]); }
216 
217  do{
218  // generate k accordimg to std::pow(k,-b)
219  trials++;
220 
221  // normalization constant
222  // cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
223  // y = G4UniformRand();
224  // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
225 
226  // generate theta accordimg to theta/(1+c*std::pow(theta,2)
227  // Normalization constant
228  cte2 = 2*c/std::log(1+c*pi2);
229 
230  y = G4UniformRand();
231  t = std::sqrt((G4Exp(2*c*y/cte2)-1)/c);
232  u = G4UniformRand();
233 
234  // point acceptance algorithm
235  //fk = std::pow(k,-b);
236  //ft = t/(1+c*t*t);
237  dmax = A*std::pow(k,-b)*t/(1+c*t*t);
238  ds = Calculatedsdkdt(k,t,Eel);
239 
240  // violation point
241  if(ds > dmax && nwarn >= 20) {
242  ++nwarn;
243  G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV
244  << " D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
245  << " results are not reliable!"
246  << G4endl;
247  if(20 == nwarn) {
248  G4cout << " WARNING in G4Generator2BN is closed" << G4endl;
249  }
250  }
251 
252  } while(u*dmax > ds || t > CLHEP::pi);
253 
254  G4double sint = std::sin(t);
255  G4double phi = twopi*G4UniformRand();
256 
257  fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
259 
260  return fLocalDirection;
261 }
void set(double x, double y, double z)
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
int G4int
Definition: G4Types.hh:78
G4double Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
static constexpr double twopi
Definition: G4SIunits.hh:76
tuple b
Definition: test.py:12
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
static constexpr double pi
Definition: SystemOfUnits.h:54
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
static constexpr double pi2
Definition: G4SIunits.hh:78

Here is the call graph for this function:

void G4Generator2BN::SetGammaCutValue ( G4double  cutValue)
inline

Definition at line 83 of file G4Generator2BN.hh.

83 {kcut = cutValue;};
void G4Generator2BN::SetInterpolationThetaIncrement ( G4double  increment)
inline

Definition at line 80 of file G4Generator2BN.hh.

80 {dtheta = increment;};

The documentation for this class was generated from the following files: