Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4HadronicProcess.cc
Go to the documentation of this file.
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 98000 2016-06-30 13:01:05Z 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"
70 
71 #include "G4AutoLock.hh"
72 #include "G4NistManager.hh"
73 
74 #include <typeinfo>
75 #include <sstream>
76 #include <iostream>
77 
78 #include <stdlib.h>
79 
80 // File-scope variable to capture environment variable at startup
81 
82 static const char* G4Hadronic_Random_File = getenv("G4HADRONIC_RANDOM_FILE");
83 
85 
87  G4ProcessType procType)
88  : G4VDiscreteProcess(processName, procType)
89 {
90  SetProcessSubType(fHadronInelastic); // Default unless subclass changes
91 
94  theInteraction = nullptr;
95  theCrossSectionDataStore = new G4CrossSectionDataStore();
96  theProcessStore = G4HadronicProcessStore::Instance();
97  theProcessStore->Register(this);
98  theInitialNumberOfInteractionLength = 0.0;
99  aScaleFactor = 1;
100  xBiasOn = false;
101  nMatWarn = 0;
102  useIntegralXS = true;
103  theLastCrossSection = 0.0;
104  G4HadronicProcess_debug_flag = false;
105  GetEnergyMomentumCheckEnvvars();
106 }
107 
109 
111  G4HadronicProcessType aHadSubType)
112  : G4VDiscreteProcess(processName, fHadronic)
113 {
114  SetProcessSubType(aHadSubType);
115 
118  theInteraction = nullptr;
119  theCrossSectionDataStore = new G4CrossSectionDataStore();
120  theProcessStore = G4HadronicProcessStore::Instance();
121  theProcessStore->Register(this);
122  theInitialNumberOfInteractionLength = 0.0;
123  aScaleFactor = 1;
124  xBiasOn = false;
125  useIntegralXS = true;
126  nMatWarn = 0;
127  theLastCrossSection = 0.0;
128  G4HadronicProcess_debug_flag = false;
129  GetEnergyMomentumCheckEnvvars();
130 }
131 
132 
134 {
135  theProcessStore->DeRegister(this);
136  delete theTotalResult;
137  delete theCrossSectionDataStore;
138 }
139 
140 void G4HadronicProcess::GetEnergyMomentumCheckEnvvars() {
141  levelsSetByProcess = false;
142 
143  epReportLevel = getenv("G4Hadronic_epReportLevel") ?
144  strtol(getenv("G4Hadronic_epReportLevel"),0,10) : 0;
145 
146  epCheckLevels.first = getenv("G4Hadronic_epCheckRelativeLevel") ?
147  strtod(getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX;
148 
149  epCheckLevels.second = getenv("G4Hadronic_epCheckAbsoluteLevel") ?
150  strtod(getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX;
151 }
152 
154 {
155  if(!a) { return; }
156  try{ theEnergyRangeManager.RegisterMe( a ); }
157  catch(G4HadronicException & aE)
158  {
160  aE.Report(ed);
161  ed << "Unrecoverable error in " << GetProcessName()
162  << " to register " << a->GetModelName() << G4endl;
163  G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException,
164  ed);
165  }
167 }
168 
169 G4double
171  const G4Element * elm,
172  const G4Material* mat)
173 {
174  G4Material* aMaterial = const_cast<G4Material*>(mat);
175  if(!aMaterial)
176  {
177  // Because NeutronHP needs a material pointer (for instance to get the
178  // temperature), we ask the Nist manager to find a simple material
179  // from the (integer) Z of the element.
180  aMaterial =
182  if(!aMaterial) {
183  ++nMatWarn;
184  static const G4int nmax = 5;
185  if(nMatWarn < nmax) {
187  ed << "Cannot compute Element x-section for " << GetProcessName()
188  << " because no material defined \n"
189  << " Please, specify material pointer or define simple material"
190  << " for Z= " << elm->GetZasInt();
191  G4Exception("G4HadronicProcess::GetElementCrossSection", "had066", JustWarning,
192  ed);
193  }
194  }
195  }
196  G4double x =
197  std::max(theCrossSectionDataStore->GetCrossSection(part, elm, aMaterial),0.0);
198  return x;
199 }
200 
202 {
203  if(getenv("G4HadronicProcess_debug")) {
204  G4HadronicProcess_debug_flag = true;
205  }
206  theProcessStore->RegisterParticle(this, &p);
207 }
208 
210 {
211  try
212  {
213  theCrossSectionDataStore->BuildPhysicsTable(p);
214  theEnergyRangeManager.BuildPhysicsTable(p);
215  }
216  catch(G4HadronicException aR)
217  {
219  aR.Report(ed);
220  ed << " hadronic initialisation fails";
221  G4Exception("G4HadronicProcess::BuildPhysicsTable", "had000",
222  FatalException,ed);
223  }
225 }
226 
229 {
230  //G4cout << "GetMeanFreePath " << aTrack.GetDefinition()->GetParticleName()
231  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
232  try
233  {
234  theLastCrossSection = aScaleFactor*
235  theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(),
236  aTrack.GetMaterial());
237  }
238  catch(G4HadronicException aR)
239  {
241  aR.Report(ed);
242  DumpState(aTrack,"GetMeanFreePath",ed);
243  ed << " Cross section is not available" << G4endl;
244  G4Exception("G4HadronicProcess::GetMeanFreePath", "had002", FatalException,
245  ed);
246  }
247  G4double res = (theLastCrossSection > 0.0) ? 1.0/theLastCrossSection : DBL_MAX;
248  //G4cout << " xsection= " << res << G4endl;
249  return res;
250 }
251 
254 {
255  //G4cout << "PostStepDoIt " << aTrack.GetDefinition()->GetParticleName()
256  // << " Ekin= " << aTrack.GetKineticEnergy() << G4endl;
257  // if primary is not Alive then do nothing
259  theTotalResult->Initialize(aTrack);
261  if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
262 
263  // Find cross section at end of step and check if <= 0
264  //
265  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
266  G4Material* aMaterial = aTrack.GetMaterial();
267 
268  // check only for charged particles
269  if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
270  G4double xs = 0.0;
271  try
272  {
273  xs = aScaleFactor*theCrossSectionDataStore->GetCrossSection(aParticle,aMaterial);
274  }
275  catch(G4HadronicException aR)
276  {
278  aR.Report(ed);
279  DumpState(aTrack,"PostStepDoIt",ed);
280  ed << " Cross section is not available" << G4endl;
281  G4Exception("G4HadronicProcess::PostStepDoIt", "had002", FatalException,ed);
282  }
283  if(xs <= 0.0 || (useIntegralXS && xs < theLastCrossSection*G4UniformRand())) {
284  // No interaction
285  return theTotalResult;
286  }
287  }
288 
289  G4Element* anElement = nullptr;
290  try
291  {
292  anElement = theCrossSectionDataStore->SampleZandA(aParticle,
293  aMaterial,
294  targetNucleus);
295  }
296  catch(G4HadronicException & aR)
297  {
299  aR.Report(ed);
300  DumpState(aTrack,"SampleZandA",ed);
301  ed << " PostStepDoIt failed on element selection" << G4endl;
302  G4Exception("G4HadronicProcess::PostStepDoIt", "had003", FatalException,
303  ed);
304  }
305 
306  // Next check for illegal track status
307  //
308  if (aTrack.GetTrackStatus() != fAlive &&
309  aTrack.GetTrackStatus() != fSuspend) {
310  if (aTrack.GetTrackStatus() == fStopAndKill ||
312  aTrack.GetTrackStatus() == fPostponeToNextEvent) {
314  ed << "G4HadronicProcess: track in unusable state - "
315  << aTrack.GetTrackStatus() << G4endl;
316  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
317  DumpState(aTrack,"PostStepDoIt",ed);
318  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
319  }
320  // No warning for fStopButAlive which is a legal status here
321  return theTotalResult;
322  }
323 
324  // Initialize the hadronic projectile from the track
325  thePro.Initialise(aTrack);
326 
327  try
328  {
329  theInteraction =
330  ChooseHadronicInteraction( thePro, targetNucleus, aMaterial, anElement );
331  }
332  catch(G4HadronicException & aE)
333  {
335  aE.Report(ed);
336  ed << "Target element "<<anElement->GetName()<<" Z= "
337  << targetNucleus.GetZ_asInt() << " A= "
338  << targetNucleus.GetA_asInt() << G4endl;
339  DumpState(aTrack,"ChooseHadronicInteraction",ed);
340  ed << " No HadronicInteraction found out" << G4endl;
341  G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException,
342  ed);
343  }
344 
345  G4HadFinalState* result = nullptr;
346  G4int reentryCount = 0;
347 
348  do
349  {
350  try
351  {
352  // Save random engine if requested for debugging
355  }
356  // Call the interaction
357  result = theInteraction->ApplyYourself( thePro, targetNucleus);
358  ++reentryCount;
359  }
360  catch(G4HadronicException aR)
361  {
363  aR.Report(ed);
364  ed << "Call for " << theInteraction->GetModelName() << G4endl;
365  ed << "Target element "<<anElement->GetName()<<" Z= "
366  << targetNucleus.GetZ_asInt()
367  << " A= " << targetNucleus.GetA_asInt() << G4endl;
368  DumpState(aTrack,"ApplyYourself",ed);
369  ed << " ApplyYourself failed" << G4endl;
370  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
371  ed);
372  }
373 
374  // Check the result for catastrophic energy non-conservation
375  CheckResult(thePro, targetNucleus, result);
376 
377  if(reentryCount>100) {
379  ed << "Call for " << theInteraction->GetModelName() << G4endl;
380  ed << "Target element "<<anElement->GetName()<<" Z= "
381  << targetNucleus.GetZ_asInt()
382  << " A= " << targetNucleus.GetA_asInt() << G4endl;
383  DumpState(aTrack,"ApplyYourself",ed);
384  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
385  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
386  ed);
387  }
388  }
389  while(!result); /* Loop checking, 30-Oct-2015, G.Folger */
390 
391  // Check whether kaon0 or anti_kaon0 are present between the secondaries:
392  // if this is the case, transform them into either kaon0S or kaon0L,
393  // with equal, 50% probability, keeping their dynamical masses (and
394  // the other kinematical properties).
395  // When this happens - very rarely - a "JustWarning" exception is thrown.
396  G4int nSec = result->GetNumberOfSecondaries();
397  if ( nSec > 0 ) {
398  for ( G4int i = 0; i < nSec; ++i ) {
399  G4DynamicParticle* dynamicParticle = result->GetSecondary(i)->GetParticle();
400  const G4ParticleDefinition* particleDefinition = dynamicParticle->GetParticleDefinition();
401  if ( particleDefinition == G4KaonZero::Definition() ||
402  particleDefinition == G4AntiKaonZero::Definition() ) {
403  G4ParticleDefinition* newParticleDefinition = G4KaonZeroShort::Definition();
404  if ( G4UniformRand() > 0.5 ) newParticleDefinition = G4KaonZeroLong::Definition();
405  G4double dynamicMass = dynamicParticle->GetMass();
406  dynamicParticle->SetDefinition( newParticleDefinition );
407  dynamicParticle->SetMass( dynamicMass );
409  ed << " Hadronic model " << theInteraction->GetModelName() << G4endl;
410  ed << " created " << particleDefinition->GetParticleName() << G4endl;
411  ed << " -> forced to be " << newParticleDefinition->GetParticleName() << G4endl;
412  G4Exception( "G4HadronicProcess::PostStepDoIt", "had007", JustWarning, ed );
413  //if ( std::abs( dynamicMass - particleDefinition->GetPDGMass() ) > 0.1*CLHEP::eV ) {
414  // G4cout << "\t dynamicMass=" << dynamicMass << " != PDGMass="
415  // << particleDefinition->GetPDGMass() << G4endl;
416  //}
417  }
418  }
419  }
420 
422 
424 
425  FillResult(result, aTrack);
426 
427  if (epReportLevel != 0) {
428  CheckEnergyMomentumConservation(aTrack, targetNucleus);
429  }
430  //G4cout << "PostStepDoIt done " << G4endl;
431  return theTotalResult;
432 }
433 
434 
435 void G4HadronicProcess::ProcessDescription(std::ostream& outFile) const
436 {
437  outFile << "The description for this process has not been written yet.\n";
438 }
439 
440 
441 G4double G4HadronicProcess::XBiasSurvivalProbability()
442 {
443  G4double result = 0;
445  G4double biasedProbability = 1.-std::exp(-nLTraversed);
446  G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
447  result = (biasedProbability-realProbability)/biasedProbability;
448  return result;
449 }
450 
451 G4double G4HadronicProcess::XBiasSecondaryWeight()
452 {
453  G4double result = 0;
455  result =
456  1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
457  return result;
458 }
459 
460 void
462 {
464 
465  G4double rotation = CLHEP::twopi*G4UniformRand();
466  G4ThreeVector it(0., 0., 1.);
467 
468  G4double efinal = aR->GetEnergyChange();
469  if(efinal < 0.0) { efinal = 0.0; }
470 
471  // check status of primary
472  if(aR->GetStatusChange() == stopAndKill) {
475 
476  // check its final energy
477  } else if(0.0 == efinal) {
480  ->GetAtRestProcessVector()->size() > 0)
483 
484  // primary is not killed apply rotation and Lorentz transformation
485  } else {
488  G4double newE = efinal + mass;
489  G4double newP = std::sqrt(efinal*(efinal + 2*mass));
490  G4ThreeVector newPV = newP*aR->GetMomentumChange();
491  G4LorentzVector newP4(newE, newPV);
492  newP4.rotate(rotation, it);
493  newP4 *= aR->GetTrafoToLab();
495  newE = newP4.e() - mass;
496  if(G4HadronicProcess_debug_flag && newE <= 0.0) {
498  DumpState(aT,"Primary has zero energy after interaction",ed);
499  G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning, ed);
500  }
501  if(newE < 0.0) { newE = 0.0; }
502  theTotalResult->ProposeEnergy( newE );
503  }
504  //G4cout << "FillResult: Efinal= " << efinal << " status= "
505  // << theTotalResult->GetTrackStatus()
506  // << " fKill= " << fStopAndKill << G4endl;
507 
508  // check secondaries: apply rotation and Lorentz transformation
509  G4int nSec = aR->GetNumberOfSecondaries();
511  G4double weight = aT.GetWeight();
512 
513  if (nSec > 0) {
514  G4double time0 = aT.GetGlobalTime();
515  for (G4int i = 0; i < nSec; ++i) {
517  theM.rotate(rotation, it);
518  theM *= aR->GetTrafoToLab();
519  aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
520 
521  // time of interaction starts from zero
522  G4double time = aR->GetSecondary(i)->GetTime();
523  if (time < 0.0) { time = 0.0; }
524 
525  // take into account global time
526  time += time0;
527 
528  G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
529  time, aT.GetPosition());
531  G4double newWeight = weight*aR->GetSecondary(i)->GetWeight();
532  // G4cout << "#### ParticleDebug "
533  // <<GetProcessName()<<" "
534  //<<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
535  // <<aScaleFactor<<" "
536  // <<XBiasSurvivalProbability()<<" "
537  // <<XBiasSecondaryWeight()<<" "
538  // <<aT.GetWeight()<<" "
539  // <<aR->GetSecondary(i)->GetWeight()<<" "
540  // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
541  // <<G4endl;
542  track->SetWeight(newWeight);
545  if (G4HadronicProcess_debug_flag) {
546  G4double e = track->GetKineticEnergy();
547  if (e <= 0.0) {
549  DumpState(aT,"Secondary has zero energy",ed);
550  ed << "Secondary " << track->GetDefinition()->GetParticleName()
551  << G4endl;
552  G4Exception("G4HadronicProcess::FillResults", "had011",
553  JustWarning,ed);
554  }
555  }
556  }
557  }
558 
559  aR->Clear();
560 }
561 
563 {
564  xBiasOn = true;
565  aScaleFactor = aScale;
566  G4String it = GetProcessName();
567  if ((it != "photonNuclear") &&
568  (it != "electronNuclear") &&
569  (it != "positronNuclear") ) {
571  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had009",
572  FatalException, ed,
573  "Cross-section biasing available only for gamma and electro nuclear reactions.");
574  }
575 
576  if (aScale < 100) {
578  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010", JustWarning,ed,
579  "Cross-section bias readjusted to be above safe limit. New value is 100");
580  aScaleFactor = 100.;
581  }
582 }
583 
585  const G4Nucleus &aNucleus,
586  G4HadFinalState * result)
587 {
588  // check for catastrophic energy non-conservation
589  // to re-sample the interaction
590 
592  G4double nuclearMass(0);
593  if (theModel) {
594 
595  // Compute final-state total energy
596  G4double finalE(0.);
597  G4int nSec = result->GetNumberOfSecondaries();
598 
599  nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
600  aNucleus.GetZ_asInt());
601  if (result->GetStatusChange() != stopAndKill) {
602  // Interaction didn't complete, returned "do nothing" state
603  // and reset nucleus or the primary survived the interaction
604  // (e.g. electro-nuclear ) => keep nucleus
605  finalE=result->GetLocalEnergyDeposit() +
606  aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
607  if( nSec == 0 ){
608  // Since there are no secondaries, there is no recoil nucleus.
609  // To check energy balance we must neglect the initial nucleus too.
610  nuclearMass=0.0;
611  }
612  }
613  for (G4int i = 0; i < nSec; i++) {
614  G4DynamicParticle *pdyn=result->GetSecondary(i)->GetParticle();
615  finalE += pdyn->GetTotalEnergy();
616  G4double mass_pdg=pdyn->GetDefinition()->GetPDGMass();
617  G4double mass_dyn=pdyn->GetMass();
618  if ( std::abs(mass_pdg - mass_dyn) > 0.1*mass_pdg + 1.*MeV){
619  result->Clear();
620  result = 0;
622  desc << "Warning: Secondary with off-shell dynamic mass detected: " << G4endl
623  << " " << pdyn->GetDefinition()->GetParticleName()
624  << ", PDG mass: " << mass_pdg << ", dynamic mass: "<< mass_dyn << G4endl
625  << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl
626  << " Process / Model: " << GetProcessName()<< " / "
627  << theModel->GetModelName() << G4endl
628  << " Primary: " << aPro.GetDefinition()->GetParticleName()
629  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
630  << " E= " << aPro.Get4Momentum().e()
631  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
632  << aNucleus.GetA_asInt() << ")" << G4endl;
633  G4Exception("G4HadronicProcess:CheckResult()", "had012",
635  // must return here.....
636  return result;
637  }
638  }
639  G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
640 
641  std::pair<G4double, G4double> checkLevels =
642  theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
643  if (std::abs(deltaE) > checkLevels.second &&
644  std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
645  // do not delete result, this is a pointer to a data member;
646  result->Clear();
647  result = 0;
649  desc << "Warning: Bad energy non-conservation detected, will "
650  << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl
651  << " Process / Model: " << GetProcessName()<< " / "
652  << theModel->GetModelName() << G4endl
653  << " Primary: " << aPro.GetDefinition()->GetParticleName()
654  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "), "
655  << " E= " << aPro.Get4Momentum().e()
656  << ", target nucleus (" << aNucleus.GetZ_asInt() << ", "
657  << aNucleus.GetA_asInt() << ")" << G4endl
658  << " E(initial - final) = " << deltaE << " MeV." << G4endl;
659  G4Exception("G4HadronicProcess:CheckResult()", "had012",
661  }
662  }
663  return result;
664 }
665 
666 void
668  const G4Nucleus& aNucleus)
669 {
670  G4int target_A=aNucleus.GetA_asInt();
671  G4int target_Z=aNucleus.GetZ_asInt();
672  G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
673  G4LorentzVector target4mom(0, 0, 0, targetMass);
674 
675  G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
676  G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
677  G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
678 
679  G4int initial_A = target_A + track_A;
680  G4int initial_Z = target_Z + track_Z;
681 
682  G4LorentzVector initial4mom = projectile4mom + target4mom;
683 
684  // Compute final-state momentum for scattering and "do nothing" results
685  G4LorentzVector final4mom;
686  G4int final_A(0), final_Z(0);
687 
689  if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
690  // Either interaction didn't complete, returned "do nothing" state
691  // or the primary survived the interaction (e.g. electro-nucleus )
692  G4Track temp(aTrack);
693 
694  // Use the final energy / momentum
697 
698  if( nSec == 0 ){
699  // Interaction didn't complete, returned "do nothing" state
700  // - or suppressed recoil (e.g. Neutron elastic )
701  final4mom = temp.GetDynamicParticle()->Get4Momentum() + target4mom;
702  final_A = initial_A;
703  final_Z = initial_Z;
704  }else{
705  // The primary remains in final state (e.g. electro-nucleus )
706  final4mom = temp.GetDynamicParticle()->Get4Momentum();
707  final_A = track_A;
708  final_Z = track_Z;
709  // Expect that the target nucleus will have interacted,
710  // and its products, including recoil, will be included in secondaries.
711  }
712  }
713  if( nSec > 0 ) {
714  G4Track* sec;
715 
716  for (G4int i = 0; i < nSec; i++) {
717  sec = theTotalResult->GetSecondary(i);
718  final4mom += sec->GetDynamicParticle()->Get4Momentum();
719  final_A += sec->GetDefinition()->GetBaryonNumber();
720  final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
721  }
722  }
723 
724  // Get level-checking information (used to cut-off relative checks)
725  G4String processName = GetProcessName();
727  G4String modelName("none");
728  if (theModel) modelName = theModel->GetModelName();
729  std::pair<G4double, G4double> checkLevels = epCheckLevels;
730  if (!levelsSetByProcess) {
731  if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
732  checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
733  checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
734  }
735 
736  // Compute absolute total-energy difference, and relative kinetic-energy
737  G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
738 
739  G4LorentzVector diff = initial4mom - final4mom;
740  G4double absolute = diff.e();
741  G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
742 
743  G4double absolute_mom = diff.vect().mag();
744  G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
745 
746  // Evaluate relative and absolute conservation
747  G4bool relPass = true;
748  G4String relResult = "pass";
749  if ( std::abs(relative) > checkLevels.first
750  || std::abs(relative_mom) > checkLevels.first) {
751  relPass = false;
752  relResult = checkRelative ? "fail" : "N/A";
753  }
754 
755  G4bool absPass = true;
756  G4String absResult = "pass";
757  if ( std::abs(absolute) > checkLevels.second
758  || std::abs(absolute_mom) > checkLevels.second ) {
759  absPass = false ;
760  absResult = "fail";
761  }
762 
763  G4bool chargePass = true;
764  G4String chargeResult = "pass";
765  if ( (initial_A-final_A)!=0
766  || (initial_Z-final_Z)!=0 ) {
767  chargePass = checkLevels.second < DBL_MAX ? false : true;
768  chargeResult = "fail";
769  }
770 
771  G4bool conservationPass = (relPass || absPass) && chargePass;
772 
773  std::stringstream Myout;
774  G4bool Myout_notempty(false);
775  // Options for level of reporting detail:
776  // 0. off
777  // 1. report only when E/p not conserved
778  // 2. report regardless of E/p conservation
779  // 3. report only when E/p not conserved, with model names, process names, and limits
780  // 4. report regardless of E/p conservation, with model names, process names, and limits
781  // negative -1.., as above, but send output to stderr
782 
783  if( std::abs(epReportLevel) == 4
784  || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
785  Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
786  Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
787  << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
788  << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
789  << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
790  << aNucleus.GetA_asInt() << ")" << G4endl;
791  Myout_notempty=true;
792  }
793  if ( std::abs(epReportLevel) == 4
794  || std::abs(epReportLevel) == 2
795  || ! conservationPass ){
796 
797  Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
798  << relative << " p/p(0)= " << relative_mom << G4endl;
799  Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
800  << absolute/MeV << " / " << absolute_mom/MeV << " 3mom: " << (diff.vect())*1./MeV << G4endl;
801  Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
802  Myout_notempty=true;
803 
804  }
805  Myout.flush();
806  if ( Myout_notempty ) {
807  if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
808  else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
809  }
810 }
811 
812 
814  const G4String& method,
816 {
817  ed << "Unrecoverable error in the method " << method << " of "
818  << GetProcessName() << G4endl;
819  ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
820  << aTrack.GetParentID()
821  << " " << aTrack.GetParticleDefinition()->GetParticleName()
822  << G4endl;
823  ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
824  << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
825  ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
826 
827  if (aTrack.GetMaterial()) {
828  ed << " material " << aTrack.GetMaterial()->GetName();
829  }
830  ed << G4endl;
831 
832  if (aTrack.GetVolume()) {
833  ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
834  << ">" << G4endl;
835  }
836 }
G4double G4ParticleHPJENDLHEData::G4double result
G4ParticleDefinition * GetDefinition() const
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4int GetParentID() const
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4int GetNumberOfSecondaries() const
void RegisterMe(G4HadronicInteraction *a)
G4Track * GetSecondary(G4int anIndex) const
void RegisterInteraction(G4HadronicProcess *, G4HadronicInteraction *)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4HadSecondary * GetSecondary(size_t i)
G4double GetTotalEnergy() const
G4LorentzRotation & GetTrafoToLab()
void BiasCrossSectionByFactor(G4double aScale)
static G4HadronicProcessStore * Instance()
const G4DynamicParticle * GetDynamicParticle() const
static G4KaonZero * Definition()
Definition: G4KaonZero.cc:53
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4int first(char) const
static G4AntiKaonZero * Definition()
const char * p
Definition: xmltok.h:285
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
const G4String & GetName() const
Definition: G4Material.hh:178
G4double GetElementCrossSection(const G4DynamicParticle *part, const G4Element *elm, const G4Material *mat=nullptr)
const G4ThreeVector & GetMomentumChange() const
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4HadProjectile thePro
G4double GetEnergyChange() const
G4ParticleDefinition * GetDefinition() const
tuple x
Definition: test.py:50
const G4ThreeVector * GetMomentumDirection() const
void CheckEnergyMomentumConservation(const G4Track &, const G4Nucleus &)
const G4String & GetModelName() const
G4double GetTime() const
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4HadronicProcessType
const G4String & GetParticleName() const
G4int GetCreatorModelType() const
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetCreatorModelIndex(G4int idx)
void RegisterMe(G4HadronicInteraction *a)
void ProposeWeight(G4double finalWeight)
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
const G4LorentzRotation & GetTrafoToLab() const
void SetSecondaryWeightByProcess(G4bool)
void DumpState(const G4Track &, const G4String &, G4ExceptionDescription &)
void SetTrafoToLab(const G4LorentzRotation &aT)
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
G4double GetKineticEnergy() const
void FillResult(G4HadFinalState *aR, const G4Track &aT)
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4double GetTotalNumberOfInteractionLengthTraversed() const
Definition: G4VProcess.hh:458
G4ParticleChange * theTotalResult
const G4ParticleDefinition * GetDefinition() const
void Register(G4HadronicProcess *)
const G4String & GetName() const
G4double GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
static constexpr double mm
Definition: SystemOfUnits.h:95
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
G4double GetKineticEnergy() const
G4double GetCrossSection(const G4DynamicParticle *, const G4Material *)
const G4ParticleDefinition * GetParticleDefinition() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
const G4int nmax
Definition: G4Step.hh:76
G4int GetTrackID() const
const G4ParticleDefinition * GetParticleDefinition() const
virtual void ProcessDescription(std::ostream &outFile) const
G4double GetGlobalTime() const
static G4KaonZeroLong * Definition()
static const char * G4Hadronic_Random_File
void RegisterParticle(G4HadronicProcess *, const G4ParticleDefinition *)
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
const G4LorentzVector & Get4Momentum() const
void DeRegister(G4HadronicProcess *)
HepLorentzVector & rotate(double, const Hep3Vector &)
const G4TouchableHandle & GetTouchableHandle() const
G4LorentzVector Get4Momentum() const
G4Material * GetMaterial() const
void Initialise(const G4Track &aT)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void Set4Momentum(const G4LorentzVector &momentum)
G4int size() const
void BuildPhysicsTable(const G4ParticleDefinition &)
G4ThreeVector GetMomentum() const
virtual void Initialize(const G4Track &)
static constexpr double GeV
const G4ThreeVector & GetMomentumDirection() const
G4HadronicInteraction * GetHadronicInteraction() const
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ProcessManager * GetProcessManager() const
void SetNumberOfSecondaries(G4int totSecondaries)
static void saveEngineStatus(const char filename[]="Config.conf")
Definition: Random.cc:275
Hep3Vector unit() const
G4DynamicParticle * GetParticle()
void BuildPhysicsTable(const G4ParticleDefinition &)
G4int GetZasInt() const
Definition: G4Element.hh:132
static G4KaonZeroShort * Definition()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void ProposeEnergy(G4double finalEnergy)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VPhysicalVolume * GetVolume() const
void AddSecondary(G4Track *aSecondary)
G4double GetWeight() const
G4HadronicInteraction * ChooseHadronicInteraction(const G4HadProjectile &aHadProjectile, G4Nucleus &aTargetNucleus, G4Material *aMaterial, G4Element *anElement)
G4double GetEnergy() const
#define G4endl
Definition: G4ios.hh:61
G4HadronicProcess(const G4String &processName="Hadronic", G4ProcessType procType=fHadronic)
static constexpr double MeV
Definition: G4SIunits.hh:214
G4TrackStatus GetTrackStatus() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4Element * SampleZandA(const G4DynamicParticle *, const G4Material *, G4Nucleus &target)
void SetKineticEnergy(const G4double aValue)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
const G4String & GetName() const
Definition: G4Element.hh:127
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4ForceCondition
G4double GetPDGCharge() const
double mag() const
#define DBL_MAX
Definition: templates.hh:83
void SetMass(G4double mass)
void Report(std::ostream &aS)
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4Material * FindSimpleMaterial(G4int Z) const
G4double GetLocalEnergyDeposit() const
void SetMomentumDirection(const G4ThreeVector &aValue)
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result)
G4GLOB_DLL std::ostream G4cerr
G4HadFinalStateStatus GetStatusChange() const
void PrintInfo(const G4ParticleDefinition *)
G4double GetTotalEnergy() const
G4ProcessType
G4double GetWeight() const