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

#include <G4RDGenerator2BN.hh>

Inheritance diagram for G4RDGenerator2BN:
Collaboration diagram for G4RDGenerator2BN:

Public Member Functions

 G4RDGenerator2BN (const G4String &name)
 
 ~G4RDGenerator2BN ()
 
G4double PolarAngle (const G4double initial_energy, const G4double final_energy, const G4int Z)
 
void PrintGeneratorInformation () const
 
void SetInterpolationThetaIncrement (G4double increment)
 
G4double GetInterpolationThetaIncrement ()
 
void SetGammaCutValue (G4double cutValue)
 
G4double GetGammaCutValue ()
 
void ConstructMajorantSurface ()
 
- Public Member Functions inherited from G4RDVBremAngularDistribution
 G4RDVBremAngularDistribution (const G4String &name)
 
virtual ~G4RDVBremAngularDistribution ()
 

Protected Member Functions

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

Detailed Description

Definition at line 60 of file G4RDGenerator2BN.hh.

Constructor & Destructor Documentation

G4RDGenerator2BN::G4RDGenerator2BN ( const G4String name)

Definition at line 154 of file G4RDGenerator2BN.cc.

155 {
156  b = 1.2;
157  index_min = -300;
158  index_max = 20;
159 
160  // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
161  kmin = 100*eV;
162  Ekmin = 250*eV;
163  kcut = 100*eV;
164 
165  // Increment Theta value for surface interpolation
166  dtheta = 0.1*rad;
167 
168  // Construct Majorant Surface to 2BN cross-section
169  // (we dont need this. Pre-calculated values are used instead due to performance issues
170  // ConstructMajorantSurface();
171 }
static constexpr double rad
Definition: G4SIunits.hh:149
static constexpr double eV
Definition: G4SIunits.hh:215
G4RDVBremAngularDistribution(const G4String &name)
G4RDGenerator2BN::~G4RDGenerator2BN ( )

Definition at line 175 of file G4RDGenerator2BN.cc.

176 {;}

Member Function Documentation

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

Definition at line 202 of file G4RDGenerator2BN.cc.

203 {
204 
205  G4double dsdkdt_value = 0.;
206  G4double Z = 1;
207  // classic radius (in cm)
208  G4double r0 = 2.82E-13;
209  //squared classic radius (in barn)
210  G4double r02 = r0*r0*1.0E+24;
211 
212  // Photon energy cannot be greater than electron kinetic energy
213  if(kout > (Eel-electron_mass_c2)){
214  dsdkdt_value = 0;
215  return dsdkdt_value;
216  }
217 
218  G4double E0 = Eel/electron_mass_c2;
219  G4double k = kout/electron_mass_c2;
220  G4double E = E0 - k;
221 
222  if(E <= 1*MeV ){ // Kinematic limit at 1 MeV !!
223  dsdkdt_value = 0;
224  return dsdkdt_value;
225  }
226 
227 
228  G4double p0 = std::sqrt(E0*E0-1);
229  G4double p = std::sqrt(E*E-1);
230  G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
231  G4double delta0 = E0 - p0*std::cos(theta);
232  G4double epsilon = std::log((E+p)/(E-p));
233  G4double Z2 = Z*Z;
234  G4double sintheta2 = std::sin(theta)*std::sin(theta);
235  G4double E02 = E0*E0;
236  G4double E2 = E*E;
237  G4double p02 = E0*E0-1;
238  G4double k2 = k*k;
239  G4double delta02 = delta0*delta0;
240  G4double delta04 = delta02* delta02;
241  G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
242  G4double Q2 = Q*Q;
243  G4double epsilonQ = std::log((Q+p)/(Q-p));
244 
245 
246  dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
247  ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
248  ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
249  ((2*(p02-k2))/((Q2*delta02))) +
250  ((4*E)/(p02*delta0)) +
251  (L/(p*p0))*(
252  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
253  ((4*E02*(E02+E2))/(p02*delta02)) +
254  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
255  (2*k*(E02+E*E0-1))/((p02*delta0))
256  ) -
257  ((4*epsilon)/(p*delta0)) +
258  ((epsilonQ)/(p*Q))*
259  (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
260  );
261 
262 
263 
264  dsdkdt_value = dsdkdt_value*std::sin(theta);
265  return dsdkdt_value;
266 }
const char * p
Definition: xmltok.h:285
static double Q[]
float electron_mass_c2
Definition: hepunit.py:274
static constexpr double MeV
Definition: G4SIunits.hh:214
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
static constexpr double L
Definition: G4SIunits.hh:124
double epsilon(double density, double temperature)

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 195 of file G4RDGenerator2BN.cc.

196 {
197  G4double Fkt_value = 0;
198  Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
199  return Fkt_value;
200 }
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 G4RDGenerator2BN::ConstructMajorantSurface ( )

Definition at line 268 of file G4RDGenerator2BN.cc.

