Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPCaptureFS.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// neutron_hp -- source file
27// J.P. Wellisch, Nov-1996
28// A prototype of the low energy neutron transport model.
29//
30// 12-April-06 Enable IC electron emissions T. Koi
31// 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION flag
32// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33// 101203 Bugzilla/Geant4 Problem 1155 Lack of residual in some case
34// 110430 Temporary solution in the case of being MF6 final state in Capture reaction (MT102)
35//
36// P. Arce, June-2014 Conversion neutron_hp to particle_hp
37//
41#include "G4SystemOfUnits.hh"
42#include "G4Gamma.hh"
43#include "G4ReactionProduct.hh"
44#include "G4Nucleus.hh"
46#include "G4Fragment.hh"
47#include "G4IonTable.hh"
50
51
53 {
54 secID = G4PhysicsModelCatalog::GetModelID( "model_NeutronHPCapture" );
55 hasXsec = false;
56 hasExactMF6 = false;
57 targetMass = 0;
58 }
59
60
62 {
63
64 if ( theResult.Get() == nullptr ) theResult.Put( new G4HadFinalState );
65 theResult.Get()->Clear();
66
67 G4int i;
68
69// prepare neutron
70 G4double eKinetic = theTrack.GetKineticEnergy();
71 const G4HadProjectile *incidentParticle = &theTrack;
72 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition() ) );
73 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
74 theNeutron.SetKineticEnergy( eKinetic );
75
76 // Prepare target
77 G4ReactionProduct theTarget;
78 G4Nucleus aNucleus;
79 G4double eps = 0.0001;
80 if (targetMass < 500*MeV) targetMass =
81 (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps) )) /
83 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
84 G4double temperature = theTrack.GetMaterial()->GetTemperature();
85 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
87
88 // Put neutron in nucleus rest system
89 theNeutron.Lorentz(theNeutron, theTarget);
90 eKinetic = theNeutron.GetKineticEnergy();
91
92 // Sample the photons
93 G4ReactionProductVector * thePhotons = 0;
94 if ( HasFSData() && !G4ParticleHPManager::GetInstance()->GetUseOnlyPhotoEvaporation() )
95 {
96 //NDL has final state data
97 if ( hasExactMF6 ) {
98 theMF6FinalState.SetTarget(theTarget);
99 theMF6FinalState.SetProjectileRP(theNeutron);
100 thePhotons = theMF6FinalState.Sample( eKinetic );
101 } else {
102 thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
103 }
104 if ( thePhotons == NULL ) {
105 throw G4HadronicException(__FILE__, __LINE__, "Final state data for photon is not properly allocated");
106 }
107 }
108 else
109 {
110 //NDL does not have final state data or forced to use PhotoEvaporation model
111 G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum();
112 G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
113 G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4);
114 G4PhotonEvaporation photonEvaporation;
115 // T. K. add
116 photonEvaporation.SetICM( TRUE );
117 G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
118 thePhotons = new G4ReactionProductVector;
119 for(auto it=products->cbegin(); it!=products->cend(); ++it)
120 {
122 // T. K. add
123 if ( (*it)->GetParticleDefinition() != 0 )
124 theOne->SetDefinition( (*it)->GetParticleDefinition() );
125 else
126 theOne->SetDefinition( G4Gamma::Gamma() ); // this definiion will be over writen
127
128 // T. K. comment out below line
129 //theOne->SetDefinition( G4Gamma::Gamma() );
131 if ( (*it)->GetMomentum().mag() > 10*MeV)
132 theOne->SetDefinition(theTable->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0) );
133
134 if ( (*it)->GetExcitationEnergy() > 1.0e-2*eV) {
135 G4double ex = (*it)->GetExcitationEnergy();
137 aPhoton->SetDefinition( G4Gamma::Gamma() );
138 aPhoton->SetMomentum( (*it)->GetMomentum().vect().unit() * ex );
139 //aPhoton->SetTotalEnergy( ex ); //will be calculated from momentum
140 thePhotons->push_back(aPhoton);
141 }
142
143 theOne->SetMomentum( (*it)->GetMomentum().vect() * ( (*it)->GetMomentum().t() - (*it)->GetExcitationEnergy() ) / (*it)->GetMomentum().t() ) ;
144 thePhotons->push_back(theOne);
145 delete *it;
146 }
147 delete products;
148 }
149
150 // Add them to the final state
151 G4int nPhotons = 0;
152 nPhotons=(G4int)thePhotons->size();
153
154///*
156//Make at least one photon
157//101203 TK
158 if ( nPhotons == 0 )
159 {
161 theOne->SetDefinition( G4Gamma::Gamma() );
162 // Bug #1745 DHW G4double theta = pi*G4UniformRand();
163 G4double costheta = 2.*G4UniformRand()-1.;
164 G4double theta = std::acos(costheta);
165 G4double phi = twopi*G4UniformRand();
166 G4double sinth = std::sin(theta);
167 G4ThreeVector direction(sinth*std::cos(phi), sinth*std::sin(phi), costheta);
168 theOne->SetMomentum(direction);
169 thePhotons->push_back(theOne);
170 ++nPhotons; // 0 -> 1
171 }
172//One photon case: energy set to Q-value
173//101203 TK
174 //if ( nPhotons == 1 )
175 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
176 {
177 G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
178
180 - G4IonTable::GetIonTable()->GetIonMass(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
181
182 thePhotons->operator[](0)->SetMomentum( Q*direction );
183 }
184//
185 }
186
187 // back to lab system
188 for(i=0; i<nPhotons; i++)
189 {
190 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1*theTarget);
191 }
192
193 // Recoil, if only one gamma
194 //if (1==nPhotons)
195 if ( nPhotons == 1 && thePhotons->operator[](0)->GetDefinition()->GetBaryonNumber() == 0 )
196 {
199 ->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
200 theOne->SetDefinition(aRecoil);
201 // Now energy;
202 // Can be done slightly better @
203 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect()
204 +theTarget.GetMomentum()
205 -thePhotons->operator[](0)->GetMomentum();
206
207 //TKDB 140520
208 //G4ThreeVector theMomUnit = aMomentum.unit();
209 //G4double aKinEnergy = theTrack.GetKineticEnergy()
210 // +theTarget.GetKineticEnergy(); // gammas come from Q-value
211 //G4double theResMass = aRecoil->GetPDGMass();
212 //G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
213 //G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
214 //G4ThreeVector theMomentum = theAbsMom*theMomUnit;
215 //theOne->SetMomentum(theMomentum);
216
217 theOne->SetMomentum(aMomentum);
218 theResult.Get()->AddSecondary(theOne, secID);
219 }
220
221 // Now fill in the gammas.
222 for(i=0; i<nPhotons; i++)
223 {
224 // back to lab system
226 theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
227 theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
228 theResult.Get()->AddSecondary(theOne, secID);
229 delete thePhotons->operator[](i);
230 }
231 delete thePhotons;
232
233//101203TK
234 G4bool residual = false;
236 ->GetIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0);
237 for ( std::size_t j = 0 ; j != theResult.Get()->GetNumberOfSecondaries() ; j++ )
238 {
239 if ( theResult.Get()->GetSecondary(j)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
240 }
241
242 if ( residual == false )
243 {
244 G4int nNonZero = 0;
245 G4LorentzVector p_photons(0,0,0,0);
246 for ( std::size_t j = 0 ; j != theResult.Get()->GetNumberOfSecondaries() ; j++ )
247 {
248 p_photons += theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum();
249 // To many 0 momentum photons -> Check PhotonDist
250 if ( theResult.Get()->GetSecondary(j)->GetParticle()->Get4Momentum().e() > 0 ) nNonZero++;
251 }
252
253 // Can we include kinetic energy here?
254 G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
255 - ( p_photons.e() + aRecoil->GetPDGMass() );
256
257//Add photons
258 if ( nPhotons - nNonZero > 0 )
259 {
260 //G4cout << "TKDB G4ParticleHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl;
261 std::vector<G4double> vRand;
262 vRand.push_back( 0.0 );
263 for ( G4int j = 0 ; j != nPhotons - nNonZero - 1 ; j++ )
264 {
265 vRand.push_back( G4UniformRand() );
266 }
267 vRand.push_back( 1.0 );
268 std::sort( vRand.begin(), vRand.end() );
269
270 std::vector<G4double> vEPhoton;
271 for ( G4int j = 0 ; j < (G4int)vRand.size() - 1 ; j++ )
272 {
273 vEPhoton.push_back( deltaE * ( vRand[j+1] - vRand[j] ) );
274 }
275 std::sort( vEPhoton.begin(), vEPhoton.end() );
276
277 for ( G4int j = 0 ; j < nPhotons - nNonZero - 1 ; j++ )
278 {
279 //Isotopic in LAB OK?
280 // Bug # 1745 DHW G4double theta = pi*G4UniformRand();
281 G4double costheta = 2.*G4UniformRand()-1.;
282 G4double theta = std::acos(costheta);
283 G4double phi = twopi*G4UniformRand();
284 G4double sinth = std::sin(theta);
285 G4double en = vEPhoton[j];
286 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*costheta);
287
288 p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
290 theOne->SetDefinition( G4Gamma::Gamma() );
291 theOne->SetMomentum( tempVector );
292 theResult.Get()->AddSecondary(theOne, secID);
293 }
294
295// Add last photon
297 theOne->SetDefinition( G4Gamma::Gamma() );
298// For better momentum conservation
299 G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
300 p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
301 theOne->SetMomentum( lastPhoton );
302 theResult.Get()->AddSecondary(theOne, secID);
303 }
304
305//Add residual
307 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
308 - p_photons.vect();
309 theOne->SetDefinition(aRecoil);
310 theOne->SetMomentum( aMomentum );
311 theResult.Get()->AddSecondary(theOne, secID);
312
313 }
314//101203TK END
315
316// clean up the primary neutron
318 return theResult.Get();
319 }
320
321#include <sstream>
323 {
324
325 //TK110430 BEGIN
326 std::stringstream ss;
327 ss << static_cast<G4int>(Z);
328 G4String sZ;
329 ss >> sZ;
330 ss.clear();
331 ss << static_cast<G4int>(A);
332 G4String sA;
333 ss >> sA;
334
335 ss.clear();
336 G4String sM;
337 if ( M > 0 )
338 {
339 ss << "m";
340 ss << M;
341 ss >> sM;
342 ss.clear();
343 }
344
345 G4String element_name = theNames.GetName( static_cast<G4int>(Z)-1 );
346 G4String filenameMF6 = dirName+"/FSMF6/"+sZ+"_"+sA+sM+"_"+element_name;
347 //std::ifstream dummyIFS(filenameMF6, std::ios::in);
348 //if ( dummyIFS.good() == true ) hasExactMF6=true;
349 std::istringstream theData(std::ios::in);
350 G4ParticleHPManager::GetInstance()->GetDataStream(filenameMF6,theData);
351
352 //TK110430 Only use MF6MT102 which has exactly same A and Z
353 //Even _nat_ do not select and there is no _nat_ case in ENDF-VII.0
354 if ( theData.good() == true ) {
355 hasExactMF6=true;
356 theMF6FinalState.Init(theData);
357 //theData.close();
358 return;
359 }
360 //TK110430 END
361
362
363 G4String tString = "/FS";
364 G4bool dbool;
365 G4ParticleHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
366
367 G4String filename = aFile.GetName();
368 SetAZMs( A, Z, M, aFile );
369 //theBaseA = A;
370 //theBaseZ = G4int(Z+.5);
371 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
372 {
373 hasAnyData = false;
374 hasFSData = false;
375 hasXsec = false;
376 return;
377 }
378 //std::ifstream theData(filename, std::ios::in);
379 //std::istringstream theData(std::ios::in);
380 theData.clear();
382 hasFSData = theFinalStatePhotons.InitMean(theData);
383 if(hasFSData)
384 {
385 targetMass = theFinalStatePhotons.GetTargetMass();
386 theFinalStatePhotons.InitAngular(theData);
387 theFinalStatePhotons.InitEnergies(theData);
388 }
389 //theData.close();
390 }
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
@ stopAndKill
CLHEP::HepLorentzVector G4LorentzVector
#define M(row, col)
std::vector< G4ReactionProduct * > G4ReactionProductVector
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
const G4double A[17]
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector unit() const
double mag() const
Hep3Vector vect() const
value_type & Get() const
Definition: G4Cache.hh:315
void Put(const value_type &val) const
Definition: G4Cache.hh:321
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
void SetMomentum(const G4ThreeVector &momentum)
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4HadSecondary * GetSecondary(size_t i)
const G4Material * GetMaterial() const
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4DynamicParticle * GetParticle()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4double GetIonMass(G4int Z, G4int A, G4int nL=0, G4int lvl=0) const
Definition: G4IonTable.cc:1517
G4double GetTemperature() const
Definition: G4Material.hh:177
static G4Neutron * Neutron()
Definition: G4Neutron.cc:103
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:118
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
G4HadFinalState * ApplyYourself(const G4HadProjectile &theTrack)
void Init(std::istream &aDataFile)
void SetProjectileRP(G4ReactionProduct &aIncidentPart)
void SetTarget(G4ReactionProduct &aTarget)
G4ReactionProductVector * Sample(G4double anEnergy)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
G4Cache< G4HadFinalState * > theResult
static G4ParticleHPManager * GetInstance()
void GetDataStream(G4String, std::istringstream &iss)
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
void InitEnergies(std::istream &aDataFile)
G4ReactionProductVector * GetPhotons(G4double anEnergy)
void InitAngular(std::istream &aDataFile)
G4bool InitMean(std::istream &aDataFile)
void SetICM(G4bool) override
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
static G4int GetModelID(const G4int modelIndex)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void SetDefinitionAndUpdateE(const G4ParticleDefinition *aParticleDefinition)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
void SetKineticEnergy(const G4double en)
#define TRUE
Definition: globals.hh:41