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