Geant4_10
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 67989 2013-03-13 10:54:03Z 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 <typeinfo>
72 #include <sstream>
73 #include <iostream>
74 
75 #include <stdlib.h>
76 
77 // File-scope variable to capture environment variable at startup
78 
79 static const char* G4Hadronic_Random_File = getenv("G4HADRONIC_RANDOM_FILE");
80 
82 
84  G4ProcessType procType)
85  : G4VDiscreteProcess(processName, procType)
86 {
87  SetProcessSubType(fHadronInelastic); // Default unless subclass changes
88 
91  theInteraction = 0;
92  theCrossSectionDataStore = new G4CrossSectionDataStore();
94  aScaleFactor = 1;
95  xBiasOn = false;
96  G4HadronicProcess_debug_flag = false;
97 
98  GetEnergyMomentumCheckEnvvars();
99 }
100 
102 
104  G4HadronicProcessType aHadSubType)
105  : G4VDiscreteProcess(processName, fHadronic)
106 {
107  SetProcessSubType(aHadSubType);
108 
111  theInteraction = 0;
112  theCrossSectionDataStore = new G4CrossSectionDataStore();
114  aScaleFactor = 1;
115  xBiasOn = false;
116  G4HadronicProcess_debug_flag = false;
117 
118  GetEnergyMomentumCheckEnvvars();
119 }
120 
121 
123 {
125  delete theTotalResult;
126  delete theCrossSectionDataStore;
127 }
128 
129 void G4HadronicProcess::GetEnergyMomentumCheckEnvvars() {
130  levelsSetByProcess = false;
131 
132  epReportLevel = getenv("G4Hadronic_epReportLevel") ?
133  strtol(getenv("G4Hadronic_epReportLevel"),0,10) : 0;
134 
135  epCheckLevels.first = getenv("G4Hadronic_epCheckRelativeLevel") ?
136  strtod(getenv("G4Hadronic_epCheckRelativeLevel"),0) : DBL_MAX;
137 
138  epCheckLevels.second = getenv("G4Hadronic_epCheckAbsoluteLevel") ?
139  strtod(getenv("G4Hadronic_epCheckAbsoluteLevel"),0) : DBL_MAX;
140 }
141 
143 {
144  if(!a) { return; }
145  try{GetManagerPointer()->RegisterMe( a );}
146  catch(G4HadronicException & aE)
147  {
149  aE.Report(ed);
150  ed << "Unrecoverable error in " << GetProcessName()
151  << " to register " << a->GetModelName() << G4endl;
152  G4Exception("G4HadronicProcess::RegisterMe", "had001", FatalException,
153  ed);
154  }
156 }
157 
159 {
160  if(getenv("G4HadronicProcess_debug")) {
161  G4HadronicProcess_debug_flag = true;
162  }
164 }
165 
167 {
168  try
169  {
170  theCrossSectionDataStore->BuildPhysicsTable(p);
171  }
172  catch(G4HadronicException aR)
173  {
175  aR.Report(ed);
176  ed << " hadronic initialisation fails" << G4endl;
177  G4Exception("G4HadronicProcess::BuildPhysicsTable", "had000",
178  FatalException,ed);
179  }
181 }
182 
185 {
186  try
187  {
188  theLastCrossSection = aScaleFactor*
189  theCrossSectionDataStore->GetCrossSection(aTrack.GetDynamicParticle(),
190  aTrack.GetMaterial());
191  }
192  catch(G4HadronicException aR)
193  {
195  aR.Report(ed);
196  DumpState(aTrack,"GetMeanFreePath",ed);
197  ed << " Cross section is not available" << G4endl;
198  G4Exception("G4HadronicProcess::GetMeanFreePath", "had002", FatalException,
199  ed);
200  }
201  G4double res = DBL_MAX;
202  if( theLastCrossSection > 0.0 ) { res = 1.0/theLastCrossSection; }
203  return res;
204 }
205 
208 {
209  // if primary is not Alive then do nothing
211  theTotalResult->Initialize(aTrack);
213  if(aTrack.GetTrackStatus() != fAlive) { return theTotalResult; }
214 
215  // Find cross section at end of step and check if <= 0
216  //
217  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
218  G4Material* aMaterial = aTrack.GetMaterial();
219 
220  G4Element* anElement = 0;
221  try
222  {
223  anElement = theCrossSectionDataStore->SampleZandA(aParticle,
224  aMaterial,
225  targetNucleus);
226  }
227  catch(G4HadronicException & aR)
228  {
230  aR.Report(ed);
231  DumpState(aTrack,"SampleZandA",ed);
232  ed << " PostStepDoIt failed on element selection" << G4endl;
233  G4Exception("G4HadronicProcess::PostStepDoIt", "had003", FatalException,
234  ed);
235  }
236 
237  // check only for charged particles
238  if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
239  if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) {
240  // No interaction
241  return theTotalResult;
242  }
243  }
244 
245  // Next check for illegal track status
246  //
247  if (aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend) {
248  if (aTrack.GetTrackStatus() == fStopAndKill ||
250  aTrack.GetTrackStatus() == fPostponeToNextEvent) {
252  ed << "G4HadronicProcess: track in unusable state - "
253  << aTrack.GetTrackStatus() << G4endl;
254  ed << "G4HadronicProcess: returning unchanged track " << G4endl;
255  DumpState(aTrack,"PostStepDoIt",ed);
256  G4Exception("G4HadronicProcess::PostStepDoIt", "had004", JustWarning, ed);
257  }
258  // No warning for fStopButAlive which is a legal status here
259  return theTotalResult;
260  }
261 
262  // Go on to regular case
263  //
264  G4double originalEnergy = aParticle->GetKineticEnergy();
265  G4double kineticEnergy = originalEnergy;
266 
267  // Get kinetic energy per nucleon for ions
268  if(aParticle->GetParticleDefinition()->GetBaryonNumber() > 1.5)
269  kineticEnergy/=aParticle->GetParticleDefinition()->GetBaryonNumber();
270 
271  try
272  {
273  theInteraction =
274  ChooseHadronicInteraction( kineticEnergy, aMaterial, anElement );
275  }
276  catch(G4HadronicException & aE)
277  {
279  aE.Report(ed);
280  ed << "Target element "<<anElement->GetName()<<" Z= "
281  << targetNucleus.GetZ_asInt() << " A= "
282  << targetNucleus.GetA_asInt() << G4endl;
283  DumpState(aTrack,"ChooseHadronicInteraction",ed);
284  ed << " No HadronicInteraction found out" << G4endl;
285  G4Exception("G4HadronicProcess::PostStepDoIt", "had005", FatalException,
286  ed);
287  }
288 
289  // Initialize the hadronic projectile from the track
290  thePro.Initialise(aTrack);
291  G4HadFinalState* result = 0;
292  G4int reentryCount = 0;
293 
294  do
295  {
296  try
297  {
298  // Save random engine if requested for debugging
299  if (G4Hadronic_Random_File) {
300  CLHEP::HepRandom::saveEngineStatus(G4Hadronic_Random_File);
301  }
302  // Call the interaction
303  result = theInteraction->ApplyYourself( thePro, targetNucleus);
304  ++reentryCount;
305  }
306  catch(G4HadronicException aR)
307  {
309  aR.Report(ed);
310  ed << "Call for " << theInteraction->GetModelName() << G4endl;
311  ed << "Target element "<<anElement->GetName()<<" Z= "
312  << targetNucleus.GetZ_asInt()
313  << " A= " << targetNucleus.GetA_asInt() << G4endl;
314  DumpState(aTrack,"ApplyYourself",ed);
315  ed << " ApplyYourself failed" << G4endl;
316  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
317  ed);
318  }
319 
320  // Check the result for catastrophic energy non-conservation
321  result = CheckResult(thePro,targetNucleus, result);
322 
323  if(reentryCount>100) {
325  ed << "Call for " << theInteraction->GetModelName() << G4endl;
326  ed << "Target element "<<anElement->GetName()<<" Z= "
327  << targetNucleus.GetZ_asInt()
328  << " A= " << targetNucleus.GetA_asInt() << G4endl;
329  DumpState(aTrack,"ApplyYourself",ed);
330  ed << " ApplyYourself does not completed after 100 attempts" << G4endl;
331  G4Exception("G4HadronicProcess::PostStepDoIt", "had006", FatalException,
332  ed);
333  }
334  }
335  while(!result);
336 
338 
340 
341  FillResult(result, aTrack);
342 
343  if (epReportLevel != 0) {
344  CheckEnergyMomentumConservation(aTrack, targetNucleus);
345  }
346  return theTotalResult;
347 }
348 
349 
351 {
352  outFile << "The description for this process has not been written yet.\n";
353 }
354 
355 
356 G4double G4HadronicProcess::XBiasSurvivalProbability()
357 {
358  G4double result = 0;
360  G4double biasedProbability = 1.-std::exp(-nLTraversed);
361  G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
362  result = (biasedProbability-realProbability)/biasedProbability;
363  return result;
364 }
365 
366 G4double G4HadronicProcess::XBiasSecondaryWeight()
367 {
368  G4double result = 0;
370  result =
371  1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
372  return result;
373 }
374 
375 void
377 {
379 
380  G4double rotation = CLHEP::twopi*G4UniformRand();
381  G4ThreeVector it(0., 0., 1.);
382 
383  G4double efinal = aR->GetEnergyChange();
384  if(efinal < 0.0) { efinal = 0.0; }
385 
386  // check status of primary
387  if(aR->GetStatusChange() == stopAndKill) {
390 
391  // check its final energy
392  } else if(0.0 == efinal) {
395  ->GetAtRestProcessVector()->size() > 0)
398 
399  // primary is not killed apply rotation and Lorentz transformation
400  } else {
403  G4double newE = efinal + mass;
404  G4double newP = std::sqrt(efinal*(efinal + 2*mass));
405  G4ThreeVector newPV = newP*aR->GetMomentumChange();
406  G4LorentzVector newP4(newE, newPV);
407  newP4.rotate(rotation, it);
408  newP4 *= aR->GetTrafoToLab();
410  newE = newP4.e() - mass;
411  if(G4HadronicProcess_debug_flag && newE <= 0.0) {
413  DumpState(aT,"Primary has zero energy after interaction",ed);
414  G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning, ed);
415  }
416  if(newE < 0.0) { newE = 0.0; }
417  theTotalResult->ProposeEnergy( newE );
418  }
419 
420  // check secondaries: apply rotation and Lorentz transformation
421  G4int nSec = aR->GetNumberOfSecondaries();
423  G4double weight = aT.GetWeight();
424 
425  if (nSec > 0) {
426  G4double time0 = aT.GetGlobalTime();
427  for (G4int i = 0; i < nSec; ++i) {
429  theM.rotate(rotation, it);
430  theM *= aR->GetTrafoToLab();
431  aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
432 
433  // time of interaction starts from zero
434  G4double time = aR->GetSecondary(i)->GetTime();
435  if (time < 0.0) { time = 0.0; }
436 
437  // take into account global time
438  time += time0;
439 
440  G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
441  time, aT.GetPosition());
442  G4double newWeight = weight*aR->GetSecondary(i)->GetWeight();
443  // G4cout << "#### ParticleDebug "
444  // <<GetProcessName()<<" "
445  // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
446  // <<aScaleFactor<<" "
447  // <<XBiasSurvivalProbability()<<" "
448  // <<XBiasSecondaryWeight()<<" "
449  // <<aT.GetWeight()<<" "
450  // <<aR->GetSecondary(i)->GetWeight()<<" "
451  // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
452  // <<G4endl;
453  track->SetWeight(newWeight);
456  if (G4HadronicProcess_debug_flag) {
457  G4double e = track->GetKineticEnergy();
458  if (e <= 0.0) {
460  DumpState(aT,"Secondary has zero energy",ed);
461  ed << "Secondary " << track->GetDefinition()->GetParticleName()
462  << G4endl;
463  G4Exception("G4HadronicProcess::FillResults", "had011", JustWarning,ed);
464  }
465  }
466  }
467  }
468 
469  aR->Clear();
470  return;
471 }
472 /*
473 void
474 G4HadronicProcess::FillTotalResult(G4HadFinalState* aR, const G4Track& aT)
475 {
476  theTotalResult->Clear();
477  theTotalResult->ProposeLocalEnergyDeposit(0.);
478  theTotalResult->Initialize(aT);
479  theTotalResult->SetSecondaryWeightByProcess(true);
480  theTotalResult->ProposeTrackStatus(fAlive);
481  G4double rotation = CLHEP::twopi*G4UniformRand();
482  G4ThreeVector it(0., 0., 1.);
483 
484  if(aR->GetStatusChange()==stopAndKill)
485  {
486  if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
487  {
488  theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
489  }
490  else
491  {
492  theTotalResult->ProposeTrackStatus(fStopAndKill);
493  theTotalResult->ProposeEnergy( 0.0 );
494  }
495  }
496  else if(aR->GetStatusChange()!=stopAndKill )
497  {
498  if(aR->GetStatusChange()==suspend)
499  {
500  theTotalResult->ProposeTrackStatus(fSuspend);
501  if(xBiasOn)
502  {
503  G4ExceptionDescription ed;
504  DumpState(aT,"FillTotalResult",ed);
505  G4Exception("G4HadronicProcess::FillTotalResult", "had007", FatalException,
506  ed,"Cannot cross-section bias a process that suspends tracks.");
507  }
508  } else if (aT.GetKineticEnergy() == 0) {
509  theTotalResult->ProposeTrackStatus(fStopButAlive);
510  }
511 
512  if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
513  {
514  theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
515  G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
516  G4double newM=aT.GetParticleDefinition()->GetPDGMass();
517  G4double newE=aR->GetEnergyChange() + newM;
518  G4double newP=std::sqrt(newE*newE - newM*newM);
519  G4DynamicParticle * aNew =
520  new G4DynamicParticle(aT.GetParticleDefinition(), newE, newP*aR->GetMomentumChange());
521  aR->AddSecondary(G4HadSecondary(aNew, newWeight));
522  }
523  else
524  {
525  G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
526  theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
527  if(aR->GetEnergyChange()>-.5)
528  {
529  theTotalResult->ProposeEnergy(aR->GetEnergyChange());
530  }
531  G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
532  newDirection*=aR->GetTrafoToLab();
533  theTotalResult->ProposeMomentumDirection(newDirection.vect());
534  }
535  }
536  else
537  {
538  G4ExceptionDescription ed;
539  ed << "Call for " << theInteraction->GetModelName() << G4endl;
540  ed << "Target Z= "
541  << targetNucleus.GetZ_asInt()
542  << " A= " << targetNucleus.GetA_asInt() << G4endl;
543  DumpState(aT,"FillTotalResult",ed);
544  G4Exception("G4HadronicProcess", "had008", FatalException,
545  "use of unsupported track-status.");
546  }
547 
548  if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
549  && theTotalResult->GetTrackStatus()==fAlive
550  && aR->GetStatusChange()==isAlive)
551  {
552  // Use for debugging: G4double newWeight = theTotalResult->GetParentWeight();
553 
554  G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
555  G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetParticleDefinition(),
556  aR->GetMomentumChange(),
557  newKE);
558  aR->AddSecondary(aNew);
559  aR->SetStatusChange(stopAndKill);
560 
561  theTotalResult->ProposeTrackStatus(fStopAndKill);
562  theTotalResult->ProposeEnergy( 0.0 );
563 
564  }
565  theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
566  theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
567 
568  if(aR->GetStatusChange() != stopAndKill)
569  {
570  G4double newM=aT.GetParticleDefinition()->GetPDGMass();
571  G4double newE=aR->GetEnergyChange() + newM;
572  G4double newP=std::sqrt(newE*newE - newM*newM);
573  G4ThreeVector newPV = newP*aR->GetMomentumChange();
574  G4LorentzVector newP4(newE, newPV);
575  newP4.rotate(rotation, it);
576  newP4*=aR->GetTrafoToLab();
577  theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
578  }
579 
580  for(G4int i=0; i<aR->GetNumberOfSecondaries(); ++i)
581  {
582  G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
583  theM.rotate(rotation, it);
584  theM*=aR->GetTrafoToLab();
585  aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
586  G4double time = aR->GetSecondary(i)->GetTime();
587  if(time<0) time = aT.GetGlobalTime();
588 
589  G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
590  time,
591  aT.GetPosition());
592 
593  G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
594  if(xBiasOn) { newWeight *= XBiasSecondaryWeight(); }
595  track->SetWeight(newWeight);
596  track->SetTouchableHandle(aT.GetTouchableHandle());
597  theTotalResult->AddSecondary(track);
598  }
599 
600  aR->Clear();
601  return;
602 }
603 */
604 
606 {
607  xBiasOn = true;
608  aScaleFactor = aScale;
609  G4String it = GetProcessName();
610  if( (it != "PhotonInelastic") &&
611  (it != "ElectroNuclear") &&
612  (it != "PositronNuclear") )
613  {
615  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had009", FatalException, ed,
616  "Cross-section biasing available only for gamma and electro nuclear reactions.");
617  }
618  if(aScale<100)
619  {
621  G4Exception("G4HadronicProcess::BiasCrossSectionByFactor", "had010", JustWarning,ed,
622  "Cross-section bias readjusted to be above safe limit. New value is 100");
623  aScaleFactor = 100.;
624  }
625 }
626 
628 {
629  // check for catastrophic energy non-conservation, to re-sample the interaction
630 
632  G4double nuclearMass(0);
633  if (theModel){
634 
635  // Compute final-state total energy
636  G4double finalE(0.);
637  G4int nSec = result->GetNumberOfSecondaries();
638 
639  nuclearMass = G4NucleiProperties::GetNuclearMass(aNucleus.GetA_asInt(),
640  aNucleus.GetZ_asInt());
641  if (result->GetStatusChange() != stopAndKill) {
642  // Interaction didn't complete, returned "do nothing" state => reset nucleus
643  // or the primary survived the interaction (e.g. electro-nuclear ) => keep nucleus
644  finalE=result->GetLocalEnergyDeposit() +
645  aPro.GetDefinition()->GetPDGMass() + result->GetEnergyChange();
646  if( nSec == 0 ){
647  // Since there are no secondaries, there is no recoil nucleus.
648  // To check energy balance we must neglect the initial nucleus too.
649  nuclearMass=0.0;
650  }
651  }
652  for (G4int i = 0; i < nSec; i++) {
653  finalE += result->GetSecondary(i)->GetParticle()->GetTotalEnergy();
654  }
655  G4double deltaE= nuclearMass + aPro.GetTotalEnergy() - finalE;
656 
657  std::pair<G4double, G4double> checkLevels = theModel->GetFatalEnergyCheckLevels(); // (relative, absolute)
658  if (std::abs(deltaE) > checkLevels.second && std::abs(deltaE) > checkLevels.first*aPro.GetKineticEnergy()){
659  // do not delete result, this is a pointer to a data member;
660  result=0;
662  desc << "Warning: Bad energy non-conservation detected, will "
663  << (epReportLevel<0 ? "abort the event" : "re-sample the interaction") << G4endl
664  << " Process / Model: " << GetProcessName()<< " / " << theModel->GetModelName() << G4endl
665  << " Primary: " << aPro.GetDefinition()->GetParticleName()
666  << " (" << aPro.GetDefinition()->GetPDGEncoding() << "),"
667  << " E= " << aPro.Get4Momentum().e()
668  << ", target nucleus (" << aNucleus.GetZ_asInt() << ","<< aNucleus.GetA_asInt() << ")" << G4endl
669  << " E(initial - final) = " << deltaE << " MeV." << G4endl;
670  G4Exception("G4HadronicProcess:CheckResult()", "had012", epReportLevel<0 ? EventMustBeAborted : JustWarning,desc);
671  }
672  }
673  return result;
674 }
675 
676 void
678  const G4Nucleus& aNucleus)
679 {
680  G4int target_A=aNucleus.GetA_asInt();
681  G4int target_Z=aNucleus.GetZ_asInt();
682  G4double targetMass = G4NucleiProperties::GetNuclearMass(target_A,target_Z);
683  G4LorentzVector target4mom(0, 0, 0, targetMass);
684 
685  G4LorentzVector projectile4mom = aTrack.GetDynamicParticle()->Get4Momentum();
686  G4int track_A = aTrack.GetDefinition()->GetBaryonNumber();
687  G4int track_Z = G4lrint(aTrack.GetDefinition()->GetPDGCharge());
688 
689  G4int initial_A = target_A + track_A;
690  G4int initial_Z = target_Z + track_Z;
691 
692  G4LorentzVector initial4mom = projectile4mom + target4mom;
693 
694  // Compute final-state momentum for scattering and "do nothing" results
695  G4LorentzVector final4mom;
696  G4int final_A(0), final_Z(0);
697 
699  if (theTotalResult->GetTrackStatus() != fStopAndKill) { // If it is Alive
700  // Either interaction didn't complete, returned "do nothing" state
701  // or the primary survived the interaction (e.g. electro-nucleus )
702  G4Track temp(aTrack);
703 
704  // Use the final energy / momentum
707 
708  if( nSec == 0 ){
709  // Interaction didn't complete, returned "do nothing" state
710  // - or suppressed recoil (e.g. Neutron elastic )
711  final4mom = temp.GetDynamicParticle()->Get4Momentum() + target4mom;
712  final_A = initial_A;
713  final_Z = initial_Z;
714  }else{
715  // The primary remains in final state (e.g. electro-nucleus )
716  final4mom = temp.GetDynamicParticle()->Get4Momentum();
717  final_A = track_A;
718  final_Z = track_Z;
719  // Expect that the target nucleus will have interacted,
720  // and its products, including recoil, will be included in secondaries.
721  }
722  }
723  if( nSec > 0 ) {
724  G4Track* sec;
725 
726  for (G4int i = 0; i < nSec; i++) {
727  sec = theTotalResult->GetSecondary(i);
728  final4mom += sec->GetDynamicParticle()->Get4Momentum();
729  final_A += sec->GetDefinition()->GetBaryonNumber();
730  final_Z += G4lrint(sec->GetDefinition()->GetPDGCharge());
731  }
732  }
733 
734  // Get level-checking information (used to cut-off relative checks)
735  G4String processName = GetProcessName();
737  G4String modelName("none");
738  if (theModel) modelName = theModel->GetModelName();
739  std::pair<G4double, G4double> checkLevels = epCheckLevels;
740  if (!levelsSetByProcess) {
741  if (theModel) checkLevels = theModel->GetEnergyMomentumCheckLevels();
742  checkLevels.first= std::min(checkLevels.first, epCheckLevels.first);
743  checkLevels.second=std::min(checkLevels.second, epCheckLevels.second);
744  }
745 
746  // Compute absolute total-energy difference, and relative kinetic-energy
747  G4bool checkRelative = (aTrack.GetKineticEnergy() > checkLevels.second);
748 
749  G4LorentzVector diff = initial4mom - final4mom;
750  G4double absolute = diff.e();
751  G4double relative = checkRelative ? absolute/aTrack.GetKineticEnergy() : 0.;
752 
753  G4double absolute_mom = diff.vect().mag();
754  G4double relative_mom = checkRelative ? absolute_mom/aTrack.GetMomentum().mag() : 0.;
755 
756  // Evaluate relative and absolute conservation
757  G4bool relPass = true;
758  G4String relResult = "pass";
759  if ( std::abs(relative) > checkLevels.first
760  || std::abs(relative_mom) > checkLevels.first) {
761  relPass = false;
762  relResult = checkRelative ? "fail" : "N/A";
763  }
764 
765  G4bool absPass = true;
766  G4String absResult = "pass";
767  if ( std::abs(absolute) > checkLevels.second
768  || std::abs(absolute_mom) > checkLevels.second ) {
769  absPass = false ;
770  absResult = "fail";
771  }
772 
773  G4bool chargePass = true;
774  G4String chargeResult = "pass";
775  if ( (initial_A-final_A)!=0
776  || (initial_Z-final_Z)!=0 ) {
777  chargePass = checkLevels.second < DBL_MAX ? false : true;
778  chargeResult = "fail";
779  }
780 
781  G4bool conservationPass = (relPass || absPass) && chargePass;
782 
783  std::stringstream Myout;
784  G4bool Myout_notempty(false);
785  // Options for level of reporting detail:
786  // 0. off
787  // 1. report only when E/p not conserved
788  // 2. report regardless of E/p conservation
789  // 3. report only when E/p not conserved, with model names, process names, and limits
790  // 4. report regardless of E/p conservation, with model names, process names, and limits
791  // negative -1.., as above, but send output to stderr
792 
793  if( std::abs(epReportLevel) == 4
794  || ( std::abs(epReportLevel) == 3 && ! conservationPass ) ){
795  Myout << " Process: " << processName << " , Model: " << modelName << G4endl;
796  Myout << " Primary: " << aTrack.GetParticleDefinition()->GetParticleName()
797  << " (" << aTrack.GetParticleDefinition()->GetPDGEncoding() << "),"
798  << " E= " << aTrack.GetDynamicParticle()->Get4Momentum().e()
799  << ", target nucleus (" << aNucleus.GetZ_asInt() << ","
800  << aNucleus.GetA_asInt() << ")" << G4endl;
801  Myout_notempty=true;
802  }
803  if ( std::abs(epReportLevel) == 4
804  || std::abs(epReportLevel) == 2
805  || ! conservationPass ){
806 
807  Myout << " "<< relResult <<" relative, limit " << checkLevels.first << ", values E/T(0) = "
808  << relative << " p/p(0)= " << relative_mom << G4endl;
809  Myout << " "<< absResult << " absolute, limit (MeV) " << checkLevels.second/MeV << ", values E / p (MeV) = "
810  << absolute/MeV << " / " << absolute_mom/MeV << G4endl;
811  Myout << " "<< chargeResult << " charge/baryon number balance " << (initial_Z-final_Z) << " / " << (initial_A-final_A) << " "<< G4endl;
812  Myout_notempty=true;
813 
814  }
815  Myout.flush();
816  if ( Myout_notempty ) {
817  if (epReportLevel > 0) G4cout << Myout.str()<< G4endl;
818  else if (epReportLevel < 0) G4cerr << Myout.str()<< G4endl;
819  }
820 }
821 
822 
824  const G4String& method,
826 {
827  ed << "Unrecoverable error in the method " << method << " of "
828  << GetProcessName() << G4endl;
829  ed << "TrackID= "<< aTrack.GetTrackID() << " ParentID= "
830  << aTrack.GetParentID()
831  << " " << aTrack.GetParticleDefinition()->GetParticleName()
832  << G4endl;
833  ed << "Ekin(GeV)= " << aTrack.GetKineticEnergy()/CLHEP::GeV
834  << "; direction= " << aTrack.GetMomentumDirection() << G4endl;
835  ed << "Position(mm)= " << aTrack.GetPosition()/CLHEP::mm << ";";
836 
837  if (aTrack.GetMaterial()) {
838  ed << " material " << aTrack.GetMaterial()->GetName();
839  }
840  ed << G4endl;
841 
842  if (aTrack.GetVolume()) {
843  ed << "PhysicalVolume <" << aTrack.GetVolume()->GetName()
844  << ">" << G4endl;
845  }
846 }
847 /*
848 G4ParticleDefinition* G4HadronicProcess::GetTargetDefinition()
849 {
850  const G4Nucleus* nuc = GetTargetNucleus();
851  G4int Z = nuc->GetZ_asInt();
852  G4int A = nuc->GetA_asInt();
853  return G4ParticleTable::GetParticleTable()->GetIon(Z,A,0*eV);
854 }
855 */
856 /* end of file */
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 *)
tuple a
Definition: test.py:11
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4HadSecondary * GetSecondary(size_t i)
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4LorentzRotation & GetTrafoToLab()
void BiasCrossSectionByFactor(G4double aScale)
static G4HadronicProcessStore * Instance()
std::ofstream outFile
Definition: GammaRayTel.cc:68
const G4DynamicParticle * GetDynamicParticle() const
G4int first(char) const
const char * p
Definition: xmltok.h:285
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
const G4String & GetName() const
Definition: G4Material.hh:176
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
G4double G4NeutronHPJENDLHEData::G4double result
virtual void PreparePhysicsTable(const G4ParticleDefinition &)
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4HadProjectile thePro
G4double GetEnergyChange() const
G4ParticleDefinition * GetDefinition() const
const G4ThreeVector * GetMomentumDirection() const
void CheckEnergyMomentumConservation(const G4Track &, const G4Nucleus &)
const G4String & GetModelName() const
double weight
Definition: plottest35.C:25
G4double GetTime() const
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
G4ProcessManager * GetProcessManager() const
int G4int
Definition: G4Types.hh:78
G4HadronicProcessType
const G4String & GetParticleName() const
G4EnergyRangeManager * GetManagerPointer()
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void RegisterMe(G4HadronicInteraction *a)
void ProposeWeight(G4double finalWeight)
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
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:87
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 *)
bool G4bool
Definition: G4Types.hh:79
virtual std::pair< G4double, G4double > GetEnergyMomentumCheckLevels() const
G4double GetKineticEnergy() const
G4HadronicInteraction * ChooseHadronicInteraction(G4double kineticEnergy, G4Material *aMaterial, G4Element *anElement)
G4HadFinalState * CheckResult(const G4HadProjectile &thePro, const G4Nucleus &targetNucleus, G4HadFinalState *result) const
G4double GetCrossSection(const G4DynamicParticle *, const G4Material *)
const G4ParticleDefinition * GetParticleDefinition() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
G4int GetTrackID() const
const G4ParticleDefinition * GetParticleDefinition() const
virtual void ProcessDescription(std::ostream &outFile) const
G4double GetGlobalTime() const
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 &)
G4double GetElementCrossSection(const G4DynamicParticle *part, const G4Element *elm, const G4Material *mat=0)
const G4ThreeVector & GetMomentumDirection() const
G4HadronicInteraction * GetHadronicInteraction() const
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
const G4ThreeVector & GetMomentumChange() const
void SetNumberOfSecondaries(G4int totSecondaries)
static void saveEngineStatus(const char filename[]="Config.conf")
Definition: Random.cc:175
Hep3Vector unit() const
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
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
G4double GetEnergy() const
#define G4endl
Definition: G4ios.hh:61
G4HadronicProcess(const G4String &processName="Hadronic", G4ProcessType procType=fHadronic)
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
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
G4ForceCondition
G4double GetPDGCharge() const
double mag() const
#define DBL_MAX
Definition: templates.hh:83
void Report(std::ostream &aS)
G4int GetNumberOfSecondaries() const
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)=0
G4double GetLocalEnergyDeposit() const
void SetMomentumDirection(const G4ThreeVector &aValue)
G4GLOB_DLL std::ostream G4cerr
G4HadFinalStateStatus GetStatusChange() const
void PrintInfo(const G4ParticleDefinition *)
const G4LorentzRotation & GetTrafoToLab() const
G4double GetTotalEnergy() const
G4ProcessType
G4double GetWeight() const