| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
27 // ------------------------------------------- <<
28 // 24 //
29 // GEANT4 Class source file << 25 // HPW to implement the choosing of an element for scattering.
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 26
52 #include "G4Types.hh" 27 #include "G4Types.hh"
53 #include "G4SystemOfUnits.hh" << 28
>> 29 #include <fstream>
>> 30 #include <strstream>
>> 31 #include <stdlib.h>
>> 32 #include "G4HadronicProcess.hh"
>> 33 #include "G4EffectiveCharge.hh"
54 #include "G4HadProjectile.hh" 34 #include "G4HadProjectile.hh"
55 #include "G4ElementVector.hh" 35 #include "G4ElementVector.hh"
56 #include "G4Track.hh" 36 #include "G4Track.hh"
57 #include "G4Step.hh" 37 #include "G4Step.hh"
58 #include "G4Element.hh" 38 #include "G4Element.hh"
59 #include "G4ParticleChange.hh" 39 #include "G4ParticleChange.hh"
>> 40 #include "G4TransportationManager.hh"
>> 41 #include "G4Navigator.hh"
60 #include "G4ProcessVector.hh" 42 #include "G4ProcessVector.hh"
61 #include "G4ProcessManager.hh" 43 #include "G4ProcessManager.hh"
62 #include "G4NucleiProperties.hh" << 44 #include "G4StableIsotopes.hh"
>> 45 #include "G4HadTmpUtil.hh"
63 46
>> 47 #include "G4HadLeadBias.hh"
64 #include "G4HadronicException.hh" 48 #include "G4HadronicException.hh"
65 #include "G4HadronicProcessStore.hh" << 49 #include "G4HadReentrentException.hh"
66 #include "G4HadronicParameters.hh" << 50 #include "G4HadronicInteractionWrapper.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 51
79 namespace << 52 #include "G4HadSignalHandler.hh"
80 { <<
81 constexpr G4double lambdaFactor = 0.8; <<
82 constexpr G4double invLambdaFactor = 1.0/lam <<
83 } <<
84 <<
85 ////////////////////////////////////////////// <<
86 53
87 G4HadronicProcess::G4HadronicProcess(const G4S << 54 //@@ add model name info, once typeinfo available #include <typeinfo.h>
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 55
103 G4HadronicProcess::~G4HadronicProcess() << 56 namespace G4HadronicProcess_local
104 { << 57 {
105 theProcessStore->DeRegister(this); << 58 extern "C" void G4HadronicProcessHandler_1(int)
106 delete theTotalResult; << 59 {
107 delete theCrossSectionDataStore; << 60 G4HadronicWhiteBoard::Instance().Dump();
108 if(isMaster) { <<
109 if (fXSpeaks != nullptr) { <<
110 for (auto const& e : *fXSpeaks ) { <<
111 delete e; <<
112 } <<
113 } 61 }
114 delete fXSpeaks; << 62 }
115 delete theEnergyOfCrossSectionMax; << 63
116 } << 64 G4IsoParticleChange * G4HadronicProcess::theIsoResult = NULL;
117 } << 65 G4IsoParticleChange * G4HadronicProcess::theOldIsoResult = NULL;
118 << 66 G4bool G4HadronicProcess::isoIsEnabled = true;
119 void G4HadronicProcess::InitialiseLocal() { << 67
120 theTotalResult = new G4ParticleChange(); << 68 void G4HadronicProcess::
121 theTotalResult->SetSecondaryWeightByProcess( << 69 EnableIsotopeProductionGlobally() {isoIsEnabled = true;}
122 theCrossSectionDataStore = new G4CrossSectio << 70
123 theProcessStore = G4HadronicProcessStore::In << 71 void G4HadronicProcess::
124 theProcessStore->Register(this); << 72 DisableIsotopeProductionGlobally() {isoIsEnabled = false;}
125 minKinEnergy = 1*CLHEP::MeV; << 73
126 << 74 G4HadronicProcess::G4HadronicProcess( const G4String &processName) :
127 G4HadronicParameters* param = G4HadronicPara << 75 G4VDiscreteProcess( processName )
128 epReportLevel = param->GetEPReportLevel(); << 76 {
129 epCheckLevels.first = param->GetEPRelativeLe << 77 ModelingState = 0;
130 epCheckLevels.second = param->GetEPAbsoluteL << 78 isoIsOnAnyway = 0;
131 << 79 theTotalResult = new G4ParticleChange();
132 unitVector.set(0.0, 0.0, 0.1); << 80 theCrossSectionDataStore = new G4CrossSectionDataStore();
133 if(G4Threading::IsWorkerThread()) { isMaster << 81 aScaleFactor = 1;
134 } << 82 xBiasOn = false;
135 << 83 if(getenv("SwitchLeadBiasOn")) theBias.push_back(new G4HadLeadBias());
136 void G4HadronicProcess::RegisterMe( G4Hadronic << 84 }
137 { << 85
138 if(nullptr == a) { return; } << 86 G4HadronicProcess::~G4HadronicProcess()
139 theEnergyRangeManager.RegisterMe( a ); << 87 {
140 G4HadronicProcessStore::Instance()->Register << 88 delete theTotalResult;
141 } << 89 std::for_each(theProductionModels.begin(),
>> 90 theProductionModels.end(),
>> 91 G4Delete());
>> 92 std::for_each(theBias.begin(),
>> 93 theBias.end(),
>> 94 G4Delete());
>> 95 if(theOldIsoResult) delete theOldIsoResult;
>> 96 // if(theIsoResult) delete theIsoResult;
>> 97 }
>> 98
>> 99 void G4HadronicProcess::RegisterMe( G4HadronicInteraction *a )
>> 100 {
>> 101 try{GetManagerPointer()->RegisterMe( a );}
>> 102 catch(G4HadronicException & aE)
>> 103 {
>> 104 aE.Report(std::cout);
>> 105 G4Exception("G4HadronicProcess", "007", FatalException,
>> 106 "Could not register G4HadronicInteraction");
>> 107 }
>> 108 }
142 109
143 G4double << 110 G4double G4HadronicProcess::
144 G4HadronicProcess::GetElementCrossSection(cons << 111 GetMeanFreePath(const G4Track &aTrack, G4double, G4ForceCondition *)
145 const G4Element * elm, << 112 {
146 const G4Material* mat) << 113 /*
147 { << 114 if(ReStarted)
148 if(nullptr == mat) <<
149 { 115 {
150 static const G4int nmax = 5; << 116 if(trackIdCache == aTrack.GetTrackId())
151 if(nMatWarn < nmax) { << 117 {
152 ++nMatWarn; << 118 theInitialNumberOfInteractionLength += G4VProcess::theNumberOfInteractionLengthLeft;
153 G4ExceptionDescription ed; << 119
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 } 120 }
>> 121 else
>> 122 {
>> 123 theInitialNumberOfInteractionLength = G4VProcess::theNumberOfInteractionLengthLeft;
>> 124 }
>> 125 trackIdCache = aTrack.GetTrackId();
>> 126 ReStarted = false;
161 } 127 }
162 return theCrossSectionDataStore->GetCrossSec << 128 */
163 } << 129 G4double sigma = 0.0;
164 << 130 try
165 void G4HadronicProcess::PreparePhysicsTable(co << 131 {
166 { << 132 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
167 if(nullptr == firstParticle) { firstParticle << 133 if( !IsApplicable(*aParticle->GetDefinition()))
168 theProcessStore->RegisterParticle(this, &p); << 134 {
>> 135 G4cout << "Unrecoverable error: "<<G4endl;
>> 136 G4ProcessManager * it = aParticle->GetDefinition()->GetProcessManager();
>> 137 G4ProcessVector * itv = it->GetProcessList();
>> 138 G4cout <<aParticle->GetDefinition()->GetParticleName()<<
>> 139 " has the following processes:"<<G4endl;
>> 140 for(G4int i=0; i<itv->size(); i++)
>> 141 {
>> 142 G4cout <<" "<<(*itv)[i]->GetProcessName()<<G4endl;
>> 143 }
>> 144 G4cout << "for kinetic energy "<<aParticle->GetKineticEnergy()<<G4endl;
>> 145 G4cout << "and material "<<aTrack.GetMaterial()->GetName()<<G4endl;
>> 146 G4Exception("G4HadronicProcess", "007", FatalException,
>> 147 std::string(this->GetProcessName()+
>> 148 " was called for "+
>> 149 aParticle->GetDefinition()->GetParticleName()).c_str() );
>> 150 }
>> 151 G4Material *aMaterial = aTrack.GetMaterial();
>> 152 G4int nElements = aMaterial->GetNumberOfElements();
>> 153 ModelingState = 1;
>> 154
>> 155 // returns the mean free path in GEANT4 internal units
>> 156
>> 157 const G4double *theAtomicNumDensityVector =
>> 158 aMaterial->GetAtomicNumDensityVector();
>> 159
>> 160 G4double aTemp = aMaterial->GetTemperature();
>> 161
>> 162 for( G4int i=0; i<nElements; ++i )
>> 163 {
>> 164 G4double xSection =
>> 165 GetMicroscopicCrossSection( aParticle, (*aMaterial->GetElementVector())[i], aTemp);
>> 166 sigma += theAtomicNumDensityVector[i] * xSection;
>> 167 }
>> 168 sigma *= aScaleFactor;
>> 169 theLastCrossSection = sigma;
>> 170 }
>> 171 catch(G4HadronicException aR)
>> 172 {
>> 173 aR.Report(G4cout);
>> 174 G4Exception("G4HadronicProcess", "007", FatalException,
>> 175 "G4HadronicProcess::GetMeanFreePath failed");
>> 176 }
>> 177 if( sigma > 0.0 )
>> 178 return 1.0/sigma;
>> 179 else
>> 180 return DBL_MAX;
169 } 181 }
170 182
171 void G4HadronicProcess::BuildPhysicsTable(cons << 183 G4double G4HadronicProcess::GetDistanceToBoundary(const G4Track & aT)
172 { << 184 {
173 if(firstParticle != &p) { return; } << 185 G4TransportationManager * aTM =
174 << 186 G4TransportationManager::GetTransportationManager();
175 theCrossSectionDataStore->BuildPhysicsTable( << 187 G4Navigator * aN = aTM->GetNavigatorForTracking();
176 theEnergyRangeManager.BuildPhysicsTable(p); << 188 G4ThreeVector pGlobalPoint = aT.GetStep()->GetPreStepPoint()->GetPosition();
177 G4HadronicParameters* param = G4HadronicPara << 189 G4ThreeVector pDirection = aT.GetMomentumDirection();
178 << 190 G4double dummy(0);
179 G4int subtype = GetProcessSubType(); << 191 G4double result = aN->ComputeStep(pGlobalPoint, pDirection, DBL_MAX, dummy);
180 if(useIntegralXS) { << 192 aN->LocateGlobalPointAndSetup(pGlobalPoint);
181 if(subtype == fHadronInelastic) { << 193 return result;
182 useIntegralXS = param->EnableIntegralIne << 194 }
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 195
205 // check particle for integral method << 196 G4Element * G4HadronicProcess::ChooseAandZ(
206 if(isMaster || nullptr == masterProcess) { << 197 const G4DynamicParticle *aParticle, const G4Material *aMaterial )
207 G4double charge = p.GetPDGCharge()/eplus; << 198 {
208 << 199 static G4bool noIsotopeWiseCrossSections=getenv("GHAD_DISABLE_ISOTOPE_WISE_CROSS_SECTIONS");
209 // select cross section shape << 200 static G4StableIsotopes theIso;
210 if(charge != 0.0 && useIntegralXS) { << 201 currentZ = 0;
211 G4double tmax = param->GetMaxEnergy(); << 202 currentN = 0;
212 currentParticle = firstParticle; << 203 const G4int numberOfElements = aMaterial->GetNumberOfElements();
213 // initialisation in the master thread << 204 const G4ElementVector *theElementVector = aMaterial->GetElementVector();
214 G4int pdg = p.GetPDGEncoding(); << 205 G4int i;
215 if (std::abs(pdg) == 211) { << 206 if( numberOfElements == 1 )
216 fXSType = fHadTwoPeaks; << 207 {
217 } else if (pdg == 321) { << 208 currentZ = G4double( ((*theElementVector)[0])->GetZ());
218 fXSType = fHadOnePeak; << 209 G4int localZ = G4lrint(currentZ);
219 } else if (pdg == -321) { << 210 if(noIsotopeWiseCrossSections)
220 fXSType = fHadDecreasing; << 211 {
221 } else if (pdg == 2212) { << 212 currentN = (*theElementVector)[0]->GetN();
222 fXSType = fHadTwoPeaks; << 213 }
223 } else if (pdg == -2212 || pdg == -10000 << 214 else
224 pdg == -1000020030 || pdg == -1000020040) << 215 {
225 fXSType = fHadDecreasing; << 216 G4double * running = new G4double[theIso.GetNumberOfIsotopes(localZ)];
226 } else if (charge > 0.0 || pdg == 11 || << 217 for (i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
227 fXSType = fHadIncreasing; << 218 {
228 } << 219 G4double fracInPercent=theIso.GetAbundance(theIso.GetFirstIsotope(localZ)+i);
229 << 220 G4double runningA=theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
230 delete theEnergyOfCrossSectionMax; << 221 running[i]=fracInPercent*pow(runningA, 2./3.);
231 theEnergyOfCrossSectionMax = nullptr; << 222 // rough approximation; to get it better, redesign getMSC to not use G4Element, see also below
232 if(fXSType == fHadTwoPeaks) { << 223 if(i!=0) running[i] += running[i-1];
233 if (fXSpeaks != nullptr) { <<
234 for (auto const& e : *fXSpeaks ) { <<
235 delete e; <<
236 } <<
237 } 224 }
238 delete fXSpeaks; << 225 G4double trial = G4UniformRand();
239 fXSpeaks = << 226 G4double sum = running[theIso.GetNumberOfIsotopes(localZ)-1];
240 G4HadXSHelper::FillPeaksStructure(this, &p << 227 for(i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
241 if(nullptr == fXSpeaks) { << 228 {
242 fXSType = fHadOnePeak; << 229 currentN = theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
>> 230 if(running[i]/sum>trial) break;
243 } 231 }
>> 232 delete [] running;
244 } 233 }
245 if(fXSType == fHadOnePeak) { << 234 targetNucleus.SetParameters(currentN, currentZ);
246 theEnergyOfCrossSectionMax = << 235 return (*theElementVector)[0];
247 G4HadXSHelper::FindCrossSectionMax(this, & << 236 }
248 if(nullptr == theEnergyOfCrossSectionMax) { << 237
249 fXSType = fHadIncreasing; << 238 const G4double *theAtomicNumberDensity = aMaterial->GetAtomicNumDensityVector();
>> 239 G4double aTemp = aMaterial->GetTemperature();
>> 240 G4double crossSectionTotal = 0;
>> 241 std::vector<G4double> runningSum;
>> 242 for( i=0; i < numberOfElements; ++i )
>> 243 {
>> 244 runningSum.push_back(theAtomicNumberDensity[i] *
>> 245 dispatch->GetMicroscopicCrossSection( aParticle, (*theElementVector)[i], aTemp));
>> 246 crossSectionTotal+=runningSum[i];
>> 247 }
>> 248
>> 249 G4double random = G4UniformRand();
>> 250 for( i=0; i < numberOfElements; ++i )
>> 251 {
>> 252 if(i!=0) runningSum[i]+=runningSum[i-1];
>> 253 if( random<=runningSum[i]/crossSectionTotal )
>> 254 {
>> 255 currentZ = G4double( ((*theElementVector)[i])->GetZ());
>> 256 G4int localZ = G4lrint(currentZ);
>> 257 if(noIsotopeWiseCrossSections)
>> 258 {
>> 259 currentN = ((*theElementVector)[i])->GetN();
>> 260 }
>> 261 else
>> 262 {
>> 263 G4double * running = new G4double[theIso.GetNumberOfIsotopes(localZ)];
>> 264 for (i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
>> 265 {
>> 266 G4double fracInPercent=theIso.GetAbundance(theIso.GetFirstIsotope(localZ)+i);
>> 267 G4double runningA=theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
>> 268 running[i]=fracInPercent*pow(runningA, 2./3.);
>> 269 if(i!=0) running[i] += running[i-1];
>> 270 }
>> 271 G4double trial = G4UniformRand();
>> 272 for(i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
>> 273 {
>> 274 currentN = theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
>> 275 if(running[i]/running[theIso.GetNumberOfIsotopes(localZ)-1]>trial) break;
>> 276 }
>> 277 delete [] running;
250 } 278 }
>> 279 targetNucleus.SetParameters(currentN, currentZ);
>> 280 return (*theElementVector)[i];
251 } 281 }
252 } 282 }
253 } else { << 283 currentZ = G4double((*theElementVector)[numberOfElements-1]->GetZ());
254 // initialisation in worker threads << 284 G4int localZ = G4lrint(currentZ);
255 fXSType = masterProcess->CrossSectionType( << 285 if(noIsotopeWiseCrossSections)
256 fXSpeaks = masterProcess->TwoPeaksXS(); << 286 {
257 theEnergyOfCrossSectionMax = masterProcess << 287 currentN = (*theElementVector)[numberOfElements-1]->GetN();
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 } 288 }
369 // No warning for fStopButAlive which is a << 289 else
370 return theTotalResult; << 290 {
371 } << 291 G4double * running = new G4double[theIso.GetNumberOfIsotopes(localZ)];
372 << 292 for (i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
373 // Initialize the hadronic projectile from t << 293 {
374 thePro.Initialise(aTrack); << 294 G4double fracInPercent=theIso.GetAbundance(theIso.GetFirstIsotope(localZ)+i);
375 << 295 G4double runningA=theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
376 theInteraction = ChooseHadronicInteraction(t << 296 running[i]=fracInPercent*pow(runningA, 2./3.);
377 a << 297 // rough approximation; to get it better, redesign getMSC to not use G4Element
378 if(nullptr == theInteraction) { << 298 if(i!=0) running[i] += running[i-1];
379 G4ExceptionDescription ed; << 299 }
380 ed << "Target element "<<anElement->GetNam << 300 G4double trial = G4UniformRand();
381 << targetNucleus.GetZ_asInt() << " A= << 301 for(i=0; i<theIso.GetNumberOfIsotopes(localZ); i++)
382 << targetNucleus.GetA_asInt() << G4endl << 302 {
383 DumpState(aTrack,"ChooseHadronicInteractio << 303 currentN = theIso.GetIsotopeNucleonCount(theIso.GetFirstIsotope(localZ)+i);
384 ed << " No HadronicInteraction found out" << 304 if(running[i]/running[theIso.GetNumberOfIsotopes(localZ)-1]>trial) break;
385 G4Exception("G4HadronicProcess::PostStepDo << 305 }
386 FatalException, ed); << 306 delete [] running;
387 return theTotalResult; << 307 }
>> 308 targetNucleus.SetParameters(currentN, currentZ);
>> 309 return (*theElementVector)[numberOfElements-1];
388 } 310 }
389 << 311
390 G4HadFinalState* result = nullptr; << 312 G4VParticleChange *G4HadronicProcess::GeneralPostStepDoIt(
391 G4int reentryCount = 0; << 313 const G4Track &aTrack, const G4Step &)
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 { 314 {
>> 315 // G4cout << theNumberOfInteractionLengthLeft<<G4endl;
>> 316 #ifndef G4HadSignalHandler_off
>> 317 G4HadSignalHandler aHandler(G4HadronicProcess_local::G4HadronicProcessHandler_1);
>> 318 #endif
>> 319 if(aTrack.GetTrackStatus() != fAlive)
>> 320 {
>> 321 G4cerr << "G4HadronicProcess: track in unusable state - "
>> 322 <<aTrack.GetTrackStatus()<<G4endl;
>> 323 G4cerr << "G4HadronicProcess: returning unchanged track "<<G4endl;
>> 324 G4Exception("G4HadronicProcess", "001", JustWarning, "bailing out");
>> 325 theTotalResult->Clear();
>> 326 theTotalResult->Initialize(aTrack);
>> 327 return theTotalResult;
>> 328 }
>> 329 const G4DynamicParticle *aParticle = aTrack.GetDynamicParticle();
>> 330 G4Material *aMaterial = aTrack.GetMaterial();
>> 331 G4double originalEnergy = aParticle->GetKineticEnergy();
>> 332 G4double kineticEnergy = originalEnergy;
>> 333 if(aParticle->GetDefinition()->GetBaryonNumber()>1.5)
>> 334 {
>> 335 kineticEnergy/=aParticle->GetDefinition()->GetBaryonNumber();
>> 336 }
>> 337 G4Element * anElement = 0;
402 try 338 try
403 { 339 {
404 // Call the interaction << 340 anElement = ChooseAandZ( aParticle, aMaterial );
405 result = theInteraction->ApplyYourself( <<
406 ++reentryCount; <<
407 } 341 }
408 catch(G4HadronicException & aR) 342 catch(G4HadronicException & aR)
409 { 343 {
410 G4ExceptionDescription ed; << 344 aR.Report(G4cout);
411 aR.Report(ed); << 345 G4cout << "Unrecoverable error for:"<<G4endl;
412 ed << "Call for " << theInteraction->Get << 346 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
413 ed << "Target element "<<anElement->GetN << 347 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
414 << targetNucleus.GetZ_asInt() << 348 G4cout << " - Particle type = "
415 << " A= " << targetNucleus.GetA_asInt() << << 349 <<aParticle->GetDefinition()->GetParticleName()<<G4endl;
416 DumpState(aTrack,"ApplyYourself",ed); << 350 G4Exception("G4HadronicProcess", "007", FatalException,
417 ed << " ApplyYourself failed" << G4endl; << 351 "GeneralPostStepDoIt failed on element selection.");
418 G4Exception("G4HadronicProcess::PostStep <<
419 ed); <<
420 } <<
421 <<
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 } 352 }
436 } << 353 try
437 while(!result); /* Loop checking, 30-Oct-20 << 354 {
438 << 355 theInteraction = ChooseHadronicInteraction( kineticEnergy,
439 // Check whether kaon0 or anti_kaon0 are pre << 356 aMaterial, anElement );
440 // if this is the case, transform them into << 357 }
441 // with equal, 50% probability, keeping thei << 358 catch(G4HadronicException & aE)
442 // the other kinematical properties). << 359 {
443 // When this happens - very rarely - a "Just << 360 aE.Report(std::cout);
444 // Because Fluka-Cern produces kaon0 and ant << 361 G4cout << "Unrecoverable error for:"<<G4endl;
445 // of warnings to max 1 per thread. << 362 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
446 G4int nSec = (G4int)result->GetNumberOfSecon << 363 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
447 if ( nSec > 0 ) { << 364 G4cout << " - Particle type = "<<aParticle->GetDefinition()->GetParticleName()<<G4endl;
448 for ( G4int i = 0; i < nSec; ++i ) { << 365 G4Exception("G4HadronicProcess", "007", FatalException,
449 auto dynamicParticle = result->GetSecond << 366 "ChooseHadronicInteraction failed.");
450 auto part = dynamicParticle->GetParticle << 367 }
451 if ( part == G4KaonZero::Definition() || << 368 G4HadProjectile thePro(aTrack);
452 part == G4AntiKaonZero::Definition( << 369
453 G4ParticleDefinition* newPart; << 370 G4HadFinalState *result = 0;
454 if ( G4UniformRand() > 0.5 ) { newPart << 371 G4int reentryCount = 0;
455 else { newPart = G4KaonZeroLong::Defin << 372 do
456 dynamicParticle->SetDefinition( newPar << 373 {
457 if ( nKaonWarn < 1 ) { << 374 try
458 ++nKaonWarn; << 375 {
459 G4ExceptionDescription ed; << 376 G4HadronicInteractionWrapper aW;
460 ed << " Hadronic model " << theInteraction << 377 result = aW.ApplyInteraction(thePro, targetNucleus, theInteraction);
461 ed << " created " << part->GetParticleName << 378 }
462 ed << " -> forced to be " << newPart->GetP << 379 catch(G4HadReentrentException aR)
463 G4Exception( "G4HadronicProcess::PostStepD << 380 {
>> 381 aR.Report(G4cout);
>> 382 G4cout << " G4HadronicProcess re-entering the ApplyYourself call for"<<G4endl;
>> 383 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
>> 384 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
>> 385 G4cout << " - Particle type = "<<aParticle->GetDefinition()->GetParticleName()<<G4endl;
>> 386 result = 0; // here would still be leaking...
>> 387 if(reentryCount>100)
>> 388 {
>> 389 G4Exception("G4HadronicProcess", "007", FatalException,
>> 390 "GetHadronicProcess: Reentering ApplyYourself too often - GeneralPostStepDoIt failed.");
464 } 391 }
>> 392 G4Exception("G4HadronicProcess", "007", FatalException,
>> 393 "GetHadronicProcess: GeneralPostStepDoIt failed (Reentering ApplyYourself not yet supported.)");
>> 394 }
>> 395 catch(G4HadronicException aR)
>> 396 {
>> 397 aR.Report(G4cout);
>> 398 G4cout << " G4HadronicProcess failed in ApplyYourself call for"<<G4endl;
>> 399 G4cout << " - Particle energy[GeV] = "<< originalEnergy/GeV<<G4endl;
>> 400 G4cout << " - Material = "<<aMaterial->GetName()<<G4endl;
>> 401 G4cout << " - Particle type = "<<aParticle->GetDefinition()->GetParticleName()<<G4endl;
>> 402 G4Exception("G4HadronicProcess", "007", FatalException,
>> 403 "GeneralPostStepDoIt failed.");
465 } 404 }
466 } 405 }
467 } << 406 while(!result);
468 << 407 if(!ModelingState && !getenv("BypassAllSafetyChecks") )
469 result->SetTrafoToLab(thePro.GetTrafoToLab() << 408 {
470 FillResult(result, aTrack); << 409 G4cout << "ERROR IN EXECUTION -- HADRONIC PROCESS STATE NOT VALID"<<G4endl;
471 << 410 G4cout << "Result will be of undefined quality."<<G4endl;
472 if (epReportLevel != 0) { << 411 }
473 CheckEnergyMomentumConservation(aTrack, ta << 412 if(result->GetStatusChange() == isAlive && thePro.GetDefinition() != aTrack.GetDefinition())
474 } << 413 {
475 //G4cout << "PostStepDoIt done nICelectrons= << 414 G4DynamicParticle * aP = const_cast<G4DynamicParticle *>(aTrack.GetDynamicParticle());
476 return theTotalResult; << 415 aP->SetDefinition(const_cast<G4ParticleDefinition *>(thePro.GetDefinition()));
477 } << 416 }
478 << 417 result->SetTrafoToLab(thePro.GetTrafoToLab());
479 void G4HadronicProcess::ProcessDescription(std << 418 for(G4int i=0; i<result->GetNumberOfSecondaries(); i++)
480 { << 419 {
481 outFile << "The description for this process << 420 G4DynamicParticle* aSecTrack = result->GetSecondary(i)->GetParticle();
482 } << 421 if(aSecTrack->GetDefinition()->GetPDGCharge()>1.5)
483 << 422 {
484 G4double G4HadronicProcess::XBiasSurvivalProba << 423 G4EffectiveCharge aCalculator;
485 { << 424 G4double charge = aCalculator.GetCharge(aMaterial, aSecTrack->GetKineticEnergy(),
486 G4double nLTraversed = GetTotalNumberOfInter << 425 aSecTrack->GetDefinition()->GetPDGMass(),
487 G4double biasedProbability = 1.-G4Exp(-nLTra << 426 aSecTrack->GetDefinition()->GetPDGCharge());
488 G4double realProbability = 1-G4Exp(-nLTraver << 427 if(getenv("GHADChargeDebug"))
489 G4double result = (biasedProbability-realPro << 428 {
490 return result; << 429 std::cout << "Recoil fractional charge is "
491 } << 430 << charge/aSecTrack->GetDefinition()->GetPDGCharge()<<" "
492 << 431 << charge <<" "<<aSecTrack->GetDefinition()->GetPDGCharge()<<std::endl;
493 G4double G4HadronicProcess::XBiasSecondaryWeig << 432 }
494 { << 433 aSecTrack->SetCharge(charge);
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 } 434 }
569 dynParticle->SetKineticEnergy(e); <<
570 dynParticle->SetMass(mass); <<
571 } 435 }
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 436
578 G4Track* track = new G4Track(dynParticle, << 437 if(getenv("HadronicDoitLogging") )
579 track->SetCreatorModelID(idModel); << 438 {
580 track->SetParentResonanceDef(aR->GetSecond << 439 G4cout << "HadronicDoitLogging "
581 track->SetParentResonanceID(aR->GetSeconda << 440 << GetProcessName() <<" "
582 G4double newWeight = fWeight*aR->GetSecond << 441 << aParticle->GetDefinition()->GetPDGEncoding()<<" "
583 track->SetWeight(newWeight); << 442 << originalEnergy<<" "
584 track->SetTouchableHandle(aT.GetTouchableH << 443 << aParticle->GetMomentum()<<" "
585 theTotalResult->AddSecondary(track); << 444 << targetNucleus.GetN()<<" "
586 } << 445 << targetNucleus.GetZ()<<" "
587 aR->Clear(); << 446 << G4endl;
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 } 447 }
694 } << 448 ClearNumberOfInteractionLengthLeft();
695 return result; << 449 if(isoIsOnAnyway!=-1)
696 } << 450 {
697 << 451 if(isoIsEnabled||isoIsOnAnyway)
698 void << 452 {
699 G4HadronicProcess::CheckEnergyMomentumConserva << 453 result = DoIsotopeCounting(result, aTrack, targetNucleus);
700 << 454 }
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 } 455 }
744 } << 456
745 if( nSec > 0 ) { << 457 G4double e=aTrack.GetKineticEnergy();
746 G4Track* sec; << 458 ModelingState = 0;
747 << 459 if(e<5*GeV)
748 for (G4int i = 0; i < nSec; i++) { << 460 {
749 sec = theTotalResult->GetSecondary(i); << 461 for(size_t i=0; i<theBias.size(); i++)
750 final4mom += sec->GetDynamicParticle()-> << 462 {
751 final_A += sec->GetDefinition()->GetBary << 463 result = theBias[i]->Bias(result);
752 final_Z += G4lrint(sec->GetDefinition()- << 464 }
753 } 465 }
>> 466 FillTotalResult(result, aTrack);
>> 467 return theTotalResult;
754 } 468 }
755 469
756 // Get level-checking information (used to c << 470 G4HadFinalState * G4HadronicProcess::
757 G4String processName = GetProcessName(); << 471 DoIsotopeCounting(G4HadFinalState * aResult,
758 G4HadronicInteraction* theModel = GetHadroni << 472 const G4Track & aTrack,
759 G4String modelName("none"); << 473 const G4Nucleus & aNucleus)
760 if (theModel) modelName = theModel->GetModel << 474 {
761 std::pair<G4double, G4double> checkLevels = << 475 // get the PC from iso-production
762 if (!levelsSetByProcess) { << 476 if(theOldIsoResult) delete theOldIsoResult;
763 if (theModel) checkLevels = theModel->GetE << 477 if(theIsoResult) delete theIsoResult;
764 checkLevels.first= std::min(checkLevels.fi << 478 theIsoResult = new G4IsoParticleChange;
765 checkLevels.second=std::min(checkLevels.se << 479 G4bool done = false;
766 } << 480 G4IsoResult * anIsoResult = NULL;
767 << 481 for(unsigned int i=0; i<theProductionModels.size(); i++)
768 // Compute absolute total-energy difference, << 482 {
769 G4bool checkRelative = (aTrack.GetKineticEne << 483 anIsoResult = theProductionModels[i]->GetIsotope(aTrack, aNucleus);
770 << 484 if(anIsoResult!=NULL)
771 G4LorentzVector diff = initial4mom - final4m << 485 {
772 G4double absolute = diff.e(); << 486 done = true;
773 G4double relative = checkRelative ? absolute << 487 break;
774 << 488 }
775 G4double absolute_mom = diff.vect().mag(); << 489 }
776 G4double relative_mom = checkRelative ? abso << 490 // if none in charge, use default iso production
777 << 491 if(!done) anIsoResult = ExtractResidualNucleus(aTrack, aNucleus, aResult);
778 // Evaluate relative and absolute conservati << 492
779 G4bool relPass = true; << 493 // Add all info explicitely and add typename from model called.
780 G4String relResult = "pass"; << 494 theIsoResult->SetIsotope(anIsoResult->GetIsotope());
781 if ( std::abs(relative) > checkLevels.first << 495 theIsoResult->SetProductionPosition(aTrack.GetPosition());
782 || std::abs(relative_mom) > checkLevels.fir << 496 theIsoResult->SetProductionTime(aTrack.GetGlobalTime());
783 relPass = false; << 497 theIsoResult->SetParentParticle(*aTrack.GetDynamicParticle());
784 relResult = checkRelative ? "fail" : "N/A" << 498 theIsoResult->SetMotherNucleus(anIsoResult->GetMotherNucleus());
785 } << 499 // theIsoResult->SetProducer(typeid(*theInteraction).name()); @@@@@@@
786 << 500 G4String aWorkaround("WaitingForTypeidToBeAvailableInCompilers"); // @@@@@ workaround for DEC.
787 G4bool absPass = true; << 501 theIsoResult->SetProducer(aWorkaround);
788 G4String absResult = "pass"; << 502
789 if ( std::abs(absolute) > checkLevels.seco << 503 delete anIsoResult;
790 || std::abs(absolute_mom) > checkLevels. << 504
791 absPass = false ; << 505 return aResult;
792 absResult = "fail"; << 506 }
793 } << 507
794 << 508 G4IsoResult * G4HadronicProcess::
795 G4bool chargePass = true; << 509 ExtractResidualNucleus(const G4Track & ,
796 G4String chargeResult = "pass"; << 510 const G4Nucleus & aNucleus,
797 if ( (initial_A-final_A)!=0 << 511 G4HadFinalState * aResult)
798 || (initial_Z-final_Z)!=0 ) { << 512 {
799 chargePass = checkLevels.second < DBL_MAX << 513 G4double A = aNucleus.GetN();
800 chargeResult = "fail"; << 514 G4double Z = aNucleus.GetZ();
801 } << 515 G4double bufferA = 0;
802 << 516 G4double bufferZ = 0;
803 G4bool conservationPass = (relPass || absPas << 517
804 << 518 // loop over aResult, and decrement A, Z accordingly
805 std::stringstream Myout; << 519 // cash the max
806 G4bool Myout_notempty(false); << 520 for(G4int i=0; i<aResult->GetNumberOfSecondaries(); i++)
807 // Options for level of reporting detail: << 521 {
808 // 0. off << 522 G4HadSecondary* aSecTrack = aResult->GetSecondary(i);
809 // 1. report only when E/p not conserved << 523 if(bufferA<aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber())
810 // 2. report regardless of E/p conservation << 524 {
811 // 3. report only when E/p not conserved, w << 525 bufferA = aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
812 // 4. report regardless of E/p conservation << 526 bufferZ = aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
813 // negative -1.., as above, but send output << 527 }
814 << 528 Z-=aSecTrack->GetParticle()->GetDefinition()->GetPDGCharge();
815 if( std::abs(epReportLevel) == 4 << 529 A-=aSecTrack->GetParticle()->GetDefinition()->GetBaryonNumber();
816 || ( std::abs(epReportLevel) == 3 && ! cons << 530 }
817 Myout << " Process: " << processName << << 531
818 Myout << " Primary: " << aTrack.GetParti << 532 // if the fragment was part of the final state, it is
819 << " (" << aTrack.GetParticleDefin << 533 // assumed to be the heaviest secondary.
820 << " E= " << aTrack.GetDynamicPar << 534 if(A<0.1)
821 << ", target nucleus (" << aNucleus.GetZ << 535 {
822 << aNucleus.GetA_asInt() << ")" << G4end << 536 A = bufferA;
823 Myout_notempty=true; << 537 Z = bufferZ;
824 } << 538 }
825 if ( std::abs(epReportLevel) == 4 << 539
826 || std::abs(epReportLevel) == 2 << 540 // prepare the IsoResult.
827 || ! conservationPass ){ << 541 char the1[100] = {""};
828 << 542 std::ostrstream ost1(the1, 100, std::ios::out);
829 Myout << " "<< relResult <<" relative << 543 ost1 <<Z<<"_"<<A<<"\0";
830 << relative << " p/p(0)= " << rel << 544 G4String * biff = new G4String(the1);
831 Myout << " "<< absResult << " absolute << 545 G4IsoResult * theResult = new G4IsoResult(*biff, aNucleus);
832 << absolute/MeV << " / " << absol << 546
833 Myout << " "<< chargeResult << " charg << 547 // cleaning up.
834 Myout_notempty=true; << 548 delete biff;
835 << 549
836 } << 550 return theResult;
837 Myout.flush(); <<
838 if ( Myout_notempty ) { <<
839 if (epReportLevel > 0) G4cout << Myo <<
840 else if (epReportLevel < 0) G4cerr << Myo <<
841 } 551 }
842 } <<
843 552
844 void G4HadronicProcess::DumpState(const G4Trac << 553 G4double G4HadronicProcess::
845 const G4String& method, << 554 XBiasSurvivalProbability()
846 G4ExceptionDescription& ed) <<
847 { 555 {
848 ed << "Unrecoverable error in the method " < << 556 G4double result = 0;
849 << GetProcessName() << G4endl; << 557 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
850 ed << "TrackID= "<< aTrack.GetTrackID() << " << 558 G4double biasedProbability = 1.-exp(-nLTraversed);
851 << aTrack.GetParentID() << 559 G4double realProbability = 1-exp(-nLTraversed/aScaleFactor);
852 << " " << aTrack.GetParticleDefinition() << 560 result = (biasedProbability-realProbability)/biasedProbability;
853 << G4endl; << 561 return result;
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 } 562 }
895 563
896 G4double << 564 G4double G4HadronicProcess::
897 G4HadronicProcess::ComputeCrossSection(const G << 565 XBiasSecondaryWeight()
898 const G4Material* mat, <<
899 const G4double kinEnergy) <<
900 { 566 {
901 auto dp = new G4DynamicParticle(part, unitVe << 567 G4double result = 0;
902 G4double xs = theCrossSectionDataStore->Comp << 568 G4double nLTraversed = GetTotalNumberOfInteractionLengthTraversed();
903 delete dp; << 569 result = 1./aScaleFactor*exp(-nLTraversed/aScaleFactor*(1-1./aScaleFactor));
904 return xs; << 570 return result;
905 } 571 }
906 572
907 void G4HadronicProcess::RecomputeXSandMFP(cons << 573 struct G4Nancheck{ bool operator()(G4double aV){return (!(aV<1))&&(!(aV>-1));}};
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 574
916 void G4HadronicProcess::UpdateCrossSectionAndM << 575 void G4HadronicProcess::FillTotalResult(G4HadFinalState * aR, const G4Track & aT)
917 { 576 {
918 if(fXSType == fHadNoIntegral) { << 577 // G4cout << "############# Entry debug "
919 DefineXSandMFP(); << 578 // <<GetProcessName()<<" "
920 << 579 // <<aT.GetDynamicParticle()->GetDefinition()->GetParticleName()<<" "
921 } else if(fXSType == fHadIncreasing) { << 580 // <<aT.GetDynamicParticle()<<" "
922 if(e*invLambdaFactor < mfpKinEnergy) { << 581 // <<aScaleFactor<<" "
923 mfpKinEnergy = e; << 582 // <<aT.GetWeight()<<" "
924 ComputeXSandMFP(); << 583 // <<G4endl;
925 } << 584 G4Nancheck go_wild;
926 << 585 theTotalResult->Clear();
927 } else if(fXSType == fHadDecreasing) { << 586 theTotalResult->SetLocalEnergyDeposit(0.);
928 if(e < mfpKinEnergy && mfpKinEnergy > minK << 587 theTotalResult->Initialize(aT);
929 G4double e1 = std::max(e*lambdaFactor, m << 588 theTotalResult->SetSecondaryWeightByProcess(true);
930 mfpKinEnergy = e1; << 589 theTotalResult->SetStatusChange(fAlive);
931 RecomputeXSandMFP(e1); << 590 G4double rotation = 2.*pi*G4UniformRand();
932 } << 591 G4ThreeVector it(0., 0., 1.);
933 << 592 /*
934 } else if(fXSType == fHadOnePeak) { << 593 if(xBiasOn)
935 G4double epeak = (*theEnergyOfCrossSection << 594 {
936 if(e <= epeak) { << 595 G4cout << "BiasDebug "<<GetProcessName()<<" "
937 if(e*invLambdaFactor < mfpKinEnergy) { << 596 <<aScaleFactor<<" "
938 mfpKinEnergy = e; << 597 <<XBiasSurvivalProbability()<<" "
939 ComputeXSandMFP(); << 598 <<XBiasSecondaryWeight()<<" "
940 } << 599 <<G4endl;
941 } else if(e < mfpKinEnergy) { << 600 }
942 G4double e1 = std::max(epeak, e*lambdaFa << 601 */
943 mfpKinEnergy = e1; << 602 // if(GetProcessName() != "LElastic") std::cout << "Debug -1 "<<aR->GetStatusChange()<<std::endl;
944 RecomputeXSandMFP(e1); << 603 if(aR->GetStatusChange()==stopAndKill)
945 } << 604 {
946 << 605 if( xBiasOn && G4UniformRand()<XBiasSurvivalProbability() )
947 } else if(fXSType == fHadTwoPeaks) { << 606 {
948 G4TwoPeaksHadXS* xs = (*fXSpeaks)[matIdx]; << 607 theTotalResult->SetWeightChange( XBiasSurvivalProbability()*aT.GetWeight() );
949 const G4double e1peak = xs->e1peak; << 608 }
950 << 609 else
951 // below the 1st peak << 610 {
952 if(e <= e1peak) { << 611 theTotalResult->SetStatusChange(fStopAndKill);
953 if(e*invLambdaFactor < mfpKinEnergy) { << 612 theTotalResult->SetEnergyChange( 0.0 );
954 mfpKinEnergy = e; << 613 }
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 } 614 }
976 return; << 615 else if(aR->GetStatusChange()!=stopAndKill )
977 } << 616 {
978 const G4double e2deep = xs->e2deep; << 617 if(aR->GetStatusChange()==suspend)
979 // above the 2nd peak, below the deep << 618 {
980 if(e <= e2deep) { << 619 theTotalResult->SetStatusChange(fSuspend);
981 if(mfpKinEnergy >= e2deep || e <= mfpKin << 620 if(xBiasOn)
982 const G4double e1 = std::max(e2peak, e << 621 {
983 mfpKinEnergy = e1; << 622 G4Exception("G4HadronicProcess", "007", FatalException,
984 RecomputeXSandMFP(e1); << 623 "Cannot cross-section bias a process that suspends tracks.");
>> 624 }
>> 625 }
>> 626 if(xBiasOn && G4UniformRand()<XBiasSurvivalProbability())
>> 627 {
>> 628 theTotalResult->SetWeightChange( XBiasSurvivalProbability()*aT.GetWeight() );
>> 629 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
>> 630 if(go_wild(aR->GetEnergyChange()))
>> 631 {
>> 632 G4Exception("G4HadronicProcess", "007", FatalException,
>> 633 "surviving track received NaN energy.");
>> 634 }
>> 635 if(go_wild(aR->GetMomentumChange().x()) ||
>> 636 go_wild(aR->GetMomentumChange().y()) ||
>> 637 go_wild(aR->GetMomentumChange().z()))
>> 638 {
>> 639 G4Exception("G4HadronicProcess", "007", FatalException,
>> 640 "surviving track received NaN momentum.");
>> 641 }
>> 642 G4double newM=aT.GetDefinition()->GetPDGMass();
>> 643 G4double newE=aR->GetEnergyChange() + newM;
>> 644 G4double newP=sqrt(newE*newE - newM*newM);
>> 645 G4DynamicParticle * aNew =
>> 646 new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
>> 647 G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
>> 648 aR->AddSecondary(theSec);
>> 649 }
>> 650 else
>> 651 {
>> 652 G4double newWeight = aR->GetWeightChange()*aT.GetWeight();
>> 653 theTotalResult->SetWeightChange(newWeight); // This is multiplicative
>> 654 if(aR->GetEnergyChange()>-.5)
>> 655 {
>> 656 if(go_wild(aR->GetEnergyChange()))
>> 657 {
>> 658 G4Exception("G4HadronicProcess", "007", FatalException,
>> 659 "track received NaN energy.");
>> 660 }
>> 661 theTotalResult->SetEnergyChange(aR->GetEnergyChange());
>> 662 }
>> 663 G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.);
>> 664 newDirection*=aR->GetTrafoToLab();
>> 665 theTotalResult->SetMomentumDirectionChange(newDirection.vect());
>> 666 }
985 } 667 }
986 return; << 668 else
987 } << 669 {
988 const G4double e3peak = xs->e3peak; << 670 G4cerr << "Track status is "<< aR->GetStatusChange()<<G4endl;
989 // above the deep, below 3d peak << 671 G4Exception("G4HadronicProcess", "007", FatalException,
990 if(e <= e3peak) { << 672 "use of unsupported track-status.");
991 if(e*invLambdaFactor < mfpKinEnergy) { << 673 }
992 mfpKinEnergy = e; << 674 if(GetProcessName() != "LElastic"
993 ComputeXSandMFP(); << 675 &&
>> 676 AlwaysKillLeadingHadron()
>> 677 &&
>> 678 theTotalResult->GetStatusChange()==fAlive
>> 679 &&
>> 680 aR->GetStatusChange()==isAlive)
>> 681 {
>> 682 G4double newWeight = theTotalResult->GetWeightChange();
>> 683 G4double newM=aT.GetDefinition()->GetPDGMass();
>> 684 G4double newE=aR->GetEnergyChange() + newM;
>> 685 G4double newP=sqrt(newE*newE - newM*newM);
>> 686 G4DynamicParticle * aNew =
>> 687 new G4DynamicParticle(aT.GetDefinition(), newE, newP*aR->GetMomentumChange());
>> 688 // std::cout << "Debug 0 "<<aR->GetNumberOfSecondaries()<<std::endl;
>> 689 //std::cout << "Debug 1 "<<aR->GetEnergyChange()<<" "<< aNew->GetTotalEnergy() <<std::endl;
>> 690 //std::cout << "Debug 2 "<<aR->GetMomentumChange()<<" "<< aNew->GetMomentum() << std::endl;
>> 691 //std::cout << "Debug 3 "<<newWeight<<std::endl;
>> 692 //std::cout << std::endl;
>> 693 G4HadSecondary * theSec = new G4HadSecondary(aNew, newWeight);
>> 694 aR->AddSecondary(theSec);
>> 695 aR->SetStatusChange(stopAndKill);
>> 696 theTotalResult->SetStatusChange(fStopAndKill);
>> 697 theTotalResult->SetEnergyChange( 0.0 );
>> 698 //std::cout << "Debug 4 "<< aR->GetNumberOfSecondaries() <<std::endl;
>> 699 }
>> 700 theTotalResult->SetLocalEnergyDeposit(aR->GetLocalEnergyDeposit());
>> 701 theTotalResult->SetNumberOfSecondaries(aR->GetNumberOfSecondaries());
>> 702
>> 703 for(G4int i=0; i<aR->GetNumberOfSecondaries(); i++)
>> 704 {
>> 705 //std::cout << "Debug 5 "<< aR->GetNumberOfSecondaries() <<std::endl;
>> 706 G4LorentzVector theM = aR->GetSecondary(i)->GetParticle()->Get4Momentum();
>> 707 theM.rotate(rotation, it);
>> 708 theM*=aR->GetTrafoToLab();
>> 709 if(go_wild(theM.e()))
>> 710 {
>> 711 G4Exception("G4HadronicProcess", "007", FatalException,
>> 712 "secondary track received NaN energy.");
>> 713 }
>> 714 if(go_wild(theM.x()) ||
>> 715 go_wild(theM.y()) ||
>> 716 go_wild(theM.z()))
>> 717 {
>> 718 G4Exception("G4HadronicProcess", "007", FatalException,
>> 719 "secondary track received NaN momentum.");
>> 720 }
>> 721 aR->GetSecondary(i)->GetParticle()->Set4Momentum(theM);
>> 722 G4double time = aR->GetSecondary(i)->GetTime();
>> 723 if(time<0) time = aT.GetGlobalTime();
>> 724 G4Track* track = new G4Track(aR->GetSecondary(i)->GetParticle(),
>> 725 aT.GetGlobalTime(),
>> 726 aT.GetPosition());
>> 727 G4double newWeight = aT.GetWeight()*aR->GetSecondary(i)->GetWeight();
>> 728 //static G4double pinelcount=0;
>> 729 if(xBiasOn) newWeight *= XBiasSecondaryWeight();
>> 730 // G4cout << "#### ParticleDebug "
>> 731 // <<GetProcessName()<<" "
>> 732 // <<aR->GetSecondary(i)->GetParticle()->GetDefinition()->GetParticleName()<<" "
>> 733 // <<aScaleFactor<<" "
>> 734 // <<XBiasSurvivalProbability()<<" "
>> 735 // <<XBiasSecondaryWeight()<<" "
>> 736 // <<aT.GetWeight()<<" "
>> 737 // <<aR->GetSecondary(i)->GetWeight()<<" "
>> 738 // <<aR->GetSecondary(i)->GetParticle()->Get4Momentum()<<" "
>> 739 // <<G4endl;
>> 740 track->SetWeight(newWeight);
>> 741 G4double trackDeb = track->GetKineticEnergy();
>> 742 if( ( trackDeb<0
>> 743 || (trackDeb>aT.GetKineticEnergy()+1*GeV) ) && getenv("GHADEnergyBalanceDebug") )
>> 744 {
>> 745 G4cout << "Debugging hadronic processes: "<<track->GetKineticEnergy()
>> 746 <<" "<<aT.GetKineticEnergy()
>> 747 <<" "<<GetProcessName()
>> 748 <<" "<<aT.GetDefinition()->GetParticleName()
>> 749 <<G4endl;
>> 750 }
>> 751 /*if(GetProcessName()=="PhotonInelastic")
>> 752 {
>> 753 if(aR->GetSecondary(i)->GetParticle()->GetDefinition()==G4Neutron::NeutronDefinition())
>> 754 {
>> 755 pinelcount+= newWeight;
>> 756 G4cout << "=======> Neutrons from gamma-nuclear "<<pinelcount<<G4endl;
>> 757 }
>> 758 }*/
>> 759
>> 760 theTotalResult->AddSecondary(track);
994 } 761 }
>> 762 aR->Clear();
995 return; 763 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 <<
1004 } else { <<
1005 DefineXSandMFP(); <<
1006 } <<
1007 } 764 }
>> 765 /* end of file */
1008 766