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