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