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*std::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*std::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*std::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.-std::exp(-nLTraversed); 851 << aTrack.GetParentID() << 559 G4double realProbability = 1-std::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*std::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->ProposeLocalEnergyDeposit(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->ProposeTrackStatus(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->ProposeParentWeight( 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->ProposeTrackStatus(fStopAndKill); 953 if(e*invLambdaFactor < mfpKinEnergy) { << 612 theTotalResult->ProposeEnergy( 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->ProposeTrackStatus(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->ProposeParentWeight( 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=std::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->ProposeParentWeight(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->ProposeEnergy(aR->GetEnergyChange()); >> 662 } >> 663 G4LorentzVector newDirection(aR->GetMomentumChange().unit(), 1.); >> 664 newDirection*=aR->GetTrafoToLab(); >> 665 theTotalResult->ProposeMomentumDirection(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->GetTrackStatus()==fAlive >> 679 && >> 680 aR->GetStatusChange()==isAlive) >> 681 { >> 682 G4double newWeight = theTotalResult->GetParentWeight(); >> 683 G4double newM=aT.GetDefinition()->GetPDGMass(); >> 684 G4double newE=aR->GetEnergyChange() + newM; >> 685 G4double newP=std::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->ProposeTrackStatus(fStopAndKill); >> 697 theTotalResult->ProposeEnergy( 0.0 ); >> 698 //std::cout << "Debug 4 "<< aR->GetNumberOfSecondaries() <<std::endl; >> 699 } >> 700 theTotalResult->ProposeLocalEnergyDeposit(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 track->SetTouchableHandle(aT.GetTouchableHandle()); >> 761 theTotalResult->AddSecondary(track); 994 } 762 } >> 763 aR->Clear(); 995 return; 764 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 } 765 } >> 766 /* end of file */ 1008 767