G4RDGenerator2BN Class Reference

#include <G4RDGenerator2BN.hh>

Inheritance 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.

References python.hepunit::eV, and python.hepunit::rad.

                                                      :G4RDVBremAngularDistribution(name)
{
  b = 1.2;
  index_min = -300;
  index_max = 20;

  // 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;

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

G4RDGenerator2BN::~G4RDGenerator2BN

(

)

Definition at line 175 of file G4RDGenerator2BN.cc.

{;}

Member Function Documentation

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

protected

Definition at line 202 of file G4RDGenerator2BN.cc.

References python.hepunit::electron_mass_c2, python.hepunit::MeV, and python.hepunit::pi.

Referenced by ConstructMajorantSurface(), and Generate2BN().

{

  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 L = 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)) +
         (L/(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;
}

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

protected

Definition at line 195 of file G4RDGenerator2BN.cc.

Referenced by ConstructMajorantSurface().

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

void G4RDGenerator2BN::ConstructMajorantSurface

(

)

Definition at line 268 of file G4RDGenerator2BN.cc.

References test::c, Calculatedsdkdt(), CalculateFkt(), python.hepunit::electron_mass_c2, G4cout, G4endl, G4InuclParticleNames::k0, python.hepunit::pi, and test::v.

{

  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 fraction = (v/100.);
    k = std::pow(10.,fraction)*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++;
  }

}

G4double G4RDGenerator2BN::Generate2BN ( G4double  Ek, G4double  k  ) const

protected

Definition at line 360 of file G4RDGenerator2BN.cc.

References test::c, Calculatedsdkdt(), python.hepunit::electron_mass_c2, G4cout, G4endl, G4UniformRand, and python.hepunit::pi2.

Referenced by PolarAngle().

{

  G4double Eel;
  G4double t;
  G4double cte2;
  G4double y, u;
  G4double fk, ft;
  G4double ds;
  G4double A2;
  G4double A, c;

  G4int trials = 0;
  G4int index;

  // find table index
  index = G4int(std::log10(Ek)*100) - index_min;
  Eel = Ek + electron_mass_c2;

  c = ctab[index];
  A = Atab[index];
  if(index < index_max){
    A2 = Atab[index+1];
    if(A2 > A) A = A2;
  }

  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((std::exp(2*c*y/cte2)-1)/c);
  u = G4UniformRand();

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

  // violation point
  if(ds > (A*fk*ft)) G4cout << "WARNING IN G4RDGenerator2BN !!!" << Ek << " " <<  (ds-A*fk*ft)/ds << G4endl;

  }while(u*A*fk*ft > ds);

  return t;

}

G4double G4RDGenerator2BN::GetGammaCutValue ( )

inline

Definition at line 81 of file G4RDGenerator2BN.hh.

{return kcut;};

G4double G4RDGenerator2BN::GetInterpolationThetaIncrement ( )

inline

Definition at line 78 of file G4RDGenerator2BN.hh.

{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.

References Generate2BN().

{

  G4double theta = 0;

  G4double k = initial_energy - final_energy;

  theta = Generate2BN(initial_energy, k);

  return theta;
}

void G4RDGenerator2BN::PrintGeneratorInformation ( ) const

virtual

Implements G4RDVBremAngularDistribution.

Definition at line 417 of file G4RDGenerator2BN.cc.

References G4cout, and G4endl.

{
  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;
} 

void G4RDGenerator2BN::SetGammaCutValue ( G4double  cutValue )

inline

Definition at line 80 of file G4RDGenerator2BN.hh.

{kcut = cutValue;};

void G4RDGenerator2BN::SetInterpolationThetaIncrement ( G4double  increment )

inline

Definition at line 77 of file G4RDGenerator2BN.hh.

{dtheta = increment;};

---

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

• G4RDGenerator2BN.hh
• G4RDGenerator2BN.cc