source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPCaptureFS.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 12.5 KB
Line 
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//
35#include "G4NeutronHPCaptureFS.hh"
36#include "G4Gamma.hh"
37#include "G4ReactionProduct.hh"
38#include "G4Nucleus.hh"
39#include "G4PhotonEvaporation.hh"
40#include "G4Fragment.hh"
41#include "G4ParticleTable.hh"
42#include "G4NeutronHPDataUsed.hh"
43
44 G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself(const G4HadProjectile & theTrack)
45 {
46
47 G4int i;
48 theResult.Clear();
49// prepare neutron
50 G4double eKinetic = theTrack.GetKineticEnergy();
51 const G4HadProjectile *incidentParticle = &theTrack;
52 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
53 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
54 theNeutron.SetKineticEnergy( eKinetic );
55
56// prepare target
57 G4ReactionProduct theTarget;
58 G4Nucleus aNucleus;
59 G4double eps = 0.0001;
60 if(targetMass<500*MeV)
61 targetMass = ( G4NucleiProperties::GetNuclearMass( static_cast<G4int>(theBaseA+eps) , static_cast<G4int>(theBaseZ+eps) )) /
62 G4Neutron::Neutron()->GetPDGMass();
63 G4ThreeVector neutronVelocity = 1./G4Neutron::Neutron()->GetPDGMass()*theNeutron.GetMomentum();
64 G4double temperature = theTrack.GetMaterial()->GetTemperature();
65 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass, neutronVelocity, temperature);
66
67// go to nucleus rest system
68 theNeutron.Lorentz(theNeutron, -1*theTarget);
69 eKinetic = theNeutron.GetKineticEnergy();
70
71// dice the photons
72
73 G4ReactionProductVector * thePhotons = 0;
74 if ( HasFSData() && !getenv ( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) )
75 {
76 thePhotons = theFinalStatePhotons.GetPhotons(eKinetic);
77 }
78 else
79 {
80 G4ThreeVector aCMSMomentum = theNeutron.GetMomentum()+theTarget.GetMomentum();
81 G4LorentzVector p4(aCMSMomentum, theTarget.GetTotalEnergy() + theNeutron.GetTotalEnergy());
82 G4Fragment nucleus(static_cast<G4int>(theBaseA+1), static_cast<G4int>(theBaseZ) ,p4);
83 G4PhotonEvaporation photonEvaporation;
84 // T. K. add
85 photonEvaporation.SetICM( TRUE );
86 G4FragmentVector* products = photonEvaporation.BreakItUp(nucleus);
87 G4FragmentVector::iterator i;
88 thePhotons = new G4ReactionProductVector;
89 for(i=products->begin(); i!=products->end(); i++)
90 {
91 G4ReactionProduct * theOne = new G4ReactionProduct;
92 // T. K. add
93 if ( (*i)->GetParticleDefinition() != 0 )
94 theOne->SetDefinition( (*i)->GetParticleDefinition() );
95 else
96 theOne->SetDefinition( G4Gamma::Gamma() ); // this definiion will be over writen
97
98 // T. K. comment out below line
99 //theOne->SetDefinition( G4Gamma::Gamma() );
100 G4ParticleTable* theTable = G4ParticleTable::GetParticleTable();
101 if((*i)->GetMomentum().mag() > 10*MeV)
102 theOne->SetDefinition(
103 theTable->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ)) );
104 theOne->SetMomentum( (*i)->GetMomentum().vect() ) ;
105 theOne->SetTotalEnergy( (*i)->GetMomentum().t() );
106 thePhotons->push_back(theOne);
107 delete *i;
108 }
109 delete products;
110 }
111
112
113
114// add them to the final state
115
116 G4int nPhotons = 0;
117 if(thePhotons!=0) nPhotons=thePhotons->size();
118 G4int nParticles = nPhotons;
119 if(1==nPhotons) nParticles = 2;
120
121
122//Make at least one photon
123//101203 TK
124 if ( nPhotons == 0 )
125 {
126 G4ReactionProduct * theOne = new G4ReactionProduct;
127 theOne->SetDefinition( G4Gamma::Gamma() );
128 G4double theta = pi*G4UniformRand();
129 G4double phi = twopi*G4UniformRand();
130 G4double sinth = std::sin(theta);
131 G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) );
132 theOne->SetMomentum( direction ) ;
133 thePhotons->push_back(theOne);
134 nPhotons++; // 0 -> 1
135 }
136//One photon case: energy set to Q-value
137//101203 TK
138 if ( nPhotons == 1 )
139 {
140 G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
141 G4double Q = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass()
142 - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass();
143 thePhotons->operator[](0)->SetMomentum( Q*direction );
144 }
145//
146
147 // back to lab system
148 for(i=0; i<nPhotons; i++)
149 {
150 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), theTarget);
151 }
152
153 // Recoil, if only one gamma
154 if (1==nPhotons)
155 {
156 G4DynamicParticle * theOne = new G4DynamicParticle;
157 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
158 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
159 theOne->SetDefinition(aRecoil);
160 // Now energy;
161 // Can be done slightly better @
162 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect()
163 +theTarget.GetMomentum()
164 -thePhotons->operator[](0)->GetMomentum();
165
166 G4ThreeVector theMomUnit = aMomentum.unit();
167 G4double aKinEnergy = theTrack.GetKineticEnergy()
168 +theTarget.GetKineticEnergy(); // gammas come from Q-value
169 G4double theResMass = aRecoil->GetPDGMass();
170 G4double theResE = aRecoil->GetPDGMass()+aKinEnergy;
171 G4double theAbsMom = std::sqrt(theResE*theResE - theResMass*theResMass);
172 G4ThreeVector theMomentum = theAbsMom*theMomUnit;
173 theOne->SetMomentum(theMomentum);
174 theResult.AddSecondary(theOne);
175 }
176
177 // Now fill in the gammas.
178 for(i=0; i<nPhotons; i++)
179 {
180 // back to lab system
181 G4DynamicParticle * theOne = new G4DynamicParticle;
182 theOne->SetDefinition(thePhotons->operator[](i)->GetDefinition());
183 theOne->SetMomentum(thePhotons->operator[](i)->GetMomentum());
184 theResult.AddSecondary(theOne);
185 delete thePhotons->operator[](i);
186 }
187 delete thePhotons;
188
189//101203TK
190 G4bool residual = false;
191 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
192 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
193 for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
194 {
195 if ( theResult.GetSecondary(i)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
196 }
197
198 if ( residual == false )
199 {
200 G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
201 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
202 G4int nNonZero = 0;
203 G4LorentzVector p_photons(0,0,0,0);
204 for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
205 {
206 p_photons += theResult.GetSecondary(i)->GetParticle()->Get4Momentum();
207 // To many 0 momentum photons -> Check PhotonDist
208 if ( theResult.GetSecondary(i)->GetParticle()->Get4Momentum() > 0 ) nNonZero++;
209 }
210
211 // Can we include kinetic energy here?
212 G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
213 - ( p_photons.e() + aRecoil->GetPDGMass() );
214
215//Add photons
216 if ( nPhotons - nNonZero > 0 )
217 {
218 //G4cout << "TKDB G4NeutronHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl;
219 std::vector<G4double> vRand;
220 vRand.push_back( 0.0 );
221 for ( G4int i = 0 ; i != nPhotons - nNonZero - 1 ; i++ )
222 {
223 vRand.push_back( G4UniformRand() );
224 }
225 vRand.push_back( 1.0 );
226 std::sort( vRand.begin(), vRand.end() );
227
228 std::vector<G4double> vEPhoton;
229 for ( G4int i = 0 ; i < (G4int)vRand.size() - 1 ; i++ )
230 {
231 vEPhoton.push_back( deltaE * ( vRand[i+1] - vRand[i] ) );
232 }
233 std::sort( vEPhoton.begin(), vEPhoton.end() );
234
235 for ( G4int i = 0 ; i < nPhotons - nNonZero - 1 ; i++ )
236 {
237 //Isotopic in LAB OK?
238 G4double theta = pi*G4UniformRand();
239 G4double phi = twopi*G4UniformRand();
240 G4double sinth = std::sin(theta);
241 G4double en = vEPhoton[i];
242 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
243
244 p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
245 G4DynamicParticle * theOne = new G4DynamicParticle;
246 theOne->SetDefinition( G4Gamma::Gamma() );
247 theOne->SetMomentum( tempVector );
248 theResult.AddSecondary(theOne);
249 }
250
251// Add last photon
252 G4DynamicParticle * theOne = new G4DynamicParticle;
253 theOne->SetDefinition( G4Gamma::Gamma() );
254// For better momentum conservation
255 G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
256 p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
257 theOne->SetMomentum( lastPhoton );
258 theResult.AddSecondary(theOne);
259 }
260
261//Add residual
262 G4DynamicParticle * theOne = new G4DynamicParticle;
263 G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
264 - p_photons.vect();
265 theOne->SetDefinition(aRecoil);
266 theOne->SetMomentum( aMomentum );
267 theResult.AddSecondary(theOne);
268
269 }
270//101203TK END
271
272// clean up the primary neutron
273 theResult.SetStatusChange(stopAndKill);
274 return &theResult;
275 }
276
277 void G4NeutronHPCaptureFS::Init (G4double A, G4double Z, G4String & dirName, G4String & )
278 {
279 G4String tString = "/FS/";
280 G4bool dbool;
281 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), dirName, tString, dbool);
282 G4String filename = aFile.GetName();
283 theBaseA = A;
284 theBaseZ = G4int(Z+.5);
285 if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
286 {
287 hasAnyData = false;
288 hasFSData = false;
289 hasXsec = false;
290 return;
291 }
292 std::ifstream theData(filename, std::ios::in);
293
294 hasFSData = theFinalStatePhotons.InitMean(theData);
295 if(hasFSData)
296 {
297 targetMass = theFinalStatePhotons.GetTargetMass();
298 theFinalStatePhotons.InitAngular(theData);
299 theFinalStatePhotons.InitEnergies(theData);
300 }
301 theData.close();
302 }
Note: See TracBrowser for help on using the repository browser.