269 {
270 
271  G4double Eel;
272  G4int vmax;
273  G4double Ek;
274  G4double k, theta, k0, theta0, thetamax;
275  G4double ds, df, dsmax, dk, dt;
276  G4double ratmin;
277  G4double ratio = 0;
278  G4double Vds, Vdf;
279  G4double A, c;
280 
281  G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
282 
283  if(kcut > kmin) kmin = kcut;
284 
285  G4int i = 0;
286  // Loop on energy index
287  for(G4int index = index_min; index < index_max; index++){
288 
289  G4double fraction = index/100.;
290  Ek = std::pow(10.,fraction);
291  Eel = Ek + electron_mass_c2;
292 
293  // find x-section maximum at k=kmin
294  dsmax = 0.;
295  thetamax = 0.;
296 
297  for(theta = 0.; theta < pi; theta = theta + dtheta){
298 
299  ds = Calculatedsdkdt(kmin, theta, Eel);
300 
301  if(ds > dsmax){
302  dsmax = ds;
303  thetamax = theta;
304  }
305  }
306 
307  // Compute surface paremeters at kmin
308  if(Ek < kmin || thetamax == 0){
309  c = 0;
310  A = 0;
311  }else{
312  c = 1/(thetamax*thetamax);
313  A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
314  }
315 
316  // look for correction factor to normalization at kmin
317  ratmin = 1.;
318 
319  // Volume under surfaces
320  Vds = 0.;
321  Vdf = 0.;
322  k0 = 0.;
323  theta0 = 0.;
324 
325  vmax = G4int(100.*std::log10(Ek/kmin));
326 
327  for(G4int v = 0; v < vmax; v++){
328  G4double fraction = (v/100.);
329  k = std::pow(10.,fraction)*kmin;
330 
331  for(theta = 0.; theta < pi; theta = theta + dtheta){
332  dk = k - k0;
333  dt = theta - theta0;
334  ds = Calculatedsdkdt(k,theta, Eel);
335  Vds = Vds + ds*dk*dt;
336  df = CalculateFkt(k,theta, A, c);
337  Vdf = Vdf + df*dk*dt;
338 
339  if(ds != 0.){
340  if(df != 0.) ratio = df/ds;
341  }
342 
343  if(ratio < ratmin && ratio != 0.){
344  ratmin = ratio;
345  }
346  }
347  }
348 
349 
350  // sampling efficiency
351  Atab[i] = A/ratmin * 1.04;
352  ctab[i] = c;
353 
354 // G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
355  i++;
356  }
357 
358 }
G4double Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
int G4int
Definition: G4Types.hh:78
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 G4RDGenerator2BN::Generate2BN ( G4double  Ek,
G4double  k 
) const
protected

Definition at line 360 of file G4RDGenerator2BN.cc.

361 {
362 
363  G4double Eel;
364  G4double t;
365  G4double cte2;
366  G4double y, u;
367  G4double fk, ft;
368  G4double ds;
369  G4double A2;
370  G4double A, c;
371 
372  G4int trials = 0;
373  G4int index;
374 
375  // find table index
376  index = G4int(std::log10(Ek)*100) - index_min;
377  Eel = Ek + electron_mass_c2;
378 
379  c = ctab[index];
380  A = Atab[index];
381  if(index < index_max){
382  A2 = Atab[index+1];
383  if(A2 > A) A = A2;
384  }
385 
386  do{
387  // generate k accordimg to std::pow(k,-b)
388  trials++;
389 
390  // normalization constant
391 // cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
392 // y = G4UniformRand();
393 // k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
394 
395  // generate theta accordimg to theta/(1+c*std::pow(theta,2)
396  // Normalization constant
397  cte2 = 2*c/std::log(1+c*pi2);
398 
399  y = G4UniformRand();
400  t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
401  u = G4UniformRand();
402 
403  // point acceptance algorithm
404  fk = std::pow(k,-b);
405  ft = t/(1+c*t*t);
406  ds = Calculatedsdkdt(k,t,Eel);
407 
408  // violation point
409  if(ds > (A*fk*ft)) G4cout << "WARNING IN G4RDGenerator2BN !!!" << Ek << " " << (ds-A*fk*ft)/ds << G4endl;
410 
411  }while(u*A*fk*ft > ds);
412 
413  return t;
414 
415 }
G4double Calculatedsdkdt(G4double kout, G4double theta, G4double Eel) const
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
float electron_mass_c2
Definition: hepunit.py:274
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
static constexpr double pi2
Definition: G4SIunits.hh:78

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4RDGenerator2BN::GetGammaCutValue ( )
inline

Definition at line 81 of file G4RDGenerator2BN.hh.

81 {return kcut;};
G4double G4RDGenerator2BN::GetInterpolationThetaIncrement ( )
inline

Definition at line 78 of file G4RDGenerator2BN.hh.

78 {return dtheta;};
G4double G4RDGenerator2BN::PolarAngle ( const G4double  initial_energy,
const G4double  final_energy,
const G4int  Z 
)
virtual

Implements G4RDVBremAngularDistribution.

Definition at line 180 of file G4RDGenerator2BN.cc.

183 {
184 
185  G4double theta = 0;
186 
187  G4double k = initial_energy - final_energy;
188 
189  theta = Generate2BN(initial_energy, k);
190 
191  return theta;
192 }
G4double Generate2BN(G4double Ek, G4double k) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4RDGenerator2BN::PrintGeneratorInformation ( ) const
virtual

Implements G4RDVBremAngularDistribution.

Definition at line 417 of file G4RDGenerator2BN.cc.

418 {
419  G4cout << "\n" << G4endl;
420  G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
421  G4cout << "\n" << G4endl;
422 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
void G4RDGenerator2BN::SetGammaCutValue ( G4double  cutValue)
inline

Definition at line 80 of file G4RDGenerator2BN.hh.

80 {kcut = cutValue;};
void G4RDGenerator2BN::SetInterpolationThetaIncrement ( G4double  increment)
inline

Definition at line 77 of file G4RDGenerator2BN.hh.

77 {dtheta = increment;};

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