Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4VEnergyLossProcess.cc,v 1.102 2007/05/08 17:52:00 vnivanch Exp $ >> 27 // GEANT4 tag $Name: geant4-08-03-patch-02 $ >> 28 // 26 // ------------------------------------------- 29 // ------------------------------------------------------------------- 27 // 30 // 28 // GEANT4 Class file 31 // GEANT4 Class file 29 // 32 // 30 // 33 // 31 // File name: G4VEnergyLossProcess 34 // File name: G4VEnergyLossProcess 32 // 35 // 33 // Author: Vladimir Ivanchenko 36 // Author: Vladimir Ivanchenko 34 // 37 // 35 // Creation date: 03.01.2002 38 // Creation date: 03.01.2002 36 // 39 // 37 // Modifications: Vladimir Ivanchenko << 40 // Modifications: 38 // 41 // >> 42 // 13-11-02 Minor fix - use normalised direction (V.Ivanchenko) >> 43 // 04-12-02 Minor change in PostStepDoIt (V.Ivanchenko) >> 44 // 23-12-02 Change interface in order to move to cut per region (V.Ivanchenko) >> 45 // 26-12-02 Secondary production moved to derived classes (V.Ivanchenko) >> 46 // 04-01-03 Fix problem of very small steps for ions (V.Ivanchenko) >> 47 // 20-01-03 Migrade to cut per region (V.Ivanchenko) >> 48 // 24-01-03 Temporarily close a control on usage of couples (V.Ivanchenko) >> 49 // 24-01-03 Make models region aware (V.Ivanchenko) >> 50 // 05-02-03 Fix compilation warnings (V.Ivanchenko) >> 51 // 06-02-03 Add control on tmax in PostStepDoIt (V.Ivanchenko) >> 52 // 13-02-03 SubCutoffProcessors defined for regions (V.Ivanchenko) >> 53 // 15-02-03 Lambda table can be scaled (V.Ivanchenko) >> 54 // 17-02-03 Fix problem of store/restore tables (V.Ivanchenko) >> 55 // 18-02-03 Add control on CutCouple usage (V.Ivanchenko) >> 56 // 26-02-03 Simplify control on GenericIons (V.Ivanchenko) >> 57 // 06-03-03 Control on GenericIons using SubType + update verbose (V.Ivanchenko) >> 58 // 10-03-03 Add Ion registration (V.Ivanchenko) >> 59 // 22-03-03 Add Initialisation of cash (V.Ivanchenko) >> 60 // 26-03-03 Remove finalRange modification (V.Ivanchenko) >> 61 // 09-04-03 Fix problem of negative range limit for non integral (V.Ivanchenko) >> 62 // 26-04-03 Fix retrieve tables (V.Ivanchenko) >> 63 // 06-05-03 Set defalt finalRange = 1 mm (V.Ivanchenko) >> 64 // 12-05-03 Update range calculations + lowKinEnergy (V.Ivanchenko) >> 65 // 13-05-03 Add calculation of precise range (V.Ivanchenko) >> 66 // 23-05-03 Remove tracking cuts (V.Ivanchenko) >> 67 // 03-06-03 Fix initialisation problem for STD ionisation (V.Ivanchenko) >> 68 // 21-07-03 Add UpdateEmModel method (V.Ivanchenko) >> 69 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko) >> 70 // 04-11-03 Add checks in RetrievePhysicsTable (V.Ivanchenko) >> 71 // 12-11-03 G4EnergyLossSTD -> G4EnergyLossProcess (V.Ivanchenko) >> 72 // 21-01-04 Migrade to G4ParticleChangeForLoss (V.Ivanchenko) >> 73 // 27-02-04 Fix problem of loss in low presure gases, cleanup precise range >> 74 // calculation, use functions ForLoss in AlongStepDoIt (V.Ivanchenko) >> 75 // 10-03-04 Fix a problem of Precise Range table (V.Ivanchenko) >> 76 // 19-03-04 Fix a problem energy below lowestKinEnergy (V.Ivanchenko) >> 77 // 31-03-04 Fix a problem of retrieve tables (V.Ivanchenko) >> 78 // 21-07-04 Check weather AtRest are active or not (V.Ivanchenko) >> 79 // 03-08-04 Add pointer of DEDX table to all processes (V.Ivanchenko) >> 80 // 06-08-04 Clear up names of member functions (V.Ivanchenko) >> 81 // 06-08-04 Clear up names of member functions (V.Ivanchenko) >> 82 // 27-08-04 Add NeedBuildTables method (V.Ivanchneko) >> 83 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko) >> 84 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko) >> 85 // 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko) >> 86 // 11-04-05 Use MaxSecondaryEnergy from a model (V.Ivanchenko) >> 87 // 25-07-05 Add extra protection PostStep for non-integral mode (V.Ivanchenko) >> 88 // 12-08-05 Integral=false; SetStepFunction(0.2, 0.1*mm) (mma) >> 89 // 18-08-05 Return back both AlongStep and PostStep from 7.0 (V.Ivanchenko) >> 90 // 02-09-05 Default StepFunction 0.2 1 mm + integral (V.Ivanchenko) >> 91 // 04-09-05 default lambdaFactor 0.8 (V.Ivanchenko) >> 92 // 05-10-05 protection against 0 energy loss added (L.Urban) >> 93 // 17-10-05 protection above has been removed (L.Urban) >> 94 // 06-01-06 reset currentCouple when StepFunction is changed (V.Ivanchenko) >> 95 // 10-01-06 PreciseRange -> CSDARange (V.Ivantchenko) >> 96 // 18-01-06 Clean up subcutoff including recalculation of presafety (VI) >> 97 // 20-01-06 Introduce G4EmTableType and reducing number of methods (VI) >> 98 // 22-03-06 Add control on warning printout AlongStep (VI) >> 99 // 23-03-06 Use isIonisation flag (V.Ivanchenko) >> 100 // 07-06-06 Do not reflect AlongStep in subcutoff regime (V.Ivanchenko) >> 101 // 14-01-07 add SetEmModel(index) and SetFluctModel() (mma) >> 102 // 16-01-07 add IonisationTable and IonisationSubTable (V.Ivanchenko) >> 103 // 16-02-07 set linLossLimit=1.e-6 (V.Ivanchenko) >> 104 // 13-03-07 use SafetyHelper instead of navigator (V.Ivanchenko) 39 // 105 // 40 // Class Description: 106 // Class Description: 41 // 107 // 42 // It is the unified energy loss process it ca 108 // It is the unified energy loss process it calculates the continuous 43 // energy loss for charged particles using a s 109 // energy loss for charged particles using a set of Energy Loss 44 // models valid for different energy regions. 110 // models valid for different energy regions. There are a possibility 45 // to create and access to dE/dx and range tab 111 // to create and access to dE/dx and range tables, or to calculate 46 // that information on fly. 112 // that information on fly. 47 // ------------------------------------------- 113 // ------------------------------------------------------------------- 48 // 114 // 49 //....oooOO0OOooo........oooOO0OOooo........oo 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 51 117 52 #include "G4VEnergyLossProcess.hh" 118 #include "G4VEnergyLossProcess.hh" 53 #include "G4PhysicalConstants.hh" << 54 #include "G4SystemOfUnits.hh" << 55 #include "G4ProcessManager.hh" << 56 #include "G4LossTableManager.hh" 119 #include "G4LossTableManager.hh" 57 #include "G4LossTableBuilder.hh" << 58 #include "G4Step.hh" 120 #include "G4Step.hh" 59 #include "G4ParticleDefinition.hh" 121 #include "G4ParticleDefinition.hh" 60 #include "G4ParticleTable.hh" << 61 #include "G4EmParameters.hh" << 62 #include "G4EmUtility.hh" << 63 #include "G4EmTableUtil.hh" << 64 #include "G4VEmModel.hh" 122 #include "G4VEmModel.hh" 65 #include "G4VEmFluctuationModel.hh" 123 #include "G4VEmFluctuationModel.hh" 66 #include "G4DataVector.hh" 124 #include "G4DataVector.hh" 67 #include "G4PhysicsLogVector.hh" 125 #include "G4PhysicsLogVector.hh" 68 #include "G4VParticleChange.hh" 126 #include "G4VParticleChange.hh" >> 127 #include "G4Gamma.hh" 69 #include "G4Electron.hh" 128 #include "G4Electron.hh" >> 129 #include "G4Positron.hh" >> 130 #include "G4Proton.hh" 70 #include "G4ProcessManager.hh" 131 #include "G4ProcessManager.hh" 71 #include "G4UnitsTable.hh" 132 #include "G4UnitsTable.hh" >> 133 #include "G4GenericIon.hh" >> 134 #include "G4ProductionCutsTable.hh" 72 #include "G4Region.hh" 135 #include "G4Region.hh" 73 #include "G4RegionStore.hh" 136 #include "G4RegionStore.hh" 74 #include "G4PhysicsTableHelper.hh" 137 #include "G4PhysicsTableHelper.hh" 75 #include "G4SafetyHelper.hh" 138 #include "G4SafetyHelper.hh" 76 #include "G4EmDataHandler.hh" << 77 #include "G4TransportationManager.hh" 139 #include "G4TransportationManager.hh" 78 #include "G4VAtomDeexcitation.hh" << 79 #include "G4VSubCutProducer.hh" << 80 #include "G4EmBiasingManager.hh" << 81 #include "G4Log.hh" << 82 #include <iostream> << 83 140 84 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 142 86 namespace << 87 { << 88 G4String tnames[7] = << 89 {"DEDX","Ionisation","DEDXnr","CSDARange", << 90 } << 91 << 92 << 93 G4VEnergyLossProcess::G4VEnergyLossProcess(con 143 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 94 G4P << 144 G4ProcessType type): G4VContinuousDiscreteProcess(name, type), 95 G4VContinuousDiscreteProcess(name, type) << 145 nSCoffRegions(0), 96 { << 146 idxSCoffRegions(0), 97 theParameters = G4EmParameters::Instance(); << 147 nProcesses(0), 98 SetVerboseLevel(1); << 148 theDEDXTable(0), 99 << 149 theDEDXSubTable(0), 100 // low energy limit << 150 theDEDXunRestrictedTable(0), 101 lowestKinEnergy = theParameters->LowestElect << 151 theIonisationTable(0), >> 152 theIonisationSubTable(0), >> 153 theRangeTableForLoss(0), >> 154 theCSDARangeTable(0), >> 155 theSecondaryRangeTable(0), >> 156 theInverseRangeTable(0), >> 157 theLambdaTable(0), >> 158 theSubLambdaTable(0), >> 159 theDEDXAtMaxEnergy(0), >> 160 theRangeAtMaxEnergy(0), >> 161 theEnergyOfCrossSectionMax(0), >> 162 theCrossSectionMax(0), >> 163 particle(0), >> 164 baseParticle(0), >> 165 secondaryParticle(0), >> 166 currentCouple(0), >> 167 nBins(90), >> 168 nBinsCSDA(70), >> 169 nWarnings(0), >> 170 linLossLimit(0.05), >> 171 minSubRange(0.1), >> 172 lambdaFactor(0.8), >> 173 mfpKinEnergy(0.0), >> 174 lossFluctuationFlag(true), >> 175 lossFluctuationArePossible(true), >> 176 rndmStepFlag(false), >> 177 tablesAreBuilt(false), >> 178 integral(true), >> 179 meanFreePath(false), >> 180 aboveCSmax(true), >> 181 isIonisation(true), >> 182 useSubCutoff(false) >> 183 { >> 184 >> 185 lowestKinEnergy = 1.*eV; >> 186 minKinEnergy = 0.1*keV; >> 187 maxKinEnergy = 100.0*GeV; >> 188 maxKinEnergyCSDA = 1.0*GeV; 102 189 103 // Size of tables << 190 pParticleChange = &fParticleChange; 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 191 112 // default linear loss limit << 192 // default dRoverRange and finalRange 113 finalRange = 1.*CLHEP::mm; << 193 SetStepFunction(0.2, 1.0*mm); >> 194 SetVerboseLevel(1); >> 195 thePositron = G4Positron::Positron(); 114 196 115 // run time objects << 116 pParticleChange = &fParticleChange; << 117 fParticleChange.SetSecondaryWeightByProcess( << 118 modelManager = new G4EmModelManager(); 197 modelManager = new G4EmModelManager(); 119 safetyHelper = G4TransportationManager::GetT << 198 (G4LossTableManager::Instance())->Register(this); 120 ->GetSafetyHelper(); << 199 emModel[0]=emModel[1]=emModel[2]=emModel[3]=emModel[4]=0; 121 aGPILSelection = CandidateForSelection; << 200 fluctModel = 0; 122 << 201 scoffRegions.clear(); 123 // initialise model << 202 scProcesses.clear(); 124 lManager = G4LossTableManager::Instance(); << 203 125 lManager->Register(this); << 204 safetyHelper = new G4SafetyHelper(); 126 isMaster = lManager->IsMaster(); << 205 127 << 206 const G4int n = 7; 128 G4LossTableBuilder* bld = lManager->GetTable << 207 vstrag = new G4PhysicsLogVector(keV, GeV, n); 129 theDensityFactor = bld->GetDensityFactors(); << 208 G4double s[n] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9}; 130 theDensityIdx = bld->GetCoupleIndexes(); << 209 for(G4int i=0; i<n; i++) {vstrag->PutValue(i, s[i]);} 131 << 132 scTracks.reserve(10); << 133 secParticles.reserve(12); << 134 emModels = new std::vector<G4VEmModel*>; << 135 } 210 } 136 211 137 //....oooOO0OOooo........oooOO0OOooo........oo 212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 138 213 139 G4VEnergyLossProcess::~G4VEnergyLossProcess() 214 G4VEnergyLossProcess::~G4VEnergyLossProcess() 140 { 215 { 141 if (isMaster) { << 216 delete vstrag; 142 if(nullptr == baseParticle) { delete theDa << 217 delete safetyHelper; 143 delete theEnergyOfCrossSectionMax; << 218 Clear(); 144 if(nullptr != fXSpeaks) { << 219 145 for(auto const & v : *fXSpeaks) { delete << 220 if ( !baseParticle ) { 146 delete fXSpeaks; << 221 if(theDEDXTable && theRangeTableForLoss) { 147 } << 222 theDEDXTable->clearAndDestroy(); >> 223 if(theDEDXSubTable) theDEDXSubTable->clearAndDestroy(); >> 224 } >> 225 if(theIonisationTable) theIonisationTable->clearAndDestroy(); >> 226 if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy(); >> 227 if(theDEDXunRestrictedTable && theCSDARangeTable) >> 228 theDEDXunRestrictedTable->clearAndDestroy(); >> 229 if(theCSDARangeTable) theCSDARangeTable->clearAndDestroy(); >> 230 if(theRangeTableForLoss) theRangeTableForLoss->clearAndDestroy(); >> 231 if(theInverseRangeTable) theInverseRangeTable->clearAndDestroy(); >> 232 if(theLambdaTable) theLambdaTable->clearAndDestroy(); >> 233 if(theSubLambdaTable) theSubLambdaTable->clearAndDestroy(); 148 } 234 } >> 235 149 delete modelManager; 236 delete modelManager; 150 delete biasManager; << 237 (G4LossTableManager::Instance())->DeRegister(this); 151 delete scoffRegions; << 152 delete emModels; << 153 lManager->DeRegister(this); << 154 } 238 } 155 239 156 //....oooOO0OOooo........oooOO0OOooo........oo 240 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 157 241 158 G4double G4VEnergyLossProcess::MinPrimaryEnerg << 242 void G4VEnergyLossProcess::Clear() 159 << 160 << 161 { 243 { 162 return cut; << 244 if(1 < verboseLevel) 163 } << 245 G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName() >> 246 << G4endl; 164 247 165 //....oooOO0OOooo........oooOO0OOooo........oo << 248 if(theDEDXAtMaxEnergy) delete [] theDEDXAtMaxEnergy; >> 249 if(theRangeAtMaxEnergy) delete [] theRangeAtMaxEnergy; >> 250 if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax; >> 251 if(theCrossSectionMax) delete [] theCrossSectionMax; >> 252 if(idxSCoffRegions) delete [] idxSCoffRegions; 166 253 167 void G4VEnergyLossProcess::AddEmModel(G4int or << 254 theDEDXAtMaxEnergy = 0; 168 G4VEmFlu << 255 theRangeAtMaxEnergy = 0; 169 const G4 << 256 theEnergyOfCrossSectionMax = 0; 170 { << 257 theCrossSectionMax = 0; 171 if(nullptr == ptr) { return; } << 258 tablesAreBuilt = false; 172 G4VEmFluctuationModel* afluc = (nullptr == f << 173 modelManager->AddEmModel(order, ptr, afluc, << 174 ptr->SetParticleChange(pParticleChange, aflu << 175 } << 176 << 177 //....oooOO0OOooo........oooOO0OOooo........oo << 178 259 179 void G4VEnergyLossProcess::SetEmModel(G4VEmMod << 260 scTracks.clear(); 180 { << 261 scProcesses.clear(); 181 if(nullptr == ptr) { return; } << 182 if(!emModels->empty()) { << 183 for(auto & em : *emModels) { if(em == ptr) << 184 } << 185 emModels->push_back(ptr); << 186 } 262 } 187 263 188 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 189 265 190 void G4VEnergyLossProcess::SetDynamicMassCharg << 266 void G4VEnergyLossProcess::PreparePhysicsTable( 191 << 267 const G4ParticleDefinition& part) 192 { 268 { 193 massRatio = massratio; << 194 logMassRatio = G4Log(massRatio); << 195 fFactor = charge2ratio*biasFactor; << 196 if(baseMat) { fFactor *= (*theDensityFactor) << 197 chargeSqRatio = charge2ratio; << 198 reduceFactor = 1.0/(fFactor*massRatio); << 199 } << 200 269 201 //....oooOO0OOooo........oooOO0OOooo........oo << 270 // Are particle defined? >> 271 if( !particle ) { >> 272 if(part.GetParticleType() == "nucleus" && >> 273 part.GetParticleSubType() == "generic") >> 274 particle = G4GenericIon::GenericIon(); >> 275 else particle = ∂ >> 276 } 202 277 203 void << 278 if(1 < verboseLevel) { 204 G4VEnergyLossProcess::PreparePhysicsTable(cons << 279 G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for " 205 { << 280 << GetProcessName() 206 particle = G4EmTableUtil::CheckIon(this, &pa << 281 << " for " << part.GetParticleName() 207 verboseLe << 282 << " local: " << particle->GetParticleName() >> 283 << G4endl; >> 284 } 208 285 209 if( particle != &part ) { << 286 G4LossTableManager* lManager = G4LossTableManager::Instance(); 210 if(!isIon) { lManager->RegisterExtraPartic << 287 211 if(1 < verboseLevel) { << 288 if(&part != particle) { 212 G4cout << "### G4VEnergyLossProcess::Pre << 289 if(part.GetParticleType() == "nucleus") lManager->RegisterIon(&part, this); 213 << " interrupted for " << GetProc << 290 else lManager->RegisterExtraParticle(&part, this); 214 << part.GetParticleName() << " is << 215 << " spline=" << spline << G4endl << 216 } << 217 return; 291 return; 218 } 292 } 219 293 220 tablesAreBuilt = false; << 294 Clear(); 221 if (GetProcessSubType() == fIonisation) { Se << 222 295 223 G4LossTableBuilder* bld = lManager->GetTable << 296 currentCouple = 0; 224 lManager->PreparePhysicsTable(&part, this); << 297 preStepLambda = 0.0; >> 298 mfpKinEnergy = DBL_MAX; >> 299 preStepMFP = DBL_MAX; >> 300 fRange = DBL_MAX; 225 301 226 // Base particle and set of models can be de 302 // Base particle and set of models can be defined here 227 InitialiseEnergyLossProcess(particle, basePa 303 InitialiseEnergyLossProcess(particle, baseParticle); 228 304 229 // parameters of the process << 305 // Tables preparation 230 if(!actLossFluc) { lossFluctuationFlag = the << 306 if (!baseParticle) { 231 useCutAsFinalRange = theParameters->UseCutAs << 307 232 if(!actMinKinEnergy) { minKinEnergy = thePar << 308 theDEDXTable = G4PhysicsTableHelper::PreparePhysicsTable(theDEDXTable); 233 if(!actMaxKinEnergy) { maxKinEnergy = thePar << 309 if (lManager->BuildCSDARange()) { 234 if(!actBinning) { nBins = theParameters->Num << 310 theDEDXunRestrictedTable = 235 maxKinEnergyCSDA = theParameters->MaxEnergyF << 311 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXunRestrictedTable); 236 nBinsCSDA = theParameters->NumberOfBinsPerDe << 312 theCSDARangeTable = 237 *G4lrint(std::log10(maxKinEnergyCSDA/minKi << 313 G4PhysicsTableHelper::PreparePhysicsTable(theCSDARangeTable); 238 if(!actLinLossLimit) { linLossLimit = thePar << 314 } 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 << 246 theParameters->DefineRegParamForLoss(this); << 247 315 248 fRangeEnergy = 0.0; << 316 theRangeTableForLoss = >> 317 G4PhysicsTableHelper::PreparePhysicsTable(theRangeTableForLoss); >> 318 theInverseRangeTable = >> 319 G4PhysicsTableHelper::PreparePhysicsTable(theInverseRangeTable); >> 320 >> 321 theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable); >> 322 if (nSCoffRegions) { >> 323 theDEDXSubTable = >> 324 G4PhysicsTableHelper::PreparePhysicsTable(theDEDXSubTable); >> 325 theSubLambdaTable = >> 326 G4PhysicsTableHelper::PreparePhysicsTable(theSubLambdaTable); >> 327 } >> 328 } 249 329 250 G4double initialCharge = particle->GetPDGCha 330 G4double initialCharge = particle->GetPDGCharge(); 251 G4double initialMass = particle->GetPDGMas 331 G4double initialMass = particle->GetPDGMass(); >> 332 chargeSquare = initialCharge*initialCharge/(eplus*eplus); >> 333 chargeSqRatio = 1.0; >> 334 massRatio = 1.0; >> 335 reduceFactor = 1.0; 252 336 253 theParameters->FillStepFunction(particle, th << 337 if (baseParticle) { 254 << 338 massRatio = (baseParticle->GetPDGMass())/initialMass; 255 // parameters for scaling from the base part << 256 if (nullptr != baseParticle) { << 257 massRatio = (baseParticle->GetPDGMass() << 258 logMassRatio = G4Log(massRatio); << 259 G4double q = initialCharge/baseParticle->G 339 G4double q = initialCharge/baseParticle->GetPDGCharge(); 260 chargeSqRatio = q*q; 340 chargeSqRatio = q*q; 261 if(chargeSqRatio > 0.0) { reduceFactor = 1 << 341 reduceFactor = 1.0/(chargeSqRatio*massRatio); 262 } 342 } 263 lowestKinEnergy = (initialMass < CLHEP::MeV) << 264 ? theParameters->LowestElectronEnergy() << 265 : theParameters->LowestMuHadEnergy(); << 266 343 267 // Tables preparation << 344 theCuts = modelManager->Initialise(particle, secondaryParticle, 268 if (isMaster && nullptr == baseParticle) { << 345 minSubRange, verboseLevel); 269 if(nullptr == theData) { theData = new G4E << 270 346 271 if(nullptr != theDEDXTable && isIonisation << 347 // Sub Cutoff Regime 272 if(nullptr != theIonisationTable && theD << 348 scProcesses.clear(); 273 theData->CleanTable(0); << 349 nProcesses = 0; 274 theDEDXTable = theIonisationTable; << 350 275 theIonisationTable = nullptr; << 351 if (nSCoffRegions>0) { >> 352 theSubCuts = modelManager->SubCutoff(); >> 353 >> 354 const G4ProductionCutsTable* theCoupleTable= >> 355 G4ProductionCutsTable::GetProductionCutsTable(); >> 356 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 357 idxSCoffRegions = new G4int[numOfCouples]; >> 358 >> 359 for (size_t j=0; j<numOfCouples; j++) { >> 360 >> 361 const G4MaterialCutsCouple* couple = >> 362 theCoupleTable->GetMaterialCutsCouple(j); >> 363 const G4ProductionCuts* pcuts = couple->GetProductionCuts(); >> 364 G4int reg = 0; >> 365 for(G4int i=0; i<nSCoffRegions; i++) { >> 366 if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1; 276 } 367 } 277 } << 368 idxSCoffRegions[j] = reg; 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 << 288 theLambdaTable = theData->MakeTable(4); << 289 if(isIonisation) { << 290 theRangeTableForLoss = theData->MakeTabl << 291 theInverseRangeTable = theData->MakeTabl << 292 } 369 } 293 } 370 } 294 371 295 // forced biasing << 372 lManager->EnergyLossProcessIsInitialised(particle, this); 296 if(nullptr != biasManager) { << 297 biasManager->Initialise(part,GetProcessNam << 298 biasFlag = false; << 299 } << 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 } << 319 } << 320 373 321 if(1 < verboseLevel) { << 374 if (1 < verboseLevel) { 322 G4cout << "G4VEnergyLossProcess::PrepearPh << 375 G4cout << "G4VEnergyLossProcess::Initialise() is done " 323 << " for " << GetProcessName() << " << 376 << " 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 377 << " massRatio= " << massRatio 331 << " reduceFactor= " << reduceFacto 378 << " reduceFactor= " << reduceFactor << G4endl; 332 if (nSCoffRegions > 0) { << 379 if (nSCoffRegions) { 333 G4cout << " SubCut secondary production << 380 G4cout << " SubCutoff Regime is ON for regions: " << G4endl; 334 for (G4int i=0; i<nSCoffRegions; ++i) { << 381 for (G4int i=0; i<nSCoffRegions; i++) { 335 const G4Region* r = (*scoffRegions)[i] << 382 const G4Region* r = scoffRegions[i]; 336 G4cout << " " << r->GetName( << 383 G4cout << " " << r->GetName() << G4endl; 337 } 384 } 338 } else if(nullptr != subcutProducer) { << 339 G4cout << " SubCut secondary production << 340 } 385 } 341 } 386 } 342 } 387 } 343 388 344 //....oooOO0OOooo........oooOO0OOooo........oo 389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 345 390 346 void G4VEnergyLossProcess::BuildPhysicsTable(c 391 void G4VEnergyLossProcess::BuildPhysicsTable(const G4ParticleDefinition& part) 347 { 392 { 348 if(1 < verboseLevel) { << 393 if(1 < verboseLevel) { 349 G4cout << "### G4VEnergyLossProcess::Build 394 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() for " 350 << GetProcessName() 395 << GetProcessName() 351 << " and particle " << part.GetPart 396 << " and particle " << part.GetParticleName() 352 << "; the first particle " << parti << 397 << "; local: " << particle->GetParticleName(); 353 if(baseParticle) { << 398 if(baseParticle) G4cout << "; base: " << baseParticle->GetParticleName(); 354 G4cout << "; base: " << baseParticle->Ge << 355 } << 356 G4cout << G4endl; 399 G4cout << G4endl; 357 G4cout << " TablesAreBuilt= " << tables << 358 << " spline=" << spline << " ptr: " << 359 } 400 } 360 401 361 if(&part == particle) { << 402 if(!tablesAreBuilt && &part == particle) 362 if(isMaster) { << 403 G4LossTableManager::Instance()->BuildPhysicsTable(particle, this); 363 lManager->BuildPhysicsTable(particle, th << 404 364 << 405 if(0 < verboseLevel && (&part == particle) && !baseParticle) 365 } else { << 406 PrintInfoDefinition(); 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 } << 388 407 389 // protection against double printout << 390 if(theParameters->IsPrintLocked()) { return; << 391 << 392 // explicitly defined printout by particle n << 393 G4String num = part.GetParticleName(); << 394 if(1 < verboseLevel || << 395 (0 < verboseLevel && (num == "e-" || << 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 } << 404 if(1 < verboseLevel) { 408 if(1 < verboseLevel) { 405 G4cout << "### G4VEnergyLossProcess::Build 409 G4cout << "### G4VEnergyLossProcess::BuildPhysicsTable() done for " 406 << GetProcessName() 410 << GetProcessName() 407 << " and particle " << part.GetPart 411 << " and particle " << part.GetParticleName(); 408 if(isIonisation) { G4cout << " isIonisati << 412 if(isIonisation) G4cout << " isIonisation flag = 1"; 409 G4cout << " baseMat=" << baseMat << G4endl << 413 G4cout << G4endl; 410 } 414 } 411 } 415 } 412 416 413 //....oooOO0OOooo........oooOO0OOooo........oo 417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 414 418 415 G4PhysicsTable* G4VEnergyLossProcess::BuildDED << 419 void G4VEnergyLossProcess::SetEmModel(G4VEmModel* p, G4int index) 416 { 420 { 417 G4PhysicsTable* table = nullptr; << 421 emModel[index] = p; 418 G4double emax = maxKinEnergy; << 419 G4int bin = nBins; << 420 << 421 if(fTotal == tType) { << 422 emax = maxKinEnergyCSDA; << 423 bin = nBinsCSDA; << 424 table = theDEDXunRestrictedTable; << 425 } else if(fRestricted == tType) { << 426 table = theDEDXTable; << 427 } else { << 428 G4cout << "G4VEnergyLossProcess::BuildDEDX << 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 } 422 } 445 423 446 //....oooOO0OOooo........oooOO0OOooo........oo 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 447 425 448 G4PhysicsTable* G4VEnergyLossProcess::BuildLam << 426 G4VEmModel* G4VEnergyLossProcess::EmModel(G4int index) 449 { 427 { 450 if(nullptr == theLambdaTable) { return theLa << 428 return emModel[index]; 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 } 429 } 464 430 465 //....oooOO0OOooo........oooOO0OOooo........oo 431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 466 432 467 void G4VEnergyLossProcess::StreamInfo(std::ost << 433 void G4VEnergyLossProcess::SetFluctModel(G4VEmFluctuationModel* p) 468 const G4ParticleDefinition& pa << 469 { 434 { 470 G4String indent = (rst ? " " : ""); << 435 fluctModel = p; 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 } << 511 } 436 } 512 437 513 //....oooOO0OOooo........oooOO0OOooo........oo 438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 514 439 515 void G4VEnergyLossProcess::ActivateSubCutoff(c << 440 G4VEmFluctuationModel* G4VEnergyLossProcess::FluctModel() 516 { 441 { 517 if(nullptr == scoffRegions) { << 442 return fluctModel; 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 } 443 } 530 444 531 //....oooOO0OOooo........oooOO0OOooo........oo 445 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 532 446 533 G4bool G4VEnergyLossProcess::IsRegionForCubcut << 447 void G4VEnergyLossProcess::AddEmModel(G4int order, G4VEmModel* p, >> 448 G4VEmFluctuationModel* fluc, >> 449 const G4Region* region) 534 { 450 { 535 if(0 == nSCoffRegions) { return true; } << 451 modelManager->AddEmModel(order, p, fluc, region); 536 const G4Region* r = aTrack.GetVolume()->GetL << 452 if(p) p->SetParticleChange(pParticleChange, fluc); 537 for(auto & reg : *scoffRegions) { << 453 if(!fluc) { 538 if(r == reg) { return true; } << 454 lossFluctuationFlag = false; >> 455 lossFluctuationArePossible = false; 539 } 456 } 540 return false; << 541 } 457 } 542 458 543 //....oooOO0OOooo........oooOO0OOooo........oo 459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 544 460 545 void G4VEnergyLossProcess::StartTracking(G4Tra << 461 void G4VEnergyLossProcess::UpdateEmModel(const G4String& nam, >> 462 G4double emin, G4double emax) 546 { 463 { 547 // reset parameters for the new track << 464 modelManager->UpdateEmModel(nam, emin, emax); 548 theNumberOfInteractionLengthLeft = -1.0; << 549 mfpKinEnergy = DBL_MAX; << 550 preStepLambda = 0.0; << 551 currentCouple = nullptr; << 552 << 553 // reset ion << 554 if(isIon) { << 555 const G4double newmass = track->GetDefinit << 556 massRatio = (nullptr == baseParticle) ? CL << 557 : baseParticle->GetPDGMass()/newmass; << 558 logMassRatio = G4Log(massRatio); << 559 } << 560 // forced biasing only for primary particles << 561 if(nullptr != biasManager) { << 562 if(0 == track->GetParentID()) { << 563 biasFlag = true; << 564 biasManager->ResetForcedInteraction(); << 565 } << 566 } << 567 } 465 } 568 466 569 //....oooOO0OOooo........oooOO0OOooo........oo 467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 570 468 571 G4double G4VEnergyLossProcess::AlongStepGetPhy << 469 void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r) 572 const G4Track& tr << 470 { 573 G4GPILSelection* << 471 G4RegionStore* regionStore = G4RegionStore::GetInstance(); 574 { << 472 if(val) { 575 G4double x = DBL_MAX; << 473 useSubCutoff = true; 576 *selection = aGPILSelection; << 474 if (!r) r = regionStore->GetRegion("DefaultRegionForTheWorld", false); 577 if(isIonisation && currentModel->IsActive(pr << 475 if (nSCoffRegions) { 578 GetScaledRangeForScaledEnergy(preStepScale << 476 for (G4int i=0; i<nSCoffRegions; i++) { 579 x = (useCutAsFinalRange) ? std::min(finalR << 477 if (r == scoffRegions[i]) return; 580 currentCouple->GetProductionCuts()->GetP << 478 } 581 x = (fRange > x) ? fRange*dRoverRange + x* << 479 } 582 : fRange; << 480 scoffRegions.push_back(r); 583 /* << 481 nSCoffRegions++; 584 G4cout<<"AlongStepGPIL: " << GetProcessN << 482 } else { 585 << " fRange=" << fRange << " finR=" << finR << 483 useSubCutoff = false; 586 */ << 587 } 484 } 588 return x; << 589 } 485 } 590 486 591 //....oooOO0OOooo........oooOO0OOooo........oo 487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 592 488 593 G4double G4VEnergyLossProcess::PostStepGetPhys << 489 G4PhysicsTable* G4VEnergyLossProcess::BuildDEDXTable(G4EmTableType tType) 594 const G4Track& tr << 595 G4double previo << 596 G4ForceCondition* << 597 { 490 { 598 // condition is set to "Not Forced" << 491 if(1 < verboseLevel) { 599 *condition = NotForced; << 492 G4cout << "G4VEnergyLossProcess::BuildDEDXTable() of type " << tType 600 G4double x = DBL_MAX; << 493 << " for " << GetProcessName() >> 494 << " and particle " << particle->GetParticleName() >> 495 << G4endl; >> 496 } >> 497 G4PhysicsTable* table = 0; >> 498 G4double emin = minKinEnergy; >> 499 G4double emax = maxKinEnergy; >> 500 G4int bin = nBins; >> 501 >> 502 if(fTotal == tType) { >> 503 emax = maxKinEnergyCSDA; >> 504 bin = nBinsCSDA; >> 505 table = theDEDXunRestrictedTable; >> 506 } else if(fRestricted == tType) { >> 507 table = theDEDXTable; >> 508 if(theIonisationTable) >> 509 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable); >> 510 } else if(fSubRestricted == tType) { >> 511 table = theDEDXSubTable; >> 512 if(theIonisationSubTable) >> 513 table = G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable); >> 514 } else { >> 515 G4cout << "G4VEnergyLossProcess::BuildDEDXTable WARNING: wrong type " >> 516 << tType << G4endl; >> 517 } 601 518 602 // initialisation of material, mass, charge, << 519 // Access to materials 603 // at the beginning of the step << 520 const G4ProductionCutsTable* theCoupleTable= 604 DefineMaterial(track.GetMaterialCutsCouple() << 521 G4ProductionCutsTable::GetProductionCutsTable(); 605 preStepKinEnergy = track.GetKineticEne << 522 size_t numOfCouples = theCoupleTable->GetTableSize(); 606 preStepScaledEnergy = preStepKinEnergy*ma << 607 SelectModel(preStepScaledEnergy); << 608 523 609 if(!currentModel->IsActive(preStepScaledEner << 524 if(1 < verboseLevel) { 610 theNumberOfInteractionLengthLeft = -1.0; << 525 G4cout << numOfCouples << " materials" 611 mfpKinEnergy = DBL_MAX; << 526 << " minKinEnergy= " << minKinEnergy 612 preStepLambda = 0.0; << 527 << " maxKinEnergy= " << maxKinEnergy 613 currentInteractionLength = DBL_MAX; << 528 << " EmTableType= " << tType 614 return x; << 529 << " table= " << table >> 530 << G4endl; 615 } 531 } >> 532 if(!table) return table; >> 533 >> 534 for(size_t i=0; i<numOfCouples; i++) { 616 535 617 // change effective charge of a charged part << 536 if(1 < verboseLevel) 618 if(isIon) { << 537 G4cout << "G4VEnergyLossProcess::BuildDEDXVector flag= " 619 const G4double q2 = currentModel->ChargeSq << 538 << table->GetFlag(i) << G4endl; 620 fFactor = q2*biasFactor; << 539 621 if(baseMat) { fFactor *= (*theDensityFacto << 540 if (table->GetFlag(i)) { 622 reduceFactor = 1.0/(fFactor*massRatio); << 541 623 if (lossFluctuationFlag) { << 542 // create physics vector and fill it 624 auto fluc = currentModel->GetModelOfFluc << 543 const G4MaterialCutsCouple* couple = 625 fluc->SetParticleAndCharge(track.GetDefi << 544 theCoupleTable->GetMaterialCutsCouple(i); >> 545 G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin); >> 546 modelManager->FillDEDXVector(aVector, couple, tType); >> 547 >> 548 // Insert vector for this material into the table >> 549 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); 626 } 550 } 627 } 551 } 628 552 629 // forced biasing only for primary particles << 553 if(1 < verboseLevel) { 630 if(biasManager) { << 554 G4cout << "G4VEnergyLossProcess::BuildDEDXTable(): table is built for " 631 if(0 == track.GetParentID() && biasFlag && << 555 << particle->GetParticleName() 632 biasManager->ForcedInteractionRegion((G << 556 << " and process " << GetProcessName() 633 return biasManager->GetStepLimit((G4int) << 557 << G4endl; 634 } << 558 // if(2 < verboseLevel) G4cout << (*table) << G4endl; 635 } 559 } 636 560 637 ComputeLambdaForScaledEnergy(preStepScaledEn << 561 return table; 638 << 562 } 639 // zero cross section << 563 640 if(preStepLambda <= 0.0) { << 564 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 641 theNumberOfInteractionLengthLeft = -1.0; << 565 642 currentInteractionLength = DBL_MAX; << 566 G4PhysicsTable* G4VEnergyLossProcess::BuildLambdaTable(G4EmTableType tType) >> 567 { >> 568 G4PhysicsTable* table = 0; >> 569 >> 570 if(fRestricted == tType) { >> 571 table = theLambdaTable; >> 572 } else if(fSubRestricted == tType) { >> 573 table = theSubLambdaTable; 643 } else { 574 } else { >> 575 G4cout << "G4VEnergyLossProcess::BuildLambdaTable WARNING: wrong type " >> 576 << tType << G4endl; >> 577 } 644 578 645 // non-zero cross section << 579 if(1 < verboseLevel) { 646 if (theNumberOfInteractionLengthLeft < 0.0 << 580 G4cout << "G4VEnergyLossProcess::BuildLambdaTable() of type " >> 581 << tType << " for process " >> 582 << GetProcessName() << " and particle " >> 583 << particle->GetParticleName() >> 584 << " EmTableType= " << tType >> 585 << " table= " << table >> 586 << G4endl; >> 587 } >> 588 if(!table) return table; >> 589 >> 590 // Access to materials >> 591 const G4ProductionCutsTable* theCoupleTable= >> 592 G4ProductionCutsTable::GetProductionCutsTable(); >> 593 size_t numOfCouples = theCoupleTable->GetTableSize(); >> 594 >> 595 for(size_t i=0; i<numOfCouples; i++) { >> 596 >> 597 if (table->GetFlag(i)) { >> 598 >> 599 // create physics vector and fill it >> 600 const G4MaterialCutsCouple* couple = >> 601 theCoupleTable->GetMaterialCutsCouple(i); >> 602 G4double cut = (*theCuts)[i]; >> 603 if(fSubRestricted == tType) cut = (*theSubCuts)[i]; >> 604 G4PhysicsVector* aVector = LambdaPhysicsVector(couple, cut); >> 605 modelManager->FillLambdaVector(aVector, couple, true, tType); 647 606 648 // beggining of tracking (or just after << 607 // Insert vector for this material into the table 649 theNumberOfInteractionLengthLeft = -G4Lo << 608 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); 650 theInitialNumberOfInteractionLength = th << 651 << 652 } else if(currentInteractionLength < DBL_M << 653 << 654 // subtract NumberOfInteractionLengthLef << 655 theNumberOfInteractionLengthLeft -= << 656 previousStepSize/currentInteractionLen << 657 << 658 theNumberOfInteractionLengthLeft = << 659 std::max(theNumberOfInteractionLengthL << 660 } << 661 << 662 // new mean free path and step limit << 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 } << 681 << 682 //....oooOO0OOooo........oooOO0OOooo........oo << 683 << 684 void << 685 G4VEnergyLossProcess::ComputeLambdaForScaledEn << 686 { << 687 // cross section increased with energy << 688 if(fXSType == fEmIncreasing) { << 689 if(e*invLambdaFactor < mfpKinEnergy) { << 690 preStepLambda = GetLambdaForScaledEnergy << 691 mfpKinEnergy = (preStepLambda > 0.0) ? e << 692 } << 693 << 694 // cross section has one peak << 695 } else if(fXSType == fEmOnePeak) { << 696 const G4double epeak = (*theEnergyOfCrossS << 697 if(e <= epeak) { << 698 if(e*invLambdaFactor < mfpKinEnergy) { << 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 << 708 // cross section has more than one peaks << 709 } else if(fXSType == fEmTwoPeaks) { << 710 G4TwoPeaksXS* xs = (*fXSpeaks)[basedCouple << 711 const G4double e1peak = xs->e1peak; << 712 << 713 // below the 1st peak << 714 if(e <= e1peak) { << 715 if(e*invLambdaFactor < mfpKinEnergy) { << 716 preStepLambda = GetLambdaForScaledEner << 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 } 609 } 765 // integral method is not used << 766 } else { << 767 preStepLambda = GetLambdaForScaledEnergy(e << 768 } 610 } >> 611 >> 612 if(1 < verboseLevel) { >> 613 G4cout << "Lambda table is built for " >> 614 << particle->GetParticleName() >> 615 << G4endl; >> 616 } >> 617 >> 618 return table; 769 } 619 } 770 620 771 //....oooOO0OOooo........oooOO0OOooo........oo 621 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 772 622 773 G4VParticleChange* G4VEnergyLossProcess::Along 623 G4VParticleChange* G4VEnergyLossProcess::AlongStepDoIt(const G4Track& track, 774 624 const G4Step& step) 775 { 625 { 776 fParticleChange.InitializeForAlongStep(track 626 fParticleChange.InitializeForAlongStep(track); 777 // The process has range table - calculate e 627 // The process has range table - calculate energy loss 778 if(!isIonisation || !currentModel->IsActive( << 628 if(!isIonisation) return &fParticleChange; 779 return &fParticleChange; << 780 } << 781 629 >> 630 // Get the actual (true) Step length 782 G4double length = step.GetStepLength(); 631 G4double length = step.GetStepLength(); >> 632 if(length == 0.0) return &fParticleChange; 783 G4double eloss = 0.0; 633 G4double eloss = 0.0; 784 << 634 785 /* << 635 /* 786 if(-1 < verboseLevel) { 636 if(-1 < verboseLevel) { 787 const G4ParticleDefinition* d = track.GetP << 637 const G4ParticleDefinition* d = track.GetDefinition(); 788 G4cout << "AlongStepDoIt for " 638 G4cout << "AlongStepDoIt for " 789 << GetProcessName() << " and partic << 639 << GetProcessName() << " and particle " 790 << " eScaled(MeV)=" << preStepScal << 640 << d->GetParticleName() 791 << " range(mm)=" << fRange/mm << " << 641 << " eScaled(MeV)= " << preStepScaledEnergy/MeV 792 << " rf=" << reduceFactor << " q^ << 642 << " range(mm)= " << fRange/mm 793 << " md=" << d->GetPDGMass() << " << 643 << " s(mm)= " << length/mm 794 << " " << track.GetMaterial()->Get << 644 << " q^2= " << chargeSqRatio >> 645 << " md= " << d->GetPDGMass() >> 646 << " status= " << track.GetTrackStatus() >> 647 << G4endl; 795 } 648 } 796 */ 649 */ 797 const G4DynamicParticle* dynParticle = track << 798 650 799 // define new weight for primary and seconda << 651 // stopping 800 G4double weight = fParticleChange.GetParentW << 652 if (length >= fRange) { 801 if(weightFlag) { << 802 weight /= biasFactor; << 803 fParticleChange.ProposeWeight(weight); << 804 } << 805 << 806 // stopping, check actual range and kinetic << 807 if (length >= fRange || preStepKinEnergy <= << 808 eloss = preStepKinEnergy; 653 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 654 822 // Short step 655 // Short step 823 eloss = length*GetDEDXForScaledEnergy(preSte << 656 } else if( length <= linLossLimit * fRange ) { 824 LogSca << 657 eloss = GetDEDXForScaledEnergy(preStepScaledEnergy)*length; 825 /* << 658 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 659 // Long step 834 if(eloss > preStepKinEnergy*linLossLimit) { << 660 } else { >> 661 G4double r = GetScaledRangeForScaledEnergy(preStepScaledEnergy); >> 662 G4double x = r - length/reduceFactor; >> 663 if(x < 0.0) { >> 664 if(0 < verboseLevel && nWarnings<0) { >> 665 G4cout << "WARNING! G4VEnergyLossProcess::AlongStepDoIt: x= " << x >> 666 << " for eScaled(MeV)= " << preStepScaledEnergy/MeV >> 667 << " step(mm)= " << length/mm >> 668 << " range(mm)= " << fRange/mm >> 669 << " for " << track.GetDefinition()->GetParticleName() >> 670 << G4endl; >> 671 nWarnings++; >> 672 } >> 673 x = 0.0; >> 674 } >> 675 eloss = (ScaledKinEnergyForLoss(r) - ScaledKinEnergyForLoss(x))/massRatio; 835 676 836 const G4double x = (fRange - length)/reduc << 837 const G4double de = preStepKinEnergy - Sca << 838 if(de > 0.0) { eloss = de; } << 839 /* 677 /* 840 if(-1 < verboseLevel) 678 if(-1 < verboseLevel) 841 G4cout << " Long STEP: rPre(mm)=" << 679 G4cout << "Long STEP: rPre(mm)= " << r/mm 842 << GetScaledRangeForScaledEnergy( << 680 << " rPost(mm)= " << x/mm 843 << " x(mm)=" << x/mm << 681 << " ePre(MeV)= " << preStepScaledEnergy/MeV 844 << " eloss(MeV)=" << eloss/MeV << 682 << " eloss(MeV)= " << eloss/MeV 845 << " rFactor=" << reduceFactor << 683 << " eloss0(MeV)= " 846 << " massRatio=" << massRatio << 684 << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV 847 << G4endl; 685 << G4endl; 848 */ 686 */ 849 } 687 } 850 688 851 /* << 689 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); >> 690 G4VEmModel* currentModel = SelectModel(preStepScaledEnergy); >> 691 /* >> 692 G4double eloss0 = eloss; 852 if(-1 < verboseLevel ) { 693 if(-1 < verboseLevel ) { 853 G4cout << "Before fluct: eloss(MeV)= " << 694 G4cout << "Before fluct: eloss(MeV)= " << eloss/MeV 854 << " e-eloss= " << preStepKinEnergy 695 << " e-eloss= " << preStepKinEnergy-eloss 855 << " step(mm)= " << length/mm << " << 696 << " step(mm)= " << length/mm 856 << " fluct= " << lossFluctuationFla << 697 << " range(mm)= " << fRange/mm >> 698 << " fluct= " << lossFluctuationFlag >> 699 << G4endl; 857 } 700 } 858 */ << 701 */ 859 702 860 const G4double cut = (*theCuts)[currentCoupl << 703 G4double cut = (*theCuts)[currentMaterialIndex]; 861 G4double esec = 0.0; 704 G4double esec = 0.0; >> 705 G4double esecdep = 0.0; 862 706 863 // Corrections, which cannot be tabulated << 707 // SubCutOff 864 if(isIon) { << 708 if(useSubCutoff) { 865 currentModel->CorrectionsAlongStep(current << 709 if(idxSCoffRegions[currentMaterialIndex]) { 866 length, << 710 867 eloss = std::max(eloss, 0.0); << 711 G4double currentMinSafety = 0.0; >> 712 G4StepPoint* prePoint = step.GetPreStepPoint(); >> 713 G4StepPoint* postPoint = step.GetPostStepPoint(); >> 714 G4double preSafety = prePoint->GetSafety(); >> 715 G4double postSafety = preSafety - length; >> 716 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1); >> 717 >> 718 // recompute safety >> 719 if(prePoint->GetStepStatus() != fGeomBoundary && >> 720 postPoint->GetStepStatus() != fGeomBoundary) { >> 721 >> 722 /* >> 723 // G4bool yes = (track.GetTrackID() == 5512); >> 724 G4bool yes = false; >> 725 if(yes) >> 726 G4cout << "G4VEnergyLoss: presafety= " << preSafety >> 727 << " rcut= " << rcut << " length= " << length >> 728 << " dir " << track.GetMomentumDirection() >> 729 << G4endl; >> 730 */ >> 731 >> 732 if(preSafety < rcut) >> 733 preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition()); >> 734 >> 735 //if(yes) { >> 736 // G4cout << "G4VEnergyLoss: newsafety= " << preSafety << G4endl; >> 737 // if(preSafety==0.0 && track.GetTrackID() == 5512 ) exit(1); >> 738 //} >> 739 if(postSafety < rcut) >> 740 postSafety = safetyHelper->ComputeSafety(postPoint->GetPosition()); >> 741 /* >> 742 if(-1 < verboseLevel) >> 743 G4cout << "Subcutoff: presafety(mm)= " << preSafety/mm >> 744 << " postsafety(mm)= " << postSafety/mm >> 745 << " rcut(mm)= " << rcut/mm >> 746 << G4endl; >> 747 */ >> 748 currentMinSafety = std::min(preSafety,postSafety); >> 749 } >> 750 >> 751 // Decide to start subcut sampling >> 752 if(currentMinSafety < rcut) { >> 753 >> 754 cut = (*theSubCuts)[currentMaterialIndex]; >> 755 eloss -= GetSubDEDXForScaledEnergy(preStepScaledEnergy)*length; >> 756 if(eloss < 0.0) eloss = 0.0; >> 757 SampleSubCutSecondaries(scTracks, step, currentModel,currentMaterialIndex); >> 758 /* >> 759 if(nProcesses) { >> 760 for(G4int i=0; i<nProcesses; i++) { >> 761 (scProcesses[i])->SampleSubCutSecondaries( >> 762 scTracks, step, (scProcesses[i])-> >> 763 SelectModelForMaterial(preStepKinEnergy, currentMaterialIndex), >> 764 currentMaterialIndex); >> 765 } >> 766 } >> 767 */ >> 768 G4int n = scTracks.size(); >> 769 if(n) { >> 770 G4ThreeVector mom = dynParticle->GetMomentum(); >> 771 fParticleChange.SetNumberOfSecondaries(n); >> 772 for(G4int i=0; i<n; i++) { >> 773 G4Track* t = scTracks[i]; >> 774 G4double ekin = t->GetKineticEnergy(); >> 775 G4double e = ekin; >> 776 if (t->GetDefinition() == thePositron) e += 2.0*electron_mass_c2; >> 777 G4bool addSec = true; >> 778 >> 779 // do not track very low-energy delta-electrons >> 780 if(theSecondaryRangeTable) { >> 781 G4bool b; >> 782 G4double rg = >> 783 ((*theSecondaryRangeTable)[currentMaterialIndex]->GetValue(e, b)); >> 784 // if(rg < currentMinSafety) { >> 785 if(rg < safetyHelper->ComputeSafety(t->GetPosition())) { >> 786 esecdep += e; >> 787 delete t; >> 788 addSec = false; >> 789 } >> 790 } >> 791 if(addSec) { >> 792 esec += e; >> 793 pParticleChange->AddSecondary(t); >> 794 mom -= t->GetMomentum(); >> 795 } >> 796 } >> 797 scTracks.clear(); >> 798 // fParticleChange.SetProposedMomentum(mom); >> 799 } >> 800 } >> 801 } 868 } 802 } 869 803 870 // Sample fluctuations if not full energy lo << 804 // Corrections, which cannot be tabulated 871 if(eloss >= preStepKinEnergy) { << 805 CorrectionsAlongStep(currentCouple, dynParticle, eloss, length); 872 eloss = preStepKinEnergy; << 873 806 874 } else if (lossFluctuationFlag) { << 807 // Sample fluctuations 875 const G4double tmax = currentModel->MaxSec << 808 if (lossFluctuationFlag && 876 const G4double tcut = std::min(cut, tmax); << 809 eloss + esec + esecdep + lowestKinEnergy < preStepKinEnergy) { 877 G4VEmFluctuationModel* fluc = currentModel << 810 878 eloss = fluc->SampleFluctuations(currentCo << 811 G4double tmax = 879 tcut, tma << 812 std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut); 880 /* << 813 eloss = currentModel->GetModelOfFluctuations()-> 881 if(-1 < verboseLevel) << 814 SampleFluctuations(currentMaterial,dynParticle,tmax,length,eloss); >> 815 >> 816 /* >> 817 if(-1 < verboseLevel) 882 G4cout << "After fluct: eloss(MeV)= " << 818 G4cout << "After fluct: eloss(MeV)= " << eloss/MeV 883 << " fluc= " << (eloss-eloss0)/Me 819 << " fluc= " << (eloss-eloss0)/MeV 884 << " ChargeSqRatio= " << chargeSq << 820 << " currentChargeSquare= " << chargeSquare 885 << " massRatio= " << massRatio << << 821 << " massRatio= " << massRatio >> 822 << " tmax= " << tmax >> 823 << G4endl; 886 */ 824 */ 887 } 825 } >> 826 // add low-energy subcutoff particles >> 827 eloss += esecdep; 888 828 889 // deexcitation << 829 // Energy balanse 890 if (useDeexcitation) { << 891 G4double esecfluo = preStepKinEnergy; << 892 G4double de = esecfluo; << 893 atomDeexcitation->AlongStepDeexcitation(sc << 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 - 830 G4double finalT = preStepKinEnergy - eloss - esec; 916 if (finalT <= lowestKinEnergy) { 831 if (finalT <= lowestKinEnergy) { 917 eloss += finalT; 832 eloss += finalT; 918 finalT = 0.0; 833 finalT = 0.0; 919 } else if(isIon) { << 920 fParticleChange.SetProposedCharge( << 921 currentModel->GetParticleCharge(track.Ge << 922 currentM << 923 } 834 } 924 eloss = std::max(eloss, 0.0); << 925 835 926 fParticleChange.SetProposedKineticEnergy(fin 836 fParticleChange.SetProposedKineticEnergy(finalT); 927 fParticleChange.ProposeLocalEnergyDeposit(el 837 fParticleChange.ProposeLocalEnergyDeposit(eloss); 928 /* << 838 >> 839 /* 929 if(-1 < verboseLevel) { 840 if(-1 < verboseLevel) { 930 G4double del = finalT + eloss + esec - pre << 931 G4cout << "Final value eloss(MeV)= " << el 841 G4cout << "Final value eloss(MeV)= " << eloss/MeV 932 << " preStepKinEnergy= " << preStep 842 << " preStepKinEnergy= " << preStepKinEnergy 933 << " postStepKinEnergy= " << finalT 843 << " postStepKinEnergy= " << finalT 934 << " de(keV)= " << del/keV << 935 << " lossFlag= " << lossFluctuation 844 << " lossFlag= " << lossFluctuationFlag 936 << " status= " << track.GetTrackSt 845 << " status= " << track.GetTrackStatus() 937 << G4endl; 846 << G4endl; 938 } 847 } 939 */ 848 */ >> 849 940 return &fParticleChange; 850 return &fParticleChange; 941 } 851 } 942 852 943 //....oooOO0OOooo........oooOO0OOooo........oo 853 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 944 854 945 void G4VEnergyLossProcess::FillSecondariesAlon << 855 void G4VEnergyLossProcess::SampleSubCutSecondaries( 946 { << 856 std::vector<G4Track*>& tracks, 947 const std::size_t n0 = scTracks.size(); << 857 const G4Step& step, 948 G4double weight = wt; << 858 G4VEmModel* model, 949 // weight may be changed by biasing manager << 859 G4int idx) 950 if(biasManager) { << 860 { 951 if(biasManager->SecondaryBiasingRegion((G4 << 861 // Fast check weather subcutoff can work 952 weight *= << 862 G4double subcut = (*theSubCuts)[idx]; 953 biasManager->ApplySecondaryBiasing(scT << 863 G4double cut = (*theCuts)[idx]; 954 } << 864 if(cut <= subcut) return; 955 } << 865 956 << 866 const G4Track* track = step.GetTrack(); 957 // fill secondaries << 867 const G4DynamicParticle* dp = track->GetDynamicParticle(); 958 const std::size_t n = scTracks.size(); << 868 G4bool b; 959 fParticleChange.SetNumberOfSecondaries((G4in << 869 G4double cross = 960 << 870 chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(dp->GetKineticEnergy(),b)); 961 for(std::size_t i=0; i<n; ++i) { << 871 G4double length = step.GetStepLength(); 962 G4Track* t = scTracks[i]; << 872 if(length*cross < 1.e-9) return; 963 if(nullptr != t) { << 873 /* 964 t->SetWeight(weight); << 874 if(-1 < verboseLevel) 965 pParticleChange->AddSecondary(t); << 875 G4cout << "<<< Subcutoff for " << GetProcessName() 966 G4int pdg = t->GetDefinition()->GetPDGEn << 876 << " cross(1/mm)= " << cross*mm << ">>>" 967 if (i < n0) { << 877 << " e(MeV)= " << preStepScaledEnergy 968 if (pdg == 22) { << 878 << " matIdx= " << currentMaterialIndex 969 t->SetCreatorModelID(gpixeID); << 879 << G4endl; 970 } else if (pdg == 11) { << 880 */ 971 t->SetCreatorModelID(epixeID); << 881 972 } else { << 882 // Sample subcutoff secondaries 973 t->SetCreatorModelID(biasID); << 883 G4StepPoint* preStepPoint = step.GetPreStepPoint(); 974 } << 884 G4ThreeVector prepoint = preStepPoint->GetPosition(); 975 } else { << 885 G4ThreeVector dr = step.GetPostStepPoint()->GetPosition() - prepoint; 976 t->SetCreatorModelID(biasID); << 886 G4double pretime = preStepPoint->GetGlobalTime(); >> 887 G4double dt = length/preStepPoint->GetVelocity(); >> 888 G4double fragment = 0.0; >> 889 >> 890 do { >> 891 G4double del = -std::log(G4UniformRand())/cross; >> 892 fragment += del/length; >> 893 if (fragment > 1.0) break; >> 894 std::vector<G4DynamicParticle*>* newp = >> 895 model->SampleSecondaries(track->GetMaterialCutsCouple(),dp,subcut,cut); >> 896 if (newp) { >> 897 >> 898 G4DynamicParticle* p; >> 899 G4int nNew = newp->size(); >> 900 for (G4int i=0; i<nNew; i++) { >> 901 >> 902 p = (*newp)[i]; >> 903 G4ThreeVector r = prepoint + fragment*dr; >> 904 G4Track* t = new G4Track(p, pretime + fragment*dt, r); >> 905 t->SetTouchableHandle(track->GetTouchableHandle()); >> 906 tracks.push_back(t); >> 907 /* >> 908 if(-1 < verboseLevel) >> 909 G4cout << "New track " << p->GetDefinition()->GetParticleName() >> 910 << " e(keV)= " << p->GetKineticEnergy()/keV >> 911 << " fragment= " << fragment >> 912 << G4endl; >> 913 */ 977 } 914 } >> 915 delete newp; 978 } 916 } 979 } << 917 } while (fragment <= 1.0); 980 scTracks.clear(); << 918 } 981 } << 982 919 983 //....oooOO0OOooo........oooOO0OOooo........oo 920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 984 921 985 G4VParticleChange* G4VEnergyLossProcess::PostS 922 G4VParticleChange* G4VEnergyLossProcess::PostStepDoIt(const G4Track& track, 986 923 const G4Step& step) 987 { 924 { 988 // clear number of interaction lengths in an << 989 theNumberOfInteractionLengthLeft = -1.0; << 990 mfpKinEnergy = DBL_MAX; << 991 << 992 fParticleChange.InitializeForPostStep(track) 925 fParticleChange.InitializeForPostStep(track); 993 const G4double finalT = track.GetKineticEner << 926 G4double finalT = track.GetKineticEnergy(); 994 << 927 if(finalT == 0.0) return &fParticleChange; 995 const G4double postStepScaledEnergy = finalT << 996 SelectModel(postStepScaledEnergy); << 997 928 998 if(!currentModel->IsActive(postStepScaledEne << 929 G4double postStepScaledEnergy = finalT*massRatio; 999 return &fParticleChange; << 1000 } << 1001 /* 930 /* 1002 if(1 < verboseLevel) { << 931 if(-1 < verboseLevel) { 1003 G4cout<<GetProcessName()<<" PostStepDoIt: << 932 G4cout << GetProcessName() >> 933 << "::PostStepDoIt: E(MeV)= " << finalT/MeV >> 934 << G4endl; 1004 } 935 } 1005 */ 936 */ 1006 // forced process - should happen only once << 937 // Integral approach 1007 if(biasFlag) { << 938 if (integral) { 1008 if(biasManager->ForcedInteractionRegion(( << 939 G4double lx = GetLambdaForScaledEnergy(postStepScaledEnergy); 1009 biasFlag = false; << 940 /* >> 941 if(preStepLambda<lx && 1 < verboseLevel && nWarnings<200) { >> 942 G4cout << "WARNING: for " << particle->GetParticleName() >> 943 << " and " << GetProcessName() >> 944 << " E(MeV)= " << finalT/MeV >> 945 << " preLambda= " << preStepLambda >> 946 << " < " << lx << " (postLambda) " >> 947 << G4endl; >> 948 nWarnings++; 1010 } 949 } >> 950 */ >> 951 if(preStepLambda*G4UniformRand() > lx) >> 952 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step); 1011 } 953 } 1012 const G4DynamicParticle* dp = track.GetDyna << 1013 954 1014 // Integral approach << 955 G4VEmModel* currentModel = SelectModel(postStepScaledEnergy); 1015 if (fXSType != fEmNoIntegral) { << 956 const G4DynamicParticle* dynParticle = track.GetDynamicParticle(); 1016 const G4double logFinalT = dp->GetLogKine << 957 G4double tcut = (*theCuts)[currentMaterialIndex]; 1017 G4double lx = GetLambdaForScaledEnergy(po << 958 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle); 1018 lo << 1019 lx = std::max(lx, 0.0); << 1020 959 1021 // if both lg and lx are zero then no int << 960 if (tcut < tmax) { 1022 if(preStepLambda*G4UniformRand() >= lx) { << 961 std::vector<G4DynamicParticle*>* newp = SecondariesPostStep( 1023 return &fParticleChange; << 962 currentModel, currentCouple, dynParticle, tcut); >> 963 >> 964 if(newp) { >> 965 G4int num = newp->size(); >> 966 fParticleChange.SetNumberOfSecondaries(num); >> 967 for (G4int i=0; i<num; i++) { >> 968 fParticleChange.AddSecondary((*newp)[i]); >> 969 } >> 970 delete newp; 1024 } 971 } >> 972 >> 973 finalT = fParticleChange.GetProposedKineticEnergy(); 1025 } 974 } 1026 975 1027 // define new weight for primary and second << 976 /* 1028 G4double weight = fParticleChange.GetParent << 977 if(-1 < verboseLevel) { 1029 if(weightFlag) { << 978 G4cout << "::PostStepDoIt: Sample secondary; Efin= " << finalT/MeV 1030 weight /= biasFactor; << 979 << " MeV; model= (" << currentModel->LowEnergyLimit() 1031 fParticleChange.ProposeWeight(weight); << 980 << ", " << currentModel->HighEnergyLimit() << ")" >> 981 << " preStepLambda= " << preStepLambda >> 982 << " dir= " << track.GetMomentumDirection() >> 983 << " status= " << track.GetTrackStatus() >> 984 << G4endl; >> 985 } >> 986 */ >> 987 if (finalT <= lowestKinEnergy) { >> 988 fParticleChange.SetProposedKineticEnergy(0.0); >> 989 return &fParticleChange; 1032 } 990 } 1033 991 1034 const G4double tcut = (*theCuts)[currentCou << 992 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step); >> 993 } 1035 994 1036 // sample secondaries << 995 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1037 secParticles.clear(); << 1038 currentModel->SampleSecondaries(&secParticl << 1039 996 1040 const G4int num0 = (G4int)secParticles.size << 997 void G4VEnergyLossProcess::PrintInfoDefinition() >> 998 { >> 999 if(0 < verboseLevel) { >> 1000 G4cout << G4endl << GetProcessName() << ": tables are built for " >> 1001 << particle->GetParticleName() >> 1002 << G4endl >> 1003 << " dE/dx and range tables from " >> 1004 << G4BestUnit(minKinEnergy,"Energy") >> 1005 << " to " << G4BestUnit(maxKinEnergy,"Energy") >> 1006 << " in " << nBins << " bins." << G4endl >> 1007 << " Lambda tables from threshold to " >> 1008 << G4BestUnit(maxKinEnergy,"Energy") >> 1009 << " in " << nBins << " bins." >> 1010 << G4endl; >> 1011 PrintInfo(); >> 1012 if(theRangeTableForLoss && isIonisation) >> 1013 G4cout << " Step function: finalRange(mm)= " << finalRange/mm >> 1014 << ", dRoverRange= " << dRoverRange >> 1015 << ", integral: " << integral >> 1016 << G4endl; >> 1017 >> 1018 if(theCSDARangeTable && isIonisation) >> 1019 G4cout << " CSDA range table up" >> 1020 << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy") >> 1021 << " in " << nBinsCSDA << " bins." << G4endl; >> 1022 >> 1023 if(nSCoffRegions>0) >> 1024 G4cout << " Subcutoff sampling in " << nSCoffRegions >> 1025 << " regions" << G4endl; >> 1026 >> 1027 if(2 < verboseLevel) { >> 1028 G4cout << "DEDXTable address= " << theDEDXTable << G4endl; >> 1029 if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl; >> 1030 G4cout << "non restricted DEDXTable address= " >> 1031 << theDEDXunRestrictedTable << G4endl; >> 1032 if(theDEDXunRestrictedTable && isIonisation) >> 1033 G4cout << (*theDEDXunRestrictedTable) << G4endl; >> 1034 if(theDEDXSubTable && isIonisation) G4cout << (*theDEDXSubTable) >> 1035 << G4endl; >> 1036 G4cout << "CSDARangeTable address= " << theCSDARangeTable >> 1037 << G4endl; >> 1038 if(theCSDARangeTable && isIonisation) G4cout << (*theCSDARangeTable) >> 1039 << G4endl; >> 1040 G4cout << "RangeTableForLoss address= " << theRangeTableForLoss >> 1041 << G4endl; >> 1042 if(theRangeTableForLoss && isIonisation) >> 1043 G4cout << (*theRangeTableForLoss) << G4endl; >> 1044 G4cout << "InverseRangeTable address= " << theInverseRangeTable >> 1045 << G4endl; >> 1046 if(theInverseRangeTable && isIonisation) >> 1047 G4cout << (*theInverseRangeTable) << G4endl; >> 1048 G4cout << "LambdaTable address= " << theLambdaTable << G4endl; >> 1049 if(theLambdaTable && isIonisation) G4cout << (*theLambdaTable) << G4endl; >> 1050 G4cout << "SubLambdaTable address= " << theSubLambdaTable << G4endl; >> 1051 if(theSubLambdaTable && isIonisation) G4cout << (*theSubLambdaTable) >> 1052 << G4endl; >> 1053 } >> 1054 } >> 1055 } >> 1056 >> 1057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1041 1058 1042 // bremsstrahlung splitting or Russian roul << 1059 void G4VEnergyLossProcess::SetDEDXTable(G4PhysicsTable* p, G4EmTableType tType) 1043 if(biasManager) { << 1060 { 1044 if(biasManager->SecondaryBiasingRegion((G << 1061 if(fTotal == tType && theDEDXunRestrictedTable != p) { 1045 G4double eloss = 0.0; << 1062 if(theDEDXunRestrictedTable) theDEDXunRestrictedTable->clearAndDestroy(); 1046 weight *= biasManager->ApplySecondaryBi << 1063 theDEDXunRestrictedTable = p; 1047 secPart << 1064 if(p) { 1048 track, << 1065 size_t n = p->length(); 1049 &fParti << 1066 G4PhysicsVector* pv = (*p)[0]; 1050 (G4int) << 1067 G4double emax = maxKinEnergyCSDA; 1051 step.Ge << 1068 G4bool b; 1052 if(eloss > 0.0) { << 1069 theDEDXAtMaxEnergy = new G4double [n]; 1053 eloss += fParticleChange.GetLocalEner << 1070 1054 fParticleChange.ProposeLocalEnergyDep << 1071 for (size_t i=0; i<n; i++) { >> 1072 pv = (*p)[i]; >> 1073 G4double dedx = pv->GetValue(emax, b); >> 1074 theDEDXAtMaxEnergy[i] = dedx; >> 1075 //G4cout << "i= " << i << " emax(MeV)= " << emax/MeV<< " dedx= " >> 1076 //<< dedx << G4endl; 1055 } 1077 } 1056 } 1078 } >> 1079 >> 1080 } else if(fRestricted == tType) { >> 1081 theDEDXTable = p; >> 1082 } else if(fSubRestricted == tType) { >> 1083 theDEDXSubTable = p; >> 1084 } else if(fIonisation == tType && theIonisationTable != p) { >> 1085 if(theIonisationTable) theIonisationTable->clearAndDestroy(); >> 1086 theIonisationTable = p; >> 1087 } else if(fSubIonisation == tType && theIonisationSubTable != p) { >> 1088 if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy(); >> 1089 theIonisationSubTable = p; 1057 } 1090 } >> 1091 } 1058 1092 1059 // save secondaries << 1093 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1060 const G4int num = (G4int)secParticles.size( << 1061 if(num > 0) { << 1062 1094 1063 fParticleChange.SetNumberOfSecondaries(nu << 1095 void G4VEnergyLossProcess::SetCSDARangeTable(G4PhysicsTable* p) 1064 G4double time = track.GetGlobalTime(); << 1096 { >> 1097 if(theCSDARangeTable != p) theCSDARangeTable = p; 1065 1098 1066 G4int n1(0), n2(0); << 1099 if(p) { 1067 if(num0 > mainSecondaries) { << 1100 size_t n = p->length(); 1068 currentModel->FillNumberOfSecondaries(n << 1101 G4PhysicsVector* pv = (*p)[0]; >> 1102 G4double emax = maxKinEnergyCSDA; >> 1103 G4bool b; >> 1104 theRangeAtMaxEnergy = new G4double [n]; >> 1105 >> 1106 for (size_t i=0; i<n; i++) { >> 1107 pv = (*p)[i]; >> 1108 G4double r2 = pv->GetValue(emax, b); >> 1109 theRangeAtMaxEnergy[i] = r2; >> 1110 //G4cout << "i= " << i << " e2(MeV)= " << emax/MeV << " r2= " >> 1111 //<< r2<< G4endl; 1069 } 1112 } >> 1113 } >> 1114 } 1070 1115 1071 for (G4int i=0; i<num; ++i) { << 1116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1072 if(nullptr != secParticles[i]) { << 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 1117 1088 //G4cout << "Secondary(post step) has << 1118 void G4VEnergyLossProcess::SetRangeTableForLoss(G4PhysicsTable* p) 1089 // << ", kenergy " << t->GetKin << 1119 { 1090 // << " time= " << time/ns << " << 1120 if(theRangeTableForLoss != p) { 1091 pParticleChange->AddSecondary(t); << 1121 theRangeTableForLoss = p; 1092 } << 1122 if(1 < verboseLevel) { >> 1123 G4cout << "### Set Range table " << p >> 1124 << " for " << particle->GetParticleName() >> 1125 << " and process " << GetProcessName() << G4endl; 1093 } 1126 } 1094 } 1127 } >> 1128 } 1095 1129 1096 if(0.0 == fParticleChange.GetProposedKineti << 1130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1097 fAlive == fParticleChange.GetTrackStatus << 1098 if(particle->GetProcessManager()->GetAtRe << 1099 { fParticleChange.ProposeTrackStatus << 1100 else { fParticleChange.ProposeTrackStatus << 1101 } << 1102 1131 1103 /* << 1132 void G4VEnergyLossProcess::SetSecondaryRangeTable(G4PhysicsTable* p) 1104 if(-1 < verboseLevel) { << 1133 { 1105 G4cout << "::PostStepDoIt: Sample seconda << 1134 theSecondaryRangeTable = p; 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 } 1135 } 1117 1136 1118 //....oooOO0OOooo........oooOO0OOooo........o 1137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1119 1138 1120 G4bool G4VEnergyLossProcess::StorePhysicsTabl << 1139 void G4VEnergyLossProcess::SetInverseRangeTable(G4PhysicsTable* p) 1121 const G4ParticleDefinition* part, cons << 1122 { 1140 { 1123 if (!isMaster || nullptr != baseParticle || << 1141 if(theInverseRangeTable != p) { 1124 for(std::size_t i=0; i<7; ++i) { << 1142 theInverseRangeTable = p; 1125 // ionisation table only for ionisation p << 1143 if(1 < verboseLevel) { 1126 if (nullptr == theData->Table(i) || (!isI << 1144 G4cout << "### Set InverseRange table " << p 1127 continue; << 1145 << " for " << particle->GetParticleName() 1128 } << 1146 << " and process " << GetProcessName() << G4endl; 1129 if (-1 < verboseLevel) { << 1130 G4cout << "G4VEnergyLossProcess::StoreP << 1131 << " " << particle->GetParticleName() << 1132 << " " << GetProcessName() << 1133 << " " << tnames[i] << " " << theDat << 1134 } << 1135 if (!G4EmTableUtil::StoreTable(this, part << 1136 dir, tnames[i], verboseLevel, asci << 1137 return false; << 1138 } 1147 } 1139 } 1148 } 1140 return true; << 1141 } 1149 } 1142 1150 1143 //....oooOO0OOooo........oooOO0OOooo........o << 1151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1144 1152 1145 G4bool << 1153 void G4VEnergyLossProcess::SetLambdaTable(G4PhysicsTable* p) 1146 G4VEnergyLossProcess::RetrievePhysicsTable(co << 1147 co << 1148 { 1154 { 1149 if (!isMaster || nullptr != baseParticle || << 1155 if(1 < verboseLevel) { 1150 for(std::size_t i=0; i<7; ++i) { << 1156 G4cout << "### Set Lambda table " << p 1151 // ionisation table only for ionisation p << 1157 << " for " << particle->GetParticleName() 1152 if (!isIonisation && 1 == i) { continue; << 1158 << " and process " << GetProcessName() << G4endl; 1153 if(!G4EmTableUtil::RetrieveTable(this, pa << 1159 } 1154 verboseL << 1160 if(theLambdaTable != p) theLambdaTable = p; 1155 return false; << 1161 tablesAreBuilt = true; >> 1162 >> 1163 if(p) { >> 1164 size_t n = p->length(); >> 1165 G4PhysicsVector* pv = (*p)[0]; >> 1166 G4double e, s, smax, emax; >> 1167 theEnergyOfCrossSectionMax = new G4double [n]; >> 1168 theCrossSectionMax = new G4double [n]; >> 1169 G4bool b; >> 1170 >> 1171 for (size_t i=0; i<n; i++) { >> 1172 pv = (*p)[i]; >> 1173 emax = DBL_MAX; >> 1174 smax = 0.0; >> 1175 if(pv) { >> 1176 size_t nb = pv->GetVectorLength(); >> 1177 emax = pv->GetLowEdgeEnergy(nb); >> 1178 for (size_t j=0; j<nb; j++) { >> 1179 e = pv->GetLowEdgeEnergy(j); >> 1180 s = pv->GetValue(e,b); >> 1181 if(s > smax) { >> 1182 smax = s; >> 1183 emax = e; >> 1184 } >> 1185 } >> 1186 } >> 1187 theEnergyOfCrossSectionMax[i] = emax; >> 1188 theCrossSectionMax[i] = smax; >> 1189 if(1 < verboseLevel) { >> 1190 G4cout << "For " << particle->GetParticleName() >> 1191 << " Max CS at i= " << i << " emax(MeV)= " << emax/MeV >> 1192 << " lambda= " << smax << G4endl; >> 1193 } 1156 } 1194 } 1157 } 1195 } 1158 return true; << 1159 } 1196 } 1160 1197 1161 //....oooOO0OOooo........oooOO0OOooo........o 1198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1162 1199 1163 G4double G4VEnergyLossProcess::GetDEDXDispers << 1200 void G4VEnergyLossProcess::SetSubLambdaTable(G4PhysicsTable* p) 1164 const G4Mat << 1165 const G4Dyn << 1166 G4dou << 1167 { 1201 { 1168 DefineMaterial(couple); << 1202 if(theSubLambdaTable != p) theSubLambdaTable = p; 1169 G4double ekin = dp->GetKineticEnergy(); << 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 } 1203 } 1178 1204 1179 //....oooOO0OOooo........oooOO0OOooo........o 1205 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1180 1206 1181 G4double << 1207 G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector( 1182 G4VEnergyLossProcess::CrossSectionPerVolume(G << 1208 const G4MaterialCutsCouple* couple, G4double cut) 1183 c << 1184 G << 1185 { 1209 { 1186 // Cross section per volume is calculated << 1210 // G4double cut = (*theCuts)[couple->GetIndex()]; 1187 DefineMaterial(couple); << 1211 // G4int nbins = nLambdaBins; 1188 G4double cross = 0.0; << 1212 G4double tmin = 1189 if (nullptr != theLambdaTable) { << 1213 std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut), 1190 cross = GetLambdaForScaledEnergy(kineticE << 1214 minKinEnergy); 1191 logKinet << 1215 if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy; 1192 } else { << 1216 G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins); 1193 SelectModel(kineticEnergy*massRatio); << 1217 return v; 1194 cross = (!baseMat) ? biasFactor : biasFac << 1195 cross *= (currentModel->CrossSectionPerVo << 1196 << 1197 } << 1198 return std::max(cross, 0.0); << 1199 } 1218 } 1200 1219 1201 //....oooOO0OOooo........oooOO0OOooo........o 1220 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1202 1221 1203 G4double G4VEnergyLossProcess::MeanFreePath(c << 1222 G4double G4VEnergyLossProcess::MicroscopicCrossSection( >> 1223 G4double kineticEnergy, const G4MaterialCutsCouple* couple) 1204 { 1224 { 1205 DefineMaterial(track.GetMaterialCutsCouple( << 1225 // Cross section per atom is calculated 1206 const G4double kinEnergy = track.GetKine << 1226 DefineMaterial(couple); 1207 const G4double logKinEnergy = track.GetDyna << 1227 G4double cross = 0.0; 1208 const G4double cs = GetLambdaForScaledEnerg << 1228 G4bool b; 1209 << 1229 if(theLambdaTable) 1210 return (0.0 < cs) ? 1.0/cs : DBL_MAX; << 1230 cross = >> 1231 ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b)/ >> 1232 currentMaterial->GetTotNbOfAtomsPerVolume(); >> 1233 >> 1234 return cross; 1211 } 1235 } 1212 1236 1213 //....oooOO0OOooo........oooOO0OOooo........o 1237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1214 1238 1215 G4double G4VEnergyLossProcess::ContinuousStep << 1239 G4double G4VEnergyLossProcess::MeanFreePath( 1216 << 1240 const G4Track& track, G4double s, G4ForceCondition* cond) 1217 << 1218 { 1241 { 1219 return AlongStepGetPhysicalInteractionLengt << 1242 return GetMeanFreePath(track, s, cond); 1220 } 1243 } 1221 1244 1222 //....oooOO0OOooo........oooOO0OOooo........o 1245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1223 1246 1224 G4double G4VEnergyLossProcess::GetMeanFreePat << 1247 G4double G4VEnergyLossProcess::ContinuousStepLimit( 1225 const G4Track& t << 1248 const G4Track& track, G4double x, G4double y, G4double& z) 1226 G4double, << 1227 G4ForceCondition << 1228 << 1229 { 1249 { 1230 *condition = NotForced; << 1250 return GetContinuousStepLimit(track, x, y, z); 1231 return MeanFreePath(track); << 1232 } 1251 } 1233 1252 1234 //....oooOO0OOooo........oooOO0OOooo........o 1253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1235 1254 1236 G4double G4VEnergyLossProcess::GetContinuousS << 1255 void G4VEnergyLossProcess::SetIntegral(G4bool val) 1237 const G4Track&, << 1238 G4double, G4double, G4double& << 1239 { 1256 { 1240 return DBL_MAX; << 1257 integral = val; 1241 } 1258 } 1242 1259 1243 //....oooOO0OOooo........oooOO0OOooo........o 1260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1244 1261 1245 G4PhysicsVector* << 1262 void G4VEnergyLossProcess::SetStepFunction(G4double v1, G4double v2) 1246 G4VEnergyLossProcess::LambdaPhysicsVector(con << 1247 G4d << 1248 { 1263 { 1249 DefineMaterial(couple); << 1264 dRoverRange = v1; 1250 G4PhysicsVector* v = (*theLambdaTable)[base << 1265 finalRange = v2; 1251 return new G4PhysicsVector(*v); << 1266 if (dRoverRange > 0.999) dRoverRange = 1.0; >> 1267 currentCouple = 0; >> 1268 mfpKinEnergy = DBL_MAX; 1252 } 1269 } 1253 1270 1254 //....oooOO0OOooo........oooOO0OOooo........o 1271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1255 1272 1256 void << 1273 void G4VEnergyLossProcess::SetParticle(const G4ParticleDefinition* p) 1257 G4VEnergyLossProcess::SetDEDXTable(G4PhysicsT << 1258 { 1274 { 1259 if(1 < verboseLevel) { << 1275 particle = p; 1260 G4cout << "### Set DEDX table " << p << " << 1261 << " " << theDEDXunRestrictedTable << << 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 } 1276 } 1280 1277 1281 //....oooOO0OOooo........oooOO0OOooo........o 1278 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1282 1279 1283 void G4VEnergyLossProcess::SetCSDARangeTable( << 1280 void G4VEnergyLossProcess::SetBaseParticle(const G4ParticleDefinition* p) 1284 { 1281 { 1285 theCSDARangeTable = p; << 1282 baseParticle = p; 1286 } 1283 } 1287 1284 1288 //....oooOO0OOooo........oooOO0OOooo........o 1285 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1289 1286 1290 void G4VEnergyLossProcess::SetRangeTableForLo << 1287 void G4VEnergyLossProcess::SetSecondaryParticle(const G4ParticleDefinition* p) 1291 { 1288 { 1292 theRangeTableForLoss = p; << 1289 secondaryParticle = p; 1293 } 1290 } 1294 1291 1295 //....oooOO0OOooo........oooOO0OOooo........o 1292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1296 1293 1297 void G4VEnergyLossProcess::SetInverseRangeTab << 1294 G4bool G4VEnergyLossProcess::StorePhysicsTable( >> 1295 const G4ParticleDefinition* part, const G4String& directory, >> 1296 G4bool ascii) 1298 { 1297 { 1299 theInverseRangeTable = p; << 1298 G4bool res = true; >> 1299 if ( baseParticle || part != particle ) return res; >> 1300 >> 1301 if ( theDEDXTable ) { >> 1302 const G4String name = GetPhysicsTableFileName(part,directory,"DEDX",ascii); >> 1303 if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false; >> 1304 } >> 1305 >> 1306 if ( theDEDXunRestrictedTable ) { >> 1307 const G4String name = >> 1308 GetPhysicsTableFileName(part,directory,"DEDXnr",ascii); >> 1309 if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false; >> 1310 } >> 1311 >> 1312 if ( theDEDXSubTable ) { >> 1313 const G4String name = >> 1314 GetPhysicsTableFileName(part,directory,"SubDEDX",ascii); >> 1315 if( !theDEDXSubTable->StorePhysicsTable(name,ascii)) res = false; >> 1316 } >> 1317 >> 1318 if ( theIonisationTable ) { >> 1319 const G4String name = >> 1320 GetPhysicsTableFileName(part,directory,"Ionisation",ascii); >> 1321 if( !theIonisationTable->StorePhysicsTable(name,ascii)) res = false; >> 1322 } >> 1323 >> 1324 if ( theIonisationSubTable ) { >> 1325 const G4String name = >> 1326 GetPhysicsTableFileName(part,directory,"SubIonisation",ascii); >> 1327 if( !theIonisationSubTable->StorePhysicsTable(name,ascii)) res = false; >> 1328 } >> 1329 >> 1330 if ( theCSDARangeTable && isIonisation ) { >> 1331 const G4String name = >> 1332 GetPhysicsTableFileName(part,directory,"CSDARange",ascii); >> 1333 if( !theCSDARangeTable->StorePhysicsTable(name,ascii)) res = false; >> 1334 } >> 1335 >> 1336 if ( theRangeTableForLoss && isIonisation ) { >> 1337 const G4String name = >> 1338 GetPhysicsTableFileName(part,directory,"Range",ascii); >> 1339 if( !theRangeTableForLoss->StorePhysicsTable(name,ascii)) res = false; >> 1340 } >> 1341 >> 1342 if ( theInverseRangeTable && isIonisation ) { >> 1343 const G4String name = >> 1344 GetPhysicsTableFileName(part,directory,"InverseRange",ascii); >> 1345 if( !theInverseRangeTable->StorePhysicsTable(name,ascii)) res = false; >> 1346 } >> 1347 >> 1348 if ( theLambdaTable && isIonisation) { >> 1349 const G4String name = >> 1350 GetPhysicsTableFileName(part,directory,"Lambda",ascii); >> 1351 if( !theLambdaTable->StorePhysicsTable(name,ascii)) res = false; >> 1352 } >> 1353 >> 1354 if ( theSubLambdaTable && isIonisation) { >> 1355 const G4String name = >> 1356 GetPhysicsTableFileName(part,directory,"SubLambda",ascii); >> 1357 if( !theSubLambdaTable->StorePhysicsTable(name,ascii)) res = false; >> 1358 } >> 1359 if ( res ) { >> 1360 if(0 < verboseLevel) { >> 1361 G4cout << "Physics tables are stored for " << particle->GetParticleName() >> 1362 << " and process " << GetProcessName() >> 1363 << " in the directory <" << directory >> 1364 << "> " << G4endl; >> 1365 } >> 1366 } else { >> 1367 G4cout << "Fail to store Physics Tables for " >> 1368 << particle->GetParticleName() >> 1369 << " and process " << GetProcessName() >> 1370 << " in the directory <" << directory >> 1371 << "> " << G4endl; >> 1372 } >> 1373 return res; 1300 } 1374 } 1301 1375 1302 //....oooOO0OOooo........oooOO0OOooo........o << 1376 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 1303 1377 1304 void G4VEnergyLossProcess::SetLambdaTable(G4P << 1378 G4bool G4VEnergyLossProcess::RetrievePhysicsTable( >> 1379 const G4ParticleDefinition* part, const G4String& directory, >> 1380 G4bool ascii) 1305 { 1381 { >> 1382 G4bool res = true; >> 1383 const G4String particleName = part->GetParticleName(); >> 1384 1306 if(1 < verboseLevel) { 1385 if(1 < verboseLevel) { 1307 G4cout << "### Set Lambda table " << p << << 1386 G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for " 1308 << " for " << particle->GetParticl << 1387 << particleName << " and process " << GetProcessName() 1309 << " and process " << GetProcessNa << 1388 << "; tables_are_built= " << tablesAreBuilt >> 1389 << G4endl; 1310 } 1390 } 1311 theLambdaTable = p; << 1391 if(particle == part) { 1312 tablesAreBuilt = true; << 1313 1392 1314 if(isMaster && nullptr != p) { << 1393 G4bool yes = true; 1315 delete theEnergyOfCrossSectionMax; << 1394 G4bool fpi = true; 1316 theEnergyOfCrossSectionMax = nullptr; << 1395 if ( !baseParticle ) { 1317 if(fEmTwoPeaks == fXSType) { << 1396 G4String filename; 1318 if(nullptr != fXSpeaks) { << 1397 1319 for(auto & ptr : *fXSpeaks) { delete ptr; } << 1398 filename = GetPhysicsTableFileName(part,directory,"DEDX",ascii); 1320 delete fXSpeaks; << 1399 yes = theDEDXTable->ExistPhysicsTable(filename); >> 1400 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1401 theDEDXTable,filename,ascii); >> 1402 if(yes) { >> 1403 if (0 < verboseLevel) { >> 1404 G4cout << "DEDX table for " << particleName >> 1405 << " is Retrieved from <" >> 1406 << filename << ">" >> 1407 << G4endl; >> 1408 } >> 1409 } else { >> 1410 fpi = false; >> 1411 if (1 < verboseLevel) { >> 1412 G4cout << "DEDX table for " << particleName << " from file <" >> 1413 << filename << "> is not Retrieved" >> 1414 << G4endl; >> 1415 } 1321 } 1416 } 1322 G4LossTableBuilder* bld = lManager->Get << 1417 1323 fXSpeaks = G4EmUtility::FillPeaksStruct << 1418 filename = GetPhysicsTableFileName(part,directory,"Range",ascii); 1324 if(nullptr == fXSpeaks) { fXSType = fEm << 1419 yes = theRangeTableForLoss->ExistPhysicsTable(filename); 1325 } << 1420 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( 1326 if(fXSType == fEmOnePeak) { << 1421 theRangeTableForLoss,filename,ascii); 1327 theEnergyOfCrossSectionMax = G4EmUtilit << 1422 if(yes) { 1328 if(nullptr == theEnergyOfCrossSectionMa << 1423 if (0 < verboseLevel) { >> 1424 G4cout << "Range table for loss for " << particleName >> 1425 << " is Retrieved from <" >> 1426 << filename << ">" >> 1427 << G4endl; >> 1428 } >> 1429 } else { >> 1430 if(fpi) { >> 1431 res = false; >> 1432 G4cout << "Range table for loss for " << particleName >> 1433 << " from file <" >> 1434 << filename << "> is not Retrieved" >> 1435 << G4endl; >> 1436 } >> 1437 } >> 1438 >> 1439 filename = GetPhysicsTableFileName(part,directory,"DEDXnr",ascii); >> 1440 yes = theDEDXunRestrictedTable->ExistPhysicsTable(filename); >> 1441 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1442 theDEDXunRestrictedTable,filename,ascii); >> 1443 if(yes) { >> 1444 if (0 < verboseLevel) { >> 1445 G4cout << "Non-restricted DEDX table for " << particleName >> 1446 << " is Retrieved from <" >> 1447 << filename << ">" >> 1448 << G4endl; >> 1449 } >> 1450 } else { >> 1451 if (1 < verboseLevel) { >> 1452 G4cout << "Non-restricted DEDX table for " << particleName >> 1453 << " from file <" >> 1454 << filename << "> is not Retrieved" >> 1455 << G4endl; >> 1456 } >> 1457 } >> 1458 >> 1459 filename = GetPhysicsTableFileName(part,directory,"CSDARange",ascii); >> 1460 yes = theCSDARangeTable->ExistPhysicsTable(filename); >> 1461 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1462 theCSDARangeTable,filename,ascii); >> 1463 if(yes) { >> 1464 if (0 < verboseLevel) { >> 1465 G4cout << "CSDA Range table for " << particleName >> 1466 << " is Retrieved from <" >> 1467 << filename << ">" >> 1468 << G4endl; >> 1469 } >> 1470 } else { >> 1471 G4cout << "CSDA Range table for loss for " << particleName >> 1472 << " does not exist" >> 1473 << G4endl; >> 1474 } >> 1475 >> 1476 filename = GetPhysicsTableFileName(part,directory,"InverseRange",ascii); >> 1477 yes = theInverseRangeTable->ExistPhysicsTable(filename); >> 1478 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1479 theInverseRangeTable,filename,ascii); >> 1480 if(yes) { >> 1481 if (0 < verboseLevel) { >> 1482 G4cout << "InverseRange table for " << particleName >> 1483 << " is Retrieved from <" >> 1484 << filename << ">" >> 1485 << G4endl; >> 1486 } >> 1487 } else { >> 1488 if(fpi) { >> 1489 res = false; >> 1490 G4cout << "InverseRange table for " << particleName >> 1491 << " from file <" >> 1492 << filename << "> is not Retrieved" >> 1493 << G4endl; >> 1494 >> 1495 } >> 1496 } >> 1497 >> 1498 filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii); >> 1499 yes = theLambdaTable->ExistPhysicsTable(filename); >> 1500 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1501 theLambdaTable,filename,ascii); >> 1502 if(yes) { >> 1503 if (0 < verboseLevel) { >> 1504 G4cout << "Lambda table for " << particleName >> 1505 << " is Retrieved from <" >> 1506 << filename << ">" >> 1507 << G4endl; >> 1508 } >> 1509 } else { >> 1510 if(fpi) { >> 1511 res = false; >> 1512 G4cout << "Lambda table for " << particleName << " from file <" >> 1513 << filename << "> is not Retrieved" >> 1514 << G4endl; >> 1515 } >> 1516 } >> 1517 >> 1518 filename = GetPhysicsTableFileName(part,directory,"SubDEDX",ascii); >> 1519 yes = theDEDXSubTable->ExistPhysicsTable(filename); >> 1520 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1521 theDEDXSubTable,filename,ascii); >> 1522 if(yes) { >> 1523 if (0 < verboseLevel) { >> 1524 G4cout << "SubDEDX table for " << particleName >> 1525 << " is Retrieved from <" >> 1526 << filename << ">" >> 1527 << G4endl; >> 1528 } >> 1529 } else { >> 1530 if(nSCoffRegions) { >> 1531 res=false; >> 1532 G4cout << "SubDEDX table for " << particleName << " from file <" >> 1533 << filename << "> is not Retrieved" >> 1534 << G4endl; >> 1535 } >> 1536 } >> 1537 >> 1538 filename = GetPhysicsTableFileName(part,directory,"SubLambda",ascii); >> 1539 yes = theSubLambdaTable->ExistPhysicsTable(filename); >> 1540 if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1541 theSubLambdaTable,filename,ascii); >> 1542 if(yes) { >> 1543 if (0 < verboseLevel) { >> 1544 G4cout << "SubLambda table for " << particleName >> 1545 << " is Retrieved from <" >> 1546 << filename << ">" >> 1547 << G4endl; >> 1548 } >> 1549 } else { >> 1550 if(nSCoffRegions) { >> 1551 res=false; >> 1552 G4cout << "SubLambda table for " << particleName << " from file <" >> 1553 << filename << "> is not Retrieved" >> 1554 << G4endl; >> 1555 } >> 1556 } >> 1557 >> 1558 filename = GetPhysicsTableFileName(part,directory,"Ionisation",ascii); >> 1559 yes = theIonisationTable->ExistPhysicsTable(filename); >> 1560 if(yes) { >> 1561 theIonisationTable = >> 1562 G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable); >> 1563 >> 1564 yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1565 theIonisationTable,filename,ascii); >> 1566 } >> 1567 if(yes) { >> 1568 if (0 < verboseLevel) { >> 1569 G4cout << "Ionisation table for " << particleName >> 1570 << " is Retrieved from <" >> 1571 << filename << ">" >> 1572 << G4endl; >> 1573 } >> 1574 } >> 1575 >> 1576 filename = GetPhysicsTableFileName(part,directory,"SubIonisation",ascii); >> 1577 yes = theIonisationSubTable->ExistPhysicsTable(filename); >> 1578 if(yes) { >> 1579 theIonisationSubTable = >> 1580 G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable); >> 1581 yes = G4PhysicsTableHelper::RetrievePhysicsTable( >> 1582 theIonisationSubTable,filename,ascii); >> 1583 } >> 1584 if(yes) { >> 1585 if (0 < verboseLevel) { >> 1586 G4cout << "SubIonisation table for " << particleName >> 1587 << " is Retrieved from <" >> 1588 << filename << ">" >> 1589 << G4endl; >> 1590 } >> 1591 } 1329 } 1592 } 1330 } 1593 } >> 1594 >> 1595 return res; 1331 } 1596 } 1332 1597 1333 //....oooOO0OOooo........oooOO0OOooo........o 1598 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1334 1599 1335 void G4VEnergyLossProcess::SetEnergyOfCrossSe << 1600 void G4VEnergyLossProcess::SetLinearLossLimit(G4double val) 1336 { 1601 { 1337 theEnergyOfCrossSectionMax = p; << 1602 linLossLimit = val; 1338 } 1603 } 1339 1604 1340 //....oooOO0OOooo........oooOO0OOooo........o 1605 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1341 1606 1342 void G4VEnergyLossProcess::SetTwoPeaksXS(std: << 1607 void G4VEnergyLossProcess::SetLossFluctuations(G4bool val) 1343 { 1608 { 1344 fXSpeaks = ptr; << 1609 if(val && !lossFluctuationArePossible) return; >> 1610 lossFluctuationFlag = val; 1345 } 1611 } 1346 1612 1347 //....oooOO0OOooo........oooOO0OOooo........o 1613 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1348 1614 1349 const G4Element* G4VEnergyLossProcess::GetCur << 1615 void G4VEnergyLossProcess::SetRandomStep(G4bool val) 1350 { 1616 { 1351 return (nullptr != currentModel) << 1617 rndmStepFlag = val; 1352 ? currentModel->GetCurrentElement(current << 1353 } 1618 } 1354 1619 1355 //....oooOO0OOooo........oooOO0OOooo........o 1620 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1356 1621 1357 void G4VEnergyLossProcess::SetCrossSectionBia << 1622 void G4VEnergyLossProcess::SetMinSubRange(G4double val) 1358 << 1359 { 1623 { 1360 if(f > 0.0) { << 1624 minSubRange = val; 1361 biasFactor = f; << 1362 weightFlag = flag; << 1363 if(1 < verboseLevel) { << 1364 G4cout << "### SetCrossSectionBiasingFa << 1365 << " process " << GetProcessName << 1366 << " biasFactor= " << f << " wei << 1367 << G4endl; << 1368 } << 1369 } << 1370 } 1625 } 1371 1626 1372 //....oooOO0OOooo........oooOO0OOooo........o 1627 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1373 1628 1374 void G4VEnergyLossProcess::ActivateForcedInte << 1629 G4bool G4VEnergyLossProcess::TablesAreBuilt() const 1375 << 1376 << 1377 { 1630 { 1378 if(nullptr == biasManager) { biasManager = << 1631 return tablesAreBuilt; 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 } 1632 } 1390 1633 1391 //....oooOO0OOooo........oooOO0OOooo........o 1634 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1392 1635 1393 void << 1636 G4int G4VEnergyLossProcess::NumberOfSubCutoffRegions() const 1394 G4VEnergyLossProcess::ActivateSecondaryBiasin << 1395 << 1396 << 1397 { 1637 { 1398 if (0.0 <= factor) { << 1638 return nSCoffRegions; 1399 // Range cut can be applied only for e- << 1639 } 1400 if(0.0 == factor && secondaryParticle != << 1401 { return; } << 1402 1640 1403 if(nullptr == biasManager) { biasManager << 1641 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1404 biasManager->ActivateSecondaryBiasing(reg << 1642 1405 if(1 < verboseLevel) { << 1643 void G4VEnergyLossProcess::SetDEDXBinning(G4int nbins) 1406 G4cout << "### ActivateSecondaryBiasing << 1644 { 1407 << " process " << GetProcessName << 1645 nBins = nbins; 1408 << " factor= " << factor << 1409 << " in G4Region <" << region << 1410 << "> energyLimit(MeV)= " << ene << 1411 << G4endl; << 1412 } << 1413 } << 1414 } 1646 } 1415 1647 1416 //....oooOO0OOooo........oooOO0OOooo........o 1648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1417 1649 1418 void G4VEnergyLossProcess::SetIonisation(G4bo << 1650 void G4VEnergyLossProcess::SetLambdaBinning(G4int nbins) 1419 { 1651 { 1420 isIonisation = val; << 1652 nBins = nbins; 1421 aGPILSelection = (val) ? CandidateForSelect << 1422 } 1653 } 1423 1654 1424 //....oooOO0OOooo........oooOO0OOooo........o 1655 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1425 1656 1426 void G4VEnergyLossProcess::SetLinearLossLimi << 1657 void G4VEnergyLossProcess::SetDEDXBinningForCSDARange(G4int nbins) 1427 { 1658 { 1428 if(0.0 < val && val < 1.0) { << 1659 nBinsCSDA = nbins; 1429 linLossLimit = val; << 1430 actLinLossLimit = true; << 1431 } else { PrintWarning("SetLinearLossLimit", << 1432 } 1660 } 1433 1661 1434 //....oooOO0OOooo........oooOO0OOooo........o 1662 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1435 1663 1436 void G4VEnergyLossProcess::SetStepFunction(G4 << 1664 G4double G4VEnergyLossProcess::MinKinEnergy() const 1437 { 1665 { 1438 if(0.0 < v1 && 0.0 < v2) { << 1666 return minKinEnergy; 1439 dRoverRange = std::min(1.0, v1); << 1440 finalRange = std::min(v2, 1.e+50); << 1441 } else { << 1442 PrintWarning("SetStepFunctionV1", v1); << 1443 PrintWarning("SetStepFunctionV2", v2); << 1444 } << 1445 } 1667 } 1446 1668 1447 //....oooOO0OOooo........oooOO0OOooo........o 1669 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1448 1670 1449 void G4VEnergyLossProcess::SetLowestEnergyLim << 1671 void G4VEnergyLossProcess::SetMinKinEnergy(G4double e) 1450 { 1672 { 1451 if(1.e-18 < val && val < 1.e+50) { lowestKi << 1673 minKinEnergy = e; 1452 else { PrintWarning("SetLowestEnergyLimit", << 1453 } 1674 } 1454 1675 1455 //....oooOO0OOooo........oooOO0OOooo........o 1676 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1456 1677 1457 void G4VEnergyLossProcess::SetDEDXBinning(G4i << 1678 void G4VEnergyLossProcess::SetMaxKinEnergy(G4double e) 1458 { 1679 { 1459 if(2 < n && n < 1000000000) { << 1680 maxKinEnergy = e; 1460 nBins = n; << 1681 if(e < maxKinEnergyCSDA) maxKinEnergyCSDA = e; 1461 actBinning = true; << 1462 } else { << 1463 G4double e = (G4double)n; << 1464 PrintWarning("SetDEDXBinning", e); << 1465 } << 1466 } 1682 } 1467 1683 1468 //....oooOO0OOooo........oooOO0OOooo........o 1684 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1469 1685 1470 void G4VEnergyLossProcess::SetMinKinEnergy(G4 << 1686 void G4VEnergyLossProcess::SetMaxKinEnergyForCSDARange(G4double e) 1471 { 1687 { 1472 if(1.e-18 < e && e < maxKinEnergy) { << 1688 maxKinEnergyCSDA = e; 1473 minKinEnergy = e; << 1474 actMinKinEnergy = true; << 1475 } else { PrintWarning("SetMinKinEnergy", e) << 1476 } 1689 } 1477 1690 1478 //....oooOO0OOooo........oooOO0OOooo........o 1691 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1479 1692 1480 void G4VEnergyLossProcess::SetMaxKinEnergy(G4 << 1693 G4double G4VEnergyLossProcess::MaxKinEnergy() const >> 1694 { >> 1695 return maxKinEnergy; >> 1696 } >> 1697 >> 1698 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1699 >> 1700 void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*) >> 1701 {} >> 1702 >> 1703 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1704 >> 1705 void G4VEnergyLossProcess::SetLambdaFactor(G4double val) 1481 { 1706 { 1482 if(minKinEnergy < e && e < 1.e+50) { << 1707 if(val > 0.0 && val <= 1.0) lambdaFactor = val; 1483 maxKinEnergy = e; << 1484 actMaxKinEnergy = true; << 1485 if(e < maxKinEnergyCSDA) { maxKinEnergyCS << 1486 } else { PrintWarning("SetMaxKinEnergy", e) << 1487 } 1708 } 1488 1709 1489 //....oooOO0OOooo........oooOO0OOooo........o 1710 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1490 1711 1491 void G4VEnergyLossProcess::PrintWarning(const << 1712 void G4VEnergyLossProcess::SetIonisation(G4bool val) 1492 { 1713 { 1493 G4String ss = "G4VEnergyLossProcess::" + ti << 1714 isIonisation = val; 1494 G4ExceptionDescription ed; << 1495 ed << "Parameter is out of range: " << val << 1496 << " it will have no effect!\n" << " Pr << 1497 << GetProcessName() << " nbins= " << nB << 1498 << " Emin(keV)= " << minKinEnergy/keV << 1499 << " Emax(GeV)= " << maxKinEnergy/GeV; << 1500 G4Exception(ss, "em0044", JustWarning, ed); << 1501 } 1715 } 1502 1716 1503 //....oooOO0OOooo........oooOO0OOooo........o 1717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 1504 1718 1505 void G4VEnergyLossProcess::ProcessDescription << 1719 G4bool G4VEnergyLossProcess::IsIonisationProcess() const 1506 { 1720 { 1507 if(nullptr != particle) { StreamInfo(out, * << 1721 return isIonisation; 1508 } 1722 } 1509 1723 1510 //....oooOO0OOooo........oooOO0OOooo........o 1724 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1725 >> 1726 void G4VEnergyLossProcess::AddCollaborativeProcess( >> 1727 G4VEnergyLossProcess* p) >> 1728 { >> 1729 scProcesses.push_back(p); >> 1730 nProcesses++; >> 1731 if (0 < verboseLevel) >> 1732 G4cout << "### The process " << p->GetProcessName() >> 1733 << " is added to the list of collaborative processes of " >> 1734 << GetProcessName() << G4endl; >> 1735 } >> 1736 >> 1737 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 1738 1511 1739