| 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 // $Id: G4PenelopeBremsstrahlungAngular.cc,v 1.7 2006/06/29 19:40:37 gunter Exp $
>> 27 // GEANT4 tag $Name: geant4-09-01-patch-03 $
>> 28 //
27 // ------------------------------------------- 29 // --------------------------------------------------------------
28 // 30 //
29 // File name: G4PenelopeBremsstrahlungAngu 31 // File name: G4PenelopeBremsstrahlungAngular
30 // 32 //
31 // Author: Luciano Pandola 33 // Author: Luciano Pandola
32 // << 34 //
33 // Creation date: November 2010 << 35 // Creation date: February 2003
34 // 36 //
35 // History: 37 // History:
36 // ----------- 38 // -----------
37 // 23 Nov 2010 L. Pandola 1st implement << 39 // 04 Feb 2003 L. Pandola 1st implementation
38 // 24 May 2011 L. Pandola Renamed (make << 40 // 19 Mar 2003 L. Pandola Bugs fixed
39 // 13 Mar 2012 L. Pandola Made a derive << 41 // 07 Nov 2003 L. Pandola Added GetAtomicNumber method for testing
40 // and update th << 42 // purposes
41 // 18 Jul 2012 L. Pandola Migrated to t <<
42 // Now returns a <<
43 // 03 Oct 2013 L. Pandola Migrated to M <<
44 // 17 Oct 2013 L. Pandola Partially rev <<
45 // thread-local <<
46 // <<
47 //-------------------------------------------- 43 //----------------------------------------------------------------
48 44
49 #include "G4PenelopeBremsstrahlungAngular.hh" 45 #include "G4PenelopeBremsstrahlungAngular.hh"
50 << 46 #include "G4PenelopeInterpolator.hh"
51 #include "globals.hh" <<
52 #include "G4PhysicalConstants.hh" <<
53 #include "G4SystemOfUnits.hh" <<
54 #include "G4PhysicsFreeVector.hh" <<
55 #include "G4PhysicsTable.hh" <<
56 #include "G4Material.hh" <<
57 #include "Randomize.hh" 47 #include "Randomize.hh"
58 #include "G4Exp.hh" << 48 #include "globals.hh"
59 <<
60 G4PenelopeBremsstrahlungAngular::G4PenelopeBre <<
61 G4VEmAngularDistribution("Penelope"), fEffec <<
62 fLorentzTables1(nullptr),fLorentzTables2(nul <<
63 49
>> 50 G4PenelopeBremsstrahlungAngular::G4PenelopeBremsstrahlungAngular (G4int Zed)
>> 51 : Zmat(Zed)
64 { 52 {
65 fDataRead = false; << 53 InterpolationTableForZ();
66 fVerbosityLevel = 0; << 54 InterpolationForK();
67 } 55 }
68 56
69 //....oooOO0OOooo........oooOO0OOooo........oo <<
70 57
71 G4PenelopeBremsstrahlungAngular::~G4PenelopeBr 58 G4PenelopeBremsstrahlungAngular::~G4PenelopeBremsstrahlungAngular()
72 { 59 {
73 ClearTables(); <<
74 } 60 }
75 61
76 //....oooOO0OOooo........oooOO0OOooo........oo << 62 G4int G4PenelopeBremsstrahlungAngular::GetAtomicNumber()
77 <<
78 void G4PenelopeBremsstrahlungAngular::Initiali <<
79 { 63 {
80 ClearTables(); << 64 return Zmat;
81 } 65 }
82 66
83 //....oooOO0OOooo........oooOO0OOooo........oo << 67 void G4PenelopeBremsstrahlungAngular::InterpolationTableForZ()
84 <<
85 void G4PenelopeBremsstrahlungAngular::ClearTab <<
86 { 68 {
87 if (fLorentzTables1) << 69 G4double pZ[NumberofZPoints] = {2.0,8.0,13.0,47.0,79.0,92.0};
88 { << 70 G4double pX[NumberofZPoints],pY[NumberofZPoints];
89 for (auto j=fLorentzTables1->begin(); j << 71 G4double QQ1[NumberofZPoints][NumberofEPoints][NumberofKPoints];
90 { << 72 G4double QQ2[NumberofZPoints][NumberofEPoints][NumberofKPoints];
91 G4PhysicsTable* tab = j->second; << 73
92 tab->clearAndDestroy(); << 74 //Read information from DataBase file
93 delete tab; << 75 char* path = getenv("G4LEDATA");
94 } <<
95 fLorentzTables1->clear(); <<
96 delete fLorentzTables1; <<
97 fLorentzTables1 = nullptr; <<
98 } <<
99 <<
100 if (fLorentzTables2) <<
101 { <<
102 for (auto j=fLorentzTables2->begin(); j <<
103 { <<
104 G4PhysicsTable* tab = j->second; <<
105 tab->clearAndDestroy(); <<
106 delete tab; <<
107 } <<
108 fLorentzTables2->clear(); <<
109 delete fLorentzTables2; <<
110 fLorentzTables2 = nullptr; <<
111 } <<
112 if (fEffectiveZSq) <<
113 { <<
114 delete fEffectiveZSq; <<
115 fEffectiveZSq = nullptr; <<
116 } <<
117 } <<
118 <<
119 //....oooOO0OOooo........oooOO0OOooo........oo <<
120 <<
121 void G4PenelopeBremsstrahlungAngular::ReadData <<
122 { <<
123 //Read information from DataBase file <<
124 const char* path = G4FindDataDir("G4LEDATA") <<
125 if (!path) 76 if (!path)
126 { 77 {
127 G4String excep = << 78 G4String excep = "G4PenelopeBremsstrahlungAngular - G4LEDATA environment variable not set!";
128 "G4PenelopeBremsstrahlungAngular - G4LEDATA << 79 G4Exception(excep);
129 G4Exception("G4PenelopeBremsstrahlungAng <<
130 "em0006",FatalException,excep); <<
131 return; <<
132 } 80 }
133 G4String pathString(path); 81 G4String pathString(path);
134 G4String pathFile = pathString + "/penelope/ << 82 G4String pathFile = pathString + "/penelope/br-ang-pen.dat";
135 std::ifstream file(pathFile); 83 std::ifstream file(pathFile);
136 << 84 std::filebuf* lsdp = file.rdbuf();
137 if (!file.is_open()) << 85
>> 86 if (!(lsdp->is_open()))
138 { 87 {
139 G4String excep = "G4PenelopeBremsstrahlu 88 G4String excep = "G4PenelopeBremsstrahlungAngular - data file " + pathFile + " not found!";
140 G4Exception("G4PenelopeBremsstrahlungAng << 89 G4Exception(excep);
141 "em0003",FatalException,excep); <<
142 return; <<
143 } 90 }
144 G4int i=0,j=0,k=0; // i=index for Z, j=index << 91 G4int i=0,j=0,k=0; // i=index for Z, j=index for E, k=index for K
145 << 92 G4double a1,a2;
146 for (k=0;k<fNumberofKPoints;k++) << 93 while(i != -1) {
147 for (i=0;i<fNumberofZPoints;i++) << 94 file >> i >> j >> k >> a1 >> a2;
148 for (j=0;j<fNumberofEPoints;j++) << 95 if (i > -1){
149 { << 96 QQ1[i][j][k]=a1;
150 G4double a1,a2; << 97 QQ2[i][j][k]=a2;
151 G4int ik1,iz1,ie1; << 98 }
152 G4double zr,er,kr; << 99 }
153 file >> iz1 >> ie1 >> ik1 >> zr >> er >> k <<
154 //check the data are correct <<
155 if ((iz1-1 == i) && (ik1-1 == k) && (ie1-1 <<
156 { <<
157 fQQ1[i][j][k]=a1; <<
158 fQQ2[i][j][k]=a2; <<
159 } <<
160 else <<
161 { <<
162 G4ExceptionDescription ed; <<
163 ed << "Corrupted data file " << pathFi <<
164 G4Exception("G4PenelopeBremsstrahlungA <<
165 "em0005",FatalException,ed); <<
166 } <<
167 } <<
168 file.close(); 100 file.close();
169 fDataRead = true; <<
170 } <<
171 <<
172 //....oooOO0OOooo........oooOO0OOooo........oo <<
173 <<
174 void G4PenelopeBremsstrahlungAngular::PrepareT <<
175 { <<
176 //Unused at the moment: the G4PenelopeBremss <<
177 //builds its own version of the tables. <<
178 101
179 //Check if data file has already been read <<
180 if (!fDataRead) <<
181 { <<
182 ReadDataFile(); <<
183 if (!fDataRead) <<
184 G4Exception("G4PenelopeBremsstrahlungAngular <<
185 "em2001",FatalException,"Unable to bui <<
186 } <<
187 <<
188 if (!fLorentzTables1) <<
189 fLorentzTables1 = new std::map<G4double, <<
190 if (!fLorentzTables2) <<
191 fLorentzTables2 = new std::map<G4double,G4 <<
192 <<
193 G4double Zmat = CalculateEffectiveZ(material <<
194 <<
195 const G4int reducedEnergyGrid=21; <<
196 //Support arrays. <<
197 G4double betas[fNumberofEPoints]; //betas fo <<
198 //tables for interpolation <<
199 G4double Q1[fNumberofEPoints][fNumberofKPoin <<
200 G4double Q2[fNumberofEPoints][fNumberofKPoin <<
201 //expanded tables for interpolation <<
202 G4double Q1E[fNumberofEPoints][reducedEnergy <<
203 G4double Q2E[fNumberofEPoints][reducedEnergy <<
204 G4double pZ[fNumberofZPoints] = {2.0,8.0,13. <<
205 102
206 G4int i=0,j=0,k=0; // i=index for Z, j=index <<
207 //Interpolation in Z 103 //Interpolation in Z
208 for (i=0;i<fNumberofEPoints;i++) << 104 for (i=0;i<NumberofEPoints;i++){
209 { << 105 for (j=0;j<NumberofKPoints;j++){
210 for (j=0;j<fNumberofKPoints;j++) << 106 for (k=0;k<NumberofZPoints;k++){
211 { << 107 pX[k]=std::log(QQ1[k][i][j]);
212 G4PhysicsFreeVector* QQ1vector = << 108 pY[k]=QQ2[k][i][j];
213 new G4PhysicsFreeVector(fNumberofZPoints << 109 }
214 G4PhysicsFreeVector* QQ2vector = << 110 G4PenelopeInterpolator* interpolator1 = new G4PenelopeInterpolator(pZ,pX,NumberofZPoints);
215 new G4PhysicsFreeVector(fNumberofZPoints << 111 Q1[i][j]=std::exp(interpolator1->CubicSplineInterpolation((G4double) Zmat));
216 << 112 delete interpolator1;
217 //fill vectors << 113 G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pZ,pY,NumberofZPoints);
218 for (k=0;k<fNumberofZPoints;k++) << 114 Q2[i][j]=interpolator2->CubicSplineInterpolation((G4double) Zmat);
219 { << 115 delete interpolator2;
220 QQ1vector->PutValues(k,pZ[k],G4Log(fQQ <<
221 QQ2vector->PutValues(k,pZ[k],fQQ2[k][i <<
222 } <<
223 //Filled: now calculate derivatives <<
224 QQ1vector->FillSecondDerivatives(); <<
225 QQ2vector->FillSecondDerivatives(); <<
226 <<
227 Q1[i][j]= G4Exp(QQ1vector->Value(Zmat)); <<
228 Q2[i][j]=QQ2vector->Value(Zmat); <<
229 delete QQ1vector; <<
230 delete QQ2vector; <<
231 } <<
232 } 116 }
233 G4double pE[fNumberofEPoints] = {1.0e-03*MeV << 117 }
234 1.0e-01*MeV,5.0e-01*MeV}; << 118
235 G4double pK[fNumberofKPoints] = {0.0,0.6,0.8 << 119
>> 120 //std::ofstream fil("matrice.dat",std::ios::app);
>> 121 //fil << "Numero atomico: " << Zmat << G4endl;
>> 122 //for (i=0;i<NumberofEPoints;i++)
>> 123 //{
>> 124 // fil << Q1[i][0] << " " << Q1[i][1] << " " << Q1[i][2] << " " << Q1[i][3] << G4endl;
>> 125 //}
>> 126 //fil.close();
>> 127
>> 128 }
>> 129
>> 130 void G4PenelopeBremsstrahlungAngular::InterpolationForK()
>> 131 {
>> 132 G4double pE[NumberofEPoints] = {1.0e-03,5.0e-03,1.0e-02,5.0e-02,1.0e-01,5.0e-01};
>> 133 G4double pK[NumberofKPoints] = {0.0,0.6,0.8,0.95};
236 G4double ppK[reducedEnergyGrid]; 134 G4double ppK[reducedEnergyGrid];
>> 135 G4double pX[NumberofKPoints];
>> 136 G4int i,j;
237 137
238 for(i=0;i<reducedEnergyGrid;i++) << 138 for(i=0;i<reducedEnergyGrid;i++){
239 ppK[i]=((G4double) i) * 0.05; 139 ppK[i]=((G4double) i) * 0.05;
>> 140 }
240 141
241 << 142 for(i=0;i<NumberofEPoints;i++){
242 for(i=0;i<fNumberofEPoints;i++) <<
243 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron 143 betas[i]=std::sqrt(pE[i]*(pE[i]+2*electron_mass_c2))/(pE[i]+electron_mass_c2);
>> 144 }
244 145
245 << 146 for (i=0;i<NumberofEPoints;i++){
246 for (i=0;i<fNumberofEPoints;i++) << 147 for (j=0;j<NumberofKPoints;j++){
247 { << 148 Q1[i][j]=Q1[i][j]/((G4double) Zmat);
248 for (j=0;j<fNumberofKPoints;j++) <<
249 Q1[i][j]=Q1[i][j]/Zmat; <<
250 } 149 }
>> 150 }
251 151
252 //Expanded table of distribution parameters 152 //Expanded table of distribution parameters
253 for (i=0;i<fNumberofEPoints;i++) << 153 for (i=0;i<NumberofEPoints;i++){
254 { << 154 for (j=0;j<NumberofKPoints;j++){
255 G4PhysicsFreeVector* Q1vector = new G4Ph << 155 pX[j]=std::log(Q1[i][j]); //logarithmic
256 G4PhysicsFreeVector* Q2vector = new G4Ph <<
257 <<
258 for (j=0;j<fNumberofKPoints;j++) <<
259 { <<
260 Q1vector->PutValues(j,pK[j],G4Log(Q1[i][j] <<
261 Q2vector->PutValues(j,pK[j],Q2[i][j]); <<
262 } <<
263 <<
264 for (j=0;j<reducedEnergyGrid;j++) <<
265 { <<
266 Q1E[i][j]=Q1vector->Value(ppK[j]); <<
267 Q2E[i][j]=Q2vector->Value(ppK[j]); <<
268 } <<
269 delete Q1vector; <<
270 delete Q2vector; <<
271 } <<
272 // <<
273 //TABLES to be stored <<
274 // <<
275 G4PhysicsTable* theTable1 = new G4PhysicsTab <<
276 G4PhysicsTable* theTable2 = new G4PhysicsTab <<
277 // the table will contain reducedEnergyGrid <<
278 // values of k, <<
279 // Each of the G4PhysicsFreeVectors has a pr <<
280 // y vs. E <<
281 // <<
282 //reserve space of the vectors. <<
283 for (j=0;j<reducedEnergyGrid;j++) <<
284 { <<
285 G4PhysicsFreeVector* thevec = new G4Phys <<
286 theTable1->push_back(thevec); <<
287 G4PhysicsFreeVector* thevec2 = new G4Phy <<
288 theTable2->push_back(thevec2); <<
289 } 156 }
290 << 157 G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(pK,pX,NumberofKPoints);
291 for (j=0;j<reducedEnergyGrid;j++) << 158 for (j=0;j<reducedEnergyGrid;j++){
292 { << 159 Q1E[i][j]=interpolator->CubicSplineInterpolation(ppK[j]);
293 G4PhysicsFreeVector* thevec = (G4Physics <<
294 G4PhysicsFreeVector* thevec2 = (G4Physic <<
295 for (i=0;i<fNumberofEPoints;i++) <<
296 { <<
297 thevec->PutValues(i,betas[i],Q1E[i][j]); <<
298 thevec2->PutValues(i,betas[i],Q2E[i][j]); <<
299 } <<
300 //Vectors are filled: calculate derivati <<
301 thevec->FillSecondDerivatives(); <<
302 thevec2->FillSecondDerivatives(); <<
303 } 160 }
304 << 161 delete interpolator;
305 if (fLorentzTables1 && fLorentzTables2) << 162 for (j=0;j<NumberofKPoints;j++){
306 { << 163 pX[j]=Q2[i][j];
307 fLorentzTables1->insert(std::make_pair(Z <<
308 fLorentzTables2->insert(std::make_pair(Z <<
309 } 164 }
310 else << 165 G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(pK,pX,NumberofKPoints);
311 { << 166 for (j=0;j<reducedEnergyGrid;j++){
312 G4ExceptionDescription ed; << 167 Q2E[i][j]=interpolator2->CubicSplineInterpolation(ppK[j]);
313 ed << "Unable to create tables of Lorent <<
314 ed << "<Z>= " << Zmat << " in G4Penelop <<
315 delete theTable1; <<
316 delete theTable2; <<
317 G4Exception("G4PenelopeBremsstrahlungAng <<
318 "em2005",FatalException,ed); <<
319 } 168 }
320 << 169 delete interpolator2;
321 return; << 170 }
322 } 171 }
323 172
324 //....oooOO0OOooo........oooOO0OOooo........oo << 173 G4double G4PenelopeBremsstrahlungAngular::ExtractCosTheta(G4double e1,G4double e2)
325 <<
326 G4ThreeVector& G4PenelopeBremsstrahlungAngular <<
327 G4double eGamma, <<
328 G4int, <<
329 const G4Material* material) <<
330 { 174 {
331 if (!material) << 175 //e1 = kinetic energy of the electron
332 { << 176 //e2 = energy of the bremsstrahlung photon
333 G4Exception("G4PenelopeBremsstrahlungAng <<
334 "em2040",FatalException,"The pointer to <<
335 return fLocalDirection; <<
336 } <<
337 <<
338 //Retrieve the effective Z <<
339 G4double Zmat = 0; <<
340 <<
341 //The model might be initialized incorrectly <<
342 //G4PenelopeBremsstrahungModel: make sure it <<
343 if (!fEffectiveZSq) <<
344 { <<
345 G4Exception("G4PenelopeBremsstrahlungAng <<
346 "em2040",JustWarning,"EffectiveZSq table <<
347 PrepareTables(material,false); <<
348 //return fLocalDirection; <<
349 } <<
350 <<
351 //found in the table: return it <<
352 if (fEffectiveZSq->count(material)) <<
353 Zmat = fEffectiveZSq->find(material)->seco <<
354 else //this can happen in unit tests or when <<
355 //models other than G4PenelopeBremsstrahun <<
356 { <<
357 G4Exception("G4PenelopeBremsstrahlungAng <<
358 "em2040",JustWarning,"Material not found <<
359 PrepareTables(material,false); <<
360 Zmat = fEffectiveZSq->find(material)->se <<
361 // return fLocalDirection; <<
362 } <<
363 <<
364 if (fVerbosityLevel > 0) <<
365 { <<
366 G4cout << "Effective <Z> for material : <<
367 " = " << Zmat << G4endl; <<
368 } <<
369 <<
370 G4double ePrimary = dp->GetKineticEnergy(); <<
371 <<
372 G4double beta = std::sqrt(ePrimary*(ePrimary <<
373 (ePrimary+electron_mass_c2); <<
374 G4double cdt = 0; <<
375 G4double sinTheta = 0; <<
376 G4double phi = 0; <<
377 177
378 //Use a pure dipole distribution for energy << 178 G4double beta = std::sqrt(e1*(e1+2*electron_mass_c2))/(e1+electron_mass_c2);
379 if (ePrimary > 500*keV) << 179
380 { <<
381 cdt = 2.0*G4UniformRand() - 1.0; <<
382 if (G4UniformRand() > 0.75) <<
383 { <<
384 if (cdt<0) <<
385 cdt = -1.0*std::pow(-cdt,1./3.); <<
386 else <<
387 cdt = std::pow(cdt,1./3.); <<
388 } <<
389 cdt = (cdt+beta)/(1.0+beta*cdt); <<
390 //Get primary kinematics <<
391 sinTheta = std::sqrt(1. - cdt*cdt); <<
392 phi = twopi * G4UniformRand(); <<
393 fLocalDirection.set(sinTheta* std::cos(p <<
394 sinTheta* std::sin(phi), <<
395 cdt); <<
396 //rotate <<
397 fLocalDirection.rotateUz(dp->GetMomentum <<
398 //return <<
399 return fLocalDirection; <<
400 } <<
401 <<
402 if (!(fLorentzTables1->count(Zmat)) || !(fLo <<
403 { <<
404 G4ExceptionDescription ed; <<
405 ed << "Unable to retrieve Lorentz tables <<
406 G4Exception("G4PenelopeBremsstrahlungAng <<
407 "em2006",FatalException,ed); <<
408 } <<
409 180
410 //retrieve actual tables <<
411 const G4PhysicsTable* theTable1 = fLorentzTa <<
412 const G4PhysicsTable* theTable2 = fLorentzTa <<
413 181
414 G4double RK=20.0*eGamma/ePrimary; << 182 G4double RK=20.0*e2/e1;
415 G4int ik=std::min((G4int) RK,19); 183 G4int ik=std::min((G4int) RK,19);
416 << 184
417 G4double P10=0,P11=0,P1=0; 185 G4double P10=0,P11=0,P1=0;
418 G4double P20=0,P21=0,P2=0; 186 G4double P20=0,P21=0,P2=0;
419 << 187 G4double pX[NumberofEPoints];
420 //First coefficient 188 //First coefficient
421 const G4PhysicsFreeVector* v1 = (G4PhysicsFr << 189 G4int i;
422 const G4PhysicsFreeVector* v2 = (G4PhysicsFr << 190 G4int j = ik;
423 P10 = v1->Value(beta); << 191 for (i=0;i<NumberofEPoints;i++){
424 P11 = v2->Value(beta); << 192 pX[i]=Q1E[i][j];
>> 193 }
>> 194 G4PenelopeInterpolator* interpolator = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
>> 195 P10=interpolator->CubicSplineInterpolation(beta);
>> 196 delete interpolator;
>> 197 j++; //(j=ik+1)
>> 198 for (i=0;i<NumberofEPoints;i++){
>> 199 pX[i]=Q1E[i][j];
>> 200 }
>> 201 G4PenelopeInterpolator* interpolator2 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
>> 202 P11=interpolator2->CubicSplineInterpolation(beta);
>> 203 delete interpolator2;
425 P1=P10+(RK-(G4double) ik)*(P11-P10); 204 P1=P10+(RK-(G4double) ik)*(P11-P10);
426 << 205
427 //Second coefficient 206 //Second coefficient
428 const G4PhysicsFreeVector* v3 = (G4PhysicsFr << 207 j = ik;
429 const G4PhysicsFreeVector* v4 = (G4PhysicsFr << 208 for (i=0;i<NumberofEPoints;i++){
430 P20=v3->Value(beta); << 209 pX[i]=Q2E[i][j];
431 P21=v4->Value(beta); << 210 }
>> 211 G4PenelopeInterpolator* interpolator3 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
>> 212 P20=interpolator3->CubicSplineInterpolation(beta);
>> 213 delete interpolator3;
>> 214 j++; //(j=ik+1)
>> 215 for (i=0;i<NumberofEPoints;i++){
>> 216 pX[i]=Q2E[i][j];
>> 217 }
>> 218 G4PenelopeInterpolator* interpolator4 = new G4PenelopeInterpolator(betas,pX,NumberofEPoints);
>> 219 P21=interpolator4->CubicSplineInterpolation(beta);
>> 220 delete interpolator4;
432 P2=P20+(RK-(G4double) ik)*(P21-P20); 221 P2=P20+(RK-(G4double) ik)*(P21-P20);
433 << 222
434 //Sampling from the Lorenz-trasformed dipole 223 //Sampling from the Lorenz-trasformed dipole distributions
435 P1=std::min(G4Exp(P1)/beta,1.0); << 224 P1=std::min(std::exp(P1)/beta,1.0);
436 G4double betap = std::min(std::max(beta*(1.0 225 G4double betap = std::min(std::max(beta*(1.0+P2/beta),0.0),0.9999);
437 << 226
438 G4double testf=0; << 227 G4double cdt=0,testf=0;
439 << 228
440 if (G4UniformRand() < P1) << 229 if (G4UniformRand() < P1){
441 { << 230 do{
442 do{ << 231 cdt = 2.0*G4UniformRand()-1.0;
443 cdt = 2.0*G4UniformRand()-1.0; << 232 testf=2.0*G4UniformRand()-(1.0+cdt*cdt);
444 testf=2.0*G4UniformRand()-(1.0+cdt*cdt); << 233 }while(testf>0);
445 }while(testf>0); << 234 }
446 } << 235 else{
447 else << 236 do{
448 { << 237 cdt = 2.0*G4UniformRand()-1.0;
449 do{ << 238 testf=G4UniformRand()-(1.0-cdt*cdt);
450 cdt = 2.0*G4UniformRand()-1.0; << 239 }while(testf>0);
451 testf=G4UniformRand()-(1.0-cdt*cdt); << 240 }
452 }while(testf>0); <<
453 } <<
454 cdt = (cdt+betap)/(1.0+betap*cdt); 241 cdt = (cdt+betap)/(1.0+betap*cdt);
455 << 242 return cdt;
456 //Get primary kinematics <<
457 sinTheta = std::sqrt(1. - cdt*cdt); <<
458 phi = twopi * G4UniformRand(); <<
459 fLocalDirection.set(sinTheta* std::cos(phi), <<
460 sinTheta* std::sin(phi), <<
461 cdt); <<
462 //rotate <<
463 fLocalDirection.rotateUz(dp->GetMomentumDire <<
464 //return <<
465 return fLocalDirection; <<
466 <<
467 } <<
468 <<
469 //....oooOO0OOooo........oooOO0OOooo........oo <<
470 <<
471 G4double G4PenelopeBremsstrahlungAngular::Calc <<
472 { <<
473 if (!fEffectiveZSq) <<
474 fEffectiveZSq = new std::map<const G4Mater <<
475 <<
476 //Already exists: return it <<
477 if (fEffectiveZSq->count(material)) <<
478 return fEffectiveZSq->find(material)->seco <<
479 <<
480 //Helper for the calculation <<
481 std::vector<G4double> *StechiometricFactors <<
482 G4int nElements = (G4int)material->GetNumber <<
483 const G4ElementVector* elementVector = mater <<
484 const G4double* fractionVector = material->G <<
485 for (G4int i=0;i<nElements;++i) <<
486 { <<
487 G4double fraction = fractionVector[i]; <<
488 G4double atomicWeigth = (*elementVector) <<
489 StechiometricFactors->push_back(fraction <<
490 } <<
491 //Find max <<
492 G4double MaxStechiometricFactor = 0.; <<
493 for (G4int i=0;i<nElements;++i) <<
494 { <<
495 if ((*StechiometricFactors)[i] > MaxStec <<
496 MaxStechiometricFactor = (*Stechiometr <<
497 } <<
498 //Normalize <<
499 for (G4int i=0;i<nElements;++i) <<
500 (*StechiometricFactors)[i] /= MaxStechiom <<
501 <<
502 G4double sumz2 = 0; <<
503 G4double sums = 0; <<
504 for (G4int i=0;i<nElements;++i) <<
505 { <<
506 G4double Z = (*elementVector)[i]->GetZ() <<
507 sumz2 += (*StechiometricFactors)[i]*Z*Z; <<
508 sums += (*StechiometricFactors)[i]; <<
509 } <<
510 delete StechiometricFactors; <<
511 <<
512 G4double ZBR = std::sqrt(sumz2/sums); <<
513 fEffectiveZSq->insert(std::make_pair(materia <<
514 <<
515 return ZBR; <<
516 } 243 }
517 244