| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class source file
30 //
31 // G4HadronicProcess
32 //
33 // original by H.P.Wellisch
34 // J.L. Chuma, TRIUMF, 10-Mar-1997
35 //
36 // Modifications:
37 // 05-Jul-2010 V.Ivanchenko cleanup commented lines
38 // 20-Jul-2011 M.Kelsey -- null-pointer checks in DumpState()
39 // 24-Sep-2011 M.Kelsey -- Use envvar G4HADRONIC_RANDOM_FILE to save random
40 // engine state before each model call
41 // 18-Oct-2011 M.Kelsey -- Handle final-state cases in conservation checks.
42 // 14-Mar-2012 G.Folger -- enhance checks for conservation of energy, etc.
43 // 28-Jul-2012 M.Maire -- add function GetTargetDefinition()
44 // 14-Sep-2012 Inherit from RestDiscrete, use subtype code (now in ctor) to
45 // configure base-class
46 // 28-Sep-2012 Restore inheritance from G4VDiscreteProcess, remove enable-flag
47 // changing, remove warning message from original ctor.
48 // 21-Aug-2019 V.Ivanchenko leave try/catch only for ApplyYourself(..), cleanup
49
50 #include "G4HadronicProcess.hh"
51
52 #include "G4Types.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4HadProjectile.hh"
55 #include "G4ElementVector.hh"
56 #include "G4Track.hh"
57 #include "G4Step.hh"
58 #include "G4Element.hh"
59 #include "G4ParticleChange.hh"
60 #include "G4ProcessVector.hh"
61 #include "G4ProcessManager.hh"
62 #include "G4NucleiProperties.hh"
63
64 #include "G4HadronicException.hh"
65 #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
75 #include <typeinfo>
76 #include <sstream>
77 #include <iostream>
78
79 namespace
80 {
81 constexpr G4double lambdaFactor = 0.8;
82 constexpr G4double invLambdaFactor = 1.0/lambdaFactor;
83 }
84
85 //////////////////////////////////////////////////////////////////
86
87 G4HadronicProcess::G4HadronicProcess(const G4String& processName,
88 G4ProcessType procType)
89 : G4VDiscreteProcess(processName, procType)
90 {
91 SetProcessSubType(fHadronInelastic); // Default unless subclass changes
92 InitialiseLocal();
93 }
94
95 G4HadronicProcess::G4HadronicProcess(const G4String& processName,
96 G4HadronicProcessType aHadSubType)
97 : G4VDiscreteProcess(processName, fHadronic)
98 {
99 SetProcessSubType(aHadSubType);
100 InitialiseLocal();
101 }
102
103 G4HadronicProcess::~G4HadronicProcess()
104 {
105 theProcessStore->DeRegister(this);
106 delete theTotalResult;
107 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 }
118
119 void G4HadronicProcess::InitialiseLocal() {
120 theTotalResult = new G4ParticleChange();
121 theTotalResult->SetSecondaryWeightByProcess(true);
122 theCrossSectionDataStore = new G4CrossSectionDataStore();
123 theProcessStore = G4HadronicProcessStore::Instance();
124 theProcessStore->Register(this);
125 minKinEnergy = 1*CLHEP::MeV;
126
127 G4HadronicParameters* param = G4HadronicParameters::Instance();
128 epReportLevel = param->GetEPReportLevel();
129 epCheckLevels.first = param->GetEPRelativeLevel();
130 epCheckLevels.second = param->GetEPAbsoluteLevel();
131
132 unitVector.set(0.0, 0.0, 0.1);
133 if(G4Threading::IsWorkerThread()) { isMaster = false; }
134 }
135
136 void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
137 {
138 if(nullptr == a) { return; }
139 theEnergyRangeManager.RegisterMe( a );
140 G4HadronicProcessStore::Instance()->RegisterInteraction(this, a);
141 }
142
143 G4double
144 G4HadronicProcess::GetElementCrossSection(const G4DynamicParticle * dp,
145 const G4Element * elm,
146 const G4Material* mat)
147 {
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 for " << GetProcessName()
155 << " because no material defined \n"
156 << " Please, specify material pointer or define simple material"
157 << " for Z= " << elm->GetZasInt();
158 G4Exception("G4HadronicProcess::GetElementCrossSection", "had066",
159 JustWarning, ed);
160 }
161 }
162 return theCrossSectionDataStore->GetCrossSection(dp, elm, mat);
163 }
164
165 void G4HadronicProcess::PreparePhysicsTable(const G4ParticleDefinition& p)
166 {
167 if(nullptr == firstParticle) { firstParticle = &p; }
168 theProcessStore->RegisterParticle(this, &p);
169 }
170
171 void G4HadronicProcess::BuildPhysicsTable(const G4ParticleDefinition& p)
172 {
173 if(firstParticle != &p) { return; }
174
175 theCrossSectionDataStore->BuildPhysicsTable(p);
176 theEnergyRangeManager.BuildPhysicsTable(p);
177 G4HadronicParameters* param = G4HadronicParameters::Instance();
178
179 G4int subtype = GetProcessSubType();
180 if(useIntegralXS) {
181 if(subtype == fHadronInelastic) {
182 useIntegralXS = param->EnableIntegralInelasticXS();
183 } else if(subtype == fHadronElastic) {
184 useIntegralXS = param->EnableIntegralElasticXS();
185 }
186 }
187 fXSType = fHadNoIntegral;
188
189 if(nullptr == masterProcess) {
190 masterProcess = dynamic_cast<const G4HadronicProcess*>(GetMasterProcess());
191 }
192 if(nullptr == masterProcess) {
193 if(1 < param->GetVerboseLevel()) {
194 G4ExceptionDescription ed;
195 ed << "G4HadronicProcess::BuildPhysicsTable: for "
196 << GetProcessName() << " for " << p.GetParticleName()
197 << " fail due to undefined pointer to the master process \n"
198 << " ThreadID= " << G4Threading::G4GetThreadId()
199 << " initialisation of worker started before master initialisation";
200 G4Exception("G4HadronicProcess::BuildPhysicsTable", "had066",
201 JustWarning, ed);
202 }
203 }
204
205 // check particle for integral method
206 if(isMaster || nullptr == masterProcess) {
207 G4double charge = p.GetPDGCharge()/eplus;
208
209 // select cross section shape
210 if(charge != 0.0 && useIntegralXS) {
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 == -1000010020 || pdg == -1000010030 ||
224 pdg == -1000020030 || pdg == -1000020040) {
225 fXSType = fHadDecreasing;
226 } else if (charge > 0.0 || pdg == 11 || pdg == 13) {
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, minKinEnergy, tmax);
241 if(nullptr == fXSpeaks) {
242 fXSType = fHadOnePeak;
243 }
244 }
245 if(fXSType == fHadOnePeak) {
246 theEnergyOfCrossSectionMax =
247 G4HadXSHelper::FindCrossSectionMax(this, &p, minKinEnergy, tmax);
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->EnergyOfCrossSectionMax();
258 }
259 if(isMaster && 1 < param->GetVerboseLevel()) {
260 G4cout << "G4HadronicProcess::BuildPhysicsTable: for "
261 << GetProcessName() << " and " << p.GetParticleName()
262 << " typeXS=" << fXSType << G4endl;
263 }
264 G4HadronicProcessStore::Instance()->PrintInfo(&p);
265 }
266
267 void G4HadronicProcess::StartTracking(G4Track* track)
268 {
269 currentMat = nullptr;
270 currentParticle = track->GetDefinition();
271 fDynParticle = track->GetDynamicParticle();
272 theNumberOfInteractionLengthLeft = -1.0;
273 }
274
275 G4double G4HadronicProcess::PostStepGetPhysicalInteractionLength(
276 const G4Track& track,
277 G4double previousStepSize,
278 G4ForceCondition* condition)
279 {
280 *condition = NotForced;
281
282 const G4Material* mat = track.GetMaterial();
283 if(mat != currentMat) {
284 currentMat = mat;
285 mfpKinEnergy = DBL_MAX;
286 matIdx = (G4int)track.GetMaterial()->GetIndex();
287 }
288 UpdateCrossSectionAndMFP(track.GetKineticEnergy());
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( G4UniformRand() );
300 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
301 } else {
302 theNumberOfInteractionLengthLeft -=
303 previousStepSize/currentInteractionLength;
304 theNumberOfInteractionLengthLeft =
305 std::max(theNumberOfInteractionLengthLeft, 0.0);
306 }
307 currentInteractionLength = theMFP;
308 return theNumberOfInteractionLengthLeft*theMFP;
309 }
310
311 G4double G4HadronicProcess::GetMeanFreePath(
312 const G4Track &aTrack, G4double,
313 G4ForceCondition*)
314 {
315 G4double xs = aScaleFactor*theCrossSectionDataStore
316 ->ComputeCrossSection(aTrack.GetDynamicParticle(),aTrack.GetMaterial());
317 return (xs > 0.0) ? 1.0/xs : DBL_MAX;
318 }
319
320 G4VParticleChange*
321 G4HadronicProcess::PostStepDoIt(const G4Track& aTrack, const G4Step&)
322 {
323 theNumberOfInteractionLengthLeft = -1.0;
324
325 //G4cout << "PostStepDoIt " << aTrack.GetDefinition()->GetParticleName()
326 // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
327 // if primary is not Alive then do nothing
328 theTotalResult->Clear();
329 theTotalResult->Initialize(aTrack);
330 fWeight = aTrack.GetWeight();
331 theTotalResult->ProposeWeight(fWeight);
332 if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
333
334 // Find cross section at end of step and check if <= 0
335 //
336 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
337 const G4Material* aMaterial = aTrack.GetMaterial();
338
339 // check only for charged particles
340 if(fXSType != fHadNoIntegral) {
341 mfpKinEnergy = DBL_MAX;
342 G4double xs = aScaleFactor*
343 theCrossSectionDataStore->ComputeCrossSection(aParticle,aMaterial);
344 //G4cout << "xs=" << xs << " xs0=" << theLastCrossSection
345 // << " " << aMaterial->GetName() << G4endl;
346 if(xs < theLastCrossSection*G4UniformRand()) {
347 // No interaction
348 return theTotalResult;
349 }
350 }
351
352 const G4Element* anElement =
353 theCrossSectionDataStore->SampleZandA(aParticle,aMaterial,targetNucleus);
354
355 // Next check for illegal track status
356 //
357 if (aTrack.GetTrackStatus() != fAlive &&
358 aTrack.GetTrackStatus() != fSuspend) {
359 if (aTrack.GetTrackStatus() == fStopAndKill ||
360 aTrack.GetTrackStatus() == fKillTrackAndSecondaries ||
361 aTrack.GetTrackStatus() == fPostponeToNextEvent) {
362 G4ExceptionDescription ed;
363 ed << "G4HadronicProcess: track in unusable state - "
364 << aTrack.GetTrackStatus() << G4endl;
365 ed << "G4HadronicProcess: returning unchanged track " << G4endl;
366 DumpState(aTrack,"PostStepDoIt",ed);
367 G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
368 }
369 // No warning for fStopButAlive which is a legal status here
370 return theTotalResult;
371 }
372
373 // Initialize the hadronic projectile from the track
374 thePro.Initialise(aTrack);
375
376 theInteraction = ChooseHadronicInteraction(thePro, targetNucleus,
377 aMaterial, anElement);
378 if(nullptr == theInteraction) {
379 G4ExceptionDescription ed;
380 ed << "Target element "<<anElement->GetName()<<" Z= "
381 << targetNucleus.GetZ_asInt() << " A= "
382 << targetNucleus.GetA_asInt() << G4endl;
383 DumpState(aTrack,"ChooseHadronicInteraction",ed);
384 ed << " No HadronicInteraction found out" << G4endl;
385 G4Exception("G4HadronicProcess::PostStepDoIt", "had005",
386 FatalException, ed);
387 return theTotalResult;
388 }
389
390 G4HadFinalState* result = nullptr;
391 G4int reentryCount = 0;
392 /*
393 G4cout << "### " << aParticle->GetDefinition()->GetParticleName()
394 << " Ekin(MeV)= " << aParticle->GetKineticEnergy()
395 << " Z= " << targetNucleus.GetZ_asInt()
396 << " A= " << targetNucleus.GetA_asInt()
397 << " by " << theInteraction->GetModelName()
398 << G4endl;
399 */
400 do
401 {
402 try
403 {
404 // Call the interaction
405 result = theInteraction->ApplyYourself( thePro, targetNucleus);
406 ++reentryCount;
407 }
408 catch(G4HadronicException & aR)
409 {
410 G4ExceptionDescription ed;
411 aR.Report(ed);
412 ed << "Call for " << theInteraction->GetModelName() << G4endl;
413 ed << "Target element "<<anElement->GetName()<<" Z= "
414 << targetNucleus.GetZ_asInt()
415 << " A= " << targetNucleus.GetA_asInt() << G4endl;
416 DumpState(aTrack,"ApplyYourself",ed);
417 ed << " ApplyYourself failed" << G4endl;
418 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
419 ed);
420 }
421
422 // Check the result for catastrophic energy non-conservation
423 result = CheckResult(thePro, targetNucleus, result);
424
425 if(reentryCount>100) {
426 G4ExceptionDescription ed;
427 ed << "Call for " << theInteraction->GetModelName() << G4endl;
428 ed << "Target element "<<anElement->GetName()<<" Z= "
429 << targetNucleus.GetZ_asInt()
430 << " A= " << targetNucleus.GetA_asInt() << G4endl;
431 DumpState(aTrack,"ApplyYourself",ed);
432 ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
433 G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
434 ed);
435 }
436 }
437 while(!result); /* Loop checking, 30-Oct-2015, G.Folger */
438
439 // Check whether kaon0 or anti_kaon0 are present between the secondaries:
440 // if this is the case, transform them into either kaon0S or kaon0L,
441 // with equal, 50% probability, keeping their dynamical masses (and
442 // the other kinematical properties).
443 // When this happens - very rarely - a "JustWarning" exception is thrown.
444 // Because Fluka-Cern produces kaon0 and anti_kaon0, we reduce the number
445 // of warnings to max 1 per thread.
446 G4int nSec = (G4int)result->GetNumberOfSecondaries();
447 if ( nSec > 0 ) {
448 for ( G4int i = 0; i < nSec; ++i ) {
449 auto dynamicParticle = result->GetSecondary(i)->GetParticle();
450 auto part = dynamicParticle->GetParticleDefinition();
451 if ( part == G4KaonZero::Definition() ||
452 part == G4AntiKaonZero::Definition() ) {
453 G4ParticleDefinition* newPart;
454 if ( G4UniformRand() > 0.5 ) { newPart = G4KaonZeroShort::Definition(); }
455 else { newPart = G4KaonZeroLong::Definition(); }
456 dynamicParticle->SetDefinition( newPart );
457 if ( nKaonWarn < 1 ) {
458 ++nKaonWarn;
459 G4ExceptionDescription ed;
460 ed << " Hadronic model " << theInteraction->GetModelName() << G4endl;
461 ed << " created " << part->GetParticleName() << G4endl;
462 ed << " -> forced to be " << newPart->GetParticleName() << G4endl;
463 G4Exception( "G4HadronicProcess::PostStepDoIt", "had007", JustWarning, ed );
464 }
465 }
466 }
467 }
468
469 result->SetTrafoToLab(thePro.GetTrafoToLab());
470 FillResult(result, aTrack);
471
472 if (epReportLevel != 0) {
473 CheckEnergyMomentumConservation(aTrack, targetNucleus);
474 }
475 //G4cout << "PostStepDoIt done nICelectrons= " << nICelectrons << G4endl;
476 return theTotalResult;
477 }
478
479 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const
480 {
481 outFile << "The description for this process has not been written yet.\n";
482 }
483
484 G4double G4HadronicProcess::XBiasSurvivalProbability()
485 {
486 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
487 G4double biasedProbability = 1.-G4Exp(-nLTraversed);
488 G4double realProbability = 1-G4Exp(-nLTraversed/aScaleFactor);
489 G4double result = (biasedProbability-realProbability)/biasedProbability;
490 return result;
491 }
492
493 G4double G4HadronicProcess::XBiasSecondaryWeight()
494 {
495 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
496 G4double result =
497 1./aScaleFactor*G4Exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
498 return result;
499 }
500
501 void
502 G4HadronicProcess::FillResult(G4HadFinalState * aR, const G4Track & aT)
503 {
504 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
505 const G4ThreeVector& dir = aT.GetMomentumDirection();
506
507 G4double efinal = std::max(aR->GetEnergyChange(), 0.0);
508
509 // check status of primary
510 if(aR->GetStatusChange() == stopAndKill) {
511 theTotalResult->ProposeTrackStatus(fStopAndKill);
512 theTotalResult->ProposeEnergy( 0.0 );
513
514 // check its final energy
515 } else if(0.0 == efinal) {
516 theTotalResult->ProposeEnergy( 0.0 );
517 if(aT.GetParticleDefinition()->GetProcessManager()
518 ->GetAtRestProcessVector()->size() > 0)
519 { theTotalResult->ProposeTrackStatus(fStopButAlive); }
520 else { theTotalResult->ProposeTrackStatus(fStopAndKill); }
521
522 // primary is not killed apply rotation and Lorentz transformation
523 } else {
524 theTotalResult->ProposeTrackStatus(fAlive);
525 G4ThreeVector newDir = aR->GetMomentumChange();
526 newDir.rotateUz(dir);
527 theTotalResult->ProposeMomentumDirection(newDir);
528 theTotalResult->ProposeEnergy(efinal);
529 }
530 //G4cout << "FillResult: Efinal= " << efinal << " status= "
531 // << theTotalResult->GetTrackStatus()
532 // << " fKill= " << fStopAndKill << G4endl;
533
534 // check secondaries
535 nICelectrons = 0;
536 G4int nSec = (G4int)aR->GetNumberOfSecondaries();
537 theTotalResult->SetNumberOfSecondaries(nSec);
538 G4double time0 = aT.GetGlobalTime();
539
540 for (G4int i = 0; i < nSec; ++i) {
541 G4DynamicParticle* dynParticle = aR->GetSecondary(i)->GetParticle();
542
543 // apply rotation
544 G4ThreeVector newDir = dynParticle->GetMomentumDirection();
545 newDir.rotateUz(dir);
546 dynParticle->SetMomentumDirection(newDir);
547
548 // check if secondary is on the mass shell
549 const G4ParticleDefinition* part = dynParticle->GetDefinition();
550 G4double mass = part->GetPDGMass();
551 G4double dmass= dynParticle->GetMass();
552 const G4double delta_mass_lim = 1.0*CLHEP::keV;
553 const G4double delta_ekin = 0.001*CLHEP::eV;
554 if(std::abs(dmass - mass) > delta_mass_lim) {
555 G4double e =
556 std::max(dynParticle->GetKineticEnergy() + dmass - mass, delta_ekin);
557 if(verboseLevel > 1) {
558 G4ExceptionDescription ed;
559 ed << "TrackID= "<< aT.GetTrackID()
560 << " " << aT.GetParticleDefinition()->GetParticleName()
561 << " Target Z= " << targetNucleus.GetZ_asInt() << " A= "
562 << targetNucleus.GetA_asInt()
563 << " Ekin(GeV)= " << aT.GetKineticEnergy()/CLHEP::GeV
564 << "\n Secondary is out of mass shell: " << part->GetParticleName()
565 << " EkinNew(MeV)= " << e
566 << " DeltaMass(MeV)= " << dmass - mass << G4endl;
567 G4Exception("G4HadronicProcess::FillResults", "had012", JustWarning, ed);
568 }
569 dynParticle->SetKineticEnergy(e);
570 dynParticle->SetMass(mass);
571 }
572 G4int idModel = aR->GetSecondary(i)->GetCreatorModelID();
573 if(part->GetPDGEncoding() == 11) { ++nICelectrons; }
574
575 // time of interaction starts from zero + global time
576 G4double time = std::max(aR->GetSecondary(i)->GetTime(), 0.0) + time0;
577
578 G4Track* track = new G4Track(dynParticle, time, aT.GetPosition());
579 track->SetCreatorModelID(idModel);
580 track->SetParentResonanceDef(aR->GetSecondary(i)->GetParentResonanceDef());
581 track->SetParentResonanceID(aR->GetSecondary(i)->GetParentResonanceID());
582 G4double newWeight = fWeight*aR->GetSecondary(i)->GetWeight();
583 track->SetWeight(newWeight);
584 track->SetTouchableHandle(aT.GetTouchableHandle());
585 theTotalResult->AddSecondary(track);
586 }
587 aR->Clear();
588 // G4cout << "FillResults done nICe= " << nICelectrons << G4endl;
589 }
590
591 void G4HadronicProcess::MultiplyCrossSectionBy(G4double factor)
592 {
593 BiasCrossSectionByFactor(factor);
594 }
595
596 void G4HadronicProcess::BiasCrossSectionByFactor(G4double aScale)
597 {
598 if (aScale <= 0.0) {
599 G4ExceptionDescription ed;
600 ed << " Wrong biasing factor " << aScale << " for " << GetProcessName();
601 G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010",
602 JustWarning, ed, "Cross-section bias is ignored");
603 } else {
604 aScaleFactor = aScale;
605 }
606 }
607
608 G4HadFinalState* G4HadronicProcess::CheckResult(const G4HadProjectile & aPro,
609 const G4Nucleus &aNucleus,
610 G4HadFinalState * result)
611 {
612 // check for catastrophic energy non-conservation
613 // to re-sample the interaction
614 G4HadronicInteraction * theModel = GetHadronicInteraction();
615 G4double nuclearMass(0);
616 if (nullptr != theModel) {
617
618 // Compute final-state total energy
619 G4double finalE(0.);
620 G4int nSec = (G4int)result->GetNumberOfSecondaries();
621
622 nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
623 aNucleus.GetZ_asInt());
624 if (result->GetStatusChange() != stopAndKill) {
625 // Interaction didn't complete, returned "do nothing" state
626 // and reset nucleus or the primary survived the interaction
627 // (e.g. electro-nuclear ) => keep nucleus
628 finalE=result->GetLocalEnergyDeposit() +
629 aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
630 if( nSec == 0 ){
631 // Since there are no secondaries, there is no recoil nucleus.
632 // To check energy balance we must neglect the initial nucleus too.
633 nuclearMass=0.0;
634 }
635 }
636 for (G4int i = 0; i < nSec; ++i) {
637 G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle();
638 finalE += pdyn->GetTotalEnergy();
639 G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass();
640 G4double mass_dyn=pdyn->GetMass();
641 if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV ) {
642 // If it is shortlived, then a difference less than 3 times the width is acceptable
643 if ( pdyn->GetDefinition()->IsShortLived() &&
644 std::abs(mass_pdg - mass_dyn) < 3.0*pdyn->GetDefinition()->GetPDGWidth() ) {
645 continue;
646 }
647 result->Clear();
648 result = nullptr;
649 G4ExceptionDescription desc;
650 desc << "Warning: Secondary with off-shell dynamic mass detected: "
651 << G4endl
652 << " " << pdyn->GetDefinition()->GetParticleName()
653 << ", PDG mass: " << mass_pdg << ", dynamic mass: "
654 << mass_dyn << G4endl
655 << (epReportLevel<0 ? "abort the event"
656 : "re-sample the interaction") << G4endl
657 << " Process / Model: " << GetProcessName()<< " / "
658 << theModel->GetModelName() << G4endl
659 << " Primary: " << aPro.GetDefinition()->GetParticleName()
660 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
661 << " E= " << aPro.Get4Momentum().e()
662 << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
663 << aNucleus.GetA_asInt() << ")" << G4endl;
664 G4Exception("G4HadronicProcess:CheckResult()", "had012",
665 epReportLevel<0 ? EventMustBeAborted : JustWarning,desc);
666 // must return here.....
667 return result;
668 }
669 }
670 G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
671
672 std::pair<G4double, G4double> checkLevels =
673 theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
674 if (std::abs(deltaE) > checkLevels.second &&
675 std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
676 // do not delete result, this is a pointer to a data member;
677 result->Clear();
678 result = nullptr;
679 G4ExceptionDescription desc;
680 desc << "Warning: Bad energy non-conservation detected, will "
681 << (epReportLevel<0 ? "abort the event"
682 : "re-sample the interaction") << G4endl
683 << " Process / Model: " << GetProcessName()<< " / "
684 << theModel->GetModelName() << G4endl
685 << " Primary: " << aPro.GetDefinition()->GetParticleName()
686 << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
687 << " E= " << aPro.Get4Momentum().e()
688 << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
689 << aNucleus.GetA_asInt() << ")" << G4endl
690 << " E(initial - final) = " << deltaE << " MeV." << G4endl;
691 G4Exception("G4HadronicProcess:CheckResult()", "had012",
692 epReportLevel<0 ? EventMustBeAborted : JustWarning,desc);
693 }
694 }
695 return result;
696 }
697
698 void
699 G4HadronicProcess::CheckEnergyMomentumConservation(const G4Track& aTrack,
700 const G4Nucleus& aNucleus)
701 {
702 G4int target_A=aNucleus.GetA_asInt();
703 G4int target_Z=aNucleus.GetZ_asInt();
704 G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
705 G4LorentzVector target4mom(0, 0, 0, targetMass
706 + nICelectrons*CLHEP::electron_mass_c2);
707
708 G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
709 G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
710 G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
711
712 G4int initial_A = target_A + track_A;
713 G4int initial_Z = target_Z + track_Z - nICelectrons;
714
715 G4LorentzVector initial4mom = projectile4mom + target4mom;
716
717 // Compute final-state momentum for scattering and "do nothing" results
718 G4LorentzVector final4mom;
719 G4int final_A(0), final_Z(0);
720
721 G4int nSec = theTotalResult->GetNumberOfSecondaries();
722 if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
723 // Either interaction didn't complete, returned "do nothing" state
724 // or the primary survived the interaction (e.g. electro-nucleus )
725
726 // Interaction didn't complete, returned "do nothing" state
727 // - or suppressed recoil (e.g. Neutron elastic )
728 final4mom = initial4mom;
729 final_A = initial_A;
730 final_Z = initial_Z;
731 if (nSec > 0) {
732 // The primary remains in final state (e.g. electro-nucleus )
733 // Use the final energy / momentum
734 const G4ThreeVector& v = *theTotalResult->GetMomentumDirection();
735 G4double ekin = theTotalResult->GetEnergy();
736 G4double mass = aTrack.GetDefinition()->GetPDGMass();
737 G4double ptot = std::sqrt(ekin*(ekin + 2*mass));
738 final4mom.set(ptot*v.x(), ptot*v.y(), ptot*v.z(), mass + ekin);
739 final_A = track_A;
740 final_Z = track_Z;
741 // Expect that the target nucleus will have interacted,
742 // and its products, including recoil, will be included in secondaries.
743 }
744 }
745 if( nSec > 0 ) {
746 G4Track* sec;
747
748 for (G4int i = 0; i < nSec; i++) {
749 sec = theTotalResult->GetSecondary(i);
750 final4mom += sec->GetDynamicParticle()->Get4Momentum();
751 final_A += sec->GetDefinition()->GetBaryonNumber();
752 final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
753 }
754 }
755
756 // Get level-checking information (used to cut-off relative checks)
757 G4String processName = GetProcessName();
758 G4HadronicInteraction* theModel = GetHadronicInteraction();
759 G4String modelName("none");
760 if (theModel) modelName = theModel->GetModelName();
761 std::pair<G4double, G4double> checkLevels = epCheckLevels;
762 if (!levelsSetByProcess) {
763 if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
764 checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
765 checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
766 }
767
768 // Compute absolute total-energy difference, and relative kinetic-energy
769 G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
770
771 G4LorentzVector diff = initial4mom - final4mom;
772 G4double absolute = diff.e();
773 G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
774
775 G4double absolute_mom = diff.vect().mag();
776 G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
777
778 // Evaluate relative and absolute conservation
779 G4bool relPass = true;
780 G4String relResult = "pass";
781 if ( std::abs(relative) > checkLevels.first
782 || std::abs(relative_mom) > checkLevels.first) {
783 relPass = false;
784 relResult = checkRelative ? "fail" : "N/A";
785 }
786
787 G4bool absPass = true;
788 G4String absResult = "pass";
789 if ( std::abs(absolute) > checkLevels.second
790 || std::abs(absolute_mom) > checkLevels.second ) {
791 absPass = false ;
792 absResult = "fail";
793 }
794
795 G4bool chargePass = true;
796 G4String chargeResult = "pass";
797 if ( (initial_A-final_A)!=0
798 || (initial_Z-final_Z)!=0 ) {
799 chargePass = checkLevels.second < DBL_MAX ? false : true;
800 chargeResult = "fail";
801 }
802
803 G4bool conservationPass = (relPass || absPass) && chargePass;
804
805 std::stringstream Myout;
806 G4bool Myout_notempty(false);
807 // Options for level of reporting detail:
808 // 0. off
809 // 1. report only when E/p not conserved
810 // 2. report regardless of E/p conservation
811 // 3. report only when E/p not conserved, with model names, process names, and limits
812 // 4. report regardless of E/p conservation, with model names, process names, and limits
813 // negative -1.., as above, but send output to stderr
814
815 if( std::abs(epReportLevel) == 4
816 || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
817 Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
818 Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
819 << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
820 << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
821 << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
822 << aNucleus.GetA_asInt() << ")" << G4endl;
823 Myout_notempty=true;
824 }
825 if ( std::abs(epReportLevel) == 4
826 || std::abs(epReportLevel) == 2
827 || ! conservationPass ){
828
829 Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
830 << relative << " p/p(0)= " << relative_mom << G4endl;
831 Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
832 << absolute/MeV << " / " << absolute_mom/MeV << " 3mom: " << (diff.vect())*1./MeV << G4endl;
833 Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
834 Myout_notempty=true;
835
836 }
837 Myout.flush();
838 if ( Myout_notempty ) {
839 if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
840 else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
841 }
842 }
843
844 void G4HadronicProcess::DumpState(const G4Track& aTrack,
845 const G4String& method,
846 G4ExceptionDescription& ed)
847 {
848 ed << "Unrecoverable error in the method " << method << " of "
849 << GetProcessName() << G4endl;
850 ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
851 << aTrack.GetParentID()
852 << " " << aTrack.GetParticleDefinition()->GetParticleName()
853 << G4endl;
854 ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
855 << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
856 ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
857
858 if (aTrack.GetMaterial()) {
859 ed << " material " << aTrack.GetMaterial()->GetName();
860 }
861 ed << G4endl;
862
863 if (aTrack.GetVolume()) {
864 ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
865 << ">" << G4endl;
866 }
867 }
868
869 void G4HadronicProcess::DumpPhysicsTable(const G4ParticleDefinition& p)
870 {
871 theCrossSectionDataStore->DumpPhysicsTable(p);
872 }
873
874 void G4HadronicProcess::AddDataSet(G4VCrossSectionDataSet * aDataSet)
875 {
876 theCrossSectionDataStore->AddDataSet(aDataSet);
877 }
878
879 std::vector<G4HadronicInteraction*>&
880 G4HadronicProcess::GetHadronicInteractionList()
881 {
882 return theEnergyRangeManager.GetHadronicInteractionList();
883 }
884
885 G4HadronicInteraction*
886 G4HadronicProcess::GetHadronicModel(const G4String& modelName)
887 {
888 std::vector<G4HadronicInteraction*>& list
889 = theEnergyRangeManager.GetHadronicInteractionList();
890 for (auto & mod : list) {
891 if (mod->GetModelName() == modelName) return mod;
892 }
893 return nullptr;
894 }
895
896 G4double
897 G4HadronicProcess::ComputeCrossSection(const G4ParticleDefinition* part,
898 const G4Material* mat,
899 const G4double kinEnergy)
900 {
901 auto dp = new G4DynamicParticle(part, unitVector, kinEnergy);
902 G4double xs = theCrossSectionDataStore->ComputeCrossSection(dp, mat);
903 delete dp;
904 return xs;
905 }
906
907 void G4HadronicProcess::RecomputeXSandMFP(const G4double kinEnergy)
908 {
909 auto dp = new G4DynamicParticle(currentParticle, unitVector, kinEnergy);
910 theLastCrossSection = aScaleFactor*
911 theCrossSectionDataStore->ComputeCrossSection(dp, currentMat);
912 theMFP = (theLastCrossSection > 0.0) ? 1.0/theLastCrossSection : DBL_MAX;
913 delete dp;
914 }
915
916 void G4HadronicProcess::UpdateCrossSectionAndMFP(const G4double e)
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 > minKinEnergy) {
929 G4double e1 = std::max(e*lambdaFactor, minKinEnergy);
930 mfpKinEnergy = e1;
931 RecomputeXSandMFP(e1);
932 }
933
934 } else if(fXSType == fHadOnePeak) {
935 G4double epeak = (*theEnergyOfCrossSectionMax)[matIdx];
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*lambdaFactor);
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 <= mfpKinEnergy) {
963 const G4double e1 = std::max(e1peak, e*lambdaFactor);
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 <= mfpKinEnergy) {
982 const G4double e1 = std::max(e2peak, e*lambdaFactor);
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*lambdaFactor);
1000 mfpKinEnergy = e1;
1001 RecomputeXSandMFP(e1);
1002 }
1003
1004 } else {
1005 DefineXSandMFP();
1006 }
1007 }
1008