G4Generator2BN Class Reference

#include <G4Generator2BN.hh>

Inheritance diagram for G4Generator2BN:

G4VEmAngularDistribution

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 ()

Protected Member Functions

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

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.

00157  : G4VEmAngularDistribution("AngularGen2BN")
00158 {
00159   b = 1.2;
00160   index_min = -300;
00161   index_max = 319;
00162 
00163   // Set parameters minimum limits Ekelectron = 250 eV and kphoton = 100 eV
00164   kmin = 100*eV;
00165   Ekmin = 250*eV;
00166   kcut = 100*eV;
00167 
00168   // Increment Theta value for surface interpolation
00169   dtheta = 0.1*rad;
00170 
00171   nwarn = 0;
00172 
00173   // Construct Majorant Surface to 2BN cross-section
00174   // (we dont need this. Pre-calculated values are used instead due to performance issues
00175   // ConstructMajorantSurface();
00176 }

G4Generator2BN::~G4Generator2BN (  )  [virtual]

Definition at line 180 of file G4Generator2BN.cc.

00181 {}

Member Function Documentation

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

Definition at line 269 of file G4Generator2BN.cc.

References G4INCL::Math::pi.

Referenced by ConstructMajorantSurface(), and SampleDirection().

00270 {
00271 
00272   G4double dsdkdt_value = 0.;
00273   G4double Z = 1;
00274   // classic radius (in cm)
00275   G4double r0 = 2.82E-13;
00276   //squared classic radius (in barn)
00277   G4double r02 = r0*r0*1.0E+24;
00278 
00279   // Photon energy cannot be greater than electron kinetic energy
00280   if(kout > (Eel-electron_mass_c2)){
00281     dsdkdt_value = 0;
00282     return dsdkdt_value;
00283   }
00284 
00285      G4double E0 = Eel/electron_mass_c2;
00286      G4double k =  kout/electron_mass_c2;
00287      G4double E =  E0 - k;
00288 
00289      if(E <= 1*MeV ){                                 // Kinematic limit at 1 MeV !!
00290        dsdkdt_value = 0;
00291        return dsdkdt_value;
00292      }
00293 
00294 
00295      G4double p0 = std::sqrt(E0*E0-1);
00296      G4double p = std::sqrt(E*E-1);
00297      G4double L = std::log((E*E0-1+p*p0)/(E*E0-1-p*p0));
00298      G4double delta0 = E0 - p0*std::cos(theta);
00299      G4double epsilon = std::log((E+p)/(E-p));
00300      G4double Z2 = Z*Z;
00301      G4double sintheta2 = std::sin(theta)*std::sin(theta);
00302      G4double E02 = E0*E0;
00303      G4double E2 = E*E;
00304      G4double p02 = E0*E0-1;
00305      G4double k2 = k*k;
00306      G4double delta02 = delta0*delta0;
00307      G4double delta04 =  delta02* delta02;
00308      G4double Q = std::sqrt(p02+k2-2*k*p0*std::cos(theta));
00309      G4double Q2 = Q*Q;
00310      G4double epsilonQ = std::log((Q+p)/(Q-p));
00311 
00312 
00313      dsdkdt_value = Z2 * (r02/(8*pi*137)) * (1/k) * (p/p0) *
00314        ( (8 * (sintheta2*(2*E02+1))/(p02*delta04)) -
00315          ((2*(5*E02+2*E*E0+3))/(p02 * delta02)) -
00316          ((2*(p02-k2))/((Q2*delta02))) +
00317          ((4*E)/(p02*delta0)) +
00318          (L/(p*p0))*(
00319                  ((4*E0*sintheta2*(3*k-p02*E))/(p02*delta04)) +
00320                  ((4*E02*(E02+E2))/(p02*delta02)) +
00321                  ((2-2*(7*E02-3*E*E0+E2))/(p02*delta02)) +
00322                  (2*k*(E02+E*E0-1))/((p02*delta0))
00323                  ) -
00324          ((4*epsilon)/(p*delta0)) +
00325          ((epsilonQ)/(p*Q))*
00326          (4/delta02-(6*k/delta0)-(2*k*(p02-k2))/(Q2*delta0))
00327          );
00328 
00329 
00330 
00331      dsdkdt_value = dsdkdt_value*std::sin(theta);
00332      return dsdkdt_value;
00333 }

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

