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