| 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),isInitialised(false),smallEnergy(4.*MeV)
48 { 48 {
49 fParticleChange = nullptr; << 49 fParticleChange = 0;
50 50
51 lowEnergyLimit = 4.0*electron_mass_c2; 51 lowEnergyLimit = 4.0*electron_mass_c2;
52 momentumThreshold_c = energyThreshold = xb = <<
53 52
54 verboseLevel= 0; 53 verboseLevel= 0;
55 // Verbosity scale for debugging purposes: 54 // Verbosity scale for debugging purposes:
56 // 0 = nothing 55 // 0 = nothing
57 // 1 = calculation of cross sections, file o 56 // 1 = calculation of cross sections, file openings...
58 // 2 = entering in methods 57 // 2 = entering in methods
59 58
60 if(verboseLevel > 0) 59 if(verboseLevel > 0)
61 { 60 {
62 G4cout << "G4BoldyshevTripletModel is cons 61 G4cout << "G4BoldyshevTripletModel is constructed " << G4endl;
63 } 62 }
64 } 63 }
65 64
66 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
67 66
68 G4BoldyshevTripletModel::~G4BoldyshevTripletMo 67 G4BoldyshevTripletModel::~G4BoldyshevTripletModel()
69 { 68 {
70 if(IsMaster()) { 69 if(IsMaster()) {
71 for(G4int i=0; i<maxZ; ++i) { 70 for(G4int i=0; i<maxZ; ++i) {
72 if(data[i]) { 71 if(data[i]) {
73 delete data[i]; 72 delete data[i];
74 data[i] = nullptr; << 73 data[i] = 0;
75 } 74 }
76 } 75 }
77 } 76 }
78 } 77 }
79 78
80 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 80
82 void G4BoldyshevTripletModel::Initialise(const << 81 void G4BoldyshevTripletModel::Initialise(
83 const G4DataVector&) << 82 const G4ParticleDefinition* particle,
>> 83 const G4DataVector& cuts)
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"
92 << G4endl; 92 << G4endl;
93 } 93 }
94 // compute values only once <<
95 energyThreshold = 1.1*electron_mass_c2; <<
96 momentumThreshold_c = std::sqrt(energyThresh <<
97 - electron_mas <<
98 G4double momentumThreshold_N = momentumThres <<
99 G4double t = 0.5*G4Log(momentumThreshold_N + <<
100 std::sqrt(momentumThreshold_N*momentumT <<
101 //G4cout << 0.5*asinh(momentumThreshold_N) < <<
102 G4double sinht = std::sinh(t); <<
103 G4double cosht = std::cosh(t); <<
104 G4double logsinht = G4Log(2.*sinht); <<
105 G4double J1 = 0.5*(t*cosht/sinht - logsinht) <<
106 G4double J2 = (-2./3.)*logsinht + t*cosht/si <<
107 + (sinht - t*cosht*cosht*cosht)/(3.*sinht* <<
108 <<
109 xb = 2.*(J1-J2)/J1; <<
110 xn = 1. - xb/6.; <<
111 94
112 if(IsMaster()) 95 if(IsMaster())
113 { 96 {
114 // Access to elements << 97
115 const char* path = G4FindDataDir("G4LEDATA << 98 // Initialise element selector
>> 99
>> 100 InitialiseElementSelectors(particle, cuts);
>> 101
>> 102 // Access to elements
>> 103
>> 104 char* path = getenv("G4LEDATA");
116 105
117 G4ProductionCutsTable* theCoupleTable = 106 G4ProductionCutsTable* theCoupleTable =
118 G4ProductionCutsTable::GetProductionCuts 107 G4ProductionCutsTable::GetProductionCutsTable();
119 108
120 G4int numOfCouples = (G4int)theCoupleTable << 109 G4int numOfCouples = theCoupleTable->GetTableSize();
121 110
122 for(G4int i=0; i<numOfCouples; ++i) 111 for(G4int i=0; i<numOfCouples; ++i)
123 { 112 {
124 const G4Material* material = 113 const G4Material* material =
125 theCoupleTable->GetMaterialCutsCouple( 114 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
126 const G4ElementVector* theElementVector 115 const G4ElementVector* theElementVector = material->GetElementVector();
127 std::size_t nelm = material->GetNumberOf << 116 G4int nelm = material->GetNumberOfElements();
128 117
129 for (std::size_t j=0; j<nelm; ++j) << 118 for (G4int j=0; j<nelm; ++j)
130 { 119 {
131 G4int Z = std::min((*theElementVector) << 120 G4int Z = (G4int)(*theElementVector)[j]->GetZ();
>> 121 if(Z < 1) { Z = 1; }
>> 122 else if(Z > maxZ) { Z = maxZ; }
132 if(!data[Z]) { ReadData(Z, path); } 123 if(!data[Z]) { ReadData(Z, path); }
133 } 124 }
134 } 125 }
135 } 126 }
136 if(!fParticleChange) { << 127 if(isInitialised) { return; }
137 fParticleChange = GetParticleChangeForGamm << 128 fParticleChange = GetParticleChangeForGamma();
138 } << 129 isInitialised = true;
>> 130 }
>> 131
>> 132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 133
>> 134 void G4BoldyshevTripletModel::InitialiseLocal(
>> 135 const G4ParticleDefinition*, G4VEmModel* masterModel)
>> 136 {
>> 137 SetElementSelectors(masterModel->GetElementSelectors());
139 } 138 }
140 139
141 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 141
143 G4double 142 G4double
144 G4BoldyshevTripletModel::MinPrimaryEnergy(cons 143 G4BoldyshevTripletModel::MinPrimaryEnergy(const G4Material*,
145 const G4ParticleDefinition*, << 144 const G4ParticleDefinition*,
146 G4double) << 145 G4double)
147 { 146 {
148 return lowEnergyLimit; 147 return lowEnergyLimit;
149 } 148 }
150 149
151 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 151
153 void G4BoldyshevTripletModel::ReadData(size_t 152 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path)
154 { 153 {
155 if (verboseLevel > 1) 154 if (verboseLevel > 1)
156 { 155 {
157 G4cout << "Calling ReadData() of G4Boldysh 156 G4cout << "Calling ReadData() of G4BoldyshevTripletModel"
158 << G4endl; 157 << G4endl;
159 } 158 }
160 159
161 if(data[Z]) { return; } 160 if(data[Z]) { return; }
162 161
163 const char* datadir = path; 162 const char* datadir = path;
164 163
165 if(!datadir) 164 if(!datadir)
166 { 165 {
167 datadir = G4FindDataDir("G4LEDATA"); << 166 datadir = getenv("G4LEDATA");
168 if(!datadir) 167 if(!datadir)
169 { 168 {
170 G4Exception("G4BoldyshevTripletModel::Re 169 G4Exception("G4BoldyshevTripletModel::ReadData()",
171 "em0006",FatalException, 170 "em0006",FatalException,
172 "Environment variable G4LEDATA not defin 171 "Environment variable G4LEDATA not defined");
173 return; 172 return;
174 } 173 }
175 } 174 }
>> 175
>> 176 //
>> 177
>> 178 data[Z] = new G4LPhysicsFreeVector();
>> 179
>> 180 //
176 181
177 data[Z] = new G4PhysicsFreeVector(0,/*spline <<
178 std::ostringstream ost; 182 std::ostringstream ost;
179 ost << datadir << "/livermore/tripdata/pp-tr << 183 ost << datadir << "livermore/tripdata/pp-trip-cs-" << Z <<".dat";
180 std::ifstream fin(ost.str().c_str()); 184 std::ifstream fin(ost.str().c_str());
181 185
182 if( !fin.is_open()) 186 if( !fin.is_open())
183 { 187 {
184 G4ExceptionDescription ed; 188 G4ExceptionDescription ed;
185 ed << "G4BoldyshevTripletModel data file < 189 ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str()
186 << "> is not opened!" << G4endl; 190 << "> is not opened!" << G4endl;
187 G4Exception("G4BoldyshevTripletModel::Read 191 G4Exception("G4BoldyshevTripletModel::ReadData()",
188 "em0003",FatalException, 192 "em0003",FatalException,
189 ed,"G4LEDATA version should be G4EMLOW6.27 193 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
190 return; 194 return;
191 } 195 }
192 196
193 else 197 else
194 { 198 {
195 199
196 if(verboseLevel > 3) { G4cout << "File " < 200 if(verboseLevel > 3) { G4cout << "File " << ost.str()
197 << " is opened by G4BoldyshevTripletMod 201 << " is opened by G4BoldyshevTripletModel" << G4endl;}
198 202
199 data[Z]->Retrieve(fin, true); 203 data[Z]->Retrieve(fin, true);
200 } 204 }
201 205
202 // Activation of spline interpolation 206 // Activation of spline interpolation
203 data[Z]->FillSecondDerivatives(); << 207 data[Z] ->SetSpline(true);
>> 208
204 } 209 }
205 210
206 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 212
208 G4double G4BoldyshevTripletModel::ComputeCross << 213 G4double
209 const G4ParticleDefinition* part, << 214 G4BoldyshevTripletModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
210 G4double GammaEnergy, G4double Z, G4double, << 215 G4double GammaEnergy,
>> 216 G4double Z, G4double,
>> 217 G4double, G4double)
211 { 218 {
212 if (verboseLevel > 1) 219 if (verboseLevel > 1)
213 { 220 {
214 G4cout << "Calling ComputeCrossSectionPerA 221 G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
215 << G4endl; 222 << G4endl;
216 } 223 }
217 224
218 if (GammaEnergy < lowEnergyLimit) { return 0 225 if (GammaEnergy < lowEnergyLimit) { return 0.0; }
219 226
220 G4double xs = 0.0; << 227 G4double xs = 0.0;
221 G4int intZ = std::max(1, std::min(G4lrint(Z) << 228
222 G4PhysicsFreeVector* pv = data[intZ]; << 229 G4int intZ=G4int(Z);
>> 230
>> 231 if(intZ < 1 || intZ > maxZ) { return xs; }
>> 232
>> 233 G4LPhysicsFreeVector* pv = data[intZ];
223 234
224 // if element was not initialised 235 // if element was not initialised
225 // do initialisation safely for MT mode 236 // do initialisation safely for MT mode
226 if(!pv) 237 if(!pv)
227 { 238 {
228 InitialiseForElement(part, intZ); << 239 InitialiseForElement(0, intZ);
229 pv = data[intZ]; 240 pv = data[intZ];
230 if(!pv) { return xs; } 241 if(!pv) { return xs; }
231 } 242 }
232 // x-section is taken from the table 243 // x-section is taken from the table
233 xs = pv->Value(GammaEnergy); 244 xs = pv->Value(GammaEnergy);
234 245
235 if(verboseLevel > 1) << 246 if(verboseLevel > 0)
236 { 247 {
237 G4cout << "*** Triplet conversion xs for << 248 G4int n = pv->GetVectorLength() - 1;
238 << GammaEnergy/MeV << " cs=" << xs/mil << 249 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
>> 250 << GammaEnergy/MeV << G4endl;
>> 251 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
>> 252 G4cout << " -> first cs value in EADL data file (iu) =" << (*pv)[0] << G4endl;
>> 253 G4cout << " -> last cs value in EADL data file (iu) =" << (*pv)[n] << G4endl;
>> 254 G4cout << "*********************************************************" << G4endl;
239 } 255 }
>> 256
240 return xs; 257 return xs;
>> 258
241 } 259 }
242 260
243 //....oooOO0OOooo........oooOO0OOooo........oo 261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 262
245 void G4BoldyshevTripletModel::SampleSecondarie 263 void G4BoldyshevTripletModel::SampleSecondaries(
246 std::vector<G 264 std::vector<G4DynamicParticle*>* fvect,
247 const G4MaterialCutsCouple* /*couple* 265 const G4MaterialCutsCouple* /*couple*/,
248 const G4DynamicParticle* aDynamicGamm 266 const G4DynamicParticle* aDynamicGamma,
249 G4double, G4double) 267 G4double, G4double)
250 { 268 {
251 269
252 // The energies of the secondary particles a << 270 // The energies of the secondary particles are sampled using // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
253 // a modified Wheeler-Lamb model (see PhysRe << 271
254 if (verboseLevel > 1) { 272 if (verboseLevel > 1) {
255 G4cout << "Calling SampleSecondaries() of 273 G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel"
256 << G4endl; 274 << G4endl;
257 } 275 }
258 276
259 G4double photonEnergy = aDynamicGamma->GetKi 277 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
260 G4ParticleMomentum photonDirection = aDynami 278 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
261 279
262 G4double epsilon; << 280 G4double epsilon ;
>> 281 // G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
263 282
264 CLHEP::HepRandomEngine* rndmEngine = G4Rando << 283
>> 284 G4double p0 = electron_mass_c2;
>> 285 G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
>> 286 G4double ener_re=0., theta_re, phi_re, phi;
>> 287
>> 288 // Calculo de theta - elecron de recoil
>> 289
>> 290 G4double energyThreshold = sqrt(2.)*electron_mass_c2; // -> momentumThreshold_N = 1
>> 291
>> 292 energyThreshold = 1.1*electron_mass_c2;
>> 293 // G4cout << energyThreshold << G4endl;
>> 294
>> 295
>> 296 G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2); // momentun in MeV/c unit
>> 297 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; // momentun in mc unit
265 298
266 // recoil electron thould be 3d particle << 299 // Calculation of recoil electron production
267 G4DynamicParticle* particle3 = nullptr; << 300
268 static const G4double costlim = std::cos(4.4 << 301
269 << 302 G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ;
270 G4double loga, f1_re, greject, cost; << 303 G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
271 G4double cosThetaMax = (energyThreshold - el << 304 G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0);
272 + electron_mass_c2*(energyThreshold + ele << 305 G4double recoilProb = G4UniformRand();
273 /momentumThreshold_c; << 306 //G4cout << "SIGMA TOT " << SigmaTot << " " << "SigmaQ " << SigmaQ << " " << SigmaQ/SigmaTot << " " << recoilProb << G4endl;
274 if (cosThetaMax > 1.) { << 307
275 //G4cout << "G4BoldyshevTripletModel::Samp << 308
276 // << cosThetaMax << G4endl; << 309 if (recoilProb >= SigmaQ/SigmaTot) // create electron recoil
277 cosThetaMax = 1.0; << 310 {
278 } << 311
>> 312 G4double cosThetaMax = ( ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
>> 313 ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
>> 314
>> 315
>> 316 if (cosThetaMax > 1) G4cout << "ERRORE " << G4endl;
279 317
280 G4double logcostm = G4Log(cosThetaMax); << 318 G4double r1;
281 do { << 319 G4double r2;
282 cost = G4Exp(logcostm*rndmEngine->flat()); << 320 G4double are, bre, loga, f1_re, greject, cost;
283 G4double are = 1./(14.*cost*cost); << 321
284 G4double bre = (1.-5.*cost*cost)/(2.*cost) << 322 do {
285 loga = G4Log((1.+ cost)/(1.- cost)); << 323 r1 = G4UniformRand();
286 f1_re = 1. - bre*loga; << 324 r2 = G4UniformRand();
287 greject = (cost < costlim) ? are*f1_re : 1 << 325 // cost = (pow(4./enern,0.5*r1)) ;
288 } while(greject < rndmEngine->flat()); << 326
>> 327 cost = pow(cosThetaMax,r1);
>> 328 theta_re = acos(cost);
>> 329 are = 1./(14.*cost*cost);
>> 330 bre = (1.-5.*cost*cost)/(2.*cost);
>> 331 loga = log((1.+ cost)/(1.- cost));
>> 332 f1_re = 1. - bre*loga;
>> 333
>> 334 if ( theta_re >= 4.47*CLHEP::pi/180.)
>> 335 {
>> 336 greject = are*f1_re;
>> 337 } else {
>> 338 greject = 1. ;
>> 339 }
>> 340 } while(greject < r2);
289 341
290 // Calculo de phi - elecron de recoil << 342 // Calculo de phi - elecron de recoil
291 G4double sint2 = (1. - cost)*(1. + cost); << 343
292 G4double fp = 1. - sint2*loga/(2.*cost) ; << 344 G4double r3, r4, rt;
293 G4double rt, phi_re; <<
294 do { <<
295 phi_re = twopi*rndmEngine->flat(); <<
296 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/t <<
297 } while(rt < rndmEngine->flat()); <<
298 <<
299 // Calculo de la energia - elecron de recoil <<
300 G4double S = electron_mass_c2*(2.* photonEn <<
301 G4double P2 = S - electron_mass_c2*electron_ <<
302 345
303 G4double D2 = 4.*S * electron_mass_c2*electr << 346 do {
304 G4double ener_re = electron_mass_c2 * (S + e << 347 r3 = G4UniformRand();
>> 348 r4 = G4UniformRand();
>> 349 phi_re = twopi*r3 ;
>> 350 G4double sint2 = 1. - cost*cost ;
>> 351 G4double fp = 1. - sint2*loga/(2.*cost) ;
>> 352 rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
>> 353
>> 354 } while(rt < r4);
>> 355
>> 356 // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
>> 357
>> 358 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
>> 359
>> 360 G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
>> 361 + (S - electron_mass_c2*electron_mass_c2)
>> 362 *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
>> 363 ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
305 364
306 if(ener_re >= energyThreshold) << 365 // New Recoil energy calculation
307 { << 366
308 G4double electronRKineEnergy = ener_re - << 367 G4double momentum_recoil = 2* (electron_mass_c2) * (std::cos(theta_re)/(std::sin(phi_re)*std::sin(phi_re)));
309 G4double sint = std::sqrt(sint2); << 368 G4double ener_recoil = sqrt( momentum_recoil*momentum_recoil + electron_mass_c2*electron_mass_c2);
310 G4ThreeVector electronRDirection (sint*s << 369 ener_re = ener_recoil;
>> 370
>> 371 // G4cout << "electron de retroceso " << ener_re << " " << theta_re << " " << phi_re << G4endl;
>> 372
>> 373 // Recoil electron creation
>> 374 G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
>> 375
>> 376 G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
>> 377
>> 378 G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
311 electronRDirection.rotateUz(photonDirect 379 electronRDirection.rotateUz(photonDirection);
312 particle3 = new G4DynamicParticle (G4Ele << 380
313 electronRDirection, << 381 G4DynamicParticle* particle3 = new G4DynamicParticle (G4Electron::Electron(),
314 electronRKineEnergy); << 382 electronRDirection,
>> 383 electronRKineEnergy);
>> 384 fvect->push_back(particle3);
>> 385
315 } 386 }
316 else 387 else
317 { 388 {
318 // deposito la energia ener_re - electr << 389 // deposito la energia ener_re - electron_mass_c2
319 // G4cout << "electron de retroceso " << << 390 // G4cout << "electron de retroceso " << ener_re << G4endl;
320 fParticleChange->ProposeLocalEnergyDepos << 391
321 ener_re = 0.0; << 392 fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
322 } 393 }
323 394
324 // Depaola (2004) suggested distribution for << 395
325 // VI: very suspect that 1 random number is << 396 // Depaola (2004) suggested distribution for e+e- energy
326 // and sampling below is not correct - s <<
327 G4double re = rndmEngine->flat(); <<
328 <<
329 G4double a = std::sqrt(16./xb - 3. - 36.*re <<
330 G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + <<
331 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0 <<
332 <<
333 G4double photonEnergy1 = photonEnergy - ener <<
334 // resto al foton la energia del electron de <<
335 G4double positronTotEnergy = std::max(epsilo <<
336 G4double electronTotEnergy = std::max(photon <<
337 <<
338 static const G4double a1 = 1.6; <<
339 static const G4double a2 = 0.5333333333; <<
340 G4double uu = -G4Log(rndmEngine->flat()*rndm <<
341 G4double u = (0.25 > rndmEngine->flat()) ? u <<
342 <<
343 G4double thetaEle = u*electron_mass_c2/elect <<
344 G4double sinte = std::sin(thetaEle); <<
345 G4double coste = std::cos(thetaEle); <<
346 <<
347 G4double thetaPos = u*electron_mass_c2/posit <<
348 G4double sintp = std::sin(thetaPos); <<
349 G4double costp = std::cos(thetaPos); <<
350 <<
351 G4double phi = twopi * rndmEngine->flat(); <<
352 G4double sinp = std::sin(phi); <<
353 G4double cosp = std::cos(phi); <<
354 <<
355 // Kinematics of the created pair: <<
356 // the electron and positron are assumed to <<
357 // distribution with respect to the Z axis a <<
358 397
359 G4double electronKineEnergy = electronTotEne << 398 // G4double t = 0.5*asinh(momentumThreshold_N);
>> 399 G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
>> 400 //G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
>> 401 G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
>> 402 G4double J2 = (-2./3.)*log(2.*sinh(t)) + t*cosh(t)/sinh(t) + (sinh(t)-t*pow(cosh(t),3))/(3.*pow(sinh(t),3));
>> 403 G4double b = 2.*(J1-J2)/J1;
>> 404
>> 405 G4double n = 1 - b/6.;
>> 406 G4double re=0.;
>> 407 re = G4UniformRand();
>> 408 G4double a = 0.;
>> 409
>> 410 G4double b1 = 16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) +
>> 411 6.*pow(b,2.)*re*n;
>> 412 a = pow((b1/b),0.5);
>> 413 G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
>> 414 epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
>> 415
>> 416 G4double photonEnergy1 = photonEnergy - ener_re ; // resto al foton la energia del electron de retro.
>> 417 positronTotEnergy = epsilon*photonEnergy1;
>> 418 electronTotEnergy = photonEnergy1 - positronTotEnergy; // temporarly
>> 419
>> 420 G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
>> 421 electron_mass_c2*electron_mass_c2) ;
>> 422 G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
>> 423 electron_mass_c2*electron_mass_c2) ;
>> 424
>> 425 thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
>> 426 thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
>> 427 phi = twopi * G4UniformRand();
>> 428
>> 429 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
>> 430 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
>> 431
>> 432 // Kinematics of the created pair:
>> 433
>> 434 // the electron and positron are assumed to have a symetric angular
>> 435 // distribution with respect to the Z axis along the parent photon
>> 436
360 437
361 G4ThreeVector electronDirection (sinte*cosp, << 438 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
>> 439
>> 440
>> 441 G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
362 electronDirection.rotateUz(photonDirection); 442 electronDirection.rotateUz(photonDirection);
363 443
364 G4DynamicParticle* particle1 = new G4Dynamic 444 G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
365 445 electronDirection,
366 electronKineEnergy); 446 electronKineEnergy);
367 447
368 G4double positronKineEnergy = positronTotEne << 448 // The e+ is always created (even with kinetic energy = 0) for further annihilation
>> 449
>> 450 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
369 451
370 G4ThreeVector positronDirection (-sintp*cosp << 452 G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
371 positronDirection.rotateUz(photonDirection); 453 positronDirection.rotateUz(photonDirection);
372 454
373 // Create G4DynamicParticle object for the p 455 // Create G4DynamicParticle object for the particle2
>> 456
374 G4DynamicParticle* particle2 = new G4Dynamic 457 G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
375 458 positronDirection, positronKineEnergy);
376 // Fill output vector 459 // Fill output vector
377 460
378 fvect->push_back(particle1); 461 fvect->push_back(particle1);
379 fvect->push_back(particle2); 462 fvect->push_back(particle2);
380 463
381 if(particle3) { fvect->push_back(particle3); <<
382 464
383 // kill incident photon << 465 // kill incident photon
384 fParticleChange->SetProposedKineticEnergy(0. 466 fParticleChange->SetProposedKineticEnergy(0.);
385 fParticleChange->ProposeTrackStatus(fStopAnd 467 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 468
>> 469
>> 470
386 } 471 }
387 472
388 //....oooOO0OOooo........oooOO0OOooo........oo 473 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
389 474
390 #include "G4AutoLock.hh" 475 #include "G4AutoLock.hh"
391 namespace { G4Mutex BoldyshevTripletModelMutex 476 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; }
392 477
393 void G4BoldyshevTripletModel::InitialiseForEle 478 void G4BoldyshevTripletModel::InitialiseForElement(
394 const G4ParticleDefinition*, G4int Z) << 479 const G4ParticleDefinition*,
>> 480 G4int Z)
395 { 481 {
396 G4AutoLock l(&BoldyshevTripletModelMutex); 482 G4AutoLock l(&BoldyshevTripletModelMutex);
397 // G4cout << "G4BoldyshevTripletModel::Initi << 483 // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= "
398 // << Z << G4endl; << 484 // << Z << G4endl;
399 if(!data[Z]) { ReadData(Z); } 485 if(!data[Z]) { ReadData(Z); }
400 l.unlock(); 486 l.unlock();
401 } 487 }
402 488
403 //....oooOO0OOooo........oooOO0OOooo........oo 489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404 490