| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 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" <<
51 27
52 #include "G4Types.hh" 28 #include "G4Types.hh"
53 #include "G4SystemOfUnits.hh" << 29
>> 30 #include <fstream>
>> 31 #include <sstream>
>> 32 #include <stdlib.h>
>> 33 #include "G4HadronicProcess.hh"
>> 34 // #include "G4EffectiveCharge.hh"
54 #include "G4HadProjectile.hh" 35 #include "G4HadProjectile.hh"
55 #include "G4ElementVector.hh" 36 #include "G4ElementVector.hh"
56 #include "G4Track.hh" 37 #include "G4Track.hh"
57 #include "G4Step.hh" 38 #include "G4Step.hh"
58 #include "G4Element.hh" 39 #include "G4Element.hh"
59 #include "G4ParticleChange.hh" 40 #include "G4ParticleChange.hh"
>> 41 #include "G4TransportationManager.hh"
>> 42 #include "G4Navigator.hh"
60 #include "G4ProcessVector.hh" 43 #include "G4ProcessVector.hh"
61 #include "G4ProcessManager.hh" 44 #include "G4ProcessManager.hh"
62 #include "G4NucleiProperties.hh" << 45 #include "G4StableIsotopes.hh"
>> 46 #include "G4HadTmpUtil.hh"
63 47
>> 48 #include "G4HadLeadBias.hh"
64 #include "G4HadronicException.hh" 49 #include "G4HadronicException.hh"
65 #include "G4HadronicProcessStore.hh" << 50 #include "G4HadReentrentException.hh"
66 #include "G4HadronicParameters.hh" << 51 #include "G4HadronicInteractionWrapper.hh"
67 #include "G4VCrossSectionDataSet.hh" << 52
68 << 53 #include "G4HadSignalHandler.hh"
69 #include "G4NistManager.hh" <<
70 #include "G4VLeadingParticleBiasing.hh" <<
71 #include "G4HadXSHelper.hh" <<
72 #include "G4Threading.hh" <<
73 #include "G4Exp.hh" <<
74 54
75 #include <typeinfo> 55 #include <typeinfo>
76 #include <sstream> <<
77 #include <iostream> <<
78 56
79 namespace << 57 namespace G4HadronicProcess_local
80 { 58 {
81 constexpr G4double lambdaFactor = 0.8; << 59 extern "C" void G4HadronicProcessHandler_1(int)
82 constexpr G4double invLambdaFactor = 1.0/lam << 60 {
83 } << 61 G4HadronicWhiteBoard::Instance().Dump();
>> 62 }
>> 63 }
84 64
85 ////////////////////////////////////////////// << 65 G4IsoParticleChange * G4HadronicProcess::theIsoResult = 0;
>> 66 G4IsoParticleChange * G4HadronicProcess::theOldIsoResult = 0;
>> 67 G4bool G4HadronicProcess::isoIsEnabled = true;
86 68
87 G4HadronicProcess::G4HadronicProcess(const G4S << 69 void G4HadronicProcess::
88 G4Process << 70 EnableIsotopeProductionGlobally() {isoIsEnabled = true;}
89 : G4VDiscreteProcess(processName, procType) <<
90 { <<
91 SetProcessSubType(fHadronInelastic); // Def <<
92 InitialiseLocal(); <<
93 } <<
94 71
95 G4HadronicProcess::G4HadronicProcess(const G4S << 72 void G4HadronicProcess::
96 G4Hadroni << 73 DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
97 : G4VDiscreteProcess(processName, fHadronic) << 74
98 { << 75 G4HadronicProcess::G4HadronicProcess( const G4String &processName,
99 SetProcessSubType(aHadSubType); << 76 G4ProcessType aType ) :
100 InitialiseLocal(); << 77 G4VDiscreteProcess( processName, aType)
>> 78 {
>> 79 ModelingState = 0;
>> 80 isoIsOnAnyway = -1;
>> 81 theTotalResult = new G4ParticleChange();
>> 82 theCrossSectionDataStore = new G4CrossSectionDataStore();
>> 83 aScaleFactor = 1;
>> 84 xBiasOn = false;
>> 85 if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
101 } 86 }
102 87
103 G4HadronicProcess::~G4HadronicProcess() 88 G4HadronicProcess::~G4HadronicProcess()
104 { << 89 {
105 theProcessStore->DeRegister(this); <<
106 delete theTotalResult; 90 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 91
132 unitVector.set(0.0, 0.0, 0.1); << 92 std::for_each(theProductionModels.begin(),
133 if(G4Threading::IsWorkerThread()) { isMaster << 93 theProductionModels.end(), G4Delete());
>> 94 std::for_each(theBias.begin(), theBias.end(), G4Delete());
>> 95
>> 96 delete theOldIsoResult; delete theIsoResult;
>> 97 delete theCrossSectionDataStore;
134 } 98 }
135 99
136 void G4HadronicProcess::RegisterMe( G4Hadronic 100 void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
137 { << 101 {
138 if(nullptr == a) { return; } << 102 try{GetManagerPointer()->RegisterMe( a );}
139 theEnergyRangeManager.RegisterMe( a ); << 103 catch(G4HadronicException & aE)
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 { 104 {
150 static const G4int nmax = 5; << 105 aE.Report(std::cout);
151 if(nMatWarn < nmax) { << 106 G4Exception("G4HadronicProcess", "007", FatalException,
152 ++nMatWarn; << 107 "Could not register G4HadronicInteraction");
153 G4ExceptionDescription ed; <<
154 ed << "Cannot compute Element x-section <<
155 << " because no material defined \n" <<
156 << " Please, specify material pointer or de <<
157 << " for Z= " << elm->GetZasInt(); <<
158 G4Exception("G4HadronicProcess::GetEleme <<
159 JustWarning, ed); <<
160 } <<
161 } 108 }
162 return theCrossSectionDataStore->GetCrossSec <<
163 } 109 }
164 110
165 void G4HadronicProcess::PreparePhysicsTable(co << 111 G4double G4HadronicProcess::
166 { << 112 GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
167 if(nullptr == firstParticle) { firstParticle << 113 {
168 theProcessStore->RegisterParticle(this, &p); << 114 G4double sigma = 0.0;
>> 115 try
>> 116 {
>> 117 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
>> 118 if( !IsApplicable(*aParticle->GetDefinition()))
>> 119 {
>> 120 G4cout << "Unrecoverable error: "<<G4endl;
>> 121 G4ProcessManager * it = aParticle->GetDefinition()->GetProcessManager();
>> 122 G4ProcessVector * itv = it->GetProcessList();
>> 123 G4cout <<aParticle->GetDefinition()->GetParticleName()<<
>> 124 " has the following processes:"<<G4endl;
>> 125 for(G4int i=0; i<itv->size(); i++)
>> 126 {
>> 127 G4cout <<" "<<(*itv)[i]->GetProcessName()<<G4endl;
>> 128 }
>> 129 G4cout << "for kinetic energy "<<aParticle->GetKineticEnergy()<<G4endl;
>> 130 G4cout << "and material "<<aTrack.GetMaterial()->GetName()<<G4endl;
>> 131 G4Exception("G4HadronicProcess", "007", FatalException,
>> 132 std::string(this->GetProcessName()+
>> 133 " was called for "+
>> 134 aParticle->GetDefinition()->GetParticleName()).c_str() );
>> 135 }
>> 136 G4Material *aMaterial = aTrack.GetMaterial();
>> 137 ModelingState = 1;
>> 138
>> 139 sigma = theCrossSectionDataStore->GetCrossSection(aParticle, aMaterial);
>> 140
>> 141 sigma *= aScaleFactor;
>> 142 theLastCrossSection = sigma;
>> 143 }
>> 144 catch(G4HadronicException aR)
>> 145 {
>> 146 aR.Report(G4cout);
>> 147 G4Exception("G4HadronicProcess", "007", FatalException,
>> 148 "G4HadronicProcess::GetMeanFreePath failed");
>> 149 }
>> 150 if( sigma > 0.0 )
>> 151 return 1.0/sigma;
>> 152 else
>> 153 return DBL_MAX;
169 } 154 }
170 155
171 void G4HadronicProcess::BuildPhysicsTable(cons << 156
>> 157 G4Element* G4HadronicProcess::ChooseAandZ(
>> 158 const G4DynamicParticle *aParticle, const G4Material *aMaterial )
172 { 159 {
173 if(firstParticle != &p) { return; } << 160 std::pair<G4double, G4double> ZA =
>> 161 theCrossSectionDataStore->SelectRandomIsotope(aParticle, aMaterial);
>> 162 G4double ZZ = ZA.first;
>> 163 G4double AA = ZA.second;
174 164
175 theCrossSectionDataStore->BuildPhysicsTable( << 165 targetNucleus.SetParameters(AA, ZZ);
176 theEnergyRangeManager.BuildPhysicsTable(p); <<
177 G4HadronicParameters* param = G4HadronicPara <<
178 <<
179 G4int subtype = GetProcessSubType(); <<
180 if(useIntegralXS) { <<
181 if(subtype == fHadronInelastic) { <<
182 useIntegralXS = param->EnableIntegralIne <<
183 } else if(subtype == fHadronElastic) { <<
184 useIntegralXS = param->EnableIntegralEla <<
185 } <<
186 } <<
187 fXSType = fHadNoIntegral; <<
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 } <<
204 166
205 // check particle for integral method << 167 const G4int numberOfElements = aMaterial->GetNumberOfElements();
206 if(isMaster || nullptr == masterProcess) { << 168 const G4ElementVector* theElementVector = aMaterial->GetElementVector();
207 G4double charge = p.GetPDGCharge()/eplus; << 169 G4Element* chosen = 0;
208 << 170 for (G4int i = 0; i < numberOfElements; i++) {
209 // select cross section shape << 171 chosen = (*theElementVector)[i];
210 if(charge != 0.0 && useIntegralXS) { << 172 if (chosen->GetZ() == ZZ) break;
211 G4double tmax = param->GetMaxEnergy(); <<
212 currentParticle = firstParticle; <<
213 // initialisation in the master thread <<
214 G4int pdg = p.GetPDGEncoding(); <<
215 if (std::abs(pdg) == 211) { <<
216 fXSType = fHadTwoPeaks; <<
217 } else if (pdg == 321) { <<
218 fXSType = fHadOnePeak; <<
219 } else if (pdg == -321) { <<
220 fXSType = fHadDecreasing; <<
221 } else if (pdg == 2212) { <<
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 } <<
252 } <<
253 } else { <<
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 } 173 }
264 G4HadronicProcessStore::Instance()->PrintInf << 174 return chosen;
265 } 175 }
266 176
267 void G4HadronicProcess::StartTracking(G4Track* <<
268 { <<
269 currentMat = nullptr; <<
270 currentParticle = track->GetDefinition(); <<
271 fDynParticle = track->GetDynamicParticle(); <<
272 theNumberOfInteractionLengthLeft = -1.0; <<
273 } <<
274 177
275 G4double G4HadronicProcess::PostStepGetPhysica << 178 struct G4Nancheck{ bool operator()(G4double aV){return (!(aV<1))&&(!(aV>-1));}};
276 const G4Track& tr << 179
277 G4double previousStepSize, << 180 G4VParticleChange *G4HadronicProcess::GeneralPostStepDoIt(
278 G4ForceCondition* << 181 const G4Track &aTrack, const G4Step &)
279 { 182 {
280 *condition = NotForced; << 183 // Debugging stuff
281 184
282 const G4Material* mat = track.GetMaterial(); << 185 bool G4HadronicProcess_debug_flag = false;
283 if(mat != currentMat) { << 186 if(getenv("G4HadronicProcess_debug")) G4HadronicProcess_debug_flag = true;
284 currentMat = mat; << 187 if(G4HadronicProcess_debug_flag)
285 mfpKinEnergy = DBL_MAX; << 188 std::cout << "@@@@ hadronic process start "<< std::endl;
286 matIdx = (G4int)track.GetMaterial()->GetIn << 189 // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
287 } << 190 #ifndef G4HadSignalHandler_off
288 UpdateCrossSectionAndMFP(track.GetKineticEne << 191 G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
289 << 192 #endif
290 // zero cross section <<
291 if(theLastCrossSection <= 0.0) { <<
292 theNumberOfInteractionLengthLeft = -1.0; <<
293 currentInteractionLength = DBL_MAX; <<
294 return DBL_MAX; <<
295 } <<
296 193
297 // non-zero cross section << 194 if(aTrack.GetTrackStatus() != fAlive && aTrack.GetTrackStatus() != fSuspend)
298 if (theNumberOfInteractionLengthLeft < 0.0) << 195 {
299 theNumberOfInteractionLengthLeft = -G4Log( << 196 G4cerr << "G4HadronicProcess: track in unusable state - "
300 theInitialNumberOfInteractionLength = theN << 197 <<aTrack.GetTrackStatus()<<G4endl;
301 } else { << 198 G4cerr << "G4HadronicProcess: returning unchanged track "<<G4endl;
302 theNumberOfInteractionLengthLeft -= << 199 G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
303 previousStepSize/currentInteractionLengt << 200 theTotalResult->Clear();
304 theNumberOfInteractionLengthLeft = << 201 theTotalResult->Initialize(aTrack);
305 std::max(theNumberOfInteractionLengthLef << 202 return theTotalResult;
306 } << 203 }
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 204
320 G4VParticleChange* << 205 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
321 G4HadronicProcess::PostStepDoIt(const G4Track& << 206 G4Material *aMaterial = aTrack.GetMaterial();
322 { << 207 G4double originalEnergy = aParticle->GetKineticEnergy();
323 theNumberOfInteractionLengthLeft = -1.0; << 208 G4double kineticEnergy = originalEnergy;
>> 209
>> 210 // More debugging
>> 211
>> 212 G4Nancheck go_wild;
>> 213 if(go_wild(originalEnergy) ||
>> 214 go_wild(aParticle->Get4Momentum().x()) ||
>> 215 go_wild(aParticle->Get4Momentum().y()) ||
>> 216 go_wild(aParticle->Get4Momentum().z()) ||
>> 217 go_wild(aParticle->Get4Momentum().t())
>> 218 )
>> 219 {
>> 220 G4Exception("G4HadronicProcess", "001", JustWarning, "NaN in input energy or momentum - bailing out.");
>> 221 theTotalResult->Clear();
>> 222 theTotalResult->Initialize(aTrack);
>> 223 return theTotalResult;
>> 224 }
324 225
325 //G4cout << "PostStepDoIt " << aTrack.GetDef << 226 // Get kinetic energy per nucleon for ions
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 227
334 // Find cross section at end of step and che << 228 if(aParticle->GetDefinition()->GetBaryonNumber() > 1.5)
335 // << 229 kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
336 const G4DynamicParticle* aParticle = aTrack. <<
337 const G4Material* aMaterial = aTrack.GetMate <<
338 230
339 // check only for charged particles << 231 G4Element* anElement = 0;
340 if(fXSType != fHadNoIntegral) { << 232 try
341 mfpKinEnergy = DBL_MAX; << 233 {
342 G4double xs = aScaleFactor* << 234 anElement = ChooseAandZ( aParticle, aMaterial );
343 theCrossSectionDataStore->ComputeCrossSe << 235 }
344 //G4cout << "xs=" << xs << " xs0=" << theL << 236 catch(G4HadronicException & aR)
345 // << " " << aMaterial->GetName() << << 237 {
346 if(xs < theLastCrossSection*G4UniformRand( << 238 aR.Report(G4cout);
347 // No interaction << 239 G4cout << "Unrecoverable error for:"<<G4endl;
348 return theTotalResult; << 240 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
349 } << 241 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
>> 242 G4cout << " - Particle type = "
>> 243 <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
>> 244 G4Exception("G4HadronicProcess", "007", FatalException,
>> 245 "GeneralPostStepDoIt failed on element selection.");
350 } 246 }
351 247
352 const G4Element* anElement = << 248 try
353 theCrossSectionDataStore->SampleZandA(aPar << 249 {
354 << 250 theInteraction = ChooseHadronicInteraction( kineticEnergy,
355 // Next check for illegal track status << 251 aMaterial, anElement );
356 // << 252 }
357 if (aTrack.GetTrackStatus() != fAlive && << 253 catch(G4HadronicException & aE)
358 aTrack.GetTrackStatus() != fSuspend) { << 254 {
359 if (aTrack.GetTrackStatus() == fStopAndKil << 255 aE.Report(std::cout);
360 aTrack.GetTrackStatus() == fKillTrackA << 256 G4cout << "Unrecoverable error for:"<<G4endl;
361 aTrack.GetTrackStatus() == fPostponeTo << 257 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
362 G4ExceptionDescription ed; << 258 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
363 ed << "G4HadronicProcess: track in unusa << 259 G4cout << " - Particle type = "
364 << aTrack.GetTrackStatus() << G4endl; << 260 << aParticle->GetDefinition()->GetParticleName()<<G4endl;
365 ed << "G4HadronicProcess: returning unch << 261 G4Exception("G4HadronicProcess", "007", FatalException,
366 DumpState(aTrack,"PostStepDoIt",ed); << 262 "ChooseHadronicInteraction failed.");
367 G4Exception("G4HadronicProcess::PostStep <<
368 } <<
369 // No warning for fStopButAlive which is a <<
370 return theTotalResult; <<
371 } 263 }
372 264
373 // Initialize the hadronic projectile from t 265 // Initialize the hadronic projectile from the track
374 thePro.Initialise(aTrack); <<
375 <<
376 theInteraction = ChooseHadronicInteraction(t <<
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 266
390 G4HadFinalState* result = nullptr; << 267 G4HadProjectile thePro(aTrack);
>> 268
>> 269 G4HadFinalState* result = 0;
391 G4int reentryCount = 0; 270 G4int reentryCount = 0;
392 /* << 271
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 272 do
401 { 273 {
402 try 274 try
403 { 275 {
404 // Call the interaction 276 // Call the interaction
405 result = theInteraction->ApplyYourself( << 277
406 ++reentryCount; << 278 G4HadronicInteractionWrapper aW;
>> 279 result = aW.ApplyInteraction(thePro, targetNucleus, theInteraction,
>> 280 GetProcessName(),
>> 281 theInteraction->GetModelName());
407 } 282 }
408 catch(G4HadronicException & aR) << 283 catch(G4HadReentrentException aR)
409 { 284 {
410 G4ExceptionDescription ed; << 285 aR.Report(G4cout);
411 aR.Report(ed); << 286 G4cout << " G4HadronicProcess re-entering the ApplyYourself call for "
412 ed << "Call for " << theInteraction->Get << 287 <<G4endl;
413 ed << "Target element "<<anElement->GetN << 288 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
414 << targetNucleus.GetZ_asInt() << 289 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
415 << " A= " << targetNucleus.GetA_asInt() << << 290 G4cout << " - Particle type = "
416 DumpState(aTrack,"ApplyYourself",ed); << 291 << aParticle->GetDefinition()->GetParticleName() << G4endl;
417 ed << " ApplyYourself failed" << G4endl; << 292 result = 0; // here would still be leaking...
418 G4Exception("G4HadronicProcess::PostStep << 293 if(reentryCount>100)
419 ed); << 294 {
420 } << 295 G4Exception("G4HadronicProcess", "007", FatalException,
421 << 296 "GetHadronicProcess: Reentering ApplyYourself too often - GeneralPostStepDoIt failed.");
422 // Check the result for catastrophic energ <<
423 result = CheckResult(thePro, targetNucleus <<
424 <<
425 if(reentryCount>100) { <<
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 } 297 }
>> 298 G4Exception("G4HadronicProcess", "007", FatalException,
>> 299 "GetHadronicProcess: GeneralPostStepDoIt failed (Reentering ApplyYourself not yet supported.)");
>> 300 }
>> 301 catch(G4HadronicException aR)
>> 302 {
>> 303 aR.Report(G4cout);
>> 304 G4cout << " G4HadronicProcess failed in ApplyYourself call for"
>> 305 << G4endl;
>> 306 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
>> 307 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
>> 308 G4cout << " - Particle type = "
>> 309 << aParticle->GetDefinition()->GetParticleName() << G4endl;
>> 310 G4Exception("G4HadronicProcess", "007", FatalException,
>> 311 "GeneralPostStepDoIt failed.");
466 } 312 }
467 } 313 }
>> 314 while(!result);
468 315
469 result->SetTrafoToLab(thePro.GetTrafoToLab() << 316 if(!ModelingState && !getenv("BypassAllSafetyChecks") )
470 FillResult(result, aTrack); << 317 {
471 << 318 G4cout << "ERROR IN EXECUTION -- HADRONIC PROCESS STATE NOT VALID"<<G4endl;
472 if (epReportLevel != 0) { << 319 G4cout << "Result will be of undefined quality."<<G4endl;
473 CheckEnergyMomentumConservation(aTrack, ta <<
474 } 320 }
475 //G4cout << "PostStepDoIt done nICelectrons= <<
476 return theTotalResult; <<
477 } <<
478 321
479 void G4HadronicProcess::ProcessDescription(std << 322 // NOT USED ?? Projectile particle has changed character during interaction
480 { << 323 if(result->GetStatusChange() == isAlive &&
481 outFile << "The description for this process << 324 thePro.GetDefinition() != aTrack.GetDefinition())
482 } << 325 {
483 << 326 G4DynamicParticle * aP =
484 G4double G4HadronicProcess::XBiasSurvivalProba << 327 const_cast<G4DynamicParticle *>(aTrack.GetDynamicParticle());
485 { << 328 aP->SetDefinition(const_cast<G4ParticleDefinition *>(thePro.GetDefinition()));
486 G4double nLTraversed = GetTotalNumberOfInter << 329 }
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 330
507 G4double efinal = std::max(aR->GetEnergyChan << 331 result->SetTrafoToLab(thePro.GetTrafoToLab());
508 332
509 // check status of primary << 333 /*
510 if(aR->GetStatusChange() == stopAndKill) { << 334 // Loop over charged ion secondaries
511 theTotalResult->ProposeTrackStatus(fStopAn <<
512 theTotalResult->ProposeEnergy( 0.0 ); <<
513 335
514 // check its final energy << 336 for(G4int i=0; i<result->GetNumberOfSecondaries(); i++)
515 } else if(0.0 == efinal) { << 337 {
516 theTotalResult->ProposeEnergy( 0.0 ); << 338 G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle();
517 if(aT.GetParticleDefinition()->GetProcessM << 339 if(aSecTrack->GetDefinition()->GetPDGCharge()>1.5)
518 ->GetAtRestProcessVector()->size() > 0) << 340 {
519 { theTotalResult->ProposeTrackStatus( << 341 G4EffectiveCharge aCalculator;
520 else { theTotalResult->ProposeTrackStatus( << 342 G4double charge =
521 << 343 aCalculator.GetCharge(aMaterial, aSecTrack->GetKineticEnergy(),
522 // primary is not killed apply rotation an << 344 aSecTrack->GetDefinition()->GetPDGMass(),
523 } else { << 345 aSecTrack->GetDefinition()->GetPDGCharge());
524 theTotalResult->ProposeTrackStatus(fAlive) << 346 if(getenv("GHADChargeDebug"))
525 G4ThreeVector newDir = aR->GetMomentumChan << 347 {
526 newDir.rotateUz(dir); << 348 std::cout << "Recoil fractional charge is "
527 theTotalResult->ProposeMomentumDirection(n << 349 << charge/aSecTrack->GetDefinition()->GetPDGCharge()<<" "
528 theTotalResult->ProposeEnergy(efinal); << 350 << charge <<" "<<aSecTrack->GetDefinition()->GetPDGCharge()<<std::endl;
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 } 351 }
569 dynParticle->SetKineticEnergy(e); << 352 aSecTrack->SetCharge(charge);
570 dynParticle->SetMass(mass); <<
571 } 353 }
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 } 354 }
587 aR->Clear(); << 355 */
588 // G4cout << "FillResults done nICe= " << nI <<
589 } <<
590 <<
591 void G4HadronicProcess::MultiplyCrossSectionBy <<
592 { <<
593 BiasCrossSectionByFactor(factor); <<
594 } <<
595 356
596 void G4HadronicProcess::BiasCrossSectionByFact << 357 if(getenv("HadronicDoitLogging") )
597 { << 358 {
598 if (aScale <= 0.0) { << 359 G4cout << "HadronicDoitLogging "
599 G4ExceptionDescription ed; << 360 << GetProcessName() <<" "
600 ed << " Wrong biasing factor " << aScale < << 361 << aParticle->GetDefinition()->GetPDGEncoding()<<" "
601 G4Exception("G4HadronicProcess::BiasCrossS << 362 << originalEnergy<<" "
602 JustWarning, ed, "Cross-sectio << 363 << aParticle->GetMomentum()<<" "
603 } else { << 364 << targetNucleus.GetN()<<" "
604 aScaleFactor = aScale; << 365 << targetNucleus.GetZ()<<" "
>> 366 << G4endl;
605 } 367 }
606 } <<
607 368
608 G4HadFinalState* G4HadronicProcess::CheckResul << 369 ClearNumberOfInteractionLengthLeft();
609 const G4Nucleus &aNucleus, << 370 if(isoIsOnAnyway!=-1)
610 G4HadFinalState * result) << 371 {
611 { << 372 if(isoIsEnabled||isoIsOnAnyway)
612 // check for catastrophic energy non-conserv << 373 {
613 // to re-sample the interaction << 374 result = DoIsotopeCounting(result, aTrack, targetNucleus);
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 } 375 }
670 G4double deltaE= nuclearMass + aPro.GetTo << 376 }
671 << 377
672 std::pair<G4double, G4double> checkLevels << 378 G4double e=aTrack.GetKineticEnergy();
673 theModel->GetFatalEnergyCheckLevels(); << 379 ModelingState = 0;
674 if (std::abs(deltaE) > checkLevels.second << 380 if(e<5*GeV)
675 std::abs(deltaE) > checkLevels.first*a << 381 {
676 // do not delete result, this is a point << 382 for(size_t i=0; i<theBias.size(); i++)
677 result->Clear(); << 383 {
678 result = nullptr; << 384 result = theBias[i]->Bias(result);
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 } 385 }
694 } 386 }
695 return result; << 387
>> 388 // Put hadronic final state particles into G4ParticleChange
>> 389
>> 390 FillTotalResult(result, aTrack);
>> 391 if(G4HadronicProcess_debug_flag)
>> 392 std::cout << "@@@@ hadronic process end "<< std::endl;
>> 393
>> 394 return theTotalResult;
696 } 395 }
697 396
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 } <<
754 } <<
755 <<
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 397
858 if (aTrack.GetMaterial()) { << 398 G4HadFinalState*
859 ed << " material " << aTrack.GetMaterial( << 399 G4HadronicProcess::DoIsotopeCounting(G4HadFinalState * aResult,
>> 400 const G4Track & aTrack,
>> 401 const G4Nucleus & aNucleus)
>> 402 {
>> 403 // get the PC from iso-production
>> 404 delete theOldIsoResult;
>> 405 theOldIsoResult = 0;
>> 406 delete theIsoResult;
>> 407 theIsoResult = new G4IsoParticleChange;
>> 408 G4bool done = false;
>> 409 G4IsoResult * anIsoResult = 0;
>> 410 for(unsigned int i=0; i<theProductionModels.size(); i++)
>> 411 {
>> 412 anIsoResult = theProductionModels[i]->GetIsotope(aTrack, aNucleus);
>> 413 if(anIsoResult!=0)
>> 414 {
>> 415 done = true;
>> 416 break;
>> 417 }
860 } 418 }
861 ed << G4endl; <<
862 419
863 if (aTrack.GetVolume()) { << 420 // If no production models active, use default iso production
864 ed << "PhysicalVolume <" << aTrack.GetVol << 421 if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult);
865 << ">" << G4endl; <<
866 } <<
867 } <<
868 422
869 void G4HadronicProcess::DumpPhysicsTable(const << 423 // Add all info explicitely and add typename from model called.
870 { << 424 theIsoResult->SetIsotope(anIsoResult->GetIsotope());
871 theCrossSectionDataStore->DumpPhysicsTable(p << 425 theIsoResult->SetProductionPosition(aTrack.GetPosition());
872 } << 426 theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
>> 427 theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
>> 428 theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
>> 429 theIsoResult->SetProducer(typeid(*theInteraction).name());
>> 430
>> 431 delete anIsoResult;
873 432
874 void G4HadronicProcess::AddDataSet(G4VCrossSec << 433 // If isotope production is enabled the GetIsotopeProductionInfo()
875 { << 434 // method must be called or else a memory leak will result
876 theCrossSectionDataStore->AddDataSet(aDataSe << 435 //
877 } << 436 // The following code will fix the memory leak, but remove the
>> 437 // isotope information:
>> 438 //
>> 439 // if(theIsoResult) {
>> 440 // delete theIsoResult;
>> 441 // theIsoResult = 0;
>> 442 // }
>> 443
>> 444 return aResult;
>> 445 }
>> 446
>> 447 G4IsoResult*
>> 448 G4HadronicProcess::ExtractResidualNucleus(const G4Track&,
>> 449 const G4Nucleus& aNucleus,
>> 450 G4HadFinalState* aResult)
>> 451 {
>> 452 G4double A = aNucleus.GetN();
>> 453 G4double Z = aNucleus.GetZ();
>> 454 G4double bufferA = 0;
>> 455 G4double bufferZ = 0;
>> 456
>> 457 // loop over aResult, and decrement A, Z accordingly
>> 458 // cash the max
>> 459 for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
>> 460 {
>> 461 G4HadSecondary* aSecTrack = aResult->GetSecondary(i);
>> 462 if(bufferA<aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber())
>> 463 {
>> 464 bufferA = aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
>> 465 bufferZ = aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
>> 466 }
>> 467 Z-=aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
>> 468 A-=aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
>> 469 }
>> 470
>> 471 // if the fragment was part of the final state, it is
>> 472 // assumed to be the heaviest secondary.
>> 473 if(A<0.1)
>> 474 {
>> 475 A = bufferA;
>> 476 Z = bufferZ;
>> 477 }
>> 478
>> 479 // prepare the IsoResult.
878 480
879 std::vector<G4HadronicInteraction*>& << 481 std::ostringstream ost1;
880 G4HadronicProcess::GetHadronicInteractionList( << 482 ost1 <<Z<<"_"<<A;
881 { << 483 G4String biff = ost1.str();
882 return theEnergyRangeManager.GetHadronicInte << 484 G4IsoResult * theResult = new G4IsoResult(biff, aNucleus);
883 } <<
884 485
885 G4HadronicInteraction* << 486 return theResult;
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 } 487 }
895 488
896 G4double << 489 G4double G4HadronicProcess::XBiasSurvivalProbability()
897 G4HadronicProcess::ComputeCrossSection(const G <<
898 const G4Material* mat, <<
899 const G4double kinEnergy) <<
900 { 490 {
901 auto dp = new G4DynamicParticle(part, unitVe << 491 G4double result = 0;
902 G4double xs = theCrossSectionDataStore->Comp << 492 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
903 delete dp; << 493 G4double biasedProbability = 1.-std::exp(-nLTraversed);
904 return xs; << 494 G4double realProbability = 1-std::exp(-nLTraversed/aScaleFactor);
>> 495 result = (biasedProbability-realProbability)/biasedProbability;
>> 496 return result;
905 } 497 }
906 498
907 void G4HadronicProcess::RecomputeXSandMFP(cons << 499 G4double G4HadronicProcess::XBiasSecondaryWeight()
908 { 500 {
909 auto dp = new G4DynamicParticle(currentParti << 501 G4double result = 0;
910 theLastCrossSection = aScaleFactor* << 502 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
911 theCrossSectionDataStore->ComputeCrossSect << 503 result =
912 theMFP = (theLastCrossSection > 0.0) ? 1.0/t << 504 1./aScaleFactor*std::exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
913 delete dp; << 505 return result;
914 } 506 }
915 507
916 void G4HadronicProcess::UpdateCrossSectionAndM << 508 void
>> 509 G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
917 { 510 {
918 if(fXSType == fHadNoIntegral) { << 511 G4Nancheck go_wild;
919 DefineXSandMFP(); << 512 theTotalResult->Clear();
>> 513 theTotalResult->ProposeLocalEnergyDeposit(0.);
>> 514 theTotalResult->Initialize(aT);
>> 515 theTotalResult->SetSecondaryWeightByProcess(true);
>> 516 theTotalResult->ProposeTrackStatus(fAlive);
>> 517 G4double rotation = 2.*pi*G4UniformRand();
>> 518 G4ThreeVector it(0., 0., 1.);
920 519
921 } else if(fXSType == fHadIncreasing) { << 520 /*
922 if(e*invLambdaFactor < mfpKinEnergy) { << 521 if(xBiasOn)
923 mfpKinEnergy = e; << 522 {
924 ComputeXSandMFP(); << 523 G4cout << "BiasDebug "<<GetProcessName()<<" "
925 } << 524 <<aScaleFactor<<" "
926 << 525 <<XBiasSurvivalProbability()<<" "
927 } else if(fXSType == fHadDecreasing) { << 526 <<XBiasSecondaryWeight()<<" "
928 if(e < mfpKinEnergy && mfpKinEnergy > minK << 527 <<G4endl;
929 G4double e1 = std::max(e*lambdaFactor, m << 528 }
930 mfpKinEnergy = e1; << 529 */
931 RecomputeXSandMFP(e1); << 530 // if(GetProcessName() != "LElastic") std::cout << "Debug -1 "<<aR->GetStatusChange()<<std::endl;
932 } << 531 if(aR->GetStatusChange()==stopAndKill)
933 << 532 {
934 } else if(fXSType == fHadOnePeak) { << 533 if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
935 G4double epeak = (*theEnergyOfCrossSection << 534 {
936 if(e <= epeak) { << 535 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
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 } 536 }
946 << 537 else
947 } else if(fXSType == fHadTwoPeaks) { << 538 {
948 G4TwoPeaksHadXS* xs = (*fXSpeaks)[matIdx]; << 539 theTotalResult->ProposeTrackStatus(fStopAndKill);
949 const G4double e1peak = xs->e1peak; << 540 theTotalResult->ProposeEnergy( 0.0 );
950 <<
951 // below the 1st peak <<
952 if(e <= e1peak) { <<
953 if(e*invLambdaFactor < mfpKinEnergy) { <<
954 mfpKinEnergy = e; <<
955 ComputeXSandMFP(); <<
956 } <<
957 return; <<
958 } 541 }
959 const G4double e1deep = xs->e1deep; << 542 }
960 // above the 1st peak, below the deep << 543 else if(aR->GetStatusChange()!=stopAndKill )
961 if(e <= e1deep) { << 544 {
962 if(mfpKinEnergy >= e1deep || e <= mfpKin << 545 if(aR->GetStatusChange()==suspend)
963 const G4double e1 = std::max(e1peak, e << 546 {
964 mfpKinEnergy = e1; << 547 theTotalResult->ProposeTrackStatus(fSuspend);
965 RecomputeXSandMFP(e1); << 548 if(xBiasOn)
>> 549 {
>> 550 G4Exception("G4HadronicProcess", "007", FatalException,
>> 551 "Cannot cross-section bias a process that suspends tracks.");
966 } 552 }
967 return; << 553 } else if (aT.GetKineticEnergy() == 0) {
>> 554 theTotalResult->ProposeTrackStatus(fStopButAlive);
968 } 555 }
969 const G4double e2peak = xs->e2peak; << 556
970 // above the deep, below 2nd peak << 557 if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
971 if(e <= e2peak) { << 558 {
972 if(e*invLambdaFactor < mfpKinEnergy) { << 559 theTotalResult->ProposeParentWeight( XBiasSurvivalProbability()*aT.GetWeight() );
973 mfpKinEnergy = e; << 560 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
974 ComputeXSandMFP(); << 561 if(go_wild(aR->GetEnergyChange()))
975 } << 562 {
976 return; << 563 G4Exception("G4HadronicProcess", "007", FatalException,
>> 564 "surviving track received NaN energy.");
>> 565 }
>> 566 if(go_wild(aR->GetMomentumChange().x()) ||
>> 567 go_wild(aR->GetMomentumChange().y()) ||
>> 568 go_wild(aR->GetMomentumChange().z()))
>> 569 {
>> 570 G4Exception("G4HadronicProcess", "007", FatalException,
>> 571 "surviving track received NaN momentum.");
>> 572 }
>> 573 G4double newM=aT.GetDefinition()->GetPDGMass();
>> 574 G4double newE=aR->GetEnergyChange() + newM;
>> 575 G4double newP=std::sqrt(newE*newE - newM*newM);
>> 576 G4DynamicParticle * aNew =
>> 577 new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
>> 578 G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
>> 579 aR->AddSecondary(theSec);
977 } 580 }
978 const G4double e2deep = xs->e2deep; << 581 else
979 // above the 2nd peak, below the deep << 582 {
980 if(e <= e2deep) { << 583 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
981 if(mfpKinEnergy >= e2deep || e <= mfpKin << 584 theTotalResult->ProposeParentWeight(newWeight); // This is multiplicative
982 const G4double e1 = std::max(e2peak, e << 585 if(aR->GetEnergyChange()>-.5)
983 mfpKinEnergy = e1; << 586 {
984 RecomputeXSandMFP(e1); << 587 if(go_wild(aR->GetEnergyChange()))
>> 588 {
>> 589 G4Exception("G4HadronicProcess", "007", FatalException,
>> 590 "track received NaN energy.");
>> 591 }
>> 592 theTotalResult->ProposeEnergy(aR->GetEnergyChange());
985 } 593 }
986 return; << 594 G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
>> 595 newDirection*=aR->GetTrafoToLab();
>> 596 theTotalResult->ProposeMomentumDirection(newDirection.vect());
987 } 597 }
988 const G4double e3peak = xs->e3peak; << 598 }
989 // above the deep, below 3d peak << 599 else
990 if(e <= e3peak) { << 600 {
991 if(e*invLambdaFactor < mfpKinEnergy) { << 601 G4cerr << "Track status is "<< aR->GetStatusChange()<<G4endl;
992 mfpKinEnergy = e; << 602 G4Exception("G4HadronicProcess", "007", FatalException,
993 ComputeXSandMFP(); << 603 "use of unsupported track-status.");
994 } << 604 }
995 return; << 605
>> 606 if(GetProcessName() != "hElastic" && GetProcessName() != "HadronElastic"
>> 607 && theTotalResult->GetTrackStatus()==fAlive
>> 608 && aR->GetStatusChange()==isAlive
>> 609 )
>> 610 {
>> 611 // Use for debugging: G4double newWeight = theTotalResult->GetParentWeight();
>> 612
>> 613 G4double newKE = std::max(DBL_MIN, aR->GetEnergyChange());
>> 614 G4DynamicParticle* aNew = new G4DynamicParticle(aT.GetDefinition(),
>> 615 aR->GetMomentumChange(),
>> 616 newKE);
>> 617 G4HadSecondary* theSec = new G4HadSecondary(aNew, 1.0);
>> 618 aR->AddSecondary(theSec);
>> 619 aR->SetStatusChange(stopAndKill);
>> 620
>> 621 theTotalResult->ProposeTrackStatus(fStopAndKill);
>> 622 theTotalResult->ProposeEnergy( 0.0 );
>> 623
>> 624 }
>> 625 theTotalResult->ProposeLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
>> 626 theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
>> 627
>> 628 if(aR->GetStatusChange() != stopAndKill)
>> 629 {
>> 630 G4double newM=aT.GetDefinition()->GetPDGMass();
>> 631 G4double newE=aR->GetEnergyChange() + newM;
>> 632 G4double newP=std::sqrt(newE*newE - newM*newM);
>> 633 G4ThreeVector newPV = newP*aR->GetMomentumChange();
>> 634 G4LorentzVector newP4(newE, newPV);
>> 635 newP4.rotate(rotation, it);
>> 636 newP4*=aR->GetTrafoToLab();
>> 637 theTotalResult->ProposeMomentumDirection(newP4.vect().unit());
>> 638 }
>> 639
>> 640 for(G4int i=0; i<aR->GetNumberOfSecondaries(); i++)
>> 641 {
>> 642 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
>> 643 theM.rotate(rotation, it);
>> 644 theM*=aR->GetTrafoToLab();
>> 645
>> 646 if(go_wild(theM.e()))
>> 647 {
>> 648 G4Exception("G4HadronicProcess", "007", FatalException,
>> 649 "secondary track received NaN energy.");
996 } 650 }
997 // above 3d peak << 651 if(go_wild(theM.x()) ||
998 if(e <= mfpKinEnergy) { << 652 go_wild(theM.y()) ||
999 const G4double e1 = std::max(e3peak, e*l << 653 go_wild(theM.z()))
1000 mfpKinEnergy = e1; << 654 {
1001 RecomputeXSandMFP(e1); << 655 G4Exception("G4HadronicProcess", "007", FatalException,
>> 656 "secondary track received NaN momentum.");
1002 } 657 }
1003 658
1004 } else { << 659 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
1005 DefineXSandMFP(); << 660 G4double time = aR->GetSecondary(i)->GetTime();
1006 } << 661 if(time<0) time = aT.GetGlobalTime();
>> 662
>> 663 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
>> 664 time,
>> 665 aT.GetPosition());
>> 666
>> 667 G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
>> 668 //static G4double pinelcount=0;
>> 669 if(xBiasOn) newWeight *= XBiasSecondaryWeight();
>> 670 // G4cout << "#### ParticleDebug "
>> 671 // <<GetProcessName()<<" "
>> 672 // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
>> 673 // <<aScaleFactor<<" "
>> 674 // <<XBiasSurvivalProbability()<<" "
>> 675 // <<XBiasSecondaryWeight()<<" "
>> 676 // <<aT.GetWeight()<<" "
>> 677 // <<aR->GetSecondary(i)->GetWeight()<<" "
>> 678 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
>> 679 // <<G4endl;
>> 680 track->SetWeight(newWeight);
>> 681 G4double trackDeb = track->GetKineticEnergy();
>> 682 if( ( trackDeb<0
>> 683 || (trackDeb>aT.GetKineticEnergy()+1*GeV) ) && getenv("GHADEnergyBalanceDebug") )
>> 684 {
>> 685 G4cout << "Debugging hadronic processes: "<<track->GetKineticEnergy()
>> 686 <<" "<<aT.GetKineticEnergy()
>> 687 <<" "<<GetProcessName()
>> 688 <<" "<<aT.GetDefinition()->GetParticleName()
>> 689 <<G4endl;
>> 690 }
>> 691 track->SetTouchableHandle(aT.GetTouchableHandle());
>> 692 theTotalResult->AddSecondary(track);
>> 693 }
>> 694
>> 695 aR->Clear();
>> 696 return;
1007 } 697 }
>> 698 /* end of file */
1008 699