Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 29 // 30 // 12-Apr-06 fix in delayed neutron and photon 31 // 07-Sep-11 M. Kelsey -- Follow change to G4H 32 // P. Arce, June-2014 Conversion neutron_hp to 33 // 34 35 #include "G4ParticleHPFissionFS.hh" 36 37 #include "G4DynamicParticleVector.hh" 38 #include "G4Exp.hh" 39 #include "G4IonTable.hh" 40 #include "G4Nucleus.hh" 41 #include "G4ParticleHPFissionERelease.hh" 42 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicsModelCatalog.hh" 44 45 G4ParticleHPFissionFS::G4ParticleHPFissionFS() 46 { 47 secID = G4PhysicsModelCatalog::GetModelID("m 48 hasXsec = false; 49 produceFissionFragments = false; 50 } 51 52 void G4ParticleHPFissionFS::Init(G4double A, G 53 const G4Strin 54 { 55 theFS.Init(A, Z, M, dirName, aFSType, projec 56 theFC.Init(A, Z, M, dirName, aFSType, projec 57 theSC.Init(A, Z, M, dirName, aFSType, projec 58 theTC.Init(A, Z, M, dirName, aFSType, projec 59 theLC.Init(A, Z, M, dirName, aFSType, projec 60 61 theFF.Init(A, Z, M, dirName, aFSType, projec 62 if (G4ParticleHPManager::GetInstance()->GetP 63 G4cout << "Fission fragment production is 64 << "Z = " << (G4int)Z << ", A = " < 65 G4cout << "As currently modeled this optio 66 "fission fragments." 67 << G4endl; 68 produceFissionFragments = true; 69 } 70 } 71 72 G4HadFinalState* G4ParticleHPFissionFS::ApplyY 73 { 74 // Because it may change by UI command 75 produceFissionFragments = G4ParticleHPManage 76 77 // prepare neutron 78 if (theResult.Get() == nullptr) theResult.Pu 79 theResult.Get()->Clear(); 80 G4double eKinetic = theTrack.GetKineticEnerg 81 const G4HadProjectile* incidentParticle = &t 82 G4ReactionProduct theNeutron( 83 const_cast<G4ParticleDefinition*>(incident 84 theNeutron.SetMomentum(incidentParticle->Get 85 theNeutron.SetKineticEnergy(eKinetic); 86 87 // prepare target 88 G4Nucleus aNucleus; 89 G4ReactionProduct theTarget; 90 G4double targetMass = theFS.GetMass(); 91 G4ThreeVector neuVelo = 92 (1. / incidentParticle->GetDefinition()->G 93 theTarget = 94 aNucleus.GetBiasedThermalNucleus(targetMas 95 theTarget.SetDefinition( 96 G4IonTable::GetIonTable()->GetIon(G4int(th 97 98 // set neutron and target in the FS classes 99 theFS.SetNeutronRP(theNeutron); 100 theFS.SetTarget(theTarget); 101 theFC.SetNeutronRP(theNeutron); 102 theFC.SetTarget(theTarget); 103 theSC.SetNeutronRP(theNeutron); 104 theSC.SetTarget(theTarget); 105 theTC.SetNeutronRP(theNeutron); 106 theTC.SetTarget(theTarget); 107 theLC.SetNeutronRP(theNeutron); 108 theLC.SetTarget(theTarget); 109 theFF.SetNeutronRP(theNeutron); 110 theFF.SetTarget(theTarget); 111 112 // boost to target rest system and decide on 113 theNeutron.Lorentz(theNeutron, -1 * theTarge 114 115 // dice the photons 116 117 G4DynamicParticleVector* thePhotons; 118 thePhotons = theFS.GetPhotons(); 119 120 // select the FS in charge 121 122 eKinetic = theNeutron.GetKineticEnergy(); 123 G4double xSec[4]; 124 xSec[0] = theFC.GetXsec(eKinetic); 125 xSec[1] = xSec[0] + theSC.GetXsec(eKinetic); 126 xSec[2] = xSec[1] + theTC.GetXsec(eKinetic); 127 xSec[3] = xSec[2] + theLC.GetXsec(eKinetic); 128 G4int it; 129 unsigned int i = 0; 130 G4double random = G4UniformRand(); 131 if (xSec[3] == 0) { 132 it = -1; 133 } 134 else { 135 for (i = 0; i < 4; i++) { 136 it = i; 137 if (random < xSec[i] / xSec[3]) break; 138 } 139 } 140 141 // dice neutron multiplicities, energies and 142 // no energy conservation on an event-to-eve 143 // also for mean, we rely on the consistancy 144 145 G4int Prompt = 0, delayed = 0, all = 0; 146 G4DynamicParticleVector* theNeutrons = nullp 147 switch (it) // check logic, and ask, if par 148 // particles @@@ 149 { 150 case 0: 151 theFS.SampleNeutronMult(all, Prompt, del 152 if (Prompt == 0 && delayed == 0) Prompt 153 theNeutrons = theFC.ApplyYourself(Prompt 154 // take 'U' into account explicitly (see 155 break; 156 case 1: 157 theFS.SampleNeutronMult(all, Prompt, del 158 if (Prompt == 0 && delayed == 0) Prompt 159 theNeutrons = theSC.ApplyYourself(Prompt 160 break; 161 case 2: 162 theFS.SampleNeutronMult(all, Prompt, del 163 if (Prompt == 0 && delayed == 0) Prompt 164 theNeutrons = theTC.ApplyYourself(Prompt 165 break; 166 case 3: 167 theFS.SampleNeutronMult(all, Prompt, del 168 if (Prompt == 0 && delayed == 0) Prompt 169 theNeutrons = theLC.ApplyYourself(Prompt 170 break; 171 default: 172 break; 173 } 174 175 // dice delayed neutrons and photons, and fa 176 // for Prompt in case channel had no FS data 177 178 if (produceFissionFragments) delayed = 0; 179 180 G4double* theDecayConstants; 181 182 if (theNeutrons != nullptr) { 183 theDecayConstants = new G4double[delayed]; 184 for (i = 0; i < theNeutrons->size(); ++i) 185 theResult.Get()->AddSecondary(theNeutron 186 } 187 delete theNeutrons; 188 189 G4DynamicParticleVector* theDelayed = null 190 theDelayed = theFS.ApplyYourself(0, delaye 191 for (i = 0; i < theDelayed->size(); i++) { 192 G4double time = -G4Log(G4UniformRand()) 193 time += theTrack.GetGlobalTime(); 194 theResult.Get()->AddSecondary(theDelayed 195 theResult.Get()->GetSecondary(theResult. 196 } 197 delete theDelayed; 198 } 199 else { 200 theFS.SampleNeutronMult(all, Prompt, delay 201 theDecayConstants = new G4double[delayed]; 202 if (Prompt == 0 && delayed == 0) Prompt = 203 theNeutrons = theFS.ApplyYourself(Prompt, 204 G4int i0; 205 for (i0 = 0; i0 < Prompt; ++i0) { 206 theResult.Get()->AddSecondary(theNeutron 207 } 208 209 for (i0 = Prompt; i0 < Prompt + delayed; + 210 // Protect against the very rare case of 211 G4double time = 0.0; 212 if (theDecayConstants[i0 - Prompt] > 1.0 213 time = -G4Log(G4UniformRand()) / theDe 214 } 215 else { 216 G4ExceptionDescription ed; 217 ed << " theDecayConstants[i0-Prompt]=" 218 << " -> cannot sample the time : 219 G4Exception("G4ParticleHPFissionFS::Ap 220 } 221 222 time += theTrack.GetGlobalTime(); 223 theResult.Get()->AddSecondary(theNeutron 224 theResult.Get()->GetSecondary(theResult. 225 } 226 delete theNeutrons; 227 } 228 delete[] theDecayConstants; 229 230 std::size_t nPhotons = 0; 231 if (thePhotons != nullptr) { 232 nPhotons = thePhotons->size(); 233 for (i = 0; i < nPhotons; ++i) { 234 theResult.Get()->AddSecondary(thePhotons 235 } 236 delete thePhotons; 237 } 238 239 // finally deal with local energy deposition 240 241 G4ParticleHPFissionERelease* theERelease = t 242 G4double eDepByFragments = theERelease->GetF 243 // theResult.SetLocalEnergyDeposit(eDepByFra 244 if (!produceFissionFragments) theResult.Get( 245 // clean up the primary neutron 246 theResult.Get()->SetStatusChange(stopAndKill 247 248 if (produceFissionFragments) { 249 G4int fragA_Z = 0; 250 G4int fragA_A = 0; 251 G4int fragA_M = 0; 252 // System is traget rest! 253 theFF.GetAFissionFragment(eKinetic, fragA_ 254 if (0 == fragA_A) { return theResult.Get() 255 G4int fragB_Z = (G4int)theBaseZ - fragA_Z; 256 G4int fragB_A = (G4int)theBaseA - fragA_A 257 258 G4IonTable* pt = G4IonTable::GetIonTable() 259 // Excitation energy is not taken into acc 260 G4ParticleDefinition* pdA = pt->GetIon(fra 261 G4ParticleDefinition* pdB = pt->GetIon(fra 262 263 // Isotropic Distribution 264 G4double phi = twopi * G4UniformRand(); 265 // Bug #1745 DHW G4double theta = pi*G4Uni 266 G4double costheta = 2. * G4UniformRand() - 267 G4double theta = std::acos(costheta); 268 G4double sinth = std::sin(theta); 269 G4ThreeVector direction(sinth * std::cos(p 270 271 // Just use ENDF value for this 272 G4double ER = eDepByFragments; 273 G4double ma = pdA->GetPDGMass(); 274 G4double mb = pdB->GetPDGMass(); 275 G4double EA = ER / (1 + ma / mb); 276 G4double EB = ER - EA; 277 auto dpA = new G4DynamicParticle(pdA, dire 278 auto dpB = new G4DynamicParticle(pdB, -dir 279 theResult.Get()->AddSecondary(dpA, secID); 280 theResult.Get()->AddSecondary(dpB, secID); 281 } 282 283 return theResult.Get(); 284 } 285