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: G4VEnergyLossProcess.cc,v 1.10 2004/03/12 11:42:57 vnivanch Exp $ >> 24 // GEANT4 tag $Name: geant4-06-01 $ >> 25 // 26 // ------------------------------------------- 26 // ------------------------------------------------------------------- 27 // 27 // 28 // GEANT4 Class file 28 // GEANT4 Class file 29 // 29 // 30 // 30 // 31 // File name: G4VEnergyLossProcess 31 // File name: G4VEnergyLossProcess 32 // 32 // 33 // Author: Vladimir Ivanchenko 33 // Author: Vladimir Ivanchenko 34 // 34 // 35 // Creation date: 03.01.2002 35 // Creation date: 03.01.2002 36 // 36 // 37 // Modifications: Vladimir Ivanchenko << 37 // Modifications: 38 // 38 // >> 39 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko) >> 40 // 04-12-02 Minor change in PostStepDoIt (V.Ivanchenko) >> 41 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko) >> 42 // 26-12-02 Secondary production moved to derived classes (V.Ivanchenko) >> 43 // 04-01-03 Fix problem of very small steps for ions (V.Ivanchenko) >> 44 // 20-01-03 Migrade to cut per region (V.Ivanchenko) >> 45 // 24-01-03 Temporarily close a control on usage of couples (V.Ivanchenko) >> 46 // 24-01-03 Make models region aware (V.Ivanchenko) >> 47 // 05-02-03 Fix compilation warnings (V.Ivanchenko) >> 48 // 06-02-03 Add control on tmax in PostStepDoIt (V.Ivanchenko) >> 49 // 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko) >> 50 // 15-02-03 Lambda table can be scaled (V.Ivanchenko) >> 51 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko) >> 52 // 18-02-03 Add control on CutCouple usage (V.Ivanchenko) >> 53 // 26-02-03 Simplify control on GenericIons (V.Ivanchenko) >> 54 // 06-03-03 Control on GenericIons using SubType + update verbose (V.Ivanchenko) >> 55 // 10-03-03 Add Ion registration (V.Ivanchenko) >> 56 // 22-03-03 Add Initialisation of cash (V.Ivanchenko) >> 57 // 26-03-03 Remove finalRange modification (V.Ivanchenko) >> 58 // 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko) >> 59 // 26-04-03 Fix retrieve tables (V.Ivanchenko) >> 60 // 06-05-03 Set defalt finalRange = 1 mm (V.Ivanchenko) >> 61 // 12-05-03 Update range calculations + lowKinEnergy (V.Ivanchenko) >> 62 // 13-05-03 Add calculation of precise range (V.Ivanchenko) >> 63 // 23-05-03 Remove tracking cuts (V.Ivanchenko) >> 64 // 03-06-03 Fix initialisation problem for STD ionisation (V.Ivanchenko) >> 65 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko) >> 66 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko) >> 67 // 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko) >> 68 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko) >> 69 // 21-01-04 Migrade to G4ParticleChangeForLoss (V.Ivanchenko) >> 70 // 27-02-04 Fix problem of loss in low presure gases, cleanup precise range >> 71 // calculation, use functions ForLoss in AlongStepDoIt (V.Ivanchenko) >> 72 // 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko) 39 // 73 // 40 // Class Description: 74 // Class Description: 41 // 75 // 42 // It is the unified energy loss process it ca 76 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 77 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 78 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 79 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 80 // that information on fly. 47 // ------------------------------------------- 81 // ------------------------------------------------------------------- 48 // 82 // 49 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 85 52 #include "G4VEnergyLossProcess.hh" 86 #include "G4VEnergyLossProcess.hh" 53 #include "G4PhysicalConstants.hh" << 54 #include "G4SystemOfUnits.hh" << 55 #include "G4ProcessManager.hh" << 56 #include "G4LossTableManager.hh" 87 #include "G4LossTableManager.hh" 57 #include "G4LossTableBuilder.hh" << 58 #include "G4Step.hh" 88 #include "G4Step.hh" 59 #include "G4ParticleDefinition.hh" 89 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" << 61 #include "G4EmParameters.hh" << 62 #include "G4EmUtility.hh" << 63 #include "G4EmTableUtil.hh" << 64 #include "G4VEmModel.hh" 90 #include "G4VEmModel.hh" 65 #include "G4VEmFluctuationModel.hh" 91 #include "G4VEmFluctuationModel.hh" 66 #include "G4DataVector.hh" 92 #include "G4DataVector.hh" >> 93 #include "G4PhysicsTable.hh" >> 94 #include "G4PhysicsVector.hh" 67 #include "G4PhysicsLogVector.hh" 95 #include "G4PhysicsLogVector.hh" 68 #include "G4VParticleChange.hh" 96 #include "G4VParticleChange.hh" >> 97 #include "G4Gamma.hh" 69 #include "G4Electron.hh" 98 #include "G4Electron.hh" >> 99 #include "G4Proton.hh" >> 100 #include "G4VSubCutoffProcessor.hh" 70 #include "G4ProcessManager.hh" 101 #include "G4ProcessManager.hh" 71 #include "G4UnitsTable.hh" 102 #include "G4UnitsTable.hh" >> 103 #include "G4GenericIon.hh" >> 104 #include "G4ProductionCutsTable.hh" 72 #include "G4Region.hh" 105 #include "G4Region.hh" 73 #include "G4RegionStore.hh" 106 #include "G4RegionStore.hh" 74 #include "G4PhysicsTableHelper.hh" << 75 #include "G4SafetyHelper.hh" << 76 #include "G4EmDataHandler.hh" << 77 #include "G4TransportationManager.hh" << 78 #include "G4VAtomDeexcitation.hh" << 79 #include "G4VSubCutProducer.hh" << 80 #include "G4EmBiasingManager.hh" << 81 #include "G4Log.hh" << 82 #include <iostream> << 83 107 84 //....oooOO0OOooo........oooOO0OOooo........oo 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 109 86 namespace << 110 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, G4ProcessType type): 87 { << 111 G4VContinuousDiscreteProcess(name, type), 88 G4String tnames[7] = << 112 nSCoffRegions(0), 89 {"DEDX","Ionisation","DEDXnr","CSDARange", << 113 idxSCoffRegions(0), 90 } << 114 theDEDXTable(0), 91 << 115 theRangeTableForLoss(0), >> 116 thePreciseRangeTable(0), >> 117 theSecondaryRangeTable(0), >> 118 theInverseRangeTable(0), >> 119 theLambdaTable(0), >> 120 theSubLambdaTable(0), >> 121 theDEDXAtMaxEnergy(0), >> 122 theRangeAtMaxEnergy(0), >> 123 particle(0), >> 124 baseParticle(0), >> 125 secondaryParticle(0), >> 126 currentCouple(0), >> 127 nDEDXBins(90), >> 128 nDEDXBinsForRange(70), >> 129 nLambdaBins(90), >> 130 linLossLimit(0.05), >> 131 minSubRange(0.1), >> 132 defaultRoverRange(0.2), >> 133 defaultIntegralRange(1.0), >> 134 lossFluctuationFlag(true), >> 135 rndmStepFlag(false), >> 136 hasRestProcess(true), >> 137 tablesAreBuilt(false), >> 138 integral(true), >> 139 meanFreePath(true) >> 140 { >> 141 >> 142 lowestKinEnergy = 1.*eV; >> 143 minKinEnergy = 0.1*keV; >> 144 maxKinEnergy = 100.0*GeV; >> 145 maxKinEnergyForRange = 1.0*GeV; 92 146 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con << 147 pParticleChange = &fParticleChange; 94 G4P << 95 G4VContinuousDiscreteProcess(name, type) << 96 { << 97 theParameters = G4EmParameters::Instance(); << 98 SetVerboseLevel(1); << 99 << 100 // low energy limit << 101 lowestKinEnergy = theParameters->LowestElect << 102 << 103 // Size of tables << 104 minKinEnergy = 0.1*CLHEP::keV; << 105 maxKinEnergy = 100.0*CLHEP::TeV; << 106 maxKinEnergyCSDA = 1.0*CLHEP::GeV; << 107 nBins = 84; << 108 nBinsCSDA = 35; << 109 << 110 invLambdaFactor = 1.0/lambdaFactor; << 111 148 112 // default linear loss limit << 149 // default dRoverRange and finalRange 113 finalRange = 1.*CLHEP::mm; << 150 SetStepFunction(defaultIntegralRange, 1.0*mm); >> 151 SetVerboseLevel(0); 114 152 115 // run time objects << 116 pParticleChange = &fParticleChange; << 117 fParticleChange.SetSecondaryWeightByProcess( << 118 modelManager = new G4EmModelManager(); 153 modelManager = new G4EmModelManager(); 119 safetyHelper = G4TransportationManager::GetT << 154 (G4LossTableManager::Instance())->Register(this); 120 ->GetSafetyHelper(); << 155 scoffProcessors.clear(); 121 aGPILSelection = CandidateForSelection; << 156 scoffRegions.clear(); 122 << 123 // initialise model << 124 lManager = G4LossTableManager::Instance(); << 125 lManager->Register(this); << 126 isMaster = lManager->IsMaster(); << 127 << 128 G4LossTableBuilder* bld = lManager->GetTable << 129 theDensityFactor = bld->GetDensityFactors(); << 130 theDensityIdx = bld->GetCoupleIndexes(); << 131 << 132 scTracks.reserve(10); << 133 secParticles.reserve(12); << 134 emModels = new std::vector<G4VEmModel*>; << 135 } 157 } 136 158 137 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 160 139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 161 G4VEnergyLossProcess::~G4VEnergyLossProcess() 140 { 162 { 141 if (isMaster) { << 163 Clear(); 142 if(nullptr == baseParticle) { delete theDa << 164 143 delete theEnergyOfCrossSectionMax; << 165 if (nSCoffRegions) { 144 if(nullptr != fXSpeaks) { << 166 for (G4int i=0; i<nSCoffRegions; i++) { 145 for(auto const & v : *fXSpeaks) { delete << 167 if (scoffProcessors[i]) { 146 delete fXSpeaks; << 168 for (G4int j=i+1; j<nSCoffRegions; j++) { >> 169 if(scoffProcessors[i] == scoffProcessors[j]) scoffProcessors[j] = 0; >> 170 } >> 171 delete scoffProcessors[i]; >> 172 } 147 } 173 } >> 174 scoffProcessors.clear(); >> 175 scoffRegions.clear(); 148 } 176 } 149 delete modelManager; 177 delete modelManager; 150 delete biasManager; << 178 (G4LossTableManager::Instance())->DeRegister(this); 151 delete scoffRegions; << 152 delete emModels; << 153 lManager->DeRegister(this); << 154 } 179 } 155 180 156 //....oooOO0OOooo........oooOO0OOooo........oo 181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 182 158 G4double G4VEnergyLossProcess::MinPrimaryEnerg << 183 void G4VEnergyLossProcess::Clear() 159 << 160 << 161 { 184 { 162 return cut; << 185 if(0 < verboseLevel) { 163 } << 186 G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName() << G4endl; 164 << 187 } 165 //....oooOO0OOooo........oooOO0OOooo........oo << 188 if ( !baseParticle ) { 166 << 189 if(theDEDXTable) theDEDXTable->clearAndDestroy(); 167 void G4VEnergyLossProcess::AddEmModel(G4int or << 190 if(thePreciseRangeTable) thePreciseRangeTable->clearAndDestroy(); 168 G4VEmFlu << 191 if(theRangeTableForLoss) theRangeTableForLoss->clearAndDestroy(); 169 const G4 << 192 if(theInverseRangeTable) theInverseRangeTable->clearAndDestroy(); 170 { << 193 if(theLambdaTable) theLambdaTable->clearAndDestroy(); 171 if(nullptr == ptr) { return; } << 194 if(theSubLambdaTable) theSubLambdaTable->clearAndDestroy(); 172 G4VEmFluctuationModel* afluc = (nullptr == f << 195 if(theDEDXAtMaxEnergy) delete [] theDEDXAtMaxEnergy; 173 modelManager->AddEmModel(order, ptr, afluc, << 196 if(theRangeAtMaxEnergy) delete [] theRangeAtMaxEnergy; 174 ptr->SetParticleChange(pParticleChange, aflu << 175 } << 176 << 177 //....oooOO0OOooo........oooOO0OOooo........oo << 178 << 179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod << 180 { << 181 if(nullptr == ptr) { return; } << 182 if(!emModels->empty()) { << 183 for(auto & em : *emModels) { if(em == ptr) << 184 } 197 } 185 emModels->push_back(ptr); << 186 } << 187 << 188 //....oooOO0OOooo........oooOO0OOooo........oo << 189 198 190 void G4VEnergyLossProcess::SetDynamicMassCharg << 199 theDEDXTable = 0; 191 << 200 thePreciseRangeTable = 0; 192 { << 201 theRangeTableForLoss = 0; 193 massRatio = massratio; << 202 theInverseRangeTable = 0; 194 logMassRatio = G4Log(massRatio); << 203 theSecondaryRangeTable = 0; 195 fFactor = charge2ratio*biasFactor; << 204 theLambdaTable = 0; 196 if(baseMat) { fFactor *= (*theDensityFactor) << 205 theSubLambdaTable = 0; 197 chargeSqRatio = charge2ratio; << 206 theDEDXAtMaxEnergy = 0; 198 reduceFactor = 1.0/(fFactor*massRatio); << 207 theRangeAtMaxEnergy = 0; >> 208 tablesAreBuilt = false; 199 } 209 } 200 210 201 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 202 212 203 void << 213 void G4VEnergyLossProcess::Initialise() 204 G4VEnergyLossProcess::PreparePhysicsTable(cons << 205 { 214 { 206 particle = G4EmTableUtil::CheckIon(this, &pa << 215 if(0 < verboseLevel) { 207 verboseLe << 216 G4cout << "G4VEnergyLossProcess::Initialise() for " 208 << 217 << GetProcessName() 209 if( particle != &part ) { << 218 << " for " << particle->GetParticleName() 210 if(!isIon) { lManager->RegisterExtraPartic << 219 << G4endl; 211 if(1 < verboseLevel) { << 212 G4cout << "### G4VEnergyLossProcess::Pre << 213 << " interrupted for " << GetProc << 214 << part.GetParticleName() << " is << 215 << " spline=" << spline << G4endl << 216 } << 217 return; << 218 } 220 } 219 221 220 tablesAreBuilt = false; << 222 Clear(); 221 if (GetProcessSubType() == fIonisation) { Se << 222 223 223 G4LossTableBuilder* bld = lManager->GetTable << 224 G4double initialCharge = particle->GetPDGCharge(); 224 lManager->PreparePhysicsTable(&part, this); << 225 G4double initialMass = particle->GetPDGMass(); 225 << 226 chargeSquare = initialCharge*initialCharge/(eplus*eplus); 226 // Base particle and set of models can be de << 227 chargeSqRatio = 1.0; 227 InitialiseEnergyLossProcess(particle, basePa << 228 massRatio = 1.0; >> 229 reduceFactor = 1.0; >> 230 >> 231 if(particle->GetProcessManager()->GetAtRestProcessVector()->size()) >> 232 hasRestProcess = true; >> 233 else hasRestProcess = false; 228 234 229 // parameters of the process << 235 if (baseParticle) { 230 if(!actLossFluc) { lossFluctuationFlag = the << 236 massRatio = (baseParticle->GetPDGMass())/initialMass; 231 useCutAsFinalRange = theParameters->UseCutAs << 237 G4double q = initialCharge/baseParticle->GetPDGCharge(); 232 if(!actMinKinEnergy) { minKinEnergy = thePar << 238 chargeSqRatio = q*q; 233 if(!actMaxKinEnergy) { maxKinEnergy = thePar << 239 reduceFactor = 1.0/(chargeSqRatio*massRatio); 234 if(!actBinning) { nBins = theParameters->Num << 240 } 235 maxKinEnergyCSDA = theParameters->MaxEnergyF << 236 nBinsCSDA = theParameters->NumberOfBinsPerDe << 237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi << 238 if(!actLinLossLimit) { linLossLimit = thePar << 239 lambdaFactor = theParameters->LambdaFactor() << 240 invLambdaFactor = 1.0/lambdaFactor; << 241 if(isMaster) { SetVerboseLevel(theParameters << 242 else { SetVerboseLevel(theParameters->Worker << 243 // integral option may be disabled << 244 if(!theParameters->Integral()) { fXSType = f << 245 241 246 theParameters->DefineRegParamForLoss(this); << 242 theCuts = modelManager->Initialise(particle, secondaryParticle, minSubRange, verboseLevel); 247 243 248 fRangeEnergy = 0.0; << 244 // Sub Cutoff Regime 249 245 250 G4double initialCharge = particle->GetPDGCha << 246 idxSCoffRegions.clear(); 251 G4double initialMass = particle->GetPDGMas << 252 247 253 theParameters->FillStepFunction(particle, th << 248 const G4ProductionCutsTable* theCoupleTable= >> 249 G4ProductionCutsTable::GetProductionCutsTable(); >> 250 size_t numOfCouples = theCoupleTable->GetTableSize(); 254 251 255 // parameters for scaling from the base part << 252 if (nSCoffRegions) { 256 if (nullptr != baseParticle) { << 253 const G4DataVector* theSubCuts = modelManager->SubCutoff(); 257 massRatio = (baseParticle->GetPDGMass() << 258 logMassRatio = G4Log(massRatio); << 259 G4double q = initialCharge/baseParticle->G << 260 chargeSqRatio = q*q; << 261 if(chargeSqRatio > 0.0) { reduceFactor = 1 << 262 } << 263 lowestKinEnergy = (initialMass < CLHEP::MeV) << 264 ? theParameters->LowestElectronEnergy() << 265 : theParameters->LowestMuHadEnergy(); << 266 << 267 // Tables preparation << 268 if (isMaster && nullptr == baseParticle) { << 269 if(nullptr == theData) { theData = new G4E << 270 << 271 if(nullptr != theDEDXTable && isIonisation << 272 if(nullptr != theIonisationTable && theD << 273 theData->CleanTable(0); << 274 theDEDXTable = theIonisationTable; << 275 theIonisationTable = nullptr; << 276 } << 277 } << 278 << 279 theDEDXTable = theData->MakeTable(theDEDXT << 280 bld->InitialiseBaseMaterials(theDEDXTable) << 281 theData->UpdateTable(theIonisationTable, 1 << 282 << 283 if (theParameters->BuildCSDARange()) { << 284 theDEDXunRestrictedTable = theData->Make << 285 if(isIonisation) { theCSDARangeTable = t << 286 } << 287 254 288 theLambdaTable = theData->MakeTable(4); << 255 for (G4int i=0; i<nSCoffRegions; i++) { 289 if(isIonisation) { << 256 scoffProcessors[i]->Initialise(particle, secondaryParticle, theCuts, theSubCuts); 290 theRangeTableForLoss = theData->MakeTabl << 291 theInverseRangeTable = theData->MakeTabl << 292 } 257 } 293 } << 258 for (size_t j=0; j<numOfCouples; j++) { 294 259 295 // forced biasing << 260 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j); 296 if(nullptr != biasManager) { << 261 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); 297 biasManager->Initialise(part,GetProcessNam << 262 G4int reg = nSCoffRegions; 298 biasFlag = false; << 263 do {reg--;} while (reg && pcuts != (scoffRegions[reg]->GetProductionCuts())); 299 } << 264 idxSCoffRegions.push_back(reg); 300 baseMat = bld->GetBaseMaterialFlag(); << 301 numberOfModels = modelManager->NumberOfModel << 302 currentModel = modelManager->GetModel(0); << 303 G4EmTableUtil::UpdateModels(this, modelManag << 304 numberOfModels, << 305 mainSecondaries, << 306 theParameters->U << 307 theCuts = modelManager->Initialise(particle, << 308 verboseLe << 309 // subcut processor << 310 if(isIonisation) { << 311 subcutProducer = lManager->SubCutProducer( << 312 } << 313 if(1 == nSCoffRegions) { << 314 if((*scoffRegions)[0]->GetName() == "Defau << 315 delete scoffRegions; << 316 scoffRegions = nullptr; << 317 nSCoffRegions = 0; << 318 } 265 } 319 } 266 } 320 267 321 if(1 < verboseLevel) { << 268 if (0 < verboseLevel) { 322 G4cout << "G4VEnergyLossProcess::PrepearPh << 269 G4cout << "G4VEnergyLossProcess::Initialise() is done " 323 << " for " << GetProcessName() << " << 270 << " chargeSqRatio= " << chargeSqRatio 324 << " isIon= " << isIon << " spline= << 325 if(baseParticle) { << 326 G4cout << "; base: " << baseParticle->Ge << 327 } << 328 G4cout << G4endl; << 329 G4cout << " chargeSqRatio= " << chargeSqRa << 330 << " massRatio= " << massRatio 271 << " massRatio= " << massRatio 331 << " reduceFactor= " << reduceFacto 272 << " reduceFactor= " << reduceFactor << G4endl; 332 if (nSCoffRegions > 0) { << 273 if (nSCoffRegions) { 333 G4cout << " SubCut secondary production << 274 G4cout << " SubCutoff Regime is ON for regions: " << G4endl; 334 for (G4int i=0; i<nSCoffRegions; ++i) { << 275 for (G4int i=0; i<nSCoffRegions; i++) { 335 const G4Region* r = (*scoffRegions)[i] << 276 const G4Region* r = scoffRegions[i]; 336 G4cout << " " << r->GetName( << 277 G4cout << " " << r->GetName() << G4endl; 337 } 278 } 338 } else if(nullptr != subcutProducer) { << 339 G4cout << " SubCut secondary production << 340 } 279 } 341 } 280 } 342 } 281 } 343 282 344 //....oooOO0OOooo........oooOO0OOooo........oo 283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 345 284 346 void G4VEnergyLossProcess::BuildPhysicsTable(c 285 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 347 { 286 { 348 if(1 < verboseLevel) { << 287 currentCouple = 0; >> 288 preStepLambda = 0.0; >> 289 if(0 < verboseLevel) { >> 290 G4cout << "========================================================" << G4endl; 349 G4cout << "### G4VEnergyLossProcess::Build 291 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for " 350 << GetProcessName() 292 << GetProcessName() 351 << " and particle " << part.GetPart 293 << " and particle " << part.GetParticleName() 352 << "; the first particle " << parti << 294 << G4endl; 353 if(baseParticle) { << 354 G4cout << "; base: " << baseParticle->Ge << 355 } << 356 G4cout << G4endl; << 357 G4cout << " TablesAreBuilt= " << tables << 358 << " spline=" << spline << " ptr: " << 359 } << 360 << 361 if(&part == particle) { << 362 if(isMaster) { << 363 lManager->BuildPhysicsTable(particle, th << 364 << 365 } else { << 366 const auto masterProcess = << 367 static_cast<const G4VEnergyLossProcess << 368 << 369 numberOfModels = modelManager->NumberOfM << 370 G4EmTableUtil::BuildLocalElossProcess(th << 371 pa << 372 tablesAreBuilt = true; << 373 baseMat = masterProcess->UseBaseMaterial << 374 lManager->LocalPhysicsTables(particle, t << 375 } << 376 << 377 // needs to be done only once << 378 safetyHelper->InitialiseHelper(); << 379 } << 380 // Added tracking cut to avoid tracking arti << 381 // and identified deexcitation flag << 382 if(isIonisation) { << 383 atomDeexcitation = lManager->AtomDeexcitat << 384 if(nullptr != atomDeexcitation) { << 385 if(atomDeexcitation->IsPIXEActive()) { u << 386 } << 387 } 295 } 388 296 389 // protection against double printout << 297 if (part.GetParticleName() != "GenericIon" && 390 if(theParameters->IsPrintLocked()) { return; << 298 part.GetParticleType() == "nucleus" && >> 299 part.GetParticleSubType() == "generic") >> 300 { >> 301 (G4LossTableManager::Instance())->RegisterIon(&part, this); >> 302 /* >> 303 G4cout << part.GetProcessManager() << " " >> 304 << (G4GenericIon::GenericIon())->GetProcessManager() >> 305 << G4endl; >> 306 */ >> 307 return; >> 308 } 391 309 392 // explicitly defined printout by particle n << 310 // Are particle defined? 393 G4String num = part.GetParticleName(); << 311 if( !particle ) { 394 if(1 < verboseLevel || << 312 particle = ∂ 395 (0 < verboseLevel && (num == "e-" || << 313 baseParticle = DefineBaseParticle(particle); 396 num == "e+" || n << 397 num == "mu-" || n << 398 num == "pi+" || n << 399 num == "kaon+" || n << 400 num == "alpha" || n << 401 num == "GenericIon" << 402 StreamInfo(G4cout, part); << 403 } 314 } 404 if(1 < verboseLevel) { << 315 405 G4cout << "### G4VEnergyLossProcess::Build << 316 // Recalculation is needed because cuts were changed or recalculation is forced 406 << GetProcessName() << 317 G4LossTableManager* lManager = G4LossTableManager::Instance(); 407 << " and particle " << part.GetPart << 318 if ( lManager->IsRecalcNeeded(particle) ) { 408 if(isIonisation) { G4cout << " isIonisati << 319 409 G4cout << " baseMat=" << baseMat << G4endl << 320 // It is responsability of the G4LossTables to build DEDX and range tables >> 321 lManager->BuildPhysicsTable(particle); >> 322 >> 323 if(!baseParticle) PrintInfoDefinition(); >> 324 >> 325 if(0 < verboseLevel) { >> 326 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for " >> 327 << GetProcessName() >> 328 << " and particle " << part.GetParticleName() >> 329 << G4endl; >> 330 } 410 } 331 } 411 } 332 } 412 333 413 //....oooOO0OOooo........oooOO0OOooo........oo 334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 414 335 415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED << 336 void G4VEnergyLossProcess::SetParticles(const G4ParticleDefinition* p1, >> 337 const G4ParticleDefinition* p2) 416 { 338 { 417 G4PhysicsTable* table = nullptr; << 339 particle = p1; 418 G4double emax = maxKinEnergy; << 340 baseParticle = p2; 419 G4int bin = nBins; << 341 G4bool yes = true; 420 << 342 if(particle && baseParticle) { 421 if(fTotal == tType) { << 343 if((particle->GetPDGMass() < MeV && baseParticle->GetPDGMass() > MeV) || 422 emax = maxKinEnergyCSDA; << 344 (particle->GetPDGMass() > MeV && baseParticle->GetPDGMass() < MeV)) yes = false; 423 bin = nBinsCSDA; << 345 424 table = theDEDXunRestrictedTable; << 346 } 425 } else if(fRestricted == tType) { << 347 if(!yes) { 426 table = theDEDXTable; << 348 G4cout << "Warning in G4VEnergyLossProcess::SetParticle: " >> 349 << particle->GetParticleName() >> 350 << " losses cannot be obtained from " >> 351 << baseParticle->GetParticleName() >> 352 << " losses" << G4endl; 427 } else { 353 } else { 428 G4cout << "G4VEnergyLossProcess::BuildDEDX << 354 DefineBaseParticle(p1); 429 << tType << G4endl; << 430 } << 431 if(1 < verboseLevel) { << 432 G4cout << "G4VEnergyLossProcess::BuildDEDX << 433 << " for " << GetProcessName() << 434 << " and " << particle->GetParticle << 435 << "spline=" << spline << G4endl; << 436 } << 437 if(nullptr == table) { return table; } << 438 << 439 G4LossTableBuilder* bld = lManager->GetTable << 440 G4EmTableUtil::BuildDEDXTable(this, particle << 441 table, minKinE << 442 verboseLevel, << 443 return table; << 444 } << 445 << 446 //....oooOO0OOooo........oooOO0OOooo........oo << 447 << 448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam << 449 { << 450 if(nullptr == theLambdaTable) { return theLa << 451 << 452 G4double scale = theParameters->MaxKinEnergy << 453 G4int nbin = << 454 theParameters->NumberOfBinsPerDecade()*G4l << 455 scale = nbin/G4Log(scale); << 456 << 457 G4LossTableBuilder* bld = lManager->GetTable << 458 G4EmTableUtil::BuildLambdaTable(this, partic << 459 bld, theLamb << 460 minKinEnergy << 461 verboseLevel << 462 return theLambdaTable; << 463 } << 464 << 465 //....oooOO0OOooo........oooOO0OOooo........oo << 466 << 467 void G4VEnergyLossProcess::StreamInfo(std::ost << 468 const G4ParticleDefinition& pa << 469 { << 470 G4String indent = (rst ? " " : ""); << 471 out << std::setprecision(6); << 472 out << G4endl << indent << GetProcessName() << 473 if (!rst) out << " for " << part.GetParticle << 474 out << " XStype:" << fXSType << 475 << " SubType=" << GetProcessSubType() < << 476 << " dE/dx and range tables from " << 477 << G4BestUnit(minKinEnergy,"Energy") << 478 << " to " << G4BestUnit(maxKinEnergy,"En << 479 << " in " << nBins << " bins" << G4endl << 480 << " Lambda tables from threshold t << 481 << G4BestUnit(maxKinEnergy,"Energy") << 482 << ", " << theParameters->NumberOfBinsPe << 483 << " bins/decade, spline: " << spline << 484 << G4endl; << 485 if(nullptr != theRangeTableForLoss && isIoni << 486 out << " StepFunction=(" << dRoverRan << 487 << finalRange/mm << " mm)" << 488 << ", integ: " << fXSType << 489 << ", fluct: " << lossFluctuationFlag << 490 << ", linLossLim= " << linLossLimit << 491 << G4endl; << 492 } << 493 StreamProcessInfo(out); << 494 modelManager->DumpModelList(out, verboseLeve << 495 if(nullptr != theCSDARangeTable && isIonisat << 496 out << " CSDA range table up" << 497 << " to " << G4BestUnit(maxKinEnergyCS << 498 << " in " << nBinsCSDA << " bins" << G << 499 } << 500 if(nSCoffRegions>0 && isIonisation) { << 501 out << " Subcutoff sampling in " << n << 502 << " regions" << G4endl; << 503 } << 504 if(2 < verboseLevel) { << 505 for(std::size_t i=0; i<7; ++i) { << 506 auto ta = theData->Table(i); << 507 out << " " << tnames[i] << " addres << 508 if(nullptr != ta) { out << *ta << G4endl << 509 } << 510 } 355 } 511 } 356 } 512 357 513 //....oooOO0OOooo........oooOO0OOooo........oo 358 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 514 359 515 void G4VEnergyLossProcess::ActivateSubCutoff(c << 360 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p, G4VEmFluctuationModel* fluc, >> 361 const G4Region* region) 516 { 362 { 517 if(nullptr == scoffRegions) { << 363 modelManager->AddEmModel(order, p, fluc, region); 518 scoffRegions = new std::vector<const G4Reg << 519 } << 520 // the region is in the list << 521 if(!scoffRegions->empty()) { << 522 for (auto & reg : *scoffRegions) { << 523 if (reg == r) { return; } << 524 } << 525 } << 526 // new region << 527 scoffRegions->push_back(r); << 528 ++nSCoffRegions; << 529 } 364 } 530 365 531 //....oooOO0OOooo........oooOO0OOooo........oo 366 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 532 367 533 G4bool G4VEnergyLossProcess::IsRegionForCubcut << 368 void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam, G4double emin, G4double emax) 534 { 369 { 535 if(0 == nSCoffRegions) { return true; } << 370 modelManager->UpdateEmModel(nam, emin, emax); 536 const G4Region* r = aTrack.GetVolume()->GetL << 537 for(auto & reg : *scoffRegions) { << 538 if(r == reg) { return true; } << 539 } << 540 return false; << 541 } 371 } 542 372 543 //....oooOO0OOooo........oooOO0OOooo........oo 373 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 544 374 545 void G4VEnergyLossProcess::StartTracking(G4Tra << 375 void G4VEnergyLossProcess::AddSubCutoffProcessor(G4VSubCutoffProcessor* p, >> 376 const G4Region* r) 546 { 377 { 547 // reset parameters for the new track << 378 if( !p ) { 548 theNumberOfInteractionLengthLeft = -1.0; << 379 G4cout << "G4VEnergyLossProcess::AddSubCutoffProcessor WARNING: no SubCutoffProcessor defined." << G4endl; 549 mfpKinEnergy = DBL_MAX; << 380 return; 550 preStepLambda = 0.0; << 381 } 551 currentCouple = nullptr; << 382 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 552 << 383 if (!r) r = regionStore->GetRegion("DefaultRegionForTheWorld", false); 553 // reset ion << 384 if (nSCoffRegions) { 554 if(isIon) { << 385 for (G4int i=0; i<nSCoffRegions; i++) { 555 const G4double newmass = track->GetDefinit << 386 if (r == scoffRegions[i]) { 556 massRatio = (nullptr == baseParticle) ? CL << 387 if ( scoffProcessors[i] ) delete scoffProcessors[i]; 557 : baseParticle->GetPDGMass()/newmass; << 388 scoffProcessors[i] = p; 558 logMassRatio = G4Log(massRatio); << 389 return; 559 } << 390 } 560 // forced biasing only for primary particles << 561 if(nullptr != biasManager) { << 562 if(0 == track->GetParentID()) { << 563 biasFlag = true; << 564 biasManager->ResetForcedInteraction(); << 565 } 391 } 566 } 392 } >> 393 scoffProcessors.push_back(p); >> 394 scoffRegions.push_back(r); >> 395 nSCoffRegions++; 567 } 396 } 568 397 569 //....oooOO0OOooo........oooOO0OOooo........oo 398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 570 399 571 G4double G4VEnergyLossProcess::AlongStepGetPhy << 400 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable() 572 const G4Track& tr << 401 { 573 G4GPILSelection* << 402 574 { << 403 if(0 < verboseLevel) { 575 G4double x = DBL_MAX; << 404 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() for " 576 *selection = aGPILSelection; << 405 << GetProcessName() 577 if(isIonisation && currentModel->IsActive(pr << 406 << " and particle " << particle->GetParticleName() 578 GetScaledRangeForScaledEnergy(preStepScale << 407 << G4endl; 579 x = (useCutAsFinalRange) ? std::min(finalR << 580 currentCouple->GetProductionCuts()->GetP << 581 x = (fRange > x) ? fRange*dRoverRange + x* << 582 : fRange; << 583 /* << 584 G4cout<<"AlongStepGPIL: " << GetProcessN << 585 << " fRange=" << fRange << " finR=" << finR << 586 */ << 587 } 408 } 588 return x; << 409 >> 410 // vectors to provide continues dE/dx >> 411 G4DataVector factor; >> 412 G4DataVector dedxLow; >> 413 G4DataVector dedxHigh; >> 414 >> 415 // Access to materials >> 416 const G4ProductionCutsTable* theCoupleTable= >> 417 G4ProductionCutsTable::GetProductionCutsTable(); >> 418 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 419 >> 420 G4PhysicsTable* theTable = new G4PhysicsTable(numOfCouples); >> 421 >> 422 if(0 < verboseLevel) { >> 423 G4cout << numOfCouples << " materials" >> 424 << " minKinEnergy= " << minKinEnergy >> 425 << " maxKinEnergy= " << maxKinEnergy >> 426 << G4endl; >> 427 } >> 428 >> 429 for(size_t i=0; i<numOfCouples; i++) { >> 430 >> 431 // create physics vector and fill it >> 432 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); >> 433 G4PhysicsVector* aVector = DEDXPhysicsVector(couple); >> 434 modelManager->FillDEDXVector(aVector, couple); >> 435 >> 436 // Insert vector for this material into the table >> 437 theTable->insert(aVector) ; >> 438 } >> 439 >> 440 if(0 < verboseLevel) { >> 441 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for " >> 442 << particle->GetParticleName() >> 443 << G4endl; >> 444 if(2 < verboseLevel) { >> 445 G4cout << *theTable << G4endl; >> 446 } >> 447 } >> 448 >> 449 return theTable; 589 } 450 } 590 451 591 //....oooOO0OOooo........oooOO0OOooo........oo 452 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 592 453 593 G4double G4VEnergyLossProcess::PostStepGetPhys << 454 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTableForPreciseRange() 594 const G4Track& tr << 455 { 595 G4double previo << 456 596 G4ForceCondition* << 457 if(0 < verboseLevel) { 597 { << 458 G4cout << "G4VEnergyLossProcess::BuildDEDXTableForPreciseRange() for " 598 // condition is set to "Not Forced" << 459 << GetProcessName() 599 *condition = NotForced; << 460 << " and particle " << particle->GetParticleName() 600 G4double x = DBL_MAX; << 461 << G4endl; 601 << 462 } 602 // initialisation of material, mass, charge, << 463 603 // at the beginning of the step << 464 // vectors to provide continues dE/dx 604 DefineMaterial(track.GetMaterialCutsCouple() << 465 G4DataVector factor; 605 preStepKinEnergy = track.GetKineticEne << 466 G4DataVector dedxLow; 606 preStepScaledEnergy = preStepKinEnergy*ma << 467 G4DataVector dedxHigh; 607 SelectModel(preStepScaledEnergy); << 468 608 << 469 // Access to materials 609 if(!currentModel->IsActive(preStepScaledEner << 470 const G4ProductionCutsTable* theCoupleTable= 610 theNumberOfInteractionLengthLeft = -1.0; << 471 G4ProductionCutsTable::GetProductionCutsTable(); 611 mfpKinEnergy = DBL_MAX; << 472 size_t numOfCouples = theCoupleTable->GetTableSize(); 612 preStepLambda = 0.0; << 473 613 currentInteractionLength = DBL_MAX; << 474 G4PhysicsTable* theTable = new G4PhysicsTable(numOfCouples); 614 return x; << 475 615 } << 476 if(0 < verboseLevel) { 616 << 477 G4cout << numOfCouples << " materials" 617 // change effective charge of a charged part << 478 << " minKinEnergy= " << minKinEnergy 618 if(isIon) { << 479 << " maxKinEnergy= " << maxKinEnergy 619 const G4double q2 = currentModel->ChargeSq << 480 << G4endl; 620 fFactor = q2*biasFactor; << 481 } 621 if(baseMat) { fFactor *= (*theDensityFacto << 482 622 reduceFactor = 1.0/(fFactor*massRatio); << 483 for(size_t i=0; i<numOfCouples; i++) { 623 if (lossFluctuationFlag) { << 484 624 auto fluc = currentModel->GetModelOfFluc << 485 // create physics vector and fill it 625 fluc->SetParticleAndCharge(track.GetDefi << 486 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); 626 } << 487 G4PhysicsVector* aVector = DEDXPhysicsVectorForPreciseRange(couple); >> 488 modelManager->FillDEDXVectorForPreciseRange(aVector, couple); >> 489 >> 490 // Insert vector for this material into the table >> 491 theTable->insert(aVector) ; 627 } 492 } 628 493 629 // forced biasing only for primary particles << 494 if(0 < verboseLevel) { 630 if(biasManager) { << 495 G4cout << "G4VEnergyLossProcess::BuildDEDXTableForPreciseRange(): table is built for " 631 if(0 == track.GetParentID() && biasFlag && << 496 << particle->GetParticleName() 632 biasManager->ForcedInteractionRegion((G << 497 << G4endl; 633 return biasManager->GetStepLimit((G4int) << 498 if(2 < verboseLevel) { >> 499 G4cout << *theTable << G4endl; 634 } 500 } 635 } 501 } 636 502 637 ComputeLambdaForScaledEnergy(preStepScaledEn << 503 return theTable; 638 << 504 } 639 // zero cross section << 640 if(preStepLambda <= 0.0) { << 641 theNumberOfInteractionLengthLeft = -1.0; << 642 currentInteractionLength = DBL_MAX; << 643 } else { << 644 505 645 // non-zero cross section << 506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 646 if (theNumberOfInteractionLengthLeft < 0.0 << 647 507 648 // beggining of tracking (or just after << 508 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable() 649 theNumberOfInteractionLengthLeft = -G4Lo << 509 { 650 theInitialNumberOfInteractionLength = th << 651 510 652 } else if(currentInteractionLength < DBL_M << 511 if(0 < verboseLevel) { >> 512 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() for process " >> 513 << GetProcessName() << " and particle " >> 514 << particle->GetParticleName() >> 515 << G4endl; >> 516 } >> 517 >> 518 // Access to materials >> 519 const G4ProductionCutsTable* theCoupleTable= >> 520 G4ProductionCutsTable::GetProductionCutsTable(); >> 521 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 522 >> 523 G4PhysicsTable* theTable = new G4PhysicsTable(numOfCouples); 653 524 654 // subtract NumberOfInteractionLengthLef << 525 for(size_t i=0; i<numOfCouples; i++) { 655 theNumberOfInteractionLengthLeft -= << 656 previousStepSize/currentInteractionLen << 657 526 658 theNumberOfInteractionLengthLeft = << 527 // create physics vector and fill it 659 std::max(theNumberOfInteractionLengthL << 528 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); >> 529 G4PhysicsVector* aVector = LambdaPhysicsVector(couple); >> 530 modelManager->FillLambdaVector(aVector, couple); >> 531 >> 532 // Insert vector for this material into the table >> 533 theTable->insert(aVector) ; >> 534 } >> 535 >> 536 if(0 < verboseLevel) { >> 537 G4cout << "Lambda table is built for " >> 538 << particle->GetParticleName() >> 539 << G4endl; >> 540 if(2 < verboseLevel) { >> 541 G4cout << *theTable << G4endl; 660 } 542 } >> 543 } 661 544 662 // new mean free path and step limit << 545 return theTable; 663 currentInteractionLength = 1.0/preStepLamb << 664 x = theNumberOfInteractionLengthLeft * cur << 665 } << 666 #ifdef G4VERBOSE << 667 if (verboseLevel>2) { << 668 G4cout << "G4VEnergyLossProcess::PostStepG << 669 G4cout << "[ " << GetProcessName() << "]" << 670 G4cout << " for " << track.GetDefinition() << 671 << " in Material " << currentMate << 672 << " Ekin(MeV)= " << preStepKinEner << 673 << " track material: " << track.Get << 674 <<G4endl; << 675 G4cout << "MeanFreePath = " << currentInte << 676 << "InteractionLength= " << x/cm << << 677 } << 678 #endif << 679 return x; << 680 } 546 } 681 547 682 //....oooOO0OOooo........oooOO0OOooo........oo 548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 683 549 684 void << 550 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaSubTable() 685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 551 { 686 { << 552 if(0 < verboseLevel) { 687 // cross section increased with energy << 553 G4cout << "G4VEnergyLossProcess::BuildLambdaSubTable() for process " 688 if(fXSType == fEmIncreasing) { << 554 << GetProcessName() << " and particle " 689 if(e*invLambdaFactor < mfpKinEnergy) { << 555 << particle->GetParticleName() << G4endl; 690 preStepLambda = GetLambdaForScaledEnergy << 556 } 691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 692 } << 693 557 694 // cross section has one peak << 558 // Access to materials 695 } else if(fXSType == fEmOnePeak) { << 559 const G4ProductionCutsTable* theCoupleTable= 696 const G4double epeak = (*theEnergyOfCrossS << 560 G4ProductionCutsTable::GetProductionCutsTable(); 697 if(e <= epeak) { << 561 size_t numOfCouples = theCoupleTable->GetTableSize(); 698 if(e*invLambdaFactor < mfpKinEnergy) { << 562 G4PhysicsTable* theTable = new G4PhysicsTable(numOfCouples); 699 preStepLambda = GetLambdaForScaledEner << 700 mfpKinEnergy = (preStepLambda > 0.0) ? << 701 } << 702 } else if(e < mfpKinEnergy) { << 703 const G4double e1 = std::max(epeak, e*la << 704 mfpKinEnergy = e1; << 705 preStepLambda = GetLambdaForScaledEnergy << 706 } << 707 563 708 // cross section has more than one peaks << 564 for(size_t i=0; i<numOfCouples; i++) { 709 } else if(fXSType == fEmTwoPeaks) { << 565 710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple << 566 // create physics vector and fill it 711 const G4double e1peak = xs->e1peak; << 567 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); 712 << 568 G4PhysicsVector* aVector = SubLambdaPhysicsVector(couple); 713 // below the 1st peak << 569 modelManager->FillSubLambdaVector(aVector, couple); 714 if(e <= e1peak) { << 570 715 if(e*invLambdaFactor < mfpKinEnergy) { << 571 // Insert vector for this material into the table 716 preStepLambda = GetLambdaForScaledEner << 572 theTable->insert(aVector) ; 717 mfpKinEnergy = (preStepLambda > 0.0) ? << 718 } << 719 return; << 720 } << 721 const G4double e1deep = xs->e1deep; << 722 // above the 1st peak, below the deep << 723 if(e <= e1deep) { << 724 if(mfpKinEnergy >= e1deep || e <= mfpKin << 725 const G4double e1 = std::max(e1peak, e << 726 mfpKinEnergy = e1; << 727 preStepLambda = GetLambdaForScaledEner << 728 } << 729 return; << 730 } << 731 const G4double e2peak = xs->e2peak; << 732 // above the deep, below 2nd peak << 733 if(e <= e2peak) { << 734 if(e*invLambdaFactor < mfpKinEnergy) { << 735 mfpKinEnergy = e; << 736 preStepLambda = GetLambdaForScaledEner << 737 } << 738 return; << 739 } << 740 const G4double e2deep = xs->e2deep; << 741 // above the 2nd peak, below the deep << 742 if(e <= e2deep) { << 743 if(mfpKinEnergy >= e2deep || e <= mfpKin << 744 const G4double e1 = std::max(e2peak, e << 745 mfpKinEnergy = e1; << 746 preStepLambda = GetLambdaForScaledEner << 747 } << 748 return; << 749 } << 750 const G4double e3peak = xs->e3peak; << 751 // above the deep, below 3d peak << 752 if(e <= e3peak) { << 753 if(e*invLambdaFactor < mfpKinEnergy) { << 754 mfpKinEnergy = e; << 755 preStepLambda = GetLambdaForScaledEner << 756 } << 757 return; << 758 } << 759 // above 3d peak << 760 if(e <= mfpKinEnergy) { << 761 const G4double e1 = std::max(e3peak, e*l << 762 mfpKinEnergy = e1; << 763 preStepLambda = GetLambdaForScaledEnergy << 764 } << 765 // integral method is not used << 766 } else { << 767 preStepLambda = GetLambdaForScaledEnergy(e << 768 } 573 } >> 574 >> 575 if(0 < verboseLevel) { >> 576 G4cout << "Table is built for " >> 577 << particle->GetParticleName() >> 578 << G4endl; >> 579 } >> 580 >> 581 return theTable; 769 } 582 } 770 583 >> 584 771 //....oooOO0OOooo........oooOO0OOooo........oo 585 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 772 586 773 G4VParticleChange* G4VEnergyLossProcess::Along 587 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track, 774 588 const G4Step& step) 775 { 589 { 776 fParticleChange.InitializeForAlongStep(track 590 fParticleChange.InitializeForAlongStep(track); 777 // The process has range table - calculate e 591 // The process has range table - calculate energy loss 778 if(!isIonisation || !currentModel->IsActive( << 592 if(!theRangeTableForLoss) return &fParticleChange; 779 return &fParticleChange; << 780 } << 781 593 >> 594 // Get the actual (true) Step length 782 G4double length = step.GetStepLength(); 595 G4double length = step.GetStepLength(); 783 G4double eloss = 0.0; 596 G4double eloss = 0.0; 784 << 597 // G4bool b; 785 /* << 598 >> 599 /* 786 if(-1 < verboseLevel) { 600 if(-1 < verboseLevel) { 787 const G4ParticleDefinition* d = track.GetP << 601 const G4ParticleDefinition* d = track.GetDefinition(); 788 G4cout << "AlongStepDoIt for " 602 G4cout << "AlongStepDoIt for " 789 << GetProcessName() << " and partic << 603 << GetProcessName() << " and particle " 790 << " eScaled(MeV)=" << preStepScal << 604 << d->GetParticleName() 791 << " range(mm)=" << fRange/mm << " << 605 << " eScaled(MeV)= " << preStepScaledEnergy/MeV 792 << " rf=" << reduceFactor << " q^ << 606 << " slim(mm)= " << fRange/mm 793 << " md=" << d->GetPDGMass() << " << 607 << " s(mm)= " << length/mm 794 << " " << track.GetMaterial()->Get << 608 << " q^2= " << chargeSqRatio 795 } << 609 << " md= " << d->GetPDGMass() 796 */ << 610 << G4endl; 797 const G4DynamicParticle* dynParticle = track << 798 << 799 // define new weight for primary and seconda << 800 G4double weight = fParticleChange.GetParentW << 801 if(weightFlag) { << 802 weight /= biasFactor; << 803 fParticleChange.ProposeWeight(weight); << 804 } 611 } 805 << 612 */ 806 // stopping, check actual range and kinetic << 613 // stopping 807 if (length >= fRange || preStepKinEnergy <= << 614 if (length >= fRange) { 808 eloss = preStepKinEnergy; 615 eloss = preStepKinEnergy; 809 if (useDeexcitation) { << 810 atomDeexcitation->AlongStepDeexcitation( << 811 << 812 if(scTracks.size() > 0) { FillSecondarie << 813 eloss = std::max(eloss, 0.0); << 814 } << 815 fParticleChange.SetProposedKineticEnergy(0 << 816 fParticleChange.ProposeLocalEnergyDeposit( << 817 return &fParticleChange; << 818 } << 819 // zero step length with non-zero range << 820 if(length <= 0.0) { return &fParticleChange; << 821 616 822 // Short step 617 // Short step 823 eloss = length*GetDEDXForScaledEnergy(preSte << 618 } else if( length <= linLossLimit * fRange ) { 824 LogSca << 619 eloss = GetDEDXForLoss(preStepKinEnergy)*length; 825 /* << 620 826 G4cout << "##### Short STEP: eloss= " << elo << 827 << " Escaled=" << preStepScaledEnergy << 828 << " R=" << fRange << 829 << " L=" << length << 830 << " fFactor=" << fFactor << " minE=" << mi << 831 << " idxBase=" << basedCoupleIndex << G4end << 832 */ << 833 // Long step 621 // Long step 834 if(eloss > preStepKinEnergy*linLossLimit) { << 622 } else { >> 623 G4double x = (GetRangeForLoss(preStepKinEnergy) - length)/reduceFactor; >> 624 G4double postStepScaledEnergy = GetKineticEnergyForLoss(x); >> 625 >> 626 if (preStepScaledEnergy > postStepScaledEnergy) >> 627 eloss = (preStepScaledEnergy - postStepScaledEnergy)/massRatio; >> 628 else >> 629 eloss = GetDEDXForLoss(preStepKinEnergy)*length; 835 630 836 const G4double x = (fRange - length)/reduc << 837 const G4double de = preStepKinEnergy - Sca << 838 if(de > 0.0) { eloss = de; } << 839 /* 631 /* 840 if(-1 < verboseLevel) << 632 if(-1 < verboseLevel) { 841 G4cout << " Long STEP: rPre(mm)=" << 633 G4bool b; 842 << GetScaledRangeForScaledEnergy( << 634 G4cout << "fRange(mm)= " << fRange/mm 843 << " x(mm)=" << x/mm << 635 << " xPost(mm)= " << x/mm 844 << " eloss(MeV)=" << eloss/MeV << 636 << " ePre(MeV)= " << preStepScaledEnergy/MeV 845 << " rFactor=" << reduceFactor << 637 << " ePost(MeV)= " << postStepScaledEnergy/MeV 846 << " massRatio=" << massRatio << 638 << " eloss(MeV)= " << eloss/MeV >> 639 << " eloss0(MeV)= " << (((*theDEDXTable)[currentMaterialIndex])-> >> 640 GetValue(preStepScaledEnergy, b))*length*chargeSqRatio 847 << G4endl; 641 << G4endl; >> 642 } 848 */ 643 */ >> 644 849 } 645 } 850 646 >> 647 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); >> 648 G4double tmax = MaxSecondaryEnergy(dynParticle); >> 649 tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]); >> 650 851 /* 651 /* 852 if(-1 < verboseLevel ) { << 652 G4double eloss0 = eloss; 853 G4cout << "Before fluct: eloss(MeV)= " << << 653 if(-1 < verboseLevel) { >> 654 G4bool b; >> 655 //G4cout << *theDEDXTable << G4endl; >> 656 G4cout << "eloss(MeV)= " << eloss/MeV >> 657 << " eloss0(MeV)= " << (((*theDEDXTable)[currentMaterialIndex])-> >> 658 GetValue(preStepScaledEnergy, b))*length*chargeSqRatio >> 659 << " r0(mm)= " << (((*theRangeTableForLoss)[currentMaterialIndex])-> >> 660 GetValue(preStepScaledEnergy, b)) >> 661 << " tmax= " << tmax 854 << " e-eloss= " << preStepKinEnergy 662 << " e-eloss= " << preStepKinEnergy-eloss 855 << " step(mm)= " << length/mm << " << 663 // << " mat= " << currentMaterial 856 << " fluct= " << lossFluctuationFla << 664 // << " matIdx= " << currentMaterialIndex >> 665 // << " preCouple= " << (step.GetPreStepPoint())->GetMaterialCutsCouple() >> 666 // << " postCouple= " << (step.GetPostStepPoint())->GetMaterialCutsCouple() >> 667 // << " rangeTable= " << theRangeTableForLoss >> 668 << G4endl; 857 } 669 } 858 */ 670 */ 859 671 860 const G4double cut = (*theCuts)[currentCoupl << 672 // Sample fluctuations 861 G4double esec = 0.0; << 673 if (lossFluctuationFlag && eloss + lowestKinEnergy <= preStepKinEnergy) { 862 674 863 // Corrections, which cannot be tabulated << 675 eloss = modelManager->SampleFluctuations(currentMaterial, dynParticle, 864 if(isIon) { << 676 tmax, length, eloss, preStepScaledEnergy, 865 currentModel->CorrectionsAlongStep(current << 677 currentMaterialIndex); 866 length, << 867 eloss = std::max(eloss, 0.0); << 868 } 678 } 869 679 870 // Sample fluctuations if not full energy lo << 680 /* 871 if(eloss >= preStepKinEnergy) { << 681 if(-1 < verboseLevel) { 872 eloss = preStepKinEnergy; << 682 G4cout << "eloss(MeV)= " << eloss/MeV 873 << 683 << " fluc= " << (eloss-eloss0)/MeV 874 } else if (lossFluctuationFlag) { << 684 << " currentChargeSquare= " << chargeSquare 875 const G4double tmax = currentModel->MaxSec << 685 << " massRatio= " << massRatio 876 const G4double tcut = std::min(cut, tmax); << 686 << G4endl; 877 G4VEmFluctuationModel* fluc = currentModel << 878 eloss = fluc->SampleFluctuations(currentCo << 879 tcut, tma << 880 /* << 881 if(-1 < verboseLevel) << 882 G4cout << "After fluct: eloss(MeV)= " << << 883 << " fluc= " << (eloss-eloss0)/Me << 884 << " ChargeSqRatio= " << chargeSq << 885 << " massRatio= " << massRatio << << 886 */ << 887 } 687 } >> 688 */ >> 689 >> 690 // Subcutoff and/or deexcitation >> 691 std::vector<G4Track*>* newp = >> 692 SecondariesAlongStep(step, tmax, eloss, preStepScaledEnergy); >> 693 >> 694 if(newp) { >> 695 >> 696 G4int n = newp->size(); >> 697 if(n > 0) { >> 698 fParticleChange.SetNumberOfSecondaries(n); >> 699 G4Track* t; >> 700 G4double e; >> 701 for (G4int i=0; i<n; i++) { >> 702 t = (*newp)[i]; >> 703 e = t->GetKineticEnergy(); >> 704 const G4ParticleDefinition* pd = t->GetDefinition(); >> 705 if (pd != G4Gamma::Gamma() && pd != G4Electron::Electron() ) e += pd->GetPDGMass(); 888 706 889 // deexcitation << 707 preStepKinEnergy -= e; 890 if (useDeexcitation) { << 708 pParticleChange->AddSecondary(t); 891 G4double esecfluo = preStepKinEnergy; << 709 } 892 G4double de = esecfluo; << 710 } 893 atomDeexcitation->AlongStepDeexcitation(sc << 711 delete newp; 894 de << 895 << 896 // sum of de-excitation energies << 897 esecfluo -= de; << 898 << 899 // subtracted from energy loss << 900 if(eloss >= esecfluo) { << 901 esec += esecfluo; << 902 eloss -= esecfluo; << 903 } else { << 904 esec += esecfluo; << 905 eloss = 0.0; << 906 } << 907 } << 908 if(nullptr != subcutProducer && IsRegionForC << 909 subcutProducer->SampleSecondaries(step, sc << 910 } << 911 // secondaries from atomic de-excitation and << 912 if(!scTracks.empty()) { FillSecondariesAlong << 913 << 914 // Energy balance << 915 G4double finalT = preStepKinEnergy - eloss - << 916 if (finalT <= lowestKinEnergy) { << 917 eloss += finalT; << 918 finalT = 0.0; << 919 } else if(isIon) { << 920 fParticleChange.SetProposedCharge( << 921 currentModel->GetParticleCharge(track.Ge << 922 currentM << 923 } 712 } 924 eloss = std::max(eloss, 0.0); << 925 713 926 fParticleChange.SetProposedKineticEnergy(fin << 714 preStepKinEnergy -= eloss; 927 fParticleChange.ProposeLocalEnergyDeposit(el << 715 928 /* << 716 /* 929 if(-1 < verboseLevel) { 717 if(-1 < verboseLevel) { 930 G4double del = finalT + eloss + esec - pre << 718 G4cout << "eloss(MeV)= " << eloss/MeV 931 G4cout << "Final value eloss(MeV)= " << el << 932 << " preStepKinEnergy= " << preStep 719 << " preStepKinEnergy= " << preStepKinEnergy 933 << " postStepKinEnergy= " << finalT << 934 << " de(keV)= " << del/keV << 935 << " lossFlag= " << lossFluctuation 720 << " lossFlag= " << lossFluctuationFlag 936 << " status= " << track.GetTrackSt << 937 << G4endl; 721 << G4endl; 938 } 722 } 939 */ << 723 */ 940 return &fParticleChange; << 941 } << 942 724 943 //....oooOO0OOooo........oooOO0OOooo........oo << 725 if (preStepKinEnergy <= lowestKinEnergy) { 944 726 945 void G4VEnergyLossProcess::FillSecondariesAlon << 727 eloss += preStepKinEnergy; 946 { << 728 preStepKinEnergy = 0.0; 947 const std::size_t n0 = scTracks.size(); << 948 G4double weight = wt; << 949 // weight may be changed by biasing manager << 950 if(biasManager) { << 951 if(biasManager->SecondaryBiasingRegion((G4 << 952 weight *= << 953 biasManager->ApplySecondaryBiasing(scT << 954 } << 955 } << 956 729 957 // fill secondaries << 730 if (hasRestProcess) fParticleChange.SetStatusChange(fStopButAlive); 958 const std::size_t n = scTracks.size(); << 731 else fParticleChange.SetStatusChange(fStopAndKill); 959 fParticleChange.SetNumberOfSecondaries((G4in << 960 << 961 for(std::size_t i=0; i<n; ++i) { << 962 G4Track* t = scTracks[i]; << 963 if(nullptr != t) { << 964 t->SetWeight(weight); << 965 pParticleChange->AddSecondary(t); << 966 G4int pdg = t->GetDefinition()->GetPDGEn << 967 if (i < n0) { << 968 if (pdg == 22) { << 969 t->SetCreatorModelID(gpixeID); << 970 } else if (pdg == 11) { << 971 t->SetCreatorModelID(epixeID); << 972 } else { << 973 t->SetCreatorModelID(biasID); << 974 } << 975 } else { << 976 t->SetCreatorModelID(biasID); << 977 } << 978 } << 979 } 732 } 980 scTracks.clear(); << 733 >> 734 fParticleChange.SetProposedKineticEnergy(preStepKinEnergy); >> 735 fParticleChange.SetLocalEnergyDeposit(eloss); >> 736 >> 737 return &fParticleChange; 981 } 738 } 982 739 983 //....oooOO0OOooo........oooOO0OOooo........oo 740 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 984 741 985 G4VParticleChange* G4VEnergyLossProcess::PostS 742 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track, 986 << 743 const G4Step& step) 987 { 744 { 988 // clear number of interaction lengths in an << 989 theNumberOfInteractionLengthLeft = -1.0; << 990 mfpKinEnergy = DBL_MAX; << 991 << 992 fParticleChange.InitializeForPostStep(track) 745 fParticleChange.InitializeForPostStep(track); 993 const G4double finalT = track.GetKineticEner << 746 G4double finalT = track.GetKineticEnergy(); >> 747 G4double postStepScaledEnergy = finalT*massRatio; 994 748 995 const G4double postStepScaledEnergy = finalT << 749 // Integral approach 996 SelectModel(postStepScaledEnergy); << 750 if (integral) { >> 751 G4bool b; >> 752 G4double postStepLambda = chargeSqRatio* >> 753 (((*theLambdaTable)[currentMaterialIndex])->GetValue(postStepScaledEnergy,b)); 997 754 998 if(!currentModel->IsActive(postStepScaledEne << 755 if(preStepLambda*G4UniformRand() > postStepLambda) 999 return &fParticleChange; << 756 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step); 1000 } 757 } >> 758 >> 759 G4VEmModel* currentModel = SelectModel(postStepScaledEnergy); >> 760 G4double tcut = (*theCuts)[currentMaterialIndex]; >> 761 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); >> 762 G4double tmax = currentModel->MaxSecondaryEnergy(dynParticle); >> 763 1001 /* 764 /* 1002 if(1 < verboseLevel) { << 765 if(0 < verboseLevel) { 1003 G4cout<<GetProcessName()<<" PostStepDoIt: << 766 const G4ParticleDefinition* pd = dynParticle->GetDefinition(); >> 767 G4cout << "G4VEnergyLossProcess::PostStepDoIt: Sample secondary; E= " << finalT/MeV >> 768 << " MeV; model= (" << currentModel->LowEnergyLimit(pd) >> 769 << ", " << currentModel->HighEnergyLimit(pd) << ")" >> 770 << G4endl; 1004 } 771 } 1005 */ 772 */ 1006 // forced process - should happen only once << 1007 if(biasFlag) { << 1008 if(biasManager->ForcedInteractionRegion(( << 1009 biasFlag = false; << 1010 } << 1011 } << 1012 const G4DynamicParticle* dp = track.GetDyna << 1013 << 1014 // Integral approach << 1015 if (fXSType != fEmNoIntegral) { << 1016 const G4double logFinalT = dp->GetLogKine << 1017 G4double lx = GetLambdaForScaledEnergy(po << 1018 lo << 1019 lx = std::max(lx, 0.0); << 1020 << 1021 // if both lg and lx are zero then no int << 1022 if(preStepLambda*G4UniformRand() >= lx) { << 1023 return &fParticleChange; << 1024 } << 1025 } << 1026 773 1027 // define new weight for primary and second << 774 if (tcut < tmax) 1028 G4double weight = fParticleChange.GetParent << 775 SecondariesPostStep(currentModel,currentCouple,dynParticle,tcut,finalT); 1029 if(weightFlag) { << 1030 weight /= biasFactor; << 1031 fParticleChange.ProposeWeight(weight); << 1032 } << 1033 << 1034 const G4double tcut = (*theCuts)[currentCou << 1035 << 1036 // sample secondaries << 1037 secParticles.clear(); << 1038 currentModel->SampleSecondaries(&secParticl << 1039 << 1040 const G4int num0 = (G4int)secParticles.size << 1041 << 1042 // bremsstrahlung splitting or Russian roul << 1043 if(biasManager) { << 1044 if(biasManager->SecondaryBiasingRegion((G << 1045 G4double eloss = 0.0; << 1046 weight *= biasManager->ApplySecondaryBi << 1047 secPart << 1048 track, << 1049 &fParti << 1050 (G4int) << 1051 step.Ge << 1052 if(eloss > 0.0) { << 1053 eloss += fParticleChange.GetLocalEner << 1054 fParticleChange.ProposeLocalEnergyDep << 1055 } << 1056 } << 1057 } << 1058 776 1059 // save secondaries << 777 if (finalT <= 0.0) { 1060 const G4int num = (G4int)secParticles.size( << 778 fParticleChange.SetProposedKineticEnergy(0.0); 1061 if(num > 0) { << 1062 << 1063 fParticleChange.SetNumberOfSecondaries(nu << 1064 G4double time = track.GetGlobalTime(); << 1065 << 1066 G4int n1(0), n2(0); << 1067 if(num0 > mainSecondaries) { << 1068 currentModel->FillNumberOfSecondaries(n << 1069 } << 1070 779 1071 for (G4int i=0; i<num; ++i) { << 780 if (hasRestProcess) fParticleChange.SetStatusChange(fStopButAlive); 1072 if(nullptr != secParticles[i]) { << 781 else fParticleChange.SetStatusChange(fStopAndKill); 1073 G4Track* t = new G4Track(secParticles << 1074 t->SetTouchableHandle(track.GetToucha << 1075 if (biasManager) { << 1076 t->SetWeight(weight * biasManager-> << 1077 } else { << 1078 t->SetWeight(weight); << 1079 } << 1080 if(i < num0) { << 1081 t->SetCreatorModelID(secID); << 1082 } else if(i < num0 + n1) { << 1083 t->SetCreatorModelID(tripletID); << 1084 } else { << 1085 t->SetCreatorModelID(biasID); << 1086 } << 1087 782 1088 //G4cout << "Secondary(post step) has << 783 return &fParticleChange; 1089 // << ", kenergy " << t->GetKin << 1090 // << " time= " << time/ns << " << 1091 pParticleChange->AddSecondary(t); << 1092 } << 1093 } << 1094 } 784 } 1095 785 1096 if(0.0 == fParticleChange.GetProposedKineti << 786 fParticleChange.SetProposedKineticEnergy(finalT); 1097 fAlive == fParticleChange.GetTrackStatus << 1098 if(particle->GetProcessManager()->GetAtRe << 1099 { fParticleChange.ProposeTrackStatus << 1100 else { fParticleChange.ProposeTrackStatus << 1101 } << 1102 787 1103 /* << 788 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step); 1104 if(-1 < verboseLevel) { << 1105 G4cout << "::PostStepDoIt: Sample seconda << 1106 << fParticleChange.GetProposedKineticEner << 1107 << " MeV; model= (" << currentMode << 1108 << ", " << currentModel->HighEner << 1109 << " preStepLambda= " << preStepL << 1110 << " dir= " << track.GetMomentumD << 1111 << " status= " << track.GetTrackS << 1112 << G4endl; << 1113 } << 1114 */ << 1115 return &fParticleChange; << 1116 } 789 } 1117 790 1118 //....oooOO0OOooo........oooOO0OOooo........o 791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1119 792 1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl << 793 void G4VEnergyLossProcess::PrintInfoDefinition() 1121 const G4ParticleDefinition* part, cons << 1122 { 794 { 1123 if (!isMaster || nullptr != baseParticle || << 795 if(-1 < verboseLevel) { 1124 for(std::size_t i=0; i<7; ++i) { << 796 G4cout << G4endl << GetProcessName() << ": tables are built for " 1125 // ionisation table only for ionisation p << 797 << particle->GetParticleName() 1126 if (nullptr == theData->Table(i) || (!isI << 798 << G4endl 1127 continue; << 799 << " dE/dx and range tables from " 1128 } << 800 << G4BestUnit(minKinEnergy,"Energy") 1129 if (-1 < verboseLevel) { << 801 << " to " << G4BestUnit(maxKinEnergy,"Energy") 1130 G4cout << "G4VEnergyLossProcess::StoreP << 802 << " in " << nDEDXBins << " bins." << G4endl 1131 << " " << particle->GetParticleName() << 803 << " Lambda tables from threshold to " 1132 << " " << GetProcessName() << 804 << G4BestUnit(maxKinEnergy,"Energy") 1133 << " " << tnames[i] << " " << theDat << 805 << " in " << nLambdaBins << " bins." >> 806 << G4endl; >> 807 if(theRangeTableForLoss) { >> 808 G4cout << " Step function: finalRange(mm)= " << finalRange/mm >> 809 << ", dRoverRange= " << dRoverRange >> 810 << ", integral: " << integral >> 811 << G4endl; 1134 } 812 } 1135 if (!G4EmTableUtil::StoreTable(this, part << 813 if(thePreciseRangeTable) { 1136 dir, tnames[i], verboseLevel, asci << 814 G4cout << " Precise range table up" 1137 return false; << 815 << " to " << G4BestUnit(maxKinEnergyForRange,"Energy") >> 816 << " in " << nDEDXBinsForRange << " bins." << G4endl; 1138 } 817 } 1139 } << 1140 return true; << 1141 } << 1142 818 1143 //....oooOO0OOooo........oooOO0OOooo........o << 819 if(2 < verboseLevel) { 1144 << 820 G4cout << "DEDXTable address= " << theDEDXTable << G4endl; 1145 G4bool << 821 if(theDEDXTable) G4cout << (*theDEDXTable) << G4endl; 1146 G4VEnergyLossProcess::RetrievePhysicsTable(co << 822 G4cout << "PreciseRangeTable address= " << thePreciseRangeTable << G4endl; 1147 co << 823 if(thePreciseRangeTable) G4cout << (*thePreciseRangeTable) << G4endl; 1148 { << 824 G4cout << "RangeTableForLoss address= " << theRangeTableForLoss << G4endl; 1149 if (!isMaster || nullptr != baseParticle || << 825 if(theRangeTableForLoss) G4cout << (*theRangeTableForLoss) << G4endl; 1150 for(std::size_t i=0; i<7; ++i) { << 826 G4cout << "InverseRangeTable address= " << theInverseRangeTable << G4endl; 1151 // ionisation table only for ionisation p << 827 if(theInverseRangeTable) G4cout << (*theInverseRangeTable) << G4endl; 1152 if (!isIonisation && 1 == i) { continue; << 828 G4cout << "LambdaTable address= " << theLambdaTable << G4endl; 1153 if(!G4EmTableUtil::RetrieveTable(this, pa << 829 if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl; 1154 verboseL << 830 G4cout << "SubLambdaTable address= " << theSubLambdaTable << G4endl; 1155 return false; << 831 if(theSubLambdaTable) G4cout << (*theSubLambdaTable) << G4endl; 1156 } 832 } 1157 } 833 } 1158 return true; << 1159 } 834 } 1160 835 1161 //....oooOO0OOooo........oooOO0OOooo........o 836 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1162 837 1163 G4double G4VEnergyLossProcess::GetDEDXDispers << 838 void G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p) 1164 const G4Mat << 1165 const G4Dyn << 1166 G4dou << 1167 { 839 { 1168 DefineMaterial(couple); << 840 if(theDEDXTable) theDEDXTable->clearAndDestroy(); 1169 G4double ekin = dp->GetKineticEnergy(); << 841 theDEDXTable = p; 1170 SelectModel(ekin*massRatio); << 1171 G4double tmax = currentModel->MaxSecondaryK << 1172 G4double tcut = std::min(tmax,(*theCuts)[cu << 1173 G4double d = 0.0; << 1174 G4VEmFluctuationModel* fm = currentModel->G << 1175 if(nullptr != fm) { d = fm->Dispersion(curr << 1176 return d; << 1177 } 842 } 1178 843 1179 //....oooOO0OOooo........oooOO0OOooo........o 844 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1180 845 1181 G4double << 846 void G4VEnergyLossProcess::SetRangeTable(G4PhysicsTable* p) 1182 G4VEnergyLossProcess::CrossSectionPerVolume(G << 1183 c << 1184 G << 1185 { 847 { 1186 // Cross section per volume is calculated << 848 if(thePreciseRangeTable) thePreciseRangeTable->clearAndDestroy(); 1187 DefineMaterial(couple); << 849 if(theDEDXAtMaxEnergy) delete [] theDEDXAtMaxEnergy; 1188 G4double cross = 0.0; << 850 if(theRangeAtMaxEnergy) delete [] theRangeAtMaxEnergy; 1189 if (nullptr != theLambdaTable) { << 851 1190 cross = GetLambdaForScaledEnergy(kineticE << 852 thePreciseRangeTable = p; 1191 logKinet << 853 if(p) { 1192 } else { << 854 size_t n = p->length(); 1193 SelectModel(kineticEnergy*massRatio); << 855 G4PhysicsVector* pv = (*p)[0]; 1194 cross = (!baseMat) ? biasFactor : biasFac << 856 G4double emax = pv->GetLowEdgeEnergy(pv->GetVectorLength()); 1195 cross *= (currentModel->CrossSectionPerVo << 857 G4bool b; 1196 << 858 theDEDXAtMaxEnergy = new G4double [n]; >> 859 theRangeAtMaxEnergy = new G4double [n]; >> 860 >> 861 for (size_t i=0; i<n; i++) { >> 862 pv = (*p)[i]; >> 863 G4double r2 = pv->GetValue(emax, b); >> 864 G4double dedx = ((*theDEDXTable)[i])->GetValue(emax,b); >> 865 theDEDXAtMaxEnergy[i] = dedx; >> 866 theRangeAtMaxEnergy[i] = r2; >> 867 //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= " << r2 >> 868 // << " dedx= " << dedx << G4endl; >> 869 } 1197 } 870 } 1198 return std::max(cross, 0.0); << 1199 } 871 } 1200 872 1201 //....oooOO0OOooo........oooOO0OOooo........o 873 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1202 874 1203 G4double G4VEnergyLossProcess::MeanFreePath(c << 875 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p) 1204 { 876 { 1205 DefineMaterial(track.GetMaterialCutsCouple( << 877 if(theRangeTableForLoss) theRangeTableForLoss->clearAndDestroy(); 1206 const G4double kinEnergy = track.GetKine << 878 theRangeTableForLoss = p; 1207 const G4double logKinEnergy = track.GetDyna << 1208 const G4double cs = GetLambdaForScaledEnerg << 1209 << 1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; << 1211 } 879 } 1212 880 1213 //....oooOO0OOooo........oooOO0OOooo........o 881 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1214 882 1215 G4double G4VEnergyLossProcess::ContinuousStep << 883 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p) 1216 << 1217 << 1218 { 884 { 1219 return AlongStepGetPhysicalInteractionLengt << 885 theSecondaryRangeTable = p; 1220 } 886 } 1221 887 1222 //....oooOO0OOooo........oooOO0OOooo........o 888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1223 889 1224 G4double G4VEnergyLossProcess::GetMeanFreePat << 890 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p) 1225 const G4Track& t << 1226 G4double, << 1227 G4ForceCondition << 1228 << 1229 { 891 { 1230 *condition = NotForced; << 892 if(theInverseRangeTable) theInverseRangeTable->clearAndDestroy(); 1231 return MeanFreePath(track); << 893 theInverseRangeTable = p; 1232 } 894 } 1233 895 1234 //....oooOO0OOooo........oooOO0OOooo........o 896 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1235 897 1236 G4double G4VEnergyLossProcess::GetContinuousS << 898 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p) 1237 const G4Track&, << 1238 G4double, G4double, G4double& << 1239 { 899 { 1240 return DBL_MAX; << 900 if(theLambdaTable) theLambdaTable->clearAndDestroy(); >> 901 theLambdaTable = p; >> 902 tablesAreBuilt = true; 1241 } 903 } 1242 904 1243 //....oooOO0OOooo........oooOO0OOooo........o 905 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1244 906 1245 G4PhysicsVector* << 907 void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p) 1246 G4VEnergyLossProcess::LambdaPhysicsVector(con << 1247 G4d << 1248 { 908 { 1249 DefineMaterial(couple); << 909 if(theSubLambdaTable) theSubLambdaTable->clearAndDestroy(); 1250 G4PhysicsVector* v = (*theLambdaTable)[base << 910 theSubLambdaTable = p; 1251 return new G4PhysicsVector(*v); << 911 if (nSCoffRegions) { >> 912 for (G4int i=0; i<nSCoffRegions; i++) { >> 913 scoffProcessors[i]->SetLambdaSubTable(theSubLambdaTable); >> 914 } >> 915 } 1252 } 916 } 1253 917 1254 //....oooOO0OOooo........oooOO0OOooo........o 918 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1255 919 1256 void << 920 G4PhysicsVector* G4VEnergyLossProcess::DEDXPhysicsVector(const G4MaterialCutsCouple*) 1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT << 1258 { 921 { 1259 if(1 < verboseLevel) { << 922 G4int nbins = nDEDXBins; 1260 G4cout << "### Set DEDX table " << p << " << 923 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins); 1261 << " " << theDEDXunRestrictedTable << << 924 return v; 1262 << " for " << particle->GetParticl << 1263 << " and process " << GetProcessNa << 1264 << " type=" << tType << " isIonisation:" << 1265 } << 1266 if(fTotal == tType) { << 1267 theDEDXunRestrictedTable = p; << 1268 } else if(fRestricted == tType) { << 1269 theDEDXTable = p; << 1270 if(isMaster && nullptr == baseParticle) { << 1271 theData->UpdateTable(theDEDXTable, 0); << 1272 } << 1273 } else if(fIsIonisation == tType) { << 1274 theIonisationTable = p; << 1275 if(isMaster && nullptr == baseParticle) { << 1276 theData->UpdateTable(theIonisationTable << 1277 } << 1278 } << 1279 } 925 } 1280 926 1281 //....oooOO0OOooo........oooOO0OOooo........o 927 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1282 928 1283 void G4VEnergyLossProcess::SetCSDARangeTable( << 929 G4PhysicsVector* G4VEnergyLossProcess::DEDXPhysicsVectorForPreciseRange( >> 930 const G4MaterialCutsCouple*) 1284 { 931 { 1285 theCSDARangeTable = p; << 932 G4int nbins = nDEDXBinsForRange; >> 933 G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergyForRange, nbins); >> 934 return v; 1286 } 935 } 1287 936 1288 //....oooOO0OOooo........oooOO0OOooo........o 937 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1289 938 1290 void G4VEnergyLossProcess::SetRangeTableForLo << 939 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(const G4MaterialCutsCouple* couple) 1291 { 940 { 1292 theRangeTableForLoss = p; << 941 G4double cut = (*theCuts)[couple->GetIndex()]; >> 942 G4int nbins = nLambdaBins; >> 943 G4double tmin = std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut), >> 944 minKinEnergy); >> 945 if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy; >> 946 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nbins); >> 947 return v; 1293 } 948 } 1294 949 1295 //....oooOO0OOooo........oooOO0OOooo........o 950 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1296 951 1297 void G4VEnergyLossProcess::SetInverseRangeTab << 952 G4PhysicsVector* G4VEnergyLossProcess::SubLambdaPhysicsVector(const G4MaterialCutsCouple* couple) 1298 { 953 { 1299 theInverseRangeTable = p; << 954 return LambdaPhysicsVector(couple); 1300 } 955 } 1301 956 1302 //....oooOO0OOooo........oooOO0OOooo........o 957 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1303 958 1304 void G4VEnergyLossProcess::SetLambdaTable(G4P << 959 G4double G4VEnergyLossProcess::MicroscopicCrossSection(G4double kineticEnergy, >> 960 const G4MaterialCutsCouple* couple) 1305 { 961 { 1306 if(1 < verboseLevel) { << 962 // Cross section per atom is calculated 1307 G4cout << "### Set Lambda table " << p << << 963 DefineMaterial(couple); 1308 << " for " << particle->GetParticl << 964 G4double cross = 0.0; 1309 << " and process " << GetProcessNa << 965 G4bool b; 1310 } << 966 if(theLambdaTable) { 1311 theLambdaTable = p; << 967 cross = (((*theLambdaTable)[currentMaterialIndex])-> 1312 tablesAreBuilt = true; << 968 GetValue(kineticEnergy, b)); 1313 969 1314 if(isMaster && nullptr != p) { << 970 cross /= currentMaterial->GetTotNbOfAtomsPerVolume(); 1315 delete theEnergyOfCrossSectionMax; << 1316 theEnergyOfCrossSectionMax = nullptr; << 1317 if(fEmTwoPeaks == fXSType) { << 1318 if(nullptr != fXSpeaks) { << 1319 for(auto & ptr : *fXSpeaks) { delete ptr; } << 1320 delete fXSpeaks; << 1321 } << 1322 G4LossTableBuilder* bld = lManager->Get << 1323 fXSpeaks = G4EmUtility::FillPeaksStruct << 1324 if(nullptr == fXSpeaks) { fXSType = fEm << 1325 } << 1326 if(fXSType == fEmOnePeak) { << 1327 theEnergyOfCrossSectionMax = G4EmUtilit << 1328 if(nullptr == theEnergyOfCrossSectionMa << 1329 } << 1330 } 971 } >> 972 >> 973 return cross; 1331 } 974 } 1332 975 1333 //....oooOO0OOooo........oooOO0OOooo........o 976 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1334 977 1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe << 978 G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track, >> 979 G4double s, >> 980 G4ForceCondition* cond) 1336 { 981 { 1337 theEnergyOfCrossSectionMax = p; << 982 return GetMeanFreePath(track, s, cond); 1338 } 983 } 1339 984 1340 //....oooOO0OOooo........oooOO0OOooo........o 985 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1341 986 1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: << 987 G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, >> 988 G4double x, G4double y, G4double& z) 1343 { 989 { 1344 fXSpeaks = ptr; << 990 return GetContinuousStepLimit(track, x, y, z); 1345 } 991 } 1346 992 1347 //....oooOO0OOooo........oooOO0OOooo........o 993 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1348 994 1349 const G4Element* G4VEnergyLossProcess::GetCur << 995 void G4VEnergyLossProcess::SetStepLimits(G4double v1, G4double v2) 1350 { 996 { 1351 return (nullptr != currentModel) << 997 dRoverRange = v1; 1352 ? currentModel->GetCurrentElement(current << 998 finalRange = v2; >> 999 if (dRoverRange > 1.0) dRoverRange = 1.0; 1353 } 1000 } 1354 1001 1355 //....oooOO0OOooo........oooOO0OOooo........o 1002 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1356 1003 1357 void G4VEnergyLossProcess::SetCrossSectionBia << 1004 void G4VEnergyLossProcess::SetIntegral(G4bool val) 1358 << 1359 { 1005 { 1360 if(f > 0.0) { << 1006 if(integral != val) { 1361 biasFactor = f; << 1007 if(val) dRoverRange = defaultIntegralRange; 1362 weightFlag = flag; << 1008 else dRoverRange = defaultRoverRange; 1363 if(1 < verboseLevel) { << 1364 G4cout << "### SetCrossSectionBiasingFa << 1365 << " process " << GetProcessName << 1366 << " biasFactor= " << f << " wei << 1367 << G4endl; << 1368 } << 1369 } 1009 } >> 1010 integral = val; 1370 } 1011 } 1371 1012 1372 //....oooOO0OOooo........oooOO0OOooo........o 1013 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1373 1014 1374 void G4VEnergyLossProcess::ActivateForcedInte << 1015 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2) 1375 << 1016 { 1376 << 1017 dRoverRange = v1; 1377 { << 1018 finalRange = v2; 1378 if(nullptr == biasManager) { biasManager = << 1019 if (dRoverRange > 0.999) dRoverRange = 1.0; 1379 if(1 < verboseLevel) { << 1380 G4cout << "### ActivateForcedInteraction: << 1381 << " process " << GetProcessName() << 1382 << " length(mm)= " << length/mm << 1383 << " in G4Region <" << region << 1384 << "> weightFlag= " << flag << 1385 << G4endl; << 1386 } << 1387 weightFlag = flag; << 1388 biasManager->ActivateForcedInteraction(leng << 1389 } 1020 } 1390 1021 1391 //....oooOO0OOooo........oooOO0OOooo........o 1022 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1392 1023 1393 void << 1024 void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p) 1394 G4VEnergyLossProcess::ActivateSecondaryBiasin << 1025 { 1395 << 1026 particle = p; 1396 << 1027 baseParticle = DefineBaseParticle(particle); 1397 { << 1398 if (0.0 <= factor) { << 1399 // Range cut can be applied only for e- << 1400 if(0.0 == factor && secondaryParticle != << 1401 { return; } << 1402 << 1403 if(nullptr == biasManager) { biasManager << 1404 biasManager->ActivateSecondaryBiasing(reg << 1405 if(1 < verboseLevel) { << 1406 G4cout << "### ActivateSecondaryBiasing << 1407 << " process " << GetProcessName << 1408 << " factor= " << factor << 1409 << " in G4Region <" << region << 1410 << "> energyLimit(MeV)= " << ene << 1411 << G4endl; << 1412 } << 1413 } << 1414 } 1028 } 1415 1029 1416 //....oooOO0OOooo........oooOO0OOooo........o 1030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1417 1031 1418 void G4VEnergyLossProcess::SetIonisation(G4bo << 1032 void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p) 1419 { 1033 { 1420 isIonisation = val; << 1034 baseParticle = p; 1421 aGPILSelection = (val) ? CandidateForSelect << 1422 } 1035 } 1423 1036 1424 //....oooOO0OOooo........oooOO0OOooo........o 1037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1425 1038 1426 void G4VEnergyLossProcess::SetLinearLossLimi << 1039 void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p) 1427 { 1040 { 1428 if(0.0 < val && val < 1.0) { << 1041 secondaryParticle = p; 1429 linLossLimit = val; << 1430 actLinLossLimit = true; << 1431 } else { PrintWarning("SetLinearLossLimit", << 1432 } 1042 } 1433 1043 1434 //....oooOO0OOooo........oooOO0OOooo........o 1044 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1435 1045 1436 void G4VEnergyLossProcess::SetStepFunction(G4 << 1046 G4bool G4VEnergyLossProcess::StorePhysicsTable(G4ParticleDefinition* part, >> 1047 const G4String& directory, >> 1048 G4bool ascii) 1437 { 1049 { 1438 if(0.0 < v1 && 0.0 < v2) { << 1050 G4bool res = true; 1439 dRoverRange = std::min(1.0, v1); << 1051 if ( baseParticle ) return res; 1440 finalRange = std::min(v2, 1.e+50); << 1052 G4bool yes = true; 1441 } else { << 1053 1442 PrintWarning("SetStepFunctionV1", v1); << 1054 if ( theDEDXTable ) { 1443 PrintWarning("SetStepFunctionV2", v2); << 1055 const G4String name = GetPhysicsTableFileName(part,directory,"DEDX",ascii); >> 1056 yes = theDEDXTable->StorePhysicsTable(name,ascii); >> 1057 if( !yes ) res = false; 1444 } 1058 } 1445 } << 1446 1059 1447 //....oooOO0OOooo........oooOO0OOooo........o << 1060 if ( thePreciseRangeTable ) { >> 1061 const G4String name = GetPhysicsTableFileName(part,directory,"PreciseRange",ascii); >> 1062 yes = thePreciseRangeTable->StorePhysicsTable(name,ascii); >> 1063 if( !yes ) res = false; >> 1064 } 1448 1065 1449 void G4VEnergyLossProcess::SetLowestEnergyLim << 1066 if ( theRangeTableForLoss ) { 1450 { << 1067 const G4String name = GetPhysicsTableFileName(part,directory,"Range",ascii); 1451 if(1.e-18 < val && val < 1.e+50) { lowestKi << 1068 yes = theRangeTableForLoss->StorePhysicsTable(name,ascii); 1452 else { PrintWarning("SetLowestEnergyLimit", << 1069 if( !yes ) res = false; 1453 } << 1070 } 1454 1071 1455 //....oooOO0OOooo........oooOO0OOooo........o << 1072 if ( theInverseRangeTable ) { >> 1073 const G4String name = GetPhysicsTableFileName(part,directory,"InverseRange",ascii); >> 1074 yes = theInverseRangeTable->StorePhysicsTable(name,ascii); >> 1075 if( !yes ) res = false; >> 1076 } 1456 1077 1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i << 1078 if ( theLambdaTable ) { 1458 { << 1079 const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii); 1459 if(2 < n && n < 1000000000) { << 1080 yes = theLambdaTable->StorePhysicsTable(name,ascii); 1460 nBins = n; << 1081 if( !yes ) res = false; 1461 actBinning = true; << 1082 } >> 1083 >> 1084 if ( theSubLambdaTable ) { >> 1085 const G4String name = GetPhysicsTableFileName(part,directory,"SubLambda",ascii); >> 1086 yes = theSubLambdaTable->StorePhysicsTable(name,ascii); >> 1087 if( !yes ) res = false; >> 1088 } >> 1089 if ( res ) { >> 1090 G4cout << "Physics tables are stored for " << particle->GetParticleName() >> 1091 << " and process " << GetProcessName() >> 1092 << " in the directory <" << directory >> 1093 << "> " << G4endl; 1462 } else { 1094 } else { 1463 G4double e = (G4double)n; << 1095 G4cout << "Fail to store Physics Tables for " << particle->GetParticleName() 1464 PrintWarning("SetDEDXBinning", e); << 1096 << " and process " << GetProcessName() 1465 } << 1097 << " in the directory <" << directory >> 1098 << "> " << G4endl; >> 1099 } >> 1100 return res; 1466 } 1101 } 1467 1102 1468 //....oooOO0OOooo........oooOO0OOooo........o << 1103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 1469 1104 1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 << 1105 G4bool G4VEnergyLossProcess::RetrievePhysicsTable(G4ParticleDefinition* part, >> 1106 const G4String& directory, >> 1107 G4bool ascii) 1471 { 1108 { 1472 if(1.e-18 < e && e < maxKinEnergy) { << 1109 G4bool res = true; 1473 minKinEnergy = e; << 1110 currentCouple = 0; 1474 actMinKinEnergy = true; << 1111 preStepLambda = 0.0; 1475 } else { PrintWarning("SetMinKinEnergy", e) << 1112 if(0 < verboseLevel) { 1476 } << 1113 G4cout << "========================================================" << G4endl; >> 1114 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for " >> 1115 << part->GetParticleName() << " and process " << GetProcessName() >> 1116 << "; tables_are_built= " << tablesAreBuilt >> 1117 << G4endl; >> 1118 } 1477 1119 1478 //....oooOO0OOooo........oooOO0OOooo........o << 1120 const G4String particleName = part->GetParticleName(); >> 1121 if( !particle ) { >> 1122 particle = part; >> 1123 baseParticle = DefineBaseParticle(particle); >> 1124 } >> 1125 >> 1126 if(particleName != "GenericIon" && >> 1127 part->GetParticleType() == "nucleus" && >> 1128 part->GetParticleSubType() == "generic") >> 1129 { >> 1130 (G4LossTableManager::Instance())->RegisterIon(part, this); >> 1131 return res; >> 1132 } >> 1133 >> 1134 if(tablesAreBuilt) return res; >> 1135 Initialise(); >> 1136 >> 1137 // Recalculation is needed because cuts were changed or recalculation is forced >> 1138 G4LossTableManager* lManager = G4LossTableManager::Instance(); >> 1139 if ( lManager->IsRecalcNeeded(particle)) { >> 1140 >> 1141 G4bool yes = true; >> 1142 G4bool fpi = true; >> 1143 if ( !baseParticle ) { >> 1144 G4PhysicsTable* table; >> 1145 G4String filename; >> 1146 const G4ProductionCutsTable* theCoupleTable= >> 1147 G4ProductionCutsTable::GetProductionCutsTable(); >> 1148 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 1149 >> 1150 filename = GetPhysicsTableFileName(part,directory,"DEDX",ascii); >> 1151 >> 1152 table = new G4PhysicsTable(numOfCouples); >> 1153 yes = table->ExistPhysicsTable(filename); >> 1154 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); >> 1155 if(yes) { >> 1156 SetDEDXTable(table); >> 1157 if (-1 < verboseLevel) { >> 1158 G4cout << "DEDX table for " << particleName << " is retrieved from <" >> 1159 << filename << ">" >> 1160 << G4endl; >> 1161 } >> 1162 } else { >> 1163 fpi = false; >> 1164 table->clearAndDestroy(); >> 1165 if (0 < verboseLevel) { >> 1166 G4cout << "DEDX table for " << particleName << " from file <" >> 1167 << filename << "> is not retrieved" >> 1168 << G4endl; >> 1169 } >> 1170 } 1479 1171 1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 << 1172 filename = GetPhysicsTableFileName(part,directory,"Range",ascii); 1481 { << 1173 table = new G4PhysicsTable(numOfCouples); 1482 if(minKinEnergy < e && e < 1.e+50) { << 1174 yes = table->ExistPhysicsTable(filename); 1483 maxKinEnergy = e; << 1175 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); 1484 actMaxKinEnergy = true; << 1176 if(yes) { 1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS << 1177 SetRangeTableForLoss(table); 1486 } else { PrintWarning("SetMaxKinEnergy", e) << 1178 if (-1 < verboseLevel) { 1487 } << 1179 G4cout << "Range table for loss for " << particleName << " is retrieved from <" >> 1180 << filename << ">" >> 1181 << G4endl; >> 1182 } >> 1183 } else { >> 1184 table->clearAndDestroy(); >> 1185 if(fpi) { >> 1186 res = false; >> 1187 G4cout << "Range table for loss for " << particleName << " from file <" >> 1188 << filename << "> is not retrieved" >> 1189 << G4endl; >> 1190 } >> 1191 } 1488 1192 1489 //....oooOO0OOooo........oooOO0OOooo........o << 1193 filename = GetPhysicsTableFileName(part,directory,"PreciseRange",ascii); >> 1194 table = new G4PhysicsTable(numOfCouples); >> 1195 yes = table->ExistPhysicsTable(filename); >> 1196 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); >> 1197 if(yes) { >> 1198 SetRangeTable(table); >> 1199 if (-1 < verboseLevel) { >> 1200 G4cout << "Precise Range table for " << particleName << " is retrieved from <" >> 1201 << filename << ">" >> 1202 << G4endl; >> 1203 } >> 1204 } else { >> 1205 table->clearAndDestroy(); >> 1206 if(fpi) { >> 1207 res = false; >> 1208 G4cout << "Precise Range table for loss for " << particleName << " does not exist" >> 1209 << G4endl; >> 1210 } >> 1211 } 1490 1212 1491 void G4VEnergyLossProcess::PrintWarning(const << 1213 filename = GetPhysicsTableFileName(part,directory,"InverseRange",ascii); 1492 { << 1214 table = new G4PhysicsTable(numOfCouples); 1493 G4String ss = "G4VEnergyLossProcess::" + ti << 1215 yes = table->ExistPhysicsTable(filename); 1494 G4ExceptionDescription ed; << 1216 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); 1495 ed << "Parameter is out of range: " << val << 1217 if(yes) { 1496 << " it will have no effect!\n" << " Pr << 1218 SetInverseRangeTable(table); 1497 << GetProcessName() << " nbins= " << nB << 1219 if (-1 < verboseLevel) { 1498 << " Emin(keV)= " << minKinEnergy/keV << 1220 G4cout << "InverseRange table for " << particleName << " is retrieved from <" 1499 << " Emax(GeV)= " << maxKinEnergy/GeV; << 1221 << filename << ">" 1500 G4Exception(ss, "em0044", JustWarning, ed); << 1222 << G4endl; 1501 } << 1223 } >> 1224 } else { >> 1225 table->clearAndDestroy(); >> 1226 if(fpi) { >> 1227 res = false; >> 1228 G4cout << "InverseRange table for " << particleName << " from file <" >> 1229 << filename << "> is not retrieved" >> 1230 << G4endl; 1502 1231 1503 //....oooOO0OOooo........oooOO0OOooo........o << 1232 } >> 1233 } 1504 1234 1505 void G4VEnergyLossProcess::ProcessDescription << 1235 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii); 1506 { << 1236 table = new G4PhysicsTable(numOfCouples); 1507 if(nullptr != particle) { StreamInfo(out, * << 1237 yes = table->ExistPhysicsTable(filename); >> 1238 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); >> 1239 if(yes) { >> 1240 SetLambdaTable(table); >> 1241 if (-1 < verboseLevel) { >> 1242 G4cout << "Lambda table for " << particleName << " is retrieved from <" >> 1243 << filename << ">" >> 1244 << G4endl; >> 1245 } >> 1246 } else { >> 1247 table->clearAndDestroy(); >> 1248 if(fpi) { >> 1249 res = false; >> 1250 G4cout << "Lambda table for " << particleName << " from file <" >> 1251 << filename << "> is not retrieved" >> 1252 << G4endl; >> 1253 } >> 1254 } >> 1255 >> 1256 filename = GetPhysicsTableFileName(part,directory,"SubLambda",ascii); >> 1257 table = new G4PhysicsTable(numOfCouples); >> 1258 yes = table->ExistPhysicsTable(filename); >> 1259 if(yes) yes = table->RetrievePhysicsTable(filename,ascii); >> 1260 if(yes) { >> 1261 SetSubLambdaTable(table); >> 1262 if (-1 < verboseLevel) { >> 1263 G4cout << "SubLambda table for " << particleName << " is retrieved from <" >> 1264 << filename << ">" >> 1265 << G4endl; >> 1266 } >> 1267 } else { >> 1268 table->clearAndDestroy(); >> 1269 if(nSCoffRegions) { >> 1270 res=false; >> 1271 G4cout << "SubLambda table for " << particleName << " from file <" >> 1272 << filename << "> is not retrieved" >> 1273 << G4endl; >> 1274 } >> 1275 } >> 1276 if(res) PrintInfoDefinition(); >> 1277 else { >> 1278 G4cout << "### BuildPhysicsTable will be requested for " << GetProcessName() >> 1279 << " for " << particleName << G4endl; >> 1280 } >> 1281 } >> 1282 tablesAreBuilt = true; >> 1283 } >> 1284 >> 1285 lManager->RetrievePhysicsTables(particle, this); >> 1286 >> 1287 return res; 1508 } 1288 } 1509 1289 >> 1290 1510 //....oooOO0OOooo........oooOO0OOooo........o 1291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1292 1511 1293