| 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 // 10.06.97 created. V. Grichine << 26 //
27 // 18.11.98 simplified public interface; new m << 27 // $Id: G4SandiaTable.cc,v 1.27 2007/06/20 09:21:45 vnivanch Exp $
28 // 31.01.01 redesign of ComputeMatSandiaMatrix << 28 // GEANT4 tag $Name: geant4-09-00-patch-01 $
29 // 16.02.01 adapted for STL. mma << 29 //
30 // 22.02.01 GetsandiaCofForMaterial(energy) re << 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
31 // 03.04.01 fnulcof returned if energy < emin << 31 //
32 // 10.07.01 Migration to STL. M. Verderi. <<
33 // 03.02.04 Update distructor V.Ivanchenko <<
34 // 05.03.04 New methods for old sorting algori 32 // 05.03.04 New methods for old sorting algorithm for PAI model. V.Grichine
35 // 26.10.11 new scheme for G4Exception (mma) << 33 // 03.02.04 Update distructor V.Ivanchenko
36 // 22.05.13 preparation of material table with << 34 // 10.07.01 Migration to STL. M. Verderi.
37 // 09.07.14 modify low limit in GetSandiaCofPe << 35 // 03.04.01 fnulcof returned if energy < emin
38 // 10.07.14 modify low limit for water. VI << 36 // 22.02.01 GetsandiaCofForMaterial(energy) return 0 below lowest interval mma
>> 37 // 16.02.01 adapted for STL. mma
>> 38 // 31.01.01 redesign of ComputeMatSandiaMatrix(). mma
>> 39 // 18.11.98 simplified public interface; new methods for materials. mma
>> 40 // 10.06.97 created. V. Grichine
>> 41 //
>> 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
39 43
40 #include "G4SandiaTable.hh" <<
41 44
>> 45 #include "G4SandiaTable.hh"
>> 46 #include "G4StaticSandiaData.hh"
42 #include "G4Material.hh" 47 #include "G4Material.hh"
43 #include "G4MaterialTable.hh" 48 #include "G4MaterialTable.hh"
44 #include "G4PhysicalConstants.hh" <<
45 #include "G4StaticSandiaData.hh" <<
46 #include "G4SystemOfUnits.hh" <<
47 <<
48 const G4double G4SandiaTable::funitc[5] = {CLH <<
49 CLHEP::cm2* CLHEP::keV* CLHEP::keV / CLHEP:: <<
50 CLHEP::cm2* CLHEP::keV* CLHEP::keV* CLHEP::k <<
51 CLHEP::cm2* CLHEP::keV* CLHEP::keV* CLHEP::k <<
52 <<
53 G4int G4SandiaTable::fCumulInterval[] = {0}; <<
54 49
>> 50 G4int G4SandiaTable::fCumulInterval[101];
>> 51 G4double G4SandiaTable::fSandiaCofPerAtom[4];
>> 52
55 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
56 54
57 G4SandiaTable::G4SandiaTable(const G4Material* << 55 G4SandiaTable::G4SandiaTable(G4int matIndex)
58 { << 56 {
59 fMatSandiaMatrix = nullptr; << 57 fMatSandiaMatrix = 0 ;
60 fMatSandiaMatrixPAI = nullptr; << 58 fPhotoAbsorptionCof = 0 ;
61 fPhotoAbsorptionCof = nullptr; <<
62 <<
63 fMatNbOfIntervals = 0; <<
64 <<
65 fMaxInterval = 0; <<
66 fVerbose = 0; <<
67 <<
68 // build the CumulInterval array <<
69 if (0 == fCumulInterval[0]) { <<
70 fCumulInterval[0] = 1; <<
71 <<
72 for (G4int Z = 1; Z < 101; ++Z) { <<
73 fCumulInterval[Z] = fCumulInterval[Z - 1 <<
74 } <<
75 } <<
76 <<
77 fMaxInterval = 0; <<
78 fSandiaCofPerAtom.resize(4, 0.0); <<
79 fLowerI1 = false; <<
80 // compute macroscopic Sandia coefs for a ma <<
81 ComputeMatSandiaMatrix(); // mma <<
82 } <<
83 <<
84 //....oooOO0OOooo........oooOO0OOooo........oo <<
85 <<
86 G4SandiaTable::~G4SandiaTable() <<
87 { <<
88 if (fMatSandiaMatrix != nullptr) { <<
89 fMatSandiaMatrix->clearAndDestroy(); <<
90 delete fMatSandiaMatrix; <<
91 } <<
92 if (fMatSandiaMatrixPAI != nullptr) { <<
93 fMatSandiaMatrixPAI->clearAndDestroy(); <<
94 delete fMatSandiaMatrixPAI; <<
95 } <<
96 <<
97 delete[] fPhotoAbsorptionCof; <<
98 } <<
99 59
100 //....oooOO0OOooo........oooOO0OOooo........oo << 60 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
>> 61 size_t numberOfMat = G4Material::GetNumberOfMaterials();
101 62
102 void G4SandiaTable::GetSandiaCofPerAtom( << 63 if ( matIndex >= 0 && matIndex < G4int(numberOfMat) )
103 G4int Z, G4double energy, std::vector<G4doub << 64 {
104 { << 65 fMaterial = (*theMaterialTable)[matIndex];
105 #ifdef G4VERBOSE <<
106 if (Z < 1 || Z > 100) { <<
107 Z = PrintErrorZ(Z, "GetSandiaCofPerAtom"); <<
108 } <<
109 if (4 > coeff.size()) { <<
110 PrintErrorV("GetSandiaCofPerAtom(): input <<
111 coeff.resize(4); <<
112 } <<
113 #endif <<
114 G4double Emin = fSandiaTable[fCumulInterval[ <<
115 // G4double Iopot = fIonizationPotentials[Z] <<
116 // if (Emin < Iopot) Emin = Iopot; <<
117 <<
118 G4int row = 0; <<
119 if (energy <= Emin) { <<
120 energy = Emin; <<
121 } 66 }
122 else { << 67 else
123 G4int interval = fNbOfIntervals[Z] - 1; << 68 {
124 row = fCumulInterval[Z - 1] + interval; << 69 G4Exception("G4SandiaTable::G4SandiaTable(G4int matIndex): wrong matIndex ");
125 // Loop checking, 07-Aug-2015, Vladimir Iv <<
126 while ((interval > 0) && (energy < fSandia <<
127 --interval; <<
128 row = fCumulInterval[Z - 1] + interval; <<
129 } <<
130 } 70 }
131 << 71 ComputeMatTable();
132 G4double AoverAvo = Z * amu / fZtoAratio[Z]; <<
133 <<
134 coeff[0] = AoverAvo * funitc[1] * fSandiaTab <<
135 coeff[1] = AoverAvo * funitc[2] * fSandiaTab <<
136 coeff[2] = AoverAvo * funitc[3] * fSandiaTab <<
137 coeff[3] = AoverAvo * funitc[4] * fSandiaTab <<
138 } 72 }
139 73
140 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
141 75
142 void G4SandiaTable::GetSandiaCofWater(G4double << 76 G4SandiaTable::G4SandiaTable(G4Material* material)
>> 77 : fMaterial(material)
143 { 78 {
144 #ifdef G4VERBOSE << 79 fMatSandiaMatrix = 0 ;
145 if (4 > coeff.size()) { << 80 fPhotoAbsorptionCof = 0 ;
146 PrintErrorV("GetSandiaCofWater: input vect << 81 //build the CumulInterval array
147 coeff.resize(4); << 82 fCumulInterval[0] = 1;
148 } << 83 for (G4int Z=1; Z<101; Z++)
149 #endif << 84 fCumulInterval[Z] = fCumulInterval[Z-1] + fNbOfIntervals[Z];
150 G4int i = 0; << 85
151 if (energy > fH2OlowerI1[0][0] * CLHEP::keV) << 86 //compute macroscopic Sandia coefs for a material
152 i = fH2OlowerInt - 1; << 87 ComputeMatSandiaMatrix();
153 for (; i > 0; --i) { << 88
154 if (energy >= fH2OlowerI1[i][0] * CLHEP: << 89 //initialisation of fnulcof
155 break; << 90 fnulcof[0] = fnulcof[1] = fnulcof[2] = fnulcof[3] = 0.;
156 } <<
157 } <<
158 } <<
159 coeff[0] = funitc[1] * fH2OlowerI1[i][1]; <<
160 coeff[1] = funitc[2] * fH2OlowerI1[i][2]; <<
161 coeff[2] = funitc[3] * fH2OlowerI1[i][3]; <<
162 coeff[3] = funitc[4] * fH2OlowerI1[i][4]; <<
163 } 91 }
164 << 92
165 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
166 94
167 G4double G4SandiaTable::GetWaterEnergyLimit() << 95 // Fake default constructor - sets only member data and allocates memory
168 { << 96 // for usage restricted to object persistency
169 return fH2OlowerI1[fH2OlowerInt - 1][0] * CL <<
170 } <<
171 <<
172 //....oooOO0OOooo........oooOO0OOooo........oo <<
173 97
174 G4double G4SandiaTable::GetWaterCofForMaterial << 98 G4SandiaTable::G4SandiaTable(__void__&)
>> 99 : fMaterial(0), fMatSandiaMatrix(0), fPhotoAbsorptionCof(0)
175 { 100 {
176 return fH2OlowerI1[i][j] * funitc[j]; <<
177 } 101 }
178 102
179 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
180 104
181 G4double G4SandiaTable::GetZtoA(G4int Z) << 105 G4SandiaTable::~G4SandiaTable()
182 { << 106 {
183 #ifdef G4VERBOSE << 107 if(fMatSandiaMatrix) {
184 if (Z < 1 || Z > 100) { << 108 // fMatSandiaMatrix->clearAndDestroy();
185 Z = PrintErrorZ(Z, "GetSandiaCofPerAtom"); << 109 delete fMatSandiaMatrix;
>> 110 }
>> 111
>> 112 if(fPhotoAbsorptionCof)
>> 113 {
>> 114 // for(G4int i = 0 ; i < fMaxInterval ; i++) delete[] fPhotoAbsorptionCof[i];
>> 115 delete [] fPhotoAbsorptionCof ;
186 } 116 }
187 #endif <<
188 return fZtoAratio[Z]; <<
189 } <<
190 <<
191 //....oooOO0OOooo........oooOO0OOooo........oo <<
192 <<
193 #ifdef G4VERBOSE <<
194 <<
195 G4int G4SandiaTable::PrintErrorZ(G4int Z, cons <<
196 { <<
197 G4String sss = "G4SandiaTable::" + ss + "()" <<
198 G4ExceptionDescription ed; <<
199 ed << "Atomic number out of range Z= " << Z <<
200 G4Exception(sss, "mat060", JustWarning, ed, <<
201 return (Z > 100) ? 100 : 1; <<
202 } <<
203 <<
204 //....oooOO0OOooo........oooOO0OOooo........oo <<
205 <<
206 void G4SandiaTable::PrintErrorV(const G4String <<
207 { <<
208 G4String sss = "G4SandiaTable::" + ss; <<
209 G4ExceptionDescription ed; <<
210 G4Exception(sss, "mat061", JustWarning, "Wro <<
211 } 117 }
212 #endif << 118
213 <<
214 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
215 120
216 void G4SandiaTable::ComputeMatSandiaMatrix() 121 void G4SandiaTable::ComputeMatSandiaMatrix()
217 { << 122 {
218 // get list of elements << 123 //get list of elements
219 const auto NbElm = (G4int)fMaterial->GetNumb << 124 const G4int NbElm = fMaterial->GetNumberOfElements();
220 const G4ElementVector* ElementVector = fMate 125 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
>> 126
>> 127 G4int* Z = new G4int[NbElm]; //Atomic number
221 128
222 auto Z = new G4int[NbElm]; // Atomic number << 129 //
223 << 130 //determine the maximum number of energy-intervals for this material
224 // determine the maximum number of energy-in << 131 //
225 G4int MaxIntervals = 0; 132 G4int MaxIntervals = 0;
226 G4int elm, z; << 133 G4int elm;
227 << 134 for (elm=0; elm<NbElm; elm++)
228 // here we compute only for a mixture, so no << 135 { Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
229 // if z is out of validity interval << 136 MaxIntervals += fNbOfIntervals[Z[elm]];
230 for (elm = 0; elm < NbElm; ++elm) { << 137 }
231 z = G4lrint((*ElementVector)[elm]->GetZ()) << 138
232 if (z < 1) { << 139 //
233 z = 1; << 140 //copy the Energy bins in a tmp1 array
234 } <<
235 else if (z > 100) { <<
236 z = 100; <<
237 } <<
238 Z[elm] = z; <<
239 MaxIntervals += fNbOfIntervals[z]; <<
240 } <<
241 <<
242 // copy the Energy bins in a tmp1 array <<
243 //(take care of the Ionization Potential of 141 //(take care of the Ionization Potential of each element)
244 auto tmp1 = new G4double[MaxIntervals]; << 142 //
>> 143 G4double* tmp1 = new G4double[MaxIntervals];
245 G4double IonizationPot; 144 G4double IonizationPot;
246 G4int interval1 = 0; << 145 G4int interval1=0;
247 << 146 for (elm=0; elm<NbElm; elm++)
248 for (elm = 0; elm < NbElm; ++elm) { << 147 {
249 z = Z[elm]; << 148 IonizationPot = GetIonizationPot(Z[elm]);
250 IonizationPot = fIonizationPotentials[z] * << 149 for (G4int row=fCumulInterval[Z[elm]-1];row<fCumulInterval[Z[elm]];row++)
251 for (G4int row = fCumulInterval[z - 1]; ro << 150 tmp1[interval1++] = std::max(fSandiaTable[row][0]*keV,IonizationPot);
252 tmp1[interval1] = std::max(fSandiaTable[ << 151 }
253 ++interval1; << 152
254 } << 153 //
255 } << 154 //sort the energies in strickly increasing values in a tmp2 array
256 // sort the energies in strickly increasing <<
257 //(eliminate redondances) 155 //(eliminate redondances)
258 << 156 //
259 auto tmp2 = new G4double[MaxIntervals]; << 157 G4double* tmp2 = new G4double[MaxIntervals];
260 G4double Emin; 158 G4double Emin;
261 G4int interval2 = 0; 159 G4int interval2 = 0;
262 << 160
263 do { 161 do {
264 Emin = DBL_MAX; << 162 Emin = DBL_MAX;
265 << 163 for (G4int i1=0; i1<MaxIntervals; i1++)
266 for (G4int i1 = 0; i1 < MaxIntervals; ++i1 << 164 if (tmp1[i1] < Emin) Emin = tmp1[i1]; //find the minimum
267 Emin = std::min(Emin, tmp1[i1]); // fin << 165 if (Emin < DBL_MAX) tmp2[interval2++] = Emin; //copy Emin in tmp2
268 } << 166 for (G4int j1=0; j1<MaxIntervals; j1++)
269 if (Emin < DBL_MAX) { << 167 if (tmp1[j1] <= Emin) tmp1[j1] = DBL_MAX; //eliminate from tmp1
270 tmp2[interval2] = Emin; << 168 } while (Emin < DBL_MAX);
271 ++interval2; << 169
272 } << 170 //
273 // copy Emin in tmp2 << 171 //create the sandia matrix for this material
274 for (G4int j1 = 0; j1 < MaxIntervals; ++j1 << 172 //
275 if (tmp1[j1] <= Emin) { <<
276 tmp1[j1] = DBL_MAX; <<
277 } // eliminate from tmp1 <<
278 } <<
279 // Loop checking, 07-Aug-2015, Vladimir Iv <<
280 } while (Emin < DBL_MAX); <<
281 <<
282 // create the sandia matrix for this materia <<
283 <<
284 fMatSandiaMatrix = new G4OrderedTable(); 173 fMatSandiaMatrix = new G4OrderedTable();
285 G4int interval; 174 G4int interval;
>> 175 for (interval=0; interval<interval2; interval++)
>> 176 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
>> 177
>> 178 //
>> 179 //ready to compute the Sandia coefs for the material
>> 180 //
>> 181 const G4double* NbOfAtomsPerVolume = fMaterial->GetVecNbOfAtomsPerVolume();
>> 182
>> 183 const G4double prec = 1.e-03*eV;
>> 184 G4double coef, oldsum(0.), newsum(0.);
>> 185 fMatNbOfIntervals = 0;
>> 186
>> 187 for (interval=0; interval<interval2; interval++)
>> 188 {
>> 189 Emin = (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[0] = tmp2[interval];
>> 190 for (G4int k=1; k<5; k++)(*(*fMatSandiaMatrix)[fMatNbOfIntervals])[k]=0.;
>> 191 newsum = 0.;
>> 192
>> 193 for (elm=0; elm<NbElm; elm++)
>> 194 {
>> 195 GetSandiaCofPerAtom(Z[elm], Emin+prec);
>> 196 for (G4int j=1; j<5; j++)
>> 197 {
>> 198 coef = NbOfAtomsPerVolume[elm]*fSandiaCofPerAtom[j-1];
>> 199 (*(*fMatSandiaMatrix)[fMatNbOfIntervals])[j] += coef;
>> 200 newsum += std::abs(coef);
>> 201 }
>> 202 }
>> 203
>> 204 //check for null or redondant intervals
>> 205 if (newsum != oldsum) { oldsum = newsum; fMatNbOfIntervals++;}
>> 206 }
>> 207 delete [] Z;
>> 208 delete [] tmp1;
>> 209 delete [] tmp2;
>> 210 }
286 211
287 for (interval = 0; interval < interval2; ++i << 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
288 fMatSandiaMatrix->push_back(new G4DataVect <<
289 } <<
290 213
291 // ready to compute the Sandia coefs for the << 214 G4double G4SandiaTable::GetSandiaCofForMaterial(G4int interval, G4int j)
>> 215 {
>> 216 assert (interval>=0 && interval<fMatNbOfIntervals && j>=0 && j<5);
>> 217 return ((*(*fMatSandiaMatrix)[interval])[j]);
>> 218 }
292 219
293 const G4double* NbOfAtomsPerVolume = fMateri << 220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
294 221
295 static const G4double prec = 1.e-03 * CLHEP: << 222 G4double* G4SandiaTable::GetSandiaCofForMaterial(G4double energy)
296 G4double coef, oldsum(0.), newsum(0.); << 223 {
297 fMatNbOfIntervals = 0; << 224 if (energy < (*(*fMatSandiaMatrix)[0])[0]) return fnulcof;
>> 225
>> 226 G4int interval = fMatNbOfIntervals - 1;
>> 227 while ((interval>0)&&(energy<(*(*fMatSandiaMatrix)[interval])[0])) interval--;
>> 228 return &((*(*fMatSandiaMatrix)[interval])[1]);
>> 229 }
298 230
299 for (interval = 0; interval < interval2; ++i << 231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... ....oooOO0OOooo....
300 Emin = (*(*fMatSandiaMatrix)[fMatNbOfInter <<
301 232
302 for (G4int k = 1; k < 5; ++k) { << 233 ////////////////////////////////////////////////////////////////////////
303 (*(*fMatSandiaMatrix)[fMatNbOfIntervals] << 234 ////////////////////////////////////////////////////////////////////////
304 } << 235 //
305 newsum = 0.; << 236 // Methods for PAI model
306 237
307 for (elm = 0; elm < NbElm; elm++) { <<
308 GetSandiaCofPerAtom(Z[elm], Emin + prec, <<
309 238
310 for (G4int j = 1; j < 5; ++j) { << 239 ///////////////////////////////////////////////////////////////////////
311 coef = NbOfAtomsPerVolume[elm] * fSand << 240 //
312 (*(*fMatSandiaMatrix)[fMatNbOfInterval << 241 // Bubble sorting of left energy interval in SandiaTable in ascening order
313 newsum += std::abs(coef); << 242 //
314 } <<
315 } <<
316 // check for null or redondant intervals <<
317 243
318 if (newsum != oldsum) { << 244 void
319 oldsum = newsum; << 245 G4SandiaTable::SandiaSort(G4double** da ,
320 ++fMatNbOfIntervals; << 246 G4int sz )
321 } << 247 {
322 } << 248 for(G4int i = 1 ;i < sz ; i++ )
323 delete[] Z; << 249 {
324 delete[] tmp1; << 250 for(G4int j = i + 1 ;j < sz ; j++ )
325 delete[] tmp2; << 251 {
326 << 252 if(da[i][0] > da[j][0])
327 if (fVerbose > 0) { << 253 {
328 G4cout << "G4SandiaTable::ComputeMatSandia << 254 SandiaSwap(da,i,j) ;
329 << 255 }
330 for (G4int i = 0; i < fMatNbOfIntervals; + << 256 }
331 G4cout << i << "\t" << GetSandiaCofForMa << 257 }
332 << GetSandiaCofForMaterial(i, 1) <<
333 << GetSandiaCofForMaterial(i, 3) <<
334 } <<
335 } <<
336 } 258 }
337 259
338 ////////////////////////////////////////////// << 260 ////////////////////////////////////////////////////////////////////////////
>> 261 //
>> 262 // SandiaIntervals
339 // 263 //
340 // Sandia matrix for PAI models based on vecto <<
341 264
342 void G4SandiaTable::ComputeMatSandiaMatrixPAI( << 265 G4int
>> 266 G4SandiaTable::SandiaIntervals(G4int Z[],
>> 267 G4int el )
343 { 268 {
344 G4int MaxIntervals = 0; << 269 G4int c,i ;
345 G4int elm, c, i, j, jj, k, k1, k2, c1, n1, z <<
346 <<
347 const auto noElm = (G4int)fMaterial->GetNumb <<
348 const G4ElementVector* ElementVector = fMate <<
349 270
350 std::vector<G4int> Z(noElm); // Atomic numb << 271 fMaxInterval = 0 ;
351 272
352 for (elm = 0; elm < noElm; ++elm) { << 273 for(i=0;i<el;i++)
353 z = G4lrint((*ElementVector)[elm]->GetZ()) << 274 {
354 if (z < 1) { << 275 fMaxInterval += fNbOfIntervals[Z[i]] ;
355 z = 1; <<
356 } <<
357 else if (z > 100) { <<
358 z = 100; <<
359 } <<
360 Z[elm] = z; <<
361 MaxIntervals += fNbOfIntervals[Z[elm]]; <<
362 } 276 }
363 fMaxInterval = MaxIntervals + 2; << 277 fMaxInterval += 2 ;
364 278
365 if (fVerbose > 0) { << 279 // G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl ;
366 G4cout << "G4SandiaTable::ComputeMatSandia <<
367 } <<
368 280
369 G4DataVector fPhotoAbsorptionCof0(fMaxInterv << 281 fPhotoAbsorptionCof = new G4double* [fMaxInterval] ;
370 G4DataVector fPhotoAbsorptionCof1(fMaxInterv <<
371 G4DataVector fPhotoAbsorptionCof2(fMaxInterv <<
372 G4DataVector fPhotoAbsorptionCof3(fMaxInterv <<
373 G4DataVector fPhotoAbsorptionCof4(fMaxInterv <<
374 282
375 for (c = 0; c < fMaxInterval; ++c) // just << 283 for(i = 0 ; i < fMaxInterval ; i++)
376 { 284 {
377 fPhotoAbsorptionCof0[c] = 0.; << 285 fPhotoAbsorptionCof[i] = new G4double[5] ;
378 fPhotoAbsorptionCof1[c] = 0.; <<
379 fPhotoAbsorptionCof2[c] = 0.; <<
380 fPhotoAbsorptionCof3[c] = 0.; <<
381 fPhotoAbsorptionCof4[c] = 0.; <<
382 } 286 }
383 c = 1; <<
384 <<
385 for (i = 0; i < noElm; ++i) { <<
386 G4double I1 = fIonizationPotentials[Z[i]] <<
387 n1 = 1; <<
388 287
389 for (j = 1; j < Z[i]; ++j) { <<
390 n1 += fNbOfIntervals[j]; <<
391 } <<
392 288
393 G4int n2 = n1 + fNbOfIntervals[Z[i]]; << 289 // for(c = 0 ; c < fIntervalLimit ; c++) // just in case
394 290
395 for (k1 = n1; k1 < n2; ++k1) { << 291 for(c = 0 ; c < fMaxInterval ; c++) // just in case
396 if (I1 > fSandiaTable[k1][0]) { << 292 {
397 continue; // no ionization for energi << 293 fPhotoAbsorptionCof[c][0] = 0. ;
398 } // ionisation potential) << 294 }
399 break; << 295 c = 1 ;
>> 296 for(i = 0 ; i < el ; i++)
>> 297 {
>> 298 G4double I1 = fIonizationPotentials[Z[i]]*keV ; // First ionization
>> 299 G4int n1 = 1 ; // potential in keV
>> 300 G4int j, c1, k1, k2 ;
>> 301 for(j = 1 ; j < Z[i] ; j++)
>> 302 {
>> 303 n1 += fNbOfIntervals[j] ;
400 } 304 }
401 G4int flag = 0; << 305 G4int n2 = n1 + fNbOfIntervals[Z[i]] ;
402 << 306
403 for (c1 = 1; c1 < c; ++c1) { << 307 for(k1 = n1 ; k1 < n2 ; k1++)
404 if (fPhotoAbsorptionCof0[c1] == I1) // << 308 {
>> 309 if(I1 > fSandiaTable[k1][0])
>> 310 {
>> 311 continue ; // no ionization for energies smaller than I1 (first
>> 312 } // ionisation potential)
>> 313 break ;
>> 314 }
>> 315 G4int flag = 0 ;
>> 316
>> 317 for(c1 = 1 ; c1 < c ; c1++)
>> 318 {
>> 319 if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
405 { 320 {
406 flag = 1; << 321 flag = 1 ;
407 break; << 322 break ;
408 } 323 }
409 } 324 }
410 if (flag == 0) { << 325 if(flag == 0)
411 fPhotoAbsorptionCof0[c] = I1; << 326 {
412 ++c; << 327 fPhotoAbsorptionCof[c][0] = I1 ;
>> 328 c++ ;
413 } 329 }
414 for (k2 = k1; k2 < n2; ++k2) { << 330 for(k2 = k1 ; k2 < n2 ; k2++)
415 flag = 0; << 331 {
416 << 332 flag = 0 ;
417 for (c1 = 1; c1 < c; ++c1) { << 333 for(c1 = 1 ; c1 < c ; c1++)
418 if (fPhotoAbsorptionCof0[c1] == fSandi << 334 {
419 flag = 1; << 335 if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
420 break; << 336 {
>> 337 flag = 1 ;
>> 338 break ;
421 } 339 }
422 } 340 }
423 if (flag == 0) { << 341 if(flag == 0)
424 fPhotoAbsorptionCof0[c] = fSandiaTable << 342 {
425 ++c; << 343 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0] ;
426 } << 344 c++ ;
427 } <<
428 } // end for(i) <<
429 // sort out <<
430 <<
431 for (i = 1; i < c; ++i) { <<
432 for (j = i + 1; j < c; ++j) { <<
433 if (fPhotoAbsorptionCof0[i] > fPhotoAbso <<
434 G4double tmp = fPhotoAbsorptionCof0[i] <<
435 fPhotoAbsorptionCof0[i] = fPhotoAbsorp <<
436 fPhotoAbsorptionCof0[j] = tmp; <<
437 } 345 }
438 } << 346 }
439 if (fVerbose > 0) { << 347 } // end for(i)
440 G4cout << i << "\t energy = " << fPhotoA << 348
441 } << 349 SandiaSort(fPhotoAbsorptionCof,c) ;
442 } << 350 fMaxInterval = c ;
443 fMaxInterval = c; << 351 return c ;
444 << 352 }
445 const G4double* fractionW = fMaterial->GetFr <<
446 353
447 if (fVerbose > 0) { << 354 ///////////////////////////////////////////////////////////////////////
448 for (i = 0; i < noElm; ++i) { << 355 //
449 G4cout << i << " = elN, fraction = " << << 356 // SandiaMixing
450 } << 357 //
451 } <<
452 <<
453 for (i = 0; i < noElm; ++i) { <<
454 n1 = 1; <<
455 G4double I1 = fIonizationPotentials[Z[i]] <<
456 <<
457 for (j = 1; j < Z[i]; ++j) { <<
458 n1 += fNbOfIntervals[j]; <<
459 } <<
460 <<
461 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1; <<
462 358
463 for (k = n1; k < n2; ++k) { << 359 G4int
464 G4double B1 = fSandiaTable[k][0]; << 360 G4SandiaTable::SandiaMixing( G4int Z[],
465 G4double B2 = fSandiaTable[k + 1][0]; << 361 const G4double fractionW[],
466 << 362 G4int el,
467 for (G4int q = 1; q < fMaxInterval - 1; << 363 G4int mi )
468 G4double E1 = fPhotoAbsorptionCof0[q]; << 364 {
469 G4double E2 = fPhotoAbsorptionCof0[q + << 365 G4int i;
470 << 366
471 if (fVerbose > 0) { << 367 for(i = 0 ; i < mi ; i++)
472 G4cout << "k = " << k << ", q = " << << 368 {
473 << ", E1 = " << E1 << ", E2 = << 369 for(G4int j = 1 ; j < 5 ; j++) fPhotoAbsorptionCof[i][j] = 0. ;
474 } << 370 }
475 if (B1 > E1 || B2 < E2 || E1 < I1) { << 371 for(i = 0 ; i < el ; i++)
476 if (fVerbose > 0) { << 372 {
477 G4cout << "continue for: B1 = " << << 373 G4int n1 = 1 ;
478 << ", E2 = " << E2 << G4end << 374 G4int j, k ;
479 } << 375 G4double I1 = fIonizationPotentials[Z[i]]*keV ;
480 continue; << 376 for(j = 1 ; j < Z[i] ; j++)
481 } << 377 {
482 fPhotoAbsorptionCof1[q] += fSandiaTabl << 378 n1 += fNbOfIntervals[j] ;
483 fPhotoAbsorptionCof2[q] += fSandiaTabl <<
484 fPhotoAbsorptionCof3[q] += fSandiaTabl <<
485 fPhotoAbsorptionCof4[q] += fSandiaTabl <<
486 } 379 }
487 } << 380 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1 ;
488 // Last interval <<
489 <<
490 fPhotoAbsorptionCof1[fMaxInterval - 1] += <<
491 fPhotoAbsorptionCof2[fMaxInterval - 1] += <<
492 fPhotoAbsorptionCof3[fMaxInterval - 1] += <<
493 fPhotoAbsorptionCof4[fMaxInterval - 1] += <<
494 } // for(i) <<
495 c = 0; // Deleting of first intervals where <<
496 <<
497 do { <<
498 ++c; <<
499 <<
500 if (fPhotoAbsorptionCof1[c] != 0.0 || fPho <<
501 fPhotoAbsorptionCof3[c] != 0.0 || fPho <<
502 { <<
503 continue; <<
504 } <<
505 <<
506 if (fVerbose > 0) { <<
507 G4cout << c << " = number with zero cofs <<
508 } <<
509 for (jj = 2; jj < fMaxInterval; ++jj) { <<
510 fPhotoAbsorptionCof0[jj - 1] = fPhotoAbs <<
511 fPhotoAbsorptionCof1[jj - 1] = fPhotoAbs <<
512 fPhotoAbsorptionCof2[jj - 1] = fPhotoAbs <<
513 fPhotoAbsorptionCof3[jj - 1] = fPhotoAbs <<
514 fPhotoAbsorptionCof4[jj - 1] = fPhotoAbs <<
515 } <<
516 --fMaxInterval; <<
517 --c; <<
518 } <<
519 // Loop checking, 07-Aug-2015, Vladimir Ivan <<
520 while (c < fMaxInterval - 1); <<
521 <<
522 if (fPhotoAbsorptionCof0[fMaxInterval - 1] = <<
523 fMaxInterval--; <<
524 } <<
525 381
526 // create the sandia matrix for this materia << 382 for(k = n1 ; k < n2 ; k++)
527 <<
528 fMatSandiaMatrixPAI = new G4OrderedTable(); <<
529 <<
530 G4double density = fMaterial->GetDensity(); <<
531 <<
532 for (i = 0; i < fMaxInterval; ++i) // -> G4 <<
533 { <<
534 fPhotoAbsorptionCof0[i + 1] *= funitc[0]; <<
535 fPhotoAbsorptionCof1[i + 1] *= funitc[1] * <<
536 fPhotoAbsorptionCof2[i + 1] *= funitc[2] * <<
537 fPhotoAbsorptionCof3[i + 1] *= funitc[3] * <<
538 fPhotoAbsorptionCof4[i + 1] *= funitc[4] * <<
539 } <<
540 if (fLowerI1) { <<
541 if (fMaterial->GetName() == "G4_WATER") { <<
542 fMaxInterval += fH2OlowerInt; <<
543 <<
544 for (i = 0; i < fMaxInterval; ++i) // i <<
545 { 383 {
546 fMatSandiaMatrixPAI->push_back(new G4D << 384 G4double B1 = fSandiaTable[k][0] ;
547 } << 385 G4double B2 = fSandiaTable[k+1][0] ;
548 for (i = 0; i < fH2OlowerInt; ++i) { << 386 for(G4int c = 1 ; c < mi-1 ; c++)
549 (*(*fMatSandiaMatrixPAI)[i])[0] = fH2O << 387 {
550 (*(*fMatSandiaMatrixPAI)[i])[1] = fH2O << 388 G4double E1 = fPhotoAbsorptionCof[c][0] ;
551 (*(*fMatSandiaMatrixPAI)[i])[2] = fH2O << 389 G4double E2 = fPhotoAbsorptionCof[c+1][0] ;
552 (*(*fMatSandiaMatrixPAI)[i])[3] = fH2O << 390 if(B1 > E1 || B2 < E2 || E1 < I1)
553 (*(*fMatSandiaMatrixPAI)[i])[4] = fH2O << 391 {
554 } << 392 continue ;
555 for (i = fH2OlowerInt; i < fMaxInterval; << 393 }
556 (*(*fMatSandiaMatrixPAI)[i])[0] = fPho << 394 for(j = 1 ; j < 5 ; j++)
557 (*(*fMatSandiaMatrixPAI)[i])[1] = fPho << 395 {
558 (*(*fMatSandiaMatrixPAI)[i])[2] = fPho << 396 fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i] ;
559 (*(*fMatSandiaMatrixPAI)[i])[3] = fPho << 397 }
560 (*(*fMatSandiaMatrixPAI)[i])[4] = fPho << 398 }
561 } << 399 }
562 } << 400 for(j = 1 ; j < 5 ; j++) // Last interval
563 } << 401 {
564 else { << 402 fPhotoAbsorptionCof[mi-1][j] += fSandiaTable[k][j]*fractionW[i] ;
565 for (i = 0; i < fMaxInterval; ++i) // ini << 403 }
>> 404 } // for(i)
>> 405 G4int c = 0 ; // Deleting of first intervals where all coefficients = 0
>> 406 do
566 { 407 {
567 fMatSandiaMatrixPAI->push_back(new G4Dat << 408 c++ ;
568 } << 409 if( fPhotoAbsorptionCof[c][1] != 0.0 ||
569 for (i = 0; i < fMaxInterval; ++i) { << 410 fPhotoAbsorptionCof[c][2] != 0.0 ||
570 (*(*fMatSandiaMatrixPAI)[i])[0] = fPhoto << 411 fPhotoAbsorptionCof[c][3] != 0.0 ||
571 (*(*fMatSandiaMatrixPAI)[i])[1] = fPhoto << 412 fPhotoAbsorptionCof[c][4] != 0.0 )
572 (*(*fMatSandiaMatrixPAI)[i])[2] = fPhoto << 413 {
573 (*(*fMatSandiaMatrixPAI)[i])[3] = fPhoto << 414 continue ;
574 (*(*fMatSandiaMatrixPAI)[i])[4] = fPhoto << 415 }
575 } << 416 for(G4int jj = 2 ; jj < mi ; jj++)
576 } << 417 {
577 // --fMaxInterval; << 418 for(G4int kk = 0 ; kk < 5 ; kk++)
578 // to avoid duplicate at 500 keV or extra ze << 419 {
579 << 420 fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk] ;
580 if (fVerbose > 0) { << 421 }
581 G4cout << "G4SandiaTable::ComputeMatSandia << 422 }
582 << G4endl; << 423 mi-- ;
583 << 424 c-- ;
584 for (i = 0; i < fMaxInterval; ++i) { << 425 }
585 G4cout << i << "\t" << GetSandiaMatTable << 426 while(c < mi - 1) ;
586 << GetSandiaMatTablePAI(i, 1) << << 427
587 << GetSandiaMatTablePAI(i, 3) << << 428 return mi ;
588 } << 429 }
589 } <<
590 return; <<
591 } <<
592 <<
593 ////////////////////////////////////////////// <<
594 // Methods for PAI model only <<
595 // <<
596 <<
597 G4SandiaTable::G4SandiaTable(G4int matIndex) <<
598 { <<
599 fMaterial = nullptr; <<
600 fMatNbOfIntervals = 0; <<
601 fMatSandiaMatrix = nullptr; <<
602 fMatSandiaMatrixPAI = nullptr; <<
603 fPhotoAbsorptionCof = nullptr; <<
604 <<
605 fMaxInterval = 0; <<
606 fVerbose = 0; <<
607 fLowerI1 = false; <<
608 <<
609 fSandiaCofPerAtom.resize(4, 0.0); <<
610 <<
611 const G4MaterialTable* theMaterialTable = G4 <<
612 auto numberOfMat = (G4int)G4Material::GetNum <<
613 <<
614 if (matIndex >= 0 && matIndex < numberOfMat) <<
615 fMaterial = (*theMaterialTable)[matIndex]; <<
616 } <<
617 else { <<
618 G4Exception( <<
619 "G4SandiaTable::G4SandiaTable(G4int matI <<
620 } <<
621 } <<
622 <<
623 ////////////////////////////////////////////// <<
624 <<
625 G4SandiaTable::G4SandiaTable() <<
626 { <<
627 fMaterial = nullptr; <<
628 fMatNbOfIntervals = 0; <<
629 fMatSandiaMatrix = nullptr; <<
630 fMatSandiaMatrixPAI = nullptr; <<
631 fPhotoAbsorptionCof = nullptr; <<
632 <<
633 fMaxInterval = 0; <<
634 fVerbose = 0; <<
635 fLowerI1 = false; <<
636 430
637 fSandiaCofPerAtom.resize(4, 0.0); <<
638 } <<
639 <<
640 ////////////////////////////////////////////// <<
641 <<
642 void G4SandiaTable::Initialize(const G4Materia <<
643 { <<
644 fMaterial = mat; <<
645 ComputeMatSandiaMatrixPAI(); <<
646 } <<
647 <<
648 ////////////////////////////////////////////// <<
649 <<
650 G4int G4SandiaTable::GetMaxInterval() const { <<
651 431
652 ////////////////////////////////////////////// <<
653 432
654 G4double** G4SandiaTable::GetPointerToCof() << 433 ////////////////////////////////////////////////////////////////////////////
655 { <<
656 if (fPhotoAbsorptionCof == nullptr) { <<
657 ComputeMatTable(); <<
658 } <<
659 return fPhotoAbsorptionCof; <<
660 } <<
661 <<
662 ////////////////////////////////////////////// <<
663 <<
664 void G4SandiaTable::SandiaSwap(G4double** da, <<
665 { <<
666 G4double tmp = da[i][0]; <<
667 da[i][0] = da[j][0]; <<
668 da[j][0] = tmp; <<
669 } <<
670 <<
671 ////////////////////////////////////////////// <<
672 <<
673 G4double G4SandiaTable::GetPhotoAbsorpCof(G4in <<
674 { <<
675 return fPhotoAbsorptionCof[i][j] * funitc[j] <<
676 } <<
677 <<
678 ////////////////////////////////////////////// <<
679 // 434 //
680 // Bubble sorting of left energy interval in S << 435 // Sandia interval and mixing calculations for materialCutsCouple constructor
681 // 436 //
682 437
683 void G4SandiaTable::SandiaSort(G4double** da, << 438 void G4SandiaTable::ComputeMatTable()
684 { 439 {
685 for (G4int i = 1; i < sz; ++i) { << 440 G4int MaxIntervals = 0;
686 for (G4int j = i + 1; j < sz; ++j) { << 441 G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n1;
687 if (da[i][0] > da[j][0]) { <<
688 SandiaSwap(da, i, j); <<
689 } <<
690 } <<
691 } <<
692 } <<
693 442
694 ////////////////////////////////////////////// << 443 const G4int noElm = fMaterial->GetNumberOfElements();
695 // << 444 const G4ElementVector* ElementVector = fMaterial->GetElementVector();
696 // SandiaIntervals << 445 G4int* Z = new G4int[noElm]; //Atomic number
697 // << 446
>> 447 for (elm = 0; elm<noElm; elm++)
>> 448 {
>> 449 Z[elm] = (G4int)(*ElementVector)[elm]->GetZ();
>> 450 MaxIntervals += fNbOfIntervals[Z[elm]];
>> 451 }
>> 452 fMaxInterval = 0 ;
698 453
699 G4int G4SandiaTable::SandiaIntervals(G4int Z[] << 454 for(i = 0; i < noElm; i++) fMaxInterval += fNbOfIntervals[Z[i]] ;
700 { << 455
701 G4int c, i, flag = 0, n1 = 1; << 456 fMaxInterval += 2 ;
702 G4int j, c1, k1, k2; <<
703 G4double I1; <<
704 fMaxInterval = 0; <<
705 457
706 for (i = 0; i < el; ++i) { << 458 // G4cout<<"fMaxInterval = "<<fMaxInterval<<G4endl ;
707 fMaxInterval += fNbOfIntervals[Z[i]]; <<
708 } <<
709 459
710 fMaxInterval += 2; << 460 fPhotoAbsorptionCof = new G4double* [fMaxInterval] ;
711 461
712 if (fVerbose > 0) { << 462 for(i = 0 ; i < fMaxInterval ; i++)
713 G4cout << "begin sanInt, fMaxInterval = " << 463 {
>> 464 fPhotoAbsorptionCof[i] = new G4double[5] ;
714 } 465 }
715 466
716 fPhotoAbsorptionCof = new G4double*[fMaxInte <<
717 467
718 for (i = 0; i < fMaxInterval; ++i) { << 468 // for(c = 0 ; c < fIntervalLimit ; c++) // just in case
719 fPhotoAbsorptionCof[i] = new G4double[5]; <<
720 } <<
721 // for(c = 0; c < fIntervalLimit; ++c) // <<
722 469
723 for (c = 0; c < fMaxInterval; ++c) { << 470 for(c = 0 ; c < fMaxInterval ; c++) // just in case
724 fPhotoAbsorptionCof[c][0] = 0.; << 471 {
>> 472 fPhotoAbsorptionCof[c][0] = 0. ;
725 } 473 }
>> 474 c = 1 ;
726 475
727 c = 1; << 476 for(i = 0 ; i < noElm ; i++)
728 << 477 {
729 for (i = 0; i < el; ++i) { << 478 G4double I1 = fIonizationPotentials[Z[i]]*keV ; // First ionization
730 I1 = fIonizationPotentials[Z[i]] * keV; / << 479 n1 = 1 ; // potential in keV
731 n1 = 1; // potential in keV << 480
732 << 481 for(j = 1 ; j < Z[i] ; j++)
733 for (j = 1; j < Z[i]; ++j) { << 482 {
734 n1 += fNbOfIntervals[j]; << 483 n1 += fNbOfIntervals[j] ;
735 } <<
736 <<
737 G4int n2 = n1 + fNbOfIntervals[Z[i]]; <<
738 <<
739 for (k1 = n1; k1 < n2; k1++) { <<
740 if (I1 > fSandiaTable[k1][0]) { <<
741 continue; // no ionization for energi <<
742 } // ionisation potential) <<
743 break; <<
744 } 484 }
745 flag = 0; << 485 G4int n2 = n1 + fNbOfIntervals[Z[i]] ;
746 << 486
747 for (c1 = 1; c1 < c; c1++) { << 487 for(k1 = n1 ; k1 < n2 ; k1++)
748 if (fPhotoAbsorptionCof[c1][0] == I1) / << 488 {
>> 489 if(I1 > fSandiaTable[k1][0])
>> 490 {
>> 491 continue ; // no ionization for energies smaller than I1 (first
>> 492 } // ionisation potential)
>> 493 break ;
>> 494 }
>> 495 G4int flag = 0 ;
>> 496
>> 497 for(c1 = 1 ; c1 < c ; c1++)
>> 498 {
>> 499 if(fPhotoAbsorptionCof[c1][0] == I1) // this value already has existed
749 { 500 {
750 flag = 1; << 501 flag = 1 ;
751 break; << 502 break ;
752 } 503 }
753 } 504 }
754 if (flag == 0) { << 505 if(flag == 0)
755 fPhotoAbsorptionCof[c][0] = I1; << 506 {
756 ++c; << 507 fPhotoAbsorptionCof[c][0] = I1 ;
>> 508 c++ ;
757 } 509 }
758 for (k2 = k1; k2 < n2; k2++) { << 510 for(k2 = k1 ; k2 < n2 ; k2++)
759 flag = 0; << 511 {
>> 512 flag = 0 ;
760 513
761 for (c1 = 1; c1 < c; c1++) { << 514 for(c1 = 1 ; c1 < c ; c1++)
762 if (fPhotoAbsorptionCof[c1][0] == fSan << 515 {
763 flag = 1; << 516 if(fPhotoAbsorptionCof[c1][0] == fSandiaTable[k2][0])
764 break; << 517 {
>> 518 flag = 1 ;
>> 519 break ;
765 } 520 }
766 } 521 }
767 if (flag == 0) { << 522 if(flag == 0)
768 fPhotoAbsorptionCof[c][0] = fSandiaTab << 523 {
769 if (fVerbose > 0) { << 524 fPhotoAbsorptionCof[c][0] = fSandiaTable[k2][0] ;
770 G4cout << "sanInt, c = " << c << ", << 525 c++ ;
771 } <<
772 ++c; <<
773 } 526 }
774 } << 527 }
775 } // end for(i) << 528 } // end for(i)
>> 529
>> 530 SandiaSort(fPhotoAbsorptionCof,c) ;
>> 531 fMaxInterval = c ;
>> 532
776 533
777 SandiaSort(fPhotoAbsorptionCof, c); << 534 // G4int
778 fMaxInterval = c; << 535 // G4SandiaTable::SandiaMixing( G4int Z[],
779 if (fVerbose > 0) { <<
780 G4cout << "end SanInt, fMaxInterval = " << <<
781 } <<
782 return c; <<
783 } <<
784 536
785 ////////////////////////////////////////////// << 537 const G4double* fractionW = fMaterial->GetFractionVector();
786 // <<
787 // SandiaMixing <<
788 // <<
789 538
790 G4int G4SandiaTable::SandiaMixing(G4int Z[], c << 539 // G4int el,
791 { << 540 // G4int mi )
792 G4int i, j, n1, k, c = 1, jj, kk; <<
793 G4double I1, B1, B2, E1, E2; <<
794 541
795 for (i = 0; i < mi; ++i) { << 542
796 for (j = 1; j < 5; ++j) { << 543
797 fPhotoAbsorptionCof[i][j] = 0.; << 544 for(i = 0 ; i < fMaxInterval ; i++)
798 } << 545 {
>> 546 for(j = 1 ; j < 5 ; j++) fPhotoAbsorptionCof[i][j] = 0.;
799 } 547 }
800 for (i = 0; i < el; ++i) { << 548 for(i = 0 ; i < noElm; i++)
801 n1 = 1; << 549 {
802 I1 = fIonizationPotentials[Z[i]] * keV; << 550 n1 = 1 ;
>> 551 G4double I1 = fIonizationPotentials[Z[i]]*keV ;
803 552
804 for (j = 1; j < Z[i]; ++j) { << 553 for(j = 1 ; j < Z[i] ; j++)
805 n1 += fNbOfIntervals[j]; << 554 {
806 } << 555 n1 += fNbOfIntervals[j] ;
>> 556 }
>> 557 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1 ;
807 558
808 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1; << 559 for(k = n1 ; k < n2 ; k++)
>> 560 {
>> 561 G4double B1 = fSandiaTable[k][0] ;
>> 562 G4double B2 = fSandiaTable[k+1][0] ;
>> 563 for(G4int c = 1 ; c < fMaxInterval-1 ; c++)
>> 564 {
>> 565 G4double E1 = fPhotoAbsorptionCof[c][0] ;
>> 566 G4double E2 = fPhotoAbsorptionCof[c+1][0] ;
>> 567 if(B1 > E1 || B2 < E2 || E1 < I1)
>> 568 {
>> 569 continue ;
>> 570 }
>> 571 for(j = 1 ; j < 5 ; j++)
>> 572 {
>> 573 fPhotoAbsorptionCof[c][j] += fSandiaTable[k][j]*fractionW[i] ;
>> 574 }
>> 575 }
>> 576 }
>> 577 for(j = 1 ; j < 5 ; j++) // Last interval
>> 578 {
>> 579 fPhotoAbsorptionCof[fMaxInterval-1][j] += fSandiaTable[k][j]*fractionW[i] ;
>> 580 }
>> 581 } // for(i)
809 582
810 for (k = n1; k < n2; ++k) { << 583 c = 0 ; // Deleting of first intervals where all coefficients = 0
811 B1 = fSandiaTable[k][0]; <<
812 B2 = fSandiaTable[k + 1][0]; <<
813 <<
814 for (c = 1; c < mi - 1; ++c) { <<
815 E1 = fPhotoAbsorptionCof[c][0]; <<
816 E2 = fPhotoAbsorptionCof[c + 1][0]; <<
817 584
818 if (B1 > E1 || B2 < E2 || E1 < I1) { << 585 do
819 continue; << 586 {
820 } << 587 c++ ;
821 588
822 for (j = 1; j < 5; ++j) { << 589 if( fPhotoAbsorptionCof[c][1] != 0.0 ||
823 fPhotoAbsorptionCof[c][j] += fSandia << 590 fPhotoAbsorptionCof[c][2] != 0.0 ||
824 if (fVerbose > 0) { << 591 fPhotoAbsorptionCof[c][3] != 0.0 ||
825 G4cout << "c=" << c << "; j=" << j << 592 fPhotoAbsorptionCof[c][4] != 0.0 ) continue ;
826 << "; frW=" << fractionW[i] << 593
827 } << 594 for(jj = 2 ; jj < fMaxInterval ; jj++)
828 } <<
829 } <<
830 } <<
831 for (j = 1; j < 5; ++j) // Last interval <<
832 { 595 {
833 fPhotoAbsorptionCof[mi - 1][j] += fSandi << 596 for(kk = 0 ; kk < 5 ; kk++)
834 if (fVerbose > 0) { << 597 {
835 G4cout << "mi-1=" << mi - 1 << "; j=" << 598 fPhotoAbsorptionCof[jj-1][kk]= fPhotoAbsorptionCof[jj][kk] ;
836 << "; frW=" << fractionW[i] << <<
837 } 599 }
838 } 600 }
839 } // for(i) << 601 fMaxInterval-- ;
840 c = 0; // Deleting of first intervals where << 602 c-- ;
841 << 603 }
842 do { << 604 while(c < fMaxInterval - 1) ;
843 ++c; << 605
844 << 606 // create the sandia matrix for this material
845 if (fPhotoAbsorptionCof[c][1] != 0.0 || fP << 607
846 fPhotoAbsorptionCof[c][3] != 0.0 || fP << 608 fMatSandiaMatrix = new G4OrderedTable();
>> 609
>> 610 for (i = 0; i < fMaxInterval; i++)
>> 611 {
>> 612 fMatSandiaMatrix->push_back(new G4DataVector(5,0.));
>> 613 }
>> 614 for (i = 0; i < fMaxInterval; i++)
>> 615 {
>> 616 for(j = 0 ; j < 5 ; j++)
847 { 617 {
848 continue; << 618 (*(*fMatSandiaMatrix)[i])[j] = fPhotoAbsorptionCof[i+1][j];
849 } << 619 }
>> 620 }
>> 621 delete [] Z;
>> 622 return ;
>> 623 }
850 624
851 for (jj = 2; jj < mi; ++jj) { <<
852 for (kk = 0; kk < 5; ++kk) { <<
853 fPhotoAbsorptionCof[jj - 1][kk] = fPho <<
854 } <<
855 } <<
856 mi--; <<
857 c--; <<
858 } <<
859 // Loop checking, 07-Aug-2015, Vladimir Ivan <<
860 while (c < mi - 1); <<
861 625
862 if (fVerbose > 0) { <<
863 G4cout << "end SanMix, mi = " << mi << G4e <<
864 } <<
865 626
866 return mi; <<
867 } <<
868 627
869 ////////////////////////////////////////////// <<
870 628
871 G4int G4SandiaTable::GetMatNbOfIntervals() con <<
872 629
873 ////////////////////////////////////////////// << 630 ///////////////////////////////////////////////////////////////////////
>> 631 //
>> 632 // GetNbOfIntervals
>> 633 //
874 634
875 G4double G4SandiaTable::GetSandiaPerAtom(G4int << 635 G4int
>> 636 G4SandiaTable::GetNbOfIntervals(G4int Z)
876 { 637 {
877 #ifdef G4VERBOSE << 638 return fNbOfIntervals[Z];
878 if (Z < 1 || Z > 100) { <<
879 Z = PrintErrorZ(Z, "GetSandiaPerAtom"); <<
880 } <<
881 if (interval < 0 || interval >= fNbOfInterva <<
882 PrintErrorV("GetSandiaPerAtom"); <<
883 interval = (interval < 0) ? 0 : fNbOfInter <<
884 } <<
885 if (j < 0 || j > 4) { <<
886 PrintErrorV("GetSandiaPerAtom"); <<
887 j = (j < 0) ? 0 : 4; <<
888 } <<
889 #endif <<
890 G4int row = fCumulInterval[Z - 1] + interval <<
891 G4double x = fSandiaTable[row][0] * CLHEP::k <<
892 if (j > 0) { <<
893 x = Z * CLHEP::amu / fZtoAratio[Z] * fSand <<
894 } <<
895 return x; <<
896 } 639 }
897 640
898 ////////////////////////////////////////////// << 641 ///////////////////////////////////////////////////////////////////////
>> 642 //
>> 643 // GetSandiaCofPerAtom(Z, interval, index)
>> 644 //
899 645
900 G4double G4SandiaTable::GetSandiaCofForMateria << 646 G4double
>> 647 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4int interval, G4int j)
901 { 648 {
902 #ifdef G4VERBOSE << 649 assert (Z>0 && Z<101 && interval>=0 && interval<fNbOfIntervals[Z]
903 if (interval < 0 || interval >= fMatNbOfInte << 650 && j>=0 && j<5);
904 PrintErrorV("GetSandiaCofForMaterial"); <<
905 interval = (interval < 0) ? 0 : fMatNbOfIn <<
906 } <<
907 if (j < 0 || j > 4) { <<
908 PrintErrorV("GetSandiaCofForMaterial"); <<
909 j = (j < 0) ? 0 : 4; <<
910 } <<
911 #endif <<
912 return ((*(*fMatSandiaMatrix)[interval])[j]) <<
913 } <<
914 <<
915 ////////////////////////////////////////////// <<
916 651
917 const G4double* G4SandiaTable::GetSandiaCofFor << 652 G4int row = fCumulInterval[Z-1] + interval;
918 { << 653 if (j==0) return fSandiaTable[row][0]*keV;
919 G4int interval = 0; << 654 G4double AoverAvo = Z*amu/fZtoAratio[Z];
920 if (energy > (*(*fMatSandiaMatrix)[0])[0]) { << 655 return AoverAvo*(fSandiaTable[row][j]*cm2*std::pow(keV,G4double(j))/g);
921 interval = fMatNbOfIntervals - 1; <<
922 // Loop checking, 07-Aug-2015, Vladimir Iv <<
923 while ((interval > 0) && (energy < (*(*fMa <<
924 --interval; <<
925 } <<
926 } <<
927 return &((*(*fMatSandiaMatrix)[interval])[1] <<
928 } 656 }
929 657
930 ////////////////////////////////////////////// << 658 ///////////////////////////////////////////////////////////////////////
>> 659 //
>> 660 // GetSandiaCofPerAtom(Z, energy)
>> 661 //
931 662
932 G4double G4SandiaTable::GetSandiaMatTable(G4in << 663 G4double*
>> 664 G4SandiaTable::GetSandiaCofPerAtom(G4int Z, G4double energy)
933 { 665 {
934 #ifdef G4VERBOSE << 666 assert (Z > 0);
935 if (interval < 0 || interval >= fMatNbOfInte << 667
936 PrintErrorV("GetSandiaCofForMaterial"); << 668 G4double Emin = fSandiaTable[fCumulInterval[Z-1]][0]*keV;
937 interval = (interval < 0) ? 0 : fMatNbOfIn << 669 G4double Iopot = fIonizationPotentials[Z]*eV;
938 } << 670 if (Iopot > Emin) Emin = Iopot;
939 if (j < 0 || j > 4) { << 671
940 PrintErrorV("GetSandiaCofForMaterial"); << 672 G4int interval = fNbOfIntervals[Z] - 1;
941 j = (j < 0) ? 0 : 4; << 673 G4int row = fCumulInterval[Z-1] + interval;
942 } << 674 while ((interval>0) && (energy<fSandiaTable[row][0]*keV))
943 #endif << 675 row = fCumulInterval[Z-1] + --interval;
944 return ((*(*fMatSandiaMatrix)[interval])[j]) << 676
>> 677 if (energy >= Emin)
>> 678 {
>> 679 G4double AoverAvo = Z*amu/fZtoAratio[Z];
>> 680
>> 681 fSandiaCofPerAtom[0]=AoverAvo*(fSandiaTable[row][1]*cm2*keV/g);
>> 682 fSandiaCofPerAtom[1]=AoverAvo*(fSandiaTable[row][2]*cm2*keV*keV/g);
>> 683 fSandiaCofPerAtom[2]=AoverAvo*(fSandiaTable[row][3]*cm2*keV*keV*keV/g);
>> 684 fSandiaCofPerAtom[3]=AoverAvo*(fSandiaTable[row][4]*cm2*keV*keV*keV*keV/g);
>> 685 }
>> 686 else
>> 687 fSandiaCofPerAtom[0] = fSandiaCofPerAtom[1] = fSandiaCofPerAtom[2] =
>> 688 fSandiaCofPerAtom[3] = 0.;
>> 689
>> 690 return fSandiaCofPerAtom;
>> 691
945 } 692 }
946 693
947 ////////////////////////////////////////////// << 694 ///////////////////////////////////////////////////////////////////////
>> 695 //
>> 696 // GetIonizationPot
>> 697 //
948 698
949 G4double G4SandiaTable::GetSandiaMatTablePAI(G << 699 G4double
>> 700 G4SandiaTable::GetIonizationPot(G4int Z)
950 { 701 {
951 #ifdef G4VERBOSE << 702 return fIonizationPotentials[Z]*eV;
952 if (interval < 0 || interval >= fMaxInterval <<
953 PrintErrorV("GetSandiaCofForMaterialPAI"); <<
954 interval = (interval < 0) ? 0 : fMaxInterv <<
955 } <<
956 if (j < 0 || j > 4) { <<
957 PrintErrorV("GetSandiaCofForMaterialPAI"); <<
958 j = (j < 0) ? 0 : 4; <<
959 } <<
960 #endif <<
961 return ((*(*fMatSandiaMatrixPAI)[interval])[ <<
962 } 703 }
963 704
964 ////////////////////////////////////////////// << 705 ///////////////////////////////////////////////////////////////////////
965 // 706 //
966 // Sandia interval and mixing calculations fo << 707 // GetZtoA
967 // 708 //
968 709
969 void G4SandiaTable::ComputeMatTable() << 710 G4double
>> 711 G4SandiaTable::GetZtoA(G4int Z)
970 { 712 {
971 G4int elm, c, i, j, jj, k, kk, k1, k2, c1, n << 713 return fZtoAratio[Z];
972 << 714 }
973 const auto noElm = (G4int)fMaterial->GetNumb <<
974 const G4ElementVector* ElementVector = fMate <<
975 auto Z = new G4int[noElm]; // Atomic number <<
976 <<
977 fMaxInterval = 0; <<
978 for (elm = 0; elm < noElm; ++elm) { <<
979 Z[elm] = (*ElementVector)[elm]->GetZasInt( <<
980 fMaxInterval += fNbOfIntervals[Z[elm]]; <<
981 } <<
982 fMaxInterval += 2; <<
983 <<
984 // G4cout<<"fMaxInterval = "<<fMaxInterval< <<
985 <<
986 fPhotoAbsorptionCof = new G4double*[fMaxInte <<
987 <<
988 for (i = 0; i < fMaxInterval; ++i) { <<
989 fPhotoAbsorptionCof[i] = new G4double[5]; <<
990 } <<
991 <<
992 // for(c = 0; c < fIntervalLimit; ++c) // <<
993 <<
994 for (c = 0; c < fMaxInterval; ++c) // just <<
995 { <<
996 fPhotoAbsorptionCof[c][0] = 0.; <<
997 } <<
998 c = 1; <<
999 <<
1000 for (i = 0; i < noElm; ++i) { <<
1001 G4double I1 = fIonizationPotentials[Z[i]] <<
1002 n1 = 1; // potential in keV <<
1003 <<
1004 for (j = 1; j < Z[i]; ++j) { <<
1005 n1 += fNbOfIntervals[j]; <<
1006 } <<
1007 G4int n2 = n1 + fNbOfIntervals[Z[i]]; <<
1008 <<
1009 for (k1 = n1; k1 < n2; ++k1) { <<
1010 if (I1 > fSandiaTable[k1][0]) { <<
1011 continue; // no ionization for energ <<
1012 } // ionisation potential) <<
1013 break; <<
1014 } <<
1015 G4int flag = 0; <<
1016 <<
1017 for (c1 = 1; c1 < c; ++c1) { <<
1018 if (fPhotoAbsorptionCof[c1][0] == I1) <<
1019 { <<
1020 flag = 1; <<
1021 break; <<
1022 } <<
1023 } <<
1024 if (flag == 0) { <<
1025 fPhotoAbsorptionCof[c][0] = I1; <<
1026 ++c; <<
1027 } <<
1028 for (k2 = k1; k2 < n2; ++k2) { <<
1029 flag = 0; <<
1030 <<
1031 for (c1 = 1; c1 < c; ++c1) { <<
1032 if (fPhotoAbsorptionCof[c1][0] == fSa <<
1033 flag = 1; <<
1034 break; <<
1035 } <<
1036 } <<
1037 if (flag == 0) { <<
1038 fPhotoAbsorptionCof[c][0] = fSandiaTa <<
1039 ++c; <<
1040 } <<
1041 } <<
1042 } // end for(i) <<
1043 <<
1044 SandiaSort(fPhotoAbsorptionCof, c); <<
1045 fMaxInterval = c; <<
1046 <<
1047 const G4double* fractionW = fMaterial->GetF <<
1048 <<
1049 for (i = 0; i < fMaxInterval; ++i) { <<
1050 for (j = 1; j < 5; ++j) { <<
1051 fPhotoAbsorptionCof[i][j] = 0.; <<
1052 } <<
1053 } <<
1054 for (i = 0; i < noElm; ++i) { <<
1055 n1 = 1; <<
1056 G4double I1 = fIonizationPotentials[Z[i]] <<
1057 <<
1058 for (j = 1; j < Z[i]; ++j) { <<
1059 n1 += fNbOfIntervals[j]; <<
1060 } <<
1061 G4int n2 = n1 + fNbOfIntervals[Z[i]] - 1; <<
1062 <<
1063 for (k = n1; k < n2; ++k) { <<
1064 G4double B1 = fSandiaTable[k][0]; <<
1065 G4double B2 = fSandiaTable[k + 1][0]; <<
1066 for (G4int q = 1; q < fMaxInterval - 1; <<
1067 G4double E1 = fPhotoAbsorptionCof[q][ <<
1068 G4double E2 = fPhotoAbsorptionCof[q + <<
1069 if (B1 > E1 || B2 < E2 || E1 < I1) { <<
1070 continue; <<
1071 } <<
1072 for (j = 1; j < 5; ++j) { <<
1073 fPhotoAbsorptionCof[q][j] += fSandi <<
1074 } <<
1075 } <<
1076 } <<
1077 for (j = 1; j < 5; ++j) // Last interval <<
1078 { <<
1079 fPhotoAbsorptionCof[fMaxInterval - 1][j <<
1080 } <<
1081 } // for(i) <<
1082 <<
1083 c = 0; // Deleting of first intervals wher <<
1084 <<
1085 do { <<
1086 ++c; <<
1087 <<
1088 if (fPhotoAbsorptionCof[c][1] != 0.0 || f <<
1089 fPhotoAbsorptionCof[c][3] != 0.0 || f <<
1090 { <<
1091 continue; <<
1092 } <<
1093 <<
1094 for (jj = 2; jj < fMaxInterval; ++jj) { <<
1095 for (kk = 0; kk < 5; ++kk) { <<
1096 fPhotoAbsorptionCof[jj - 1][kk] = fPh <<
1097 } <<
1098 } <<
1099 --fMaxInterval; <<
1100 --c; <<
1101 } <<
1102 // Loop checking, 07-Aug-2015, Vladimir Iva <<
1103 while (c < fMaxInterval - 1); <<
1104 <<
1105 // create the sandia matrix for this materi <<
1106 <<
1107 --fMaxInterval; // vmg 20.11.10 <<
1108 <<
1109 fMatSandiaMatrix = new G4OrderedTable(); <<
1110 715
1111 for (i = 0; i < fMaxInterval; ++i) { << 716 // G4SandiaTable class -- end of implementation file
1112 fMatSandiaMatrix->push_back(new G4DataVec << 717 //
1113 } << 718 ////////////////////////////////////////////////////////////////////////////
1114 for (i = 0; i < fMaxInterval; ++i) { <<
1115 for (j = 0; j < 5; ++j) { <<
1116 (*(*fMatSandiaMatrix)[i])[j] = fPhotoAb <<
1117 } <<
1118 } <<
1119 fMatNbOfIntervals = fMaxInterval; <<
1120 719
1121 if (fVerbose > 0) { <<
1122 G4cout << "vmg, G4SandiaTable::ComputeMat <<
1123 720
1124 for (i = 0; i < fMaxInterval; ++i) { <<
1125 G4cout << i << "\t" << GetSandiaCofForM <<
1126 << this->GetSandiaCofForMaterial <<
1127 << "\t" << this->GetSandiaCofFor <<
1128 << this->GetSandiaCofForMaterial <<
1129 } <<
1130 } <<
1131 delete[] Z; <<
1132 return; <<
1133 } <<
1134 721