Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ParticleHPFissionFS.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-Apr-06 fix in delayed neutron and photon emission without FS data by T. Koi
31 // 07-Sep-11 M. Kelsey -- Follow change to G4HadFinalState interface
32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 //
34 
35 #include "G4Exp.hh"
36 #include "G4ParticleHPFissionFS.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4Nucleus.hh"
41 #include "G4IonTable.hh"
42 
44  {
45  //G4cout << "G4ParticleHPFissionFS::Init " << A << " " << Z << " " << M << G4endl;
46  theFS.Init(A, Z, M, dirName, aFSType, projectile);
47  theFC.Init(A, Z, M, dirName, aFSType, projectile);
48  theSC.Init(A, Z, M, dirName, aFSType, projectile);
49  theTC.Init(A, Z, M, dirName, aFSType, projectile);
50  theLC.Init(A, Z, M, dirName, aFSType, projectile);
51 
52  theFF.Init(A, Z, M, dirName, aFSType, projectile);
53  if ( G4ParticleHPManager::GetInstance()->GetProduceFissionFragments() && theFF.HasFSData() )
54  {
55  G4cout << "Fission fragment production is now activated in HP package for "
56  << "Z = " << (G4int)Z
57  << ", A = " << (G4int)A
58  //<< "M = " << M
59  << G4endl;
60  G4cout << "As currently modeled this option precludes production of delayed neutrons from fission fragments." << G4endl;
61  produceFissionFragments = true;
62  }
63  }
65  {
66 
67  //Because it may change by UI command
69 
70  //G4cout << "G4ParticleHPFissionFS::ApplyYourself " << G4endl;
71 // prepare neutron
72 
73  if ( theResult.Get() == NULL ) theResult.Put( new G4HadFinalState );
74  theResult.Get()->Clear();
75  G4double eKinetic = theTrack.GetKineticEnergy();
76  const G4HadProjectile *incidentParticle = &theTrack;
77  G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
78  theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
79  theNeutron.SetKineticEnergy( eKinetic );
80 
81 // prepare target
82  G4Nucleus aNucleus;
83  G4ReactionProduct theTarget;
84  G4double targetMass = theFS.GetMass();
85  G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
86  theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
87  theTarget.SetDefinition( G4IonTable::GetIonTable()->GetIon( G4int(theBaseZ), G4int(theBaseA) , 0.0 ) ); //TESTPHP
88 // set neutron and target in the FS classes
89  theFS.SetNeutronRP(theNeutron);
90  theFS.SetTarget(theTarget);
91  theFC.SetNeutronRP(theNeutron);
92  theFC.SetTarget(theTarget);
93  theSC.SetNeutronRP(theNeutron);
94  theSC.SetTarget(theTarget);
95  theTC.SetNeutronRP(theNeutron);
96  theTC.SetTarget(theTarget);
97  theLC.SetNeutronRP(theNeutron);
98  theLC.SetTarget(theTarget);
99 
100 
101  theFF.SetNeutronRP(theNeutron);
102  theFF.SetTarget(theTarget);
103 
104 //TKWORK 120531
105 //G4cout << theTarget.GetDefinition() << G4endl; this should be NULL
106 //G4cout << "Z = " << theBaseZ << ", A = " << theBaseA << ", M = " << theBaseM << G4endl;
107 // theNDLDataZ,A,M should be filled in each FS (theFS, theFC, theSC, theTC, theLC and theFF)
109 
110 // boost to target rest system and decide on channel.
111  theNeutron.Lorentz(theNeutron, -1*theTarget);
112 
113 // dice the photons
114 
115  G4DynamicParticleVector * thePhotons;
116  thePhotons = theFS.GetPhotons();
117 
118 // select the FS in charge
119 
120  eKinetic = theNeutron.GetKineticEnergy();
121  G4double xSec[4];
122  xSec[0] = theFC.GetXsec(eKinetic);
123  xSec[1] = xSec[0]+theSC.GetXsec(eKinetic);
124  xSec[2] = xSec[1]+theTC.GetXsec(eKinetic);
125  xSec[3] = xSec[2]+theLC.GetXsec(eKinetic);
126  G4int it;
127  unsigned int i=0;
128  G4double random = G4UniformRand();
129  if(xSec[3]==0)
130  {
131  it=-1;
132  }
133  else
134  {
135  for(i=0; i<4; i++)
136  {
137  it =i;
138  if(random<xSec[i]/xSec[3]) break;
139  }
140  }
141 
142 // dice neutron multiplicities, energies and momenta in Lab. @@
143 // no energy conservation on an event-to-event basis. we rely on the data to be ok. @@
144 // also for mean, we rely on the consistancy of the data. @@
145 
146  G4int Prompt=0, delayed=0, all=0;
147  G4DynamicParticleVector * theNeutrons = 0;
148  switch(it) // check logic, and ask, if partials can be assumed to correspond to individual particles @@@
149  {
150  case 0:
151  theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 0);
152  if(Prompt==0&&delayed==0) Prompt=all;
153  theNeutrons = theFC.ApplyYourself(Prompt); // delayed always in FS
154  // take 'U' into account explicitely (see 5.4) in the sampling of energy @@@@
155  break;
156  case 1:
157  theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 1);
158  if(Prompt==0&&delayed==0) Prompt=all;
159  theNeutrons = theSC.ApplyYourself(Prompt); // delayed always in FS, off done in FSFissionFS
160  break;
161  case 2:
162  theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 2);
163  if(Prompt==0&&delayed==0) Prompt=all;
164  theNeutrons = theTC.ApplyYourself(Prompt); // delayed always in FS
165  break;
166  case 3:
167  theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 3);
168  if(Prompt==0&&delayed==0) Prompt=all;
169  theNeutrons = theLC.ApplyYourself(Prompt); // delayed always in FS
170  break;
171  default:
172  break;
173  }
174 
175 // dice delayed neutrons and photons, and fallback
176 // for Prompt in case channel had no FS data; add all paricles to FS.
177 
178  //TKWORK120531
179  if ( produceFissionFragments ) delayed=0;
180 
181  G4double * theDecayConstants;
182 
183  if( theNeutrons != 0)
184  {
185  theDecayConstants = new G4double[delayed];
186  //
187  //110527TKDB Unused codes, Detected by gcc4.6 compiler
188  //G4int nPhotons = 0;
189  //if(thePhotons!=0) nPhotons = thePhotons->size();
190  for(i=0; i<theNeutrons->size(); i++)
191  {
192  theResult.Get()->AddSecondary(theNeutrons->operator[](i));
193  }
194  delete theNeutrons;
195 
196  G4DynamicParticleVector * theDelayed = 0;
197 // G4cout << "delayed" << G4endl;
198  theDelayed = theFS.ApplyYourself(0, delayed, theDecayConstants);
199  for(i=0; i<theDelayed->size(); i++)
200  {
201  G4double time = -G4Log(G4UniformRand())/theDecayConstants[i];
202  time += theTrack.GetGlobalTime();
203  theResult.Get()->AddSecondary(theDelayed->operator[](i));
204  theResult.Get()->GetSecondary(theResult.Get()->GetNumberOfSecondaries()-1)->SetTime(time);
205  }
206  delete theDelayed;
207  }
208  else
209  {
210 // cout << " all = "<<all<<G4endl;
211  theFS.SampleNeutronMult(all, Prompt, delayed, eKinetic, 0);
212  theDecayConstants = new G4double[delayed];
213  if(Prompt==0&&delayed==0) Prompt=all;
214  theNeutrons = theFS.ApplyYourself(Prompt, delayed, theDecayConstants);
215  //110527TKDB Unused codes, Detected by gcc4.6 compiler
216  //G4int nPhotons = 0;
217  //if(thePhotons!=0) nPhotons = thePhotons->size();
218  G4int i0;
219  for(i0=0; i0<Prompt; i0++)
220  {
221  theResult.Get()->AddSecondary(theNeutrons->operator[](i0));
222  }
223 
224 //G4cout << "delayed" << G4endl;
225  for(i0=Prompt; i0<Prompt+delayed; i0++)
226  {
227  G4double time = -G4Log(G4UniformRand())/theDecayConstants[i0-Prompt];
228  time += theTrack.GetGlobalTime();
229  theResult.Get()->AddSecondary(theNeutrons->operator[](i0));
230  theResult.Get()->GetSecondary(theResult.Get()->GetNumberOfSecondaries()-1)->SetTime(time);
231  }
232  delete theNeutrons;
233  }
234  delete [] theDecayConstants;
235 // cout << "all delayed "<<delayed<<G4endl;
236  unsigned int nPhotons = 0;
237  if(thePhotons!=0)
238  {
239  nPhotons = thePhotons->size();
240  for(i=0; i<nPhotons; i++)
241  {
242  theResult.Get()->AddSecondary(thePhotons->operator[](i));
243  }
244  delete thePhotons;
245  }
246 
247 // finally deal with local energy depositions.
248 // G4cout <<"Number of secondaries = "<<theResult.GetNumberOfSecondaries()<< G4endl;
249 // G4cout <<"Number of photons = "<<nPhotons<<G4endl;
250 // G4cout <<"Number of Prompt = "<<Prompt<<G4endl;
251 // G4cout <<"Number of delayed = "<<delayed<<G4endl;
252 
253  G4ParticleHPFissionERelease * theERelease = theFS.GetEnergyRelease();
254  G4double eDepByFragments = theERelease->GetFragmentKinetic();
255  //theResult.SetLocalEnergyDeposit(eDepByFragments);
256  if ( !produceFissionFragments ) theResult.Get()->SetLocalEnergyDeposit(eDepByFragments);
257 // cout << "local energy deposit" << eDepByFragments<<G4endl;
258 // clean up the primary neutron
260  //G4cout << "Prompt = " << Prompt << ", Delayed = " << delayed << ", All= " << all << G4endl;
261  //G4cout << "local energy deposit " << eDepByFragments/MeV << "MeV " << G4endl;
262 
263  //TKWORK120531
264  if ( produceFissionFragments )
265  {
266  G4int fragA_Z=0;
267  G4int fragA_A=0;
268  G4int fragA_M=0;
269  // System is traget rest!
270  theFF.GetAFissionFragment(eKinetic,fragA_Z,fragA_A,fragA_M);
271  G4int fragB_Z=(G4int)theBaseZ-fragA_Z;
272  G4int fragB_A=(G4int)theBaseA-fragA_A-Prompt;
273  //fragA_M ignored
274  //G4int fragB_M=theBaseM-fragA_M;
275  //G4cout << fragA_Z << " " << fragA_A << " " << fragA_M << G4endl;
276  //G4cout << fragB_Z << " " << fragB_A << G4endl;
277 
279  //Excitation energy is not taken into account
280  G4ParticleDefinition* pdA = pt->GetIon( fragA_Z , fragA_A , 0.0 );
281  G4ParticleDefinition* pdB = pt->GetIon( fragB_Z , fragB_A , 0.0 );
282 
283  //Isotropic Distribution
284  G4double phi = twopi*G4UniformRand();
285  G4double theta = pi*G4UniformRand();
286  G4double sinth = std::sin(theta);
287  G4ThreeVector direction (sinth*std::cos(phi) , sinth*std::sin(phi), std::cos(theta) );
288 
289  // Just use ENDF value for this
290  G4double ER = eDepByFragments;
291  G4double ma = pdA->GetPDGMass();
292  G4double mb = pdB->GetPDGMass();
293  G4double EA = ER / ( 1 + ma/mb);
294  G4double EB = ER - EA;
295  G4DynamicParticle* dpA = new G4DynamicParticle( pdA , direction , EA);
296  G4DynamicParticle* dpB = new G4DynamicParticle( pdB , -direction , EB);
297  theResult.Get()->AddSecondary(dpA);
298  theResult.Get()->AddSecondary(dpB);
299  }
300  //TKWORK 120531 END
301 
302  return theResult.Get();
303  }
static G4ParticleHPManager * GetInstance()
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
G4Cache< G4HadFinalState * > theResult
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
G4HadSecondary * GetSecondary(size_t i)
G4DynamicParticleVector * ApplyYourself(G4int NNeutrons)
void SetKineticEnergy(const G4double en)
void SetMomentum(const G4double x, const G4double y, const G4double z)
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
void SetNeutronRP(const G4ReactionProduct &aNeutron)
value_type & Get() const
Definition: G4Cache.hh:282
void SampleNeutronMult(G4int &all, G4int &Prompt, G4int &delayed, G4double energy, G4int off)
void SetNeutronRP(const G4ReactionProduct &aNeutron)
void SetTarget(const G4ReactionProduct &aTarget)
int G4int
Definition: G4Types.hh:78
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetStatusChange(G4HadFinalStateStatus aS)
void GetAFissionFragment(G4double, G4int &, G4int &, G4int &)
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:113
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
virtual G4double GetXsec(G4double anEnergy)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
G4double GetKineticEnergy() const
G4ParticleHPFissionERelease * GetEnergyRelease()
G4double GetGlobalTime() const
G4DynamicParticleVector * ApplyYourself(G4int Prompt, G4int delayed, G4double *decayconst)
std::vector< G4DynamicParticle * > G4DynamicParticleVector
G4DynamicParticleVector * ApplyYourself(G4int NNeutrons)
G4DynamicParticleVector * GetPhotons()
const G4LorentzVector & Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.hh:78
G4double GetKineticEnergy() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
void SetTarget(const G4ReactionProduct &aTarget)
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetPDGMass() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
G4ThreeVector GetMomentum() const
G4double GetTemperature() const
Definition: G4Material.hh:182
void SetLocalEnergyDeposit(G4double aE)
#define G4endl
Definition: G4ios.hh:61
G4DynamicParticleVector * ApplyYourself(G4int NNeutrons)
const G4Material * GetMaterial() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
static constexpr double pi
Definition: G4SIunits.hh:75
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4DynamicParticleVector * ApplyYourself(G4int nNeutrons)
G4HadFinalState * ApplyYourself(const G4HadProjectile &theTrack)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *projectile)
double G4double
Definition: G4Types.hh:76
void Put(const value_type &val) const
Definition: G4Cache.hh:286
G4int GetNumberOfSecondaries() const
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)