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: G4HadronicProcess.cc 93817 2015-11-02 11:33:26Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class source file 30 // GEANT4 Class source file 30 // 31 // 31 // G4HadronicProcess 32 // G4HadronicProcess 32 // 33 // 33 // original by H.P.Wellisch 34 // original by H.P.Wellisch 34 // J.L. Chuma, TRIUMF, 10-Mar-1997 35 // J.L. Chuma, TRIUMF, 10-Mar-1997 35 // 36 // 36 // Modifications: 37 // Modifications: 37 // 05-Jul-2010 V.Ivanchenko cleanup commented 38 // 05-Jul-2010 V.Ivanchenko cleanup commented lines 38 // 20-Jul-2011 M.Kelsey -- null-pointer checks 39 // 20-Jul-2011 M.Kelsey -- null-pointer checks in DumpState() 39 // 24-Sep-2011 M.Kelsey -- Use envvar G4HADRON 40 // 24-Sep-2011 M.Kelsey -- Use envvar G4HADRONIC_RANDOM_FILE to save random 40 // engine state before each model call 41 // engine state before each model call 41 // 18-Oct-2011 M.Kelsey -- Handle final-state 42 // 18-Oct-2011 M.Kelsey -- Handle final-state cases in conservation checks. 42 // 14-Mar-2012 G.Folger -- enhance checks for 43 // 14-Mar-2012 G.Folger -- enhance checks for conservation of energy, etc. 43 // 28-Jul-2012 M.Maire -- add function GetTar 44 // 28-Jul-2012 M.Maire -- add function GetTargetDefinition() 44 // 14-Sep-2012 Inherit from RestDiscrete, use 45 // 14-Sep-2012 Inherit from RestDiscrete, use subtype code (now in ctor) to 45 // configure base-class 46 // configure base-class 46 // 28-Sep-2012 Restore inheritance from G4VDis 47 // 28-Sep-2012 Restore inheritance from G4VDiscreteProcess, remove enable-flag 47 // changing, remove warning message from or 48 // changing, remove warning message from original ctor. 48 // 21-Aug-2019 V.Ivanchenko leave try/catch on << 49 49 50 #include "G4HadronicProcess.hh" 50 #include "G4HadronicProcess.hh" 51 51 52 #include "G4Types.hh" 52 #include "G4Types.hh" 53 #include "G4SystemOfUnits.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4HadProjectile.hh" 54 #include "G4HadProjectile.hh" 55 #include "G4ElementVector.hh" 55 #include "G4ElementVector.hh" 56 #include "G4Track.hh" 56 #include "G4Track.hh" 57 #include "G4Step.hh" 57 #include "G4Step.hh" 58 #include "G4Element.hh" 58 #include "G4Element.hh" 59 #include "G4ParticleChange.hh" 59 #include "G4ParticleChange.hh" >> 60 #include "G4TransportationManager.hh" >> 61 #include "G4Navigator.hh" 60 #include "G4ProcessVector.hh" 62 #include "G4ProcessVector.hh" 61 #include "G4ProcessManager.hh" 63 #include "G4ProcessManager.hh" >> 64 #include "G4StableIsotopes.hh" >> 65 #include "G4HadTmpUtil.hh" 62 #include "G4NucleiProperties.hh" 66 #include "G4NucleiProperties.hh" 63 67 64 #include "G4HadronicException.hh" 68 #include "G4HadronicException.hh" 65 #include "G4HadronicProcessStore.hh" 69 #include "G4HadronicProcessStore.hh" 66 #include "G4HadronicParameters.hh" << 67 #include "G4VCrossSectionDataSet.hh" << 68 70 >> 71 #include "G4AutoLock.hh" 69 #include "G4NistManager.hh" 72 #include "G4NistManager.hh" 70 #include "G4VLeadingParticleBiasing.hh" << 71 #include "G4HadXSHelper.hh" << 72 #include "G4Threading.hh" << 73 #include "G4Exp.hh" << 74 73 75 #include <typeinfo> 74 #include <typeinfo> 76 #include <sstream> 75 #include <sstream> 77 #include <iostream> 76 #include <iostream> 78 77 79 namespace << 78 #include <stdlib.h> 80 { << 79 81 constexpr G4double lambdaFactor = 0.8; << 80 // File-scope variable to capture environment variable at startup 82 constexpr G4double invLambdaFactor = 1.0/lam << 81 83 } << 82 static const char* G4Hadronic_Random_File = getenv("G4HADRONIC_RANDOM_FILE"); 84 83 85 ////////////////////////////////////////////// 84 ////////////////////////////////////////////////////////////////// 86 85 87 G4HadronicProcess::G4HadronicProcess(const G4S 86 G4HadronicProcess::G4HadronicProcess(const G4String& processName, 88 G4Process 87 G4ProcessType procType) 89 : G4VDiscreteProcess(processName, procType) 88 : G4VDiscreteProcess(processName, procType) 90 { 89 { 91 SetProcessSubType(fHadronInelastic); // Def 90 SetProcessSubType(fHadronInelastic); // Default unless subclass changes 92 InitialiseLocal(); << 91 >> 92 theTotalResult = new G4ParticleChange(); >> 93 theTotalResult->SetSecondaryWeightByProcess(true); >> 94 theInteraction = 0; >> 95 theCrossSectionDataStore = new G4CrossSectionDataStore(); >> 96 theProcessStore = G4HadronicProcessStore::Instance(); >> 97 theProcessStore->Register(this); >> 98 theInitialNumberOfInteractionLength = 0.0; >> 99 aScaleFactor = 1; >> 100 xBiasOn = false; >> 101 theLastCrossSection = 0.0; >> 102 G4HadronicProcess_debug_flag = false; >> 103 GetEnergyMomentumCheckEnvvars(); 93 } 104 } 94 105 >> 106 ////////////////////////////////////////////////////////////////// >> 107 95 G4HadronicProcess::G4HadronicProcess(const G4S 108 G4HadronicProcess::G4HadronicProcess(const G4String& processName, 96 G4Hadroni 109 G4HadronicProcessType aHadSubType) 97 : G4VDiscreteProcess(processName, fHadronic) 110 : G4VDiscreteProcess(processName, fHadronic) 98 { 111 { 99 SetProcessSubType(aHadSubType); 112 SetProcessSubType(aHadSubType); 100 InitialiseLocal(); << 113 >> 114 theTotalResult = new G4ParticleChange(); >> 115 theTotalResult->SetSecondaryWeightByProcess(true); >> 116 theInteraction = 0; >> 117 theCrossSectionDataStore = new G4CrossSectionDataStore(); >> 118 theProcessStore = G4HadronicProcessStore::Instance(); >> 119 theProcessStore->Register(this); >> 120 theInitialNumberOfInteractionLength = 0.0; >> 121 aScaleFactor = 1; >> 122 xBiasOn = false; >> 123 theLastCrossSection = 0.0; >> 124 G4HadronicProcess_debug_flag = false; >> 125 GetEnergyMomentumCheckEnvvars(); 101 } 126 } 102 127 >> 128 103 G4HadronicProcess::~G4HadronicProcess() 129 G4HadronicProcess::~G4HadronicProcess() 104 { 130 { 105 theProcessStore->DeRegister(this); 131 theProcessStore->DeRegister(this); 106 delete theTotalResult; 132 delete theTotalResult; 107 delete theCrossSectionDataStore; 133 delete theCrossSectionDataStore; 108 if(isMaster) { << 109 if (fXSpeaks != nullptr) { << 110 for (auto const& e : *fXSpeaks ) { << 111 delete e; << 112 } << 113 } << 114 delete fXSpeaks; << 115 delete theEnergyOfCrossSectionMax; << 116 } << 117 } 134 } 118 135 119 void G4HadronicProcess::InitialiseLocal() { << 136 void G4HadronicProcess::GetEnergyMomentumCheckEnvvars() { 120 theTotalResult = new G4ParticleChange(); << 137 levelsSetByProcess = false; 121 theTotalResult->SetSecondaryWeightByProcess( << 138 122 theCrossSectionDataStore = new G4CrossSectio << 139 epReportLevel = getenv("G4Hadronic_epReportLevel") ? 123 theProcessStore = G4HadronicProcessStore::In << 140 strtol(getenv("G4Hadronic_epReportLevel"),0,10) : 0; 124 theProcessStore->Register(this); << 125 minKinEnergy = 1*CLHEP::MeV; << 126 141 127 G4HadronicParameters* param = G4HadronicPara << 142 epCheckLevels.first = getenv("G4Hadronic_epCheckRelativeLevel") ? 128 epReportLevel = param->GetEPReportLevel(); << 143 strtod(getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX; 129 epCheckLevels.first = param->GetEPRelativeLe << 130 epCheckLevels.second = param->GetEPAbsoluteL << 131 144 132 unitVector.set(0.0, 0.0, 0.1); << 145 epCheckLevels.second = getenv("G4Hadronic_epCheckAbsoluteLevel") ? 133 if(G4Threading::IsWorkerThread()) { isMaster << 146 strtod(getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX; 134 } 147 } 135 148 136 void G4HadronicProcess::RegisterMe( G4Hadronic 149 void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a ) 137 { 150 { 138 if(nullptr == a) { return; } << 151 if(!a) { return; } 139 theEnergyRangeManager.RegisterMe( a ); << 152 try{ theEnergyRangeManager.RegisterMe( a ); } >> 153 catch(G4HadronicException & aE) >> 154 { >> 155 G4ExceptionDescription ed; >> 156 aE.Report(ed); >> 157 ed << "Unrecoverable error in " << GetProcessName() >> 158 << " to register " << a->GetModelName() << G4endl; >> 159 G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException, >> 160 ed); >> 161 } 140 G4HadronicProcessStore::Instance()->Register 162 G4HadronicProcessStore::Instance()->RegisterInteraction(this, a); 141 } 163 } 142 164 143 G4double << 165 G4double G4HadronicProcess::GetElementCrossSection(const G4DynamicParticle * part, 144 G4HadronicProcess::GetElementCrossSection(cons << 166 const G4Element * elm, 145 const G4Element * elm, << 167 const G4Material* mat) 146 const G4Material* mat) << 168 { 147 { << 169 G4Material* aMaterial = const_cast<G4Material*>(mat); 148 if(nullptr == mat) << 170 if(! mat) 149 { << 171 { 150 static const G4int nmax = 5; << 172 // Because NeutronHP needs a material pointer (for instance to get the 151 if(nMatWarn < nmax) { << 173 // temperature), we ask the Nist manager to build or find a simple material 152 ++nMatWarn; << 174 // from the (integer) Z of the element. 153 G4ExceptionDescription ed; << 175 // Note that repeated calls to this method are not producing multiple copies 154 ed << "Cannot compute Element x-section << 176 // of the same material. But it needs to be protected against race conditions 155 << " because no material defined \n" << 177 // between different threads. 156 << " Please, specify material pointer or de << 178 aMaterial = InitialiseMaterial(G4int(elm->GetZ())); 157 << " for Z= " << elm->GetZasInt(); << 158 G4Exception("G4HadronicProcess::GetEleme << 159 JustWarning, ed); << 160 } << 161 } 179 } 162 return theCrossSectionDataStore->GetCrossSec << 180 G4double x = theCrossSectionDataStore->GetCrossSection(part, elm, aMaterial); >> 181 if(x < 0.0) { x = 0.0; } >> 182 return x; >> 183 } >> 184 >> 185 >> 186 namespace { G4Mutex hadronicProcessMutex = G4MUTEX_INITIALIZER; } >> 187 G4Material* G4HadronicProcess::InitialiseMaterial(G4int Z) >> 188 { >> 189 G4AutoLock l(&hadronicProcessMutex); >> 190 return G4NistManager::Instance()->FindOrBuildSimpleMaterial(Z); 163 } 191 } 164 192 >> 193 165 void G4HadronicProcess::PreparePhysicsTable(co 194 void G4HadronicProcess::PreparePhysicsTable(const G4ParticleDefinition& p) 166 { 195 { 167 if(nullptr == firstParticle) { firstParticle << 196 if(getenv("G4HadronicProcess_debug")) { >> 197 G4HadronicProcess_debug_flag = true; >> 198 } 168 theProcessStore->RegisterParticle(this, &p); 199 theProcessStore->RegisterParticle(this, &p); 169 } 200 } 170 201 171 void G4HadronicProcess::BuildPhysicsTable(cons 202 void G4HadronicProcess::BuildPhysicsTable(const G4ParticleDefinition& p) 172 { 203 { 173 if(firstParticle != &p) { return; } << 204 try 174 << 205 { 175 theCrossSectionDataStore->BuildPhysicsTable( << 206 theCrossSectionDataStore->BuildPhysicsTable(p); 176 theEnergyRangeManager.BuildPhysicsTable(p); << 207 theEnergyRangeManager.BuildPhysicsTable(p); 177 G4HadronicParameters* param = G4HadronicPara << 178 << 179 G4int subtype = GetProcessSubType(); << 180 if(useIntegralXS) { << 181 if(subtype == fHadronInelastic) { << 182 useIntegralXS = param->EnableIntegralIne << 183 } else if(subtype == fHadronElastic) { << 184 useIntegralXS = param->EnableIntegralEla << 185 } << 186 } << 187 fXSType = fHadNoIntegral; << 188 << 189 if(nullptr == masterProcess) { << 190 masterProcess = dynamic_cast<const G4Hadro << 191 } << 192 if(nullptr == masterProcess) { << 193 if(1 < param->GetVerboseLevel()) { << 194 G4ExceptionDescription ed; << 195 ed << "G4HadronicProcess::BuildPhysicsTa << 196 << GetProcessName() << " for " << p.GetPart << 197 << " fail due to undefined pointer to the m << 198 << " ThreadID= " << G4Threading::G4GetThre << 199 << " initialisation of worker started befo << 200 G4Exception("G4HadronicProcess::BuildPhy << 201 JustWarning, ed); << 202 } << 203 } 208 } 204 << 209 catch(G4HadronicException aR) 205 // check particle for integral method << 210 { 206 if(isMaster || nullptr == masterProcess) { << 211 G4ExceptionDescription ed; 207 G4double charge = p.GetPDGCharge()/eplus; << 212 aR.Report(ed); 208 << 213 ed << " hadronic initialisation fails" << G4endl; 209 // select cross section shape << 214 G4Exception("G4HadronicProcess::BuildPhysicsTable", "had000", 210 if(charge != 0.0 && useIntegralXS) { << 215 FatalException,ed); 211 G4double tmax = param->GetMaxEnergy(); << 212 currentParticle = firstParticle; << 213 // initialisation in the master thread << 214 G4int pdg = p.GetPDGEncoding(); << 215 if (std::abs(pdg) == 211) { << 216 fXSType = fHadTwoPeaks; << 217 } else if (pdg == 321) { << 218 fXSType = fHadOnePeak; << 219 } else if (pdg == -321) { << 220 fXSType = fHadDecreasing; << 221 } else if (pdg == 2212) { << 222 fXSType = fHadTwoPeaks; << 223 } else if (pdg == -2212 || pdg == -10000 << 224 pdg == -1000020030 || pdg == -1000020040) << 225 fXSType = fHadDecreasing; << 226 } else if (charge > 0.0 || pdg == 11 || << 227 fXSType = fHadIncreasing; << 228 } << 229 << 230 delete theEnergyOfCrossSectionMax; << 231 theEnergyOfCrossSectionMax = nullptr; << 232 if(fXSType == fHadTwoPeaks) { << 233 if (fXSpeaks != nullptr) { << 234 for (auto const& e : *fXSpeaks ) { << 235 delete e; << 236 } << 237 } << 238 delete fXSpeaks; << 239 fXSpeaks = << 240 G4HadXSHelper::FillPeaksStructure(this, &p << 241 if(nullptr == fXSpeaks) { << 242 fXSType = fHadOnePeak; << 243 } << 244 } << 245 if(fXSType == fHadOnePeak) { << 246 theEnergyOfCrossSectionMax = << 247 G4HadXSHelper::FindCrossSectionMax(this, & << 248 if(nullptr == theEnergyOfCrossSectionMax) { << 249 fXSType = fHadIncreasing; << 250 } << 251 } << 252 } << 253 } else { << 254 // initialisation in worker threads << 255 fXSType = masterProcess->CrossSectionType( << 256 fXSpeaks = masterProcess->TwoPeaksXS(); << 257 theEnergyOfCrossSectionMax = masterProcess << 258 } << 259 if(isMaster && 1 < param->GetVerboseLevel()) << 260 G4cout << "G4HadronicProcess::BuildPhysics << 261 << GetProcessName() << " and " << p.GetPa << 262 << " typeXS=" << fXSType << G4endl; << 263 } 216 } 264 G4HadronicProcessStore::Instance()->PrintInf 217 G4HadronicProcessStore::Instance()->PrintInfo(&p); 265 } 218 } 266 219 267 void G4HadronicProcess::StartTracking(G4Track* << 220 G4double G4HadronicProcess:: 268 { << 221 GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) 269 currentMat = nullptr; << 270 currentParticle = track->GetDefinition(); << 271 fDynParticle = track->GetDynamicParticle(); << 272 theNumberOfInteractionLengthLeft = -1.0; << 273 } << 274 << 275 G4double G4HadronicProcess::PostStepGetPhysica << 276 const G4Track& tr << 277 G4double previousStepSize, << 278 G4ForceCondition* << 279 { 222 { 280 *condition = NotForced; << 223 //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName() 281 << 224 // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl; 282 const G4Material* mat = track.GetMaterial(); << 225 try 283 if(mat != currentMat) { << 226 { 284 currentMat = mat; << 227 theLastCrossSection = aScaleFactor* 285 mfpKinEnergy = DBL_MAX; << 228 theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(), 286 matIdx = (G4int)track.GetMaterial()->GetIn << 229 aTrack.GetMaterial()); 287 } << 288 UpdateCrossSectionAndMFP(track.GetKineticEne << 289 << 290 // zero cross section << 291 if(theLastCrossSection <= 0.0) { << 292 theNumberOfInteractionLengthLeft = -1.0; << 293 currentInteractionLength = DBL_MAX; << 294 return DBL_MAX; << 295 } << 296 << 297 // non-zero cross section << 298 if (theNumberOfInteractionLengthLeft < 0.0) << 299 theNumberOfInteractionLengthLeft = -G4Log( << 300 theInitialNumberOfInteractionLength = theN << 301 } else { << 302 theNumberOfInteractionLengthLeft -= << 303 previousStepSize/currentInteractionLengt << 304 theNumberOfInteractionLengthLeft = << 305 std::max(theNumberOfInteractionLengthLef << 306 } 230 } 307 currentInteractionLength = theMFP; << 231 catch(G4HadronicException aR) 308 return theNumberOfInteractionLengthLeft*theM << 232 { 309 } << 233 G4ExceptionDescription ed; 310 << 234 aR.Report(ed); 311 G4double G4HadronicProcess::GetMeanFreePath( << 235 DumpState(aTrack,"GetMeanFreePath",ed); 312 const G4Track &aTr << 236 ed << " Cross section is not available" << G4endl; 313 G4ForceCondition*) << 237 G4Exception("G4HadronicProcess::GetMeanFreePath", "had002", FatalException, 314 { << 238 ed); 315 G4double xs = aScaleFactor*theCrossSectionDa << 239 } 316 ->ComputeCrossSection(aTrack.GetDynamicPa << 240 G4double res = DBL_MAX; 317 return (xs > 0.0) ? 1.0/xs : DBL_MAX; << 241 if( theLastCrossSection > 0.0 ) { res = 1.0/theLastCrossSection; } >> 242 //G4cout << " xsection= " << res << G4endl; >> 243 return res; 318 } 244 } 319 245 320 G4VParticleChange* 246 G4VParticleChange* 321 G4HadronicProcess::PostStepDoIt(const G4Track& 247 G4HadronicProcess::PostStepDoIt(const G4Track& aTrack, const G4Step&) 322 { 248 { 323 theNumberOfInteractionLengthLeft = -1.0; << 324 << 325 //G4cout << "PostStepDoIt " << aTrack.GetDef 249 //G4cout << "PostStepDoIt " << aTrack.GetDefinition()->GetParticleName() 326 // << " Ekin= " << aTrack.GetKineticEnergy 250 // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl; 327 // if primary is not Alive then do nothing 251 // if primary is not Alive then do nothing 328 theTotalResult->Clear(); 252 theTotalResult->Clear(); 329 theTotalResult->Initialize(aTrack); 253 theTotalResult->Initialize(aTrack); 330 fWeight = aTrack.GetWeight(); << 254 theTotalResult->ProposeWeight(aTrack.GetWeight()); 331 theTotalResult->ProposeWeight(fWeight); << 332 if(aTrack.GetTrackStatus() != fAlive) { retu 255 if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; } 333 256 334 // Find cross section at end of step and che 257 // Find cross section at end of step and check if <= 0 335 // 258 // 336 const G4DynamicParticle* aParticle = aTrack. 259 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 337 const G4Material* aMaterial = aTrack.GetMate << 260 G4Material* aMaterial = aTrack.GetMaterial(); >> 261 >> 262 G4Element* anElement = 0; >> 263 try >> 264 { >> 265 anElement = theCrossSectionDataStore->SampleZandA(aParticle, >> 266 aMaterial, >> 267 targetNucleus); >> 268 } >> 269 catch(G4HadronicException & aR) >> 270 { >> 271 G4ExceptionDescription ed; >> 272 aR.Report(ed); >> 273 DumpState(aTrack,"SampleZandA",ed); >> 274 ed << " PostStepDoIt failed on element selection" << G4endl; >> 275 G4Exception("G4HadronicProcess::PostStepDoIt", "had003", FatalException, >> 276 ed); >> 277 } 338 278 339 // check only for charged particles 279 // check only for charged particles 340 if(fXSType != fHadNoIntegral) { << 280 if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) { 341 mfpKinEnergy = DBL_MAX; << 281 if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) { 342 G4double xs = aScaleFactor* << 343 theCrossSectionDataStore->ComputeCrossSe << 344 //G4cout << "xs=" << xs << " xs0=" << theL << 345 // << " " << aMaterial->GetName() << << 346 if(xs < theLastCrossSection*G4UniformRand( << 347 // No interaction 282 // No interaction 348 return theTotalResult; 283 return theTotalResult; 349 } 284 } 350 } 285 } 351 286 352 const G4Element* anElement = << 353 theCrossSectionDataStore->SampleZandA(aPar << 354 << 355 // Next check for illegal track status 287 // Next check for illegal track status 356 // 288 // 357 if (aTrack.GetTrackStatus() != fAlive && 289 if (aTrack.GetTrackStatus() != fAlive && 358 aTrack.GetTrackStatus() != fSuspend) { 290 aTrack.GetTrackStatus() != fSuspend) { 359 if (aTrack.GetTrackStatus() == fStopAndKil 291 if (aTrack.GetTrackStatus() == fStopAndKill || 360 aTrack.GetTrackStatus() == fKillTrackA 292 aTrack.GetTrackStatus() == fKillTrackAndSecondaries || 361 aTrack.GetTrackStatus() == fPostponeTo 293 aTrack.GetTrackStatus() == fPostponeToNextEvent) { 362 G4ExceptionDescription ed; 294 G4ExceptionDescription ed; 363 ed << "G4HadronicProcess: track in unusa 295 ed << "G4HadronicProcess: track in unusable state - " 364 << aTrack.GetTrackStatus() << G4endl; 296 << aTrack.GetTrackStatus() << G4endl; 365 ed << "G4HadronicProcess: returning unch 297 ed << "G4HadronicProcess: returning unchanged track " << G4endl; 366 DumpState(aTrack,"PostStepDoIt",ed); 298 DumpState(aTrack,"PostStepDoIt",ed); 367 G4Exception("G4HadronicProcess::PostStep 299 G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed); 368 } 300 } 369 // No warning for fStopButAlive which is a 301 // No warning for fStopButAlive which is a legal status here 370 return theTotalResult; 302 return theTotalResult; 371 } 303 } 372 304 373 // Initialize the hadronic projectile from t 305 // Initialize the hadronic projectile from the track 374 thePro.Initialise(aTrack); 306 thePro.Initialise(aTrack); 375 307 376 theInteraction = ChooseHadronicInteraction(t << 308 try 377 a << 309 { 378 if(nullptr == theInteraction) { << 310 theInteraction = >> 311 ChooseHadronicInteraction( thePro, targetNucleus, aMaterial, anElement ); >> 312 } >> 313 catch(G4HadronicException & aE) >> 314 { 379 G4ExceptionDescription ed; 315 G4ExceptionDescription ed; >> 316 aE.Report(ed); 380 ed << "Target element "<<anElement->GetNam 317 ed << "Target element "<<anElement->GetName()<<" Z= " 381 << targetNucleus.GetZ_asInt() << " A= 318 << targetNucleus.GetZ_asInt() << " A= " 382 << targetNucleus.GetA_asInt() << G4endl 319 << targetNucleus.GetA_asInt() << G4endl; 383 DumpState(aTrack,"ChooseHadronicInteractio 320 DumpState(aTrack,"ChooseHadronicInteraction",ed); 384 ed << " No HadronicInteraction found out" 321 ed << " No HadronicInteraction found out" << G4endl; 385 G4Exception("G4HadronicProcess::PostStepDo << 322 G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException, 386 FatalException, ed); << 323 ed); 387 return theTotalResult; << 388 } 324 } 389 325 390 G4HadFinalState* result = nullptr; << 326 G4HadFinalState* result = 0; 391 G4int reentryCount = 0; 327 G4int reentryCount = 0; 392 /* << 328 393 G4cout << "### " << aParticle->GetDefinition << 394 << " Ekin(MeV)= " << aParticle->GetKinetic << 395 << " Z= " << targetNucleus.GetZ_asInt() << 396 << " A= " << targetNucleus.GetA_asInt() << 397 << " by " << theInteraction->GetModelName( << 398 << G4endl; << 399 */ << 400 do 329 do 401 { 330 { 402 try 331 try 403 { 332 { >> 333 // Save random engine if requested for debugging >> 334 if (G4Hadronic_Random_File) { >> 335 CLHEP::HepRandom::saveEngineStatus(G4Hadronic_Random_File); >> 336 } 404 // Call the interaction 337 // Call the interaction 405 result = theInteraction->ApplyYourself( 338 result = theInteraction->ApplyYourself( thePro, targetNucleus); 406 ++reentryCount; 339 ++reentryCount; 407 } 340 } 408 catch(G4HadronicException & aR) << 341 catch(G4HadronicException aR) 409 { 342 { 410 G4ExceptionDescription ed; 343 G4ExceptionDescription ed; 411 aR.Report(ed); 344 aR.Report(ed); 412 ed << "Call for " << theInteraction->Get 345 ed << "Call for " << theInteraction->GetModelName() << G4endl; 413 ed << "Target element "<<anElement->GetN 346 ed << "Target element "<<anElement->GetName()<<" Z= " 414 << targetNucleus.GetZ_asInt() 347 << targetNucleus.GetZ_asInt() 415 << " A= " << targetNucleus.GetA_asInt() << 348 << " A= " << targetNucleus.GetA_asInt() << G4endl; 416 DumpState(aTrack,"ApplyYourself",ed); 349 DumpState(aTrack,"ApplyYourself",ed); 417 ed << " ApplyYourself failed" << G4endl; 350 ed << " ApplyYourself failed" << G4endl; 418 G4Exception("G4HadronicProcess::PostStep 351 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException, 419 ed); 352 ed); 420 } 353 } 421 354 422 // Check the result for catastrophic energ 355 // Check the result for catastrophic energy non-conservation 423 result = CheckResult(thePro, targetNucleus << 356 CheckResult(thePro, targetNucleus, result); 424 357 425 if(reentryCount>100) { 358 if(reentryCount>100) { 426 G4ExceptionDescription ed; 359 G4ExceptionDescription ed; 427 ed << "Call for " << theInteraction->Get 360 ed << "Call for " << theInteraction->GetModelName() << G4endl; 428 ed << "Target element "<<anElement->GetN 361 ed << "Target element "<<anElement->GetName()<<" Z= " 429 << targetNucleus.GetZ_asInt() 362 << targetNucleus.GetZ_asInt() 430 << " A= " << targetNucleus.GetA_asInt() << 363 << " A= " << targetNucleus.GetA_asInt() << G4endl; 431 DumpState(aTrack,"ApplyYourself",ed); 364 DumpState(aTrack,"ApplyYourself",ed); 432 ed << " ApplyYourself does not completed 365 ed << " ApplyYourself does not completed after 100 attempts" << G4endl; 433 G4Exception("G4HadronicProcess::PostStep 366 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException, 434 ed); 367 ed); 435 } 368 } 436 } 369 } 437 while(!result); /* Loop checking, 30-Oct-20 370 while(!result); /* Loop checking, 30-Oct-2015, G.Folger */ 438 371 439 // Check whether kaon0 or anti_kaon0 are pre << 440 // if this is the case, transform them into << 441 // with equal, 50% probability, keeping thei << 442 // the other kinematical properties). << 443 // When this happens - very rarely - a "Just << 444 // Because Fluka-Cern produces kaon0 and ant << 445 // of warnings to max 1 per thread. << 446 G4int nSec = (G4int)result->GetNumberOfSecon << 447 if ( nSec > 0 ) { << 448 for ( G4int i = 0; i < nSec; ++i ) { << 449 auto dynamicParticle = result->GetSecond << 450 auto part = dynamicParticle->GetParticle << 451 if ( part == G4KaonZero::Definition() || << 452 part == G4AntiKaonZero::Definition( << 453 G4ParticleDefinition* newPart; << 454 if ( G4UniformRand() > 0.5 ) { newPart << 455 else { newPart = G4KaonZeroLong::Defin << 456 dynamicParticle->SetDefinition( newPar << 457 if ( nKaonWarn < 1 ) { << 458 ++nKaonWarn; << 459 G4ExceptionDescription ed; << 460 ed << " Hadronic model " << theInteraction << 461 ed << " created " << part->GetParticleName << 462 ed << " -> forced to be " << newPart->GetP << 463 G4Exception( "G4HadronicProcess::PostStepD << 464 } << 465 } << 466 } << 467 } << 468 << 469 result->SetTrafoToLab(thePro.GetTrafoToLab() 372 result->SetTrafoToLab(thePro.GetTrafoToLab()); >> 373 >> 374 ClearNumberOfInteractionLengthLeft(); >> 375 470 FillResult(result, aTrack); 376 FillResult(result, aTrack); 471 377 472 if (epReportLevel != 0) { 378 if (epReportLevel != 0) { 473 CheckEnergyMomentumConservation(aTrack, ta 379 CheckEnergyMomentumConservation(aTrack, targetNucleus); 474 } 380 } 475 //G4cout << "PostStepDoIt done nICelectrons= << 381 //G4cout << "PostStepDoIt done " << G4endl; 476 return theTotalResult; 382 return theTotalResult; 477 } 383 } 478 384 >> 385 479 void G4HadronicProcess::ProcessDescription(std 386 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const 480 { 387 { 481 outFile << "The description for this process 388 outFile << "The description for this process has not been written yet.\n"; 482 } 389 } 483 390 >> 391 484 G4double G4HadronicProcess::XBiasSurvivalProba 392 G4double G4HadronicProcess::XBiasSurvivalProbability() 485 { 393 { >> 394 G4double result = 0; 486 G4double nLTraversed = GetTotalNumberOfInter 395 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed(); 487 G4double biasedProbability = 1.-G4Exp(-nLTra << 396 G4double biasedProbability = 1.-std::exp(-nLTraversed); 488 G4double realProbability = 1-G4Exp(-nLTraver << 397 G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor); 489 G4double result = (biasedProbability-realPro << 398 result = (biasedProbability-realProbability)/biasedProbability; 490 return result; 399 return result; 491 } 400 } 492 401 493 G4double G4HadronicProcess::XBiasSecondaryWeig 402 G4double G4HadronicProcess::XBiasSecondaryWeight() 494 { 403 { >> 404 G4double result = 0; 495 G4double nLTraversed = GetTotalNumberOfInter 405 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed(); 496 G4double result = << 406 result = 497 1./aScaleFactor*G4Exp(-nLTraversed/aScale << 407 1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor)); 498 return result; 408 return result; 499 } 409 } 500 410 501 void 411 void 502 G4HadronicProcess::FillResult(G4HadFinalState 412 G4HadronicProcess::FillResult(G4HadFinalState * aR, const G4Track & aT) 503 { 413 { 504 theTotalResult->ProposeLocalEnergyDeposit(aR 414 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit()); 505 const G4ThreeVector& dir = aT.GetMomentumDir << 506 415 507 G4double efinal = std::max(aR->GetEnergyChan << 416 G4double rotation = CLHEP::twopi*G4UniformRand(); >> 417 G4ThreeVector it(0., 0., 1.); >> 418 >> 419 G4double efinal = aR->GetEnergyChange(); >> 420 if(efinal < 0.0) { efinal = 0.0; } 508 421 509 // check status of primary 422 // check status of primary 510 if(aR->GetStatusChange() == stopAndKill) { 423 if(aR->GetStatusChange() == stopAndKill) { 511 theTotalResult->ProposeTrackStatus(fStopAn 424 theTotalResult->ProposeTrackStatus(fStopAndKill); 512 theTotalResult->ProposeEnergy( 0.0 ); 425 theTotalResult->ProposeEnergy( 0.0 ); 513 426 514 // check its final energy 427 // check its final energy 515 } else if(0.0 == efinal) { 428 } else if(0.0 == efinal) { 516 theTotalResult->ProposeEnergy( 0.0 ); 429 theTotalResult->ProposeEnergy( 0.0 ); 517 if(aT.GetParticleDefinition()->GetProcessM 430 if(aT.GetParticleDefinition()->GetProcessManager() 518 ->GetAtRestProcessVector()->size() > 0) 431 ->GetAtRestProcessVector()->size() > 0) 519 { theTotalResult->ProposeTrackStatus( 432 { theTotalResult->ProposeTrackStatus(fStopButAlive); } 520 else { theTotalResult->ProposeTrackStatus( 433 else { theTotalResult->ProposeTrackStatus(fStopAndKill); } 521 434 522 // primary is not killed apply rotation an 435 // primary is not killed apply rotation and Lorentz transformation 523 } else { 436 } else { 524 theTotalResult->ProposeTrackStatus(fAlive) 437 theTotalResult->ProposeTrackStatus(fAlive); 525 G4ThreeVector newDir = aR->GetMomentumChan << 438 G4double mass = aT.GetParticleDefinition()->GetPDGMass(); 526 newDir.rotateUz(dir); << 439 G4double newE = efinal + mass; 527 theTotalResult->ProposeMomentumDirection(n << 440 G4double newP = std::sqrt(efinal*(efinal + 2*mass)); 528 theTotalResult->ProposeEnergy(efinal); << 441 G4ThreeVector newPV = newP*aR->GetMomentumChange(); >> 442 G4LorentzVector newP4(newE, newPV); >> 443 newP4.rotate(rotation, it); >> 444 newP4 *= aR->GetTrafoToLab(); >> 445 theTotalResult->ProposeMomentumDirection(newP4.vect().unit()); >> 446 newE = newP4.e() - mass; >> 447 if(G4HadronicProcess_debug_flag && newE <= 0.0) { >> 448 G4ExceptionDescription ed; >> 449 DumpState(aT,"Primary has zero energy after interaction",ed); >> 450 G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning, ed); >> 451 } >> 452 if(newE < 0.0) { newE = 0.0; } >> 453 theTotalResult->ProposeEnergy( newE ); 529 } 454 } 530 //G4cout << "FillResult: Efinal= " << efinal 455 //G4cout << "FillResult: Efinal= " << efinal << " status= " 531 // << theTotalResult->GetTrackStatus() 456 // << theTotalResult->GetTrackStatus() 532 // << " fKill= " << fStopAndKill << G4end 457 // << " fKill= " << fStopAndKill << G4endl; 533 << 534 // check secondaries << 535 nICelectrons = 0; << 536 G4int nSec = (G4int)aR->GetNumberOfSecondari << 537 theTotalResult->SetNumberOfSecondaries(nSec) << 538 G4double time0 = aT.GetGlobalTime(); << 539 458 540 for (G4int i = 0; i < nSec; ++i) { << 459 // check secondaries: apply rotation and Lorentz transformation 541 G4DynamicParticle* dynParticle = aR->GetSe << 460 G4int nSec = aR->GetNumberOfSecondaries(); >> 461 theTotalResult->SetNumberOfSecondaries(nSec); >> 462 G4double weight = aT.GetWeight(); 542 463 543 // apply rotation << 464 if (nSec > 0) { 544 G4ThreeVector newDir = dynParticle->GetMom << 465 G4double time0 = aT.GetGlobalTime(); 545 newDir.rotateUz(dir); << 466 for (G4int i = 0; i < nSec; ++i) { 546 dynParticle->SetMomentumDirection(newDir); << 467 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum(); 547 << 468 theM.rotate(rotation, it); 548 // check if secondary is on the mass shell << 469 theM *= aR->GetTrafoToLab(); 549 const G4ParticleDefinition* part = dynPart << 470 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM); 550 G4double mass = part->GetPDGMass(); << 471 551 G4double dmass= dynParticle->GetMass(); << 472 // time of interaction starts from zero 552 const G4double delta_mass_lim = 1.0*CLHEP: << 473 G4double time = aR->GetSecondary(i)->GetTime(); 553 const G4double delta_ekin = 0.001*CLHEP::e << 474 if (time < 0.0) { time = 0.0; } 554 if(std::abs(dmass - mass) > delta_mass_lim << 475 555 G4double e = << 476 // take into account global time 556 std::max(dynParticle->GetKineticEnergy << 477 time += time0; 557 if(verboseLevel > 1) { << 478 558 G4ExceptionDescription ed; << 479 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(), 559 ed << "TrackID= "<< aT.GetTrackID() << 480 time, aT.GetPosition()); 560 << " " << aT.GetParticleDefinition()->Ge << 481 track->SetCreatorModelIndex(aR->GetSecondary(i)->GetCreatorModelType()); 561 << " Target Z= " << targetNucleus.GetZ_as << 482 G4double newWeight = weight*aR->GetSecondary(i)->GetWeight(); 562 << targetNucleus.GetA_asInt() << 483 // G4cout << "#### ParticleDebug " 563 << " Ekin(GeV)= " << aT.GetKineticEnergy( << 484 // <<GetProcessName()<<" " 564 << "\n Secondary is out of mass shell: " << 485 //<<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" " 565 << " EkinNew(MeV)= " << e << 486 // <<aScaleFactor<<" " 566 << " DeltaMass(MeV)= " << dmass - mass << << 487 // <<XBiasSurvivalProbability()<<" " 567 G4Exception("G4HadronicProcess::FillResults" << 488 // <<XBiasSecondaryWeight()<<" " >> 489 // <<aT.GetWeight()<<" " >> 490 // <<aR->GetSecondary(i)->GetWeight()<<" " >> 491 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" " >> 492 // <<G4endl; >> 493 track->SetWeight(newWeight); >> 494 track->SetTouchableHandle(aT.GetTouchableHandle()); >> 495 theTotalResult->AddSecondary(track); >> 496 if (G4HadronicProcess_debug_flag) { >> 497 G4double e = track->GetKineticEnergy(); >> 498 if (e <= 0.0) { >> 499 G4ExceptionDescription ed; >> 500 DumpState(aT,"Secondary has zero energy",ed); >> 501 ed << "Secondary " << track->GetDefinition()->GetParticleName() >> 502 << G4endl; >> 503 G4Exception("G4HadronicProcess::FillResults", "had011", >> 504 JustWarning,ed); >> 505 } 568 } 506 } 569 dynParticle->SetKineticEnergy(e); << 570 dynParticle->SetMass(mass); << 571 } 507 } 572 G4int idModel = aR->GetSecondary(i)->GetCr << 573 if(part->GetPDGEncoding() == 11) { ++nICel << 574 << 575 // time of interaction starts from zero + << 576 G4double time = std::max(aR->GetSecondary( << 577 << 578 G4Track* track = new G4Track(dynParticle, << 579 track->SetCreatorModelID(idModel); << 580 track->SetParentResonanceDef(aR->GetSecond << 581 track->SetParentResonanceID(aR->GetSeconda << 582 G4double newWeight = fWeight*aR->GetSecond << 583 track->SetWeight(newWeight); << 584 track->SetTouchableHandle(aT.GetTouchableH << 585 theTotalResult->AddSecondary(track); << 586 } 508 } 587 aR->Clear(); << 588 // G4cout << "FillResults done nICe= " << nI << 589 } << 590 509 591 void G4HadronicProcess::MultiplyCrossSectionBy << 510 aR->Clear(); 592 { << 593 BiasCrossSectionByFactor(factor); << 594 } 511 } 595 512 596 void G4HadronicProcess::BiasCrossSectionByFact 513 void G4HadronicProcess::BiasCrossSectionByFactor(G4double aScale) 597 { 514 { 598 if (aScale <= 0.0) { << 515 xBiasOn = true; >> 516 aScaleFactor = aScale; >> 517 G4String it = GetProcessName(); >> 518 if ((it != "photonNuclear") && >> 519 (it != "electronNuclear") && >> 520 (it != "positronNuclear") ) { 599 G4ExceptionDescription ed; 521 G4ExceptionDescription ed; 600 ed << " Wrong biasing factor " << aScale < << 522 G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had009", 601 G4Exception("G4HadronicProcess::BiasCrossS << 523 FatalException, ed, 602 JustWarning, ed, "Cross-sectio << 524 "Cross-section biasing available only for gamma and electro nuclear reactions."); 603 } else { << 525 } 604 aScaleFactor = aScale; << 526 >> 527 if (aScale < 100) { >> 528 G4ExceptionDescription ed; >> 529 G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010", JustWarning,ed, >> 530 "Cross-section bias readjusted to be above safe limit. New value is 100"); >> 531 aScaleFactor = 100.; 605 } 532 } 606 } 533 } 607 534 608 G4HadFinalState* G4HadronicProcess::CheckResul 535 G4HadFinalState* G4HadronicProcess::CheckResult(const G4HadProjectile & aPro, 609 const G4Nucleus &aNucleus, 536 const G4Nucleus &aNucleus, 610 G4HadFinalState * result) 537 G4HadFinalState * result) 611 { 538 { 612 // check for catastrophic energy non-conserv 539 // check for catastrophic energy non-conservation 613 // to re-sample the interaction 540 // to re-sample the interaction >> 541 614 G4HadronicInteraction * theModel = GetHadron 542 G4HadronicInteraction * theModel = GetHadronicInteraction(); 615 G4double nuclearMass(0); 543 G4double nuclearMass(0); 616 if (nullptr != theModel) { << 544 if (theModel) { 617 545 618 // Compute final-state total energy 546 // Compute final-state total energy 619 G4double finalE(0.); 547 G4double finalE(0.); 620 G4int nSec = (G4int)result->GetNumberOfSec << 548 G4int nSec = result->GetNumberOfSecondaries(); 621 549 622 nuclearMass = G4NucleiProperties::GetNucle 550 nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(), 623 aNucleus.GetZ_asInt()); 551 aNucleus.GetZ_asInt()); 624 if (result->GetStatusChange() != stopAndKi 552 if (result->GetStatusChange() != stopAndKill) { 625 // Interaction didn't complete, returned 553 // Interaction didn't complete, returned "do nothing" state 626 // and reset nucleus or the primary surv 554 // and reset nucleus or the primary survived the interaction 627 // (e.g. electro-nuclear ) => keep nucl 555 // (e.g. electro-nuclear ) => keep nucleus 628 finalE=result->GetLocalEnergyDeposit() + 556 finalE=result->GetLocalEnergyDeposit() + 629 aPro.GetDefinition()->GetPDGMass( 557 aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange(); 630 if( nSec == 0 ){ 558 if( nSec == 0 ){ 631 // Since there are no secondaries, th 559 // Since there are no secondaries, there is no recoil nucleus. 632 // To check energy balance we must ne 560 // To check energy balance we must neglect the initial nucleus too. 633 nuclearMass=0.0; 561 nuclearMass=0.0; 634 } 562 } 635 } 563 } 636 for (G4int i = 0; i < nSec; ++i) { << 564 for (G4int i = 0; i < nSec; i++) { 637 G4DynamicParticle *pdyn=result->GetSecon 565 G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle(); 638 finalE += pdyn->GetTotalEnergy(); 566 finalE += pdyn->GetTotalEnergy(); 639 G4double mass_pdg=pdyn->GetDefinition()- 567 G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass(); 640 G4double mass_dyn=pdyn->GetMass(); 568 G4double mass_dyn=pdyn->GetMass(); 641 if ( std::abs(mass_pdg - mass_dyn) > 0.1 << 569 if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV){ 642 // If it is shortlived, then a differe << 570 result->Clear(); 643 if ( pdyn->GetDefinition()->IsShortLiv << 571 result = 0; 644 std::abs(mass_pdg - mass_dyn) < 3 << 572 G4ExceptionDescription desc; 645 continue; << 573 desc << "Warning: Secondary with off-shell dynamic mass detected: " << G4endl 646 } << 574 << " " << pdyn->GetDefinition()->GetParticleName() 647 result->Clear(); << 575 << ", PDG mass: " << mass_pdg << ", dynamic mass: "<< mass_dyn << G4endl 648 result = nullptr; << 576 << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl 649 G4ExceptionDescription desc; << 577 << " Process / Model: " << GetProcessName()<< " / " 650 desc << "Warning: Secondary with off-shell d << 578 << theModel->GetModelName() << G4endl 651 << G4endl << 579 << " Primary: " << aPro.GetDefinition()->GetParticleName() 652 << " " << pdyn->GetDefinition()->GetPar << 580 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), " 653 << ", PDG mass: " << mass_pdg << ", dyn << 581 << " E= " << aPro.Get4Momentum().e() 654 << mass_dyn << G4endl << 582 << ", target nucleus (" << aNucleus.GetZ_asInt() << ", " 655 << (epReportLevel<0 ? "abort the event" << 583 << aNucleus.GetA_asInt() << ")" << G4endl; 656 : "re-sample the interaction") << G4endl << 584 G4Exception("G4HadronicProcess:CheckResult()", "had012", 657 << " Process / Model: " << GetProcessN << 585 epReportLevel<0 ? EventMustBeAborted : JustWarning,desc); 658 << theModel->GetModelName() << G4endl << 586 // must return here..... 659 << " Primary: " << aPro.GetDefinition() << 587 return result; 660 << " (" << aPro.GetDefinition()->GetPDG << 661 << " E= " << aPro.Get4Momentum().e() << 662 << ", target nucleus (" << aNucleus.Get << 663 << aNucleus.GetA_asInt() << ")" << G4en << 664 G4Exception("G4HadronicProcess:CheckResult() << 665 epReportLevel<0 ? EventMustBeAborted : << 666 // must return here..... << 667 return result; << 668 } 588 } 669 } 589 } 670 G4double deltaE= nuclearMass + aPro.GetTo 590 G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE; 671 591 672 std::pair<G4double, G4double> checkLevels 592 std::pair<G4double, G4double> checkLevels = 673 theModel->GetFatalEnergyCheckLevels(); 593 theModel->GetFatalEnergyCheckLevels(); // (relative, absolute) 674 if (std::abs(deltaE) > checkLevels.second 594 if (std::abs(deltaE) > checkLevels.second && 675 std::abs(deltaE) > checkLevels.first*a 595 std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){ 676 // do not delete result, this is a point 596 // do not delete result, this is a pointer to a data member; 677 result->Clear(); 597 result->Clear(); 678 result = nullptr; << 598 result = 0; 679 G4ExceptionDescription desc; 599 G4ExceptionDescription desc; 680 desc << "Warning: Bad energy non-conserv 600 desc << "Warning: Bad energy non-conservation detected, will " 681 << (epReportLevel<0 ? "abort the event" << 601 << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl 682 : "re-sample the interaction") << G4 << 683 << " Process / Model: " << GetProcessNam 602 << " Process / Model: " << GetProcessName()<< " / " 684 << theModel->GetModelName() << G4endl 603 << theModel->GetModelName() << G4endl 685 << " Primary: " << aPro.GetDefinition()-> 604 << " Primary: " << aPro.GetDefinition()->GetParticleName() 686 << " (" << aPro.GetDefinition()->GetPDGEn 605 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), " 687 << " E= " << aPro.Get4Momentum().e() 606 << " E= " << aPro.Get4Momentum().e() 688 << ", target nucleus (" << aNucleus.GetZ_ 607 << ", target nucleus (" << aNucleus.GetZ_asInt() << ", " 689 << aNucleus.GetA_asInt() << ")" << G4endl 608 << aNucleus.GetA_asInt() << ")" << G4endl 690 << " E(initial - final) = " << deltaE << 609 << " E(initial - final) = " << deltaE << " MeV." << G4endl; 691 G4Exception("G4HadronicProcess:CheckResu 610 G4Exception("G4HadronicProcess:CheckResult()", "had012", 692 epReportLevel<0 ? EventMustBeAborted : J 611 epReportLevel<0 ? EventMustBeAborted : JustWarning,desc); 693 } 612 } 694 } 613 } 695 return result; 614 return result; 696 } 615 } 697 616 698 void 617 void 699 G4HadronicProcess::CheckEnergyMomentumConserva 618 G4HadronicProcess::CheckEnergyMomentumConservation(const G4Track& aTrack, 700 619 const G4Nucleus& aNucleus) 701 { 620 { 702 G4int target_A=aNucleus.GetA_asInt(); 621 G4int target_A=aNucleus.GetA_asInt(); 703 G4int target_Z=aNucleus.GetZ_asInt(); 622 G4int target_Z=aNucleus.GetZ_asInt(); 704 G4double targetMass = G4NucleiProperties::Ge 623 G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z); 705 G4LorentzVector target4mom(0, 0, 0, targetMa << 624 G4LorentzVector target4mom(0, 0, 0, targetMass); 706 + nICelectrons*CLHEP::electron_mass << 707 625 708 G4LorentzVector projectile4mom = aTrack.GetD 626 G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum(); 709 G4int track_A = aTrack.GetDefinition()->GetB 627 G4int track_A = aTrack.GetDefinition()->GetBaryonNumber(); 710 G4int track_Z = G4lrint(aTrack.GetDefinition 628 G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge()); 711 629 712 G4int initial_A = target_A + track_A; 630 G4int initial_A = target_A + track_A; 713 G4int initial_Z = target_Z + track_Z - nICel << 631 G4int initial_Z = target_Z + track_Z; 714 632 715 G4LorentzVector initial4mom = projectile4mom 633 G4LorentzVector initial4mom = projectile4mom + target4mom; 716 634 717 // Compute final-state momentum for scatteri 635 // Compute final-state momentum for scattering and "do nothing" results 718 G4LorentzVector final4mom; 636 G4LorentzVector final4mom; 719 G4int final_A(0), final_Z(0); 637 G4int final_A(0), final_Z(0); 720 638 721 G4int nSec = theTotalResult->GetNumberOfSeco 639 G4int nSec = theTotalResult->GetNumberOfSecondaries(); 722 if (theTotalResult->GetTrackStatus() != fSto 640 if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive 723 // Either interaction didn't complete, ret << 641 // Either interaction didn't complete, returned "do nothing" state 724 // or the primary survived the interac << 642 // or the primary survived the interaction (e.g. electro-nucleus ) 725 << 643 G4Track temp(aTrack); 726 // Interaction didn't complete, returned " << 644 727 // - or suppressed recoil (e.g. Neutron << 645 // Use the final energy / momentum 728 final4mom = initial4mom; << 646 temp.SetMomentumDirection(*theTotalResult->GetMomentumDirection()); 729 final_A = initial_A; << 647 temp.SetKineticEnergy(theTotalResult->GetEnergy()); 730 final_Z = initial_Z; << 648 731 if (nSec > 0) { << 649 if( nSec == 0 ){ 732 // The primary remains in final state (e << 650 // Interaction didn't complete, returned "do nothing" state 733 // Use the final energy / momentum << 651 // - or suppressed recoil (e.g. Neutron elastic ) 734 const G4ThreeVector& v = *theTotalResult << 652 final4mom = temp.GetDynamicParticle()->Get4Momentum() + target4mom; 735 G4double ekin = theTotalResult->GetEnerg << 653 final_A = initial_A; 736 G4double mass = aTrack.GetDefinition()-> << 654 final_Z = initial_Z; 737 G4double ptot = std::sqrt(ekin*(ekin + 2 << 655 }else{ 738 final4mom.set(ptot*v.x(), ptot*v.y(), pt << 656 // The primary remains in final state (e.g. electro-nucleus ) 739 final_A = track_A; << 657 final4mom = temp.GetDynamicParticle()->Get4Momentum(); 740 final_Z = track_Z; << 658 final_A = track_A; 741 // Expect that the target nucleus will h << 659 final_Z = track_Z; 742 // and its products, including recoil, << 660 // Expect that the target nucleus will have interacted, 743 } << 661 // and its products, including recoil, will be included in secondaries. >> 662 } 744 } 663 } 745 if( nSec > 0 ) { 664 if( nSec > 0 ) { 746 G4Track* sec; 665 G4Track* sec; 747 666 748 for (G4int i = 0; i < nSec; i++) { 667 for (G4int i = 0; i < nSec; i++) { 749 sec = theTotalResult->GetSecondary(i); 668 sec = theTotalResult->GetSecondary(i); 750 final4mom += sec->GetDynamicParticle()-> 669 final4mom += sec->GetDynamicParticle()->Get4Momentum(); 751 final_A += sec->GetDefinition()->GetBary 670 final_A += sec->GetDefinition()->GetBaryonNumber(); 752 final_Z += G4lrint(sec->GetDefinition()- 671 final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge()); 753 } 672 } 754 } 673 } 755 674 756 // Get level-checking information (used to c 675 // Get level-checking information (used to cut-off relative checks) 757 G4String processName = GetProcessName(); 676 G4String processName = GetProcessName(); 758 G4HadronicInteraction* theModel = GetHadroni 677 G4HadronicInteraction* theModel = GetHadronicInteraction(); 759 G4String modelName("none"); 678 G4String modelName("none"); 760 if (theModel) modelName = theModel->GetModel 679 if (theModel) modelName = theModel->GetModelName(); 761 std::pair<G4double, G4double> checkLevels = 680 std::pair<G4double, G4double> checkLevels = epCheckLevels; 762 if (!levelsSetByProcess) { 681 if (!levelsSetByProcess) { 763 if (theModel) checkLevels = theModel->GetE 682 if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels(); 764 checkLevels.first= std::min(checkLevels.fi 683 checkLevels.first= std::min(checkLevels.first, epCheckLevels.first); 765 checkLevels.second=std::min(checkLevels.se 684 checkLevels.second=std::min(checkLevels.second, epCheckLevels.second); 766 } 685 } 767 686 768 // Compute absolute total-energy difference, 687 // Compute absolute total-energy difference, and relative kinetic-energy 769 G4bool checkRelative = (aTrack.GetKineticEne 688 G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second); 770 689 771 G4LorentzVector diff = initial4mom - final4m 690 G4LorentzVector diff = initial4mom - final4mom; 772 G4double absolute = diff.e(); 691 G4double absolute = diff.e(); 773 G4double relative = checkRelative ? absolute 692 G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.; 774 693 775 G4double absolute_mom = diff.vect().mag(); 694 G4double absolute_mom = diff.vect().mag(); 776 G4double relative_mom = checkRelative ? abso 695 G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.; 777 696 778 // Evaluate relative and absolute conservati 697 // Evaluate relative and absolute conservation 779 G4bool relPass = true; 698 G4bool relPass = true; 780 G4String relResult = "pass"; 699 G4String relResult = "pass"; 781 if ( std::abs(relative) > checkLevels.first 700 if ( std::abs(relative) > checkLevels.first 782 || std::abs(relative_mom) > checkLevels.fir 701 || std::abs(relative_mom) > checkLevels.first) { 783 relPass = false; 702 relPass = false; 784 relResult = checkRelative ? "fail" : "N/A" 703 relResult = checkRelative ? "fail" : "N/A"; 785 } 704 } 786 705 787 G4bool absPass = true; 706 G4bool absPass = true; 788 G4String absResult = "pass"; 707 G4String absResult = "pass"; 789 if ( std::abs(absolute) > checkLevels.seco 708 if ( std::abs(absolute) > checkLevels.second 790 || std::abs(absolute_mom) > checkLevels. 709 || std::abs(absolute_mom) > checkLevels.second ) { 791 absPass = false ; 710 absPass = false ; 792 absResult = "fail"; 711 absResult = "fail"; 793 } 712 } 794 713 795 G4bool chargePass = true; 714 G4bool chargePass = true; 796 G4String chargeResult = "pass"; 715 G4String chargeResult = "pass"; 797 if ( (initial_A-final_A)!=0 716 if ( (initial_A-final_A)!=0 798 || (initial_Z-final_Z)!=0 ) { 717 || (initial_Z-final_Z)!=0 ) { 799 chargePass = checkLevels.second < DBL_MAX 718 chargePass = checkLevels.second < DBL_MAX ? false : true; 800 chargeResult = "fail"; 719 chargeResult = "fail"; 801 } 720 } 802 721 803 G4bool conservationPass = (relPass || absPas 722 G4bool conservationPass = (relPass || absPass) && chargePass; 804 723 805 std::stringstream Myout; 724 std::stringstream Myout; 806 G4bool Myout_notempty(false); 725 G4bool Myout_notempty(false); 807 // Options for level of reporting detail: 726 // Options for level of reporting detail: 808 // 0. off 727 // 0. off 809 // 1. report only when E/p not conserved 728 // 1. report only when E/p not conserved 810 // 2. report regardless of E/p conservation 729 // 2. report regardless of E/p conservation 811 // 3. report only when E/p not conserved, w 730 // 3. report only when E/p not conserved, with model names, process names, and limits 812 // 4. report regardless of E/p conservation 731 // 4. report regardless of E/p conservation, with model names, process names, and limits 813 // negative -1.., as above, but send output 732 // negative -1.., as above, but send output to stderr 814 733 815 if( std::abs(epReportLevel) == 4 734 if( std::abs(epReportLevel) == 4 816 || ( std::abs(epReportLevel) == 3 && ! cons 735 || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){ 817 Myout << " Process: " << processName << 736 Myout << " Process: " << processName << " , Model: " << modelName << G4endl; 818 Myout << " Primary: " << aTrack.GetParti 737 Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName() 819 << " (" << aTrack.GetParticleDefin 738 << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << ")," 820 << " E= " << aTrack.GetDynamicPar 739 << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e() 821 << ", target nucleus (" << aNucleus.GetZ 740 << ", target nucleus (" << aNucleus.GetZ_asInt() << "," 822 << aNucleus.GetA_asInt() << ")" << G4end 741 << aNucleus.GetA_asInt() << ")" << G4endl; 823 Myout_notempty=true; 742 Myout_notempty=true; 824 } 743 } 825 if ( std::abs(epReportLevel) == 4 744 if ( std::abs(epReportLevel) == 4 826 || std::abs(epReportLevel) == 2 745 || std::abs(epReportLevel) == 2 827 || ! conservationPass ){ 746 || ! conservationPass ){ 828 747 829 Myout << " "<< relResult <<" relative 748 Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = " 830 << relative << " p/p(0)= " << rel 749 << relative << " p/p(0)= " << relative_mom << G4endl; 831 Myout << " "<< absResult << " absolute 750 Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = " 832 << absolute/MeV << " / " << absol 751 << absolute/MeV << " / " << absolute_mom/MeV << " 3mom: " << (diff.vect())*1./MeV << G4endl; 833 Myout << " "<< chargeResult << " charg 752 Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl; 834 Myout_notempty=true; 753 Myout_notempty=true; 835 754 836 } 755 } 837 Myout.flush(); 756 Myout.flush(); 838 if ( Myout_notempty ) { 757 if ( Myout_notempty ) { 839 if (epReportLevel > 0) G4cout << Myo 758 if (epReportLevel > 0) G4cout << Myout.str()<< G4endl; 840 else if (epReportLevel < 0) G4cerr << Myo 759 else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl; 841 } 760 } 842 } 761 } 843 762 >> 763 844 void G4HadronicProcess::DumpState(const G4Trac 764 void G4HadronicProcess::DumpState(const G4Track& aTrack, 845 const G4String& method, 765 const G4String& method, 846 G4ExceptionDescription& ed) 766 G4ExceptionDescription& ed) 847 { 767 { 848 ed << "Unrecoverable error in the method " < 768 ed << "Unrecoverable error in the method " << method << " of " 849 << GetProcessName() << G4endl; 769 << GetProcessName() << G4endl; 850 ed << "TrackID= "<< aTrack.GetTrackID() << " 770 ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= " 851 << aTrack.GetParentID() 771 << aTrack.GetParentID() 852 << " " << aTrack.GetParticleDefinition() 772 << " " << aTrack.GetParticleDefinition()->GetParticleName() 853 << G4endl; 773 << G4endl; 854 ed << "Ekin(GeV)= " << aTrack.GetKineticEner 774 ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV 855 << "; direction= " << aTrack.GetMomentum 775 << "; direction= " << aTrack.GetMomentumDirection() << G4endl; 856 ed << "Position(mm)= " << aTrack.GetPosition 776 ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";"; 857 777 858 if (aTrack.GetMaterial()) { 778 if (aTrack.GetMaterial()) { 859 ed << " material " << aTrack.GetMaterial( 779 ed << " material " << aTrack.GetMaterial()->GetName(); 860 } 780 } 861 ed << G4endl; 781 ed << G4endl; 862 782 863 if (aTrack.GetVolume()) { 783 if (aTrack.GetVolume()) { 864 ed << "PhysicalVolume <" << aTrack.GetVol 784 ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName() 865 << ">" << G4endl; 785 << ">" << G4endl; 866 } 786 } 867 } << 868 << 869 void G4HadronicProcess::DumpPhysicsTable(const << 870 { << 871 theCrossSectionDataStore->DumpPhysicsTable(p << 872 } << 873 << 874 void G4HadronicProcess::AddDataSet(G4VCrossSec << 875 { << 876 theCrossSectionDataStore->AddDataSet(aDataSe << 877 } << 878 << 879 std::vector<G4HadronicInteraction*>& << 880 G4HadronicProcess::GetHadronicInteractionList( << 881 { << 882 return theEnergyRangeManager.GetHadronicInte << 883 } << 884 << 885 G4HadronicInteraction* << 886 G4HadronicProcess::GetHadronicModel(const G4St << 887 { << 888 std::vector<G4HadronicInteraction*>& list << 889 = theEnergyRangeManager.GetHadronicInt << 890 for (auto & mod : list) { << 891 if (mod->GetModelName() == modelName) retu << 892 } << 893 return nullptr; << 894 } << 895 << 896 G4double << 897 G4HadronicProcess::ComputeCrossSection(const G << 898 const G4Material* mat, << 899 const G4double kinEnergy) << 900 { << 901 auto dp = new G4DynamicParticle(part, unitVe << 902 G4double xs = theCrossSectionDataStore->Comp << 903 delete dp; << 904 return xs; << 905 } << 906 << 907 void G4HadronicProcess::RecomputeXSandMFP(cons << 908 { << 909 auto dp = new G4DynamicParticle(currentParti << 910 theLastCrossSection = aScaleFactor* << 911 theCrossSectionDataStore->ComputeCrossSect << 912 theMFP = (theLastCrossSection > 0.0) ? 1.0/t << 913 delete dp; << 914 } << 915 << 916 void G4HadronicProcess::UpdateCrossSectionAndM << 917 { << 918 if(fXSType == fHadNoIntegral) { << 919 DefineXSandMFP(); << 920 << 921 } else if(fXSType == fHadIncreasing) { << 922 if(e*invLambdaFactor < mfpKinEnergy) { << 923 mfpKinEnergy = e; << 924 ComputeXSandMFP(); << 925 } << 926 << 927 } else if(fXSType == fHadDecreasing) { << 928 if(e < mfpKinEnergy && mfpKinEnergy > minK << 929 G4double e1 = std::max(e*lambdaFactor, m << 930 mfpKinEnergy = e1; << 931 RecomputeXSandMFP(e1); << 932 } << 933 << 934 } else if(fXSType == fHadOnePeak) { << 935 G4double epeak = (*theEnergyOfCrossSection << 936 if(e <= epeak) { << 937 if(e*invLambdaFactor < mfpKinEnergy) { << 938 mfpKinEnergy = e; << 939 ComputeXSandMFP(); << 940 } << 941 } else if(e < mfpKinEnergy) { << 942 G4double e1 = std::max(epeak, e*lambdaFa << 943 mfpKinEnergy = e1; << 944 RecomputeXSandMFP(e1); << 945 } << 946 << 947 } else if(fXSType == fHadTwoPeaks) { << 948 G4TwoPeaksHadXS* xs = (*fXSpeaks)[matIdx]; << 949 const G4double e1peak = xs->e1peak; << 950 << 951 // below the 1st peak << 952 if(e <= e1peak) { << 953 if(e*invLambdaFactor < mfpKinEnergy) { << 954 mfpKinEnergy = e; << 955 ComputeXSandMFP(); << 956 } << 957 return; << 958 } << 959 const G4double e1deep = xs->e1deep; << 960 // above the 1st peak, below the deep << 961 if(e <= e1deep) { << 962 if(mfpKinEnergy >= e1deep || e <= mfpKin << 963 const G4double e1 = std::max(e1peak, e << 964 mfpKinEnergy = e1; << 965 RecomputeXSandMFP(e1); << 966 } << 967 return; << 968 } << 969 const G4double e2peak = xs->e2peak; << 970 // above the deep, below 2nd peak << 971 if(e <= e2peak) { << 972 if(e*invLambdaFactor < mfpKinEnergy) { << 973 mfpKinEnergy = e; << 974 ComputeXSandMFP(); << 975 } << 976 return; << 977 } << 978 const G4double e2deep = xs->e2deep; << 979 // above the 2nd peak, below the deep << 980 if(e <= e2deep) { << 981 if(mfpKinEnergy >= e2deep || e <= mfpKin << 982 const G4double e1 = std::max(e2peak, e << 983 mfpKinEnergy = e1; << 984 RecomputeXSandMFP(e1); << 985 } << 986 return; << 987 } << 988 const G4double e3peak = xs->e3peak; << 989 // above the deep, below 3d peak << 990 if(e <= e3peak) { << 991 if(e*invLambdaFactor < mfpKinEnergy) { << 992 mfpKinEnergy = e; << 993 ComputeXSandMFP(); << 994 } << 995 return; << 996 } << 997 // above 3d peak << 998 if(e <= mfpKinEnergy) { << 999 const G4double e1 = std::max(e3peak, e*l << 1000 mfpKinEnergy = e1; << 1001 RecomputeXSandMFP(e1); << 1002 } << 1003 << 1004 } else { << 1005 DefineXSandMFP(); << 1006 } << 1007 } 787 } 1008 788