Definition at line 262 of file G4Generator2BN.cc.

Referenced by ConstructMajorantSurface().

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

void G4Generator2BN::ConstructMajorantSurface (  ) 

Definition at line 335 of file G4Generator2BN.cc.

References Calculatedsdkdt(), CalculateFkt(), G4cout, G4endl, G4InuclParticleNames::k0, and G4INCL::Math::pi.

00336 {
00337 
00338   G4double Eel;
00339   G4int vmax;
00340   G4double Ek;
00341   G4double k, theta, k0, theta0, thetamax;
00342   G4double ds, df, dsmax, dk, dt;
00343   G4double ratmin;
00344   G4double ratio = 0;
00345   G4double Vds, Vdf;
00346   G4double A, c;
00347 
00348   G4cout << "**** Constructing Majorant Surface for 2BN Distribution ****" << G4endl;
00349 
00350   if(kcut > kmin) kmin = kcut;
00351 
00352   G4int i = 0;
00353   // Loop on energy index
00354   for(G4int index = index_min; index < index_max; index++){
00355 
00356   G4double fraction = index/100.;
00357   Ek = std::pow(10.,fraction);
00358   Eel = Ek + electron_mass_c2;
00359 
00360   // find x-section maximum at k=kmin
00361   dsmax = 0.;
00362   thetamax = 0.;
00363 
00364   for(theta = 0.; theta < pi; theta = theta + dtheta){
00365 
00366     ds = Calculatedsdkdt(kmin, theta, Eel);
00367 
00368     if(ds > dsmax){
00369       dsmax = ds;
00370       thetamax = theta;
00371     }
00372   }
00373 
00374   // Compute surface paremeters at kmin
00375   if(Ek < kmin || thetamax == 0){
00376     c = 0;
00377     A = 0;
00378   }else{
00379     c = 1/(thetamax*thetamax);
00380     A = 2*std::sqrt(c)*dsmax/(std::pow(kmin,-b));
00381   }
00382 
00383   // look for correction factor to normalization at kmin 
00384   ratmin = 1.;
00385 
00386   // Volume under surfaces
00387   Vds = 0.;
00388   Vdf = 0.;
00389   k0 = 0.;
00390   theta0 = 0.;
00391 
00392   vmax = G4int(100.*std::log10(Ek/kmin));
00393 
00394   for(G4int v = 0; v < vmax; v++){
00395     G4double fractionLocal = (v/100.);
00396     k = std::pow(10.,fractionLocal)*kmin;
00397 
00398     for(theta = 0.; theta < pi; theta = theta + dtheta){
00399       dk = k - k0;
00400       dt = theta - theta0;
00401       ds = Calculatedsdkdt(k,theta, Eel);
00402       Vds = Vds + ds*dk*dt;
00403       df = CalculateFkt(k,theta, A, c);
00404       Vdf = Vdf + df*dk*dt;
00405 
00406       if(ds != 0.){
00407         if(df != 0.) ratio = df/ds;
00408       }
00409 
00410       if(ratio < ratmin && ratio != 0.){
00411         ratmin = ratio;
00412       }
00413     }
00414   }
00415 
00416 
00417   // sampling efficiency
00418   Atab[i] = A/ratmin * 1.04;
00419   ctab[i] = c;
00420 
00421 //  G4cout << Ek << " " << i << " " << index << " " << Atab[i] << " " << ctab[i] << " " << G4endl;
00422   i++;
00423   }
00424 
00425 }

G4double G4Generator2BN::GetGammaCutValue (  )  [inline]

Definition at line 84 of file G4Generator2BN.hh.

00084 {return kcut;};

