ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4XrayRayleighModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4XrayRayleighModel.cc
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 //
27 //
28 // Author: Vladimir Grichine
29 //
30 // History:
31 //
32 // 14.10.12 V.Grichine, update of xsc and angular distribution
33 // 25.05.2011 first implementation
34 
35 #include "G4XrayRayleighModel.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 
40 
41 using namespace std;
42 
44 
46 
48 
50  const G4String& nam)
51  :G4VEmModel(nam),isInitialised(false)
52 {
53  fParticleChange = 0;
54  lowEnergyLimit = 250*eV;
55  highEnergyLimit = 10.*MeV;
56  fFormFactor = 0.0;
57 
58  // SetLowEnergyLimit(lowEnergyLimit);
60  //
61  verboseLevel= 0;
62  // Verbosity scale:
63  // 0 = nothing
64  // 1 = warning for energy non-conservation
65  // 2 = details of energy budget
66  // 3 = calculation of cross sections, file openings, sampling of atoms
67  // 4 = entering in methods
68 
69  if(verboseLevel > 0)
70  {
71  G4cout << "Xray Rayleigh is constructed " << G4endl
72  << "Energy range: "
73  << lowEnergyLimit / eV << " eV - "
74  << highEnergyLimit / MeV << " MeV"
75  << G4endl;
76  }
77 }
78 
80 
82 {
83 
84 }
85 
87 
89  const G4DataVector& cuts)
90 {
91  if (verboseLevel > 3)
92  {
93  G4cout << "Calling G4XrayRayleighModel::Initialise()" << G4endl;
94  }
95 
96  InitialiseElementSelectors(particle,cuts);
97 
98 
99  if(isInitialised) return;
101  isInitialised = true;
102 
103 }
104 
106 
108  const G4ParticleDefinition*,
109  G4double gammaEnergy,
112 {
113  if (verboseLevel > 3)
114  {
115  G4cout << "Calling CrossSectionPerAtom() of G4XrayRayleighModel" << G4endl;
116  }
117  if (gammaEnergy < lowEnergyLimit || gammaEnergy > highEnergyLimit)
118  {
119  return 0.0;
120  }
121  G4double k = gammaEnergy/hbarc;
122  k *= Bohr_radius;
123  G4double p0 = 0.680654;
124  G4double p1 = -0.0224188;
125  G4double lnZ = std::log(Z);
126 
127  G4double lna = p0 + p1*lnZ;
128 
129  G4double alpha = std::exp(lna);
130 
131  G4double fo = std::pow(k, alpha);
132 
133  p0 = 3.68455;
134  p1 = -0.464806;
135  lna = p0 + p1*lnZ;
136 
137  fo *= 0.01*std::exp(lna);
138 
139  fFormFactor = fo;
140 
141  G4double b = 1. + 2.*fo;
142  G4double b2 = b*b;
143  G4double b3 = b*b2;
144 
145  G4double xsc = fCofR*Z*Z/b3;
146  xsc *= fo*fo + (1. + fo)*(1. + fo);
147 
148 
149  return xsc;
150 
151 }
152 
154 
155 void G4XrayRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
156  const G4MaterialCutsCouple* couple,
157  const G4DynamicParticle* aDPGamma,
158  G4double,
159  G4double)
160 {
161  if ( verboseLevel > 3)
162  {
163  G4cout << "Calling SampleSecondaries() of G4XrayRayleighModel" << G4endl;
164  }
165  G4double photonEnergy0 = aDPGamma->GetKineticEnergy();
166 
167  G4ParticleMomentum photonDirection0 = aDPGamma->GetMomentumDirection();
168 
169 
170  // Sample the angle of the scattered photon
171  // according to 1 + cosTheta*cosTheta distribution
172 
173  G4double cosDipole, cosTheta, sinTheta;
174  G4double c, delta, cofA, signc = 1., a, power = 1./3.;
175 
176  c = 4. - 8.*G4UniformRand();
177  a = c;
178 
179  if( c < 0. )
180  {
181  signc = -1.;
182  a = -c;
183  }
184  delta = std::sqrt(a*a+4.);
185  delta += a;
186  delta *= 0.5;
187  cofA = -signc*std::pow(delta, power);
188  cosDipole = cofA - 1./cofA;
189 
190  // select atom
191  const G4Element* elm = SelectTargetAtom(couple, aDPGamma->GetParticleDefinition(),
192  photonEnergy0,aDPGamma->GetLogKineticEnergy());
193  G4double Z = elm->GetZ();
194 
195  G4double k = photonEnergy0/hbarc;
196  k *= Bohr_radius;
197  G4double p0 = 0.680654;
198  G4double p1 = -0.0224188;
199  G4double lnZ = std::log(Z);
200 
201  G4double lna = p0 + p1*lnZ;
202 
203  G4double alpha = std::exp(lna);
204 
205  G4double fo = std::pow(k, alpha);
206 
207  p0 = 3.68455;
208  p1 = -0.464806;
209  lna = p0 + p1*lnZ;
210 
211  fo *= 0.01*pi*std::exp(lna);
212 
213 
214  G4double beta = fo/(1 + fo);
215 
216  cosTheta = (cosDipole + beta)/(1. + cosDipole*beta);
217 
218 
219  if( cosTheta > 1.) cosTheta = 1.;
220  if( cosTheta < -1.) cosTheta = -1.;
221 
222  sinTheta = std::sqrt( (1. - cosTheta)*(1. + cosTheta) );
223 
224  // Scattered photon angles. ( Z - axis along the parent photon)
225 
227  G4double dirX = sinTheta*std::cos(phi);
228  G4double dirY = sinTheta*std::sin(phi);
229  G4double dirZ = cosTheta;
230 
231  // Update G4VParticleChange for the scattered photon
232 
233  G4ThreeVector photonDirection1(dirX, dirY, dirZ);
234  photonDirection1.rotateUz(photonDirection0);
235  fParticleChange->ProposeMomentumDirection(photonDirection1);
236 
237  fParticleChange->SetProposedKineticEnergy(photonEnergy0);
238 }
239 
240