Geant4-11
Public Member Functions | Protected Member Functions | Protected Attributes | Private Attributes | Static Private Attributes
G4Generator2BN Class Reference

#include <G4Generator2BN.hh>

Inheritance diagram for G4Generator2BN:
G4VEmAngularDistribution

Public Member Functions

void ConstructMajorantSurface ()
 
 G4Generator2BN (const G4Generator2BN &)=delete
 
 G4Generator2BN (const G4String &name="")
 
G4double GetGammaCutValue ()
 
G4double GetInterpolationThetaIncrement ()
 
const G4StringGetName () const
 
G4Generator2BNoperator= (const G4Generator2BN &right)=delete
 
void PrintGeneratorInformation () const override
 
G4ThreeVectorSampleDirection (const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=nullptr) override
 
virtual G4ThreeVectorSampleDirectionForShell (const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
 
virtual void SamplePairDirections (const G4DynamicParticle *dp, G4double elecKinEnergy, G4double posiKinEnergy, G4ThreeVector &dirElectron, G4ThreeVector &dirPositron, G4int Z=0, const G4Material *mat=nullptr)
 
void SetGammaCutValue (G4double cutValue)
 
void SetInterpolationThetaIncrement (G4double increment)
 
virtual ~G4Generator2BN ()
 

Protected Member Functions

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

Protected Attributes

G4ThreeVector fLocalDirection
 
G4bool fPolarisation
 

Private Attributes

G4double b
 
G4double dtheta
 
G4double Ekmin
 
G4Generator2BS fGenerator2BS
 
G4String fName
 
G4int index_max
 
G4int index_min
 
G4double kcut
 
G4double kmin
 
G4int nwarn
 

Static Private Attributes

static G4double Atab [320]
 
static G4double ctab [320]
 

Detailed Description

Definition at line 61 of file G4Generator2BN.hh.

Constructor & Destructor Documentation

◆ G4Generator2BN() [1/2]

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

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:129
static constexpr double eV
Definition: G4SIunits.hh:201
G4VEmAngularDistribution(const G4String &name)

References b, dtheta, Ekmin, eV, index_max, index_min, kcut, kmin, nwarn, and rad.

◆ ~G4Generator2BN()

virtual G4Generator2BN::~G4Generator2BN ( )
inlinevirtual

Definition at line 67 of file G4Generator2BN.hh.

67{;};

◆ G4Generator2BN() [2/2]

G4Generator2BN::G4Generator2BN ( const G4Generator2BN )
delete

Member Function Documentation

◆ Calculatedsdkdt()

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

Definition at line 268 of file G4Generator2BN.cc.

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

References source.hepunit::electron_mass_c2, epsilon(), LL, MeV, pi, Q, and Z.

Referenced by ConstructMajorantSurface(), and SampleDirection().

◆ CalculateFkt()

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

Definition at line 259 of file G4Generator2BN.cc.

260{
261 G4double Fkt_value = 0;
262 Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
263 return Fkt_value;
264}
const G4double A[17]

References A, and b.

Referenced by ConstructMajorantSurface().

◆ ConstructMajorantSurface()

void G4Generator2BN::ConstructMajorantSurface ( )

Definition at line 332 of file G4Generator2BN.cc.

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

References A, Atab, b, Calculatedsdkdt(), CalculateFkt(), ctab, anonymous_namespace{G4QuasiElRatios.cc}::ds, dtheta, source.hepunit::electron_mass_c2, G4cout, G4endl, index_max, index_min, G4InuclParticleNames::k0, kcut, kmin, and pi.

◆ GetGammaCutValue()

G4double G4Generator2BN::GetGammaCutValue ( )
inline

Definition at line 80 of file G4Generator2BN.hh.

80{return kcut;};

References kcut.

◆ GetInterpolationThetaIncrement()

G4double G4Generator2BN::GetInterpolationThetaIncrement ( )
inline

Definition at line 77 of file G4Generator2BN.hh.

77{return dtheta;};

References dtheta.

◆ GetName()

const G4String & G4VEmAngularDistribution::GetName ( ) const
inlineinherited

Definition at line 111 of file G4VEmAngularDistribution.hh.

112{
113 return fName;
114}

References G4VEmAngularDistribution::fName.

◆ operator=()

G4Generator2BN & G4Generator2BN::operator= ( const G4Generator2BN right)
delete

◆ PrintGeneratorInformation()

void G4Generator2BN::PrintGeneratorInformation ( ) const
overridevirtual

Reimplemented from G4VEmAngularDistribution.

Definition at line 423 of file G4Generator2BN.cc.

424{
425 G4cout << "\n" << G4endl;
426 G4cout << "Bremsstrahlung Angular Generator is 2BN Generator from 2BN Koch & Motz distribution (Rev Mod Phys 31(4), 920 (1959))" << G4endl;
427 G4cout << "\n" << G4endl;
428}

References G4cout, and G4endl.

◆ SampleDirection()

G4ThreeVector & G4Generator2BN::SampleDirection ( const G4DynamicParticle dp,
G4double  out_energy,
G4int  Z,
const G4Material mat = nullptr 
)
overridevirtual

Implements G4VEmAngularDistribution.

Definition at line 181 of file G4Generator2BN.cc.

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

References A, Atab, b, Calculatedsdkdt(), ctab, anonymous_namespace{G4QuasiElRatios.cc}::ds, fGenerator2BS, G4VEmAngularDistribution::fLocalDirection, G4cout, G4endl, G4Exp(), G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetTotalEnergy(), index_max, index_min, G4INCL::Math::max(), MeV, nwarn, CLHEP::pi, pi2, CLHEP::Hep3Vector::rotateUz(), G4Generator2BS::SampleDirection(), CLHEP::Hep3Vector::set(), and twopi.

◆ SampleDirectionForShell()

G4ThreeVector & G4VEmAngularDistribution::SampleDirectionForShell ( const G4DynamicParticle dp,
G4double  finalTotalEnergy,
G4int  Z,
G4int  shellID,
const G4Material mat 
)
virtualinherited

◆ SamplePairDirections()

void G4VEmAngularDistribution::SamplePairDirections ( const G4DynamicParticle dp,
G4double  elecKinEnergy,
G4double  posiKinEnergy,
G4ThreeVector dirElectron,
G4ThreeVector dirPositron,
G4int  Z = 0,
const G4Material mat = nullptr 
)
virtualinherited

◆ SetGammaCutValue()

void G4Generator2BN::SetGammaCutValue ( G4double  cutValue)
inline

Definition at line 79 of file G4Generator2BN.hh.

79{kcut = cutValue;};

References kcut.

◆ SetInterpolationThetaIncrement()

void G4Generator2BN::SetInterpolationThetaIncrement ( G4double  increment)
inline

Definition at line 76 of file G4Generator2BN.hh.

76{dtheta = increment;};

References dtheta.

Field Documentation

◆ Atab

G4double G4Generator2BN::Atab
staticprivate

Definition at line 94 of file G4Generator2BN.hh.

Referenced by ConstructMajorantSurface(), and SampleDirection().

◆ b

G4double G4Generator2BN::b
private

◆ ctab

G4double G4Generator2BN::ctab
staticprivate

Definition at line 95 of file G4Generator2BN.hh.

Referenced by ConstructMajorantSurface(), and SampleDirection().

◆ dtheta

G4double G4Generator2BN::dtheta
private

◆ Ekmin

G4double G4Generator2BN::Ekmin
private

Definition at line 98 of file G4Generator2BN.hh.

Referenced by G4Generator2BN().

◆ fGenerator2BS

G4Generator2BS G4Generator2BN::fGenerator2BS
private

Definition at line 93 of file G4Generator2BN.hh.

Referenced by SampleDirection().

◆ fLocalDirection

G4ThreeVector G4VEmAngularDistribution::fLocalDirection
protectedinherited

◆ fName

G4String G4VEmAngularDistribution::fName
privateinherited

Definition at line 108 of file G4VEmAngularDistribution.hh.

Referenced by G4VEmAngularDistribution::GetName().

◆ fPolarisation

G4bool G4VEmAngularDistribution::fPolarisation
protectedinherited

◆ index_max

G4int G4Generator2BN::index_max
private

Definition at line 102 of file G4Generator2BN.hh.

Referenced by ConstructMajorantSurface(), G4Generator2BN(), and SampleDirection().

◆ index_min

G4int G4Generator2BN::index_min
private

Definition at line 102 of file G4Generator2BN.hh.

Referenced by ConstructMajorantSurface(), G4Generator2BN(), and SampleDirection().

◆ kcut

G4double G4Generator2BN::kcut
private

◆ kmin

G4double G4Generator2BN::kmin
private

Definition at line 98 of file G4Generator2BN.hh.

Referenced by ConstructMajorantSurface(), and G4Generator2BN().

◆ nwarn

G4int G4Generator2BN::nwarn
private

Definition at line 103 of file G4Generator2BN.hh.

Referenced by G4Generator2BN(), and SampleDirection().


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