G4Generator2BN Class Reference

#include <G4Generator2BN.hh>

Inheritance diagram for G4Generator2BN:

Collaboration diagram for G4Generator2BN:

Public Member Functions
	G4Generator2BN (const G4String &name="")
virtual	~G4Generator2BN ()
virtual G4ThreeVector &	SampleDirection (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 G4ThreeVector &	SampleDirectionForShell (const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
const G4String &	GetName () const

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

Private Member Functions
G4Generator2BN &	operator= (const G4Generator2BN &right)
	G4Generator2BN (const G4Generator2BN &)

Private Attributes
G4Generator2BS	fGenerator2BS
G4double	b
G4int	index_min
G4int	index_max
G4double	kmin
G4double	Ekmin
G4double	dtheta
G4double	kcut
G4int	nwarn

Static Private Attributes
static G4double	Atab [320]
static G4double	ctab [320]

Additional Inherited Members
Protected Attributes inherited from G4VEmAngularDistribution
G4ThreeVector	fLocalDirection

Detailed Description

Definition at line 62 of file G4Generator2BN.hh.

Constructor & Destructor Documentation

G4Generator2BN() [1/2]

G4Generator2BN::G4Generator2BN

(

const G4String &

name = ""

)

Definition at line 157 of file G4Generator2BN.cc.

 : G4VEmAngularDistribution("AngularGen2BN")
{
  b = 1.2;
  index_min = -300;
  index_max = 319;

  // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
  kmin = 100*eV;
  Ekmin = 250*eV;
  kcut = 100*eV;

  // Increment Theta value for surface interpolation
  dtheta = 0.1*rad;

  nwarn = 0;

  // Construct Majorant Surface to 2BN cross-section
  // (we dont need this. Pre-calculated values are used instead due to performance issues
  // ConstructMajorantSurface();
}

Here is the caller graph for this function:

~G4Generator2BN()

G4Generator2BN::~G4Generator2BN ( )

virtual

Definition at line 181 of file G4Generator2BN.cc.

{}

G4Generator2BN() [2/2]

G4Generator2BN::G4Generator2BN ( const G4Generator2BN &  )

private

Member Function Documentation

Calculatedsdkdt()

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

protected

Definition at line 270 of file G4Generator2BN.cc.

{

  G4double dsdkdt_value = 0.;
  G4double Z = 1;
  // classic radius (in cm)
  G4double r0 = 2.82E-13;
  //squared classic radius (in barn)
  G4double r02 = r0*r0*1.0E+24;

  // Photon energy cannot be greater than electron kinetic energy
  if(kout > (Eel-electron_mass_c2)){
    dsdkdt_value = 0;
    return dsdkdt_value;
  }

     G4double E0 = Eel/electron_mass_c2;
     G4double k =  kout/electron_mass_c2;
     G4double E =  E0 - k;

     if(E <= 1*MeV ){                                 // Kinematic limit at 1 MeV !!
       dsdkdt_value = 0;
       return dsdkdt_value;
     }


     G4double p0 = std::sqrt(E0*E0-1);
     G4double p = std::sqrt(E*E-1);
     G4double LL = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
     G4double delta0 = E0 - p0*std::cos(theta);
     G4double epsilon = std::log((E+p)/(E-p));
     G4double Z2 = Z*Z;
     G4double sintheta2 = std::sin(theta)*std::sin(theta);
     G4double E02 = E0*E0;
     G4double E2 = E*E;
     G4double p02 = E0*E0-1;
     G4double k2 = k*k;
     G4double delta02 = delta0*delta0;
     G4double delta04 =  delta02* delta02;
     G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
     G4double Q2 = Q*Q;
     G4double epsilonQ = std::log((Q+p)/(Q-p));


     dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
       ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
         ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
         ((2*(p02-k2))/((Q2*delta02))) +
         ((4*E)/(p02*delta0)) +
         (LL/(p*p0))*(
                 ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
                 ((4*E02*(E02+E2))/(p02*delta02)) +
                 ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
                 (2*k*(E02+E*E0-1))/((p02*delta0))
                 ) -
         ((4*epsilon)/(p*delta0)) +
         ((epsilonQ)/(p*Q))*
         (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
         );



     dsdkdt_value = dsdkdt_value*std::sin(theta);
     return dsdkdt_value;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateFkt()

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

protected

Definition at line 263 of file G4Generator2BN.cc.

{
  G4double Fkt_value = 0;
  Fkt_value = A*std::pow(k,-b)*theta/(1+c*theta*theta);
  return Fkt_value;
}

Here is the caller graph for this function:

ConstructMajorantSurface()

void G4Generator2BN::ConstructMajorantSurface

(

)

Definition at line 336 of file G4Generator2BN.cc.

{

  G4double Eel;
  G4int vmax;
  G4double Ek;
  G4double k, theta, k0, theta0, thetamax;
  G4double ds, df, dsmax, dk, dt;
  G4double ratmin;
  G4double ratio = 0;
  G4double Vds, Vdf;
  G4double A, c;

  G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;

  if(kcut > kmin) kmin = kcut;

  G4int i = 0;
  // Loop on energy index
  for(G4int index = index_min; index < index_max; index++){

  G4double fraction = index/100.;
  Ek = std::pow(10.,fraction);
  Eel = Ek + electron_mass_c2;

  // find x-section maximum at k=kmin
  dsmax = 0.;
  thetamax = 0.;

  for(theta = 0.; theta < pi; theta = theta + dtheta){

    ds = Calculatedsdkdt(kmin, theta, Eel);

    if(ds > dsmax){
      dsmax = ds;
      thetamax = theta;
    }
  }

  // Compute surface paremeters at kmin
  if(Ek < kmin || thetamax == 0){
    c = 0;
    A = 0;
  }else{
    c = 1/(thetamax*thetamax);
    A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
  }

  // look for correction factor to normalization at kmin
  ratmin = 1.;

  // Volume under surfaces
  Vds = 0.;
  Vdf = 0.;
  k0 = 0.;
  theta0 = 0.;

  vmax = G4int(100.*std::log10(Ek/kmin));

  for(G4int v = 0; v < vmax; v++){
    G4double fractionLocal = (v/100.);
    k = std::pow(10.,fractionLocal)*kmin;

    for(theta = 0.; theta < pi; theta = theta + dtheta){
      dk = k - k0;
      dt = theta - theta0;
      ds = Calculatedsdkdt(k,theta, Eel);
      Vds = Vds + ds*dk*dt;
      df = CalculateFkt(k,theta, A, c);
      Vdf = Vdf + df*dk*dt;

      if(ds != 0.){
    if(df != 0.) ratio = df/ds;
      }

      if(ratio < ratmin && ratio != 0.){
    ratmin = ratio;
      }
    }
  }


  // sampling efficiency
  Atab[i] = A/ratmin * 1.04;
  ctab[i] = c;

//  G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
  i++;
  }

}

Here is the call graph for this function:

Here is the caller graph for this function:

GetGammaCutValue()

G4double G4Generator2BN::GetGammaCutValue ( )

inline

Definition at line 84 of file G4Generator2BN.hh.

{return kcut;};

Here is the call graph for this function:

GetInterpolationThetaIncrement()

G4double G4Generator2BN::GetInterpolationThetaIncrement ( )

inline

Definition at line 81 of file G4Generator2BN.hh.

{return dtheta;};

operator=()

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

private

Here is the caller graph for this function:

PrintGeneratorInformation()

void G4Generator2BN::PrintGeneratorInformation

(

)

const

Definition at line 428 of file G4Generator2BN.cc.

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

SampleDirection()

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.

{
  G4double Ek  = dp->GetKineticEnergy();
  G4double Eel = dp->GetTotalEnergy();
  if(Eel > 2*MeV) {
    return fGenerator2BS.SampleDirection(dp, out_energy, 0);
  }

  G4double k   = Eel - out_energy;

  G4double t;
  G4double cte2;
  G4double y, u;
  G4double ds;
  G4double dmax;

  G4int trials = 0;

  // find table index
  G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
  if(index > index_max) { index = index_max; }
  else if(index < 0) { index = 0; }

  //kmax = Ek;
  //kmin2 = kcut;

  G4double c = ctab[index];
  G4double A = Atab[index];
  if(index < index_max) { A = std::max(A, Atab[index+1]); }

  do{
    // generate k accordimg to std::pow(k,-b)
    trials++;

    // normalization constant
    //  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
    //  y = G4UniformRand();
    //  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));

    // generate theta accordimg to theta/(1+c*std::pow(theta,2)
    // Normalization constant
    cte2 = 2*c/std::log(1+c*pi2);

    y = G4UniformRand();
    t = std::sqrt((G4Exp(2*c*y/cte2)-1)/c);
    u = G4UniformRand();

    // point acceptance algorithm
    //fk = std::pow(k,-b);
    //ft = t/(1+c*t*t);
    dmax = A*std::pow(k,-b)*t/(1+c*t*t);
    ds = Calculatedsdkdt(k,t,Eel);

    // violation point
    if(ds > dmax && nwarn >= 20) {
      ++nwarn;
      G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV
         << "  D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
         << "  results are not reliable!"
         << G4endl;
      if(20 == nwarn) {
    G4cout << "   WARNING in G4Generator2BN is closed" << G4endl;
      }
    }

  } while(u*dmax > ds || t > CLHEP::pi);

  G4double sint = std::sin(t);
  G4double phi  = twopi*G4UniformRand();

  fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
  fLocalDirection.rotateUz(dp->GetMomentumDirection());

  return fLocalDirection;
}

Here is the call graph for this function:

SetGammaCutValue()

void G4Generator2BN::SetGammaCutValue ( G4double  cutValue )

inline

Definition at line 83 of file G4Generator2BN.hh.

{kcut = cutValue;};

SetInterpolationThetaIncrement()

void G4Generator2BN::SetInterpolationThetaIncrement ( G4double  increment )

inline

Definition at line 80 of file G4Generator2BN.hh.

{dtheta = increment;};

Member Data Documentation

Atab

G4double G4Generator2BN::Atab

staticprivate

Definition at line 106 of file G4Generator2BN.hh.

b

G4double G4Generator2BN::b

private

Definition at line 101 of file G4Generator2BN.hh.

ctab

G4double G4Generator2BN::ctab

staticprivate

Definition at line 107 of file G4Generator2BN.hh.

dtheta

G4double G4Generator2BN::dtheta

private

Definition at line 104 of file G4Generator2BN.hh.

Ekmin

G4double G4Generator2BN::Ekmin

private

Definition at line 103 of file G4Generator2BN.hh.

fGenerator2BS

G4Generator2BS G4Generator2BN::fGenerator2BS

private

Definition at line 99 of file G4Generator2BN.hh.

index_max

G4int G4Generator2BN::index_max

private

Definition at line 102 of file G4Generator2BN.hh.

index_min

G4int G4Generator2BN::index_min

private

Definition at line 102 of file G4Generator2BN.hh.

kcut

G4double G4Generator2BN::kcut

private

Definition at line 105 of file G4Generator2BN.hh.

kmin

G4double G4Generator2BN::kmin

private

Definition at line 103 of file G4Generator2BN.hh.

nwarn

G4int G4Generator2BN::nwarn

private

Definition at line 109 of file G4Generator2BN.hh.

---

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

• G4Generator2BN.hh
• G4Generator2BN.cc