Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeBremsstrahlungAngular Class Reference

#include <G4PenelopeBremsstrahlungAngular.hh>

Inheritance diagram for G4PenelopeBremsstrahlungAngular:
Collaboration diagram for G4PenelopeBremsstrahlungAngular:

Public Member Functions

 G4PenelopeBremsstrahlungAngular ()
 
 ~G4PenelopeBremsstrahlungAngular ()
 
G4double PolarAngle (const G4double initial_energy, const G4double final_energy, const G4int Z)
 
G4ThreeVectorSampleDirection (const G4DynamicParticle *dp, G4double out_energy, G4int Z, const G4Material *mat=0)
 Samples the direction of the outgoing photon (in global coordinates). More...
 
void SetVerbosityLevel (G4int vl)
 Set/Get Verbosity level. More...
 
G4int GetVerbosityLevel ()
 
void Initialize ()
 
void PrepareTables (const G4Material *material, G4bool isMaster)
 Reserved for Master Model. More...
 
- Public Member Functions inherited from G4VEmAngularDistribution
 G4VEmAngularDistribution (const G4String &name)
 
virtual ~G4VEmAngularDistribution ()
 
virtual G4ThreeVectorSampleDirectionForShell (const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
 
const G4StringGetName () const
 

Additional Inherited Members

- Protected Attributes inherited from G4VEmAngularDistribution
G4ThreeVector fLocalDirection
 

Detailed Description

Definition at line 56 of file G4PenelopeBremsstrahlungAngular.hh.

Constructor & Destructor Documentation

G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular ( )

Definition at line 61 of file G4PenelopeBremsstrahlungAngular.cc.

61  :
62  G4VEmAngularDistribution("Penelope"), theEffectiveZSq(0),
63  theLorentzTables1(0),theLorentzTables2(0)
64 
65 {
66  dataRead = false;
67  verbosityLevel = 0;
68 }
G4VEmAngularDistribution(const G4String &name)
G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular ( )

Definition at line 72 of file G4PenelopeBremsstrahlungAngular.cc.

73 {
74  ClearTables();
75 }

Member Function Documentation

G4int G4PenelopeBremsstrahlungAngular::GetVerbosityLevel ( )
inline

Definition at line 78 of file G4PenelopeBremsstrahlungAngular.hh.

78 {return verbosityLevel;};
void G4PenelopeBremsstrahlungAngular::Initialize ( )

Reserved for Master Model The Initialize() method forces the cleaning of tables

Definition at line 79 of file G4PenelopeBremsstrahlungAngular.cc.

80 {
81  ClearTables();
82 }

Here is the caller graph for this function:

G4double G4PenelopeBremsstrahlungAngular::PolarAngle ( const G4double  initial_energy,
const G4double  final_energy,
const G4int  Z 
)

Old interface, backwards compatibility. Will not work in this case it will produce a G4Exception().

Definition at line 467 of file G4PenelopeBremsstrahlungAngular.cc.

470 {
471  G4cout << "WARNING: G4PenelopeBremsstrahlungAngular() does NOT support PolarAngle()" << G4endl;
472  G4cout << "Please use the alternative interface SampleDirection()" << G4endl;
473  G4Exception("G4PenelopeBremsstrahlungAngular::PolarAngle()",
474  "em0005",FatalException,"Unsupported interface");
475  return 0;
476 }
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

void G4PenelopeBremsstrahlungAngular::PrepareTables ( const G4Material material,
G4bool  isMaster 
)

Reserved for Master Model.

Definition at line 173 of file G4PenelopeBremsstrahlungAngular.cc.

174 {
175  //Unused at the moment: the G4PenelopeBremsstrahlungAngular is thread-local, so each worker
176  //builds its own version of the tables.
177  /*
178  if (!isMaster)
179  //Should not be here!
180  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareTables()",
181  "em0100",FatalException,"Worker thread in this method");
182  */
183 
184  //Check if data file has already been read
185  if (!dataRead)
186  {
187  ReadDataFile();
188  if (!dataRead)
189  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
190  "em2001",FatalException,"Unable to build interpolation table");
191  }
192 
193  if (!theLorentzTables1)
194  theLorentzTables1 = new std::map<G4double,G4PhysicsTable*>;
195  if (!theLorentzTables2)
196  theLorentzTables2 = new std::map<G4double,G4PhysicsTable*>;
197 
198  G4double Zmat = CalculateEffectiveZ(material);
199 
200  const G4int reducedEnergyGrid=21;
201  //Support arrays.
202  G4double betas[NumberofEPoints]; //betas for interpolation
203  //tables for interpolation
204  G4double Q1[NumberofEPoints][NumberofKPoints];
205  G4double Q2[NumberofEPoints][NumberofKPoints];
206  //expanded tables for interpolation
207  G4double Q1E[NumberofEPoints][reducedEnergyGrid];
208  G4double Q2E[NumberofEPoints][reducedEnergyGrid];
209  G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
210 
211  G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
212  //Interpolation in Z
213  for (i=0;i<NumberofEPoints;i++)
214  {
215  for (j=0;j<NumberofKPoints;j++)
216  {
217  G4PhysicsFreeVector* QQ1vector = new G4PhysicsFreeVector(NumberofZPoints);
218  G4PhysicsFreeVector* QQ2vector = new G4PhysicsFreeVector(NumberofZPoints);
219 
220  //fill vectors
221  for (k=0;k<NumberofZPoints;k++)
222  {
223  QQ1vector->PutValue(k,pZ[k],std::log(QQ1[k][i][j]));
224  QQ2vector->PutValue(k,pZ[k],QQ2[k][i][j]);
225  }
226 
227  QQ1vector->SetSpline(true);
228  QQ2vector->SetSpline(true);
229 
230  Q1[i][j]= G4Exp(QQ1vector->Value(Zmat));
231  Q2[i][j]=QQ2vector->Value(Zmat);
232  delete QQ1vector;
233  delete QQ2vector;
234  }
235  }
236  G4double pE[NumberofEPoints] = {1.0e-03*MeV,5.0e-03*MeV,1.0e-02*MeV,5.0e-02*MeV,
237  1.0e-01*MeV,5.0e-01*MeV};
238  G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
239  G4double ppK[reducedEnergyGrid];
240 
241  for(i=0;i<reducedEnergyGrid;i++)
242  ppK[i]=((G4double) i) * 0.05;
243 
244 
245  for(i=0;i<NumberofEPoints;i++)
246  betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
247 
248 
249  for (i=0;i<NumberofEPoints;i++)
250  {
251  for (j=0;j<NumberofKPoints;j++)
252  Q1[i][j]=Q1[i][j]/Zmat;
253  }
254 
255  //Expanded table of distribution parameters
256  for (i=0;i<NumberofEPoints;i++)
257  {
258  G4PhysicsFreeVector* Q1vector = new G4PhysicsFreeVector(NumberofKPoints);
259  G4PhysicsFreeVector* Q2vector = new G4PhysicsFreeVector(NumberofKPoints);
260 
261  for (j=0;j<NumberofKPoints;j++)
262  {
263  Q1vector->PutValue(j,pK[j],std::log(Q1[i][j])); //logarithmic
264  Q2vector->PutValue(j,pK[j],Q2[i][j]);
265  }
266 
267  for (j=0;j<reducedEnergyGrid;j++)
268  {
269  Q1E[i][j]=Q1vector->Value(ppK[j]);
270  Q2E[i][j]=Q2vector->Value(ppK[j]);
271  }
272  delete Q1vector;
273  delete Q2vector;
274  }
275  //
276  //TABLES to be stored
277  //
278  G4PhysicsTable* theTable1 = new G4PhysicsTable();
279  G4PhysicsTable* theTable2 = new G4PhysicsTable();
280  // the table will contain reducedEnergyGrid G4PhysicsFreeVectors with different
281  // values of k,
282  // Each of the G4PhysicsFreeVectors has a profile of
283  // y vs. E
284  //
285  //reserve space of the vectors.
286  for (j=0;j<reducedEnergyGrid;j++)
287  {
288  G4PhysicsFreeVector* thevec = new G4PhysicsFreeVector(NumberofEPoints);
289  theTable1->push_back(thevec);
290  G4PhysicsFreeVector* thevec2 = new G4PhysicsFreeVector(NumberofEPoints);
291  theTable2->push_back(thevec2);
292  }
293 
294  for (j=0;j<reducedEnergyGrid;j++)
295  {
296  G4PhysicsFreeVector* thevec = (G4PhysicsFreeVector*) (*theTable1)[j];
297  G4PhysicsFreeVector* thevec2 = (G4PhysicsFreeVector*) (*theTable2)[j];
298  for (i=0;i<NumberofEPoints;i++)
299  {
300  thevec->PutValue(i,betas[i],Q1E[i][j]);
301  thevec2->PutValue(i,betas[i],Q2E[i][j]);
302  }
303  thevec->SetSpline(true);
304  thevec2->SetSpline(true);
305  }
306 
307  if (theLorentzTables1 && theLorentzTables2)
308  {
309  theLorentzTables1->insert(std::make_pair(Zmat,theTable1));
310  theLorentzTables2->insert(std::make_pair(Zmat,theTable2));
311  }
312  else
313  {
315  ed << "Unable to create tables of Lorentz coefficients for " << G4endl;
316  ed << "<Z>= " << Zmat << " in G4PenelopeBremsstrahlungAngular" << G4endl;
317  delete theTable1;
318  delete theTable2;
319  G4Exception("G4PenelopeBremsstrahlungAngular::PrepareInterpolationTables()",
320  "em2005",FatalException,ed);
321  }
322 
323  return;
324 }
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
void push_back(G4PhysicsVector *)
void PutValue(size_t index, G4double energy, G4double dataValue)
int G4int
Definition: G4Types.hh:78
void SetSpline(G4bool)
float electron_mass_c2
Definition: hepunit.py:274
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Samples the direction of the outgoing photon (in global coordinates).

