G4Generator2BN Class Reference

#include <G4Generator2BN.hh>

Inheritance 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

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 156 of file G4Generator2BN.cc.

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

 : 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();
}

G4Generator2BN::~G4Generator2BN ( )

virtual

Definition at line 180 of file G4Generator2BN.cc.

{}

Member Function Documentation

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

protected

Definition at line 269 of file G4Generator2BN.cc.

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

Referenced by ConstructMajorantSurface(), and SampleDirection().

{

  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 G4Generator2BN::CalculateFkt ( G4double  k, G4double  theta, G4double  A, G4double  c  ) const

protected

Definition at line 262 of file G4Generator2BN.cc.

Referenced by ConstructMajorantSurface().

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

void G4Generator2BN::ConstructMajorantSurface

(

)

Definition at line 335 of file G4Generator2BN.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 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++;
  }

}

G4double G4Generator2BN::GetGammaCutValue ( )

inline

Definition at line 84 of file G4Generator2BN.hh.

{return kcut;};

G4double G4Generator2BN::GetInterpolationThetaIncrement ( )

inline

Definition at line 81 of file G4Generator2BN.hh.

{return dtheta;};

void G4Generator2BN::PrintGeneratorInformation

(

)

const

Definition at line 427 of file G4Generator2BN.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;
} 

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

virtual

Implements G4VEmAngularDistribution.

Definition at line 183 of file G4Generator2BN.cc.

References test::c, Calculatedsdkdt(), G4VEmAngularDistribution::fLocalDirection, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetTotalEnergy(), G4INCL::Math::max(), python.hepunit::MeV, python.hepunit::pi2, CLHEP::Hep3Vector::rotateUz(), G4Generator2BS::SampleDirection(), CLHEP::Hep3Vector::set(), and python.hepunit::twopi.

{
  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((std::exp(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;
}

void G4Generator2BN::SetGammaCutValue ( G4double  cutValue )

inline

Definition at line 83 of file G4Generator2BN.hh.

{kcut = cutValue;};

void G4Generator2BN::SetInterpolationThetaIncrement ( G4double  increment )

inline

Definition at line 80 of file G4Generator2BN.hh.

{dtheta = increment;};

---

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

• G4Generator2BN.hh
• G4Generator2BN.cc