Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4HadronicProcess.cc 110727 2018-06-11 06:08:11Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class source file
31 //
32 // G4HadronicProcess
33 //
34 // original by H.P.Wellisch
35 // J.L. Chuma, TRIUMF, 10-Mar-1997
36 //
37 // Modifications:
38 // 05-Jul-2010 V.Ivanchenko cleanup commented lines
39 // 20-Jul-2011 M.Kelsey -- null-pointer checks in DumpState()
40 // 24-Sep-2011 M.Kelsey -- Use envvar G4HADRONIC_RANDOM_FILE to save random
41 // engine state before each model call
42 // 18-Oct-2011 M.Kelsey -- Handle final-state cases in conservation checks.
43 // 14-Mar-2012 G.Folger -- enhance checks for conservation of energy, etc.
44 // 28-Jul-2012 M.Maire -- add function GetTargetDefinition()
45 // 14-Sep-2012 Inherit from RestDiscrete, use subtype code (now in ctor) to
46 // configure base-class
47 // 28-Sep-2012 Restore inheritance from G4VDiscreteProcess, remove enable-flag
48 // changing, remove warning message from original ctor.
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"
61 #include "G4Navigator.hh"
62 #include "G4ProcessVector.hh"
63 #include "G4ProcessManager.hh"
64 #include "G4StableIsotopes.hh"
65 #include "G4HadTmpUtil.hh"
66 #include "G4NucleiProperties.hh"
67 
68 #include "G4HadronicException.hh"
71 
72 #include "G4AutoLock.hh"
73 #include "G4NistManager.hh"
74 #include "G4PhysicsModelCatalog.hh"
76 #include "G4Exp.hh"
77 
78 #include <typeinfo>
79 #include <sstream>
80 #include <iostream>
81 
82 #include <stdlib.h>
83 
84 // File-scope variable to capture environment variable at startup
85 
86 static const char* G4Hadronic_Random_File = getenv("G4HADRONIC_RANDOM_FILE");
87 
89 
91  G4ProcessType procType)
92  : G4VDiscreteProcess(processName, procType)
93 {
94  SetProcessSubType(fHadronInelastic); // Default unless subclass changes
96 }
97 
99 
101  G4HadronicProcessType aHadSubType)
102  : G4VDiscreteProcess(processName, fHadronic)
103 {
104  SetProcessSubType(aHadSubType);
105  InitialiseLocal();
106 }
107 
109 {
111  delete theTotalResult;
113 }
114 
118  theInteraction = nullptr;
121  theProcessStore->Register(this);
123  aScaleFactor = 1.0;
124  fWeight = 1.0;
125  xBiasOn = false;
126  nMatWarn = 0;
127  useIntegralXS = true;
128  theLastCrossSection = 0.0;
129  nICelectrons = 0;
130  idxIC = -1;
133 }
134 
136  levelsSetByProcess = false;
137 
138  epReportLevel = getenv("G4Hadronic_epReportLevel") ?
139  strtol(getenv("G4Hadronic_epReportLevel"),0,10) : 0;
140 
141  epCheckLevels.first = getenv("G4Hadronic_epCheckRelativeLevel") ?
142  strtod(getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX;
143 
144  epCheckLevels.second = getenv("G4Hadronic_epCheckAbsoluteLevel") ?
145  strtod(getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX;
146 }
147 
149 {
150  if(!a) { return; }
151  try{ theEnergyRangeManager.RegisterMe( a ); }
152  catch(G4HadronicException & aE)
153  {
155  aE.Report(ed);
156  ed << "Unrecoverable error in " << GetProcessName()
157  << " to register " << a->GetModelName() << G4endl;
158  G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException,
159  ed);
160  }
162 }
163 
164 G4double
166  const G4Element * elm,
167  const G4Material* mat)
168 {
169  if(!mat)
170  {
171  ++nMatWarn;
172  static const G4int nmax = 5;
173  if(nMatWarn < nmax) {
175  ed << "Cannot compute Element x-section for " << GetProcessName()
176  << " because no material defined \n"
177  << " Please, specify material pointer or define simple material"
178  << " for Z= " << elm->GetZasInt();
179  G4Exception("G4HadronicProcess::GetElementCrossSection", "had066",
180  JustWarning, ed);
181  }
182  }
183  return
184  std::max(theCrossSectionDataStore->GetCrossSection(part, elm, mat),0.0);
185 }
186 
188 {
189  if(getenv("G4HadronicProcess_debug")) {
191  }
193 }
194 
196 {
197  try
198  {
201  }
202  catch(G4HadronicException & aR)
203  {
205  aR.Report(ed);
206  ed << " hadronic initialisation fails";
207  G4Exception("G4HadronicProcess::BuildPhysicsTable", "had000",
208  FatalException,ed);
209  }
211 }
212 
215 {
216  //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName()
217  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
218  try
219  {
222  aTrack.GetMaterial());
223  }
224  catch(G4HadronicException & aR)
225  {
227  aR.Report(ed);
228  DumpState(aTrack,"GetMeanFreePath",ed);
229  ed << " Cross section is not available" << G4endl;
230  G4Exception("G4HadronicProcess::GetMeanFreePath", "had002", FatalException,
231  ed);
232  }
234  //G4cout << " xsection= " << theLastCrossSection << G4endl;
235  return res;
236 }
237 
240 {
241  //G4cout << "PostStepDoIt " << aTrack.GetDefinition()->GetParticleName()
242  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
243  // if primary is not Alive then do nothing
245  theTotalResult->Initialize(aTrack);
246  fWeight = aTrack.GetWeight();
248  if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
249 
250  // Find cross section at end of step and check if <= 0
251  //
252  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
253  const G4Material* aMaterial = aTrack.GetMaterial();
254 
255  // check only for charged particles
256  if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
257  G4double xs = 0.0;
258  try
259  {
260  xs = aScaleFactor*
261  theCrossSectionDataStore->ComputeCrossSection(aParticle,aMaterial);
262  }
263  catch(G4HadronicException & aR)
264  {
266  aR.Report(ed);
267  DumpState(aTrack,"PostStepDoIt",ed);
268  ed << " Cross section is not available" << G4endl;
269  G4Exception("G4HadronicProcess::PostStepDoIt","had002",FatalException,ed);
270  }
271  if(xs <= 0.0 || xs < theLastCrossSection*G4UniformRand()) {
272  // No interaction
273  return theTotalResult;
274  }
275  }
276 
277  const G4Element* anElement = nullptr;
278  try
279  {
280  anElement = theCrossSectionDataStore->SampleZandA(aParticle, aMaterial,
281  targetNucleus);
282  }
283  catch(G4HadronicException & aR)
284  {
286  aR.Report(ed);
287  DumpState(aTrack,"SampleZandA",ed);
288  ed << " PostStepDoIt failed on element selection" << G4endl;
289  G4Exception("G4HadronicProcess::PostStepDoIt", "had003", FatalException,
290  ed);
291  }
292 
293  // Next check for illegal track status
294  //
295  if (aTrack.GetTrackStatus() != fAlive &&
296  aTrack.GetTrackStatus() != fSuspend) {
297  if (aTrack.GetTrackStatus() == fStopAndKill ||
299  aTrack.GetTrackStatus() == fPostponeToNextEvent) {
301  ed << "G4HadronicProcess: track in unusable state - "
302  << aTrack.GetTrackStatus() << G4endl;
303  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
304  DumpState(aTrack,"PostStepDoIt",ed);
305  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
306  }
307  // No warning for fStopButAlive which is a legal status here
308  return theTotalResult;
309  }
310 
311  // Initialize the hadronic projectile from the track
312  thePro.Initialise(aTrack);
313 
314  try
315  {
317  ChooseHadronicInteraction( thePro, targetNucleus, aMaterial, anElement );
318  }
319  catch(G4HadronicException & aE)
320  {
322  aE.Report(ed);
323  ed << "Target element "<<anElement->GetName()<<" Z= "
324  << targetNucleus.GetZ_asInt() << " A= "
326  DumpState(aTrack,"ChooseHadronicInteraction",ed);
327  ed << " No HadronicInteraction found out" << G4endl;
328  G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException,
329  ed);
330  }
331 
332  G4HadFinalState* result = nullptr;
333  G4int reentryCount = 0;
334 
335  do
336  {
337  try
338  {
339  // Save random engine if requested for debugging
342  }
343  // Call the interaction
345  ++reentryCount;
346  }
347  catch(G4HadronicException & aR)
348  {
350  aR.Report(ed);
351  ed << "Call for " << theInteraction->GetModelName() << G4endl;
352  ed << "Target element "<<anElement->GetName()<<" Z= "
354  << " A= " << targetNucleus.GetA_asInt() << G4endl;
355  DumpState(aTrack,"ApplyYourself",ed);
356  ed << " ApplyYourself failed" << G4endl;
357  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
358  ed);
359  }
360 
361  // Check the result for catastrophic energy non-conservation
362  result = CheckResult(thePro, targetNucleus, result);
363 
364  if(reentryCount>100) {
366  ed << "Call for " << theInteraction->GetModelName() << G4endl;
367  ed << "Target element "<<anElement->GetName()<<" Z= "
369  << " A= " << targetNucleus.GetA_asInt() << G4endl;
370  DumpState(aTrack,"ApplyYourself",ed);
371  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
372  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
373  ed);
374  }
375  }
376  while(!result); /* Loop checking, 30-Oct-2015, G.Folger */
377 
378  // Check whether kaon0 or anti_kaon0 are present between the secondaries:
379  // if this is the case, transform them into either kaon0S or kaon0L,
380  // with equal, 50% probability, keeping their dynamical masses (and
381  // the other kinematical properties).
382  // When this happens - very rarely - a "JustWarning" exception is thrown.
383  G4int nSec = result->GetNumberOfSecondaries();
384  if ( nSec > 0 ) {
385  for ( G4int i = 0; i < nSec; ++i ) {
386  G4DynamicParticle* dynamicParticle = result->GetSecondary(i)->GetParticle();
387  const G4ParticleDefinition* particleDefinition = dynamicParticle->GetParticleDefinition();
388  if ( particleDefinition == G4KaonZero::Definition() ||
389  particleDefinition == G4AntiKaonZero::Definition() ) {
390  G4ParticleDefinition* newParticleDefinition = G4KaonZeroShort::Definition();
391  if ( G4UniformRand() > 0.5 ) newParticleDefinition = G4KaonZeroLong::Definition();
392  G4double dynamicMass = dynamicParticle->GetMass();
393  dynamicParticle->SetDefinition( newParticleDefinition );
394  dynamicParticle->SetMass( dynamicMass );
396  ed << " Hadronic model " << theInteraction->GetModelName() << G4endl;
397  ed << " created " << particleDefinition->GetParticleName() << G4endl;
398  ed << " -> forced to be " << newParticleDefinition->GetParticleName() << G4endl;
399  G4Exception( "G4HadronicProcess::PostStepDoIt", "had007", JustWarning, ed );
400  //if ( std::abs( dynamicMass - particleDefinition->GetPDGMass() ) > 0.1*CLHEP::eV ) {
401  // G4cout << "\t dynamicMass=" << dynamicMass << " != PDGMass="
402  // << particleDefinition->GetPDGMass() << G4endl;
403  //}
404  }
405  }
406  }
407 
409 
411 
412  FillResult(result, aTrack);
413 
414  if (epReportLevel != 0) {
416  }
417  //G4cout << "PostStepDoIt done " << G4endl;
418  return theTotalResult;
419 }
420 
421 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const
422 {
423  outFile << "The description for this process has not been written yet.\n";
424 }
425 
426 
428 {
430  G4double biasedProbability = 1.-G4Exp(-nLTraversed);
431  G4double realProbability = 1-G4Exp(-nLTraversed/aScaleFactor);
432  G4double result = (biasedProbability-realProbability)/biasedProbability;
433  return result;
434 }
435 
437 {
439  G4double result =
440  1./aScaleFactor*G4Exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
441  return result;
442 }
443 
444 void
446 {
448 
449  G4double rotation = CLHEP::twopi*G4UniformRand();
450  G4ThreeVector it(0., 0., 1.);
451 
452  G4double efinal = aR->GetEnergyChange();
453  if(efinal < 0.0) { efinal = 0.0; }
454 
455  // check status of primary
456  if(aR->GetStatusChange() == stopAndKill) {
459 
460  // check its final energy
461  } else if(0.0 == efinal) {
464  ->GetAtRestProcessVector()->size() > 0)
467 
468  // primary is not killed apply rotation and Lorentz transformation
469  } else {
472  G4double newE = efinal + mass;
473  G4double newP = std::sqrt(efinal*(efinal + 2*mass));
474  G4ThreeVector newPV = newP*aR->GetMomentumChange();
475  G4LorentzVector newP4(newE, newPV);
476  newP4.rotate(rotation, it);
477  newP4 *= aR->GetTrafoToLab();
479  newE = newP4.e() - mass;
480  if(G4HadronicProcess_debug_flag && newE <= 0.0) {
482  DumpState(aT,"Primary has zero energy after interaction",ed);
483  G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning, ed);
484  }
485  if(newE < 0.0) { newE = 0.0; }
486  theTotalResult->ProposeEnergy( newE );
487  }
488  //G4cout << "FillResult: Efinal= " << efinal << " status= "
489  // << theTotalResult->GetTrackStatus()
490  // << " fKill= " << fStopAndKill << G4endl;
491 
492  // check secondaries: apply rotation and Lorentz transformation
493  nICelectrons = 0;
494  if(idxIC == -1) {
495  G4int idx = G4PhysicsModelCatalog::GetIndex("e-InternalConvertion");
496  idxIC = -1 == idx ? -2 : idx;
497  }
498  G4int nSec = aR->GetNumberOfSecondaries();
500 
501  if (nSec > 0) {
502  G4double time0 = aT.GetGlobalTime();
503  for (G4int i = 0; i < nSec; ++i) {
505  theM.rotate(rotation, it);
506  theM *= aR->GetTrafoToLab();
507  aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
508 
509  if(idxIC == aR->GetSecondary(i)->GetCreatorModelType())
510  { ++nICelectrons; }
511 
512  // time of interaction starts from zero
513  G4double time = aR->GetSecondary(i)->GetTime();
514  if (time < 0.0) { time = 0.0; }
515 
516  // take into account global time
517  time += time0;
518 
519  G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
520  time, aT.GetPosition());
522  G4double newWeight = fWeight*aR->GetSecondary(i)->GetWeight();
523  // G4cout << "#### ParticleDebug "
524  // <<GetProcessName()<<" "
525  //<<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()
526  //<<" "<<aScaleFactor<<" "
527  // <<XBiasSurvivalProbability()<<" "
528  // <<XBiasSecondaryWeight()<<" "
529  // <<aT.GetWeight()<<" "
530  // <<aR->GetSecondary(i)->GetWeight()<<" "
531  // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
532  // <<G4endl;
533  track->SetWeight(newWeight);
537  G4double e = track->GetKineticEnergy();
538  if (e <= 0.0) {
540  DumpState(aT,"Secondary has zero energy",ed);
541  ed << "Secondary " << track->GetDefinition()->GetParticleName()
542  << G4endl;
543  G4Exception("G4HadronicProcess::FillResults", "had011",
544  JustWarning,ed);
545  }
546  }
547  }
548  }
549  aR->Clear();
550 }
551 
553 {
554  BiasCrossSectionByFactor(factor);
555 }
556 
558 {
559  if (aScale <= 0.0) {
561  ed << " Wrong biasing factor " << aScale << " for " << GetProcessName();
562  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010",
563  JustWarning, ed, "Cross-section bias is ignored");
564  } else {
565  xBiasOn = true;
566  aScaleFactor = aScale;
567  }
568 }
569 
571  const G4Nucleus &aNucleus,
573 {
574  // check for catastrophic energy non-conservation
575  // to re-sample the interaction
576 
578  G4double nuclearMass(0);
579  if (theModel) {
580 
581  // Compute final-state total energy
582  G4double finalE(0.);
583  G4int nSec = result->GetNumberOfSecondaries();
584 
585  nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
586  aNucleus.GetZ_asInt());
587  if (result->GetStatusChange() != stopAndKill) {
588  // Interaction didn't complete, returned "do nothing" state
589  // and reset nucleus or the primary survived the interaction
590  // (e.g. electro-nuclear ) => keep nucleus
591  finalE=result->GetLocalEnergyDeposit() +
592  aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
593  if( nSec == 0 ){
594  // Since there are no secondaries, there is no recoil nucleus.
595  // To check energy balance we must neglect the initial nucleus too.
596  nuclearMass=0.0;
597  }
598  }
599  for (G4int i = 0; i < nSec; i++) {
600  G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle();
601  finalE += pdyn->GetTotalEnergy();
602  G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass();
603  G4double mass_dyn=pdyn->GetMass();
604  if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV){
605  result->Clear();
606  result = nullptr;
608  desc << "Warning: Secondary with off-shell dynamic mass detected: "
609  << G4endl
610  << " " << pdyn->GetDefinition()->GetParticleName()
611  << ", PDG mass: " << mass_pdg << ", dynamic mass: "
612  << mass_dyn << G4endl
613  << (epReportLevel<0 ? "abort the event"
614  : "re-sample the interaction") << G4endl
615  << " Process / Model: " << GetProcessName()<< " / "
616  << theModel->GetModelName() << G4endl
617  << " Primary: " << aPro.GetDefinition()->GetParticleName()
618  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
619  << " E= " << aPro.Get4Momentum().e()
620  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
621  << aNucleus.GetA_asInt() << ")" << G4endl;
622  G4Exception("G4HadronicProcess:CheckResult()", "had012",
624  // must return here.....
625  return result;
626  }
627  }
628  G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
629 
630  std::pair<G4double, G4double> checkLevels =
631  theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
632  if (std::abs(deltaE) > checkLevels.second &&
633  std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
634  // do not delete result, this is a pointer to a data member;
635  result->Clear();
636  result = nullptr;
638  desc << "Warning: Bad energy non-conservation detected, will "
639  << (epReportLevel<0 ? "abort the event"
640  : "re-sample the interaction") << G4endl
641  << " Process / Model: " << GetProcessName()<< " / "
642  << theModel->GetModelName() << G4endl
643  << " Primary: " << aPro.GetDefinition()->GetParticleName()
644  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
645  << " E= " << aPro.Get4Momentum().e()
646  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
647  << aNucleus.GetA_asInt() << ")" << G4endl
648  << " E(initial - final) = " << deltaE << " MeV." << G4endl;
649  G4Exception("G4HadronicProcess:CheckResult()", "had012",
651  }
652  }
653  return result;
654 }
655 
656 void
658  const G4Nucleus& aNucleus)
659 {
660  G4int target_A=aNucleus.GetA_asInt();
661  G4int target_Z=aNucleus.GetZ_asInt();
662  G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
663  G4LorentzVector target4mom(0, 0, 0, targetMass
665 
666  G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
667  G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
668  G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
669 
670  G4int initial_A = target_A + track_A;
671  G4int initial_Z = target_Z + track_Z - nICelectrons;
672 
673  G4LorentzVector initial4mom = projectile4mom + target4mom;
674 
675  // Compute final-state momentum for scattering and "do nothing" results
676  G4LorentzVector final4mom;
677  G4int final_A(0), final_Z(0);
678 
680  if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
681  // Either interaction didn't complete, returned "do nothing" state
682  // or the primary survived the interaction (e.g. electro-nucleus )
683  G4Track temp(aTrack);
684 
685  // Use the final energy / momentum
688 
689  if( nSec == 0 ){
690  // Interaction didn't complete, returned "do nothing" state
691  // - or suppressed recoil (e.g. Neutron elastic )
692  final4mom = temp.GetDynamicParticle()->Get4Momentum() + target4mom;
693  final_A = initial_A;
694  final_Z = initial_Z;
695  }else{
696  // The primary remains in final state (e.g. electro-nucleus )
697  final4mom = temp.GetDynamicParticle()->Get4Momentum();
698  final_A = track_A;
699  final_Z = track_Z;
700  // Expect that the target nucleus will have interacted,
701  // and its products, including recoil, will be included in secondaries.
702  }
703  }
704  if( nSec > 0 ) {
705  G4Track* sec;
706 
707  for (G4int i = 0; i < nSec; i++) {
708  sec = theTotalResult->GetSecondary(i);
709  final4mom += sec->GetDynamicParticle()->Get4Momentum();
710  final_A += sec->GetDefinition()->GetBaryonNumber();
711  final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
712  }
713  }
714 
715  // Get level-checking information (used to cut-off relative checks)
716  G4String processName = GetProcessName();
718  G4String modelName("none");
719  if (theModel) modelName = theModel->GetModelName();
720  std::pair<G4double, G4double> checkLevels = epCheckLevels;
721  if (!levelsSetByProcess) {
722  if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
723  checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
724  checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
725  }
726 
727  // Compute absolute total-energy difference, and relative kinetic-energy
728  G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
729 
730  G4LorentzVector diff = initial4mom - final4mom;
731  G4double absolute = diff.e();
732  G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
733 
734  G4double absolute_mom = diff.vect().mag();
735  G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
736 
737  // Evaluate relative and absolute conservation
738  G4bool relPass = true;
739  G4String relResult = "pass";
740  if ( std::abs(relative) > checkLevels.first
741  || std::abs(relative_mom) > checkLevels.first) {
742  relPass = false;
743  relResult = checkRelative ? "fail" : "N/A";
744  }
745 
746  G4bool absPass = true;
747  G4String absResult = "pass";
748  if ( std::abs(absolute) > checkLevels.second
749  || std::abs(absolute_mom) > checkLevels.second ) {
750  absPass = false ;
751  absResult = "fail";
752  }
753 
754  G4bool chargePass = true;
755  G4String chargeResult = "pass";
756  if ( (initial_A-final_A)!=0
757  || (initial_Z-final_Z)!=0 ) {
758  chargePass = checkLevels.second < DBL_MAX ? false : true;
759  chargeResult = "fail";
760  }
761 
762  G4bool conservationPass = (relPass || absPass) && chargePass;
763 
764  std::stringstream Myout;
765  G4bool Myout_notempty(false);
766  // Options for level of reporting detail:
767  // 0. off
768  // 1. report only when E/p not conserved
769  // 2. report regardless of E/p conservation
770  // 3. report only when E/p not conserved, with model names, process names, and limits
771  // 4. report regardless of E/p conservation, with model names, process names, and limits
772  // negative -1.., as above, but send output to stderr
773 
774  if( std::abs(epReportLevel) == 4
775  || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
776  Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
777  Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
778  << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
779  << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
780  << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
781  << aNucleus.GetA_asInt() << ")" << G4endl;
782  Myout_notempty=true;
783  }
784  if ( std::abs(epReportLevel) == 4
785  || std::abs(epReportLevel) == 2
786  || ! conservationPass ){
787 
788  Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
789  << relative << " p/p(0)= " << relative_mom << G4endl;
790  Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
791  << absolute/MeV << " / " << absolute_mom/MeV << " 3mom: " << (diff.vect())*1./MeV << G4endl;
792  Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
793  Myout_notempty=true;
794 
795  }
796  Myout.flush();
797  if ( Myout_notempty ) {
798  if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
799  else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
800  }
801 }
802 
804  const G4String& method,
806 {
807  ed << "Unrecoverable error in the method " << method << " of "
808  << GetProcessName() << G4endl;
809  ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
810  << aTrack.GetParentID()
811  << " " << aTrack.GetParticleDefinition()->GetParticleName()
812  << G4endl;
813  ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
814  << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
815  ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
816 
817  if (aTrack.GetMaterial()) {
818  ed << " material " << aTrack.GetMaterial()->GetName();
819  }
820  ed << G4endl;
821 
822  if (aTrack.GetVolume()) {
823  ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
824  << ">" << G4endl;
825  }
826 }
827 
829 {
831 }
832 
834 {
836 }
837 
838 std::vector<G4HadronicInteraction*>&
840 {
842 }
void SetMass(G4double mass)
G4double GetKineticEnergy() const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4ParticleChange * theTotalResult
void SetTrafoToLab(const G4LorentzRotation &aT)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int GetCreatorModelType() const
static G4AntiKaonZero * Definition()
G4double GetWeight() const
G4int GetNumberOfSecondaries() const
G4LorentzRotation & GetTrafoToLab()
G4HadProjectile thePro
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:45
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
const G4ThreeVector & GetMomentumChange() const
static constexpr double MeV
Definition: G4SIunits.hh:214
void SetWeight(G4double aValue)
G4CrossSectionDataStore * theCrossSectionDataStore
const G4String & GetName() const
Definition: G4Element.hh:127
G4HadSecondary * GetSecondary(size_t i)
void ProposeWeight(G4double finalWeight)
G4HadronicInteraction * GetHadronicInteraction() const
void SetCreatorModelIndex(G4int idx)
void RegisterMe(G4HadronicInteraction *a)
#define G4endl
Definition: G4ios.hh:61
const char * p
Definition: xmltok.h:285
G4EnergyRangeManager theEnergyRangeManager
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
void DumpPhysicsTable(const G4ParticleDefinition &)
G4double ComputeCrossSection(const G4DynamicParticle *, const G4Material *)
G4HadFinalStateStatus GetStatusChange() const
void RegisterParticle(G4HadronicProcess *, const G4ParticleDefinition *)
void Register(G4HadronicProcess *)
const G4String & GetParticleName() const
G4int GetTrackID() const
const G4TouchableHandle & GetTouchableHandle() const
G4double GetPDGCharge() const
void RegisterMe(G4HadronicInteraction *a)
void SetSecondaryWeightByProcess(G4bool)
G4HadronicProcessType
static G4KaonZeroLong * Definition()
void RegisterInteraction(G4HadronicProcess *, G4HadronicInteraction *)
static void saveEngineStatus(const char filename[]="Config.conf")
Definition: Random.cc:275
const G4ParticleDefinition * GetParticleDefinition() const
void PreparePhysicsTable(const G4ParticleDefinition &) override
std::vector< G4HadronicInteraction * > & GetHadronicInteractionList()
void ProcessDescription(std::ostream &outFile) const override
G4double GetTotalEnergy() const
G4double XBiasSurvivalProbability()
G4ParticleDefinition * GetDefinition() const
G4double GetPDGMass() const
void AddSecondary(G4Track *aSecondary)
G4double GetElementCrossSection(const G4DynamicParticle *part, const G4Element *elm, const G4Material *mat=nullptr)
G4double GetWeight() const
const G4int nmax
G4HadronicProcessStore * theProcessStore
static G4int GetIndex(const G4String &)
const G4String & GetName() const
Definition: G4Material.hh:179
static G4KaonZero * Definition()
Definition: G4KaonZero.cc:53
void Initialise(const G4Track &aT)
G4ThreeVector GetMomentum() const
void BuildPhysicsTable(const G4ParticleDefinition &) override
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
TString part[npart]
Definition: Style.C:32
G4double GetTime() const
G4double theInitialNumberOfInteractionLength
G4ParticleDefinition * GetDefinition() const
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *) override
G4double GetTotalNumberOfInteractionLengthTraversed() const
Definition: G4VProcess.hh:461
void GetEnergyMomentumCheckEnvvars()
virtual void Initialize(const G4Track &)
void DumpPhysicsTable(const G4ParticleDefinition &p)
G4ProcessType
static constexpr double electron_mass_c2
std::pair< G4double, G4double > epCheckLevels
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
const G4ThreeVector * GetMomentumDirection() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
#define G4UniformRand()
Definition: Randomize.hh:53
G4double GetGlobalTime() const
void Set4Momentum(const G4LorentzVector &momentum)
G4LorentzVector Get4Momentum() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
const G4String & GetProcessName() const
Definition: G4VProcess.hh:411
const G4ThreeVector & GetPosition() const
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result)
void BiasCrossSectionByFactor(G4double aScale)
static G4double GetNuclearMass(const G4double A, const G4double Z)
void SetKineticEnergy(const G4double aValue)
Definition: G4Step.hh:76
void ProposeEnergy(G4double finalEnergy)
Float_t mat
void BuildPhysicsTable(const G4ParticleDefinition &)
void AddDataSet(G4VCrossSectionDataSet *)
void MultiplyCrossSectionBy(G4double factor)
G4Material * GetMaterial() const
void FillResult(G4HadFinalState *aR, const G4Track &aT)
G4GLOB_DLL std::ostream G4cerr
const G4LorentzVector & Get4Momentum() const
G4int GetNumberOfSecondaries() const
G4int GetZasInt() const
Definition: G4Element.hh:132
~G4HadronicProcess() override
void PrintInfo(const G4ParticleDefinition *)
static G4HadronicProcessStore * Instance()
Hep3Vector unit() const
G4double G4ParticleHPJENDLHEData::G4double result
G4double GetKineticEnergy() const
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4double GetLocalEnergyDeposit() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
const G4ParticleDefinition * GetDefinition() const
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
void CheckEnergyMomentumConservation(const G4Track &, const G4Nucleus &)
double mag() const
void AddDataSet(G4VCrossSectionDataSet *aDataSet)
int G4lrint(double ad)
Definition: templates.hh:151
int G4int
Definition: G4Types.hh:78
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4double XBiasSecondaryWeight()
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4ProcessManager * GetProcessManager() const
static G4KaonZeroShort * Definition()
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void BuildPhysicsTable(const G4ParticleDefinition &)
G4TrackStatus GetTrackStatus() const
G4DynamicParticle * GetParticle()
static constexpr double mm
Definition: SystemOfUnits.h:95
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:435
G4ForceCondition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4int size() const
G4GLOB_DLL std::ostream G4cout
G4Track * GetSecondary(G4int anIndex) const
G4HadronicProcess(const G4String &processName="Hadronic", G4ProcessType procType=fHadronic)
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:450
void DeRegister(G4HadronicProcess *)
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Hep3Vector vect() const
void Report(std::ostream &aS)
G4double GetCrossSection(const G4DynamicParticle *, const G4Material *)
G4double GetMass() const
G4double GetTotalEnergy() const
G4double GetEnergy() const
const G4ParticleDefinition * GetParticleDefinition() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double GetEnergyChange() const
G4int first(char) const
G4TrackStatus GetTrackStatus() const
const G4ThreeVector & GetMomentumDirection() const
#define DBL_MAX
Definition: templates.hh:83
G4HadronicInteraction * theInteraction
void ProposeTrackStatus(G4TrackStatus status)
const G4LorentzRotation & GetTrafoToLab() const
static constexpr double GeV
const G4DynamicParticle * GetDynamicParticle() const
G4HadronicInteraction * ChooseHadronicInteraction(const G4HadProjectile &aHadProjectile, G4Nucleus &aTargetNucleus, const G4Material *aMaterial, const G4Element *anElement)
void SetMomentumDirection(const G4ThreeVector &aValue)
static const char * G4Hadronic_Random_File
std::vector< G4HadronicInteraction * > & GetHadronicInteractionList()
const G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
const G4String & GetName() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VPhysicalVolume * GetVolume() const
G4int GetParentID() const
const G4String & GetModelName() const
HepLorentzVector & rotate(double, const Hep3Vector &)
void SetNumberOfSecondaries(G4int totSecondaries)