ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronicProcess.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4HadronicProcess.cc
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"
67 
68 #include "G4NistManager.hh"
69 #include "G4PhysicsModelCatalog.hh"
71 #include "G4Exp.hh"
72 
73 #include <typeinfo>
74 #include <sstream>
75 #include <iostream>
76 
77 // File-scope variable to capture environment variable at startup
78 
79 static const char* G4Hadronic_Random_File = std::getenv("G4HADRONIC_RANDOM_FILE");
80 
82 
84  G4ProcessType procType)
85  : G4VDiscreteProcess(processName, procType)
86 {
87  SetProcessSubType(fHadronInelastic); // Default unless subclass changes
89 }
90 
92 
94  G4HadronicProcessType aHadSubType)
95  : G4VDiscreteProcess(processName, fHadronic)
96 {
97  SetProcessSubType(aHadSubType);
99 }
100 
102 {
104  delete theTotalResult;
106 }
107 
111  theInteraction = nullptr;
114  theProcessStore->Register(this);
116  aScaleFactor = 1.0;
117  fWeight = 1.0;
118  xBiasOn = false;
119  nMatWarn = 0;
120  useIntegralXS = true;
121  theLastCrossSection = 0.0;
122  nICelectrons = 0;
123  idxIC = -1;
126 }
127 
129  levelsSetByProcess = false;
130 
131  epReportLevel = std::getenv("G4Hadronic_epReportLevel") ?
132  std::strtol(std::getenv("G4Hadronic_epReportLevel"),0,10) : 0;
133 
134  epCheckLevels.first = std::getenv("G4Hadronic_epCheckRelativeLevel") ?
135  std::strtod(std::getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX;
136 
137  epCheckLevels.second = std::getenv("G4Hadronic_epCheckAbsoluteLevel") ?
138  std::strtod(std::getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX;
139 }
140 
142 {
143  if(!a) { return; }
144  try{ theEnergyRangeManager.RegisterMe( a ); }
145  catch(G4HadronicException & aE)
146  {
148  aE.Report(ed);
149  ed << "Unrecoverable error in " << GetProcessName()
150  << " to register " << a->GetModelName() << G4endl;
151  G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException,
152  ed);
153  }
155 }
156 
157 G4double
159  const G4Element * elm,
160  const G4Material* mat)
161 {
162  if(!mat)
163  {
164  ++nMatWarn;
165  static const G4int nmax = 5;
166  if(nMatWarn < nmax) {
168  ed << "Cannot compute Element x-section for " << GetProcessName()
169  << " because no material defined \n"
170  << " Please, specify material pointer or define simple material"
171  << " for Z= " << elm->GetZasInt();
172  G4Exception("G4HadronicProcess::GetElementCrossSection", "had066",
173  JustWarning, ed);
174  }
175  }
176  return
177  std::max(theCrossSectionDataStore->GetCrossSection(part, elm, mat),0.0);
178 }
179 
181 {
182  if(std::getenv("G4HadronicProcess_debug")) {
184  }
186 }
187 
189 {
193 }
194 
197 {
198  //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName()
199  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
203  //G4cout << " xsection= " << theLastCrossSection << G4endl;
204  return res;
205 }
206 
209 {
210  //G4cout << "PostStepDoIt " << aTrack.GetDefinition()->GetParticleName()
211  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
212  // if primary is not Alive then do nothing
214  theTotalResult->Initialize(aTrack);
215  fWeight = aTrack.GetWeight();
217  if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
218 
219  // Find cross section at end of step and check if <= 0
220  //
221  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
222  const G4Material* aMaterial = aTrack.GetMaterial();
223 
224  // check only for charged particles
225  if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
226  G4double xs = aScaleFactor*
227  theCrossSectionDataStore->ComputeCrossSection(aParticle,aMaterial);
228  if(xs <= 0.0 || xs < theLastCrossSection*G4UniformRand()) {
229  // No interaction
230  return theTotalResult;
231  }
232  }
233 
234  const G4Element* anElement =
235  theCrossSectionDataStore->SampleZandA(aParticle,aMaterial,targetNucleus);
236 
237  // Next check for illegal track status
238  //
239  if (aTrack.GetTrackStatus() != fAlive &&
240  aTrack.GetTrackStatus() != fSuspend) {
241  if (aTrack.GetTrackStatus() == fStopAndKill ||
243  aTrack.GetTrackStatus() == fPostponeToNextEvent) {
245  ed << "G4HadronicProcess: track in unusable state - "
246  << aTrack.GetTrackStatus() << G4endl;
247  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
248  DumpState(aTrack,"PostStepDoIt",ed);
249  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
250  }
251  // No warning for fStopButAlive which is a legal status here
252  return theTotalResult;
253  }
254 
255  // Initialize the hadronic projectile from the track
256  thePro.Initialise(aTrack);
257 
259  aMaterial, anElement);
260  if(!theInteraction) {
262  ed << "Target element "<<anElement->GetName()<<" Z= "
263  << targetNucleus.GetZ_asInt() << " A= "
265  DumpState(aTrack,"ChooseHadronicInteraction",ed);
266  ed << " No HadronicInteraction found out" << G4endl;
267  G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException, ed);
268  return theTotalResult;
269  }
270 
271  G4HadFinalState* result = nullptr;
272  G4int reentryCount = 0;
273  /*
274  G4cout << "### " << aParticle->GetDefinition()->GetParticleName()
275  << " Ekin(MeV)= " << aParticle->GetKineticEnergy()
276  << " Z= " << targetNucleus.GetZ_asInt()
277  << " A= " << targetNucleus.GetA_asInt()
278  << " by " << theInteraction->GetModelName()
279  << G4endl;
280  */
281  do
282  {
283  try
284  {
285  // Save random engine if requested for debugging
288  }
289  // Call the interaction
291  ++reentryCount;
292  }
293  catch(G4HadronicException & aR)
294  {
296  aR.Report(ed);
297  ed << "Call for " << theInteraction->GetModelName() << G4endl;
298  ed << "Target element "<<anElement->GetName()<<" Z= "
300  << " A= " << targetNucleus.GetA_asInt() << G4endl;
301  DumpState(aTrack,"ApplyYourself",ed);
302  ed << " ApplyYourself failed" << G4endl;
303  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
304  ed);
305  }
306 
307  // Check the result for catastrophic energy non-conservation
308  result = CheckResult(thePro, targetNucleus, result);
309 
310  if(reentryCount>100) {
312  ed << "Call for " << theInteraction->GetModelName() << G4endl;
313  ed << "Target element "<<anElement->GetName()<<" Z= "
315  << " A= " << targetNucleus.GetA_asInt() << G4endl;
316  DumpState(aTrack,"ApplyYourself",ed);
317  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
318  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
319  ed);
320  }
321  }
322  while(!result); /* Loop checking, 30-Oct-2015, G.Folger */
323 
324  // Check whether kaon0 or anti_kaon0 are present between the secondaries:
325  // if this is the case, transform them into either kaon0S or kaon0L,
326  // with equal, 50% probability, keeping their dynamical masses (and
327  // the other kinematical properties).
328  // When this happens - very rarely - a "JustWarning" exception is thrown.
329  G4int nSec = result->GetNumberOfSecondaries();
330  if ( nSec > 0 ) {
331  for ( G4int i = 0; i < nSec; ++i ) {
332  G4DynamicParticle* dynamicParticle = result->GetSecondary(i)->GetParticle();
333  const G4ParticleDefinition* particleDefinition =
334  dynamicParticle->GetParticleDefinition();
335  if ( particleDefinition == G4KaonZero::Definition() ||
336  particleDefinition == G4AntiKaonZero::Definition() ) {
337  G4ParticleDefinition* newPart;
338  if( G4UniformRand() > 0.5 ) { newPart = G4KaonZeroShort::Definition(); }
339  else { newPart = G4KaonZeroLong::Definition(); }
340  dynamicParticle->SetDefinition( newPart );
342  ed << " Hadronic model " << theInteraction->GetModelName() << G4endl;
343  ed << " created " << particleDefinition->GetParticleName() << G4endl;
344  ed << " -> forced to be " << newPart->GetParticleName() << G4endl;
345  G4Exception( "G4HadronicProcess::PostStepDoIt", "had007", JustWarning, ed );
346  }
347  }
348  }
349 
351 
353 
354  FillResult(result, aTrack);
355 
356  if (epReportLevel != 0) {
358  }
359  //G4cout << "PostStepDoIt done nICelectrons= " << nICelectrons << G4endl;
360  return theTotalResult;
361 }
362 
363 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const
364 {
365  outFile << "The description for this process has not been written yet.\n";
366 }
367 
368 
370 {
372  G4double biasedProbability = 1.-G4Exp(-nLTraversed);
373  G4double realProbability = 1-G4Exp(-nLTraversed/aScaleFactor);
374  G4double result = (biasedProbability-realProbability)/biasedProbability;
375  return result;
376 }
377 
379 {
381  G4double result =
382  1./aScaleFactor*G4Exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
383  return result;
384 }
385 
386 void
388 {
390  const G4ThreeVector& dir = aT.GetMomentumDirection();
391 
392  G4double efinal = std::max(aR->GetEnergyChange(), 0.0);
393 
394  // check status of primary
395  if(aR->GetStatusChange() == stopAndKill) {
398 
399  // check its final energy
400  } else if(0.0 == efinal) {
403  ->GetAtRestProcessVector()->size() > 0)
406 
407  // primary is not killed apply rotation and Lorentz transformation
408  } else {
410  G4ThreeVector newDir = aR->GetMomentumChange();
411  newDir.rotateUz(dir);
413  theTotalResult->ProposeEnergy(efinal);
414  }
415  //G4cout << "FillResult: Efinal= " << efinal << " status= "
416  // << theTotalResult->GetTrackStatus()
417  // << " fKill= " << fStopAndKill << G4endl;
418 
419  // check secondaries
420  nICelectrons = 0;
421  if(idxIC == -1) {
422  G4int idx = G4PhysicsModelCatalog::GetIndex("e-InternalConvertion");
423  idxIC = -1 == idx ? -2 : idx;
424  }
425  G4int nSec = aR->GetNumberOfSecondaries();
427  G4double time0 = aT.GetGlobalTime();
428 
429  for (G4int i = 0; i < nSec; ++i) {
430  G4DynamicParticle* dynParticle = aR->GetSecondary(i)->GetParticle();
431 
432  // apply rotation
433  G4ThreeVector newDir = dynParticle->GetMomentumDirection();
434  newDir.rotateUz(dir);
435  dynParticle->SetMomentumDirection(newDir);
436 
437  // check if secondary is on the mass shell
438  const G4ParticleDefinition* part = dynParticle->GetDefinition();
439  G4double mass = part->GetPDGMass();
440  G4double dmass= dynParticle->GetMass();
441  const G4double delta_mass_lim = 1.0*CLHEP::keV;
442  const G4double delta_ekin = 0.001*CLHEP::eV;
443  if(std::abs(dmass - mass) > delta_mass_lim) {
444  G4double e = std::max(dynParticle->GetKineticEnergy() + dmass - mass, delta_ekin);
447  ed << "TrackID= "<< aT.GetTrackID()
448  << " " << aT.GetParticleDefinition()->GetParticleName()
449  << " Target Z= " << targetNucleus.GetZ_asInt() << " A= "
451  << " Ekin(GeV)= " << aT.GetKineticEnergy()/CLHEP::GeV
452  << "\n Secondary is out of mass shell: " << part->GetParticleName()
453  << " EkinNew(MeV)= " << e
454  << " DeltaMass(MeV)= " << dmass - mass << G4endl;
455  G4Exception("G4HadronicProcess::FillResults", "had012", JustWarning, ed);
456  }
457  dynParticle->SetKineticEnergy(e);
458  dynParticle->SetMass(mass);
459  }
460  G4int idxModel = aR->GetSecondary(i)->GetCreatorModelType();
461  //if(idxIC == idxModel) { ++nICelectrons; }
462  if(part->GetPDGEncoding() == 11) { ++nICelectrons; }
463 
464  // time of interaction starts from zero + global time
465  G4double time = std::max(aR->GetSecondary(i)->GetTime(), 0.0) + time0;
466 
467  G4Track* track = new G4Track(dynParticle, time, aT.GetPosition());
468  track->SetCreatorModelIndex(idxModel);
469  G4double newWeight = fWeight*aR->GetSecondary(i)->GetWeight();
470  track->SetWeight(newWeight);
474  G4double e = dynParticle->GetKineticEnergy();
475  if (e == 0.0) {
477  DumpState(aT,"Secondary has zero energy",ed);
478  ed << "Secondary " << part->GetParticleName()
479  << G4endl;
480  G4Exception("G4HadronicProcess::FillResults", "had011",
481  JustWarning,ed);
482  }
483  }
484  }
485  aR->Clear();
486  // G4cout << "FillResults done nICe= " << nICelectrons << G4endl;
487 }
488 
490 {
491  BiasCrossSectionByFactor(factor);
492 }
493 
495 {
496  if (aScale <= 0.0) {
498  ed << " Wrong biasing factor " << aScale << " for " << GetProcessName();
499  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010",
500  JustWarning, ed, "Cross-section bias is ignored");
501  } else {
502  xBiasOn = true;
503  aScaleFactor = aScale;
504  }
505 }
506 
508  const G4Nucleus &aNucleus,
509  G4HadFinalState * result)
510 {
511  // check for catastrophic energy non-conservation
512  // to re-sample the interaction
513 
515  G4double nuclearMass(0);
516  if (theModel) {
517 
518  // Compute final-state total energy
519  G4double finalE(0.);
520  G4int nSec = result->GetNumberOfSecondaries();
521 
522  nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
523  aNucleus.GetZ_asInt());
524  if (result->GetStatusChange() != stopAndKill) {
525  // Interaction didn't complete, returned "do nothing" state
526  // and reset nucleus or the primary survived the interaction
527  // (e.g. electro-nuclear ) => keep nucleus
528  finalE=result->GetLocalEnergyDeposit() +
529  aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
530  if( nSec == 0 ){
531  // Since there are no secondaries, there is no recoil nucleus.
532  // To check energy balance we must neglect the initial nucleus too.
533  nuclearMass=0.0;
534  }
535  }
536  for (G4int i = 0; i < nSec; i++) {
537  G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle();
538  finalE += pdyn->GetTotalEnergy();
539  G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass();
540  G4double mass_dyn=pdyn->GetMass();
541  if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV ) {
542  // If it is shortlived, then a difference less than 3 times the width is acceptable
543  if ( pdyn->GetDefinition()->IsShortLived() &&
544  std::abs(mass_pdg - mass_dyn) < 3.0*pdyn->GetDefinition()->GetPDGWidth() ) {
545  continue;
546  }
547  result->Clear();
548  result = nullptr;
550  desc << "Warning: Secondary with off-shell dynamic mass detected: "
551  << G4endl
552  << " " << pdyn->GetDefinition()->GetParticleName()
553  << ", PDG mass: " << mass_pdg << ", dynamic mass: "
554  << mass_dyn << G4endl
555  << (epReportLevel<0 ? "abort the event"
556  : "re-sample the interaction") << G4endl
557  << " Process / Model: " << GetProcessName()<< " / "
558  << theModel->GetModelName() << G4endl
559  << " Primary: " << aPro.GetDefinition()->GetParticleName()
560  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
561  << " E= " << aPro.Get4Momentum().e()
562  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
563  << aNucleus.GetA_asInt() << ")" << G4endl;
564  G4Exception("G4HadronicProcess:CheckResult()", "had012",
566  // must return here.....
567  return result;
568  }
569  }
570  G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
571 
572  std::pair<G4double, G4double> checkLevels =
573  theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
574  if (std::abs(deltaE) > checkLevels.second &&
575  std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
576  // do not delete result, this is a pointer to a data member;
577  result->Clear();
578  result = nullptr;
580  desc << "Warning: Bad energy non-conservation detected, will "
581  << (epReportLevel<0 ? "abort the event"
582  : "re-sample the interaction") << G4endl
583  << " Process / Model: " << GetProcessName()<< " / "
584  << theModel->GetModelName() << G4endl
585  << " Primary: " << aPro.GetDefinition()->GetParticleName()
586  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
587  << " E= " << aPro.Get4Momentum().e()
588  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
589  << aNucleus.GetA_asInt() << ")" << G4endl
590  << " E(initial - final) = " << deltaE << " MeV." << G4endl;
591  G4Exception("G4HadronicProcess:CheckResult()", "had012",
593  }
594  }
595  return result;
596 }
597 
598 void
600  const G4Nucleus& aNucleus)
601 {
602  G4int target_A=aNucleus.GetA_asInt();
603  G4int target_Z=aNucleus.GetZ_asInt();
604  G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
605  G4LorentzVector target4mom(0, 0, 0, targetMass
607 
608  G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
609  G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
610  G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
611 
612  G4int initial_A = target_A + track_A;
613  G4int initial_Z = target_Z + track_Z - nICelectrons;
614 
615  G4LorentzVector initial4mom = projectile4mom + target4mom;
616 
617  // Compute final-state momentum for scattering and "do nothing" results
618  G4LorentzVector final4mom;
619  G4int final_A(0), final_Z(0);
620 
622  if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
623  // Either interaction didn't complete, returned "do nothing" state
624  // or the primary survived the interaction (e.g. electro-nucleus )
625 
626  // Interaction didn't complete, returned "do nothing" state
627  // - or suppressed recoil (e.g. Neutron elastic )
628  final4mom = initial4mom;
629  final_A = initial_A;
630  final_Z = initial_Z;
631  if (nSec > 0 && aTrack.GetDynamicParticle()) {
632  // The primary remains in final state (e.g. electro-nucleus )
633  G4Track temp(aTrack);
634 
635  // Use the final energy / momentum
638  final4mom = temp.GetDynamicParticle()->Get4Momentum();
639  final_A = track_A;
640  final_Z = track_Z;
641  // Expect that the target nucleus will have interacted,
642  // and its products, including recoil, will be included in secondaries.
643  }
644  }
645  if( nSec > 0 ) {
646  G4Track* sec;
647 
648  for (G4int i = 0; i < nSec; i++) {
649  sec = theTotalResult->GetSecondary(i);
650  final4mom += sec->GetDynamicParticle()->Get4Momentum();
651  final_A += sec->GetDefinition()->GetBaryonNumber();
652  final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
653  }
654  }
655 
656  // Get level-checking information (used to cut-off relative checks)
657  G4String processName = GetProcessName();
659  G4String modelName("none");
660  if (theModel) modelName = theModel->GetModelName();
661  std::pair<G4double, G4double> checkLevels = epCheckLevels;
662  if (!levelsSetByProcess) {
663  if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
664  checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
665  checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
666  }
667 
668  // Compute absolute total-energy difference, and relative kinetic-energy
669  G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
670 
671  G4LorentzVector diff = initial4mom - final4mom;
672  G4double absolute = diff.e();
673  G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
674 
675  G4double absolute_mom = diff.vect().mag();
676  G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
677 
678  // Evaluate relative and absolute conservation
679  G4bool relPass = true;
680  G4String relResult = "pass";
681  if ( std::abs(relative) > checkLevels.first
682  || std::abs(relative_mom) > checkLevels.first) {
683  relPass = false;
684  relResult = checkRelative ? "fail" : "N/A";
685  }
686 
687  G4bool absPass = true;
688  G4String absResult = "pass";
689  if ( std::abs(absolute) > checkLevels.second
690  || std::abs(absolute_mom) > checkLevels.second ) {
691  absPass = false ;
692  absResult = "fail";
693  }
694 
695  G4bool chargePass = true;
696  G4String chargeResult = "pass";
697  if ( (initial_A-final_A)!=0
698  || (initial_Z-final_Z)!=0 ) {
699  chargePass = checkLevels.second < DBL_MAX ? false : true;
700  chargeResult = "fail";
701  }
702 
703  G4bool conservationPass = (relPass || absPass) && chargePass;
704 
705  std::stringstream Myout;
706  G4bool Myout_notempty(false);
707  // Options for level of reporting detail:
708  // 0. off
709  // 1. report only when E/p not conserved
710  // 2. report regardless of E/p conservation
711  // 3. report only when E/p not conserved, with model names, process names, and limits
712  // 4. report regardless of E/p conservation, with model names, process names, and limits
713  // negative -1.., as above, but send output to stderr
714 
715  if( std::abs(epReportLevel) == 4
716  || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
717  Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
718  Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
719  << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
720  << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
721  << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
722  << aNucleus.GetA_asInt() << ")" << G4endl;
723  Myout_notempty=true;
724  }
725  if ( std::abs(epReportLevel) == 4
726  || std::abs(epReportLevel) == 2
727  || ! conservationPass ){
728 
729  Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
730  << relative << " p/p(0)= " << relative_mom << G4endl;
731  Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
732  << absolute/MeV << " / " << absolute_mom/MeV << " 3mom: " << (diff.vect())*1./MeV << G4endl;
733  Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
734  Myout_notempty=true;
735 
736  }
737  Myout.flush();
738  if ( Myout_notempty ) {
739  if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
740  else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
741  }
742 }
743 
745  const G4String& method,
747 {
748  ed << "Unrecoverable error in the method " << method << " of "
749  << GetProcessName() << G4endl;
750  ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
751  << aTrack.GetParentID()
752  << " " << aTrack.GetParticleDefinition()->GetParticleName()
753  << G4endl;
754  ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
755  << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
756  ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
757 
758  if (aTrack.GetMaterial()) {
759  ed << " material " << aTrack.GetMaterial()->GetName();
760  }
761  ed << G4endl;
762 
763  if (aTrack.GetVolume()) {
764  ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
765  << ">" << G4endl;
766  }
767 }
768 
770 {
772 }
773 
775 {
777 }
778 
779 std::vector<G4HadronicInteraction*>&
781 {
783 }
784 
787 {
788  std::vector<G4HadronicInteraction*>& list
790  for (size_t li=0; li<list.size(); li++) {
791  if (list[li]->GetModelName() == modelName) return list[li];
792  }
793  return nullptr;
794 }