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 // 26 // 27 // Geant4 class G4EmUtility 27 // Geant4 class G4EmUtility 28 // 28 // 29 // Author V.Ivanchenko 14.03.2022 29 // Author V.Ivanchenko 14.03.2022 30 // 30 // 31 31 32 #include "G4EmUtility.hh" 32 #include "G4EmUtility.hh" 33 #include "G4RegionStore.hh" 33 #include "G4RegionStore.hh" 34 #include "G4ProductionCutsTable.hh" 34 #include "G4ProductionCutsTable.hh" 35 #include "G4VEmProcess.hh" 35 #include "G4VEmProcess.hh" 36 #include "G4EmParameters.hh" 36 #include "G4EmParameters.hh" 37 #include "G4PhysicsVector.hh" 37 #include "G4PhysicsVector.hh" 38 #include "Randomize.hh" 38 #include "Randomize.hh" 39 #include "G4Log.hh" 39 #include "G4Log.hh" 40 #include "G4Exp.hh" 40 #include "G4Exp.hh" 41 41 42 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 43 43 44 static const G4double g4log10 = G4Log(10.); 44 static const G4double g4log10 = G4Log(10.); 45 45 46 const G4Region* 46 const G4Region* 47 G4EmUtility::FindRegion(const G4String& region 47 G4EmUtility::FindRegion(const G4String& regionName, const G4int verbose) 48 { 48 { 49 const G4Region* reg = nullptr; 49 const G4Region* reg = nullptr; 50 G4RegionStore* regStore = G4RegionStore::Get 50 G4RegionStore* regStore = G4RegionStore::GetInstance(); 51 G4String r = regionName; 51 G4String r = regionName; 52 if(r == "") { r = "DefaultRegionForTheWorld" 52 if(r == "") { r = "DefaultRegionForTheWorld"; } 53 reg = regStore->GetRegion(r, true); 53 reg = regStore->GetRegion(r, true); 54 if(nullptr == reg && verbose > 0) { 54 if(nullptr == reg && verbose > 0) { 55 G4cout << "### G4EmUtility WARNING: fails 55 G4cout << "### G4EmUtility WARNING: fails to find a region <" 56 << r << G4endl; 56 << r << G4endl; 57 } else if(verbose > 1) { 57 } else if(verbose > 1) { 58 G4cout << "### G4EmUtility finds out G4Reg 58 G4cout << "### G4EmUtility finds out G4Region <" << r << ">" 59 << G4endl; 59 << G4endl; 60 } 60 } 61 return reg; 61 return reg; 62 } 62 } 63 63 64 //....oooOO0OOooo........oooOO0OOooo........oo 64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 65 65 66 const G4Element* G4EmUtility::SampleRandomElem 66 const G4Element* G4EmUtility::SampleRandomElement(const G4Material* mat) 67 { 67 { 68 const G4Element* elm = mat->GetElement(0); 68 const G4Element* elm = mat->GetElement(0); 69 std::size_t nElements = mat->GetNumberOfElem 69 std::size_t nElements = mat->GetNumberOfElements(); 70 if(1 < nElements) { 70 if(1 < nElements) { 71 G4double x = mat->GetTotNbOfElectPerVolume 71 G4double x = mat->GetTotNbOfElectPerVolume()*G4UniformRand(); 72 const G4double* y = mat->GetVecNbOfAtomsPe 72 const G4double* y = mat->GetVecNbOfAtomsPerVolume(); 73 for(std::size_t i=0; i<nElements; ++i) { 73 for(std::size_t i=0; i<nElements; ++i) { 74 elm = mat->GetElement((G4int)i); 74 elm = mat->GetElement((G4int)i); 75 x -= y[i]*elm->GetZ(); 75 x -= y[i]*elm->GetZ(); 76 if(x <= 0.0) { break; } 76 if(x <= 0.0) { break; } 77 } 77 } 78 } 78 } 79 return elm; 79 return elm; 80 } 80 } 81 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 83 83 84 const G4Isotope* G4EmUtility::SampleRandomIsot 84 const G4Isotope* G4EmUtility::SampleRandomIsotope(const G4Element* elm) 85 { 85 { 86 const std::size_t ni = elm->GetNumberOfIsoto 86 const std::size_t ni = elm->GetNumberOfIsotopes(); 87 const G4Isotope* iso = elm->GetIsotope(0); 87 const G4Isotope* iso = elm->GetIsotope(0); 88 if(ni > 1) { 88 if(ni > 1) { 89 const G4double* ab = elm->GetRelativeAbund 89 const G4double* ab = elm->GetRelativeAbundanceVector(); 90 G4double x = G4UniformRand(); 90 G4double x = G4UniformRand(); 91 for(std::size_t idx=0; idx<ni; ++idx) { 91 for(std::size_t idx=0; idx<ni; ++idx) { 92 x -= ab[idx]; 92 x -= ab[idx]; 93 if (x <= 0.0) { 93 if (x <= 0.0) { 94 iso = elm->GetIsotope((G4int)idx); 94 iso = elm->GetIsotope((G4int)idx); 95 break; 95 break; 96 } 96 } 97 } 97 } 98 } 98 } 99 return iso; 99 return iso; 100 } 100 } 101 101 102 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 103 103 104 std::vector<G4double>* G4EmUtility::FindCrossS 104 std::vector<G4double>* G4EmUtility::FindCrossSectionMax(G4PhysicsTable* p) 105 { 105 { 106 std::vector<G4double>* ptr = nullptr; 106 std::vector<G4double>* ptr = nullptr; 107 if(nullptr == p) { return ptr; } 107 if(nullptr == p) { return ptr; } 108 108 109 const std::size_t n = p->length(); 109 const std::size_t n = p->length(); 110 ptr = new std::vector<G4double>; 110 ptr = new std::vector<G4double>; 111 ptr->resize(n, DBL_MAX); 111 ptr->resize(n, DBL_MAX); 112 112 113 G4bool isPeak = false; 113 G4bool isPeak = false; 114 G4double e, ss, ee, xs; 114 G4double e, ss, ee, xs; 115 115 116 // first loop on existing vectors 116 // first loop on existing vectors 117 for (std::size_t i=0; i<n; ++i) { 117 for (std::size_t i=0; i<n; ++i) { 118 const G4PhysicsVector* pv = (*p)[i]; 118 const G4PhysicsVector* pv = (*p)[i]; 119 xs = ee = 0.0; 119 xs = ee = 0.0; 120 if(nullptr != pv) { 120 if(nullptr != pv) { 121 G4int nb = (G4int)pv->GetVectorLength(); 121 G4int nb = (G4int)pv->GetVectorLength(); 122 for (G4int j=0; j<nb; ++j) { 122 for (G4int j=0; j<nb; ++j) { 123 e = pv->Energy(j); 123 e = pv->Energy(j); 124 ss = (*pv)(j); 124 ss = (*pv)(j); 125 if(ss >= xs) { 125 if(ss >= xs) { 126 xs = ss; 126 xs = ss; 127 ee = e; 127 ee = e; 128 continue; 128 continue; 129 } else { 129 } else { 130 isPeak = true; 130 isPeak = true; 131 (*ptr)[i] = ee; 131 (*ptr)[i] = ee; 132 break; 132 break; 133 } 133 } 134 } 134 } 135 } 135 } 136 } 136 } 137 137 138 // there is no peak for any material 138 // there is no peak for any material 139 if(!isPeak) { 139 if(!isPeak) { 140 delete ptr; 140 delete ptr; 141 ptr = nullptr; 141 ptr = nullptr; 142 } 142 } 143 return ptr; 143 return ptr; 144 } 144 } 145 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 147 147 148 std::vector<G4double>* 148 std::vector<G4double>* 149 G4EmUtility::FindCrossSectionMax(G4VDiscretePr 149 G4EmUtility::FindCrossSectionMax(G4VDiscreteProcess* p, 150 const G4Parti 150 const G4ParticleDefinition* part) 151 { 151 { 152 std::vector<G4double>* ptr = nullptr; 152 std::vector<G4double>* ptr = nullptr; 153 if (nullptr == p || nullptr == part) { retur 153 if (nullptr == p || nullptr == part) { return ptr; } 154 154 155 G4EmParameters* theParameters = G4EmParamete 155 G4EmParameters* theParameters = G4EmParameters::Instance(); 156 const G4double tmin = theParameters->MinKinE 156 const G4double tmin = theParameters->MinKinEnergy(); 157 const G4double tmax = theParameters->MaxKinE 157 const G4double tmax = theParameters->MaxKinEnergy(); 158 const G4double ee = G4Log(tmax/tmin); 158 const G4double ee = G4Log(tmax/tmin); 159 const G4double scale = theParameters->Number 159 const G4double scale = theParameters->NumberOfBinsPerDecade()/g4log10; 160 G4int nbin = static_cast<G4int>(ee*scale); 160 G4int nbin = static_cast<G4int>(ee*scale); 161 nbin = std::max(nbin, 4); 161 nbin = std::max(nbin, 4); 162 G4double x = G4Exp(ee/(G4double)nbin); 162 G4double x = G4Exp(ee/(G4double)nbin); 163 163 164 const G4ProductionCutsTable* theCoupleTable= 164 const G4ProductionCutsTable* theCoupleTable= 165 G4ProductionCutsTable::GetProductionCu 165 G4ProductionCutsTable::GetProductionCutsTable(); 166 std::size_t n = theCoupleTable->GetTableSize 166 std::size_t n = theCoupleTable->GetTableSize(); 167 ptr = new std::vector<G4double>; 167 ptr = new std::vector<G4double>; 168 ptr->resize(n, DBL_MAX); 168 ptr->resize(n, DBL_MAX); 169 169 170 G4bool isPeak = false; 170 G4bool isPeak = false; 171 171 172 // first loop on existing vectors 172 // first loop on existing vectors 173 const G4int nn = static_cast<G4int>(n); 173 const G4int nn = static_cast<G4int>(n); 174 for (G4int i=0; i<nn; ++i) { 174 for (G4int i=0; i<nn; ++i) { 175 G4double sm = 0.0; 175 G4double sm = 0.0; 176 G4double em = 0.0; 176 G4double em = 0.0; 177 G4double e = tmin; 177 G4double e = tmin; 178 for (G4int j=0; j<=nbin; ++j) { 178 for (G4int j=0; j<=nbin; ++j) { 179 G4double sig = p->GetCrossSection(e, the 179 G4double sig = p->GetCrossSection(e, theCoupleTable->GetMaterialCutsCouple(i)); 180 if (sig >= sm) { 180 if (sig >= sm) { 181 em = e; 181 em = e; 182 sm = sig; 182 sm = sig; 183 e = (j+1 < nbin) ? e*x : tmax; 183 e = (j+1 < nbin) ? e*x : tmax; 184 } else { 184 } else { 185 isPeak = true; 185 isPeak = true; 186 (*ptr)[i] = em; 186 (*ptr)[i] = em; 187 break; 187 break; 188 } 188 } 189 } 189 } 190 } 190 } 191 // there is no peak for any couple 191 // there is no peak for any couple 192 if (!isPeak) { 192 if (!isPeak) { 193 delete ptr; 193 delete ptr; 194 ptr = nullptr; 194 ptr = nullptr; 195 } 195 } 196 return ptr; 196 return ptr; 197 } 197 } 198 198 199 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 200 200 201 std::vector<G4TwoPeaksXS*>* 201 std::vector<G4TwoPeaksXS*>* 202 G4EmUtility::FillPeaksStructure(G4PhysicsTable 202 G4EmUtility::FillPeaksStructure(G4PhysicsTable* p, G4LossTableBuilder* bld) 203 { 203 { 204 std::vector<G4TwoPeaksXS*>* ptr = nullptr; 204 std::vector<G4TwoPeaksXS*>* ptr = nullptr; 205 if(nullptr == p) { return ptr; } 205 if(nullptr == p) { return ptr; } 206 206 207 const G4int n = (G4int)p->length(); 207 const G4int n = (G4int)p->length(); 208 ptr = new std::vector<G4TwoPeaksXS*>; 208 ptr = new std::vector<G4TwoPeaksXS*>; 209 ptr->resize(n, nullptr); 209 ptr->resize(n, nullptr); 210 210 211 G4double e, ss, xs, ee; 211 G4double e, ss, xs, ee; 212 G4double e1peak, e1deep, e2peak, e2deep, e3p 212 G4double e1peak, e1deep, e2peak, e2deep, e3peak; 213 G4bool isDeep = false; 213 G4bool isDeep = false; 214 214 215 // first loop on existing vectors 215 // first loop on existing vectors 216 for (G4int i=0; i<n; ++i) { 216 for (G4int i=0; i<n; ++i) { 217 const G4PhysicsVector* pv = (*p)[i]; 217 const G4PhysicsVector* pv = (*p)[i]; 218 ee = xs = 0.0; 218 ee = xs = 0.0; 219 e1peak = e1deep = e2peak = e2deep = e3peak 219 e1peak = e1deep = e2peak = e2deep = e3peak = DBL_MAX; 220 if(nullptr != pv) { 220 if(nullptr != pv) { 221 G4int nb = (G4int)pv->GetVectorLength(); 221 G4int nb = (G4int)pv->GetVectorLength(); 222 for (G4int j=0; j<nb; ++j) { 222 for (G4int j=0; j<nb; ++j) { 223 e = pv->Energy(j); 223 e = pv->Energy(j); 224 ss = (*pv)(j); 224 ss = (*pv)(j); 225 // find out 1st peak 225 // find out 1st peak 226 if(e1peak == DBL_MAX) { 226 if(e1peak == DBL_MAX) { 227 if(ss >= xs) { 227 if(ss >= xs) { 228 xs = ss; 228 xs = ss; 229 ee = e; 229 ee = e; 230 continue; 230 continue; 231 } else { 231 } else { 232 e1peak = ee; 232 e1peak = ee; 233 } 233 } 234 } 234 } 235 // find out the deep 235 // find out the deep 236 if(e1deep == DBL_MAX) { 236 if(e1deep == DBL_MAX) { 237 if(ss <= xs) { 237 if(ss <= xs) { 238 xs = ss; 238 xs = ss; 239 ee = e; 239 ee = e; 240 continue; 240 continue; 241 } else { 241 } else { 242 e1deep = ee; 242 e1deep = ee; 243 isDeep = true; 243 isDeep = true; 244 } 244 } 245 } 245 } 246 // find out 2nd peak 246 // find out 2nd peak 247 if(e2peak == DBL_MAX) { 247 if(e2peak == DBL_MAX) { 248 if(ss >= xs) { 248 if(ss >= xs) { 249 xs = ss; 249 xs = ss; 250 ee = e; 250 ee = e; 251 continue; 251 continue; 252 } else { 252 } else { 253 e2peak = ee; 253 e2peak = ee; 254 } 254 } 255 } 255 } 256 if(e2deep == DBL_MAX) { 256 if(e2deep == DBL_MAX) { 257 if(ss <= xs) { 257 if(ss <= xs) { 258 xs = ss; 258 xs = ss; 259 ee = e; 259 ee = e; 260 continue; 260 continue; 261 } else { 261 } else { 262 e2deep = ee; 262 e2deep = ee; 263 break; 263 break; 264 } 264 } 265 } 265 } 266 // find out 3d peak 266 // find out 3d peak 267 if(e3peak == DBL_MAX) { 267 if(e3peak == DBL_MAX) { 268 if(ss >= xs) { 268 if(ss >= xs) { 269 xs = ss; 269 xs = ss; 270 ee = e; 270 ee = e; 271 continue; 271 continue; 272 } else { 272 } else { 273 e3peak = ee; 273 e3peak = ee; 274 } 274 } 275 } 275 } 276 } 276 } 277 } 277 } 278 G4TwoPeaksXS* x = (*ptr)[i]; 278 G4TwoPeaksXS* x = (*ptr)[i]; 279 if(nullptr == x) { 279 if(nullptr == x) { 280 x = new G4TwoPeaksXS(); 280 x = new G4TwoPeaksXS(); 281 (*ptr)[i] = x; 281 (*ptr)[i] = x; 282 } 282 } 283 x->e1peak = e1peak; 283 x->e1peak = e1peak; 284 x->e1deep = e1deep; 284 x->e1deep = e1deep; 285 x->e2peak = e2peak; 285 x->e2peak = e2peak; 286 x->e2deep = e2deep; 286 x->e2deep = e2deep; 287 x->e3peak = e3peak; 287 x->e3peak = e3peak; 288 } 288 } 289 // case of no 1st peak in all vectors 289 // case of no 1st peak in all vectors 290 if(!isDeep) { 290 if(!isDeep) { 291 delete ptr; 291 delete ptr; 292 ptr = nullptr; 292 ptr = nullptr; 293 return ptr; 293 return ptr; 294 } 294 } 295 // check base particles 295 // check base particles 296 if(!bld->GetBaseMaterialFlag()) { return ptr 296 if(!bld->GetBaseMaterialFlag()) { return ptr; } 297 297 298 auto theDensityIdx = bld->GetCoupleIndexes() 298 auto theDensityIdx = bld->GetCoupleIndexes(); 299 // second loop using base materials 299 // second loop using base materials 300 for (G4int i=0; i<n; ++i) { 300 for (G4int i=0; i<n; ++i) { 301 const G4PhysicsVector* pv = (*p)[i]; 301 const G4PhysicsVector* pv = (*p)[i]; 302 if (nullptr == pv) { 302 if (nullptr == pv) { 303 G4int j = (*theDensityIdx)[i]; 303 G4int j = (*theDensityIdx)[i]; 304 if(j == i) { continue; } 304 if(j == i) { continue; } 305 G4TwoPeaksXS* x = (*ptr)[i]; 305 G4TwoPeaksXS* x = (*ptr)[i]; 306 G4TwoPeaksXS* y = (*ptr)[j]; 306 G4TwoPeaksXS* y = (*ptr)[j]; 307 if(nullptr == x) { 307 if(nullptr == x) { 308 x = new G4TwoPeaksXS(); 308 x = new G4TwoPeaksXS(); 309 (*ptr)[i] = x; 309 (*ptr)[i] = x; 310 } 310 } 311 x->e1peak = y->e1peak; 311 x->e1peak = y->e1peak; 312 x->e1deep = y->e1deep; 312 x->e1deep = y->e1deep; 313 x->e2peak = y->e2peak; 313 x->e2peak = y->e2peak; 314 x->e2deep = y->e2deep; 314 x->e2deep = y->e2deep; 315 x->e3peak = y->e3peak; 315 x->e3peak = y->e3peak; 316 } 316 } 317 } 317 } 318 return ptr; 318 return ptr; 319 } 319 } 320 320 321 //....oooOO0OOooo........oooOO0OOooo........oo 321 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 322 322 323 void G4EmUtility::InitialiseElementSelectors(G 323 void G4EmUtility::InitialiseElementSelectors(G4VEmModel* mod, 324 const G4ParticleDefinition* par 324 const G4ParticleDefinition* part, 325 const G4DataVector& cuts, 325 const G4DataVector& cuts, 326 c 326 const G4double elow, 327 c 327 const G4double ehigh) 328 { 328 { 329 // using spline for element selectors should 329 // using spline for element selectors should be investigated in details 330 // because small number of points may provid 330 // because small number of points may provide biased results 331 // large number of points requires significa 331 // large number of points requires significant increase of memory 332 G4bool spline = false; 332 G4bool spline = false; 333 333 334 G4int nbinsPerDec = G4EmParameters::Instance 334 G4int nbinsPerDec = G4EmParameters::Instance()->NumberOfBinsPerDecade(); 335 335 336 G4ProductionCutsTable* theCoupleTable= 336 G4ProductionCutsTable* theCoupleTable= 337 G4ProductionCutsTable::GetProductionCutsTa 337 G4ProductionCutsTable::GetProductionCutsTable(); 338 std::size_t numOfCouples = theCoupleTable->G 338 std::size_t numOfCouples = theCoupleTable->GetTableSize(); 339 339 340 // prepare vector 340 // prepare vector 341 auto elmSelectors = mod->GetElementSelectors 341 auto elmSelectors = mod->GetElementSelectors(); 342 if(nullptr == elmSelectors) { 342 if(nullptr == elmSelectors) { 343 elmSelectors = new std::vector<G4EmElement 343 elmSelectors = new std::vector<G4EmElementSelector*>; 344 } 344 } 345 std::size_t nSelectors = elmSelectors->size( 345 std::size_t nSelectors = elmSelectors->size(); 346 if(numOfCouples > nSelectors) { 346 if(numOfCouples > nSelectors) { 347 for(std::size_t i=nSelectors; i<numOfCoupl 347 for(std::size_t i=nSelectors; i<numOfCouples; ++i) { 348 elmSelectors->push_back(nullptr); 348 elmSelectors->push_back(nullptr); 349 } 349 } 350 nSelectors = numOfCouples; 350 nSelectors = numOfCouples; 351 } 351 } 352 352 353 // initialise vector 353 // initialise vector 354 for(std::size_t i=0; i<numOfCouples; ++i) { 354 for(std::size_t i=0; i<numOfCouples; ++i) { 355 355 356 // no need in element selectors for infini 356 // no need in element selectors for infinite cuts 357 if(cuts[i] == DBL_MAX) { continue; } 357 if(cuts[i] == DBL_MAX) { continue; } 358 358 359 auto couple = theCoupleTable->GetMaterialC 359 auto couple = theCoupleTable->GetMaterialCutsCouple((G4int)i); 360 auto mat = couple->GetMaterial(); 360 auto mat = couple->GetMaterial(); 361 mod->SetCurrentCouple(couple); 361 mod->SetCurrentCouple(couple); 362 362 363 // selector already exist then delete 363 // selector already exist then delete 364 delete (*elmSelectors)[i]; 364 delete (*elmSelectors)[i]; 365 365 366 G4double emin = std::max(elow, mod->MinPri 366 G4double emin = std::max(elow, mod->MinPrimaryEnergy(mat, part, cuts[i])); 367 G4double emax = std::max(ehigh, 10*emin); 367 G4double emax = std::max(ehigh, 10*emin); 368 static const G4double invlog106 = 1.0/(6*G 368 static const G4double invlog106 = 1.0/(6*G4Log(10.)); 369 G4int nbins = G4lrint(nbinsPerDec*G4Log(em 369 G4int nbins = G4lrint(nbinsPerDec*G4Log(emax/emin)*invlog106); 370 nbins = std::max(nbins, 3); 370 nbins = std::max(nbins, 3); 371 371 372 (*elmSelectors)[i] = new G4EmElementSelect 372 (*elmSelectors)[i] = new G4EmElementSelector(mod,mat,nbins, 373 emin,emax,spline); 373 emin,emax,spline); 374 ((*elmSelectors)[i])->Initialise(part, cut 374 ((*elmSelectors)[i])->Initialise(part, cuts[i]); 375 /* 375 /* 376 G4cout << "G4VEmModel::InitialiseElmSele 376 G4cout << "G4VEmModel::InitialiseElmSelectors i= " << i 377 << " " << part->GetParticleName 377 << " " << part->GetParticleName() 378 << " for " << mod->GetName() << " 378 << " for " << mod->GetName() << " cut= " << cuts[i] 379 << " " << (*elmSelectors)[i] << 379 << " " << (*elmSelectors)[i] << G4endl; 380 ((*elmSelectors)[i])->Dump(part); 380 ((*elmSelectors)[i])->Dump(part); 381 */ 381 */ 382 } 382 } 383 mod->SetElementSelectors(elmSelectors); 383 mod->SetElementSelectors(elmSelectors); 384 } 384 } 385 385 386 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 387 387