| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class file 29 // GEANT4 Class file
30 // 30 //
31 // 31 //
32 // File name: G4BetheHeitler5DModel.cc 32 // File name: G4BetheHeitler5DModel.cc
33 // 33 //
34 // Authors: 34 // Authors:
35 // Igor Semeniouk and Denis Bernard, 35 // Igor Semeniouk and Denis Bernard,
36 // LLR, Ecole polytechnique & CNRS/IN2P3, 9112 36 // LLR, Ecole polytechnique & CNRS/IN2P3, 91128 Palaiseau, France
37 // 37 //
38 // Acknowledgement of the support of the Frenc 38 // Acknowledgement of the support of the French National Research Agency
39 // (ANR-13-BS05-0002). 39 // (ANR-13-BS05-0002).
40 // 40 //
41 // Reference: Nucl. Instrum. Meth. A 899 (2018 41 // Reference: Nucl. Instrum. Meth. A 899 (2018) 85 (arXiv:1802.08253 [hep-ph])
42 // Nucl. Instrum. Meth., A 936 (201 42 // Nucl. Instrum. Meth., A 936 (2019) 290
43 // 43 //
44 // Class Description: 44 // Class Description:
45 // 45 //
46 // Generates the conversion of a high-energy p 46 // Generates the conversion of a high-energy photon to an e+e- pair, either in the field of an
47 // atomic electron (triplet) or nucleus (nucle 47 // atomic electron (triplet) or nucleus (nuclear).
48 // Samples the five-dimensional (5D) different 48 // Samples the five-dimensional (5D) differential cross-section analytical expression:
49 // . Non polarized conversion: 49 // . Non polarized conversion:
50 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lon 50 // H.A. Bethe, W. Heitler, Proc. R. Soc. Lond. Ser. A 146 (1934) 83.
51 // . Polarized conversion: 51 // . Polarized conversion:
52 // T. H. Berlin and L. Madansky, Phys. Rev. 52 // T. H. Berlin and L. Madansky, Phys. Rev. 78 (1950) 623,
53 // M. M. May, Phys. Rev. 84 (1951) 265, 53 // M. M. May, Phys. Rev. 84 (1951) 265,
54 // J. M. Jauch and F. Rohrlich, The theory o 54 // J. M. Jauch and F. Rohrlich, The theory of photons and electrons, 1976.
55 // 55 //
56 // All the above expressions are named "Bethe- 56 // All the above expressions are named "Bethe-Heitler" here.
57 // 57 //
58 // Bethe & Heitler, put in Feynman diagram par 58 // Bethe & Heitler, put in Feynman diagram parlance, compute only the two dominant diagrams of
59 // the first order Born development, which is 59 // the first order Born development, which is an excellent approximation for nuclear conversion
60 // and for high-energy triplet conversion. 60 // and for high-energy triplet conversion.
61 // 61 //
62 // Only the linear polarisation of the incomin 62 // Only the linear polarisation of the incoming photon takes part in these expressions.
63 // The circular polarisation of the incoming p 63 // The circular polarisation of the incoming photon does not (take part) and no polarisation
64 // is transfered to the final leptons. 64 // is transfered to the final leptons.
65 // 65 //
66 // In case conversion takes place in the field 66 // In case conversion takes place in the field of an isolated nucleus or electron, the bare
67 // Bethe-Heitler expression is used. 67 // Bethe-Heitler expression is used.
68 // 68 //
69 // In case the nucleus or the electron are par 69 // In case the nucleus or the electron are part of an atom, the screening of the target field
70 // by the other electrons of the atom is descr 70 // by the other electrons of the atom is described by a simple form factor, function of q2:
71 // . nuclear: N.F. Mott, H.S.W. Massey, The Th 71 // . nuclear: N.F. Mott, H.S.W. Massey, The Theory of Atomic Collisions, 1934.
72 // . triplet: J.A. Wheeler and W.E. Lamb, Phys 72 // . triplet: J.A. Wheeler and W.E. Lamb, Phys. Rev. 55 (1939) 858.
73 // 73 //
74 // The nuclear form factor that affects the pr 74 // The nuclear form factor that affects the probability of very large-q2 events, is not considered.
75 // 75 //
76 // In principle the code is valid from thresho 76 // In principle the code is valid from threshold, that is from 2 * m_e c^2 for nuclear and from
77 // 4 * m_e c^2 for triplet, up to infinity, wh 77 // 4 * m_e c^2 for triplet, up to infinity, while in pratice the divergence of the differential
78 // cross section at small q2 and, at high-ener 78 // cross section at small q2 and, at high-energy, at small polar angle, make it break down at
79 // some point that depends on machine precisio 79 // some point that depends on machine precision.
80 // 80 //
81 // Very-high-energy (above a few tens of TeV) 81 // Very-high-energy (above a few tens of TeV) LPM suppression effects in the normalized differential
82 // cross-section are not considered. 82 // cross-section are not considered.
83 // 83 //
84 // The 5D differential cross section is sample 84 // The 5D differential cross section is sampled without any high-energy nor small
85 // angle approximation(s). 85 // angle approximation(s).
86 // The generation is strictly energy-momentum 86 // The generation is strictly energy-momentum conserving when all particles in the final state
87 // are taken into account, that is, including 87 // are taken into account, that is, including the recoiling target.
88 // (In contrast with the BH expressions taken 88 // (In contrast with the BH expressions taken at face values, for which the electron energy is
89 // taken to be EMinus = GammaEnergy - EPlus) 89 // taken to be EMinus = GammaEnergy - EPlus)
90 // 90 //
91 // Tests include the examination of 1D distrib 91 // Tests include the examination of 1D distributions: see TestEm15
92 // 92 //
93 // Total cross sections are not computed (we i 93 // Total cross sections are not computed (we inherit from other classes).
94 // We just convert a photon on a target when a 94 // We just convert a photon on a target when asked to do so.
95 // 95 //
96 // Pure nuclear, pure triplet and 1/Z triplet/ 96 // Pure nuclear, pure triplet and 1/Z triplet/nuclear mixture can be generated.
97 // 97 //
98 // ------------------------------------------- 98 // -------------------------------------------------------------------
99 99
100 #include "G4BetheHeitler5DModel.hh" 100 #include "G4BetheHeitler5DModel.hh"
101 #include "G4EmParameters.hh" 101 #include "G4EmParameters.hh"
102 102
103 #include "G4PhysicalConstants.hh" 103 #include "G4PhysicalConstants.hh"
104 #include "G4SystemOfUnits.hh" 104 #include "G4SystemOfUnits.hh"
105 #include "G4Electron.hh" 105 #include "G4Electron.hh"
106 #include "G4Positron.hh" 106 #include "G4Positron.hh"
107 #include "G4Gamma.hh" 107 #include "G4Gamma.hh"
>> 108 #include "G4MuonPlus.hh"
>> 109 #include "G4MuonMinus.hh"
108 #include "G4IonTable.hh" 110 #include "G4IonTable.hh"
109 #include "G4NucleiProperties.hh" 111 #include "G4NucleiProperties.hh"
110 112
111 #include "Randomize.hh" 113 #include "Randomize.hh"
112 #include "G4ParticleChangeForGamma.hh" 114 #include "G4ParticleChangeForGamma.hh"
113 #include "G4Pow.hh" 115 #include "G4Pow.hh"
114 #include "G4Log.hh" 116 #include "G4Log.hh"
115 #include "G4Exp.hh" 117 #include "G4Exp.hh"
116 118
117 #include "G4LorentzVector.hh" 119 #include "G4LorentzVector.hh"
118 #include "G4ThreeVector.hh" 120 #include "G4ThreeVector.hh"
119 #include "G4RotationMatrix.hh" 121 #include "G4RotationMatrix.hh"
120 122
121 #include <cassert> 123 #include <cassert>
122 124
>> 125 // // Q : Use enum G4EmProcessSubType hire ?
>> 126 // enum G45DConversionMode
>> 127 // {
>> 128 // kEPair, kMuPair
>> 129 // };
>> 130
123 const G4int kEPair = 0; 131 const G4int kEPair = 0;
124 const G4int kMuPair = 1; 132 const G4int kMuPair = 1;
125 133
126 134
127 //....oooOO0OOooo........oooOO0OOooo........oo 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128 136
129 G4BetheHeitler5DModel::G4BetheHeitler5DModel(c 137 G4BetheHeitler5DModel::G4BetheHeitler5DModel(const G4ParticleDefinition* pd,
130 c 138 const G4String& nam)
131 : G4PairProductionRelModel(pd, nam), << 139 : G4PairProductionRelModel(pd, nam),fVerbose(1),fConversionType(0),
>> 140 iraw(false),
132 fLepton1(G4Electron::Definition()),fLepton 141 fLepton1(G4Electron::Definition()),fLepton2(G4Positron::Definition()),
133 fTheMuPlus(nullptr),fTheMuMinus(nullptr), <<
134 fVerbose(1), <<
135 fConversionType(0), <<
136 fConvMode(kEPair), 142 fConvMode(kEPair),
137 iraw(false) << 143 fTheMuPlus(G4MuonPlus::Definition()),fTheMuMinus(G4MuonMinus::Definition())
138 { 144 {
139 theIonTable = G4IonTable::GetIonTable(); 145 theIonTable = G4IonTable::GetIonTable();
140 //Q: Do we need this on Model 146 //Q: Do we need this on Model
141 SetLowEnergyLimit(2*fTheElectron->GetPDGMass 147 SetLowEnergyLimit(2*fTheElectron->GetPDGMass());
142 } 148 }
143 149
144 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145 151
146 G4BetheHeitler5DModel::~G4BetheHeitler5DModel( << 152 G4BetheHeitler5DModel::~G4BetheHeitler5DModel()
>> 153 {}
147 154
148 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
149 156
150 void G4BetheHeitler5DModel::Initialise(const G 157 void G4BetheHeitler5DModel::Initialise(const G4ParticleDefinition* part,
151 const G4DataVector& vec) 158 const G4DataVector& vec)
152 { 159 {
153 G4PairProductionRelModel::Initialise(part, v 160 G4PairProductionRelModel::Initialise(part, vec);
154 161
155 G4EmParameters* theManager = G4EmParameters: 162 G4EmParameters* theManager = G4EmParameters::Instance();
156 // place to initialise model parameters 163 // place to initialise model parameters
157 // Verbosity levels: ( Can redefine as neede 164 // Verbosity levels: ( Can redefine as needed, but some consideration )
158 // 0 = nothing 165 // 0 = nothing
159 // > 2 print results 166 // > 2 print results
160 // > 3 print rejection warning from transfor 167 // > 3 print rejection warning from transformation (fix bug from gammaray .. )
161 // > 4 print photon direction & polarisation 168 // > 4 print photon direction & polarisation
162 fVerbose = theManager->Verbose(); 169 fVerbose = theManager->Verbose();
163 fConversionType = theManager->GetConversionT << 170 fConversionType = theManager->GetConversionType();
164 //////////////////////////////////////////// 171 //////////////////////////////////////////////////////////////
165 // iraw : 172 // iraw :
166 // true : isolated electron or nucleus 173 // true : isolated electron or nucleus.
167 // false : inside atom -> screening for 174 // false : inside atom -> screening form factor
168 iraw = theManager->OnIsolated(); 175 iraw = theManager->OnIsolated();
169 // G4cout << "BH5DModel::Initialise verbose 176 // G4cout << "BH5DModel::Initialise verbose " << fVerbose
170 // << " isolated " << iraw << " ctype "<< 177 // << " isolated " << iraw << " ctype "<< fConversionType << G4endl;
171 178
172 //Q: Do we need this on Model 179 //Q: Do we need this on Model
173 // The Leptons defined via SetLeptonPair(..) 180 // The Leptons defined via SetLeptonPair(..) method
174 SetLowEnergyLimit(2*CLHEP::electron_mass_c2) 181 SetLowEnergyLimit(2*CLHEP::electron_mass_c2);
>> 182
>> 183 if (fConvMode == kEPair) {
>> 184 assert(fLepton1->GetPDGEncoding() == fTheElectron->GetPDGEncoding()) ;
>> 185 if (fVerbose > 3)
>> 186 G4cout << "BH5DModel::Initialise conversion to e+ e-" << G4endl;
>> 187 }
>> 188
>> 189 if (fConvMode == kMuPair) {
>> 190 assert(fLepton1->GetPDGEncoding() == fTheMuMinus->GetPDGEncoding()) ;
>> 191 if (fVerbose > 3)
>> 192 G4cout << "BH5DModel::Initialise conversion to mu+ mu-" << G4endl;
>> 193 }
175 } 194 }
176 195
177 //....oooOO0OOooo........oooOO0OOooo........oo 196 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
178 197
179 void G4BetheHeitler5DModel::SetLeptonPair(cons 198 void G4BetheHeitler5DModel::SetLeptonPair(const G4ParticleDefinition* p1,
180 const G4ParticleDefinition* p2) 199 const G4ParticleDefinition* p2)
181 { 200 {
182 G4int pdg1 = p1->GetPDGEncoding(); << 201 // Lepton1 - nagative charged particle
183 G4int pdg2 = p2->GetPDGEncoding(); << 202 if ( p1->GetPDGEncoding() < 0 ){
184 G4int pdg = std::abs(pdg1); << 203 if ( p1->GetPDGEncoding() ==
185 if ( pdg1 != -pdg2 || (pdg != 11 && pdg != 1 << 204 G4Positron::Definition()->GetPDGEncoding() ) {
186 G4ExceptionDescription ed; <<
187 ed << " Wrong pair of leptons: " << p1->Ge <<
188 << " and " << p1->GetParticleName(); <<
189 G4Exception("G4BetheHeitler5DModel::SetLep <<
190 FatalErrorInArgument, ed, ""); <<
191 } else { <<
192 if ( pdg == 11 ) { <<
193 SetConversionMode(kEPair); 205 SetConversionMode(kEPair);
194 if( pdg1 == 11 ) { << 206 fLepton1 = p2;
195 fLepton1 = p1; << 207 fLepton2 = p1;
196 fLepton2 = p2; << 208 // if (fVerbose)
197 } else { <<
198 fLepton1 = p2; <<
199 fLepton2 = p1; <<
200 } <<
201 if (fVerbose > 0) <<
202 G4cout << "G4BetheHeitler5DModel::SetLeptonP 209 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
203 << G4endl; 210 << G4endl;
>> 211 } else if ( p1->GetPDGEncoding() ==
>> 212 G4MuonPlus::Definition()->GetPDGEncoding() ) {
>> 213 SetConversionMode(kMuPair);
>> 214 fLepton1 = p2;
>> 215 fLepton2 = p1;
>> 216 // if (fVerbose)
>> 217 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
>> 218 << G4endl;
204 } else { 219 } else {
>> 220 // Exception
>> 221 G4ExceptionDescription ed;
>> 222 ed << "Model not applicable to particle(s) "
>> 223 << p1->GetParticleName() << ", "
>> 224 << p2->GetParticleName();
>> 225 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
>> 226 FatalException, ed);
>> 227 }
>> 228 } else {
>> 229 if ( p1->GetPDGEncoding() ==
>> 230 G4Electron::Definition()->GetPDGEncoding() ) {
>> 231 SetConversionMode(kEPair);
>> 232 fLepton1 = p1;
>> 233 fLepton2 = p2;
>> 234 // if (fVerbose)
>> 235 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to e+ e-"
>> 236 << G4endl;
>> 237 } else if ( p1->GetPDGEncoding() ==
>> 238 G4MuonMinus::Definition()->GetPDGEncoding() ) {
205 SetConversionMode(kMuPair); 239 SetConversionMode(kMuPair);
206 if( pdg1 == 13 ) { << 240 fLepton1 = p1;
207 fLepton1 = p1; << 241 fLepton2 = p2;
208 fLepton2 = p2; << 242 // if (fVerbose)
209 } else { <<
210 fLepton1 = p2; <<
211 fLepton2 = p1; <<
212 } <<
213 fTheMuPlus = fLepton2; <<
214 fTheMuMinus= fLepton1; <<
215 if (fVerbose > 0) <<
216 G4cout << "G4BetheHeitler5DModel::SetLeptonP 243 G4cout << "G4BetheHeitler5DModel::SetLeptonPair conversion to mu+ mu-"
217 << G4endl; 244 << G4endl;
>> 245 } else {
>> 246 // Exception
>> 247 G4ExceptionDescription ed;
>> 248 ed << "Model not applicable to particle(s) "
>> 249 << p1->GetParticleName() << ", "
>> 250 << p2->GetParticleName();
>> 251 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0002",
>> 252 FatalException, ed);
218 } 253 }
219 } 254 }
>> 255 if ( fLepton1->GetPDGEncoding() != fLepton2->GetAntiPDGEncoding() ) {
>> 256 G4Exception("G4BetheHeitler5DModel::SetLeptonPair","em0007",
>> 257 FatalErrorInArgument, "pair must be particle, antiparticle ");
>> 258 G4cerr << "BH5DModel::SetLeptonPair BAD paricle/anti particle pair"
>> 259 << fLepton1->GetParticleName() << ", "
>> 260 << fLepton2->GetParticleName() << G4endl;
>> 261 }
220 } 262 }
221 263
222 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
223 265
224 G4double G4BetheHeitler5DModel::MaxDiffCrossSe 266 G4double G4BetheHeitler5DModel::MaxDiffCrossSection(const G4double* par,
225 267 G4double Z,
226 268 G4double e,
227 269 G4double loge) const
228 { 270 {
229 const G4double Q = e/par[9]; 271 const G4double Q = e/par[9];
230 return par[0] * G4Exp((par[2]+loge*par[4])*l 272 return par[0] * G4Exp((par[2]+loge*par[4])*loge)
231 / (par[1]+ G4Exp(par[3]*loge)+G4Exp(p 273 / (par[1]+ G4Exp(par[3]*loge)+G4Exp(par[5]*loge))
232 * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/ 274 * (1+par[7]*G4Exp(par[8]*G4Log(Z))*Q/(1+Q));
233 } 275 }
234 276
235 //....oooOO0OOooo........oooOO0OOooo........oo 277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236 278
237 void 279 void
238 G4BetheHeitler5DModel::SampleSecondaries(std:: 280 G4BetheHeitler5DModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
239 const 281 const G4MaterialCutsCouple* couple,
240 const 282 const G4DynamicParticle* aDynamicGamma,
241 G4dou 283 G4double, G4double)
242 { 284 {
243 // MeV 285 // MeV
244 static const G4double ElectronMass = CLHEP 286 static const G4double ElectronMass = CLHEP::electron_mass_c2;
245 287
246 const G4double LeptonMass = fLepton1->GetPDG 288 const G4double LeptonMass = fLepton1->GetPDGMass();
247 const G4double LeptonMass2 = LeptonMass*Lep 289 const G4double LeptonMass2 = LeptonMass*LeptonMass;
248 290
249 static const G4double alpha0 = CLHEP 291 static const G4double alpha0 = CLHEP::fine_structure_const;
250 // mm 292 // mm
251 static const G4double r0 = CLHEP 293 static const G4double r0 = CLHEP::classic_electr_radius;
252 // mbarn 294 // mbarn
253 static const G4double r02 = r0*r0 295 static const G4double r02 = r0*r0*1.e+25;
254 static const G4double twoPi = CLHEP 296 static const G4double twoPi = CLHEP::twopi;
255 static const G4double factor = alpha 297 static const G4double factor = alpha0 * r02 / (twoPi*twoPi);
256 // static const G4double factor1 = p 298 // static const G4double factor1 = pow((6.0 * pi),(1.0/3.0))/(8.*alpha0*ElectronMass);
257 static const G4double factor1 = 2.661 299 static const G4double factor1 = 2.66134007899/(8.*alpha0*ElectronMass);
258 // 300 //
259 G4double PairInvMassMin = 2.*LeptonMass; 301 G4double PairInvMassMin = 2.*LeptonMass;
260 G4double TrThreshold = 2.0 * ( (LeptonMass2 302 G4double TrThreshold = 2.0 * ( (LeptonMass2)/ElectronMass + LeptonMass);
261 303
262 // 304 //
263 static const G4double nu[2][10] = { 305 static const G4double nu[2][10] = {
264 //electron 306 //electron
265 { 0.0227436, 0.0582046, 3.0322675, 2.82750 307 { 0.0227436, 0.0582046, 3.0322675, 2.8275065, -0.0034004,
266 1.1212766, 1.8989468, 68.3492750, 0.0211 308 1.1212766, 1.8989468, 68.3492750, 0.0211186, 14.4},
267 //muon 309 //muon
268 {0.67810E-06, 0.86037E+05, 2.0008395, 1.67 310 {0.67810E-06, 0.86037E+05, 2.0008395, 1.6739719, -0.0057279,
269 1.4222, 0.0, 263230.0, 0.0521, 51.1338} 311 1.4222, 0.0, 263230.0, 0.0521, 51.1338}
270 }; 312 };
271 static const G4double tr[2][10] = { 313 static const G4double tr[2][10] = {
272 //electron 314 //electron
273 { 0.0332350, 4.3942537, 2.8515925, 2.6351 315 { 0.0332350, 4.3942537, 2.8515925, 2.6351695, -0.0031510,
274 1.5737305, 1.8104647, 20.6434021, -0.027 316 1.5737305, 1.8104647, 20.6434021, -0.0272586, 28.9},
275 //muon 317 //muon
276 {0.10382E-03, 0.14408E+17, 4.1368679, 3.26 318 {0.10382E-03, 0.14408E+17, 4.1368679, 3.2662121, -0.0163091,
277 0.0000, 0.0, 0.0, 0.0000, 1.0000} 319 0.0000, 0.0, 0.0, 0.0000, 1.0000}
278 }; 320 };
279 // 321 //
280 static const G4double para[2][3][2] = { 322 static const G4double para[2][3][2] = {
281 //electron 323 //electron
282 { {11., -16.},{-1.17, -2.95},{-2., -0.5} } 324 { {11., -16.},{-1.17, -2.95},{-2., -0.5} },
283 //muon 325 //muon
284 { {17.5, 1.},{-1.17, -2.95},{2., 6.} } 326 { {17.5, 1.},{-1.17, -2.95},{2., 6.} }
285 }; 327 };
286 // 328 //
287 static const G4double correctionIndex = 1.4; 329 static const G4double correctionIndex = 1.4;
288 // 330 //
289 const G4double GammaEnergy = aDynamicGamma- 331 const G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
290 // Protection, Will not be true tot cross se 332 // Protection, Will not be true tot cross section = 0
291 if ( GammaEnergy <= PairInvMassMin) { return 333 if ( GammaEnergy <= PairInvMassMin) { return; }
292 334
293 const G4double GammaEnergy2 = GammaEnergy*Ga 335 const G4double GammaEnergy2 = GammaEnergy*GammaEnergy;
294 336
295 //////////////////////////////////////////// 337 //////////////////////////////////////////////////////////////
296 const G4ParticleMomentum GammaDirection = 338 const G4ParticleMomentum GammaDirection =
297 aDynamicGamma->GetMomentumDirection(); 339 aDynamicGamma->GetMomentumDirection();
298 G4ThreeVector GammaPolarization = aDynamicGa 340 G4ThreeVector GammaPolarization = aDynamicGamma->GetPolarization();
299 341
300 // The protection polarization perpendicular 342 // The protection polarization perpendicular to the direction vector,
301 // as it done in G4LivermorePolarizedGammaCo 343 // as it done in G4LivermorePolarizedGammaConversionModel,
302 // assuming Direction is unitary vector 344 // assuming Direction is unitary vector
303 // (projection to plane) p_proj = p - (p o 345 // (projection to plane) p_proj = p - (p o d)/(d o d) x d
304 if ( GammaPolarization.howOrthogonal(GammaDi 346 if ( GammaPolarization.howOrthogonal(GammaDirection) != 0) {
305 GammaPolarization -= GammaPolarization.dot 347 GammaPolarization -= GammaPolarization.dot(GammaDirection) * GammaDirection;
306 } 348 }
307 // End of Protection 349 // End of Protection
308 // 350 //
309 const G4double GammaPolarizationMag = GammaP 351 const G4double GammaPolarizationMag = GammaPolarization.mag();
310 << 352
311 //////////////////////////////////////////// 353 //////////////////////////////////////////////////////////////
312 // target element 354 // target element
313 // select randomly one element constituting 355 // select randomly one element constituting the material
314 const G4Element* anElement = SelectTargetAt 356 const G4Element* anElement = SelectTargetAtom(couple, fTheGamma, GammaEnergy,
315 aDyna 357 aDynamicGamma->GetLogKineticEnergy() );
316 // Atomic number 358 // Atomic number
317 const G4int Z = anElement->GetZasInt() 359 const G4int Z = anElement->GetZasInt();
318 const G4int A = SelectIsotopeNumber(an 360 const G4int A = SelectIsotopeNumber(anElement);
319 const G4double iZ13 = 1./anElement->GetIonis 361 const G4double iZ13 = 1./anElement->GetIonisation()->GetZ3();
320 const G4double targetMass = G4NucleiProperti 362 const G4double targetMass = G4NucleiProperties::GetNuclearMass(A, Z);
321 363
322 const G4double NuThreshold = 2.0 * ( (Lept 364 const G4double NuThreshold = 2.0 * ( (LeptonMass2)/targetMass + LeptonMass);
323 // No conversion possible below nuclear thre 365 // No conversion possible below nuclear threshold
324 if ( GammaEnergy <= NuThreshold) { return; } 366 if ( GammaEnergy <= NuThreshold) { return; }
325 367
326 CLHEP::HepRandomEngine* rndmEngine = G4Rando 368 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
327 369
328 // itriplet : true -- triplet, false -- nucl 370 // itriplet : true -- triplet, false -- nuclear.
329 G4bool itriplet = false; 371 G4bool itriplet = false;
330 if (fConversionType == 1) { 372 if (fConversionType == 1) {
331 itriplet = false; 373 itriplet = false;
332 } else if (fConversionType == 2) { 374 } else if (fConversionType == 2) {
333 itriplet = true; 375 itriplet = true;
334 if ( GammaEnergy <= TrThreshold ) return; 376 if ( GammaEnergy <= TrThreshold ) return;
335 } else if ( GammaEnergy > TrThreshold ) { 377 } else if ( GammaEnergy > TrThreshold ) {
336 // choose triplet or nuclear from a triple 378 // choose triplet or nuclear from a triplet/nuclear=1/Z
337 // total cross section ratio. 379 // total cross section ratio.
338 // approximate at low energies ! 380 // approximate at low energies !
339 if(rndmEngine->flat()*(Z+1) < 1.) { 381 if(rndmEngine->flat()*(Z+1) < 1.) {
340 itriplet = true; 382 itriplet = true;
341 } 383 }
342 } 384 }
343 385
344 // 386 //
345 const G4double RecoilMass = itriplet ? Elec 387 const G4double RecoilMass = itriplet ? ElectronMass : targetMass;
346 const G4double RecoilMass2 = RecoilMass*Reco 388 const G4double RecoilMass2 = RecoilMass*RecoilMass;
347 const G4double sCMS = 2.*RecoilMass*G 389 const G4double sCMS = 2.*RecoilMass*GammaEnergy + RecoilMass2;
348 const G4double sCMSPlusRM2 = sCMS + RecoilMa 390 const G4double sCMSPlusRM2 = sCMS + RecoilMass2;
349 const G4double sqrts = std::sqrt(sCMS) 391 const G4double sqrts = std::sqrt(sCMS);
350 const G4double isqrts2 = 1./(2.*sqrts); 392 const G4double isqrts2 = 1./(2.*sqrts);
351 // 393 //
352 const G4double PairInvMassMax = sqrts-Reco 394 const G4double PairInvMassMax = sqrts-RecoilMass;
353 const G4double PairInvMassRange = PairInvMas 395 const G4double PairInvMassRange = PairInvMassMax/PairInvMassMin;
354 const G4double lnPairInvMassRange = G4Log(Pa 396 const G4double lnPairInvMassRange = G4Log(PairInvMassRange);
355 397
356 // initial state. Defines z axis of "0" fram 398 // initial state. Defines z axis of "0" frame as along photon propagation.
357 // Since CMS(0., 0., GammaEnergy, GammaEnerg 399 // Since CMS(0., 0., GammaEnergy, GammaEnergy+RecoilMass) set some constants
358 const G4double betaCMS = G4LorentzVector(0.0 400 const G4double betaCMS = G4LorentzVector(0.0,0.0,GammaEnergy,GammaEnergy+RecoilMass).beta();
359 401
360 // maximum value of pdf 402 // maximum value of pdf
361 const G4double EffectiveZ = iraw ? 0.5 : Z; 403 const G4double EffectiveZ = iraw ? 0.5 : Z;
362 const G4double Threshold = itriplet ? TrThr 404 const G4double Threshold = itriplet ? TrThreshold : NuThreshold;
363 const G4double AvailableEnergy = GammaEne 405 const G4double AvailableEnergy = GammaEnergy - Threshold;
364 const G4double LogAvailableEnergy = G4Log(Av 406 const G4double LogAvailableEnergy = G4Log(AvailableEnergy);
365 // 407 //
366 const G4double MaxDiffCross = itriplet 408 const G4double MaxDiffCross = itriplet
367 ? MaxDiffCrossSection(tr[fConvMode], 409 ? MaxDiffCrossSection(tr[fConvMode],
368 EffectiveZ, AvailableEnergy, LogAvaila 410 EffectiveZ, AvailableEnergy, LogAvailableEnergy)
369 : MaxDiffCrossSection(nu[fConvMode], 411 : MaxDiffCrossSection(nu[fConvMode],
370 EffectiveZ, AvailableEnergy, LogAva 412 EffectiveZ, AvailableEnergy, LogAvailableEnergy);
371 // 413 //
372 // 50% safety marging factor 414 // 50% safety marging factor
373 const G4double ymax = 1.5 * MaxDiffCross; 415 const G4double ymax = 1.5 * MaxDiffCross;
374 // x1 bounds 416 // x1 bounds
375 const G4double xu1 = (LogAvailableEnergy > 417 const G4double xu1 = (LogAvailableEnergy > para[fConvMode][2][0])
376 ? para[fConvMode][0][0] + 418 ? para[fConvMode][0][0] +
377 para[fConvMode][1][0]*LogAvailableEner 419 para[fConvMode][1][0]*LogAvailableEnergy
378 : para[fConvMode][0][0] 420 : para[fConvMode][0][0] +
379 para[fConvMode][2][0]*para[fConvMode][ 421 para[fConvMode][2][0]*para[fConvMode][1][0];
380 const G4double xl1 = (LogAvailableEnergy > 422 const G4double xl1 = (LogAvailableEnergy > para[fConvMode][2][1])
381 ? para[fConvMode][0][1] 423 ? para[fConvMode][0][1] +
382 para[fConvMode][1][1]*LogAvailableEner 424 para[fConvMode][1][1]*LogAvailableEnergy
383 : para[fConvMode][0][1] 425 : para[fConvMode][0][1] +
384 para[fConvMode][2][1]*para[fConvMode][ 426 para[fConvMode][2][1]*para[fConvMode][1][1];
385 // 427 //
386 G4LorentzVector Recoil; 428 G4LorentzVector Recoil;
387 G4LorentzVector LeptonPlus; 429 G4LorentzVector LeptonPlus;
388 G4LorentzVector LeptonMinus; 430 G4LorentzVector LeptonMinus;
389 G4double pdf = 0.; 431 G4double pdf = 0.;
390 432
391 G4double rndmv6[6] = {0.0}; << 433 G4double rndmv6[6];
392 const G4double corrFac = 1.0/(correctionInde <<
393 const G4double expLowLim = -20.; <<
394 const G4double logLowLim = G4Exp(expLowLim/c <<
395 G4double z0, z1, z2, x0, x1; <<
396 G4double betheheitler, sinTheta, cosTheta, d <<
397 // START Sampling 434 // START Sampling
398 do { 435 do {
399 436
400 rndmEngine->flatArray(6, rndmv6); 437 rndmEngine->flatArray(6, rndmv6);
401 438
402 ////////////////////////////////////////// 439 //////////////////////////////////////////////////
403 // pdf pow(x,c) with c = 1.4 440 // pdf pow(x,c) with c = 1.4
404 // integral y = pow(x,(c+1))/(c+1) @ x = 1 441 // integral y = pow(x,(c+1))/(c+1) @ x = 1 => y = 1 /(1+c)
405 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 442 // invCdf exp( log(y /* *( c + 1.0 )/ (c + 1.0 ) */ ) /( c + 1.0) )
406 ////////////////////////////////////////// 443 //////////////////////////////////////////////////
>> 444 const G4double X1 =
>> 445 G4Exp(G4Log(rndmv6[0])/(correctionIndex + 1.0));
407 446
408 z0 = (rndmv6[0] > logLowLim) ? G4Log(rndmv << 447 const G4double x0 = G4Exp(xl1 + (xu1 - xl1)*rndmv6[1]);
409 G4double X1 = (z0 > expLowLim) ? G4Exp(z0) << 448 const G4double dum0 = 1./(1.+x0);
410 z1 = xl1 + (xu1 - xl1)*rndmv6[1]; << 449 const G4double cosTheta = (x0-1.)*dum0;
411 if (z1 > expLowLim) { << 450 const G4double sinTheta = std::sqrt(4.*x0)*dum0;
412 x0 = G4Exp(z1); << 451
413 dum0 = 1.0/(1.0 + x0); << 452 const G4double PairInvMass = PairInvMassMin*G4Exp(X1*X1*lnPairInvMassRange);
414 x1 = dum0*x0; <<
415 cosTheta = -1.0 + 2.0*x1; <<
416 sinTheta = 2*std::sqrt(x1*(1.0 - x1)); <<
417 } else { <<
418 x0 = 0.0; <<
419 dum0 = 1.0; <<
420 cosTheta = -1.0; <<
421 sinTheta = 0.0; <<
422 } <<
423 453
424 z2 = X1*X1*lnPairInvMassRange; << 454 // G4double rndmv3[3];
425 const G4double PairInvMass = PairInvMassMi << 455 // rndmEngine->flatArray(3, rndmv3);
426 456
427 // cos and sin theta-lepton 457 // cos and sin theta-lepton
428 const G4double cosThetaLept = std::cos(pi* 458 const G4double cosThetaLept = std::cos(pi*rndmv6[2]);
429 // sin(ThetaLept) is always in [0,+1] if T 459 // sin(ThetaLept) is always in [0,+1] if ThetaLept is in [0,pi]
430 const G4double sinThetaLept = std::sqrt((1 460 const G4double sinThetaLept = std::sqrt((1.-cosThetaLept)*(1.+cosThetaLept));
431 // cos and sin phi-lepton 461 // cos and sin phi-lepton
432 const G4double cosPhiLept = std::cos(two 462 const G4double cosPhiLept = std::cos(twoPi*rndmv6[3]-pi);
433 // sin(PhiLept) is in [-1,0] if PhiLept in 463 // sin(PhiLept) is in [-1,0] if PhiLept in [-pi,0) and
434 // is in [0,+1] if PhiLept in 464 // is in [0,+1] if PhiLept in [0,+pi]
435 const G4double sinPhiLept = std::copysig 465 const G4double sinPhiLept = std::copysign(std::sqrt((1.-cosPhiLept)*(1.+cosPhiLept)),rndmv6[3]-0.5);
436 // cos and sin phi 466 // cos and sin phi
437 const G4double cosPhi = std::cos(two 467 const G4double cosPhi = std::cos(twoPi*rndmv6[4]-pi);
438 const G4double sinPhi = std::copysig << 468 const G4double sinPhi = std::copysign(std::sqrt((1.-cosPhi)*(1.+cosPhi)),rndmv6[4]-0.5);
439 469
440 ////////////////////////////////////////// 470 //////////////////////////////////////////////////
441 // frames: 471 // frames:
442 // 3 : the laboratory Lorentz frame, Geant 472 // 3 : the laboratory Lorentz frame, Geant4 axes definition
443 // 0 : the laboratory Lorentz frame, axes 473 // 0 : the laboratory Lorentz frame, axes along photon direction and polarisation
444 // 1 : the center-of-mass Lorentz frame 474 // 1 : the center-of-mass Lorentz frame
445 // 2 : the pair Lorentz frame 475 // 2 : the pair Lorentz frame
446 ////////////////////////////////////////// 476 //////////////////////////////////////////////////
447 477
448 // in the center-of-mass frame 478 // in the center-of-mass frame
449 479
450 const G4double RecEnergyCMS = (sCMSPlusRM 480 const G4double RecEnergyCMS = (sCMSPlusRM2-PairInvMass*PairInvMass)*isqrts2;
451 const G4double LeptonEnergy2 = PairInvMass 481 const G4double LeptonEnergy2 = PairInvMass*0.5;
452 482
453 // New way of calucaltion thePRecoil to av 483 // New way of calucaltion thePRecoil to avoid underflow
454 G4double abp = std::max((2.0*GammaEnergy*R 484 G4double abp = std::max((2.0*GammaEnergy*RecoilMass -
455 PairInvMass*PairInvMass + 2.0*PairI 485 PairInvMass*PairInvMass + 2.0*PairInvMass*RecoilMass)*
456 (2.0*GammaEnergy*R 486 (2.0*GammaEnergy*RecoilMass -
457 PairInvMass*PairInvMass - 2.0*PairI 487 PairInvMass*PairInvMass - 2.0*PairInvMass*RecoilMass),0.0);
458 488
459 G4double thePRecoil = std::sqrt(abp) * isq 489 G4double thePRecoil = std::sqrt(abp) * isqrts2;
460 490
461 // back to the center-of-mass frame 491 // back to the center-of-mass frame
462 Recoil.set( thePRecoil*sinTheta*cosPhi, 492 Recoil.set( thePRecoil*sinTheta*cosPhi,
463 thePRecoil*sinTheta*sinPhi, 493 thePRecoil*sinTheta*sinPhi,
464 thePRecoil*cosTheta, 494 thePRecoil*cosTheta,
465 RecEnergyCMS); 495 RecEnergyCMS);
466 496
467 // in the pair frame 497 // in the pair frame
468 const G4double thePLepton = std::sqrt( 498 const G4double thePLepton = std::sqrt( (LeptonEnergy2-LeptonMass)
469 * 499 *(LeptonEnergy2+LeptonMass));
470 500
471 LeptonPlus.set(thePLepton*sinThetaLept*cos 501 LeptonPlus.set(thePLepton*sinThetaLept*cosPhiLept,
472 thePLepton*sinThetaLept*sinPhiLept, 502 thePLepton*sinThetaLept*sinPhiLept,
473 thePLepton*cosThetaLept, 503 thePLepton*cosThetaLept,
474 LeptonEnergy2); 504 LeptonEnergy2);
475 505
476 LeptonMinus.set(-LeptonPlus.x(), 506 LeptonMinus.set(-LeptonPlus.x(),
477 -LeptonPlus.y(), 507 -LeptonPlus.y(),
478 -LeptonPlus.z(), 508 -LeptonPlus.z(),
479 LeptonEnergy2); 509 LeptonEnergy2);
480 510
481 511
482 // Normalisation of final state phase spac 512 // Normalisation of final state phase space:
483 // Section 47 of Particle Data Group, Chin 513 // Section 47 of Particle Data Group, Chin. Phys. C, 40, 100001 (2016)
484 // const G4double Norme = Recoil1.vect( 514 // const G4double Norme = Recoil1.vect().mag() * LeptonPlus2.vect().mag();
485 const G4double Norme = Recoil.vect().mag() 515 const G4double Norme = Recoil.vect().mag() * LeptonPlus.vect().mag();
486 516
487 // e+, e- to CMS frame from pair frame 517 // e+, e- to CMS frame from pair frame
488 518
489 // boost vector from Pair to CMS 519 // boost vector from Pair to CMS
490 const G4ThreeVector pair2cms = 520 const G4ThreeVector pair2cms =
491 G4LorentzVector( -Recoil.x(), -Recoil.y(), 521 G4LorentzVector( -Recoil.x(), -Recoil.y(), -Recoil.z(),
492 sqrts-RecEnergyCMS).boostVector(); 522 sqrts-RecEnergyCMS).boostVector();
493 523
494 LeptonPlus.boost(pair2cms); 524 LeptonPlus.boost(pair2cms);
495 LeptonMinus.boost(pair2cms); 525 LeptonMinus.boost(pair2cms);
496 526
497 // back to the laboratory frame (make use 527 // back to the laboratory frame (make use of the CMS(0,0,Eg,Eg+RM)) form
498 528
499 Recoil.boostZ(betaCMS); 529 Recoil.boostZ(betaCMS);
500 LeptonPlus.boostZ(betaCMS); 530 LeptonPlus.boostZ(betaCMS);
501 LeptonMinus.boostZ(betaCMS); 531 LeptonMinus.boostZ(betaCMS);
502 532
503 // Jacobian factors 533 // Jacobian factors
504 const G4double Jacob0 = x0*dum0*dum0; 534 const G4double Jacob0 = x0*dum0*dum0;
505 const G4double Jacob1 = 2.*X1*lnPairInvMas 535 const G4double Jacob1 = 2.*X1*lnPairInvMassRange*PairInvMass;
506 const G4double Jacob2 = std::abs(sinThetaL 536 const G4double Jacob2 = std::abs(sinThetaLept);
507 537
508 // there is no probability to have a lepto <<
509 // X and Y components of momentum may be z <<
510 const G4double EPlus = LeptonPlus.t(); 538 const G4double EPlus = LeptonPlus.t();
511 const G4double PPlus = LeptonPlus.vect().m 539 const G4double PPlus = LeptonPlus.vect().mag();
512 const G4double pPX = LeptonPlus.x(); << 540 const G4double sinThetaPlus = LeptonPlus.vect().perp()/PPlus;
513 const G4double pPY = LeptonPlus.y(); << 541 const G4double cosThetaPlus = LeptonPlus.vect().cosTheta();
514 const G4double pPZ = LeptonPlus.z(); << 542
515 G4double sinPhiPlus = 1.0; << 543 const G4double pPX = LeptonPlus.x();
516 G4double cosPhiPlus = 0.0; << 544 const G4double pPY = LeptonPlus.y();
517 G4double sinThetaPlus = 0.0; << 545 const G4double dum1 = 1./std::sqrt( pPX*pPX + pPY*pPY );
518 G4double cosThetaPlus = pPZ/PPlus; << 546 const G4double cosPhiPlus = pPX*dum1;
519 if (cosThetaPlus < 1.0 && cosThetaPlus > - << 547 const G4double sinPhiPlus = pPY*dum1;
520 sinThetaPlus = std::sqrt((1.0 - cosTheta <<
521 sinPhiPlus = pPY/(PPlus*sinThetaPlus); <<
522 cosPhiPlus = pPX/(PPlus*sinThetaPlus); <<
523 } <<
524 548
525 // denominators: 549 // denominators:
526 // the two cancelling leading terms for fo 550 // the two cancelling leading terms for forward emission at high energy, removed
527 const G4double elMassCTP = LeptonMass*cosT 551 const G4double elMassCTP = LeptonMass*cosThetaPlus;
528 const G4double ePlusSTP = EPlus*sinThetaP 552 const G4double ePlusSTP = EPlus*sinThetaPlus;
529 const G4double DPlus = (elMassCTP*elMa 553 const G4double DPlus = (elMassCTP*elMassCTP + ePlusSTP*ePlusSTP)
530 /(EPlus + PPlus* 554 /(EPlus + PPlus*cosThetaPlus);
531 555
532 // there is no probability to have a lepto <<
533 // X and Y components of momentum may be z <<
534 const G4double EMinus = LeptonMinus.t(); 556 const G4double EMinus = LeptonMinus.t();
535 const G4double PMinus = LeptonMinus.vect() 557 const G4double PMinus = LeptonMinus.vect().mag();
536 const G4double ePX = LeptonMinus.x(); << 558 const G4double sinThetaMinus = LeptonMinus.vect().perp()/PMinus;
537 const G4double ePY = LeptonMinus.y(); << 559 const G4double cosThetaMinus = LeptonMinus.vect().cosTheta();
538 const G4double ePZ = LeptonMinus.z(); << 560
539 G4double sinPhiMinus = 0.0; << 561 const G4double ePX = LeptonMinus.x();
540 G4double cosPhiMinus = 1.0; << 562 const G4double ePY = LeptonMinus.y();
541 G4double sinThetaMinus = 0.0; << 563 const G4double dum2 = 1./std::sqrt( ePX*ePX + ePY*ePY );
542 G4double cosThetaMinus = ePZ/PMinus; << 564 const G4double cosPhiMinus = ePX*dum2;
543 if (cosThetaMinus < 1.0 && cosThetaMinus > << 565 const G4double sinPhiMinus = ePY*dum2;
544 sinThetaMinus = std::sqrt((1.0 - cosThet <<
545 sinPhiMinus = ePY/(PMinus*sinThetaMinus) <<
546 cosPhiMinus = ePX/(PMinus*sinThetaMinus) <<
547 } <<
548 566
549 const G4double elMassCTM = LeptonMass*cosT 567 const G4double elMassCTM = LeptonMass*cosThetaMinus;
550 const G4double eMinSTM = EMinus*sinTheta 568 const G4double eMinSTM = EMinus*sinThetaMinus;
551 const G4double DMinus = (elMassCTM*elMa 569 const G4double DMinus = (elMassCTM*elMassCTM + eMinSTM*eMinSTM)
552 /(EMinus + PMinu 570 /(EMinus + PMinus*cosThetaMinus);
553 571
554 // cos(phiMinus-PhiPlus) 572 // cos(phiMinus-PhiPlus)
555 const G4double cosdPhi = cosPhiPlus*cosPhi 573 const G4double cosdPhi = cosPhiPlus*cosPhiMinus + sinPhiPlus*sinPhiMinus;
556 const G4double PRec = Recoil.vect().mag 574 const G4double PRec = Recoil.vect().mag();
557 const G4double q2 = PRec*PRec; 575 const G4double q2 = PRec*PRec;
558 576
559 const G4double BigPhi = -LeptonMass2 / (G 577 const G4double BigPhi = -LeptonMass2 / (GammaEnergy*GammaEnergy2 * q2*q2);
560 578
561 G4double FormFactor = 1.; 579 G4double FormFactor = 1.;
562 if (!iraw) { 580 if (!iraw) {
563 if (itriplet) { 581 if (itriplet) {
564 const G4double qun = factor1*iZ13*iZ13; 582 const G4double qun = factor1*iZ13*iZ13;
565 const G4double nun = qun * PRec; 583 const G4double nun = qun * PRec;
566 if (nun < 1.) { 584 if (nun < 1.) {
567 FormFactor = (nun < 0.01) ? (13.8-5 585 FormFactor = (nun < 0.01) ? (13.8-55.4*std::sqrt(nun))*nun
568 : std::sq 586 : std::sqrt(1-(nun-1)*(nun-1));
569 } // else FormFactor = 1 by default 587 } // else FormFactor = 1 by default
570 } else { 588 } else {
571 const G4double dum3 = 217.*PRec*iZ13; 589 const G4double dum3 = 217.*PRec*iZ13;
572 const G4double AFF = 1./(1. + dum3*dum3); 590 const G4double AFF = 1./(1. + dum3*dum3);
573 FormFactor = (1.-AFF)*(1-AFF); 591 FormFactor = (1.-AFF)*(1-AFF);
574 } 592 }
575 } // else FormFactor = 1 by default 593 } // else FormFactor = 1 by default
576 594
>> 595 G4double betheheitler;
577 if (GammaPolarizationMag==0.) { 596 if (GammaPolarizationMag==0.) {
578 const G4double pPlusSTP = PPlus*sinThe 597 const G4double pPlusSTP = PPlus*sinThetaPlus;
579 const G4double pMinusSTM = PMinus*sinTh 598 const G4double pMinusSTM = PMinus*sinThetaMinus;
580 const G4double pPlusSTPperDP = pPlusSTP 599 const G4double pPlusSTPperDP = pPlusSTP/DPlus;
581 const G4double pMinusSTMperDM = pMinusST 600 const G4double pMinusSTMperDM = pMinusSTM/DMinus;
582 const G4double dunpol = BigPhi*( 601 const G4double dunpol = BigPhi*(
583 pPlusSTPperDP *pPlusSTPperDP 602 pPlusSTPperDP *pPlusSTPperDP *(4.*EMinus*EMinus-q2)
584 + pMinusSTMperDM*pMinusSTMperD 603 + pMinusSTMperDM*pMinusSTMperDM*(4.*EPlus*EPlus - q2)
585 + 2.*pPlusSTPperDP*pMinusSTMpe 604 + 2.*pPlusSTPperDP*pMinusSTMperDM*cosdPhi
586 *(4.*EPlus*EMinus + q2 - 2 605 *(4.*EPlus*EMinus + q2 - 2.*GammaEnergy2)
587 - 2.*GammaEnergy2*(pPlusSTP*pP 606 - 2.*GammaEnergy2*(pPlusSTP*pPlusSTP+pMinusSTM*pMinusSTM)/(DMinus*DPlus));
588 betheheitler = dunpol * factor; 607 betheheitler = dunpol * factor;
589 } else { 608 } else {
590 const G4double pPlusSTP = PPlus*sinThet 609 const G4double pPlusSTP = PPlus*sinThetaPlus;
591 const G4double pMinusSTM = PMinus*sinThe 610 const G4double pMinusSTM = PMinus*sinThetaMinus;
592 const G4double pPlusSTPCPPperDP = pPlus 611 const G4double pPlusSTPCPPperDP = pPlusSTP*cosPhiPlus/DPlus;
593 const G4double pMinusSTMCPMperDM = pMinu 612 const G4double pMinusSTMCPMperDM = pMinusSTM*cosPhiMinus/DMinus;
594 const G4double caa = 2.*(EPlus*pMinusSTM 613 const G4double caa = 2.*(EPlus*pMinusSTMCPMperDM+EMinus*pPlusSTPCPPperDP);
595 const G4double cbb = pMinusSTMCPMperDM-p 614 const G4double cbb = pMinusSTMCPMperDM-pPlusSTPCPPperDP;
596 const G4double ccc = (pPlusSTP*pPlusSTP 615 const G4double ccc = (pPlusSTP*pPlusSTP + pMinusSTM*pMinusSTM
597 +2.*pPlusSTP*pMinusS 616 +2.*pPlusSTP*pMinusSTM*cosdPhi)/ (DMinus*DPlus);
598 const G4double dtot= 2.*BigPhi*( caa*caa 617 const G4double dtot= 2.*BigPhi*( caa*caa - q2*cbb*cbb - GammaEnergy2*ccc);
599 betheheitler = dtot * factor; 618 betheheitler = dtot * factor;
600 } 619 }
601 // 620 //
602 const G4double cross = Norme * Jacob0 * J 621 const G4double cross = Norme * Jacob0 * Jacob1 * Jacob2 * betheheitler
603 * FormFactor * Recoi 622 * FormFactor * RecoilMass / sqrts;
604 pdf = cross * (xu1 - xl1) / G4Exp(correcti 623 pdf = cross * (xu1 - xl1) / G4Exp(correctionIndex*G4Log(X1)); // cond1;
605 } while ( pdf < ymax * rndmv6[5] ); 624 } while ( pdf < ymax * rndmv6[5] );
606 // END of Sampling 625 // END of Sampling
607 << 626
608 if ( fVerbose > 2 ) { 627 if ( fVerbose > 2 ) {
609 G4double recul = std::sqrt(Recoil.x()*Reco 628 G4double recul = std::sqrt(Recoil.x()*Recoil.x()+Recoil.y()*Recoil.y()
610 +Recoil.z()*Reco 629 +Recoil.z()*Recoil.z());
611 G4cout << "BetheHeitler5DModel GammaEnergy 630 G4cout << "BetheHeitler5DModel GammaEnergy= " << GammaEnergy
612 << " PDF= " << pdf << " ymax= " << ymax 631 << " PDF= " << pdf << " ymax= " << ymax
613 << " recul= " << recul << G4endl; 632 << " recul= " << recul << G4endl;
614 } 633 }
615 634
616 // back to Geant4 system 635 // back to Geant4 system
617 636
618 if ( fVerbose > 4 ) { 637 if ( fVerbose > 4 ) {
619 G4cout << "BetheHeitler5DModel GammaDirect 638 G4cout << "BetheHeitler5DModel GammaDirection " << GammaDirection << G4endl;
620 G4cout << "BetheHeitler5DModel GammaPolari 639 G4cout << "BetheHeitler5DModel GammaPolarization " << GammaPolarization << G4endl;
621 G4cout << "BetheHeitler5DModel GammaEnergy 640 G4cout << "BetheHeitler5DModel GammaEnergy " << GammaEnergy << G4endl;
622 G4cout << "BetheHeitler5DModel Conv " 641 G4cout << "BetheHeitler5DModel Conv "
623 << (itriplet ? "triplet" : "nucl") << G4e 642 << (itriplet ? "triplet" : "nucl") << G4endl;
624 } 643 }
625 644
626 if (GammaPolarizationMag == 0.0) { 645 if (GammaPolarizationMag == 0.0) {
627 // set polarization axis orthohonal to dir 646 // set polarization axis orthohonal to direction
628 GammaPolarization = GammaDirection.orthogo 647 GammaPolarization = GammaDirection.orthogonal().unit();
629 } else { 648 } else {
630 // GammaPolarization not a unit vector 649 // GammaPolarization not a unit vector
631 GammaPolarization /= GammaPolarizationMag; 650 GammaPolarization /= GammaPolarizationMag;
632 } 651 }
633 652
634 // The unit norm vector that is orthogonal t 653 // The unit norm vector that is orthogonal to the two others
635 G4ThreeVector yGrec = GammaDirection.cross(G 654 G4ThreeVector yGrec = GammaDirection.cross(GammaPolarization);
636 655
637 // rotation from gamma ref. sys. to World 656 // rotation from gamma ref. sys. to World
638 G4RotationMatrix GtoW(GammaPolarization,yGre 657 G4RotationMatrix GtoW(GammaPolarization,yGrec,GammaDirection);
639 658
640 Recoil.transform(GtoW); 659 Recoil.transform(GtoW);
641 LeptonPlus.transform(GtoW); 660 LeptonPlus.transform(GtoW);
642 LeptonMinus.transform(GtoW); 661 LeptonMinus.transform(GtoW);
643 662
644 if ( fVerbose > 2 ) { 663 if ( fVerbose > 2 ) {
645 G4cout << "BetheHeitler5DModel Recoil " << 664 G4cout << "BetheHeitler5DModel Recoil " << Recoil.x() << " " << Recoil.y() << " " << Recoil.z()
646 << " " << Recoil.t() << " " << G4endl; 665 << " " << Recoil.t() << " " << G4endl;
647 G4cout << "BetheHeitler5DModel LeptonPlus 666 G4cout << "BetheHeitler5DModel LeptonPlus " << LeptonPlus.x() << " " << LeptonPlus.y() << " "
648 << LeptonPlus.z() << " " << LeptonPlus.t( 667 << LeptonPlus.z() << " " << LeptonPlus.t() << " " << G4endl;
649 G4cout << "BetheHeitler5DModel LeptonMinus 668 G4cout << "BetheHeitler5DModel LeptonMinus " << LeptonMinus.x() << " " << LeptonMinus.y() << " "
650 << LeptonMinus.z() << " " << LeptonMinus. 669 << LeptonMinus.z() << " " << LeptonMinus.t() << " " << G4endl;
651 } 670 }
652 671
653 // Create secondaries 672 // Create secondaries
654 auto aParticle1 = new G4DynamicParticle(fLep << 673 G4DynamicParticle* aParticle1 = new G4DynamicParticle(fLepton1,LeptonMinus);
655 auto aParticle2 = new G4DynamicParticle(fLep << 674 G4DynamicParticle* aParticle2 = new G4DynamicParticle(fLepton2,LeptonPlus);
656 675
657 // create G4DynamicParticle object for the p 676 // create G4DynamicParticle object for the particle3 ( recoil )
658 G4ParticleDefinition* RecoilPart; 677 G4ParticleDefinition* RecoilPart;
659 if (itriplet) { 678 if (itriplet) {
660 // triplet 679 // triplet
661 RecoilPart = fTheElectron; 680 RecoilPart = fTheElectron;
662 } else{ 681 } else{
663 RecoilPart = theIonTable->GetIon(Z, A, 0); 682 RecoilPart = theIonTable->GetIon(Z, A, 0);
664 } 683 }
665 auto aParticle3 = new G4DynamicParticle(Reco << 684 G4DynamicParticle* aParticle3 = new G4DynamicParticle(RecoilPart,Recoil);
666 685
667 // Fill output vector 686 // Fill output vector
668 fvect->push_back(aParticle1); 687 fvect->push_back(aParticle1);
669 fvect->push_back(aParticle2); 688 fvect->push_back(aParticle2);
670 fvect->push_back(aParticle3); 689 fvect->push_back(aParticle3);
671 690
672 // kill incident photon 691 // kill incident photon
673 fParticleChange->SetProposedKineticEnergy(0. 692 fParticleChange->SetProposedKineticEnergy(0.);
674 fParticleChange->ProposeTrackStatus(fStopAnd 693 fParticleChange->ProposeTrackStatus(fStopAndKill);
675 } 694 }
676 695
677 //....oooOO0OOooo........oooOO0OOooo........oo 696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
678 697