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