| 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$
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 <<
69 #include "G4NistManager.hh" <<
70 #include "G4VLeadingParticleBiasing.hh" <<
71 #include "G4HadXSHelper.hh" <<
72 #include "G4Threading.hh" <<
73 #include "G4Exp.hh" <<
74 70
75 #include <typeinfo> 71 #include <typeinfo>
76 #include <sstream> 72 #include <sstream>
77 #include <iostream> 73 #include <iostream>
78 74
79 namespace << 75 #include <stdlib.h>
80 { << 76
81 constexpr G4double lambdaFactor = 0.8; << 77 // File-scope variable to capture environment variable at startup
82 constexpr G4double invLambdaFactor = 1.0/lam << 78
83 } << 79 static const char* G4Hadronic_Random_File = getenv("G4HADRONIC_RANDOM_FILE");
84 80
85 ////////////////////////////////////////////// 81 //////////////////////////////////////////////////////////////////
86 82
87 G4HadronicProcess::G4HadronicProcess(const G4S 83 G4HadronicProcess::G4HadronicProcess(const G4String& processName,
88 G4Process 84 G4ProcessType procType)
89 : G4VDiscreteProcess(processName, procType) 85 : G4VDiscreteProcess(processName, procType)
90 { 86 {
91 SetProcessSubType(fHadronInelastic); // Def 87 SetProcessSubType(fHadronInelastic); // Default unless subclass changes
92 InitialiseLocal(); << 88
>> 89 theTotalResult = new G4ParticleChange();
>> 90 theTotalResult->SetSecondaryWeightByProcess(true);
>> 91 theInteraction = 0;
>> 92 theCrossSectionDataStore = new G4CrossSectionDataStore();
>> 93 G4HadronicProcessStore::Instance()->Register(this);
>> 94 aScaleFactor = 1;
>> 95 xBiasOn = false;
>> 96 G4HadronicProcess_debug_flag = false;
>> 97
>> 98 GetEnergyMomentumCheckEnvvars();
93 } 99 }
94 100
>> 101 //////////////////////////////////////////////////////////////////
>> 102
95 G4HadronicProcess::G4HadronicProcess(const G4S 103 G4HadronicProcess::G4HadronicProcess(const G4String& processName,
96 G4Hadroni 104 G4HadronicProcessType aHadSubType)
97 : G4VDiscreteProcess(processName, fHadronic) 105 : G4VDiscreteProcess(processName, fHadronic)
98 { 106 {
99 SetProcessSubType(aHadSubType); 107 SetProcessSubType(aHadSubType);
100 InitialiseLocal(); << 108
>> 109 theTotalResult = new G4ParticleChange();
>> 110 theTotalResult->SetSecondaryWeightByProcess(true);
>> 111 theInteraction = 0;
>> 112 theCrossSectionDataStore = new G4CrossSectionDataStore();
>> 113 G4HadronicProcessStore::Instance()->Register(this);
>> 114 aScaleFactor = 1;
>> 115 xBiasOn = false;
>> 116 G4HadronicProcess_debug_flag = false;
>> 117
>> 118 GetEnergyMomentumCheckEnvvars();
101 } 119 }
102 120
>> 121
103 G4HadronicProcess::~G4HadronicProcess() 122 G4HadronicProcess::~G4HadronicProcess()
104 { 123 {
105 theProcessStore->DeRegister(this); << 124 G4HadronicProcessStore::Instance()->DeRegister(this);
106 delete theTotalResult; 125 delete theTotalResult;
107 delete theCrossSectionDataStore; 126 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 } 127 }
118 128
119 void G4HadronicProcess::InitialiseLocal() { << 129 void G4HadronicProcess::GetEnergyMomentumCheckEnvvars() {
120 theTotalResult = new G4ParticleChange(); << 130 levelsSetByProcess = false;
121 theTotalResult->SetSecondaryWeightByProcess( << 131
122 theCrossSectionDataStore = new G4CrossSectio << 132 epReportLevel = getenv("G4Hadronic_epReportLevel") ?
123 theProcessStore = G4HadronicProcessStore::In << 133 strtol(getenv("G4Hadronic_epReportLevel"),0,10) : 0;
124 theProcessStore->Register(this); << 134
125 minKinEnergy = 1*CLHEP::MeV; << 135 epCheckLevels.first = getenv("G4Hadronic_epCheckRelativeLevel") ?
126 << 136 strtod(getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX;
127 G4HadronicParameters* param = G4HadronicPara <<
128 epReportLevel = param->GetEPReportLevel(); <<
129 epCheckLevels.first = param->GetEPRelativeLe <<
130 epCheckLevels.second = param->GetEPAbsoluteL <<
131 137
132 unitVector.set(0.0, 0.0, 0.1); << 138 epCheckLevels.second = getenv("G4Hadronic_epCheckAbsoluteLevel") ?
133 if(G4Threading::IsWorkerThread()) { isMaster << 139 strtod(getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX;
134 } 140 }
135 141
136 void G4HadronicProcess::RegisterMe( G4Hadronic 142 void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
137 { 143 {
138 if(nullptr == a) { return; } << 144 if(!a) { return; }
139 theEnergyRangeManager.RegisterMe( a ); << 145 try{GetManagerPointer()->RegisterMe( a );}
140 G4HadronicProcessStore::Instance()->Register << 146 catch(G4HadronicException & aE)
141 } << 147 {
142 << 148 G4ExceptionDescription ed;
143 G4double << 149 aE.Report(ed);
144 G4HadronicProcess::GetElementCrossSection(cons << 150 ed << "Unrecoverable error in " << GetProcessName()
145 const G4Element * elm, << 151 << " to register " << a->GetModelName() << G4endl;
146 const G4Material* mat) << 152 G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException,
147 { << 153 ed);
148 if(nullptr == mat) <<
149 { <<
150 static const G4int nmax = 5; <<
151 if(nMatWarn < nmax) { <<
152 ++nMatWarn; <<
153 G4ExceptionDescription ed; <<
154 ed << "Cannot compute Element x-section <<
155 << " because no material defined \n" <<
156 << " Please, specify material pointer or de <<
157 << " for Z= " << elm->GetZasInt(); <<
158 G4Exception("G4HadronicProcess::GetEleme <<
159 JustWarning, ed); <<
160 } <<
161 } 154 }
162 return theCrossSectionDataStore->GetCrossSec << 155 G4HadronicProcessStore::Instance()->RegisterInteraction(this, a);
163 } 156 }
164 157
165 void G4HadronicProcess::PreparePhysicsTable(co 158 void G4HadronicProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
166 { 159 {
167 if(nullptr == firstParticle) { firstParticle << 160 if(getenv("G4HadronicProcess_debug")) {
168 theProcessStore->RegisterParticle(this, &p); << 161 G4HadronicProcess_debug_flag = true;
>> 162 }
>> 163 G4HadronicProcessStore::Instance()->RegisterParticle(this, &p);
169 } 164 }
170 165
171 void G4HadronicProcess::BuildPhysicsTable(cons 166 void G4HadronicProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
172 { 167 {
173 if(firstParticle != &p) { return; } << 168 try
174 << 169 {
175 theCrossSectionDataStore->BuildPhysicsTable( << 170 theCrossSectionDataStore->BuildPhysicsTable(p);
176 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 } 171 }
204 << 172 catch(G4HadronicException aR)
205 // check particle for integral method << 173 {
206 if(isMaster || nullptr == masterProcess) { << 174 G4ExceptionDescription ed;
207 G4double charge = p.GetPDGCharge()/eplus; << 175 aR.Report(ed);
208 << 176 ed << " hadronic initialisation fails" << G4endl;
209 // select cross section shape << 177 G4Exception("G4HadronicProcess::BuildPhysicsTable", "had000",
210 if(charge != 0.0 && useIntegralXS) { << 178 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 } 179 }
264 G4HadronicProcessStore::Instance()->PrintInf 180 G4HadronicProcessStore::Instance()->PrintInfo(&p);
265 } 181 }
266 182
267 void G4HadronicProcess::StartTracking(G4Track* << 183 G4double G4HadronicProcess::
268 { << 184 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 { 185 {
280 *condition = NotForced; << 186 try
281 << 187 {
282 const G4Material* mat = track.GetMaterial(); << 188 theLastCrossSection = aScaleFactor*
283 if(mat != currentMat) { << 189 theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(),
284 currentMat = mat; << 190 aTrack.GetMaterial());
285 mfpKinEnergy = DBL_MAX; <<
286 matIdx = (G4int)track.GetMaterial()->GetIn <<
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 } 191 }
307 currentInteractionLength = theMFP; << 192 catch(G4HadronicException aR)
308 return theNumberOfInteractionLengthLeft*theM << 193 {
309 } << 194 G4ExceptionDescription ed;
310 << 195 aR.Report(ed);
311 G4double G4HadronicProcess::GetMeanFreePath( << 196 DumpState(aTrack,"GetMeanFreePath",ed);
312 const G4Track &aTr << 197 ed << " Cross section is not available" << G4endl;
313 G4ForceCondition*) << 198 G4Exception("G4HadronicProcess::GetMeanFreePath", "had002", FatalException,
314 { << 199 ed);
315 G4double xs = aScaleFactor*theCrossSectionDa << 200 }
316 ->ComputeCrossSection(aTrack.GetDynamicPa << 201 G4double res = DBL_MAX;
317 return (xs > 0.0) ? 1.0/xs : DBL_MAX; << 202 if( theLastCrossSection > 0.0 ) { res = 1.0/theLastCrossSection; }
>> 203 return res;
318 } 204 }
319 205
320 G4VParticleChange* 206 G4VParticleChange*
321 G4HadronicProcess::PostStepDoIt(const G4Track& 207 G4HadronicProcess::PostStepDoIt(const G4Track& aTrack, const G4Step&)
322 { 208 {
323 theNumberOfInteractionLengthLeft = -1.0; <<
324 <<
325 //G4cout << "PostStepDoIt " << aTrack.GetDef <<
326 // << " Ekin= " << aTrack.GetKineticEnergy <<
327 // if primary is not Alive then do nothing 209 // if primary is not Alive then do nothing
328 theTotalResult->Clear(); 210 theTotalResult->Clear();
329 theTotalResult->Initialize(aTrack); 211 theTotalResult->Initialize(aTrack);
330 fWeight = aTrack.GetWeight(); << 212 theTotalResult->ProposeWeight(aTrack.GetWeight());
331 theTotalResult->ProposeWeight(fWeight); <<
332 if(aTrack.GetTrackStatus() != fAlive) { retu 213 if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
333 214
334 // Find cross section at end of step and che 215 // Find cross section at end of step and check if <= 0
335 // 216 //
336 const G4DynamicParticle* aParticle = aTrack. 217 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
337 const G4Material* aMaterial = aTrack.GetMate << 218 G4Material* aMaterial = aTrack.GetMaterial();
>> 219
>> 220 G4Element* anElement = 0;
>> 221 try
>> 222 {
>> 223 anElement = theCrossSectionDataStore->SampleZandA(aParticle,
>> 224 aMaterial,
>> 225 targetNucleus);
>> 226 }
>> 227 catch(G4HadronicException & aR)
>> 228 {
>> 229 G4ExceptionDescription ed;
>> 230 aR.Report(ed);
>> 231 DumpState(aTrack,"SampleZandA",ed);
>> 232 ed << " PostStepDoIt failed on element selection" << G4endl;
>> 233 G4Exception("G4HadronicProcess::PostStepDoIt", "had003", FatalException,
>> 234 ed);
>> 235 }
338 236
339 // check only for charged particles 237 // check only for charged particles
340 if(fXSType != fHadNoIntegral) { << 238 if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
341 mfpKinEnergy = DBL_MAX; << 239 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 240 // No interaction
348 return theTotalResult; 241 return theTotalResult;
349 } 242 }
350 } 243 }
351 244
352 const G4Element* anElement = <<
353 theCrossSectionDataStore->SampleZandA(aPar <<
354 <<
355 // Next check for illegal track status 245 // Next check for illegal track status
356 // 246 //
357 if (aTrack.GetTrackStatus() != fAlive && << 247 if (aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
358 aTrack.GetTrackStatus() != fSuspend) { <<
359 if (aTrack.GetTrackStatus() == fStopAndKil 248 if (aTrack.GetTrackStatus() == fStopAndKill ||
360 aTrack.GetTrackStatus() == fKillTrackA 249 aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
361 aTrack.GetTrackStatus() == fPostponeTo 250 aTrack.GetTrackStatus() == fPostponeToNextEvent) {
362 G4ExceptionDescription ed; 251 G4ExceptionDescription ed;
363 ed << "G4HadronicProcess: track in unusa 252 ed << "G4HadronicProcess: track in unusable state - "
364 << aTrack.GetTrackStatus() << G4endl; 253 << aTrack.GetTrackStatus() << G4endl;
365 ed << "G4HadronicProcess: returning unch 254 ed << "G4HadronicProcess: returning unchanged track " << G4endl;
366 DumpState(aTrack,"PostStepDoIt",ed); 255 DumpState(aTrack,"PostStepDoIt",ed);
367 G4Exception("G4HadronicProcess::PostStep 256 G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
368 } 257 }
369 // No warning for fStopButAlive which is a 258 // No warning for fStopButAlive which is a legal status here
370 return theTotalResult; 259 return theTotalResult;
371 } 260 }
372 261
373 // Initialize the hadronic projectile from t << 262 // Go on to regular case
374 thePro.Initialise(aTrack); << 263 //
>> 264 G4double originalEnergy = aParticle->GetKineticEnergy();
>> 265 G4double kineticEnergy = originalEnergy;
>> 266
>> 267 // Get kinetic energy per nucleon for ions
>> 268 if(aParticle->GetParticleDefinition()->GetBaryonNumber() > 1.5)
>> 269 kineticEnergy/=aParticle->GetParticleDefinition()->GetBaryonNumber();
375 270
376 theInteraction = ChooseHadronicInteraction(t << 271 try
377 a << 272 {
378 if(nullptr == theInteraction) { << 273 theInteraction =
>> 274 ChooseHadronicInteraction( kineticEnergy, aMaterial, anElement );
>> 275 }
>> 276 catch(G4HadronicException & aE)
>> 277 {
379 G4ExceptionDescription ed; 278 G4ExceptionDescription ed;
>> 279 aE.Report(ed);
380 ed << "Target element "<<anElement->GetNam 280 ed << "Target element "<<anElement->GetName()<<" Z= "
381 << targetNucleus.GetZ_asInt() << " A= 281 << targetNucleus.GetZ_asInt() << " A= "
382 << targetNucleus.GetA_asInt() << G4endl 282 << targetNucleus.GetA_asInt() << G4endl;
383 DumpState(aTrack,"ChooseHadronicInteractio 283 DumpState(aTrack,"ChooseHadronicInteraction",ed);
384 ed << " No HadronicInteraction found out" 284 ed << " No HadronicInteraction found out" << G4endl;
385 G4Exception("G4HadronicProcess::PostStepDo << 285 G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException,
386 FatalException, ed); << 286 ed);
387 return theTotalResult; <<
388 } 287 }
389 288
390 G4HadFinalState* result = nullptr; << 289 // Initialize the hadronic projectile from the track
>> 290 thePro.Initialise(aTrack);
>> 291 G4HadFinalState* result = 0;
391 G4int reentryCount = 0; 292 G4int reentryCount = 0;
392 /* << 293
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 294 do
401 { 295 {
402 try 296 try
403 { 297 {
>> 298 // Save random engine if requested for debugging
>> 299 if (G4Hadronic_Random_File) {
>> 300 CLHEP::HepRandom::saveEngineStatus(G4Hadronic_Random_File);
>> 301 }
404 // Call the interaction 302 // Call the interaction
405 result = theInteraction->ApplyYourself( 303 result = theInteraction->ApplyYourself( thePro, targetNucleus);
406 ++reentryCount; 304 ++reentryCount;
407 } 305 }
408 catch(G4HadronicException & aR) << 306 catch(G4HadronicException aR)
409 { 307 {
410 G4ExceptionDescription ed; 308 G4ExceptionDescription ed;
411 aR.Report(ed); 309 aR.Report(ed);
412 ed << "Call for " << theInteraction->Get 310 ed << "Call for " << theInteraction->GetModelName() << G4endl;
413 ed << "Target element "<<anElement->GetN 311 ed << "Target element "<<anElement->GetName()<<" Z= "
414 << targetNucleus.GetZ_asInt() 312 << targetNucleus.GetZ_asInt()
415 << " A= " << targetNucleus.GetA_asInt() << 313 << " A= " << targetNucleus.GetA_asInt() << G4endl;
416 DumpState(aTrack,"ApplyYourself",ed); 314 DumpState(aTrack,"ApplyYourself",ed);
417 ed << " ApplyYourself failed" << G4endl; 315 ed << " ApplyYourself failed" << G4endl;
418 G4Exception("G4HadronicProcess::PostStep 316 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
419 ed); 317 ed);
420 } 318 }
421 319
422 // Check the result for catastrophic energ 320 // Check the result for catastrophic energy non-conservation
423 result = CheckResult(thePro, targetNucleus << 321 result = CheckResult(thePro,targetNucleus, result);
424 322
425 if(reentryCount>100) { 323 if(reentryCount>100) {
426 G4ExceptionDescription ed; 324 G4ExceptionDescription ed;
427 ed << "Call for " << theInteraction->Get 325 ed << "Call for " << theInteraction->GetModelName() << G4endl;
428 ed << "Target element "<<anElement->GetN 326 ed << "Target element "<<anElement->GetName()<<" Z= "
429 << targetNucleus.GetZ_asInt() 327 << targetNucleus.GetZ_asInt()
430 << " A= " << targetNucleus.GetA_asInt() << 328 << " A= " << targetNucleus.GetA_asInt() << G4endl;
431 DumpState(aTrack,"ApplyYourself",ed); 329 DumpState(aTrack,"ApplyYourself",ed);
432 ed << " ApplyYourself does not completed 330 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
433 G4Exception("G4HadronicProcess::PostStep 331 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
434 ed); 332 ed);
435 } 333 }
436 } 334 }
437 while(!result); /* Loop checking, 30-Oct-20 << 335 while(!result);
438 <<
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 336
469 result->SetTrafoToLab(thePro.GetTrafoToLab() 337 result->SetTrafoToLab(thePro.GetTrafoToLab());
>> 338
>> 339 ClearNumberOfInteractionLengthLeft();
>> 340
470 FillResult(result, aTrack); 341 FillResult(result, aTrack);
471 342
472 if (epReportLevel != 0) { 343 if (epReportLevel != 0) {
473 CheckEnergyMomentumConservation(aTrack, ta 344 CheckEnergyMomentumConservation(aTrack, targetNucleus);
474 } 345 }
475 //G4cout << "PostStepDoIt done nICelectrons= <<
476 return theTotalResult; 346 return theTotalResult;
477 } 347 }
478 348
>> 349
479 void G4HadronicProcess::ProcessDescription(std 350 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const
480 { 351 {
481 outFile << "The description for this process 352 outFile << "The description for this process has not been written yet.\n";
482 } 353 }
483 354
>> 355
484 G4double G4HadronicProcess::XBiasSurvivalProba 356 G4double G4HadronicProcess::XBiasSurvivalProbability()
485 { 357 {
>> 358 G4double result = 0;
486 G4double nLTraversed = GetTotalNumberOfInter 359 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
487 G4double biasedProbability = 1.-G4Exp(-nLTra << 360 G4double biasedProbability = 1.-std::exp(-nLTraversed);
488 G4double realProbability = 1-G4Exp(-nLTraver << 361 G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
489 G4double result = (biasedProbability-realPro << 362 result = (biasedProbability-realProbability)/biasedProbability;
490 return result; 363 return result;
491 } 364 }
492 365
493 G4double G4HadronicProcess::XBiasSecondaryWeig 366 G4double G4HadronicProcess::XBiasSecondaryWeight()
494 { 367 {
>> 368 G4double result = 0;
495 G4double nLTraversed = GetTotalNumberOfInter 369 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
496 G4double result = << 370 result =
497 1./aScaleFactor*G4Exp(-nLTraversed/aScale << 371 1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
498 return result; 372 return result;
499 } 373 }
500 374
501 void 375 void
502 G4HadronicProcess::FillResult(G4HadFinalState 376 G4HadronicProcess::FillResult(G4HadFinalState * aR, const G4Track & aT)
503 { 377 {
504 theTotalResult->ProposeLocalEnergyDeposit(aR 378 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
505 const G4ThreeVector& dir = aT.GetMomentumDir <<
506 379
507 G4double efinal = std::max(aR->GetEnergyChan << 380 G4double rotation = CLHEP::twopi*G4UniformRand();
>> 381 G4ThreeVector it(0., 0., 1.);
>> 382
>> 383 G4double efinal = aR->GetEnergyChange();
>> 384 if(efinal < 0.0) { efinal = 0.0; }
508 385
509 // check status of primary 386 // check status of primary
510 if(aR->GetStatusChange() == stopAndKill) { 387 if(aR->GetStatusChange() == stopAndKill) {
511 theTotalResult->ProposeTrackStatus(fStopAn 388 theTotalResult->ProposeTrackStatus(fStopAndKill);
512 theTotalResult->ProposeEnergy( 0.0 ); 389 theTotalResult->ProposeEnergy( 0.0 );
513 390
514 // check its final energy 391 // check its final energy
515 } else if(0.0 == efinal) { 392 } else if(0.0 == efinal) {
516 theTotalResult->ProposeEnergy( 0.0 ); 393 theTotalResult->ProposeEnergy( 0.0 );
517 if(aT.GetParticleDefinition()->GetProcessM 394 if(aT.GetParticleDefinition()->GetProcessManager()
518 ->GetAtRestProcessVector()->size() > 0) 395 ->GetAtRestProcessVector()->size() > 0)
519 { theTotalResult->ProposeTrackStatus( << 396 { aParticleChange.ProposeTrackStatus(fStopButAlive); }
520 else { theTotalResult->ProposeTrackStatus( << 397 else { aParticleChange.ProposeTrackStatus(fStopAndKill); }
521 398
522 // primary is not killed apply rotation an 399 // primary is not killed apply rotation and Lorentz transformation
523 } else { 400 } else {
524 theTotalResult->ProposeTrackStatus(fAlive) 401 theTotalResult->ProposeTrackStatus(fAlive);
525 G4ThreeVector newDir = aR->GetMomentumChan << 402 G4double mass = aT.GetParticleDefinition()->GetPDGMass();
526 newDir.rotateUz(dir); << 403 G4double newE = efinal + mass;
527 theTotalResult->ProposeMomentumDirection(n << 404 G4double newP = std::sqrt(efinal*(efinal + 2*mass));
528 theTotalResult->ProposeEnergy(efinal); << 405 G4ThreeVector newPV = newP*aR->GetMomentumChange();
529 } << 406 G4LorentzVector newP4(newE, newPV);
530 //G4cout << "FillResult: Efinal= " << efinal << 407 newP4.rotate(rotation, it);
531 // << theTotalResult->GetTrackStatus() << 408 newP4 *= aR->GetTrafoToLab();
532 // << " fKill= " << fStopAndKill << G4end << 409 theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
533 << 410 newE = newP4.e() - mass;
534 // check secondaries << 411 if(G4HadronicProcess_debug_flag && newE <= 0.0) {
535 nICelectrons = 0; << 412 G4ExceptionDescription ed;
536 G4int nSec = (G4int)aR->GetNumberOfSecondari << 413 DumpState(aT,"Primary has zero energy after interaction",ed);
>> 414 G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning, ed);
>> 415 }
>> 416 if(newE < 0.0) { newE = 0.0; }
>> 417 theTotalResult->ProposeEnergy( newE );
>> 418 }
>> 419
>> 420 // check secondaries: apply rotation and Lorentz transformation
>> 421 G4int nSec = aR->GetNumberOfSecondaries();
537 theTotalResult->SetNumberOfSecondaries(nSec) 422 theTotalResult->SetNumberOfSecondaries(nSec);
538 G4double time0 = aT.GetGlobalTime(); << 423 G4double weight = aT.GetWeight();
>> 424
>> 425 if (nSec > 0) {
>> 426 G4double time0 = aT.GetGlobalTime();
>> 427 for (G4int i = 0; i < nSec; ++i) {
>> 428 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
>> 429 theM.rotate(rotation, it);
>> 430 theM *= aR->GetTrafoToLab();
>> 431 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
>> 432
>> 433 // time of interaction starts from zero
>> 434 G4double time = aR->GetSecondary(i)->GetTime();
>> 435 if (time < 0.0) { time = 0.0; }
>> 436
>> 437 // take into account global time
>> 438 time += time0;
>> 439
>> 440 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
>> 441 time, aT.GetPosition());
>> 442 G4double newWeight = weight*aR->GetSecondary(i)->GetWeight();
>> 443 // G4cout << "#### ParticleDebug "
>> 444 // <<GetProcessName()<<" "
>> 445 // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
>> 446 // <<aScaleFactor<<" "
>> 447 // <<XBiasSurvivalProbability()<<" "
>> 448 // <<XBiasSecondaryWeight()<<" "
>> 449 // <<aT.GetWeight()<<" "
>> 450 // <<aR->GetSecondary(i)->GetWeight()<<" "
>> 451 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
>> 452 // <<G4endl;
>> 453 track->SetWeight(newWeight);
>> 454 track->SetTouchableHandle(aT.GetTouchableHandle());
>> 455 theTotalResult->AddSecondary(track);
>> 456 if (G4HadronicProcess_debug_flag) {
>> 457 G4double e = track->GetKineticEnergy();
>> 458 if (e <= 0.0) {
>> 459 G4ExceptionDescription ed;
>> 460 DumpState(aT,"Secondary has zero energy",ed);
>> 461 ed << "Secondary " << track->GetDefinition()->GetParticleName()
>> 462 << G4endl;
>> 463 G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning,ed);
>> 464 }
>> 465 }
>> 466 }
>> 467 }
539 468
540 for (G4int i = 0; i < nSec; ++i) { << 469 aR->Clear();
541 G4DynamicParticle* dynParticle = aR->GetSe << 470 return;
>> 471 }
>> 472 /*
>> 473 void
>> 474 G4HadronicProcess::FillTotalResult(G4HadFinalState* aR, const G4Track& aT)
>> 475 {
>> 476 theTotalResult->Clear();
>> 477 theTotalResult->ProposeLocalEnergyDeposit(0.);
>> 478 theTotalResult->Initialize(aT);
>> 479 theTotalResult->SetSecondaryWeightByProcess(true);
>> 480 theTotalResult->ProposeTrackStatus(fAlive);
>> 481 G4double rotation = CLHEP::twopi*G4UniformRand();
>> 482 G4ThreeVector it(0., 0., 1.);
542 483
543 // apply rotation << 484 if(aR->GetStatusChange()==stopAndKill)
544 G4ThreeVector newDir = dynParticle->GetMom << 485 {
545 newDir.rotateUz(dir); << 486 if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
546 dynParticle->SetMomentumDirection(newDir); << 487 {
547 << 488 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
548 // check if secondary is on the mass shell << 489 }
549 const G4ParticleDefinition* part = dynPart << 490 else
550 G4double mass = part->GetPDGMass(); << 491 {
551 G4double dmass= dynParticle->GetMass(); << 492 theTotalResult->ProposeTrackStatus(fStopAndKill);
552 const G4double delta_mass_lim = 1.0*CLHEP: << 493 theTotalResult->ProposeEnergy( 0.0 );
553 const G4double delta_ekin = 0.001*CLHEP::e << 494 }
554 if(std::abs(dmass - mass) > delta_mass_lim << 495 }
555 G4double e = << 496 else if(aR->GetStatusChange()!=stopAndKill )
556 std::max(dynParticle->GetKineticEnergy << 497 {
557 if(verboseLevel > 1) { << 498 if(aR->GetStatusChange()==suspend)
>> 499 {
>> 500 theTotalResult->ProposeTrackStatus(fSuspend);
>> 501 if(xBiasOn)
>> 502 {
558 G4ExceptionDescription ed; 503 G4ExceptionDescription ed;
559 ed << "TrackID= "<< aT.GetTrackID() << 504 DumpState(aT,"FillTotalResult",ed);
560 << " " << aT.GetParticleDefinition()->Ge << 505 G4Exception("G4HadronicProcess::FillTotalResult", "had007", FatalException,
561 << " Target Z= " << targetNucleus.GetZ_as << 506 ed,"Cannot cross-section bias a process that suspends tracks.");
562 << targetNucleus.GetA_asInt() <<
563 << " Ekin(GeV)= " << aT.GetKineticEnergy( <<
564 << "\n Secondary is out of mass shell: " <<
565 << " EkinNew(MeV)= " << e <<
566 << " DeltaMass(MeV)= " << dmass - mass << <<
567 G4Exception("G4HadronicProcess::FillResults" <<
568 } 507 }
569 dynParticle->SetKineticEnergy(e); << 508 } else if (aT.GetKineticEnergy() == 0) {
570 dynParticle->SetMass(mass); << 509 theTotalResult->ProposeTrackStatus(fStopButAlive);
>> 510 }
>> 511
>> 512 if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
>> 513 {
>> 514 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
>> 515 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
>> 516 G4double newM=aT.GetParticleDefinition()->GetPDGMass();
>> 517 G4double newE=aR->GetEnergyChange() + newM;
>> 518 G4double newP=std::sqrt(newE*newE - newM*newM);
>> 519 G4DynamicParticle * aNew =
>> 520 new G4DynamicParticle(aT.GetParticleDefinition(), newE, newP*aR->GetMomentumChange());
>> 521 aR->AddSecondary(G4HadSecondary(aNew, newWeight));
571 } 522 }
572 G4int idModel = aR->GetSecondary(i)->GetCr << 523 else
573 if(part->GetPDGEncoding() == 11) { ++nICel << 524 {
574 << 525 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
575 // time of interaction starts from zero + << 526 theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
576 G4double time = std::max(aR->GetSecondary( << 527 if(aR->GetEnergyChange()>-.5)
577 << 528 {
578 G4Track* track = new G4Track(dynParticle, << 529 theTotalResult->ProposeEnergy(aR->GetEnergyChange());
579 track->SetCreatorModelID(idModel); << 530 }
580 track->SetParentResonanceDef(aR->GetSecond << 531 G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
581 track->SetParentResonanceID(aR->GetSeconda << 532 newDirection*=aR->GetTrafoToLab();
582 G4double newWeight = fWeight*aR->GetSecond << 533 theTotalResult->ProposeMomentumDirection(newDirection.vect());
>> 534 }
>> 535 }
>> 536 else
>> 537 {
>> 538 G4ExceptionDescription ed;
>> 539 ed << "Call for " << theInteraction->GetModelName() << G4endl;
>> 540 ed << "Target Z= "
>> 541 << targetNucleus.GetZ_asInt()
>> 542 << " A= " << targetNucleus.GetA_asInt() << G4endl;
>> 543 DumpState(aT,"FillTotalResult",ed);
>> 544 G4Exception("G4HadronicProcess", "had008", FatalException,
>> 545 "use of unsupported track-status.");
>> 546 }
>> 547
>> 548 if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
>> 549 && theTotalResult->GetTrackStatus()==fAlive
>> 550 && aR->GetStatusChange()==isAlive)
>> 551 {
>> 552 // Use for debugging: G4double newWeight = theTotalResult->GetParentWeight();
>> 553
>> 554 G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
>> 555 G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetParticleDefinition(),
>> 556 aR->GetMomentumChange(),
>> 557 newKE);
>> 558 aR->AddSecondary(aNew);
>> 559 aR->SetStatusChange(stopAndKill);
>> 560
>> 561 theTotalResult->ProposeTrackStatus(fStopAndKill);
>> 562 theTotalResult->ProposeEnergy( 0.0 );
>> 563
>> 564 }
>> 565 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
>> 566 theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
>> 567
>> 568 if(aR->GetStatusChange() != stopAndKill)
>> 569 {
>> 570 G4double newM=aT.GetParticleDefinition()->GetPDGMass();
>> 571 G4double newE=aR->GetEnergyChange() + newM;
>> 572 G4double newP=std::sqrt(newE*newE - newM*newM);
>> 573 G4ThreeVector newPV = newP*aR->GetMomentumChange();
>> 574 G4LorentzVector newP4(newE, newPV);
>> 575 newP4.rotate(rotation, it);
>> 576 newP4*=aR->GetTrafoToLab();
>> 577 theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
>> 578 }
>> 579
>> 580 for(G4int i=0; i<aR->GetNumberOfSecondaries(); ++i)
>> 581 {
>> 582 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
>> 583 theM.rotate(rotation, it);
>> 584 theM*=aR->GetTrafoToLab();
>> 585 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
>> 586 G4double time = aR->GetSecondary(i)->GetTime();
>> 587 if(time<0) time = aT.GetGlobalTime();
>> 588
>> 589 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
>> 590 time,
>> 591 aT.GetPosition());
>> 592
>> 593 G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
>> 594 if(xBiasOn) { newWeight *= XBiasSecondaryWeight(); }
583 track->SetWeight(newWeight); 595 track->SetWeight(newWeight);
584 track->SetTouchableHandle(aT.GetTouchableH 596 track->SetTouchableHandle(aT.GetTouchableHandle());
585 theTotalResult->AddSecondary(track); 597 theTotalResult->AddSecondary(track);
586 } 598 }
587 aR->Clear(); <<
588 // G4cout << "FillResults done nICe= " << nI <<
589 } <<
590 599
591 void G4HadronicProcess::MultiplyCrossSectionBy << 600 aR->Clear();
592 { << 601 return;
593 BiasCrossSectionByFactor(factor); <<
594 } 602 }
>> 603 */
595 604
596 void G4HadronicProcess::BiasCrossSectionByFact 605 void G4HadronicProcess::BiasCrossSectionByFactor(G4double aScale)
597 { 606 {
598 if (aScale <= 0.0) { << 607 xBiasOn = true;
599 G4ExceptionDescription ed; << 608 aScaleFactor = aScale;
600 ed << " Wrong biasing factor " << aScale < << 609 G4String it = GetProcessName();
601 G4Exception("G4HadronicProcess::BiasCrossS << 610 if( (it != "PhotonInelastic") &&
602 JustWarning, ed, "Cross-sectio << 611 (it != "ElectroNuclear") &&
603 } else { << 612 (it != "PositronNuclear") )
604 aScaleFactor = aScale; << 613 {
605 } << 614 G4ExceptionDescription ed;
606 } << 615 G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had009", FatalException, ed,
607 << 616 "Cross-section biasing available only for gamma and electro nuclear reactions.");
608 G4HadFinalState* G4HadronicProcess::CheckResul <<
609 const G4Nucleus &aNucleus, <<
610 G4HadFinalState * result) <<
611 { <<
612 // check for catastrophic energy non-conserv <<
613 // to re-sample the interaction <<
614 G4HadronicInteraction * theModel = GetHadron <<
615 G4double nuclearMass(0); <<
616 if (nullptr != theModel) { <<
617 <<
618 // Compute final-state total energy <<
619 G4double finalE(0.); <<
620 G4int nSec = (G4int)result->GetNumberOfSec <<
621 <<
622 nuclearMass = G4NucleiProperties::GetNucle <<
623 aNucleus.GetZ_asInt()); <<
624 if (result->GetStatusChange() != stopAndKi <<
625 // Interaction didn't complete, returned <<
626 // and reset nucleus or the primary surv <<
627 // (e.g. electro-nuclear ) => keep nucl <<
628 finalE=result->GetLocalEnergyDeposit() + <<
629 aPro.GetDefinition()->GetPDGMass( <<
630 if( nSec == 0 ){ <<
631 // Since there are no secondaries, th <<
632 // To check energy balance we must ne <<
633 nuclearMass=0.0; <<
634 } <<
635 } 617 }
636 for (G4int i = 0; i < nSec; ++i) { << 618 if(aScale<100)
637 G4DynamicParticle *pdyn=result->GetSecon << 619 {
638 finalE += pdyn->GetTotalEnergy(); << 620 G4ExceptionDescription ed;
639 G4double mass_pdg=pdyn->GetDefinition()- << 621 G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010", JustWarning,ed,
640 G4double mass_dyn=pdyn->GetMass(); << 622 "Cross-section bias readjusted to be above safe limit. New value is 100");
641 if ( std::abs(mass_pdg - mass_dyn) > 0.1 << 623 aScaleFactor = 100.;
642 // If it is shortlived, then a differe << 624 }
643 if ( pdyn->GetDefinition()->IsShortLiv << 625 }
644 std::abs(mass_pdg - mass_dyn) < 3 << 626
645 continue; << 627 G4HadFinalState* G4HadronicProcess::CheckResult(const G4HadProjectile & aPro,const G4Nucleus &aNucleus, G4HadFinalState * result) const
646 } << 628 {
647 result->Clear(); << 629 // check for catastrophic energy non-conservation, to re-sample the interaction
648 result = nullptr; << 630
649 G4ExceptionDescription desc; << 631 G4HadronicInteraction * theModel = GetHadronicInteraction();
650 desc << "Warning: Secondary with off-shell d << 632 G4double nuclearMass(0);
651 << G4endl << 633 if (theModel){
652 << " " << pdyn->GetDefinition()->GetPar << 634
653 << ", PDG mass: " << mass_pdg << ", dyn << 635 // Compute final-state total energy
654 << mass_dyn << G4endl << 636 G4double finalE(0.);
655 << (epReportLevel<0 ? "abort the event" << 637 G4int nSec = result->GetNumberOfSecondaries();
656 : "re-sample the interaction") << G4endl << 638
657 << " Process / Model: " << GetProcessN << 639 nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
658 << theModel->GetModelName() << G4endl << 640 aNucleus.GetZ_asInt());
659 << " Primary: " << aPro.GetDefinition() << 641 if (result->GetStatusChange() != stopAndKill) {
660 << " (" << aPro.GetDefinition()->GetPDG << 642 // Interaction didn't complete, returned "do nothing" state => reset nucleus
661 << " E= " << aPro.Get4Momentum().e() << 643 // or the primary survived the interaction (e.g. electro-nuclear ) => keep nucleus
662 << ", target nucleus (" << aNucleus.Get << 644 finalE=result->GetLocalEnergyDeposit() +
663 << aNucleus.GetA_asInt() << ")" << G4en << 645 aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
664 G4Exception("G4HadronicProcess:CheckResult() << 646 if( nSec == 0 ){
665 epReportLevel<0 ? EventMustBeAborted : << 647 // Since there are no secondaries, there is no recoil nucleus.
666 // must return here..... << 648 // To check energy balance we must neglect the initial nucleus too.
667 return result; << 649 nuclearMass=0.0;
>> 650 }
>> 651 }
>> 652 for (G4int i = 0; i < nSec; i++) {
>> 653 G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle();
>> 654 finalE += pdyn->GetTotalEnergy();
>> 655 // check consistency of mass of dynamic particle with PGD value
>> 656 G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass();
>> 657 G4double mass_dyn=pdyn->GetMass();
>> 658 if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV){
>> 659 result->Clear();
>> 660 result = 0;
>> 661 G4ExceptionDescription desc;
>> 662 desc << "Warning: Secondary with off-shell dynamic mass detected: " << G4endl
>> 663 << " " << pdyn->GetDefinition()->GetParticleName()
>> 664 << ", PDG mass: " << mass_pdg << ", dynamic mass: "<< mass_dyn << G4endl
>> 665 << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl
>> 666 << " Process / Model: " << GetProcessName()<< " / "
>> 667 << theModel->GetModelName() << G4endl
>> 668 << " Primary: " << aPro.GetDefinition()->GetParticleName()
>> 669 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
>> 670 << " E= " << aPro.Get4Momentum().e()
>> 671 << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
>> 672 << aNucleus.GetA_asInt() << ")" << G4endl;
>> 673 G4Exception("G4HadronicProcess:CheckResult()", "had012",
>> 674 epReportLevel<0 ? EventMustBeAborted : JustWarning,desc);
>> 675 // must return here.....
>> 676 return result;
>> 677 }
>> 678 }
>> 679 G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
>> 680
>> 681 std::pair<G4double, G4double> checkLevels = theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
>> 682 if (std::abs(deltaE) > checkLevels.second && std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
>> 683 // do not delete result, this is a pointer to a data member;
>> 684 result->Clear();
>> 685 result=0;
>> 686 G4ExceptionDescription desc;
>> 687 desc << "Warning: Bad energy non-conservation detected, will "
>> 688 << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl
>> 689 << " Process / Model: " << GetProcessName()<< " / " << theModel->GetModelName() << G4endl
>> 690 << " Primary: " << aPro.GetDefinition()->GetParticleName()
>> 691 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "),"
>> 692 << " E= " << aPro.Get4Momentum().e()
>> 693 << ", target nucleus (" << aNucleus.GetZ_asInt() << ","<< aNucleus.GetA_asInt() << ")" << G4endl
>> 694 << " E(initial - final) = " << deltaE << " MeV." << G4endl;
>> 695 G4Exception("G4HadronicProcess:CheckResult()", "had012", epReportLevel<0 ? EventMustBeAborted : JustWarning,desc);
668 } 696 }
669 } << 697 }
670 G4double deltaE= nuclearMass + aPro.GetTo << 698 return result;
671 <<
672 std::pair<G4double, G4double> checkLevels <<
673 theModel->GetFatalEnergyCheckLevels(); <<
674 if (std::abs(deltaE) > checkLevels.second <<
675 std::abs(deltaE) > checkLevels.first*a <<
676 // do not delete result, this is a point <<
677 result->Clear(); <<
678 result = nullptr; <<
679 G4ExceptionDescription desc; <<
680 desc << "Warning: Bad energy non-conserv <<
681 << (epReportLevel<0 ? "abort the event" <<
682 : "re-sample the interaction") << G4 <<
683 << " Process / Model: " << GetProcessNam <<
684 << theModel->GetModelName() << G4endl <<
685 << " Primary: " << aPro.GetDefinition()-> <<
686 << " (" << aPro.GetDefinition()->GetPDGEn <<
687 << " E= " << aPro.Get4Momentum().e() <<
688 << ", target nucleus (" << aNucleus.GetZ_ <<
689 << aNucleus.GetA_asInt() << ")" << G4endl <<
690 << " E(initial - final) = " << deltaE << <<
691 G4Exception("G4HadronicProcess:CheckResu <<
692 epReportLevel<0 ? EventMustBeAborted : J <<
693 } <<
694 } <<
695 return result; <<
696 } 699 }
697 700
698 void 701 void
699 G4HadronicProcess::CheckEnergyMomentumConserva 702 G4HadronicProcess::CheckEnergyMomentumConservation(const G4Track& aTrack,
700 703 const G4Nucleus& aNucleus)
701 { 704 {
702 G4int target_A=aNucleus.GetA_asInt(); 705 G4int target_A=aNucleus.GetA_asInt();
703 G4int target_Z=aNucleus.GetZ_asInt(); 706 G4int target_Z=aNucleus.GetZ_asInt();
704 G4double targetMass = G4NucleiProperties::Ge 707 G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
705 G4LorentzVector target4mom(0, 0, 0, targetMa << 708 G4LorentzVector target4mom(0, 0, 0, targetMass);
706 + nICelectrons*CLHEP::electron_mass <<
707 709
708 G4LorentzVector projectile4mom = aTrack.GetD 710 G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
709 G4int track_A = aTrack.GetDefinition()->GetB 711 G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
710 G4int track_Z = G4lrint(aTrack.GetDefinition 712 G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
711 713
712 G4int initial_A = target_A + track_A; 714 G4int initial_A = target_A + track_A;
713 G4int initial_Z = target_Z + track_Z - nICel << 715 G4int initial_Z = target_Z + track_Z;
714 716
715 G4LorentzVector initial4mom = projectile4mom 717 G4LorentzVector initial4mom = projectile4mom + target4mom;
716 718
717 // Compute final-state momentum for scatteri 719 // Compute final-state momentum for scattering and "do nothing" results
718 G4LorentzVector final4mom; 720 G4LorentzVector final4mom;
719 G4int final_A(0), final_Z(0); 721 G4int final_A(0), final_Z(0);
720 722
721 G4int nSec = theTotalResult->GetNumberOfSeco 723 G4int nSec = theTotalResult->GetNumberOfSecondaries();
722 if (theTotalResult->GetTrackStatus() != fSto 724 if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
723 // Either interaction didn't complete, ret << 725 // Either interaction didn't complete, returned "do nothing" state
724 // or the primary survived the interac << 726 // or the primary survived the interaction (e.g. electro-nucleus )
725 << 727 G4Track temp(aTrack);
726 // Interaction didn't complete, returned " << 728
727 // - or suppressed recoil (e.g. Neutron << 729 // Use the final energy / momentum
728 final4mom = initial4mom; << 730 temp.SetMomentumDirection(*theTotalResult->GetMomentumDirection());
729 final_A = initial_A; << 731 temp.SetKineticEnergy(theTotalResult->GetEnergy());
730 final_Z = initial_Z; << 732
731 if (nSec > 0) { << 733 if( nSec == 0 ){
732 // The primary remains in final state (e << 734 // Interaction didn't complete, returned "do nothing" state
733 // Use the final energy / momentum << 735 // - or suppressed recoil (e.g. Neutron elastic )
734 const G4ThreeVector& v = *theTotalResult << 736 final4mom = temp.GetDynamicParticle()->Get4Momentum() + target4mom;
735 G4double ekin = theTotalResult->GetEnerg << 737 final_A = initial_A;
736 G4double mass = aTrack.GetDefinition()-> << 738 final_Z = initial_Z;
737 G4double ptot = std::sqrt(ekin*(ekin + 2 << 739 }else{
738 final4mom.set(ptot*v.x(), ptot*v.y(), pt << 740 // The primary remains in final state (e.g. electro-nucleus )
739 final_A = track_A; << 741 final4mom = temp.GetDynamicParticle()->Get4Momentum();
740 final_Z = track_Z; << 742 final_A = track_A;
741 // Expect that the target nucleus will h << 743 final_Z = track_Z;
742 // and its products, including recoil, << 744 // Expect that the target nucleus will have interacted,
743 } << 745 // and its products, including recoil, will be included in secondaries.
>> 746 }
744 } 747 }
745 if( nSec > 0 ) { 748 if( nSec > 0 ) {
746 G4Track* sec; 749 G4Track* sec;
747 750
748 for (G4int i = 0; i < nSec; i++) { 751 for (G4int i = 0; i < nSec; i++) {
749 sec = theTotalResult->GetSecondary(i); 752 sec = theTotalResult->GetSecondary(i);
750 final4mom += sec->GetDynamicParticle()-> 753 final4mom += sec->GetDynamicParticle()->Get4Momentum();
751 final_A += sec->GetDefinition()->GetBary 754 final_A += sec->GetDefinition()->GetBaryonNumber();
752 final_Z += G4lrint(sec->GetDefinition()- 755 final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
753 } 756 }
754 } 757 }
755 758
756 // Get level-checking information (used to c 759 // Get level-checking information (used to cut-off relative checks)
757 G4String processName = GetProcessName(); 760 G4String processName = GetProcessName();
758 G4HadronicInteraction* theModel = GetHadroni 761 G4HadronicInteraction* theModel = GetHadronicInteraction();
759 G4String modelName("none"); 762 G4String modelName("none");
760 if (theModel) modelName = theModel->GetModel 763 if (theModel) modelName = theModel->GetModelName();
761 std::pair<G4double, G4double> checkLevels = 764 std::pair<G4double, G4double> checkLevels = epCheckLevels;
762 if (!levelsSetByProcess) { 765 if (!levelsSetByProcess) {
763 if (theModel) checkLevels = theModel->GetE 766 if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
764 checkLevels.first= std::min(checkLevels.fi 767 checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
765 checkLevels.second=std::min(checkLevels.se 768 checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
766 } 769 }
767 770
768 // Compute absolute total-energy difference, 771 // Compute absolute total-energy difference, and relative kinetic-energy
769 G4bool checkRelative = (aTrack.GetKineticEne 772 G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
770 773
771 G4LorentzVector diff = initial4mom - final4m 774 G4LorentzVector diff = initial4mom - final4mom;
772 G4double absolute = diff.e(); 775 G4double absolute = diff.e();
773 G4double relative = checkRelative ? absolute 776 G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
774 777
775 G4double absolute_mom = diff.vect().mag(); 778 G4double absolute_mom = diff.vect().mag();
776 G4double relative_mom = checkRelative ? abso 779 G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
777 780
778 // Evaluate relative and absolute conservati 781 // Evaluate relative and absolute conservation
779 G4bool relPass = true; 782 G4bool relPass = true;
780 G4String relResult = "pass"; 783 G4String relResult = "pass";
781 if ( std::abs(relative) > checkLevels.first 784 if ( std::abs(relative) > checkLevels.first
782 || std::abs(relative_mom) > checkLevels.fir 785 || std::abs(relative_mom) > checkLevels.first) {
783 relPass = false; 786 relPass = false;
784 relResult = checkRelative ? "fail" : "N/A" 787 relResult = checkRelative ? "fail" : "N/A";
785 } 788 }
786 789
787 G4bool absPass = true; 790 G4bool absPass = true;
788 G4String absResult = "pass"; 791 G4String absResult = "pass";
789 if ( std::abs(absolute) > checkLevels.seco 792 if ( std::abs(absolute) > checkLevels.second
790 || std::abs(absolute_mom) > checkLevels. 793 || std::abs(absolute_mom) > checkLevels.second ) {
791 absPass = false ; 794 absPass = false ;
792 absResult = "fail"; 795 absResult = "fail";
793 } 796 }
794 797
795 G4bool chargePass = true; 798 G4bool chargePass = true;
796 G4String chargeResult = "pass"; 799 G4String chargeResult = "pass";
797 if ( (initial_A-final_A)!=0 800 if ( (initial_A-final_A)!=0
798 || (initial_Z-final_Z)!=0 ) { 801 || (initial_Z-final_Z)!=0 ) {
799 chargePass = checkLevels.second < DBL_MAX 802 chargePass = checkLevels.second < DBL_MAX ? false : true;
800 chargeResult = "fail"; 803 chargeResult = "fail";
801 } 804 }
802 805
803 G4bool conservationPass = (relPass || absPas 806 G4bool conservationPass = (relPass || absPass) && chargePass;
804 807
805 std::stringstream Myout; 808 std::stringstream Myout;
806 G4bool Myout_notempty(false); 809 G4bool Myout_notempty(false);
807 // Options for level of reporting detail: 810 // Options for level of reporting detail:
808 // 0. off 811 // 0. off
809 // 1. report only when E/p not conserved 812 // 1. report only when E/p not conserved
810 // 2. report regardless of E/p conservation 813 // 2. report regardless of E/p conservation
811 // 3. report only when E/p not conserved, w 814 // 3. report only when E/p not conserved, with model names, process names, and limits
812 // 4. report regardless of E/p conservation 815 // 4. report regardless of E/p conservation, with model names, process names, and limits
813 // negative -1.., as above, but send output 816 // negative -1.., as above, but send output to stderr
814 817
815 if( std::abs(epReportLevel) == 4 818 if( std::abs(epReportLevel) == 4
816 || ( std::abs(epReportLevel) == 3 && ! cons 819 || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
817 Myout << " Process: " << processName << 820 Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
818 Myout << " Primary: " << aTrack.GetParti 821 Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
819 << " (" << aTrack.GetParticleDefin 822 << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
820 << " E= " << aTrack.GetDynamicPar 823 << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
821 << ", target nucleus (" << aNucleus.GetZ 824 << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
822 << aNucleus.GetA_asInt() << ")" << G4end 825 << aNucleus.GetA_asInt() << ")" << G4endl;
823 Myout_notempty=true; 826 Myout_notempty=true;
824 } 827 }
825 if ( std::abs(epReportLevel) == 4 828 if ( std::abs(epReportLevel) == 4
826 || std::abs(epReportLevel) == 2 829 || std::abs(epReportLevel) == 2
827 || ! conservationPass ){ 830 || ! conservationPass ){
828 831
829 Myout << " "<< relResult <<" relative 832 Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
830 << relative << " p/p(0)= " << rel 833 << relative << " p/p(0)= " << relative_mom << G4endl;
831 Myout << " "<< absResult << " absolute 834 Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
832 << absolute/MeV << " / " << absol << 835 << absolute/MeV << " / " << absolute_mom/MeV << G4endl;
833 Myout << " "<< chargeResult << " charg 836 Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
834 Myout_notempty=true; 837 Myout_notempty=true;
835 838
836 } 839 }
837 Myout.flush(); 840 Myout.flush();
838 if ( Myout_notempty ) { 841 if ( Myout_notempty ) {
839 if (epReportLevel > 0) G4cout << Myo 842 if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
840 else if (epReportLevel < 0) G4cerr << Myo 843 else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
841 } 844 }
842 } 845 }
843 846
>> 847
844 void G4HadronicProcess::DumpState(const G4Trac 848 void G4HadronicProcess::DumpState(const G4Track& aTrack,
845 const G4String& method, 849 const G4String& method,
846 G4ExceptionDescription& ed) 850 G4ExceptionDescription& ed)
847 { 851 {
848 ed << "Unrecoverable error in the method " < 852 ed << "Unrecoverable error in the method " << method << " of "
849 << GetProcessName() << G4endl; 853 << GetProcessName() << G4endl;
850 ed << "TrackID= "<< aTrack.GetTrackID() << " 854 ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
851 << aTrack.GetParentID() 855 << aTrack.GetParentID()
852 << " " << aTrack.GetParticleDefinition() 856 << " " << aTrack.GetParticleDefinition()->GetParticleName()
853 << G4endl; 857 << G4endl;
854 ed << "Ekin(GeV)= " << aTrack.GetKineticEner 858 ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
855 << "; direction= " << aTrack.GetMomentum 859 << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
856 ed << "Position(mm)= " << aTrack.GetPosition 860 ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
857 861
858 if (aTrack.GetMaterial()) { 862 if (aTrack.GetMaterial()) {
859 ed << " material " << aTrack.GetMaterial( 863 ed << " material " << aTrack.GetMaterial()->GetName();
860 } 864 }
861 ed << G4endl; 865 ed << G4endl;
862 866
863 if (aTrack.GetVolume()) { 867 if (aTrack.GetVolume()) {
864 ed << "PhysicalVolume <" << aTrack.GetVol 868 ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
865 << ">" << G4endl; 869 << ">" << G4endl;
866 } 870 }
867 } 871 }
868 << 872 /*
869 void G4HadronicProcess::DumpPhysicsTable(const << 873 G4ParticleDefinition* G4HadronicProcess::GetTargetDefinition()
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 { 874 {
909 auto dp = new G4DynamicParticle(currentParti << 875 const G4Nucleus* nuc = GetTargetNucleus();
910 theLastCrossSection = aScaleFactor* << 876 G4int Z = nuc->GetZ_asInt();
911 theCrossSectionDataStore->ComputeCrossSect << 877 G4int A = nuc->GetA_asInt();
912 theMFP = (theLastCrossSection > 0.0) ? 1.0/t << 878 return G4ParticleTable::GetParticleTable()->GetIon(Z,A,0*eV);
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 } 879 }
>> 880 */
>> 881 /* end of file */
1008 882