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