Geant4_10
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 66241 2012-12-13 18:34:42Z gunter $
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 
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 using namespace std;
50 
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52 
54  const G4String& nam)
55  :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
56  crossSectionHandler(0),formFactorData(0)
57 {
58  lowEnergyLimit = 250 * eV;
59  highEnergyLimit = 100 * GeV;
60 
61  //SetLowEnergyLimit(lowEnergyLimit);
62  SetHighEnergyLimit(highEnergyLimit);
63  //
64  verboseLevel= 0;
65  // Verbosity scale:
66  // 0 = nothing
67  // 1 = warning for energy non-conservation
68  // 2 = details of energy budget
69  // 3 = calculation of cross sections, file openings, sampling of atoms
70  // 4 = entering in methods
71 
72  if(verboseLevel > 0) {
73  G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
74  << "Energy range: "
75  << lowEnergyLimit / eV << " eV - "
76  << highEnergyLimit / GeV << " GeV"
77  << G4endl;
78  }
79 }
80 
81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
82 
84 {
85  if (crossSectionHandler) delete crossSectionHandler;
86  if (formFactorData) delete formFactorData;
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 
92  const G4DataVector& cuts)
93 {
94 // Rayleigh process: The Quantum Theory of Radiation
95 // W. Heitler, Oxford at the Clarendon Press, Oxford (1954)
96 // Scattering function: A simple model of photon transport
97 // D.E. Cullen, Nucl. Instr. Meth. in Phys. Res. B 101 (1995) 499-510
98 // Polarization of the outcoming photon: Beam test of a prototype detector array for the PoGO astronomical hard X-ray/soft gamma-ray polarimeter
99 // T. Mizuno et al., Nucl. Instr. Meth. in Phys. Res. A 540 (2005) 158-168
100 
101  if (verboseLevel > 3)
102  G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
103 
104  if (crossSectionHandler)
105  {
106  crossSectionHandler->Clear();
107  delete crossSectionHandler;
108  }
109 
110  // Read data files for all materials
111 
112  crossSectionHandler = new G4CrossSectionHandler;
113  crossSectionHandler->Clear();
114  G4String crossSectionFile = "rayl/re-cs-";
115  crossSectionHandler->LoadData(crossSectionFile);
116 
117  G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
118  G4String formFactorFile = "rayl/re-ff-";
119  formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
120  formFactorData->LoadData(formFactorFile);
121 
122  InitialiseElementSelectors(particle,cuts);
123 
124  //
125  if (verboseLevel > 2)
126  G4cout << "Loaded cross section files for Livermore Polarized Rayleigh model" << G4endl;
127 
128  InitialiseElementSelectors(particle,cuts);
129 
130  if (verboseLevel > 0) {
131  G4cout << "Livermore Polarized Rayleigh model is initialized " << G4endl
132  << "Energy range: "
133  << LowEnergyLimit() / eV << " eV - "
134  << HighEnergyLimit() / GeV << " GeV"
135  << G4endl;
136  }
137 
138  if(isInitialised) return;
140  isInitialised = true;
141 }
142 
143 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
144 
146  const G4ParticleDefinition*,
147  G4double GammaEnergy,
150 {
151  if (verboseLevel > 3)
152  G4cout << "Calling CrossSectionPerAtom() of G4LivermorePolarizedRayleighModel" << G4endl;
153 
154  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
155 
156  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
157  return cs;
158 }
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
162 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
163  const G4MaterialCutsCouple* couple,
164  const G4DynamicParticle* aDynamicGamma,
165  G4double,
166  G4double)
167 {
168  if (verboseLevel > 3)
169  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl;
170 
171  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
172 
173  if (photonEnergy0 <= lowEnergyLimit)
174  {
178  return ;
179  }
180 
181  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
182 
183  // Select randomly one element in the current material
184  // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
185  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
186  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
187  G4int Z = (G4int)elm->GetZ();
188 
189  G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);
190  G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta);
191  G4double beta=GeneratePolarizationAngle();
192 
193  // incomingPhoton reference frame:
194  // z = versor parallel to the incomingPhotonDirection
195  // x = versor parallel to the incomingPhotonPolarization
196  // y = defined as z^x
197 
198  // outgoingPhoton reference frame:
199  // z' = versor parallel to the outgoingPhotonDirection
200  // x' = defined as x-x*z'z' normalized
201  // y' = defined as z'^x'
202 
203  G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit());
204  G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma));
205  G4ThreeVector y(z.cross(x));
206 
207  // z' = std::cos(phi)*std::sin(theta) x + std::sin(phi)*std::sin(theta) y + std::cos(theta) z
208  G4double xDir;
209  G4double yDir;
210  G4double zDir;
211  zDir=outcomingPhotonCosTheta;
212  xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta);
213  yDir=xDir;
214  xDir*=std::cos(outcomingPhotonPhi);
215  yDir*=std::sin(outcomingPhotonPhi);
216 
217  G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit());
218  G4ThreeVector xPrime(x.perpPart(zPrime).unit());
219  G4ThreeVector yPrime(zPrime.cross(xPrime));
220 
221  // outgoingPhotonPolarization is directed as x' std::cos(beta) + y' std::sin(beta)
222  G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta));
223 
225  fParticleChange->ProposePolarization(outcomingPhotonPolarization);
226  fParticleChange->SetProposedKineticEnergy(photonEnergy0);
227 
228 }
229 
230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231 
232 G4double G4LivermorePolarizedRayleighModel::GenerateCosTheta(G4double incomingPhotonEnergy, G4int zAtom) const
233 {
234  // d sigma k0
235  // --------- = r0^2 * pi * F^2(x, Z) * ( 2 - sin^2 theta) * std::sin (theta), x = ---- std::sin(theta/2)
236  // d theta hc
237 
238  // d sigma k0 1 - y
239  // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + y^2), x = ---- std::sqrt ( ------- ), y = std::cos(theta)
240  // d y hc 2
241 
242  // Z
243  // F(x, Z) ~ --------
244  // a + bx
245  //
246  // The time to exit from the outer loop grows as ~ k0
247  // On pcgeant2 the time is ~ 1 s for k0 ~ 1 MeV on the oxygen element. A 100 GeV
248  // event will take ~ 10 hours.
249  //
250  // On the avarage the inner loop does 1.5 iterations before exiting
251 
252  const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light);
253  //const G4VEMDataSet * formFactorData = GetScatterFunctionData();
254 
255  G4double cosTheta;
256  G4double fCosTheta;
257  G4double x;
258  G4double fValue;
259 
260  if (incomingPhotonEnergy > 5.*MeV)
261  {
262  cosTheta = 1.;
263  }
264  else
265  {
266  do
267  {
268  do
269  {
270  cosTheta = 2.*G4UniformRand()-1.;
271  fCosTheta = (1.+cosTheta*cosTheta)/2.;
272  }
273  while (fCosTheta < G4UniformRand());
274 
275  x = xFactor*std::sqrt((1.-cosTheta)/2.);
276 
277  if (x > 1.e+005)
278  fValue = formFactorData->FindValue(x, zAtom-1);
279  else
280  fValue = formFactorData->FindValue(0., zAtom-1);
281 
282  fValue/=zAtom;
283  fValue*=fValue;
284  }
285  while(fValue < G4UniformRand());
286  }
287 
288  return cosTheta;
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292 
293 G4double G4LivermorePolarizedRayleighModel::GeneratePhi(G4double cosTheta) const
294 {
295  // d sigma
296  // --------- = alpha * ( 1 - sin^2 (theta) * cos^2 (phi) )
297  // d phi
298 
299  // On the average the loop takes no more than 2 iterations before exiting
300 
301  G4double phi;
302  G4double cosPhi;
303  G4double phiProbability;
304  G4double sin2Theta;
305 
306  sin2Theta=1.-cosTheta*cosTheta;
307 
308  do
309  {
310  phi = twopi * G4UniformRand();
311  cosPhi = std::cos(phi);
312  phiProbability= 1. - sin2Theta*cosPhi*cosPhi;
313  }
314  while (phiProbability < G4UniformRand());
315 
316  return phi;
317 }
318 
319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
320 
321 G4double G4LivermorePolarizedRayleighModel::GeneratePolarizationAngle(void) const
322 {
323  // Rayleigh polarization is always on the x' direction
324 
325  return 0;
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329 
330 G4ThreeVector G4LivermorePolarizedRayleighModel::GetPhotonPolarization(const G4DynamicParticle& photon)
331 {
332 
333 // SI - From G4VLowEnergyDiscretePhotonProcess.cc
334 
335  G4ThreeVector photonMomentumDirection;
336  G4ThreeVector photonPolarization;
337 
338  photonPolarization = photon.GetPolarization();
339  photonMomentumDirection = photon.GetMomentumDirection();
340 
341  if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.)
342  {
343  // if |photonPolarization|==0. or |photonPolarization * photonDirection0| > 1e-6 * |photonPolarization ^ photonDirection0|
344  // then polarization is choosen randomly.
345 
346  G4ThreeVector e1(photonMomentumDirection.orthogonal().unit());
347  G4ThreeVector e2(photonMomentumDirection.cross(e1).unit());
348 
349  G4double angle(G4UniformRand() * twopi);
350 
351  e1*=std::cos(angle);
352  e2*=std::sin(angle);
353 
354  photonPolarization=e1+e2;
355  }
356  else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.)
357  {
358  // if |photonPolarization * photonDirection0| != 0.
359  // then polarization is made orthonormal;
360 
361  photonPolarization=photonPolarization.perpPart(photonMomentumDirection);
362  }
363 
364  return photonPolarization.unit();
365 }
366 
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
G4double GetKineticEnergy() const
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:135
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
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
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
Definition: SpaceVector.cc:237
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
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)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
Double_t y
Definition: plot.C:279
G4double FindValue(G4int Z, G4double e) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
Float_t Z
Definition: plot.C:39
const G4ThreeVector & GetMomentumDirection() const
void ProposePolarization(const G4ThreeVector &dir)
Hep3Vector perpPart() const
Hep3Vector unit() const
Hep3Vector orthogonal() const
const G4ThreeVector & GetPolarization() const
void LoadData(const G4String &dataFile)
tuple z
Definition: test.py:28
virtual G4bool LoadData(const G4String &fileName)=0
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
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:510
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121