| 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 // Author: Sebastien Incerti 26 // Author: Sebastien Incerti
27 // 22 January 2012 27 // 22 January 2012
28 // on base of G4BoldyshevTripletModel 28 // on base of G4BoldyshevTripletModel (original version)
29 // and G4LivermoreRayleighModel (MT ve 29 // and G4LivermoreRayleighModel (MT version)
30 30
31 #include "G4BoldyshevTripletModel.hh" 31 #include "G4BoldyshevTripletModel.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4Log.hh" 34 #include "G4Log.hh"
35 #include "G4Exp.hh" 35 #include "G4Exp.hh"
36 36
37 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 38
39 using namespace std; 39 using namespace std;
40 40
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 42
43 const G4int G4BoldyshevTripletModel::maxZ; << 43 G4int G4BoldyshevTripletModel::maxZ = 99;
44 G4PhysicsFreeVector* G4BoldyshevTripletModel:: << 44 G4LPhysicsFreeVector* G4BoldyshevTripletModel::data[] = {0};
45 45
46 G4BoldyshevTripletModel::G4BoldyshevTripletMod 46 G4BoldyshevTripletModel::G4BoldyshevTripletModel(const G4ParticleDefinition*, const G4String& nam)
47 :G4VEmModel(nam),smallEnergy(4.*MeV) 47 :G4VEmModel(nam),smallEnergy(4.*MeV)
48 { 48 {
49 fParticleChange = nullptr; 49 fParticleChange = nullptr;
50 50
51 lowEnergyLimit = 4.0*electron_mass_c2; 51 lowEnergyLimit = 4.0*electron_mass_c2;
52 momentumThreshold_c = energyThreshold = xb = 52 momentumThreshold_c = energyThreshold = xb = xn = lowEnergyLimit;
53 53
54 verboseLevel= 0; 54 verboseLevel= 0;
55 // Verbosity scale for debugging purposes: 55 // Verbosity scale for debugging purposes:
56 // 0 = nothing 56 // 0 = nothing
57 // 1 = calculation of cross sections, file o 57 // 1 = calculation of cross sections, file openings...
58 // 2 = entering in methods 58 // 2 = entering in methods
59 59
60 if(verboseLevel > 0) 60 if(verboseLevel > 0)
61 { 61 {
62 G4cout << "G4BoldyshevTripletModel is cons 62 G4cout << "G4BoldyshevTripletModel is constructed " << G4endl;
63 } 63 }
64 } 64 }
65 65
66 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 67
68 G4BoldyshevTripletModel::~G4BoldyshevTripletMo 68 G4BoldyshevTripletModel::~G4BoldyshevTripletModel()
69 { 69 {
70 if(IsMaster()) { 70 if(IsMaster()) {
71 for(G4int i=0; i<maxZ; ++i) { 71 for(G4int i=0; i<maxZ; ++i) {
72 if(data[i]) { 72 if(data[i]) {
73 delete data[i]; 73 delete data[i];
74 data[i] = nullptr; 74 data[i] = nullptr;
75 } 75 }
76 } 76 }
77 } 77 }
78 } 78 }
79 79
80 //....oooOO0OOooo........oooOO0OOooo........oo 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 81
82 void G4BoldyshevTripletModel::Initialise(const 82 void G4BoldyshevTripletModel::Initialise(const G4ParticleDefinition*,
83 const G4DataVector&) 83 const G4DataVector&)
84 { 84 {
85 if (verboseLevel > 1) 85 if (verboseLevel > 1)
86 { 86 {
87 G4cout << "Calling Initialise() of G4Boldy 87 G4cout << "Calling Initialise() of G4BoldyshevTripletModel."
88 << G4endl 88 << G4endl
89 << "Energy range: " 89 << "Energy range: "
90 << LowEnergyLimit() / MeV << " MeV - " 90 << LowEnergyLimit() / MeV << " MeV - "
91 << HighEnergyLimit() / GeV << " GeV isMas 91 << HighEnergyLimit() / GeV << " GeV isMaster: " << IsMaster()
92 << G4endl; 92 << G4endl;
93 } 93 }
94 // compute values only once 94 // compute values only once
95 energyThreshold = 1.1*electron_mass_c2; 95 energyThreshold = 1.1*electron_mass_c2;
96 momentumThreshold_c = std::sqrt(energyThresh 96 momentumThreshold_c = std::sqrt(energyThreshold * energyThreshold
97 - electron_mas 97 - electron_mass_c2*electron_mass_c2);
98 G4double momentumThreshold_N = momentumThres 98 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2;
99 G4double t = 0.5*G4Log(momentumThreshold_N + 99 G4double t = 0.5*G4Log(momentumThreshold_N +
100 std::sqrt(momentumThreshold_N*momentumT 100 std::sqrt(momentumThreshold_N*momentumThreshold_N + 1.0));
101 //G4cout << 0.5*asinh(momentumThreshold_N) < 101 //G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
102 G4double sinht = std::sinh(t); 102 G4double sinht = std::sinh(t);
103 G4double cosht = std::cosh(t); 103 G4double cosht = std::cosh(t);
104 G4double logsinht = G4Log(2.*sinht); 104 G4double logsinht = G4Log(2.*sinht);
105 G4double J1 = 0.5*(t*cosht/sinht - logsinht) 105 G4double J1 = 0.5*(t*cosht/sinht - logsinht);
106 G4double J2 = (-2./3.)*logsinht + t*cosht/si 106 G4double J2 = (-2./3.)*logsinht + t*cosht/sinht
107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht* 107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht*sinht*sinht);
108 108
109 xb = 2.*(J1-J2)/J1; 109 xb = 2.*(J1-J2)/J1;
110 xn = 1. - xb/6.; 110 xn = 1. - xb/6.;
111 111
112 if(IsMaster()) 112 if(IsMaster())
113 { 113 {
114 // Access to elements 114 // Access to elements
115 const char* path = G4FindDataDir("G4LEDATA << 115 char* path = std::getenv("G4LEDATA");
116 116
117 G4ProductionCutsTable* theCoupleTable = 117 G4ProductionCutsTable* theCoupleTable =
118 G4ProductionCutsTable::GetProductionCuts 118 G4ProductionCutsTable::GetProductionCutsTable();
119 119
120 G4int numOfCouples = (G4int)theCoupleTable << 120 G4int numOfCouples = theCoupleTable->GetTableSize();
121 121
122 for(G4int i=0; i<numOfCouples; ++i) 122 for(G4int i=0; i<numOfCouples; ++i)
123 { 123 {
124 const G4Material* material = 124 const G4Material* material =
125 theCoupleTable->GetMaterialCutsCouple( 125 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
126 const G4ElementVector* theElementVector 126 const G4ElementVector* theElementVector = material->GetElementVector();
127 std::size_t nelm = material->GetNumberOf << 127 G4int nelm = material->GetNumberOfElements();
128 128
129 for (std::size_t j=0; j<nelm; ++j) << 129 for (G4int j=0; j<nelm; ++j)
130 { 130 {
131 G4int Z = std::min((*theElementVector) 131 G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
132 if(!data[Z]) { ReadData(Z, path); } 132 if(!data[Z]) { ReadData(Z, path); }
133 } 133 }
134 } 134 }
135 } 135 }
136 if(!fParticleChange) { 136 if(!fParticleChange) {
137 fParticleChange = GetParticleChangeForGamm 137 fParticleChange = GetParticleChangeForGamma();
138 } 138 }
139 } 139 }
140 140
141 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 142
143 G4double 143 G4double
144 G4BoldyshevTripletModel::MinPrimaryEnergy(cons 144 G4BoldyshevTripletModel::MinPrimaryEnergy(const G4Material*,
145 const G4ParticleDefinition*, 145 const G4ParticleDefinition*,
146 G4double) 146 G4double)
147 { 147 {
148 return lowEnergyLimit; 148 return lowEnergyLimit;
149 } 149 }
150 150
151 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 152
153 void G4BoldyshevTripletModel::ReadData(size_t 153 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path)
154 { 154 {
155 if (verboseLevel > 1) 155 if (verboseLevel > 1)
156 { 156 {
157 G4cout << "Calling ReadData() of G4Boldysh 157 G4cout << "Calling ReadData() of G4BoldyshevTripletModel"
158 << G4endl; 158 << G4endl;
159 } 159 }
160 160
161 if(data[Z]) { return; } 161 if(data[Z]) { return; }
162 162
163 const char* datadir = path; 163 const char* datadir = path;
164 164
165 if(!datadir) 165 if(!datadir)
166 { 166 {
167 datadir = G4FindDataDir("G4LEDATA"); << 167 datadir = std::getenv("G4LEDATA");
168 if(!datadir) 168 if(!datadir)
169 { 169 {
170 G4Exception("G4BoldyshevTripletModel::Re 170 G4Exception("G4BoldyshevTripletModel::ReadData()",
171 "em0006",FatalException, 171 "em0006",FatalException,
172 "Environment variable G4LEDATA not defin 172 "Environment variable G4LEDATA not defined");
173 return; 173 return;
174 } 174 }
175 } 175 }
176 176
177 data[Z] = new G4PhysicsFreeVector(0,/*spline << 177 data[Z] = new G4LPhysicsFreeVector();
178 std::ostringstream ost; 178 std::ostringstream ost;
179 ost << datadir << "/livermore/tripdata/pp-tr 179 ost << datadir << "/livermore/tripdata/pp-trip-cs-" << Z <<".dat";
180 std::ifstream fin(ost.str().c_str()); 180 std::ifstream fin(ost.str().c_str());
181 181
182 if( !fin.is_open()) 182 if( !fin.is_open())
183 { 183 {
184 G4ExceptionDescription ed; 184 G4ExceptionDescription ed;
185 ed << "G4BoldyshevTripletModel data file < 185 ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str()
186 << "> is not opened!" << G4endl; 186 << "> is not opened!" << G4endl;
187 G4Exception("G4BoldyshevTripletModel::Read 187 G4Exception("G4BoldyshevTripletModel::ReadData()",
188 "em0003",FatalException, 188 "em0003",FatalException,
189 ed,"G4LEDATA version should be G4EMLOW6.27 189 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
190 return; 190 return;
191 } 191 }
192 192
193 else 193 else
194 { 194 {
195 195
196 if(verboseLevel > 3) { G4cout << "File " < 196 if(verboseLevel > 3) { G4cout << "File " << ost.str()
197 << " is opened by G4BoldyshevTripletMod 197 << " is opened by G4BoldyshevTripletModel" << G4endl;}
198 198
199 data[Z]->Retrieve(fin, true); 199 data[Z]->Retrieve(fin, true);
200 } 200 }
201 201
202 // Activation of spline interpolation 202 // Activation of spline interpolation
203 data[Z]->FillSecondDerivatives(); << 203 data[Z]->SetSpline(true);
204 } 204 }
205 205
206 //....oooOO0OOooo........oooOO0OOooo........oo 206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 207
208 G4double G4BoldyshevTripletModel::ComputeCross 208 G4double G4BoldyshevTripletModel::ComputeCrossSectionPerAtom(
209 const G4ParticleDefinition* part, 209 const G4ParticleDefinition* part,
210 G4double GammaEnergy, G4double Z, G4double, 210 G4double GammaEnergy, G4double Z, G4double, G4double, G4double)
211 { 211 {
212 if (verboseLevel > 1) 212 if (verboseLevel > 1)
213 { 213 {
214 G4cout << "Calling ComputeCrossSectionPerA 214 G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
215 << G4endl; 215 << G4endl;
216 } 216 }
217 217
218 if (GammaEnergy < lowEnergyLimit) { return 0 218 if (GammaEnergy < lowEnergyLimit) { return 0.0; }
219 219
220 G4double xs = 0.0; 220 G4double xs = 0.0;
221 G4int intZ = std::max(1, std::min(G4lrint(Z) 221 G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ));
222 G4PhysicsFreeVector* pv = data[intZ]; << 222 G4LPhysicsFreeVector* pv = data[intZ];
223 223
224 // if element was not initialised 224 // if element was not initialised
225 // do initialisation safely for MT mode 225 // do initialisation safely for MT mode
226 if(!pv) 226 if(!pv)
227 { 227 {
228 InitialiseForElement(part, intZ); 228 InitialiseForElement(part, intZ);
229 pv = data[intZ]; 229 pv = data[intZ];
230 if(!pv) { return xs; } 230 if(!pv) { return xs; }
231 } 231 }
232 // x-section is taken from the table 232 // x-section is taken from the table
233 xs = pv->Value(GammaEnergy); 233 xs = pv->Value(GammaEnergy);
234 234
235 if(verboseLevel > 1) 235 if(verboseLevel > 1)
236 { 236 {
237 G4cout << "*** Triplet conversion xs for 237 G4cout << "*** Triplet conversion xs for Z=" << Z << " at energy E(MeV)="
238 << GammaEnergy/MeV << " cs=" << xs/mil 238 << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl;
239 } 239 }
240 return xs; 240 return xs;
241 } 241 }
242 242
243 //....oooOO0OOooo........oooOO0OOooo........oo 243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 244
245 void G4BoldyshevTripletModel::SampleSecondarie 245 void G4BoldyshevTripletModel::SampleSecondaries(
246 std::vector<G 246 std::vector<G4DynamicParticle*>* fvect,
247 const G4MaterialCutsCouple* /*couple* 247 const G4MaterialCutsCouple* /*couple*/,
248 const G4DynamicParticle* aDynamicGamm 248 const G4DynamicParticle* aDynamicGamma,
249 G4double, G4double) 249 G4double, G4double)
250 { 250 {
251 251
252 // The energies of the secondary particles a 252 // The energies of the secondary particles are sampled using
253 // a modified Wheeler-Lamb model (see PhysRe 253 // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
254 if (verboseLevel > 1) { 254 if (verboseLevel > 1) {
255 G4cout << "Calling SampleSecondaries() of 255 G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel"
256 << G4endl; 256 << G4endl;
257 } 257 }
258 258
259 G4double photonEnergy = aDynamicGamma->GetKi 259 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
260 G4ParticleMomentum photonDirection = aDynami 260 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
261 261
262 G4double epsilon; 262 G4double epsilon;
263 263
264 CLHEP::HepRandomEngine* rndmEngine = G4Rando 264 CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
265 265
266 // recoil electron thould be 3d particle 266 // recoil electron thould be 3d particle
267 G4DynamicParticle* particle3 = nullptr; 267 G4DynamicParticle* particle3 = nullptr;
268 static const G4double costlim = std::cos(4.4 268 static const G4double costlim = std::cos(4.47*CLHEP::pi/180.);
269 269
270 G4double loga, f1_re, greject, cost; 270 G4double loga, f1_re, greject, cost;
271 G4double cosThetaMax = (energyThreshold - el 271 G4double cosThetaMax = (energyThreshold - electron_mass_c2
272 + electron_mass_c2*(energyThreshold + ele 272 + electron_mass_c2*(energyThreshold + electron_mass_c2)/photonEnergy )
273 /momentumThreshold_c; 273 /momentumThreshold_c;
274 if (cosThetaMax > 1.) { 274 if (cosThetaMax > 1.) {
275 //G4cout << "G4BoldyshevTripletModel::Samp 275 //G4cout << "G4BoldyshevTripletModel::SampleSecondaries: ERROR cosThetaMax= "
276 // << cosThetaMax << G4endl; 276 // << cosThetaMax << G4endl;
277 cosThetaMax = 1.0; 277 cosThetaMax = 1.0;
278 } 278 }
279 279
280 G4double logcostm = G4Log(cosThetaMax); 280 G4double logcostm = G4Log(cosThetaMax);
>> 281 G4int nn = 0;
281 do { 282 do {
282 cost = G4Exp(logcostm*rndmEngine->flat()); 283 cost = G4Exp(logcostm*rndmEngine->flat());
283 G4double are = 1./(14.*cost*cost); 284 G4double are = 1./(14.*cost*cost);
284 G4double bre = (1.-5.*cost*cost)/(2.*cost) 285 G4double bre = (1.-5.*cost*cost)/(2.*cost);
285 loga = G4Log((1.+ cost)/(1.- cost)); 286 loga = G4Log((1.+ cost)/(1.- cost));
286 f1_re = 1. - bre*loga; 287 f1_re = 1. - bre*loga;
287 greject = (cost < costlim) ? are*f1_re : 1 288 greject = (cost < costlim) ? are*f1_re : 1.0;
>> 289 // G4cout << nn << ". step of the 1st loop greject= " << greject << G4endl;
>> 290 ++nn;
288 } while(greject < rndmEngine->flat()); 291 } while(greject < rndmEngine->flat());
289 292
290 // Calculo de phi - elecron de recoil 293 // Calculo de phi - elecron de recoil
291 G4double sint2 = (1. - cost)*(1. + cost); 294 G4double sint2 = (1. - cost)*(1. + cost);
292 G4double fp = 1. - sint2*loga/(2.*cost) ; 295 G4double fp = 1. - sint2*loga/(2.*cost) ;
293 G4double rt, phi_re; 296 G4double rt, phi_re;
>> 297 nn = 0;
294 do { 298 do {
295 phi_re = twopi*rndmEngine->flat(); 299 phi_re = twopi*rndmEngine->flat();
296 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/t 300 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
>> 301 //G4cout << nn << ". step of the 2nd loop greject= " << rt << G4endl;
>> 302 ++nn;
297 } while(rt < rndmEngine->flat()); 303 } while(rt < rndmEngine->flat());
298 304
299 // Calculo de la energia - elecron de recoil 305 // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
300 G4double S = electron_mass_c2*(2.* photonEn 306 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
301 G4double P2 = S - electron_mass_c2*electron_ 307 G4double P2 = S - electron_mass_c2*electron_mass_c2;
302 308
303 G4double D2 = 4.*S * electron_mass_c2*electr 309 G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
304 G4double ener_re = electron_mass_c2 * (S + e 310 G4double ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
305 311
306 if(ener_re >= energyThreshold) 312 if(ener_re >= energyThreshold)
307 { 313 {
308 G4double electronRKineEnergy = ener_re - 314 G4double electronRKineEnergy = ener_re - electron_mass_c2;
309 G4double sint = std::sqrt(sint2); 315 G4double sint = std::sqrt(sint2);
310 G4ThreeVector electronRDirection (sint*s 316 G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
311 electronRDirection.rotateUz(photonDirect 317 electronRDirection.rotateUz(photonDirection);
312 particle3 = new G4DynamicParticle (G4Ele 318 particle3 = new G4DynamicParticle (G4Electron::Electron(),
313 electronRDirection, 319 electronRDirection,
314 electronRKineEnergy); 320 electronRKineEnergy);
315 } 321 }
316 else 322 else
317 { 323 {
318 // deposito la energia ener_re - electr 324 // deposito la energia ener_re - electron_mass_c2
319 // G4cout << "electron de retroceso " << 325 // G4cout << "electron de retroceso " << ener_re << G4endl;
320 fParticleChange->ProposeLocalEnergyDepos 326 fParticleChange->ProposeLocalEnergyDeposit(std::max(0.0, ener_re - electron_mass_c2));
321 ener_re = 0.0; 327 ener_re = 0.0;
322 } 328 }
323 329
324 // Depaola (2004) suggested distribution for 330 // Depaola (2004) suggested distribution for e+e- energy
325 // VI: very suspect that 1 random number is 331 // VI: very suspect that 1 random number is not enough
326 // and sampling below is not correct - s 332 // and sampling below is not correct - should be fixed
327 G4double re = rndmEngine->flat(); 333 G4double re = rndmEngine->flat();
328 334
329 G4double a = std::sqrt(16./xb - 3. - 36.*re 335 G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
330 G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + 336 G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + xb + 2*a)*xb*xb)/3.);
331 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0 337 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
332 338
333 G4double photonEnergy1 = photonEnergy - ener 339 G4double photonEnergy1 = photonEnergy - ener_re ;
334 // resto al foton la energia del electron de 340 // resto al foton la energia del electron de retro.
335 G4double positronTotEnergy = std::max(epsilo 341 G4double positronTotEnergy = std::max(epsilon*photonEnergy1, electron_mass_c2);
336 G4double electronTotEnergy = std::max(photon 342 G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
337 343
338 static const G4double a1 = 1.6; 344 static const G4double a1 = 1.6;
339 static const G4double a2 = 0.5333333333; 345 static const G4double a2 = 0.5333333333;
340 G4double uu = -G4Log(rndmEngine->flat()*rndm 346 G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat());
341 G4double u = (0.25 > rndmEngine->flat()) ? u 347 G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2;
342 348
343 G4double thetaEle = u*electron_mass_c2/elect 349 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
344 G4double sinte = std::sin(thetaEle); 350 G4double sinte = std::sin(thetaEle);
345 G4double coste = std::cos(thetaEle); 351 G4double coste = std::cos(thetaEle);
346 352
347 G4double thetaPos = u*electron_mass_c2/posit 353 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
348 G4double sintp = std::sin(thetaPos); 354 G4double sintp = std::sin(thetaPos);
349 G4double costp = std::cos(thetaPos); 355 G4double costp = std::cos(thetaPos);
350 356
351 G4double phi = twopi * rndmEngine->flat(); 357 G4double phi = twopi * rndmEngine->flat();
352 G4double sinp = std::sin(phi); 358 G4double sinp = std::sin(phi);
353 G4double cosp = std::cos(phi); 359 G4double cosp = std::cos(phi);
354 360
355 // Kinematics of the created pair: 361 // Kinematics of the created pair:
356 // the electron and positron are assumed to 362 // the electron and positron are assumed to have a symetric angular
357 // distribution with respect to the Z axis a 363 // distribution with respect to the Z axis along the parent photon
358 364
359 G4double electronKineEnergy = electronTotEne 365 G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
360 366
361 G4ThreeVector electronDirection (sinte*cosp, 367 G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
362 electronDirection.rotateUz(photonDirection); 368 electronDirection.rotateUz(photonDirection);
363 369
364 G4DynamicParticle* particle1 = new G4Dynamic 370 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
365 371 electronDirection,
366 electronKineEnergy); 372 electronKineEnergy);
367 373
368 G4double positronKineEnergy = positronTotEne 374 G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
369 375
370 G4ThreeVector positronDirection (-sintp*cosp 376 G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
371 positronDirection.rotateUz(photonDirection); 377 positronDirection.rotateUz(photonDirection);
372 378
373 // Create G4DynamicParticle object for the p 379 // Create G4DynamicParticle object for the particle2
374 G4DynamicParticle* particle2 = new G4Dynamic 380 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
375 381 positronDirection, positronKineEnergy);
376 // Fill output vector 382 // Fill output vector
377 383
378 fvect->push_back(particle1); 384 fvect->push_back(particle1);
379 fvect->push_back(particle2); 385 fvect->push_back(particle2);
380 386
381 if(particle3) { fvect->push_back(particle3); 387 if(particle3) { fvect->push_back(particle3); }
382 388
383 // kill incident photon 389 // kill incident photon
384 fParticleChange->SetProposedKineticEnergy(0. 390 fParticleChange->SetProposedKineticEnergy(0.);
385 fParticleChange->ProposeTrackStatus(fStopAnd 391 fParticleChange->ProposeTrackStatus(fStopAndKill);
386 } 392 }
387 393
388 //....oooOO0OOooo........oooOO0OOooo........oo 394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
389 395
390 #include "G4AutoLock.hh" 396 #include "G4AutoLock.hh"
391 namespace { G4Mutex BoldyshevTripletModelMutex 397 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; }
392 398
393 void G4BoldyshevTripletModel::InitialiseForEle 399 void G4BoldyshevTripletModel::InitialiseForElement(
394 const G4ParticleDefinition*, G4int Z) 400 const G4ParticleDefinition*, G4int Z)
395 { 401 {
396 G4AutoLock l(&BoldyshevTripletModelMutex); 402 G4AutoLock l(&BoldyshevTripletModelMutex);
397 // G4cout << "G4BoldyshevTripletModel::Initi 403 // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= "
398 // << Z << G4endl; 404 // << Z << G4endl;
399 if(!data[Z]) { ReadData(Z); } 405 if(!data[Z]) { ReadData(Z); }
400 l.unlock(); 406 l.unlock();
401 } 407 }
402 408
403 //....oooOO0OOooo........oooOO0OOooo........oo 409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 410