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 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4EmModelManager 32 // File name: G4EmModelManager 33 // 33 // 34 // Author: Vladimir Ivanchenko 34 // Author: Vladimir Ivanchenko 35 // 35 // 36 // Creation date: 07.05.2002 36 // Creation date: 07.05.2002 37 // 37 // 38 // Modifications: V.Ivanchenko 38 // Modifications: V.Ivanchenko 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // 41 // 42 // It is the unified energy loss process it ca 42 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 43 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 44 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 45 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 46 // that information on fly. 47 // ------------------------------------------- 47 // ------------------------------------------------------------------- 48 // 48 // 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 51 52 #include "G4EmModelManager.hh" 52 #include "G4EmModelManager.hh" 53 #include "G4SystemOfUnits.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4PhysicsTable.hh" 54 #include "G4PhysicsTable.hh" 55 #include "G4PhysicsVector.hh" 55 #include "G4PhysicsVector.hh" 56 #include "G4VMscModel.hh" 56 #include "G4VMscModel.hh" 57 57 58 #include "G4Step.hh" << 59 #include "G4ParticleDefinition.hh" << 60 #include "G4PhysicsVector.hh" << 61 #include "G4MaterialCutsCouple.hh" << 62 #include "G4ProductionCutsTable.hh" << 63 #include "G4RegionStore.hh" << 64 #include "G4Gamma.hh" << 65 #include "G4Electron.hh" << 66 #include "G4Positron.hh" << 67 #include "G4UnitsTable.hh" << 68 #include "G4DataVector.hh" << 69 << 70 //....oooOO0OOooo........oooOO0OOooo........oo 58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 59 72 G4RegionModels::G4RegionModels(G4int nMod, std 60 G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx, 73 G4DataVector& l 61 G4DataVector& lowE, const G4Region* reg) 74 { 62 { 75 nModelsForRegion = nMod; 63 nModelsForRegion = nMod; 76 theListOfModelIndexes = new G4int [nModelsFo 64 theListOfModelIndexes = new G4int [nModelsForRegion]; 77 lowKineticEnergy = new G4double [nModel 65 lowKineticEnergy = new G4double [nModelsForRegion+1]; 78 for (G4int i=0; i<nModelsForRegion; ++i) { 66 for (G4int i=0; i<nModelsForRegion; ++i) { 79 theListOfModelIndexes[i] = indx[i]; 67 theListOfModelIndexes[i] = indx[i]; 80 lowKineticEnergy[i] = lowE[i]; 68 lowKineticEnergy[i] = lowE[i]; 81 } 69 } 82 lowKineticEnergy[nModelsForRegion] = lowE[nM 70 lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion]; 83 theRegion = reg; 71 theRegion = reg; 84 } 72 } 85 73 86 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 87 75 88 G4RegionModels::~G4RegionModels() 76 G4RegionModels::~G4RegionModels() 89 { 77 { 90 delete [] theListOfModelIndexes; 78 delete [] theListOfModelIndexes; 91 delete [] lowKineticEnergy; 79 delete [] lowKineticEnergy; 92 } 80 } 93 81 94 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 83 >> 84 #include "G4Step.hh" >> 85 #include "G4ParticleDefinition.hh" >> 86 #include "G4PhysicsVector.hh" >> 87 #include "G4MaterialCutsCouple.hh" >> 88 #include "G4ProductionCutsTable.hh" >> 89 #include "G4RegionStore.hh" >> 90 #include "G4Gamma.hh" >> 91 #include "G4Electron.hh" >> 92 #include "G4Positron.hh" >> 93 #include "G4UnitsTable.hh" >> 94 #include "G4DataVector.hh" >> 95 95 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 96 97 97 G4EmModelManager::G4EmModelManager() << 98 G4EmModelManager::G4EmModelManager(): >> 99 nEmModels(0), >> 100 nRegions(0), >> 101 particle(0), >> 102 verboseLevel(0) 98 { 103 { >> 104 maxSubCutInRange = 0.7*mm; 99 models.reserve(4); 105 models.reserve(4); 100 flucModels.reserve(4); 106 flucModels.reserve(4); 101 regions.reserve(4); 107 regions.reserve(4); 102 orderOfModels.reserve(4); 108 orderOfModels.reserve(4); 103 isUsed.reserve(4); 109 isUsed.reserve(4); >> 110 severalModels = true; >> 111 fluoFlag = false; >> 112 currRegionModel = nullptr; >> 113 currModel = nullptr; >> 114 theCuts = nullptr; >> 115 theCutsNew = nullptr; >> 116 theSubCuts = nullptr; 104 } 117 } 105 118 106 //....oooOO0OOooo........oooOO0OOooo........oo 119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 120 108 G4EmModelManager::~G4EmModelManager() 121 G4EmModelManager::~G4EmModelManager() 109 { 122 { 110 verboseLevel = 0; // no verbosity at destruc 123 verboseLevel = 0; // no verbosity at destruction 111 Clear(); 124 Clear(); 112 delete theCutsNew; 125 delete theCutsNew; >> 126 delete theSubCuts; 113 } 127 } 114 128 115 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 130 117 void G4EmModelManager::Clear() 131 void G4EmModelManager::Clear() 118 { 132 { 119 if(1 < verboseLevel) { 133 if(1 < verboseLevel) { 120 G4cout << "G4EmModelManager::Clear()" << G 134 G4cout << "G4EmModelManager::Clear()" << G4endl; 121 } 135 } 122 std::size_t n = setOfRegionModels.size(); << 136 size_t n = setOfRegionModels.size(); 123 for(std::size_t i=0; i<n; ++i) { << 137 if(n > 0) { 124 delete setOfRegionModels[i]; << 138 for(size_t i=0; i<n; ++i) { 125 setOfRegionModels[i] = nullptr; << 139 delete setOfRegionModels[i]; >> 140 setOfRegionModels[i] = nullptr; >> 141 } 126 } 142 } 127 } 143 } 128 144 129 //....oooOO0OOooo........oooOO0OOooo........oo 145 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 130 146 131 void G4EmModelManager::AddEmModel(G4int num, G 147 void G4EmModelManager::AddEmModel(G4int num, G4VEmModel* p, 132 G4VEmFluctua 148 G4VEmFluctuationModel* fm, const G4Region* r) 133 { 149 { 134 if(nullptr == p) { << 150 if(!p) { 135 G4cout << "G4EmModelManager::AddEmModel WA 151 G4cout << "G4EmModelManager::AddEmModel WARNING: no model defined." 136 << G4endl; 152 << G4endl; 137 return; 153 return; 138 } 154 } 139 models.push_back(p); 155 models.push_back(p); 140 flucModels.push_back(fm); 156 flucModels.push_back(fm); 141 regions.push_back(r); 157 regions.push_back(r); 142 orderOfModels.push_back(num); 158 orderOfModels.push_back(num); 143 isUsed.push_back(0); 159 isUsed.push_back(0); 144 p->DefineForRegion(r); 160 p->DefineForRegion(r); 145 ++nEmModels; 161 ++nEmModels; 146 } 162 } 147 163 148 //....oooOO0OOooo........oooOO0OOooo........oo 164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 149 165 150 G4VEmModel* G4EmModelManager::GetModel(G4int i << 166 void G4EmModelManager::UpdateEmModel(const G4String& nam, >> 167 G4double emin, G4double emax) >> 168 { >> 169 if (nEmModels > 0) { >> 170 for(G4int i=0; i<nEmModels; ++i) { >> 171 if(nam == models[i]->GetName()) { >> 172 models[i]->SetLowEnergyLimit(emin); >> 173 models[i]->SetHighEnergyLimit(emax); >> 174 break; >> 175 } >> 176 } >> 177 } >> 178 G4cout << "G4EmModelManager::UpdateEmModel WARNING: no model <" >> 179 << nam << "> is found out" >> 180 << G4endl; >> 181 } >> 182 >> 183 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 184 >> 185 G4VEmModel* G4EmModelManager::GetModel(G4int i, G4bool ver) 151 { 186 { 152 G4VEmModel* model = nullptr; 187 G4VEmModel* model = nullptr; 153 if(idx >= 0 && idx < nEmModels) { model = mo << 188 if(i < nEmModels) { model = models[i]; } 154 else if(verboseLevel > 0 && ver) { 189 else if(verboseLevel > 0 && ver) { 155 G4cout << "G4EmModelManager::GetModel WARN 190 G4cout << "G4EmModelManager::GetModel WARNING: " 156 << "index " << idx << " is wrong Nm << 191 << "index " << i << " is wrong Nmodels= " 157 << nEmModels; 192 << nEmModels; 158 if(nullptr != particle) { << 193 if(particle) { G4cout << " for " << particle->GetParticleName(); } 159 G4cout << " for " << particle->GetPartic << 160 } << 161 G4cout<< G4endl; 194 G4cout<< G4endl; 162 } 195 } 163 return model; 196 return model; 164 } 197 } 165 198 166 //....oooOO0OOooo........oooOO0OOooo........oo 199 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 200 168 G4VEmModel* G4EmModelManager::GetRegionModel(G << 201 G4VEmModel* G4EmModelManager::GetRegionModel(G4int k, size_t idx) 169 { 202 { 170 G4RegionModels* rm = setOfRegionModels[idxOf 203 G4RegionModels* rm = setOfRegionModels[idxOfRegionModels[idx]]; 171 return (k < rm->NumberOfModels()) ? models[r << 204 G4VEmModel* mod = models[rm->ModelIndex(k)]; >> 205 return mod; 172 } 206 } 173 207 174 //....oooOO0OOooo........oooOO0OOooo........oo 208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 175 209 176 G4int G4EmModelManager::NumberOfRegionModels(s << 210 G4int G4EmModelManager::NumberOfRegionModels(size_t idx) const 177 { 211 { 178 G4RegionModels* rm = setOfRegionModels[idxOf 212 G4RegionModels* rm = setOfRegionModels[idxOfRegionModels[idx]]; 179 return rm->NumberOfModels(); 213 return rm->NumberOfModels(); 180 } 214 } 181 215 182 //....oooOO0OOooo........oooOO0OOooo........oo 216 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 183 217 184 const G4DataVector* 218 const G4DataVector* 185 G4EmModelManager::Initialise(const G4ParticleD 219 G4EmModelManager::Initialise(const G4ParticleDefinition* p, 186 const G4ParticleD 220 const G4ParticleDefinition* secondaryParticle, 187 G4int verb) << 221 G4double minSubRange, >> 222 G4int val) 188 { 223 { 189 verboseLevel = verb; << 224 verboseLevel = val; >> 225 G4String partname = p->GetParticleName(); 190 if(1 < verboseLevel) { 226 if(1 < verboseLevel) { 191 G4cout << "G4EmModelManager::Initialise() 227 G4cout << "G4EmModelManager::Initialise() for " 192 << p->GetParticleName() << " Nmode << 228 << partname << " Nmodels= " << nEmModels << G4endl; 193 } 229 } 194 // Are models defined? 230 // Are models defined? 195 if(nEmModels < 1) { 231 if(nEmModels < 1) { 196 G4ExceptionDescription ed; 232 G4ExceptionDescription ed; 197 ed << "No models found out for " << p->Get 233 ed << "No models found out for " << p->GetParticleName() 198 << " !"; 234 << " !"; 199 G4Exception("G4EmModelManager::Initialise" 235 G4Exception("G4EmModelManager::Initialise","em0002", 200 FatalException, ed); 236 FatalException, ed); 201 } 237 } 202 238 203 particle = p; 239 particle = p; 204 Clear(); // needed if run is not first 240 Clear(); // needed if run is not first 205 241 206 G4RegionStore* regionStore = G4RegionStore:: 242 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 207 const G4Region* world = 243 const G4Region* world = 208 regionStore->GetRegion("DefaultRegionForTh 244 regionStore->GetRegion("DefaultRegionForTheWorld", false); 209 245 210 // Identify the list of regions with differe 246 // Identify the list of regions with different set of models 211 nRegions = 1; 247 nRegions = 1; 212 std::vector<const G4Region*> setr; 248 std::vector<const G4Region*> setr; 213 setr.push_back(world); 249 setr.push_back(world); 214 G4bool isWorld = false; 250 G4bool isWorld = false; 215 251 216 for (G4int ii=0; ii<nEmModels; ++ii) { 252 for (G4int ii=0; ii<nEmModels; ++ii) { 217 const G4Region* r = regions[ii]; 253 const G4Region* r = regions[ii]; 218 if ( r == nullptr || r == world) { << 254 if ( r == 0 || r == world) { 219 isWorld = true; 255 isWorld = true; 220 regions[ii] = world; 256 regions[ii] = world; 221 } else { 257 } else { 222 G4bool newRegion = true; 258 G4bool newRegion = true; 223 if (nRegions>1) { 259 if (nRegions>1) { 224 for (G4int j=1; j<nRegions; ++j) { 260 for (G4int j=1; j<nRegions; ++j) { 225 if ( r == setr[j] ) { newRegion = fa 261 if ( r == setr[j] ) { newRegion = false; } 226 } 262 } 227 } 263 } 228 if (newRegion) { 264 if (newRegion) { 229 setr.push_back(r); 265 setr.push_back(r); 230 ++nRegions; << 266 nRegions++; 231 } 267 } 232 } 268 } 233 } 269 } 234 // Are models defined? 270 // Are models defined? 235 if(!isWorld) { 271 if(!isWorld) { 236 G4ExceptionDescription ed; 272 G4ExceptionDescription ed; 237 ed << "No models defined for the World vol 273 ed << "No models defined for the World volume for " 238 << p->GetParticleName() << " !"; 274 << p->GetParticleName() << " !"; 239 G4Exception("G4EmModelManager::Initialise" 275 G4Exception("G4EmModelManager::Initialise","em0002", 240 FatalException, ed); 276 FatalException, ed); 241 } 277 } 242 278 243 G4ProductionCutsTable* theCoupleTable= 279 G4ProductionCutsTable* theCoupleTable= 244 G4ProductionCutsTable::GetProductionCutsTa 280 G4ProductionCutsTable::GetProductionCutsTable(); 245 std::size_t numOfCouples = theCoupleTable->G << 281 size_t numOfCouples = theCoupleTable->GetTableSize(); 246 282 247 // prepare vectors, shortcut for the case of 283 // prepare vectors, shortcut for the case of only 1 model 248 // or only one region 284 // or only one region 249 if(nRegions > 1 && nEmModels > 1) { 285 if(nRegions > 1 && nEmModels > 1) { 250 idxOfRegionModels.resize(numOfCouples,0); 286 idxOfRegionModels.resize(numOfCouples,0); 251 setOfRegionModels.resize((std::size_t)nReg << 287 setOfRegionModels.resize((size_t)nRegions,0); 252 } else { 288 } else { 253 idxOfRegionModels.resize(1,0); 289 idxOfRegionModels.resize(1,0); 254 setOfRegionModels.resize(1,nullptr); << 290 setOfRegionModels.resize(1,0); 255 } 291 } 256 292 257 std::vector<G4int> modelAtRegion(nEmModel 293 std::vector<G4int> modelAtRegion(nEmModels); 258 std::vector<G4int> modelOrd(nEmModels); 294 std::vector<G4int> modelOrd(nEmModels); 259 G4DataVector eLow(nEmModels+1); 295 G4DataVector eLow(nEmModels+1); 260 G4DataVector eHigh(nEmModels); 296 G4DataVector eHigh(nEmModels); 261 297 262 if(1 < verboseLevel) { 298 if(1 < verboseLevel) { 263 G4cout << " Nregions= " << nRegions 299 G4cout << " Nregions= " << nRegions 264 << " Nmodels= " << nEmModels << G4 300 << " Nmodels= " << nEmModels << G4endl; 265 } 301 } 266 302 267 // Order models for regions 303 // Order models for regions 268 for (G4int reg=0; reg<nRegions; ++reg) { 304 for (G4int reg=0; reg<nRegions; ++reg) { 269 const G4Region* region = setr[reg]; 305 const G4Region* region = setr[reg]; 270 G4int n = 0; 306 G4int n = 0; 271 307 272 for (G4int ii=0; ii<nEmModels; ++ii) { 308 for (G4int ii=0; ii<nEmModels; ++ii) { 273 309 274 G4VEmModel* model = models[ii]; 310 G4VEmModel* model = models[ii]; 275 if ( region == regions[ii] ) { 311 if ( region == regions[ii] ) { 276 312 277 G4double tmin = model->LowEnergyLimit( 313 G4double tmin = model->LowEnergyLimit(); 278 G4double tmax = model->HighEnergyLimit 314 G4double tmax = model->HighEnergyLimit(); 279 G4int ord = orderOfModels[ii]; 315 G4int ord = orderOfModels[ii]; 280 G4bool push = true; 316 G4bool push = true; 281 G4bool insert = false; 317 G4bool insert = false; 282 G4int idx = n; 318 G4int idx = n; 283 319 284 if(1 < verboseLevel) { 320 if(1 < verboseLevel) { 285 G4cout << "Model #" << ii 321 G4cout << "Model #" << ii 286 << " <" << model->GetName() < 322 << " <" << model->GetName() << "> for region <"; 287 if (region) G4cout << region->GetNam 323 if (region) G4cout << region->GetName(); 288 G4cout << "> " 324 G4cout << "> " 289 << " tmin(MeV)= " << tmin/MeV 325 << " tmin(MeV)= " << tmin/MeV 290 << "; tmax(MeV)= " << tmax/Me 326 << "; tmax(MeV)= " << tmax/MeV 291 << "; order= " << ord 327 << "; order= " << ord 292 << "; tminAct= " << model->LowEnergyActiv 328 << "; tminAct= " << model->LowEnergyActivationLimit()/MeV 293 << "; tmaxAct= " << model->HighEnergyActi 329 << "; tmaxAct= " << model->HighEnergyActivationLimit()/MeV 294 << G4endl; 330 << G4endl; 295 } 331 } 296 332 297 static const G4double limitdelta = 0.01*eV; 333 static const G4double limitdelta = 0.01*eV; 298 if(n > 0) { 334 if(n > 0) { 299 335 300 // extend energy range to previous m 336 // extend energy range to previous models 301 tmin = std::min(tmin, eHigh[n-1]); 337 tmin = std::min(tmin, eHigh[n-1]); 302 tmax = std::max(tmax, eLow[0]); 338 tmax = std::max(tmax, eLow[0]); 303 //G4cout << "tmin= " << tmin << " t 339 //G4cout << "tmin= " << tmin << " tmax= " 304 // << tmax << " ord= " << 340 // << tmax << " ord= " << ord <<G4endl; 305 // empty energy range 341 // empty energy range 306 if( tmax - tmin <= limitdelta) { pus 342 if( tmax - tmin <= limitdelta) { push = false; } 307 // low-energy model 343 // low-energy model 308 else if (tmax == eLow[0]) { 344 else if (tmax == eLow[0]) { 309 push = false; 345 push = false; 310 insert = true; 346 insert = true; 311 idx = 0; 347 idx = 0; 312 // resolve intersections 348 // resolve intersections 313 } else if(tmin < eHigh[n-1]) { 349 } else if(tmin < eHigh[n-1]) { 314 // compare order 350 // compare order 315 for(G4int k=0; k<n; ++k) { 351 for(G4int k=0; k<n; ++k) { 316 // new model has higher order pa 352 // new model has higher order parameter, 317 // so, its application area may be red 353 // so, its application area may be reduced 318 // to avoid intersections 354 // to avoid intersections 319 if(ord >= modelOrd[k]) { 355 if(ord >= modelOrd[k]) { 320 if(tmin < eHigh[k] && tmin >= 356 if(tmin < eHigh[k] && tmin >= eLow[k]) { tmin = eHigh[k]; } 321 if(tmax <= eHigh[k] && tmax > 357 if(tmax <= eHigh[k] && tmax > eLow[k]) { tmax = eLow[k]; } 322 if(tmax > eHigh[k] && tmin < e 358 if(tmax > eHigh[k] && tmin < eLow[k]) { 323 if(tmax - eHigh[k] > eLow[k] 359 if(tmax - eHigh[k] > eLow[k] - tmin) { tmin = eHigh[k]; } 324 else { tmax = eLow[k]; } 360 else { tmax = eLow[k]; } 325 } 361 } 326 if( tmax - tmin <= limitdelta) 362 if( tmax - tmin <= limitdelta) { 327 push = false; 363 push = false; 328 break; 364 break; 329 } 365 } 330 } 366 } 331 } 367 } 332 // this model has lower order parameter 368 // this model has lower order parameter than possible 333 // other models, with which there may be 369 // other models, with which there may be intersections 334 // so, appliction area of such models ma 370 // so, appliction area of such models may be reduced 335 371 336 // insert below the first model 372 // insert below the first model 337 if (tmax <= eLow[0]) { 373 if (tmax <= eLow[0]) { 338 push = false; 374 push = false; 339 insert = true; 375 insert = true; 340 idx = 0; 376 idx = 0; 341 // resolve intersections 377 // resolve intersections 342 } else if(tmin < eHigh[n-1]) { 378 } else if(tmin < eHigh[n-1]) { 343 // last energy interval 379 // last energy interval 344 if(tmin > eLow[n-1] && tmax >= eHigh[n 380 if(tmin > eLow[n-1] && tmax >= eHigh[n-1]) { 345 eHigh[n-1] = tmin; 381 eHigh[n-1] = tmin; 346 // first energy interval 382 // first energy interval 347 } else if(tmin <= eLow[0] && tmax < eH 383 } else if(tmin <= eLow[0] && tmax < eHigh[0]) { 348 eLow[0] = tmax; 384 eLow[0] = tmax; 349 push = false; 385 push = false; 350 insert = true; 386 insert = true; 351 idx = 0; 387 idx = 0; 352 // loop over all models 388 // loop over all models 353 } else { 389 } else { 354 for(G4int k=n-1; k>=0; --k) { 390 for(G4int k=n-1; k>=0; --k) { 355 if(tmin <= eLow[k] && tmax >= eHigh[k]) 391 if(tmin <= eLow[k] && tmax >= eHigh[k]) { 356 // full overlap exclude previous model 392 // full overlap exclude previous model 357 isUsed[modelAtRegion[k]] = 0; 393 isUsed[modelAtRegion[k]] = 0; 358 idx = k; 394 idx = k; 359 if(k < n-1) { 395 if(k < n-1) { 360 // shift upper models and change ind 396 // shift upper models and change index 361 for(G4int kk=k; kk<n-1; ++kk) { 397 for(G4int kk=k; kk<n-1; ++kk) { 362 modelAtRegion[kk] = modelAtRegion[kk+1]; 398 modelAtRegion[kk] = modelAtRegion[kk+1]; 363 modelOrd[kk] = modelOrd[kk+1]; 399 modelOrd[kk] = modelOrd[kk+1]; 364 eLow[kk] = eLow[kk+1]; 400 eLow[kk] = eLow[kk+1]; 365 eHigh[kk] = eHigh[kk+1]; 401 eHigh[kk] = eHigh[kk+1]; 366 } 402 } 367 ++k; 403 ++k; 368 } 404 } 369 --n; 405 --n; 370 } else { 406 } else { 371 // partially reduce previous model are 407 // partially reduce previous model area 372 if(tmin <= eLow[k] && tmax > eLow[k]) 408 if(tmin <= eLow[k] && tmax > eLow[k]) { 373 eLow[k] = tmax; 409 eLow[k] = tmax; 374 idx = k; 410 idx = k; 375 insert = true; 411 insert = true; 376 push = false; 412 push = false; 377 } else if(tmin < eHigh[k] && tmax >= e 413 } else if(tmin < eHigh[k] && tmax >= eHigh[k]) { 378 eHigh[k] = tmin; 414 eHigh[k] = tmin; 379 idx = k + 1; 415 idx = k + 1; 380 if(idx < n) { 416 if(idx < n) { 381 insert = true; 417 insert = true; 382 push = false; 418 push = false; 383 } 419 } 384 } else if(tmin > eLow[k] && tmax < eHi 420 } else if(tmin > eLow[k] && tmax < eHigh[k]) { 385 if(eHigh[k] - tmax > tmin - eLow[k]) 421 if(eHigh[k] - tmax > tmin - eLow[k]) { 386 eLow[k] = tmax; 422 eLow[k] = tmax; 387 idx = k; 423 idx = k; 388 insert = true; 424 insert = true; 389 push = false; 425 push = false; 390 } else { 426 } else { 391 eHigh[k] = tmin; 427 eHigh[k] = tmin; 392 idx = k + 1; 428 idx = k + 1; 393 if(idx < n) { 429 if(idx < n) { 394 insert = true; 430 insert = true; 395 push = false; 431 push = false; 396 } 432 } 397 } 433 } 398 } 434 } 399 } 435 } 400 } 436 } 401 } 437 } 402 } 438 } 403 } 439 } 404 } 440 } 405 // provide space for the new model 441 // provide space for the new model 406 if(insert) { 442 if(insert) { 407 for(G4int k=n-1; k>=idx; --k) { 443 for(G4int k=n-1; k>=idx; --k) { 408 modelAtRegion[k+1] = modelAtRegion 444 modelAtRegion[k+1] = modelAtRegion[k]; 409 modelOrd[k+1] = modelOrd[k]; 445 modelOrd[k+1] = modelOrd[k]; 410 eLow[k+1] = eLow[k]; 446 eLow[k+1] = eLow[k]; 411 eHigh[k+1] = eHigh[k]; 447 eHigh[k+1] = eHigh[k]; 412 } 448 } 413 } 449 } 414 //G4cout << "push= " << push << " inse 450 //G4cout << "push= " << push << " insert= " << insert 415 // << " idx= " << idx <<G4endl; 451 // << " idx= " << idx <<G4endl; 416 // the model is added 452 // the model is added 417 if (push || insert) { 453 if (push || insert) { 418 ++n; 454 ++n; 419 modelAtRegion[idx] = ii; 455 modelAtRegion[idx] = ii; 420 modelOrd[idx] = ord; 456 modelOrd[idx] = ord; 421 eLow[idx] = tmin; 457 eLow[idx] = tmin; 422 eHigh[idx] = tmax; 458 eHigh[idx] = tmax; 423 isUsed[ii] = 1; 459 isUsed[ii] = 1; 424 } 460 } 425 // exclude models with zero energy range 461 // exclude models with zero energy range 426 for(G4int k=n-1; k>=0; --k) { 462 for(G4int k=n-1; k>=0; --k) { 427 if(eHigh[k] - eLow[k] <= limitdelta) { 463 if(eHigh[k] - eLow[k] <= limitdelta) { 428 isUsed[modelAtRegion[k]] = 0; 464 isUsed[modelAtRegion[k]] = 0; 429 if(k < n-1) { 465 if(k < n-1) { 430 for(G4int kk=k; kk<n-1; ++kk) { 466 for(G4int kk=k; kk<n-1; ++kk) { 431 modelAtRegion[kk] = modelAtRegion[kk+1]; 467 modelAtRegion[kk] = modelAtRegion[kk+1]; 432 modelOrd[kk] = modelOrd[kk+1]; 468 modelOrd[kk] = modelOrd[kk+1]; 433 eLow[kk] = eLow[kk+1]; 469 eLow[kk] = eLow[kk+1]; 434 eHigh[kk] = eHigh[kk+1]; 470 eHigh[kk] = eHigh[kk+1]; 435 } 471 } 436 } 472 } 437 --n; 473 --n; 438 } 474 } 439 } 475 } 440 } 476 } 441 } 477 } 442 eLow[0] = 0.0; 478 eLow[0] = 0.0; 443 eLow[n] = eHigh[n-1]; 479 eLow[n] = eHigh[n-1]; 444 480 445 if(1 < verboseLevel) { 481 if(1 < verboseLevel) { 446 G4cout << "### New G4RegionModels set wi 482 G4cout << "### New G4RegionModels set with " << n 447 << " models for region <"; 483 << " models for region <"; 448 if (region) { G4cout << region->GetName( 484 if (region) { G4cout << region->GetName(); } 449 G4cout << "> Elow(MeV)= "; 485 G4cout << "> Elow(MeV)= "; 450 for(G4int iii=0; iii<=n; ++iii) {G4cout 486 for(G4int iii=0; iii<=n; ++iii) {G4cout << eLow[iii]/MeV << " ";} 451 G4cout << G4endl; 487 G4cout << G4endl; 452 } 488 } 453 auto rm = new G4RegionModels(n, modelAtReg << 489 G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region); 454 setOfRegionModels[reg] = rm; 490 setOfRegionModels[reg] = rm; 455 // shortcut 491 // shortcut 456 if(1 == nEmModels) { break; } 492 if(1 == nEmModels) { break; } 457 } 493 } 458 494 459 currRegionModel = setOfRegionModels[0]; 495 currRegionModel = setOfRegionModels[0]; 460 currModel = models[0]; 496 currModel = models[0]; 461 497 462 // Access to materials and build cuts 498 // Access to materials and build cuts 463 std::size_t idx = 1; << 499 size_t idx = 1; 464 if(nullptr != secondaryParticle) { << 500 if(secondaryParticle) { 465 if( secondaryParticle == G4Gamma::Gamma() 501 if( secondaryParticle == G4Gamma::Gamma() ) { idx = 0; } 466 else if( secondaryParticle == G4Electron:: 502 else if( secondaryParticle == G4Electron::Electron()) { idx = 1; } 467 else if( secondaryParticle == G4Positron:: 503 else if( secondaryParticle == G4Positron::Positron()) { idx = 2; } 468 else { idx = 3; } 504 else { idx = 3; } 469 } 505 } 470 506 471 theCuts = 507 theCuts = 472 static_cast<const G4DataVector*>(theCouple 508 static_cast<const G4DataVector*>(theCoupleTable->GetEnergyCutsVector(idx)); 473 509 474 // for the second run the check on cuts shou 510 // for the second run the check on cuts should be repeated 475 if(nullptr != theCutsNew) { *theCutsNew = *t << 511 if(theCutsNew) { *theCutsNew = *theCuts; } >> 512 >> 513 if(minSubRange < 1.0) { >> 514 if( !theSubCuts ) { theSubCuts = new G4DataVector(); } >> 515 theSubCuts->resize(numOfCouples,DBL_MAX); >> 516 } 476 517 477 // G4cout << "========Start define cuts" << 518 // G4cout << "========Start define cuts" << G4endl; 478 // define cut values 519 // define cut values 479 for(std::size_t i=0; i<numOfCouples; ++i) { << 520 for(size_t i=0; i<numOfCouples; ++i) { 480 521 481 const G4MaterialCutsCouple* couple = 522 const G4MaterialCutsCouple* couple = 482 theCoupleTable->GetMaterialCutsCouple((G << 523 theCoupleTable->GetMaterialCutsCouple(i); 483 const G4Material* material = couple->GetMa 524 const G4Material* material = couple->GetMaterial(); 484 const G4ProductionCuts* pcuts = couple->Ge 525 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 485 526 486 G4int reg = 0; 527 G4int reg = 0; 487 if(nRegions > 1 && nEmModels > 1) { 528 if(nRegions > 1 && nEmModels > 1) { 488 reg = nRegions; 529 reg = nRegions; 489 // Loop checking, 03-Aug-2015, Vladimir 530 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 490 do {--reg;} while (reg>0 && pcuts != (se 531 do {--reg;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts())); 491 idxOfRegionModels[i] = reg; 532 idxOfRegionModels[i] = reg; 492 } 533 } 493 if(1 < verboseLevel) { 534 if(1 < verboseLevel) { 494 G4cout << "G4EmModelManager::Initialise( 535 G4cout << "G4EmModelManager::Initialise() for " 495 << material->GetName() 536 << material->GetName() 496 << " indexOfCouple= " << i 537 << " indexOfCouple= " << i 497 << " indexOfRegion= " << reg 538 << " indexOfRegion= " << reg 498 << G4endl; 539 << G4endl; 499 } 540 } 500 541 501 G4double cut = (*theCuts)[i]; 542 G4double cut = (*theCuts)[i]; 502 if(nullptr != secondaryParticle) { << 543 if(secondaryParticle) { >> 544 >> 545 // compute subcut >> 546 if( cut < DBL_MAX && minSubRange < 1.0) { >> 547 G4double subcut = minSubRange*cut; >> 548 G4double rcut = std::min(minSubRange*pcuts->GetProductionCut(idx), >> 549 maxSubCutInRange); >> 550 G4double tcutmax = >> 551 theCoupleTable->ConvertRangeToEnergy(secondaryParticle, >> 552 material,rcut); >> 553 if(tcutmax < subcut) { subcut = tcutmax; } >> 554 (*theSubCuts)[i] = subcut; >> 555 } 503 556 504 // note that idxOfRegionModels[] not alw 557 // note that idxOfRegionModels[] not always filled 505 G4int inn = 0; 558 G4int inn = 0; 506 G4int nnm = 1; 559 G4int nnm = 1; 507 if(nRegions > 1 && nEmModels > 1) { 560 if(nRegions > 1 && nEmModels > 1) { 508 inn = idxOfRegionModels[i]; 561 inn = idxOfRegionModels[i]; 509 } 562 } 510 // check cuts and introduce upper limits 563 // check cuts and introduce upper limits 511 //G4cout << "idx= " << i << " cut(keV)= 564 //G4cout << "idx= " << i << " cut(keV)= " << cut/keV << G4endl; 512 currRegionModel = setOfRegionModels[inn] 565 currRegionModel = setOfRegionModels[inn]; 513 nnm = currRegionModel->NumberOfModels(); 566 nnm = currRegionModel->NumberOfModels(); 514 567 515 //G4cout << "idx= " << i << " Nmod= " << 568 //G4cout << "idx= " << i << " Nmod= " << nnm << G4endl; 516 569 517 for(G4int jj=0; jj<nnm; ++jj) { 570 for(G4int jj=0; jj<nnm; ++jj) { 518 //G4cout << "jj= " << jj << " modidx= 571 //G4cout << "jj= " << jj << " modidx= " 519 // << currRegionModel->ModelInde 572 // << currRegionModel->ModelIndex(jj) << G4endl; 520 currModel = models[currRegionModel->Mo 573 currModel = models[currRegionModel->ModelIndex(jj)]; 521 G4double cutlim = currModel->MinEnergy 574 G4double cutlim = currModel->MinEnergyCut(particle,couple); 522 if(cutlim > cut) { 575 if(cutlim > cut) { 523 if(nullptr == theCutsNew) { theCutsN << 576 if(!theCutsNew) { theCutsNew = new G4DataVector(*theCuts); } 524 (*theCutsNew)[i] = cutlim; 577 (*theCutsNew)[i] = cutlim; 525 /* 578 /* 526 G4cout << "### " << partname << " en 579 G4cout << "### " << partname << " energy loss model in " 527 << material->GetName() 580 << material->GetName() 528 << " Cut was changed from " 581 << " Cut was changed from " << cut/keV << " keV to " 529 << cutlim/keV << " keV " << " 582 << cutlim/keV << " keV " << " due to " 530 << currModel->GetName() << G4 583 << currModel->GetName() << G4endl; 531 */ 584 */ 532 } 585 } 533 } 586 } 534 } 587 } 535 } 588 } 536 if(nullptr != theCutsNew) { theCuts = theCut << 589 if(theCutsNew) { theCuts = theCutsNew; } 537 590 538 // initialize models 591 // initialize models 539 G4int nn = 0; 592 G4int nn = 0; 540 severalModels = true; 593 severalModels = true; 541 for(G4int jj=0; jj<nEmModels; ++jj) { 594 for(G4int jj=0; jj<nEmModels; ++jj) { 542 if(1 == isUsed[jj]) { 595 if(1 == isUsed[jj]) { 543 ++nn; 596 ++nn; 544 currModel = models[jj]; 597 currModel = models[jj]; 545 currModel->Initialise(particle, *theCuts 598 currModel->Initialise(particle, *theCuts); 546 if(nullptr != flucModels[jj]) { flucMode << 599 if(flucModels[jj]) { flucModels[jj]->InitialiseMe(particle); } 547 } 600 } 548 } 601 } 549 if(1 == nn) { severalModels = false; } 602 if(1 == nn) { severalModels = false; } 550 603 551 if(1 < verboseLevel) { 604 if(1 < verboseLevel) { 552 G4cout << "G4EmModelManager for " << parti << 605 G4cout << "G4EmModelManager for " << partname 553 << " is initialised; nRegions= " < 606 << " is initialised; nRegions= " << nRegions 554 << " severalModels: " << severalMod 607 << " severalModels: " << severalModels 555 << G4endl; 608 << G4endl; 556 } 609 } >> 610 557 return theCuts; 611 return theCuts; 558 } 612 } 559 613 560 //....oooOO0OOooo........oooOO0OOooo........oo 614 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 561 615 562 void G4EmModelManager::FillDEDXVector(G4Physic 616 void G4EmModelManager::FillDEDXVector(G4PhysicsVector* aVector, 563 const G4 617 const G4MaterialCutsCouple* couple, 564 G4EmTabl 618 G4EmTableType tType) 565 { 619 { 566 std::size_t i = couple->GetIndex(); << 620 size_t i = couple->GetIndex(); 567 G4double cut = (fTotal == tType) ? DBL_MAX << 621 G4double cut = (*theCuts)[i]; >> 622 G4double emin = 0.0; >> 623 >> 624 if(fTotal == tType) { cut = DBL_MAX; } >> 625 else if(fSubRestricted == tType) { >> 626 emin = cut; >> 627 if(theSubCuts) { emin = (*theSubCuts)[i]; } >> 628 } 568 629 569 if(1 < verboseLevel) { 630 if(1 < verboseLevel) { 570 G4cout << "G4EmModelManager::FillDEDXVecto 631 G4cout << "G4EmModelManager::FillDEDXVector() for " 571 << couple->GetMaterial()->GetName() 632 << couple->GetMaterial()->GetName() 572 << " cut(MeV)= " << cut 633 << " cut(MeV)= " << cut >> 634 << " emin(MeV)= " << emin 573 << " Type " << tType 635 << " Type " << tType 574 << " for " << particle->GetParticl 636 << " for " << particle->GetParticleName() 575 << G4endl; 637 << G4endl; 576 } 638 } 577 639 578 G4int reg = 0; 640 G4int reg = 0; 579 if(nRegions > 1 && nEmModels > 1) { reg = id 641 if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; } 580 const G4RegionModels* regModels = setOfRegio 642 const G4RegionModels* regModels = setOfRegionModels[reg]; 581 G4int nmod = regModels->NumberOfModels(); 643 G4int nmod = regModels->NumberOfModels(); 582 644 583 // Calculate energy losses vector 645 // Calculate energy losses vector 584 std::size_t totBinsLoss = aVector->GetVector << 646 >> 647 //G4cout << "nmod= " << nmod << G4endl; >> 648 size_t totBinsLoss = aVector->GetVectorLength(); 585 G4double del = 0.0; 649 G4double del = 0.0; 586 G4int k0 = 0; 650 G4int k0 = 0; 587 651 588 for(std::size_t j=0; j<totBinsLoss; ++j) { << 652 for(size_t j=0; j<totBinsLoss; ++j) { >> 653 589 G4double e = aVector->Energy(j); 654 G4double e = aVector->Energy(j); 590 655 591 // Choose a model of energy losses 656 // Choose a model of energy losses 592 G4int k = 0; 657 G4int k = 0; 593 if (nmod > 1) { 658 if (nmod > 1) { 594 k = nmod; 659 k = nmod; 595 // Loop checking, 03-Aug-2015, Vladimir 660 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 596 do {--k;} while (k>0 && e <= regModels-> 661 do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k)); 597 //G4cout << "k= " << k << G4endl; 662 //G4cout << "k= " << k << G4endl; 598 if(k > 0 && k != k0) { 663 if(k > 0 && k != k0) { 599 k0 = k; 664 k0 = k; 600 G4double elow = regModels->LowEdgeEner 665 G4double elow = regModels->LowEdgeEnergy(k); 601 G4double dedx1 = << 666 G4double dedx1 = ComputeDEDX(models[regModels->ModelIndex(k-1)], 602 models[regModels->ModelIndex(k-1)]->Comput << 667 couple,elow,cut,emin); 603 G4double dedx2 = << 668 G4double dedx2 = ComputeDEDX(models[regModels->ModelIndex(k)], 604 models[regModels->ModelIndex(k)]->ComputeD << 669 couple,elow,cut,emin); 605 del = (dedx2 > 0.0) ? (dedx1/dedx2 - 1 << 670 del = 0.0; >> 671 if(dedx2 > 0.0) { del = (dedx1/dedx2 - 1.0)*elow; } 606 //G4cout << "elow= " << elow 672 //G4cout << "elow= " << elow 607 // << " dedx1= " << dedx1 << " d 673 // << " dedx1= " << dedx1 << " dedx2= " << dedx2 << G4endl; 608 } 674 } 609 } 675 } 610 G4double dedx = (1.0 + del/e)* << 676 G4double dedx = 611 models[regModels->ModelIndex(k)]->ComputeD << 677 ComputeDEDX(models[regModels->ModelIndex(k)],couple,e,cut,emin); >> 678 dedx *= (1.0 + del/e); 612 679 613 if(2 < verboseLevel) { 680 if(2 < verboseLevel) { 614 G4cout << "Material= " << couple->GetMat 681 G4cout << "Material= " << couple->GetMaterial()->GetName() 615 << " E(MeV)= " << e/MeV 682 << " E(MeV)= " << e/MeV 616 << " dEdx(MeV/mm)= " << dedx*mm/ 683 << " dEdx(MeV/mm)= " << dedx*mm/MeV 617 << " del= " << del*mm/MeV<< " k= 684 << " del= " << del*mm/MeV<< " k= " << k 618 << " modelIdx= " << regModels->Mo 685 << " modelIdx= " << regModels->ModelIndex(k) 619 << G4endl; 686 << G4endl; 620 } 687 } 621 dedx = std::max(dedx, 0.0); << 688 if(dedx < 0.0) { dedx = 0.0; } 622 aVector->PutValue(j, dedx); 689 aVector->PutValue(j, dedx); 623 } 690 } 624 } 691 } 625 692 626 //....oooOO0OOooo........oooOO0OOooo........oo 693 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 627 694 628 void G4EmModelManager::FillLambdaVector(G4Phys 695 void G4EmModelManager::FillLambdaVector(G4PhysicsVector* aVector, 629 const 696 const G4MaterialCutsCouple* couple, 630 G4bool 697 G4bool startFromNull, 631 G4EmTa 698 G4EmTableType tType) 632 { 699 { 633 std::size_t i = couple->GetIndex(); << 700 size_t i = couple->GetIndex(); 634 G4double cut = (*theCuts)[i]; 701 G4double cut = (*theCuts)[i]; 635 G4double tmax = DBL_MAX; 702 G4double tmax = DBL_MAX; >> 703 if (fSubRestricted == tType) { >> 704 tmax = cut; >> 705 if(theSubCuts) { cut = (*theSubCuts)[i]; } >> 706 } 636 707 637 G4int reg = 0; 708 G4int reg = 0; 638 if(nRegions > 1 && nEmModels > 1) { reg = i 709 if(nRegions > 1 && nEmModels > 1) { reg = idxOfRegionModels[i]; } 639 const G4RegionModels* regModels = setOfRegio 710 const G4RegionModels* regModels = setOfRegionModels[reg]; 640 G4int nmod = regModels->NumberOfModels(); 711 G4int nmod = regModels->NumberOfModels(); 641 if(1 < verboseLevel) { 712 if(1 < verboseLevel) { 642 G4cout << "G4EmModelManager::FillLambdaVec 713 G4cout << "G4EmModelManager::FillLambdaVector() for " 643 << particle->GetParticleName() 714 << particle->GetParticleName() 644 << " in " << couple->GetMaterial()- 715 << " in " << couple->GetMaterial()->GetName() 645 << " Emin(MeV)= " << aVector->Energ 716 << " Emin(MeV)= " << aVector->Energy(0) 646 << " Emax(MeV)= " << aVector->GetMa 717 << " Emax(MeV)= " << aVector->GetMaxEnergy() 647 << " cut= " << cut 718 << " cut= " << cut 648 << " Type " << tType 719 << " Type " << tType 649 << " nmod= " << nmod 720 << " nmod= " << nmod >> 721 << " theSubCuts " << theSubCuts 650 << G4endl; 722 << G4endl; 651 } 723 } 652 724 653 // Calculate lambda vector 725 // Calculate lambda vector 654 std::size_t totBinsLambda = aVector->GetVect << 726 size_t totBinsLambda = aVector->GetVectorLength(); 655 G4double del = 0.0; 727 G4double del = 0.0; 656 G4int k0 = 0; 728 G4int k0 = 0; 657 G4int k = 0; 729 G4int k = 0; 658 G4VEmModel* mod = models[regModels->ModelInd 730 G4VEmModel* mod = models[regModels->ModelIndex(0)]; 659 for(std::size_t j=0; j<totBinsLambda; ++j) { << 731 for(size_t j=0; j<totBinsLambda; ++j) { 660 732 661 G4double e = aVector->Energy(j); 733 G4double e = aVector->Energy(j); 662 734 663 // Choose a model 735 // Choose a model 664 if (nmod > 1) { 736 if (nmod > 1) { 665 k = nmod; 737 k = nmod; 666 // Loop checking, 03-Aug-2015, Vladimir 738 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko 667 do {--k;} while (k>0 && e <= regModels-> 739 do {--k;} while (k>0 && e <= regModels->LowEdgeEnergy(k)); 668 if(k > 0 && k != k0) { 740 if(k > 0 && k != k0) { 669 k0 = k; 741 k0 = k; 670 G4double elow = regModels->LowEdgeEner 742 G4double elow = regModels->LowEdgeEnergy(k); 671 G4VEmModel* mod1 = models[regModels->M 743 G4VEmModel* mod1 = models[regModels->ModelIndex(k-1)]; 672 G4double xs1 = mod1->CrossSection(cou 744 G4double xs1 = mod1->CrossSection(couple,particle,elow,cut,tmax); 673 mod = models[regModels->ModelIndex(k)] 745 mod = models[regModels->ModelIndex(k)]; 674 G4double xs2 = mod->CrossSection(coupl 746 G4double xs2 = mod->CrossSection(couple,particle,elow,cut,tmax); 675 del = (xs2 > 0.0) ? (xs1/xs2 - 1.0)*el << 747 del = 0.0; >> 748 if(xs2 > 0.0) { del = (xs1/xs2 - 1.0)*elow; } 676 //G4cout << "New model k=" << k << " E 749 //G4cout << "New model k=" << k << " E(MeV)= " << e/MeV 677 // << " Elow(MeV)= " << elow/MeV 750 // << " Elow(MeV)= " << elow/MeV << " del= " << del << G4endl; 678 } 751 } 679 } 752 } 680 G4double cross = (1.0 + del/e)*mod->CrossS << 753 G4double cross = mod->CrossSection(couple,particle,e,cut,tmax); >> 754 cross *= (1.0 + del/e); 681 if(fIsCrossSectionPrim == tType) { cross * 755 if(fIsCrossSectionPrim == tType) { cross *= e; } 682 756 683 if(j==0 && startFromNull) { cross = 0.0; } 757 if(j==0 && startFromNull) { cross = 0.0; } 684 758 685 if(2 < verboseLevel) { 759 if(2 < verboseLevel) { 686 G4cout << "FillLambdaVector: " << j << " 760 G4cout << "FillLambdaVector: " << j << ". e(MeV)= " << e/MeV 687 << " cross(1/mm)= " << cross*mm 761 << " cross(1/mm)= " << cross*mm 688 << " del= " << del*mm << " k= " < 762 << " del= " << del*mm << " k= " << k 689 << " modelIdx= " << regModels->Mo 763 << " modelIdx= " << regModels->ModelIndex(k) 690 << G4endl; 764 << G4endl; 691 } 765 } 692 cross = std::max(cross, 0.0); 766 cross = std::max(cross, 0.0); 693 aVector->PutValue(j, cross); 767 aVector->PutValue(j, cross); 694 } 768 } 695 } 769 } 696 770 697 //....oooOO0OOooo........oooOO0OOooo........oo 771 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 698 772 699 void G4EmModelManager::DumpModelList(std::ostr 773 void G4EmModelManager::DumpModelList(std::ostream& out, G4int verb) 700 { 774 { 701 if(verb == 0) { return; } 775 if(verb == 0) { return; } 702 for(G4int i=0; i<nRegions; ++i) { 776 for(G4int i=0; i<nRegions; ++i) { 703 G4RegionModels* r = setOfRegionModels[i]; 777 G4RegionModels* r = setOfRegionModels[i]; 704 const G4Region* reg = r->Region(); 778 const G4Region* reg = r->Region(); 705 G4int n = r->NumberOfModels(); 779 G4int n = r->NumberOfModels(); 706 if(n > 0) { 780 if(n > 0) { 707 out << " ===== EM models for the G4 781 out << " ===== EM models for the G4Region " << reg->GetName() 708 << " ======" << G4endl; 782 << " ======" << G4endl; 709 for(G4int j=0; j<n; ++j) { 783 for(G4int j=0; j<n; ++j) { 710 G4VEmModel* model = models[r->ModelInd 784 G4VEmModel* model = models[r->ModelIndex(j)]; 711 G4double emin = 785 G4double emin = 712 std::max(r->LowEdgeEnergy(j),model-> 786 std::max(r->LowEdgeEnergy(j),model->LowEnergyActivationLimit()); 713 G4double emax = 787 G4double emax = 714 std::min(r->LowEdgeEnergy(j+1),model 788 std::min(r->LowEdgeEnergy(j+1),model->HighEnergyActivationLimit()); 715 if(emax > emin) { 789 if(emax > emin) { 716 out << std::setw(20); 790 out << std::setw(20); 717 out << model->GetName() << " : Emin=" 791 out << model->GetName() << " : Emin=" 718 << std::setw(5) << G4BestUnit(emin,"En 792 << std::setw(5) << G4BestUnit(emin,"Energy") 719 << " Emax=" 793 << " Emax=" 720 << std::setw(5) << G4BestUnit(emax,"En 794 << std::setw(5) << G4BestUnit(emax,"Energy"); 721 G4PhysicsTable* table = model->GetCrossSec 795 G4PhysicsTable* table = model->GetCrossSectionTable(); 722 if(table) { 796 if(table) { 723 std::size_t kk = table->size(); << 797 size_t kk = table->size(); 724 for(std::size_t k=0; k<kk; ++k) { << 798 for(size_t k=0; k<kk; ++k) { 725 const G4PhysicsVector* v = (*table)[k] 799 const G4PhysicsVector* v = (*table)[k]; 726 if(v) { 800 if(v) { 727 G4int nn = G4int(v->GetVectorLength() - 1) << 801 G4int nn = v->GetVectorLength() - 1; 728 out << " Nbins=" << nn << " " 802 out << " Nbins=" << nn << " " 729 << std::setw(3) << G4BestUnit(v->Energ 803 << std::setw(3) << G4BestUnit(v->Energy(0),"Energy") 730 << " - " 804 << " - " 731 << std::setw(3) << G4BestUnit(v->Energ 805 << std::setw(3) << G4BestUnit(v->Energy(nn),"Energy"); 732 break; 806 break; 733 } 807 } 734 } 808 } 735 } 809 } 736 G4VEmAngularDistribution* an = model->GetA 810 G4VEmAngularDistribution* an = model->GetAngularDistribution(); 737 if(an) { out << " " << an->GetName(); } 811 if(an) { out << " " << an->GetName(); } 738 if(fluoFlag && model->DeexcitationFlag()) 812 if(fluoFlag && model->DeexcitationFlag()) { 739 out << " Fluo"; 813 out << " Fluo"; 740 } 814 } 741 out << G4endl; 815 out << G4endl; 742 auto msc = dynamic_cast<G4VMscModel* << 816 G4VMscModel* msc = dynamic_cast<G4VMscModel*>(model); 743 if(msc != nullptr) msc->DumpParamete 817 if(msc != nullptr) msc->DumpParameters(out); 744 } 818 } 745 } 819 } 746 } 820 } 747 if(1 == nEmModels) { break; } 821 if(1 == nEmModels) { break; } 748 } 822 } 749 if(theCutsNew) { 823 if(theCutsNew) { 750 out << " ===== Limit on energy thresh 824 out << " ===== Limit on energy threshold has been applied " << G4endl; 751 } 825 } 752 } 826 } 753 827 754 //....oooOO0OOooo........oooOO0OOooo........oo 828 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 755 829