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