Implements G4VEmAngularDistribution.

Definition at line 328 of file G4PenelopeBremsstrahlungAngular.cc.

332 {
333  if (!material)
334  {
335  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
336  "em2040",FatalException,"The pointer to G4Material* is nullptr");
337  return fLocalDirection;
338  }
339 
340  //Retrieve the effective Z
341  G4double Zmat = 0;
342 
343  if (!theEffectiveZSq)
344  {
345  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
346  "em2040",FatalException,"EffectiveZ table not available");
347  return fLocalDirection;
348  }
349 
350  //found in the table: return it
351  if (theEffectiveZSq->count(material))
352  Zmat = theEffectiveZSq->find(material)->second;
353  else
354  {
355  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
356  "em2040",FatalException,"Material not found in the effectiveZ table");
357  return fLocalDirection;
358  }
359 
360  if (verbosityLevel > 0)
361  {
362  G4cout << "Effective <Z> for material : " << material->GetName() <<
363  " = " << Zmat << G4endl;
364  }
365 
366  G4double ePrimary = dp->GetKineticEnergy();
367 
368  G4double beta = std::sqrt(ePrimary*(ePrimary+2*electron_mass_c2))/
369  (ePrimary+electron_mass_c2);
370  G4double cdt = 0;
371  G4double sinTheta = 0;
372  G4double phi = 0;
373 
374  //Use a pure dipole distribution for energy above 500 keV
375  if (ePrimary > 500*keV)
376  {
377  cdt = 2.0*G4UniformRand() - 1.0;
378  if (G4UniformRand() > 0.75)
379  {
380  if (cdt<0)
381  cdt = -1.0*std::pow(-cdt,1./3.);
382  else
383  cdt = std::pow(cdt,1./3.);
384  }
385  cdt = (cdt+beta)/(1.0+beta*cdt);
386  //Get primary kinematics
387  sinTheta = std::sqrt(1. - cdt*cdt);
388  phi = twopi * G4UniformRand();
389  fLocalDirection.set(sinTheta* std::cos(phi),
390  sinTheta* std::sin(phi),
391  cdt);
392  //rotate
394  //return
395  return fLocalDirection;
396  }
397 
398  if (!(theLorentzTables1->count(Zmat)) || !(theLorentzTables2->count(Zmat)))
399  {
401  ed << "Unable to retrieve Lorentz tables for Z= " << Zmat << G4endl;
402  G4Exception("G4PenelopeBremsstrahlungAngular::SampleDirection()",
403  "em2006",FatalException,ed);
404  }
405 
406  //retrieve actual tables
407  const G4PhysicsTable* theTable1 = theLorentzTables1->find(Zmat)->second;
408  const G4PhysicsTable* theTable2 = theLorentzTables2->find(Zmat)->second;
409 
410  G4double RK=20.0*eGamma/ePrimary;
411  G4int ik=std::min((G4int) RK,19);
412 
413  G4double P10=0,P11=0,P1=0;
414  G4double P20=0,P21=0,P2=0;
415 
416  //First coefficient
417  const G4PhysicsFreeVector* v1 = (G4PhysicsFreeVector*) (*theTable1)[ik];
418  const G4PhysicsFreeVector* v2 = (G4PhysicsFreeVector*) (*theTable1)[ik+1];
419  P10 = v1->Value(beta);
420  P11 = v2->Value(beta);
421  P1=P10+(RK-(G4double) ik)*(P11-P10);
422 
423  //Second coefficient
424  const G4PhysicsFreeVector* v3 = (G4PhysicsFreeVector*) (*theTable2)[ik];
425  const G4PhysicsFreeVector* v4 = (G4PhysicsFreeVector*) (*theTable2)[ik+1];
426  P20=v3->Value(beta);
427  P21=v4->Value(beta);
428  P2=P20+(RK-(G4double) ik)*(P21-P20);
429 
430  //Sampling from the Lorenz-trasformed dipole distributions
431  P1=std::min(G4Exp(P1)/beta,1.0);
432  G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
433 
434  G4double testf=0;
435 
436  if (G4UniformRand() < P1)
437  {
438  do{
439  cdt = 2.0*G4UniformRand()-1.0;
440  testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
441  }while(testf>0);
442  }
443  else
444  {
445  do{
446  cdt = 2.0*G4UniformRand()-1.0;
447  testf=G4UniformRand()-(1.0-cdt*cdt);
448  }while(testf>0);
449  }
450  cdt = (cdt+betap)/(1.0+betap*cdt);
451 
452  //Get primary kinematics
453  sinTheta = std::sqrt(1. - cdt*cdt);
454  phi = twopi * G4UniformRand();
455  fLocalDirection.set(sinTheta* std::cos(phi),
456  sinTheta* std::sin(phi),
457  cdt);
458  //rotate
460  //return
461  return fLocalDirection;
462 
463 }
void set(double x, double y, double z)
static const G4double * P1[nN]
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
static const G4double P20[nE]
static const G4double P10[nE]
static const G4double P11[nE]
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
string material
Definition: eplot.py:19
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
float electron_mass_c2
Definition: hepunit.py:274
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static const G4double * P2[nN]
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static constexpr double keV
Definition: G4SIunits.hh:216
static const G4double P21[nE]

Here is the call graph for this function:

Here is the caller graph for this function:

void G4PenelopeBremsstrahlungAngular::SetVerbosityLevel ( G4int  vl)
inline

Set/Get Verbosity level.

Definition at line 77 of file G4PenelopeBremsstrahlungAngular.hh.

77 {verbosityLevel = vl;};

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