ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPFissionFS.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ParticleHPFissionFS.cc
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;
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;
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  // Protect against the very rare case of division by zero
228  G4double time = 0.0;
229  if ( theDecayConstants[i0-Prompt] > 1.0e-30 ) {
230  time = -G4Log(G4UniformRand())/theDecayConstants[i0-Prompt];
231  } else {
233  ed << " theDecayConstants[i0-Prompt]=" << theDecayConstants[i0-Prompt]
234  << " -> cannot sample the time : set it to 0.0 !" << G4endl;
235  G4Exception( "G4ParticleHPFissionFS::ApplyYourself ", "HAD_FISSIONHP_001", JustWarning, ed );
236  }
237 
238  time += theTrack.GetGlobalTime();
239  theResult.Get()->AddSecondary(theNeutrons->operator[](i0));
240  theResult.Get()->GetSecondary(theResult.Get()->GetNumberOfSecondaries()-1)->SetTime(time);
241  }
242  delete theNeutrons;
243  }
244  delete [] theDecayConstants;
245 // cout << "all delayed "<<delayed<<G4endl;
246  unsigned int nPhotons = 0;
247  if(thePhotons!=0)
248  {
249  nPhotons = thePhotons->size();
250  for(i=0; i<nPhotons; i++)
251  {
252  theResult.Get()->AddSecondary(thePhotons->operator[](i));
253  }
254  delete thePhotons;
255  }
256 
257 // finally deal with local energy depositions.
258 // G4cout <<"Number of secondaries = "<<theResult.GetNumberOfSecondaries()<< G4endl;
259 // G4cout <<"Number of photons = "<<nPhotons<<G4endl;
260 // G4cout <<"Number of Prompt = "<<Prompt<<G4endl;
261 // G4cout <<"Number of delayed = "<<delayed<<G4endl;
262 
264  G4double eDepByFragments = theERelease->GetFragmentKinetic();
265  //theResult.SetLocalEnergyDeposit(eDepByFragments);
266  if ( !produceFissionFragments ) theResult.Get()->SetLocalEnergyDeposit(eDepByFragments);
267 // cout << "local energy deposit" << eDepByFragments<<G4endl;
268 // clean up the primary neutron
270  //G4cout << "Prompt = " << Prompt << ", Delayed = " << delayed << ", All= " << all << G4endl;
271  //G4cout << "local energy deposit " << eDepByFragments/MeV << "MeV " << G4endl;
272 
273  //TKWORK120531
275  {
276  G4int fragA_Z=0;
277  G4int fragA_A=0;
278  G4int fragA_M=0;
279  // System is traget rest!
280  theFF.GetAFissionFragment(eKinetic,fragA_Z,fragA_A,fragA_M);
281  G4int fragB_Z=(G4int)theBaseZ-fragA_Z;
282  G4int fragB_A=(G4int)theBaseA-fragA_A-Prompt;
283  //fragA_M ignored
284  //G4int fragB_M=theBaseM-fragA_M;
285  //G4cout << fragA_Z << " " << fragA_A << " " << fragA_M << G4endl;
286  //G4cout << fragB_Z << " " << fragB_A << G4endl;
287 
289  //Excitation energy is not taken into account
290  G4ParticleDefinition* pdA = pt->GetIon( fragA_Z , fragA_A , 0.0 );
291  G4ParticleDefinition* pdB = pt->GetIon( fragB_Z , fragB_A , 0.0 );
292 
293  //Isotropic Distribution
295  // Bug #1745 DHW G4double theta = pi*G4UniformRand();
296  G4double costheta = 2.*G4UniformRand()-1.;
297  G4double theta = std::acos(costheta);
298  G4double sinth = std::sin(theta);
299  G4ThreeVector direction(sinth*std::cos(phi), sinth*std::sin(phi), costheta);
300 
301  // Just use ENDF value for this
302  G4double ER = eDepByFragments;
303  G4double ma = pdA->GetPDGMass();
304  G4double mb = pdB->GetPDGMass();
305  G4double EA = ER / ( 1 + ma/mb);
306  G4double EB = ER - EA;
307  G4DynamicParticle* dpA = new G4DynamicParticle( pdA , direction , EA);
308  G4DynamicParticle* dpB = new G4DynamicParticle( pdB , -direction , EB);
309  theResult.Get()->AddSecondary(dpA);
310  theResult.Get()->AddSecondary(dpB);
311  }
312  //TKWORK 120531 END
313 
314  return theResult.Get();
315  }