Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/management/src/G4HadronicProcess.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/management/src/G4HadronicProcess.cc (Version 11.3.0) and /processes/hadronic/management/src/G4HadronicProcess.cc (Version 1.0)


                                                   >>   1 // This code implementation is the intellectual property of
                                                   >>   2 // the GEANT4 collaboration.
  1 //                                                  3 //
  2 // ******************************************* <<   4 // By copying, distributing or modifying the Program (or any work
  3 // * License and Disclaimer                    <<   5 // based on the Program) you indicate your acceptance of this statement,
  4 // *                                           <<   6 // and all its terms.
  5 // * The  Geant4 software  is  copyright of th <<   7 //
  6 // * the Geant4 Collaboration.  It is provided <<   8 // $Id: G4HadronicProcess.cc,v 1.7.2.1 1999/12/07 20:51:31 gunter Exp $
  7 // * conditions of the Geant4 Software License <<   9 // GEANT4 tag $Name: geant4-01-00 $
  8 // * LICENSE and available at  http://cern.ch/ <<  10 //
  9 // * include a list of copyright holders.      <<  11  // HPW to implement the choosing of an element for scattering.
 10 // *                                           <<  12 #include <fstream.h>
 11 // * Neither the authors of this software syst <<  13 #ifdef WIN32
 12 // * institutes,nor the agencies providing fin <<  14   #include <strstrea.h>
 13 // * work  make  any representation or  warran <<  15 #else
 14 // * regarding  this  software system or assum <<  16   #include <strstream.h>
 15 // * use.  Please see the license in the file  <<  17 #endif
 16 // * for the full disclaimer and the limitatio <<  18 #include <stdlib.h>
 17 // *                                           << 
 18 // * This  code  implementation is the result  << 
 19 // * technical work of the GEANT4 collaboratio << 
 20 // * By using,  copying,  modifying or  distri << 
 21 // * any work based  on the software)  you  ag << 
 22 // * use  in  resulting  scientific  publicati << 
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // ******************************************* << 
 25 //                                             << 
 26 //                                             << 
 27 // ------------------------------------------- << 
 28 //                                             << 
 29 // GEANT4 Class source file                    << 
 30 //                                             << 
 31 // G4HadronicProcess                           << 
 32 //                                             << 
 33 // original by H.P.Wellisch                    << 
 34 // J.L. Chuma, TRIUMF, 10-Mar-1997             << 
 35 //                                             << 
 36 // Modifications:                              << 
 37 // 05-Jul-2010 V.Ivanchenko cleanup commented  << 
 38 // 20-Jul-2011 M.Kelsey -- null-pointer checks << 
 39 // 24-Sep-2011 M.Kelsey -- Use envvar G4HADRON << 
 40 //    engine state before each model call      << 
 41 // 18-Oct-2011 M.Kelsey -- Handle final-state  << 
 42 // 14-Mar-2012 G.Folger -- enhance checks for  << 
 43 // 28-Jul-2012 M.Maire  -- add function GetTar << 
 44 // 14-Sep-2012 Inherit from RestDiscrete, use  << 
 45 //    configure base-class                     << 
 46 // 28-Sep-2012 Restore inheritance from G4VDis << 
 47 //    changing, remove warning message from or << 
 48 // 21-Aug-2019 V.Ivanchenko leave try/catch on << 
 49                                                << 
 50 #include "G4HadronicProcess.hh"                    19 #include "G4HadronicProcess.hh"
 51                                                <<  20 //@@ add model name info, once typeinfo available #include <typeinfo.h>
 52 #include "G4Types.hh"                          <<  21  
 53 #include "G4SystemOfUnits.hh"                  <<  22  G4IsoParticleChange * G4HadronicProcess::theIsoResult = NULL;
 54 #include "G4HadProjectile.hh"                  <<  23  G4IsoParticleChange * G4HadronicProcess::theOldIsoResult = NULL;
 55 #include "G4ElementVector.hh"                  <<  24  G4bool G4HadronicProcess::isoIsEnabled = true;
 56 #include "G4Track.hh"                          <<  25  void G4HadronicProcess::EnableIsotopeProductionGlobally()  {isoIsEnabled = true;}
 57 #include "G4Step.hh"                           <<  26  void G4HadronicProcess::DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
 58 #include "G4Element.hh"                        <<  27  
 59 #include "G4ParticleChange.hh"                 <<  28  G4Element * G4HadronicProcess::ChooseAandZ(
 60 #include "G4ProcessVector.hh"                  <<  29   const G4DynamicParticle *aParticle, const G4Material *aMaterial )
 61 #include "G4ProcessManager.hh"                 << 
 62 #include "G4NucleiProperties.hh"               << 
 63                                                << 
 64 #include "G4HadronicException.hh"              << 
 65 #include "G4HadronicProcessStore.hh"           << 
 66 #include "G4HadronicParameters.hh"             << 
 67 #include "G4VCrossSectionDataSet.hh"           << 
 68                                                << 
 69 #include "G4NistManager.hh"                    << 
 70 #include "G4VLeadingParticleBiasing.hh"        << 
 71 #include "G4HadXSHelper.hh"                    << 
 72 #include "G4Threading.hh"                      << 
 73 #include "G4Exp.hh"                            << 
 74                                                << 
 75 #include <typeinfo>                            << 
 76 #include <sstream>                             << 
 77 #include <iostream>                            << 
 78                                                << 
 79 namespace                                      << 
 80 {                                              << 
 81   constexpr G4double lambdaFactor = 0.8;       << 
 82   constexpr G4double invLambdaFactor = 1.0/lam << 
 83 }                                              << 
 84                                                << 
 85 ////////////////////////////////////////////// << 
 86                                                << 
 87 G4HadronicProcess::G4HadronicProcess(const G4S << 
 88                                      G4Process << 
 89  : G4VDiscreteProcess(processName, procType)   << 
 90 {                                              << 
 91   SetProcessSubType(fHadronInelastic);  // Def << 
 92   InitialiseLocal();                           << 
 93 }                                              << 
 94                                                << 
 95 G4HadronicProcess::G4HadronicProcess(const G4S << 
 96                                      G4Hadroni << 
 97  : G4VDiscreteProcess(processName, fHadronic)  << 
 98 {                                              << 
 99   SetProcessSubType(aHadSubType);              << 
100   InitialiseLocal();                           << 
101 }                                              << 
102                                                << 
103 G4HadronicProcess::~G4HadronicProcess()        << 
104 {                                              << 
105   theProcessStore->DeRegister(this);           << 
106   delete theTotalResult;                       << 
107   delete theCrossSectionDataStore;             << 
108   if(isMaster) {                               << 
109     if (fXSpeaks != nullptr) {                 << 
110       for (auto const& e : *fXSpeaks ) {       << 
111         delete e;                              << 
112       }                                        << 
113     }                                          << 
114     delete fXSpeaks;                           << 
115     delete theEnergyOfCrossSectionMax;         << 
116   }                                            << 
117 }                                              << 
118                                                << 
119 void G4HadronicProcess::InitialiseLocal() {    << 
120   theTotalResult = new G4ParticleChange();     << 
121   theTotalResult->SetSecondaryWeightByProcess( << 
122   theCrossSectionDataStore = new G4CrossSectio << 
123   theProcessStore = G4HadronicProcessStore::In << 
124   theProcessStore->Register(this);             << 
125   minKinEnergy = 1*CLHEP::MeV;                 << 
126                                                << 
127   G4HadronicParameters* param = G4HadronicPara << 
128   epReportLevel = param->GetEPReportLevel();   << 
129   epCheckLevels.first = param->GetEPRelativeLe << 
130   epCheckLevels.second = param->GetEPAbsoluteL << 
131                                                << 
132   unitVector.set(0.0, 0.0, 0.1);               << 
133   if(G4Threading::IsWorkerThread()) { isMaster << 
134 }                                              << 
135                                                << 
136 void G4HadronicProcess::RegisterMe( G4Hadronic << 
137 {                                              << 
138   if(nullptr == a) { return; }                 << 
139   theEnergyRangeManager.RegisterMe( a );       << 
140   G4HadronicProcessStore::Instance()->Register << 
141 }                                              << 
142                                                << 
143 G4double                                       << 
144 G4HadronicProcess::GetElementCrossSection(cons << 
145             const G4Element * elm,             << 
146             const G4Material* mat)             << 
147 {                                              << 
148   if(nullptr == mat)                           << 
149   {                                                30   {
150     static const G4int nmax = 5;               <<  31     currentZ = 0;
151     if(nMatWarn < nmax) {                      <<  32     currentN = 0;
152       ++nMatWarn;                              <<  33     const G4int numberOfElements = aMaterial->GetNumberOfElements();
153       G4ExceptionDescription ed;               <<  34     const G4ElementVector *theElementVector = aMaterial->GetElementVector();
154       ed << "Cannot compute Element x-section  <<  35     
155    << " because no material defined \n"        <<  36     if( numberOfElements == 1 ) 
156    << " Please, specify material pointer or de <<  37     {
157    << " for Z= " << elm->GetZasInt();          <<  38       currentZ = G4double((*theElementVector)(0)->GetZ());
158       G4Exception("G4HadronicProcess::GetEleme <<  39       currentN = (*theElementVector)(0)->GetN();
159       JustWarning, ed);                        <<  40       targetNucleus.SetParameters(currentN, currentZ);
160     }                                          <<  41       return (*theElementVector)(0);
161   }                                            <<  42     }
162   return theCrossSectionDataStore->GetCrossSec <<  43     
163 }                                              <<  44     const G4double *theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
164                                                <<  45     G4double crossSectionTotal = 0;
165 void G4HadronicProcess::PreparePhysicsTable(co <<  46     G4int i;
166 {                                              <<  47     for( i=0; i < numberOfElements; ++i )
167   if(nullptr == firstParticle) { firstParticle <<  48       crossSectionTotal += theAtomicNumberDensity[i] *
168   theProcessStore->RegisterParticle(this, &p); <<  49         dispatch->GetMicroscopicCrossSection( aParticle, (*theElementVector)(i) );
169 }                                              <<  50     
170                                                <<  51     G4double crossSectionSum= 0.;
171 void G4HadronicProcess::BuildPhysicsTable(cons <<  52     G4double random = G4UniformRand()*crossSectionTotal;
172 {                                              <<  53     for( i=0; i < numberOfElements; ++i )
173   if(firstParticle != &p) { return; }          <<  54     { 
174                                                <<  55       crossSectionSum += theAtomicNumberDensity[i] *
175   theCrossSectionDataStore->BuildPhysicsTable( <<  56         dispatch->GetMicroscopicCrossSection( aParticle, (*theElementVector)(i) );
176   theEnergyRangeManager.BuildPhysicsTable(p);  <<  57       if( random<=crossSectionSum )
177   G4HadronicParameters* param = G4HadronicPara <<  58       {
178                                                <<  59         currentZ = G4double((*theElementVector)(i)->GetZ());
179   G4int subtype = GetProcessSubType();         <<  60         currentN = (*theElementVector)(i)->GetN();
180   if(useIntegralXS) {                          <<  61         targetNucleus.SetParameters(currentN, currentZ);
181     if(subtype == fHadronInelastic) {          <<  62         return (*theElementVector)(i);
182       useIntegralXS = param->EnableIntegralIne <<  63       }
183     } else if(subtype == fHadronElastic) {     <<  64     }
184       useIntegralXS = param->EnableIntegralEla <<  65     currentZ = G4double((*theElementVector)(numberOfElements-1)->GetZ());
185     }                                          <<  66     currentN = (*theElementVector)(numberOfElements-1)->GetN();
186   }                                            <<  67     targetNucleus.SetParameters(currentN, currentZ);
187   fXSType = fHadNoIntegral;                    <<  68     return (*theElementVector)(numberOfElements-1);
188                                                << 
189   if(nullptr == masterProcess) {               << 
190     masterProcess = dynamic_cast<const G4Hadro << 
191   }                                            << 
192   if(nullptr == masterProcess) {               << 
193     if(1 < param->GetVerboseLevel()) {         << 
194       G4ExceptionDescription ed;               << 
195       ed << "G4HadronicProcess::BuildPhysicsTa << 
196    << GetProcessName() << " for " << p.GetPart << 
197    << " fail due to undefined pointer to the m << 
198    << "  ThreadID= " << G4Threading::G4GetThre << 
199    << "  initialisation of worker started befo << 
200       G4Exception("G4HadronicProcess::BuildPhy << 
201       JustWarning, ed);                        << 
202     }                                          << 
203   }                                                69   }
204                                                <<  70  
205   // check particle for integral method        <<  71  G4VParticleChange *G4HadronicProcess::GeneralPostStepDoIt(
206   if(isMaster || nullptr == masterProcess) {   <<  72   const G4Track &aTrack, const G4Step &aStep )
207     G4double charge = p.GetPDGCharge()/eplus;  <<  73   {
208                                                <<  74     const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
209     // select cross section shape              <<  75     G4Material *aMaterial = aTrack.GetMaterial();
210     if(charge != 0.0 && useIntegralXS) {       <<  76     G4double kineticEnergy = aParticle->GetKineticEnergy();
211       G4double tmax = param->GetMaxEnergy();   <<  77     G4Element * anElement = ChooseAandZ( aParticle, aMaterial );
212       currentParticle = firstParticle;         <<  78     theInteraction = ChooseHadronicInteraction( kineticEnergy,
213       // initialisation in the master thread   <<  79                                                           aMaterial, anElement );
214       G4int pdg = p.GetPDGEncoding();          <<  80     G4VParticleChange *result =
215       if (std::abs(pdg) == 211) {              <<  81       theInteraction->ApplyYourself( aTrack, targetNucleus);
216   fXSType = fHadTwoPeaks;                      <<  82     ResetNumberOfInteractionLengthLeft();
217       } else if (pdg == 321) {                 <<  83     if(isoIsOnAnyway!=-1)
218   fXSType = fHadOnePeak;                       <<  84     {
219       } else if (pdg == -321) {                <<  85       if(isoIsEnabled||isoIsOnAnyway)
220   fXSType = fHadDecreasing;                    <<  86       {
221       } else if (pdg == 2212) {                <<  87         result = DoIsotopeCounting(result, aTrack, targetNucleus);
222   fXSType = fHadTwoPeaks;                      << 
223       } else if (pdg == -2212 || pdg == -10000 << 
224      pdg == -1000020030 || pdg == -1000020040) << 
225   fXSType = fHadDecreasing;                    << 
226       } else if (charge > 0.0 || pdg == 11 ||  << 
227   fXSType = fHadIncreasing;                    << 
228       }                                        << 
229                                                << 
230       delete theEnergyOfCrossSectionMax;       << 
231       theEnergyOfCrossSectionMax = nullptr;    << 
232       if(fXSType == fHadTwoPeaks) {            << 
233         if (fXSpeaks != nullptr) {             << 
234           for (auto const& e : *fXSpeaks ) {   << 
235             delete e;                          << 
236           }                                    << 
237         }                                      << 
238   delete fXSpeaks;                             << 
239   fXSpeaks =                                   << 
240     G4HadXSHelper::FillPeaksStructure(this, &p << 
241   if(nullptr == fXSpeaks) {                    << 
242     fXSType = fHadOnePeak;                     << 
243   }                                            << 
244       }                                        << 
245       if(fXSType == fHadOnePeak) {             << 
246   theEnergyOfCrossSectionMax =                 << 
247     G4HadXSHelper::FindCrossSectionMax(this, & << 
248   if(nullptr == theEnergyOfCrossSectionMax) {  << 
249     fXSType = fHadIncreasing;                  << 
250   }                                            << 
251       }                                            88       }
252     }                                              89     }
253   } else {                                     <<  90     return result;
254     // initialisation in worker threads        << 
255     fXSType = masterProcess->CrossSectionType( << 
256     fXSpeaks = masterProcess->TwoPeaksXS();    << 
257     theEnergyOfCrossSectionMax = masterProcess << 
258   }                                            << 
259   if(isMaster && 1 < param->GetVerboseLevel()) << 
260     G4cout << "G4HadronicProcess::BuildPhysics << 
261      << GetProcessName() << " and " << p.GetPa << 
262      << " typeXS=" << fXSType << G4endl;       << 
263   }                                            << 
264   G4HadronicProcessStore::Instance()->PrintInf << 
265 }                                              << 
266                                                << 
267 void G4HadronicProcess::StartTracking(G4Track* << 
268 {                                              << 
269   currentMat = nullptr;                        << 
270   currentParticle = track->GetDefinition();    << 
271   fDynParticle = track->GetDynamicParticle();  << 
272   theNumberOfInteractionLengthLeft = -1.0;     << 
273 }                                              << 
274                                                << 
275 G4double G4HadronicProcess::PostStepGetPhysica << 
276                              const G4Track& tr << 
277            G4double previousStepSize,          << 
278                              G4ForceCondition* << 
279 {                                              << 
280   *condition = NotForced;                      << 
281                                                << 
282   const G4Material* mat = track.GetMaterial(); << 
283   if(mat != currentMat) {                      << 
284     currentMat = mat;                          << 
285     mfpKinEnergy = DBL_MAX;                    << 
286     matIdx = (G4int)track.GetMaterial()->GetIn << 
287   }                                            << 
288   UpdateCrossSectionAndMFP(track.GetKineticEne << 
289                                                << 
290   // zero cross section                        << 
291   if(theLastCrossSection <= 0.0) {             << 
292     theNumberOfInteractionLengthLeft = -1.0;   << 
293     currentInteractionLength = DBL_MAX;        << 
294     return DBL_MAX;                            << 
295   }                                            << 
296                                                << 
297   // non-zero cross section                    << 
298   if (theNumberOfInteractionLengthLeft < 0.0)  << 
299     theNumberOfInteractionLengthLeft = -G4Log( << 
300     theInitialNumberOfInteractionLength = theN << 
301   } else {                                     << 
302     theNumberOfInteractionLengthLeft -=        << 
303       previousStepSize/currentInteractionLengt << 
304     theNumberOfInteractionLengthLeft =         << 
305       std::max(theNumberOfInteractionLengthLef << 
306   }                                            << 
307   currentInteractionLength = theMFP;           << 
308   return theNumberOfInteractionLengthLeft*theM << 
309 }                                              << 
310                                                << 
311 G4double G4HadronicProcess::GetMeanFreePath(   << 
312                             const G4Track &aTr << 
313                             G4ForceCondition*) << 
314 {                                              << 
315   G4double xs = aScaleFactor*theCrossSectionDa << 
316      ->ComputeCrossSection(aTrack.GetDynamicPa << 
317   return (xs > 0.0) ? 1.0/xs : DBL_MAX;        << 
318 }                                              << 
319                                                << 
320 G4VParticleChange*                             << 
321 G4HadronicProcess::PostStepDoIt(const G4Track& << 
322 {                                              << 
323   theNumberOfInteractionLengthLeft = -1.0;     << 
324                                                << 
325   //G4cout << "PostStepDoIt " << aTrack.GetDef << 
326   //   << " Ekin= " << aTrack.GetKineticEnergy << 
327   // if primary is not Alive then do nothing   << 
328   theTotalResult->Clear();                     << 
329   theTotalResult->Initialize(aTrack);          << 
330   fWeight = aTrack.GetWeight();                << 
331   theTotalResult->ProposeWeight(fWeight);      << 
332   if(aTrack.GetTrackStatus() != fAlive) { retu << 
333                                                << 
334   // Find cross section at end of step and che << 
335   //                                           << 
336   const G4DynamicParticle* aParticle = aTrack. << 
337   const G4Material* aMaterial = aTrack.GetMate << 
338                                                << 
339   // check only for charged particles          << 
340   if(fXSType != fHadNoIntegral) {              << 
341     mfpKinEnergy = DBL_MAX;                    << 
342     G4double xs = aScaleFactor*                << 
343       theCrossSectionDataStore->ComputeCrossSe << 
344     //G4cout << "xs=" << xs << " xs0=" << theL << 
345     //     << "  " << aMaterial->GetName() <<  << 
346     if(xs < theLastCrossSection*G4UniformRand( << 
347       // No interaction                        << 
348       return theTotalResult;                   << 
349     }                                          << 
350   }                                            << 
351                                                << 
352   const G4Element* anElement =                 << 
353     theCrossSectionDataStore->SampleZandA(aPar << 
354                                                << 
355   // Next check for illegal track status       << 
356   //                                           << 
357   if (aTrack.GetTrackStatus() != fAlive &&     << 
358       aTrack.GetTrackStatus() != fSuspend) {   << 
359     if (aTrack.GetTrackStatus() == fStopAndKil << 
360         aTrack.GetTrackStatus() == fKillTrackA << 
361         aTrack.GetTrackStatus() == fPostponeTo << 
362       G4ExceptionDescription ed;               << 
363       ed << "G4HadronicProcess: track in unusa << 
364    << aTrack.GetTrackStatus() << G4endl;       << 
365       ed << "G4HadronicProcess: returning unch << 
366       DumpState(aTrack,"PostStepDoIt",ed);     << 
367       G4Exception("G4HadronicProcess::PostStep << 
368     }                                          << 
369     // No warning for fStopButAlive which is a << 
370     return theTotalResult;                     << 
371   }                                                91   }
372                                                    92 
373   // Initialize the hadronic projectile from t <<  93   G4VParticleChange * G4HadronicProcess::
374   thePro.Initialise(aTrack);                   <<  94   DoIsotopeCounting(G4VParticleChange * aResult,
375                                                <<  95                     const G4Track & aTrack,
376   theInteraction = ChooseHadronicInteraction(t <<  96                     const G4Nucleus & aNucleus)
377                                              a << 
378   if(nullptr == theInteraction) {              << 
379     G4ExceptionDescription ed;                 << 
380     ed << "Target element "<<anElement->GetNam << 
381        << targetNucleus.GetZ_asInt() << "  A=  << 
382        << targetNucleus.GetA_asInt() << G4endl << 
383     DumpState(aTrack,"ChooseHadronicInteractio << 
384     ed << " No HadronicInteraction found out"  << 
385     G4Exception("G4HadronicProcess::PostStepDo << 
386                 FatalException, ed);           << 
387     return theTotalResult;                     << 
388   }                                            << 
389                                                << 
390   G4HadFinalState* result = nullptr;           << 
391   G4int reentryCount = 0;                      << 
392   /*                                           << 
393   G4cout << "### " << aParticle->GetDefinition << 
394    << "  Ekin(MeV)= " << aParticle->GetKinetic << 
395    << "  Z= " << targetNucleus.GetZ_asInt()    << 
396    << "  A= " << targetNucleus.GetA_asInt()    << 
397    << "  by " << theInteraction->GetModelName( << 
398    << G4endl;                                  << 
399   */                                           << 
400   do                                           << 
401   {                                                97   {
402     try                                        <<  98     // get the PC from iso-production
                                                   >>  99     if(theOldIsoResult) delete theOldIsoResult;
                                                   >> 100     if(theIsoResult) delete theIsoResult;
                                                   >> 101     theIsoResult = new G4IsoParticleChange;
                                                   >> 102     G4bool done = false;
                                                   >> 103     G4IsoResult * anIsoResult = NULL;
                                                   >> 104     for(G4int i=0; i<theProductionModels.length(); i++)
403     {                                             105     {
404       // Call the interaction                  << 106       anIsoResult = theProductionModels(i)->GetIsotope(aTrack, aNucleus);
405       result = theInteraction->ApplyYourself(  << 107       if(anIsoResult!=NULL)
406       ++reentryCount;                          << 108       {
407     }                                          << 109         done = true;
408     catch(G4HadronicException & aR)            << 110         break;
                                                   >> 111       }
                                                   >> 112     }
                                                   >> 113     // if none in charge, use default iso production
                                                   >> 114     if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult); 
                                                   >> 115         
                                                   >> 116     // Add all info explicitely and add typename from model called.
                                                   >> 117     theIsoResult->SetIsotope(anIsoResult->GetIsotope());
                                                   >> 118     theIsoResult->SetProductionPosition(aTrack.GetPosition());
                                                   >> 119     theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
                                                   >> 120     theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
                                                   >> 121     theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
                                                   >> 122 //    theIsoResult->SetProducer(typeid(*theInteraction).name()); @@@@@@@
                                                   >> 123     G4String aWorkaround("WaitingForTypeidToBeAvailableInCompilers"); // @@@@@  workaround for DEC.
                                                   >> 124     theIsoResult->SetProducer(aWorkaround);
                                                   >> 125     
                                                   >> 126     delete anIsoResult;
                                                   >> 127 
                                                   >> 128     return aResult;
                                                   >> 129   }
                                                   >> 130   
                                                   >> 131   G4IsoResult * G4HadronicProcess::
                                                   >> 132   ExtractResidualNucleus(const G4Track & aTrack,
                                                   >> 133                          const G4Nucleus & aNucleus,
                                                   >> 134                          G4VParticleChange * aResult)
                                                   >> 135   {
                                                   >> 136     G4double A = aNucleus.GetN();
                                                   >> 137     G4double Z = aNucleus.GetZ();
                                                   >> 138     G4double bufferA = 0;
                                                   >> 139     G4double bufferZ = 0;
                                                   >> 140     
                                                   >> 141     // loop over aResult, and decrement A, Z accordingly
                                                   >> 142     // cash the max
                                                   >> 143     for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
409     {                                             144     {
410       G4ExceptionDescription ed;               << 145       G4Track* aSecTrack = aResult->GetSecondary(i);
411       aR.Report(ed);                           << 146       if(bufferA<aSecTrack->GetDefinition()->GetBaryonNumber())
412       ed << "Call for " << theInteraction->Get << 147       {
413       ed << "Target element "<<anElement->GetN << 148         bufferA = aSecTrack->GetDefinition()->GetBaryonNumber();
414    << targetNucleus.GetZ_asInt()               << 149         bufferZ = aSecTrack->GetDefinition()->GetPDGCharge();
415    << "  A= " << targetNucleus.GetA_asInt() << << 150       }
416       DumpState(aTrack,"ApplyYourself",ed);    << 151       Z-=aSecTrack->GetDefinition()->GetPDGCharge();
417       ed << " ApplyYourself failed" << G4endl; << 152       A-=aSecTrack->GetDefinition()->GetBaryonNumber();
418       G4Exception("G4HadronicProcess::PostStep << 153     }
419       ed);                                     << 154     
420     }                                          << 155     // if the fragment was part of the final state, it is 
421                                                << 156     // assumed to be the heaviest secondary.
422     // Check the result for catastrophic energ << 157     if(A<0.1)
423     result = CheckResult(thePro, targetNucleus << 158     {
424                                                << 159        A = bufferA;
425     if(reentryCount>100) {                     << 160        Z = bufferZ;
426       G4ExceptionDescription ed;               << 
427       ed << "Call for " << theInteraction->Get << 
428       ed << "Target element "<<anElement->GetN << 
429    << targetNucleus.GetZ_asInt()               << 
430    << "  A= " << targetNucleus.GetA_asInt() << << 
431       DumpState(aTrack,"ApplyYourself",ed);    << 
432       ed << " ApplyYourself does not completed << 
433       G4Exception("G4HadronicProcess::PostStep << 
434       ed);                                     << 
435     }                                          << 
436   }                                            << 
437   while(!result);  /* Loop checking, 30-Oct-20 << 
438                                                << 
439   // Check whether kaon0 or anti_kaon0 are pre << 
440   // if this is the case, transform them into  << 
441   // with equal, 50% probability, keeping thei << 
442   // the other kinematical properties).        << 
443   // When this happens - very rarely - a "Just << 
444   // Because Fluka-Cern produces kaon0 and ant << 
445   // of warnings to max 1 per thread.          << 
446   G4int nSec = (G4int)result->GetNumberOfSecon << 
447   if ( nSec > 0 ) {                            << 
448     for ( G4int i = 0; i < nSec; ++i ) {       << 
449       auto dynamicParticle = result->GetSecond << 
450       auto part = dynamicParticle->GetParticle << 
451       if ( part == G4KaonZero::Definition() || << 
452            part == G4AntiKaonZero::Definition( << 
453         G4ParticleDefinition* newPart;         << 
454         if ( G4UniformRand() > 0.5 ) { newPart << 
455         else { newPart = G4KaonZeroLong::Defin << 
456         dynamicParticle->SetDefinition( newPar << 
457   if ( nKaonWarn < 1 ) {                       << 
458     ++nKaonWarn;                               << 
459     G4ExceptionDescription ed;                 << 
460     ed << " Hadronic model " << theInteraction << 
461     ed << " created " << part->GetParticleName << 
462     ed << " -> forced to be " << newPart->GetP << 
463     G4Exception( "G4HadronicProcess::PostStepD << 
464   }                                            << 
465       }                                        << 
466     }                                          << 
467   }                                            << 
468                                                << 
469   result->SetTrafoToLab(thePro.GetTrafoToLab() << 
470   FillResult(result, aTrack);                  << 
471                                                << 
472   if (epReportLevel != 0) {                    << 
473     CheckEnergyMomentumConservation(aTrack, ta << 
474   }                                            << 
475   //G4cout << "PostStepDoIt done nICelectrons= << 
476   return theTotalResult;                       << 
477 }                                              << 
478                                                << 
479 void G4HadronicProcess::ProcessDescription(std << 
480 {                                              << 
481   outFile << "The description for this process << 
482 }                                              << 
483                                                << 
484 G4double G4HadronicProcess::XBiasSurvivalProba << 
485 {                                              << 
486   G4double nLTraversed = GetTotalNumberOfInter << 
487   G4double biasedProbability = 1.-G4Exp(-nLTra << 
488   G4double realProbability = 1-G4Exp(-nLTraver << 
489   G4double result = (biasedProbability-realPro << 
490   return result;                               << 
491 }                                              << 
492                                                << 
493 G4double G4HadronicProcess::XBiasSecondaryWeig << 
494 {                                              << 
495   G4double nLTraversed = GetTotalNumberOfInter << 
496   G4double result =                            << 
497      1./aScaleFactor*G4Exp(-nLTraversed/aScale << 
498   return result;                               << 
499 }                                              << 
500                                                << 
501 void                                           << 
502 G4HadronicProcess::FillResult(G4HadFinalState  << 
503 {                                              << 
504   theTotalResult->ProposeLocalEnergyDeposit(aR << 
505   const G4ThreeVector& dir = aT.GetMomentumDir << 
506                                                << 
507   G4double efinal = std::max(aR->GetEnergyChan << 
508                                                << 
509   // check status of primary                   << 
510   if(aR->GetStatusChange() == stopAndKill) {   << 
511     theTotalResult->ProposeTrackStatus(fStopAn << 
512     theTotalResult->ProposeEnergy( 0.0 );      << 
513                                                << 
514     // check its final energy                  << 
515   } else if(0.0 == efinal) {                   << 
516     theTotalResult->ProposeEnergy( 0.0 );      << 
517     if(aT.GetParticleDefinition()->GetProcessM << 
518        ->GetAtRestProcessVector()->size() > 0) << 
519          { theTotalResult->ProposeTrackStatus( << 
520     else { theTotalResult->ProposeTrackStatus( << 
521                                                << 
522     // primary is not killed apply rotation an << 
523   } else  {                                    << 
524     theTotalResult->ProposeTrackStatus(fAlive) << 
525     G4ThreeVector newDir = aR->GetMomentumChan << 
526     newDir.rotateUz(dir);                      << 
527     theTotalResult->ProposeMomentumDirection(n << 
528     theTotalResult->ProposeEnergy(efinal);     << 
529   }                                            << 
530   //G4cout << "FillResult: Efinal= " << efinal << 
531   //   << theTotalResult->GetTrackStatus()     << 
532   //   << "  fKill= " << fStopAndKill << G4end << 
533                                                << 
534   // check secondaries                         << 
535   nICelectrons = 0;                            << 
536   G4int nSec = (G4int)aR->GetNumberOfSecondari << 
537   theTotalResult->SetNumberOfSecondaries(nSec) << 
538   G4double time0 = aT.GetGlobalTime();         << 
539                                                << 
540   for (G4int i = 0; i < nSec; ++i) {           << 
541     G4DynamicParticle* dynParticle = aR->GetSe << 
542                                                << 
543     // apply rotation                          << 
544     G4ThreeVector newDir = dynParticle->GetMom << 
545     newDir.rotateUz(dir);                      << 
546     dynParticle->SetMomentumDirection(newDir); << 
547                                                << 
548     // check if secondary is on the mass shell << 
549     const G4ParticleDefinition* part = dynPart << 
550     G4double mass = part->GetPDGMass();        << 
551     G4double dmass= dynParticle->GetMass();    << 
552     const G4double delta_mass_lim = 1.0*CLHEP: << 
553     const G4double delta_ekin = 0.001*CLHEP::e << 
554     if(std::abs(dmass - mass) > delta_mass_lim << 
555       G4double e =                             << 
556         std::max(dynParticle->GetKineticEnergy << 
557       if(verboseLevel > 1) {                   << 
558   G4ExceptionDescription ed;                   << 
559   ed << "TrackID= "<< aT.GetTrackID()          << 
560      << "  " << aT.GetParticleDefinition()->Ge << 
561      << " Target Z= " << targetNucleus.GetZ_as << 
562      << targetNucleus.GetA_asInt()             << 
563      << " Ekin(GeV)= " << aT.GetKineticEnergy( << 
564      << "\n Secondary is out of mass shell: "  << 
565      << "  EkinNew(MeV)= " << e                << 
566      << " DeltaMass(MeV)= " << dmass - mass << << 
567   G4Exception("G4HadronicProcess::FillResults" << 
568       }                                        << 
569       dynParticle->SetKineticEnergy(e);        << 
570       dynParticle->SetMass(mass);              << 
571     }                                          << 
572     G4int idModel = aR->GetSecondary(i)->GetCr << 
573     if(part->GetPDGEncoding() == 11) { ++nICel << 
574                                                << 
575     // time of interaction starts from zero +  << 
576     G4double time = std::max(aR->GetSecondary( << 
577                                                << 
578     G4Track* track = new G4Track(dynParticle,  << 
579     track->SetCreatorModelID(idModel);         << 
580     track->SetParentResonanceDef(aR->GetSecond << 
581     track->SetParentResonanceID(aR->GetSeconda << 
582     G4double newWeight = fWeight*aR->GetSecond << 
583     track->SetWeight(newWeight);               << 
584     track->SetTouchableHandle(aT.GetTouchableH << 
585     theTotalResult->AddSecondary(track);       << 
586   }                                            << 
587   aR->Clear();                                 << 
588   // G4cout << "FillResults done nICe= " << nI << 
589 }                                              << 
590                                                << 
591 void G4HadronicProcess::MultiplyCrossSectionBy << 
592 {                                              << 
593   BiasCrossSectionByFactor(factor);            << 
594 }                                              << 
595                                                << 
596 void G4HadronicProcess::BiasCrossSectionByFact << 
597 {                                              << 
598   if (aScale <= 0.0) {                         << 
599     G4ExceptionDescription ed;                 << 
600     ed << " Wrong biasing factor " << aScale < << 
601     G4Exception("G4HadronicProcess::BiasCrossS << 
602                 JustWarning, ed, "Cross-sectio << 
603   } else {                                     << 
604     aScaleFactor = aScale;                     << 
605   }                                            << 
606 }                                              << 
607                                                << 
608 G4HadFinalState* G4HadronicProcess::CheckResul << 
609             const G4Nucleus &aNucleus,         << 
610             G4HadFinalState * result)          << 
611 {                                              << 
612   // check for catastrophic energy non-conserv << 
613   // to re-sample the interaction              << 
614   G4HadronicInteraction * theModel = GetHadron << 
615   G4double nuclearMass(0);                     << 
616   if (nullptr != theModel) {                   << 
617                                                << 
618     // Compute final-state total energy        << 
619     G4double finalE(0.);                       << 
620     G4int nSec = (G4int)result->GetNumberOfSec << 
621                                                << 
622     nuclearMass = G4NucleiProperties::GetNucle << 
623                  aNucleus.GetZ_asInt());       << 
624     if (result->GetStatusChange() != stopAndKi << 
625       // Interaction didn't complete, returned << 
626       // and reset nucleus or the primary surv << 
627       // (e.g. electro-nuclear ) => keep  nucl << 
628       finalE=result->GetLocalEnergyDeposit() + << 
629              aPro.GetDefinition()->GetPDGMass( << 
630       if( nSec == 0 ){                         << 
631          // Since there are no secondaries, th << 
632          // To check energy balance we must ne << 
633         nuclearMass=0.0;                       << 
634       }                                        << 
635     }                                          << 
636     for (G4int i = 0; i < nSec; ++i) {         << 
637       G4DynamicParticle *pdyn=result->GetSecon << 
638       finalE += pdyn->GetTotalEnergy();        << 
639       G4double mass_pdg=pdyn->GetDefinition()- << 
640       G4double mass_dyn=pdyn->GetMass();       << 
641       if ( std::abs(mass_pdg - mass_dyn) > 0.1 << 
642         // If it is shortlived, then a differe << 
643         if ( pdyn->GetDefinition()->IsShortLiv << 
644              std::abs(mass_pdg - mass_dyn) < 3 << 
645           continue;                            << 
646         }                                      << 
647   result->Clear();                             << 
648   result = nullptr;                            << 
649   G4ExceptionDescription desc;                 << 
650   desc << "Warning: Secondary with off-shell d << 
651        << G4endl                               << 
652        << " " << pdyn->GetDefinition()->GetPar << 
653        << ", PDG mass: " << mass_pdg << ", dyn << 
654        << mass_dyn << G4endl                   << 
655        << (epReportLevel<0 ? "abort the event" << 
656      : "re-sample the interaction") << G4endl  << 
657        << " Process / Model: " <<  GetProcessN << 
658        << theModel->GetModelName() << G4endl   << 
659        << " Primary: " << aPro.GetDefinition() << 
660        << " (" << aPro.GetDefinition()->GetPDG << 
661        << " E= " <<  aPro.Get4Momentum().e()   << 
662        << ", target nucleus (" << aNucleus.Get << 
663        << aNucleus.GetA_asInt() << ")" << G4en << 
664   G4Exception("G4HadronicProcess:CheckResult() << 
665         epReportLevel<0 ? EventMustBeAborted : << 
666   // must return here.....                     << 
667   return result;                               << 
668       }                                        << 
669     }                                          << 
670     G4double deltaE= nuclearMass +  aPro.GetTo << 
671                                                << 
672     std::pair<G4double, G4double> checkLevels  << 
673       theModel->GetFatalEnergyCheckLevels();   << 
674     if (std::abs(deltaE) > checkLevels.second  << 
675         std::abs(deltaE) > checkLevels.first*a << 
676       // do not delete result, this is a point << 
677       result->Clear();                         << 
678       result = nullptr;                        << 
679       G4ExceptionDescription desc;             << 
680       desc << "Warning: Bad energy non-conserv << 
681      << (epReportLevel<0 ? "abort the event"   << 
682          :  "re-sample the interaction") << G4 << 
683      << " Process / Model: " <<  GetProcessNam << 
684      << theModel->GetModelName() << G4endl     << 
685      << " Primary: " << aPro.GetDefinition()-> << 
686      << " (" << aPro.GetDefinition()->GetPDGEn << 
687      << " E= " <<  aPro.Get4Momentum().e()     << 
688      << ", target nucleus (" << aNucleus.GetZ_ << 
689      << aNucleus.GetA_asInt() << ")" << G4endl << 
690      << " E(initial - final) = " << deltaE <<  << 
691       G4Exception("G4HadronicProcess:CheckResu << 
692       epReportLevel<0 ? EventMustBeAborted : J << 
693     }                                          << 
694   }                                            << 
695   return result;                               << 
696 }                                              << 
697                                                << 
698 void                                           << 
699 G4HadronicProcess::CheckEnergyMomentumConserva << 
700                                                << 
701 {                                              << 
702   G4int target_A=aNucleus.GetA_asInt();        << 
703   G4int target_Z=aNucleus.GetZ_asInt();        << 
704   G4double targetMass = G4NucleiProperties::Ge << 
705   G4LorentzVector target4mom(0, 0, 0, targetMa << 
706            + nICelectrons*CLHEP::electron_mass << 
707                                                << 
708   G4LorentzVector projectile4mom = aTrack.GetD << 
709   G4int track_A = aTrack.GetDefinition()->GetB << 
710   G4int track_Z = G4lrint(aTrack.GetDefinition << 
711                                                << 
712   G4int initial_A = target_A + track_A;        << 
713   G4int initial_Z = target_Z + track_Z - nICel << 
714                                                << 
715   G4LorentzVector initial4mom = projectile4mom << 
716                                                << 
717   // Compute final-state momentum for scatteri << 
718   G4LorentzVector final4mom;                   << 
719   G4int final_A(0), final_Z(0);                << 
720                                                << 
721   G4int nSec = theTotalResult->GetNumberOfSeco << 
722   if (theTotalResult->GetTrackStatus() != fSto << 
723     // Either interaction didn't complete, ret << 
724     //  or    the primary survived the interac << 
725                                                << 
726     // Interaction didn't complete, returned " << 
727     //   - or suppressed recoil  (e.g. Neutron << 
728     final4mom = initial4mom;                   << 
729     final_A = initial_A;                       << 
730     final_Z = initial_Z;                       << 
731     if (nSec > 0) {                            << 
732       // The primary remains in final state (e << 
733       // Use the final energy / momentum       << 
734       const G4ThreeVector& v = *theTotalResult << 
735       G4double ekin = theTotalResult->GetEnerg << 
736       G4double mass = aTrack.GetDefinition()-> << 
737       G4double ptot = std::sqrt(ekin*(ekin + 2 << 
738       final4mom.set(ptot*v.x(), ptot*v.y(), pt << 
739       final_A = track_A;                       << 
740       final_Z = track_Z;                       << 
741       // Expect that the target nucleus will h << 
742       //  and its products, including recoil,  << 
743     }                                          << 
744   }                                            << 
745   if( nSec > 0 ) {                             << 
746     G4Track* sec;                              << 
747                                                << 
748     for (G4int i = 0; i < nSec; i++) {         << 
749       sec = theTotalResult->GetSecondary(i);   << 
750       final4mom += sec->GetDynamicParticle()-> << 
751       final_A += sec->GetDefinition()->GetBary << 
752       final_Z += G4lrint(sec->GetDefinition()- << 
753     }                                             161     }
                                                   >> 162     
                                                   >> 163     // prepare the IsoResult.
                                                   >> 164     char the1[100] = {""};
                                                   >> 165     ostrstream ost1(the1, 100, ios::out);
                                                   >> 166     ost1 <<Z<<"_"<<A;
                                                   >> 167     G4String * biff = new G4String(the1);
                                                   >> 168     G4IsoResult * theResult = new G4IsoResult(*biff, aNucleus);
                                                   >> 169     
                                                   >> 170     // cleaning up.
                                                   >> 171     delete biff;
                                                   >> 172     
                                                   >> 173     return theResult;
754   }                                               174   }
755                                                   175 
756   // Get level-checking information (used to c << 
757   G4String processName = GetProcessName();     << 
758   G4HadronicInteraction* theModel = GetHadroni << 
759   G4String modelName("none");                  << 
760   if (theModel) modelName = theModel->GetModel << 
761   std::pair<G4double, G4double> checkLevels =  << 
762   if (!levelsSetByProcess) {                   << 
763     if (theModel) checkLevels = theModel->GetE << 
764     checkLevels.first= std::min(checkLevels.fi << 
765     checkLevels.second=std::min(checkLevels.se << 
766   }                                            << 
767                                                << 
768   // Compute absolute total-energy difference, << 
769   G4bool checkRelative = (aTrack.GetKineticEne << 
770                                                << 
771   G4LorentzVector diff = initial4mom - final4m << 
772   G4double absolute = diff.e();                << 
773   G4double relative = checkRelative ? absolute << 
774                                                << 
775   G4double absolute_mom = diff.vect().mag();   << 
776   G4double relative_mom = checkRelative ? abso << 
777                                                << 
778   // Evaluate relative and absolute conservati << 
779   G4bool relPass = true;                       << 
780   G4String relResult = "pass";                 << 
781   if (  std::abs(relative) > checkLevels.first << 
782    || std::abs(relative_mom) > checkLevels.fir << 
783     relPass = false;                           << 
784     relResult = checkRelative ? "fail" : "N/A" << 
785   }                                            << 
786                                                << 
787   G4bool absPass = true;                       << 
788   G4String absResult = "pass";                 << 
789   if (   std::abs(absolute) > checkLevels.seco << 
790       || std::abs(absolute_mom) > checkLevels. << 
791     absPass = false ;                          << 
792     absResult = "fail";                        << 
793   }                                            << 
794                                                << 
795   G4bool chargePass = true;                    << 
796   G4String chargeResult = "pass";              << 
797   if (   (initial_A-final_A)!=0                << 
798       || (initial_Z-final_Z)!=0 ) {            << 
799     chargePass = checkLevels.second < DBL_MAX  << 
800     chargeResult = "fail";                     << 
801    }                                           << 
802                                                << 
803   G4bool conservationPass = (relPass || absPas << 
804                                                << 
805   std::stringstream Myout;                     << 
806   G4bool Myout_notempty(false);                << 
807   // Options for level of reporting detail:    << 
808   //  0. off                                   << 
809   //  1. report only when E/p not conserved    << 
810   //  2. report regardless of E/p conservation << 
811   //  3. report only when E/p not conserved, w << 
812   //  4. report regardless of E/p conservation << 
813   //  negative -1.., as above, but send output << 
814                                                << 
815   if(   std::abs(epReportLevel) == 4           << 
816   ||  ( std::abs(epReportLevel) == 3 && ! cons << 
817       Myout << " Process: " << processName <<  << 
818       Myout << " Primary: " << aTrack.GetParti << 
819             << " (" << aTrack.GetParticleDefin << 
820             << " E= " <<  aTrack.GetDynamicPar << 
821       << ", target nucleus (" << aNucleus.GetZ << 
822       << aNucleus.GetA_asInt() << ")" << G4end << 
823       Myout_notempty=true;                     << 
824   }                                            << 
825   if (  std::abs(epReportLevel) == 4           << 
826    || std::abs(epReportLevel) == 2             << 
827    || ! conservationPass ){                    << 
828                                                << 
829       Myout << "   "<< relResult  <<" relative << 
830              << relative << " p/p(0)= " << rel << 
831       Myout << "   "<< absResult << " absolute << 
832              << absolute/MeV << " / " << absol << 
833       Myout << "   "<< chargeResult << " charg << 
834       Myout_notempty=true;                     << 
835                                                << 
836   }                                            << 
837   Myout.flush();                               << 
838   if ( Myout_notempty ) {                      << 
839      if (epReportLevel > 0)      G4cout << Myo << 
840      else if (epReportLevel < 0) G4cerr << Myo << 
841   }                                            << 
842 }                                              << 
843                                                << 
844 void G4HadronicProcess::DumpState(const G4Trac << 
845           const G4String& method,              << 
846           G4ExceptionDescription& ed)          << 
847 {                                              << 
848   ed << "Unrecoverable error in the method " < << 
849      << GetProcessName() << G4endl;            << 
850   ed << "TrackID= "<< aTrack.GetTrackID() << " << 
851      << aTrack.GetParentID()                   << 
852      << "  " << aTrack.GetParticleDefinition() << 
853      << G4endl;                                << 
854   ed << "Ekin(GeV)= " << aTrack.GetKineticEner << 
855      << ";  direction= " << aTrack.GetMomentum << 
856   ed << "Position(mm)= " << aTrack.GetPosition << 
857                                                << 
858   if (aTrack.GetMaterial()) {                  << 
859     ed << "  material " << aTrack.GetMaterial( << 
860   }                                            << 
861   ed << G4endl;                                << 
862                                                << 
863   if (aTrack.GetVolume()) {                    << 
864     ed << "PhysicalVolume  <" << aTrack.GetVol << 
865        << ">" << G4endl;                       << 
866   }                                            << 
867 }                                              << 
868                                                << 
869 void G4HadronicProcess::DumpPhysicsTable(const << 
870 {                                              << 
871   theCrossSectionDataStore->DumpPhysicsTable(p << 
872 }                                              << 
873                                                << 
874 void G4HadronicProcess::AddDataSet(G4VCrossSec << 
875 {                                              << 
876   theCrossSectionDataStore->AddDataSet(aDataSe << 
877 }                                              << 
878                                                << 
879 std::vector<G4HadronicInteraction*>&           << 
880 G4HadronicProcess::GetHadronicInteractionList( << 
881 {                                              << 
882   return theEnergyRangeManager.GetHadronicInte << 
883 }                                              << 
884                                                << 
885 G4HadronicInteraction*                         << 
886 G4HadronicProcess::GetHadronicModel(const G4St << 
887 {                                              << 
888   std::vector<G4HadronicInteraction*>& list    << 
889         = theEnergyRangeManager.GetHadronicInt << 
890   for (auto & mod : list) {                    << 
891     if (mod->GetModelName() == modelName) retu << 
892   }                                            << 
893   return nullptr;                              << 
894 }                                              << 
895                                                << 
896 G4double                                       << 
897 G4HadronicProcess::ComputeCrossSection(const G << 
898                const G4Material* mat,          << 
899                const G4double kinEnergy)       << 
900 {                                              << 
901   auto dp = new G4DynamicParticle(part, unitVe << 
902   G4double xs = theCrossSectionDataStore->Comp << 
903   delete dp;                                   << 
904   return xs;                                   << 
905 }                                              << 
906                                                << 
907 void G4HadronicProcess::RecomputeXSandMFP(cons << 
908 {                                              << 
909   auto dp = new G4DynamicParticle(currentParti << 
910   theLastCrossSection = aScaleFactor*          << 
911     theCrossSectionDataStore->ComputeCrossSect << 
912   theMFP = (theLastCrossSection > 0.0) ? 1.0/t << 
913   delete dp;                                   << 
914 }                                              << 
915                                                << 
916 void G4HadronicProcess::UpdateCrossSectionAndM << 
917 {                                              << 
918   if(fXSType == fHadNoIntegral) {              << 
919     DefineXSandMFP();                          << 
920                                                << 
921   } else if(fXSType == fHadIncreasing) {       << 
922     if(e*invLambdaFactor < mfpKinEnergy) {     << 
923       mfpKinEnergy = e;                        << 
924       ComputeXSandMFP();                       << 
925     }                                          << 
926                                                << 
927   } else if(fXSType == fHadDecreasing) {       << 
928     if(e < mfpKinEnergy && mfpKinEnergy > minK << 
929       G4double e1 = std::max(e*lambdaFactor, m << 
930       mfpKinEnergy = e1;                       << 
931       RecomputeXSandMFP(e1);                   << 
932     }                                          << 
933                                                << 
934   } else if(fXSType == fHadOnePeak) {          << 
935     G4double epeak = (*theEnergyOfCrossSection << 
936     if(e <= epeak) {                           << 
937       if(e*invLambdaFactor < mfpKinEnergy) {   << 
938         mfpKinEnergy = e;                      << 
939   ComputeXSandMFP();                           << 
940       }                                        << 
941     } else if(e < mfpKinEnergy) {              << 
942       G4double e1 = std::max(epeak, e*lambdaFa << 
943       mfpKinEnergy = e1;                       << 
944       RecomputeXSandMFP(e1);                   << 
945     }                                          << 
946                                                << 
947   } else if(fXSType == fHadTwoPeaks) {         << 
948     G4TwoPeaksHadXS* xs = (*fXSpeaks)[matIdx]; << 
949     const G4double e1peak = xs->e1peak;        << 
950                                                << 
951     // below the 1st peak                      << 
952     if(e <= e1peak) {                          << 
953       if(e*invLambdaFactor < mfpKinEnergy) {   << 
954         mfpKinEnergy = e;                      << 
955   ComputeXSandMFP();                           << 
956       }                                        << 
957       return;                                  << 
958     }                                          << 
959     const G4double e1deep = xs->e1deep;        << 
960     // above the 1st peak, below the deep      << 
961     if(e <= e1deep) {                          << 
962       if(mfpKinEnergy >= e1deep || e <= mfpKin << 
963         const G4double e1 = std::max(e1peak, e << 
964         mfpKinEnergy = e1;                     << 
965   RecomputeXSandMFP(e1);                       << 
966       }                                        << 
967       return;                                  << 
968     }                                          << 
969     const G4double e2peak = xs->e2peak;        << 
970     // above the deep, below 2nd peak          << 
971     if(e <= e2peak) {                          << 
972       if(e*invLambdaFactor < mfpKinEnergy) {   << 
973         mfpKinEnergy = e;                      << 
974   ComputeXSandMFP();                           << 
975       }                                        << 
976       return;                                  << 
977     }                                          << 
978     const G4double e2deep = xs->e2deep;        << 
979     // above the 2nd peak, below the deep      << 
980     if(e <= e2deep) {                          << 
981       if(mfpKinEnergy >= e2deep || e <= mfpKin << 
982         const G4double e1 = std::max(e2peak, e << 
983         mfpKinEnergy = e1;                     << 
984   RecomputeXSandMFP(e1);                       << 
985       }                                        << 
986       return;                                  << 
987     }                                          << 
988     const G4double e3peak = xs->e3peak;        << 
989     // above the deep, below 3d peak           << 
990     if(e <= e3peak) {                          << 
991       if(e*invLambdaFactor < mfpKinEnergy) {   << 
992         mfpKinEnergy = e;                      << 
993   ComputeXSandMFP();                           << 
994       }                                        << 
995       return;                                  << 
996     }                                          << 
997     // above 3d peak                           << 
998     if(e <= mfpKinEnergy) {                    << 
999       const G4double e1 = std::max(e3peak, e*l << 
1000       mfpKinEnergy = e1;                      << 
1001       RecomputeXSandMFP(e1);                  << 
1002     }                                         << 
1003                                                  176 
1004   } else {                                    << 177  /* end of file */
1005     DefineXSandMFP();                         << 
1006   }                                           << 
1007 }                                             << 
1008                                                  178