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