G4double G4Generator2BN::GetInterpolationThetaIncrement (  )  [inline]

Definition at line 81 of file G4Generator2BN.hh.

00081 {return dtheta;};

void G4Generator2BN::PrintGeneratorInformation (  )  const

Definition at line 427 of file G4Generator2BN.cc.

References G4cout, and G4endl.

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

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 Calculatedsdkdt(), G4VEmAngularDistribution::fLocalDirection, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4DynamicParticle::GetTotalEnergy(), G4INCL::Math::pi, and G4Generator2BS::SampleDirection().

00187 {
00188   G4double Ek  = dp->GetKineticEnergy();
00189   G4double Eel = dp->GetTotalEnergy();
00190   if(Eel > 2*MeV) {
00191     return fGenerator2BS.SampleDirection(dp, out_energy, 0);
00192   }
00193 
00194   G4double k   = Eel - out_energy;
00195 
00196   G4double t;
00197   G4double cte2;
00198   G4double y, u;
00199   G4double ds;
00200   G4double dmax;
00201 
00202   G4int trials = 0;
00203 
00204   // find table index
00205   G4int index = G4int(std::log10(Ek/MeV)*100) - index_min;
00206   if(index > index_max) { index = index_max; }
00207   else if(index < 0) { index = 0; }
00208 
00209   //kmax = Ek;
00210   //kmin2 = kcut;
00211 
00212   G4double c = ctab[index];
00213   G4double A = Atab[index];
00214   if(index < index_max) { A = std::max(A, Atab[index+1]); }
00215 
00216   do{
00217     // generate k accordimg to std::pow(k,-b)
00218     trials++;
00219 
00220     // normalization constant 
00221     //  cte1 = (1-b)/(std::pow(kmax,(1-b))-std::pow(kmin2,(1-b)));
00222     //  y = G4UniformRand();
00223     //  k = std::pow(((1-b)/cte1*y+std::pow(kmin2,(1-b))),(1/(1-b)));
00224 
00225     // generate theta accordimg to theta/(1+c*std::pow(theta,2)
00226     // Normalization constant
00227     cte2 = 2*c/std::log(1+c*pi2);
00228 
00229     y = G4UniformRand();
00230     t = std::sqrt((std::exp(2*c*y/cte2)-1)/c);
00231     u = G4UniformRand();
00232 
00233     // point acceptance algorithm
00234     //fk = std::pow(k,-b);
00235     //ft = t/(1+c*t*t);
00236     dmax = A*std::pow(k,-b)*t/(1+c*t*t);
00237     ds = Calculatedsdkdt(k,t,Eel);
00238 
00239     // violation point
00240     if(ds > dmax && nwarn >= 20) {
00241       ++nwarn;
00242       G4cout << "### WARNING in G4Generator2BN: Ekin(MeV)= " << Ek/MeV 
00243              << "  D(Ekin,k)/Dmax-1= " << (ds/dmax - 1)
00244              << "  results are not reliable!" 
00245              << G4endl;
00246       if(20 == nwarn) { 
00247         G4cout << "   WARNING in G4Generator2BN is closed" << G4endl; 
00248       }
00249     }
00250 
00251   } while(u*dmax > ds || t > CLHEP::pi);
00252 
00253   G4double sint = std::sin(t);
00254   G4double phi  = twopi*G4UniformRand(); 
00255 
00256   fLocalDirection.set(sint*std::cos(phi), sint*std::sin(phi),std::cos(t));
00257   fLocalDirection.rotateUz(dp->GetMomentumDirection());
00258 
00259   return fLocalDirection;
00260 }

void G4Generator2BN::SetGammaCutValue ( G4double  cutValue  )  [inline]

Definition at line 83 of file G4Generator2BN.hh.

00083 {kcut = cutValue;};

void G4Generator2BN::SetInterpolationThetaIncrement ( G4double  increment  )  [inline]

Definition at line 80 of file G4Generator2BN.hh.

00080 {dtheta = increment;};

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:04 2013 for Geant4 by  doxygen 1.4.7