Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermorePolarizedRayleighModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4LivermorePolarizedRayleighModel.cc 93810 2015-11-02 11:27:56Z gcosmo $
27 //
28 // Author: Sebastien Incerti
29 // 30 October 2008
30 // on base of G4LowEnergyPolarizedRayleigh developed by R. Capra
31 //
32 // History:
33 // --------
34 // 02 May 2009 S Incerti as V. Ivanchenko proposed in G4LivermoreRayleighModel.cc
35 //
36 // Cleanup initialisation and generation of secondaries:
37 // - apply internal high-energy limit only in constructor
38 // - do not apply low-energy limit (default is 0)
39 // - remove GetMeanFreePath method and table
40 // - remove initialisation of element selector
41 // - use G4ElementSelector
42 
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4LogLogInterpolation.hh"
47 #include "G4CompositeEMDataSet.hh"
48 
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 
51 using namespace std;
52 
53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54 
55 G4int G4LivermorePolarizedRayleighModel::maxZ = 100;
56 G4LPhysicsFreeVector* G4LivermorePolarizedRayleighModel::dataCS[] = {0};
57 G4VEMDataSet* G4LivermorePolarizedRayleighModel::formFactorData = 0;
58 
60  const G4String& nam)
61  :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
62 {
63  fParticleChange =0;
64  lowEnergyLimit = 250 * eV;
65  //SetLowEnergyLimit(lowEnergyLimit);
66  //SetHighEnergyLimit(highEnergyLimit);
67  //
68  verboseLevel= 0;
69  // Verbosity scale:
70  // 0 = nothing
71  // 1 = warning for energy non-conservation
72  // 2 = details of energy budget
73  // 3 = calculation of cross sections, file openings, sampling of atoms
74  // 4 = entering in methods
75 
76  if(verboseLevel > 0) {
77  G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
78  << "Energy range: "
79  << LowEnergyLimit() / eV << " eV - "
80  << HighEnergyLimit() / GeV << " GeV"
81  << G4endl;
82  }
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86 
88 {
89  if(IsMaster()) {
90  for(G4int i=0; i<maxZ; ++i) {
91  if(dataCS[i]) {
92  delete dataCS[i];
93  dataCS[i] = 0;
94  }
95  }
96  delete formFactorData;
97  formFactorData = 0;
98 
99  }
100 }
101 
102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
103 
105  const G4DataVector& cuts)
106 {
107 // Rayleigh process: The Quantum Theory of Radiation
108 // W. Heitler, Oxford at the Clarendon Press, Oxford (1954)
109 // Scattering function: A simple model of photon transport
110 // D.E. Cullen, Nucl. Instr. Meth. in Phys. Res. B 101 (1995) 499-510
111 // Polarization of the outcoming photon: Beam test of a prototype detector array for the PoGO astronomical hard X-ray/soft gamma-ray polarimeter
112 // T. Mizuno et al., Nucl. Instr. Meth. in Phys. Res. A 540 (2005) 158-168
113 
114  if (verboseLevel > 3)
115  G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
116 
117 
118  if(IsMaster()) {
119 
120  // Form Factor
121 
122  G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
123  G4String formFactorFile = "rayl/re-ff-";
124  formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
125  formFactorData->LoadData(formFactorFile);
126 
127  // Initialise element selector
128  InitialiseElementSelectors(particle, cuts);
129 
130  // Access to elements
131  char* path = getenv("G4LEDATA");
132  G4ProductionCutsTable* theCoupleTable =
134  G4int numOfCouples = theCoupleTable->GetTableSize();
135 
136  for(G4int i=0; i<numOfCouples; ++i)
137  {
138  const G4MaterialCutsCouple* couple =
139  theCoupleTable->GetMaterialCutsCouple(i);
140  const G4Material* material = couple->GetMaterial();
141  const G4ElementVector* theElementVector = material->GetElementVector();
142  G4int nelm = material->GetNumberOfElements();
143 
144  for (G4int j=0; j<nelm; ++j)
145  {
146  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
147  if(Z < 1) { Z = 1; }
148  else if(Z > maxZ) { Z = maxZ; }
149  if( (!dataCS[Z]) ) { ReadData(Z, path); }
150  }
151  }
152  }
153 
154  if(isInitialised) { return; }
155  fParticleChange = GetParticleChangeForGamma();
156  isInitialised = true;
157 
158 }
159 
160 
161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
162 
164  G4VEmModel* masterModel)
165 {
167 }
168 
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170 
171 void G4LivermorePolarizedRayleighModel::ReadData(size_t Z, const char* path)
172 {
173  if (verboseLevel > 1)
174  {
175  G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
176  << G4endl;
177  }
178 
179  if(dataCS[Z]) { return; }
180 
181  const char* datadir = path;
182 
183  if(!datadir)
184  {
185  datadir = getenv("G4LEDATA");
186  if(!datadir)
187  {
188  G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
190  "Environment variable G4LEDATA not defined");
191  return;
192  }
193  }
194 
195  //
196 
197  dataCS[Z] = new G4LPhysicsFreeVector();
198 
199  // Activation of spline interpolation
200  //dataCS[Z] ->SetSpline(true);
201 
202  std::ostringstream ostCS;
203  ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
204  std::ifstream finCS(ostCS.str().c_str());
205 
206  if( !finCS .is_open() )
207  {
209  ed << "G4LivermorePolarizedRayleighModel data file <" << ostCS.str().c_str()
210  << "> is not opened!" << G4endl;
211  G4Exception("G4LivermorePolarizedRayleighModel::ReadData()","em0003",FatalException,
212  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
213  return;
214  }
215  else
216  {
217  if(verboseLevel > 3) {
218  G4cout << "File " << ostCS.str()
219  << " is opened by G4LivermoreRayleighModel" << G4endl;
220  }
221  dataCS[Z]->Retrieve(finCS, true);
222  }
223  }
224 
225  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 
228  const G4ParticleDefinition*,
229  G4double GammaEnergy,
230  G4double Z, G4double,
232  {
233  if (verboseLevel > 1)
234  {
235  G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()"
236  << G4endl;
237  }
238 
239  if(GammaEnergy < lowEnergyLimit) { return 0.0; }
240 
241  G4double xs = 0.0;
242 
243  G4int intZ = G4lrint(Z);
244 
245  if(intZ < 1 || intZ > maxZ) { return xs; }
246 
247  G4LPhysicsFreeVector* pv = dataCS[intZ];
248 
249  // if element was not initialised
250  // do initialisation safely for MT mode
251  if(!pv) {
252  InitialiseForElement(0, intZ);
253  pv = dataCS[intZ];
254  if(!pv) { return xs; }
255  }
256 
257  G4int n = pv->GetVectorLength() - 1;
258  G4double e = GammaEnergy/MeV;
259  if(e >= pv->Energy(n)) {
260  xs = (*pv)[n]/(e*e);
261  } else if(e >= pv->Energy(0)) {
262  xs = pv->Value(e)/(e*e);
263  }
264 
265  /* if(verboseLevel > 0)
266  {
267  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
268  << e << G4endl;
269  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
270  G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
271  << G4endl;
272  G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
273  << G4endl;
274  G4cout << "*********************************************************"
275  << G4endl;
276  }
277  */
278 
279  return xs;
280  }
281 
282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
283 
284 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
285  const G4MaterialCutsCouple* couple,
286  const G4DynamicParticle* aDynamicGamma,
287  G4double,
288  G4double)
289 {
290  if (verboseLevel > 3)
291  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl;
292 
293  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
294 
295  if (photonEnergy0 <= lowEnergyLimit)
296  {
297  fParticleChange->ProposeTrackStatus(fStopAndKill);
298  fParticleChange->SetProposedKineticEnergy(0.);
299  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
300  return ;
301  }
302 
303  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
304 
305  // Select randomly one element in the current material
306  // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
307  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
308  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
309  G4int Z = (G4int)elm->GetZ();
310 
311  G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);
312  G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta);
313  G4double beta=GeneratePolarizationAngle();
314 
315  // incomingPhoton reference frame:
316  // z = versor parallel to the incomingPhotonDirection
317  // x = versor parallel to the incomingPhotonPolarization
318  // y = defined as z^x
319 
320  // outgoingPhoton reference frame:
321  // z' = versor parallel to the outgoingPhotonDirection
322  // x' = defined as x-x*z'z' normalized
323  // y' = defined as z'^x'
324 
325  G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit());
326  G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma));
327  G4ThreeVector y(z.cross(x));
328 
329  // z' = std::cos(phi)*std::sin(theta) x + std::sin(phi)*std::sin(theta) y + std::cos(theta) z
330  G4double xDir;
331  G4double yDir;
332  G4double zDir;
333  zDir=outcomingPhotonCosTheta;
334  xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta);
335  yDir=xDir;
336  xDir*=std::cos(outcomingPhotonPhi);
337  yDir*=std::sin(outcomingPhotonPhi);
338 
339  G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit());
340  G4ThreeVector xPrime(x.perpPart(zPrime).unit());
341  G4ThreeVector yPrime(zPrime.cross(xPrime));
342 
343  // outgoingPhotonPolarization is directed as x' std::cos(beta) + y' std::sin(beta)
344  G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta));
345 
346  fParticleChange->ProposeMomentumDirection(zPrime);
347  fParticleChange->ProposePolarization(outcomingPhotonPolarization);
348  fParticleChange->SetProposedKineticEnergy(photonEnergy0);
349 
350 }
351 
352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
353 
354 G4double G4LivermorePolarizedRayleighModel::GenerateCosTheta(G4double incomingPhotonEnergy, G4int zAtom) const
355 {
356  // d sigma k0
357  // --------- = r0^2 * pi * F^2(x, Z) * ( 2 - sin^2 theta) * std::sin (theta), x = ---- std::sin(theta/2)
358  // d theta hc
359 
360  // d sigma k0 1 - y
361  // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + y^2), x = ---- std::sqrt ( ------- ), y = std::cos(theta)
362  // d y hc 2
363 
364  // Z
365  // F(x, Z) ~ --------
366  // a + bx
367  //
368  // The time to exit from the outer loop grows as ~ k0
369  // On pcgeant2 the time is ~ 1 s for k0 ~ 1 MeV on the oxygen element. A 100 GeV
370  // event will take ~ 10 hours.
371  //
372  // On the avarage the inner loop does 1.5 iterations before exiting
373 
374  const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light);
375  //const G4VEMDataSet * formFactorData = GetScatterFunctionData();
376 
377  G4double cosTheta;
378  G4double fCosTheta;
379  G4double x;
380  G4double fValue;
381 
382  if (incomingPhotonEnergy > 5.*MeV)
383  {
384  cosTheta = 1.;
385  }
386  else
387  {
388  do
389  {
390  do
391  {
392  cosTheta = 2.*G4UniformRand()-1.;
393  fCosTheta = (1.+cosTheta*cosTheta)/2.;
394  }
395  while (fCosTheta < G4UniformRand());
396 
397  x = xFactor*std::sqrt((1.-cosTheta)/2.);
398 
399  if (x > 1.e+005)
400  fValue = formFactorData->FindValue(x, zAtom-1);
401  else
402  fValue = formFactorData->FindValue(0., zAtom-1);
403 
404  fValue/=zAtom;
405  fValue*=fValue;
406  }
407  while(fValue < G4UniformRand());
408  }
409 
410  return cosTheta;
411 }
412 
413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414 
415 G4double G4LivermorePolarizedRayleighModel::GeneratePhi(G4double cosTheta) const
416 {
417  // d sigma
418  // --------- = alpha * ( 1 - sin^2 (theta) * cos^2 (phi) )
419  // d phi
420 
421  // On the average the loop takes no more than 2 iterations before exiting
422 
423  G4double phi;
424  G4double cosPhi;
425  G4double phiProbability;
426  G4double sin2Theta;
427 
428  sin2Theta=1.-cosTheta*cosTheta;
429 
430  do
431  {
432  phi = twopi * G4UniformRand();
433  cosPhi = std::cos(phi);
434  phiProbability= 1. - sin2Theta*cosPhi*cosPhi;
435  }
436  while (phiProbability < G4UniformRand());
437 
438  return phi;
439 }
440 
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
442 
443 G4double G4LivermorePolarizedRayleighModel::GeneratePolarizationAngle(void) const
444 {
445  // Rayleigh polarization is always on the x' direction
446 
447  return 0;
448 }
449 
450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
451 
452 G4ThreeVector G4LivermorePolarizedRayleighModel::GetPhotonPolarization(const G4DynamicParticle& photon)
453 {
454 
455 // SI - From G4VLowEnergyDiscretePhotonProcess.cc
456 
457  G4ThreeVector photonMomentumDirection;
458  G4ThreeVector photonPolarization;
459 
460  photonPolarization = photon.GetPolarization();
461  photonMomentumDirection = photon.GetMomentumDirection();
462 
463  if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.)
464  {
465  // if |photonPolarization|==0. or |photonPolarization * photonDirection0| > 1e-6 * |photonPolarization ^ photonDirection0|
466  // then polarization is choosen randomly.
467 
468  G4ThreeVector e1(photonMomentumDirection.orthogonal().unit());
469  G4ThreeVector e2(photonMomentumDirection.cross(e1).unit());
470 
472 
473  e1*=std::cos(angle);
474  e2*=std::sin(angle);
475 
476  photonPolarization=e1+e2;
477  }
478  else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.)
479  {
480  // if |photonPolarization * photonDirection0| != 0.
481  // then polarization is made orthonormal;
482 
483  photonPolarization=photonPolarization.perpPart(photonMomentumDirection);
484  }
485 
486  return photonPolarization.unit();
487 }
488 
489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
490 
491  #include "G4AutoLock.hh"
492  namespace { G4Mutex LivermorePolarizedRayleighModelMutex = G4MUTEX_INITIALIZER; }
493 
495  G4int Z)
496  {
497  G4AutoLock l(&LivermorePolarizedRayleighModelMutex);
498  // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
499  // << Z << G4endl;
500  if(!dataCS[Z]) { ReadData(Z); }
501  l.unlock();
502  }
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:146
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
float h_Planck
Definition: hepunit.py:263
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4double GetZ() const
Definition: G4Element.hh:131
static G4double angle[DIM]
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
double howOrthogonal(const Hep3Vector &v) const
Definition: SpaceVector.cc:219
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
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
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
const G4ThreeVector & GetMomentumDirection() const
static constexpr double cm
Definition: G4SIunits.hh:119
void ProposePolarization(const G4ThreeVector &dir)
static constexpr double eV
Definition: G4SIunits.hh:215
Hep3Vector perpPart() const
const G4int n
G4double Energy(size_t index) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:809
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
static G4ProductionCutsTable * GetProductionCutsTable()
G4int G4Mutex
Definition: G4Threading.hh:173
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:817
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
Hep3Vector unit() const
Hep3Vector orthogonal() const
const G4ThreeVector & GetPolarization() const
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
static constexpr double GeV
Definition: G4SIunits.hh:217
tuple z
Definition: test.py:28
virtual G4bool LoadData(const G4String &fileName)=0
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4LivermorePolarizedRayleighModel(const G4ParticleDefinition *p=0, const G4String &nam="LivermorePolarizedRayleigh")
double mag() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
float c_light
Definition: hepunit.py:257
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const