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