source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPElasticFS.cc@ 1199

Last change on this file since 1199 was 962, checked in by garnier, 17 years ago

update processes
File size: 14.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// 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
31// is added by T. KOI
32// 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi
33//
34#include "G4NeutronHPElasticFS.hh"
35#include "G4ReactionProduct.hh"
36#include "G4Nucleus.hh"
37#include "G4Proton.hh"
38#include "G4Deuteron.hh"
39#include "G4Triton.hh"
40#include "G4Alpha.hh"
41#include "G4ThreeVector.hh"
42#include "G4LorentzVector.hh"
43#include "G4ParticleTable.hh"
44#include "G4NeutronHPDataUsed.hh"
45
46 void G4NeutronHPElasticFS::Init (G4double A, G4double Z, G4String & dirName, G4String & )
47 {
48 G4String tString = "/FS/";
49 G4bool dbool;
50 G4NeutronHPDataUsed aFile = theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), dirName, tString, dbool);
51 G4String filename = aFile.GetName();
52 theBaseA = aFile.GetA();
53 theBaseZ = aFile.GetZ();
54 if(!dbool)
55 {
56 hasAnyData = false;
57 hasFSData = false;
58 hasXsec = false;
59 return;
60 }
61 std::ifstream theData(filename, std::ios::in);
62 theData >> repFlag >> targetMass >> frameFlag;
63 if(repFlag==1)
64 {
65 G4int nEnergy;
66 theData >> nEnergy;
67 theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
68 theCoefficients->InitInterpolation(theData);
69 G4double temp, energy;
70 G4int tempdep, nLegendre;
71 G4int i, ii;
72 for (i=0; i<nEnergy; i++)
73 {
74 theData >> temp >> energy >> tempdep >> nLegendre;
75 energy *=eV;
76 theCoefficients->Init(i, energy, nLegendre);
77 theCoefficients->SetTemperature(i, temp);
78 G4double coeff=0;
79 for(ii=0; ii<nLegendre; ii++)
80 {
81 // load legendre coefficients.
82 theData >> coeff;
83 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
84 }
85 }
86 }
87 else if (repFlag==2)
88 {
89 G4int nEnergy;
90 theData >> nEnergy;
91 theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
92 theProbArray->InitInterpolation(theData);
93 G4double temp, energy;
94 G4int tempdep, nPoints;
95 for(G4int i=0; i<nEnergy; i++)
96 {
97 theData >> temp >> energy >> tempdep >> nPoints;
98 energy *= eV;
99 theProbArray->InitInterpolation(i, theData);
100 theProbArray->SetT(i, temp);
101 theProbArray->SetX(i, energy);
102 G4double prob, costh;
103 for(G4int ii=0; ii<nPoints; ii++)
104 {
105 // fill probability arrays.
106 theData >> costh >> prob;
107 theProbArray->SetX(i, ii, costh);
108 theProbArray->SetY(i, ii, prob);
109 }
110 }
111 }
112 else if ( repFlag==3 )
113 {
114 G4int nEnergy_Legendre;
115 theData >> nEnergy_Legendre;
116 theCoefficients = new G4NeutronHPLegendreStore( nEnergy_Legendre );
117 theCoefficients->InitInterpolation( theData );
118 G4double temp, energy;
119 G4int tempdep, nLegendre;
120 G4int i, ii;
121 for ( i = 0 ; i < nEnergy_Legendre ; i++ )
122 {
123 theData >> temp >> energy >> tempdep >> nLegendre;
124 energy *=eV;
125 theCoefficients->Init( i , energy , nLegendre );
126 theCoefficients->SetTemperature( i , temp );
127 G4double coeff = 0;
128 for ( ii = 0 ; ii < nLegendre ; ii++ )
129 {
130 // load legendre coefficients.
131 theData >> coeff;
132 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
133 }
134 }
135
136 tE_of_repFlag3 = energy;
137
138 G4int nEnergy_Prob;
139 theData >> nEnergy_Prob;
140 theProbArray = new G4NeutronHPPartial( nEnergy_Prob , nEnergy_Prob );
141 theProbArray->InitInterpolation( theData );
142 G4int nPoints;
143 for ( G4int i=0 ; i < nEnergy_Prob ; i++ )
144 {
145 theData >> temp >> energy >> tempdep >> nPoints;
146
147 energy *= eV;
148
149// consistensy check
150 if ( i == 0 )
151 if ( energy != tE_of_repFlag3 )
152 G4cout << "Warning Trangition Energy of repFlag3 is not consistent." << G4endl;
153
154 theProbArray->InitInterpolation( i , theData );
155 theProbArray->SetT( i , temp );
156 theProbArray->SetX( i , energy );
157 G4double prob, costh;
158 for( G4int ii = 0 ; ii < nPoints ; ii++ )
159 {
160 // fill probability arrays.
161 theData >> costh >> prob;
162 theProbArray->SetX( i , ii , costh );
163 theProbArray->SetY( i , ii , prob );
164 }
165 }
166 }
167 else if (repFlag==0)
168 {
169 theData >> frameFlag;
170 }
171 else
172 {
173 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
174 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
175 }
176 theData.close();
177 }
178 G4HadFinalState * G4NeutronHPElasticFS::ApplyYourself(const G4HadProjectile & theTrack)
179 {
180// G4cout << "G4NeutronHPElasticFS::ApplyYourself+"<<G4endl;
181 theResult.Clear();
182 G4double eKinetic = theTrack.GetKineticEnergy();
183 const G4HadProjectile *incidentParticle = &theTrack;
184 G4ReactionProduct theNeutron( const_cast<G4ParticleDefinition *>(incidentParticle->GetDefinition()) );
185 theNeutron.SetMomentum( incidentParticle->Get4Momentum().vect() );
186 theNeutron.SetKineticEnergy( eKinetic );
187// G4cout << "G4NeutronHPElasticFS::ApplyYourself++"<<eKinetic<<" "<<G4endl;
188// G4cout << "CMSVALUES 0 "<<theNeutron.GetTotalMomentum()<<G4endl;
189
190 G4ReactionProduct theTarget;
191 G4Nucleus aNucleus;
192 G4ThreeVector neuVelo = (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
193 theTarget = aNucleus.GetBiasedThermalNucleus( targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
194// G4cout << "Nucleus-test"<<" "<<targetMass<<" ";
195// G4cout << theTarget.GetMomentum().x()<<" ";
196// G4cout << theTarget.GetMomentum().y()<<" ";
197// G4cout << theTarget.GetMomentum().z()<<G4endl;
198
199 // neutron and target defined as reaction products.
200
201// prepare lorentz-transformation to Lab.
202
203 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
204 G4double nEnergy = theNeutron.GetTotalEnergy();
205 G4ThreeVector the3Target = theTarget.GetMomentum();
206// cout << "@@@" << the3Target<<G4endl;
207 G4double tEnergy = theTarget.GetTotalEnergy();
208 G4ReactionProduct theCMS;
209 G4double totE = nEnergy+tEnergy;
210 G4ThreeVector the3CMS = the3Target+the3Neutron;
211 theCMS.SetMomentum(the3CMS);
212 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
213 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
214 theCMS.SetMass(sqrts);
215 theCMS.SetTotalEnergy(totE);
216
217 // data come as fcn of n-energy in nuclear rest frame
218 G4ReactionProduct boosted;
219 boosted.Lorentz(theNeutron, theTarget);
220 eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
221 G4double cosTh = -2;
222 if(repFlag == 1)
223 {
224 cosTh = theCoefficients->SampleElastic(eKinetic);
225 }
226
227 else if (repFlag==2)
228 {
229 cosTh = theProbArray->Sample(eKinetic);
230 }
231 else if (repFlag==3)
232 {
233 if ( eKinetic <= tE_of_repFlag3 )
234 {
235 cosTh = theCoefficients->SampleElastic(eKinetic);
236 }
237 else
238 {
239 cosTh = theProbArray->Sample(eKinetic);
240 }
241 }
242 else if (repFlag==0)
243 {
244 cosTh = 2.*G4UniformRand()-1.;
245 }
246 else
247 {
248 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
249 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::Init -- unusable number for repFlag");
250 }
251 if(cosTh<-1.1) { return 0; }
252 G4double phi = twopi*G4UniformRand();
253 G4double theta = std::acos(cosTh);
254 G4double sinth = std::sin(theta);
255 if (frameFlag == 1) // final state data given in target rest frame.
256 {
257 // we have the scattering angle, now we need the energy, then do the
258 // boosting.
259 // relativistic elastic scattering energy angular correlation:
260 theNeutron.Lorentz(theNeutron, theTarget);
261 G4double e0 = theNeutron.GetTotalEnergy();
262 G4double p0 = theNeutron.GetTotalMomentum();
263 G4double mN = theNeutron.GetMass();
264 G4double mT = theTarget.GetMass();
265 G4double eE = e0+mT;
266 G4double ap = (mT+eE)*(mT-eE) + (p0+mN)*(p0-mN);
267 G4double a = 4*(eE+p0*cosTh)*(eE-p0*cosTh);
268 G4double b = 4*ap*p0*cosTh;
269 G4double c = (2.*eE*mN-ap)*(2.*eE*mN+ap);
270 G4double en = (-b+std::sqrt(b*b - 4*a*c) )/(2*a);
271 G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
272 theNeutron.SetMomentum(tempVector);
273 theNeutron.SetTotalEnergy(std::sqrt(en*en+theNeutron.GetMass()*theNeutron.GetMass()));
274 // first to lab
275 theNeutron.Lorentz(theNeutron, -1.*theTarget);
276 // now to CMS
277 theNeutron.Lorentz(theNeutron, theCMS);
278 theTarget.SetMomentum(-theNeutron.GetMomentum());
279 theTarget.SetTotalEnergy(theNeutron.GetTotalEnergy());
280 // and back to lab
281 theNeutron.Lorentz(theNeutron, -1.*theCMS);
282 theTarget.Lorentz(theTarget, -1.*theCMS);
283 }
284 else if (frameFlag == 2) // CMS
285 {
286 theNeutron.Lorentz(theNeutron, theCMS);
287 theTarget.Lorentz(theTarget, theCMS);
288 G4double en = theNeutron.GetTotalMomentum(); // already in CMS.
289 G4ThreeVector cmsMom=theNeutron.GetMomentum(); // for neutron direction in CMS
290 G4double cms_theta=cmsMom.theta();
291 G4double cms_phi=cmsMom.phi();
292 G4ThreeVector tempVector;
293 tempVector.setX(std::cos(theta)*std::sin(cms_theta)*std::cos(cms_phi)
294 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::cos(cms_phi)
295 -std::sin(theta)*std::sin(phi)*std::sin(cms_phi) );
296 tempVector.setY(std::cos(theta)*std::sin(cms_theta)*std::sin(cms_phi)
297 +std::sin(theta)*std::cos(phi)*std::cos(cms_theta)*std::sin(cms_phi)
298 +std::sin(theta)*std::sin(phi)*std::cos(cms_phi) );
299 tempVector.setZ(std::cos(theta)*std::cos(cms_theta)
300 -std::sin(theta)*std::cos(phi)*std::sin(cms_theta) );
301 tempVector *= en;
302 theNeutron.SetMomentum(tempVector);
303 theTarget.SetMomentum(-tempVector);
304 G4double tP = theTarget.GetTotalMomentum();
305 G4double tM = theTarget.GetMass();
306 theTarget.SetTotalEnergy(std::sqrt((tP+tM)*(tP+tM)-2.*tP*tM));
307
308/*
309For debug purpose.
310Same transformation G4ReactionProduct.Lorentz() by 4vectors
311{
312G4LorentzVector n4p = G4LorentzVector ( theNeutron.GetMomentum() , theNeutron.GetKineticEnergy() + theNeutron.GetMass() );
313G4cout << "before " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
314G4LorentzVector cm4p = G4LorentzVector ( theCMS.GetMomentum() , theCMS.GetKineticEnergy() + theCMS.GetMass() );
315n4p.boost( cm4p.boostVector() );
316G4cout << cm4p/eV << G4endl;
317G4cout << "after " << ( n4p.e() - n4p.m() ) / eV<< G4endl;
318}
319*/
320
321 theNeutron.Lorentz(theNeutron, -1.*theCMS);
322//080904 Add Protection for very low energy (1e-6eV) scattering
323 if ( theNeutron.GetKineticEnergy() < 0 )
324 {
325 theNeutron.SetMomentum( G4ThreeVector(0) );
326 theNeutron.SetTotalEnergy ( theNeutron.GetMass() );
327 }
328
329 theTarget.Lorentz(theTarget, -1.*theCMS);
330//080904 Add Protection for very low energy (1e-6eV) scattering
331 if ( theTarget.GetKineticEnergy() < 0 )
332 {
333 theTarget.SetMomentum( G4ThreeVector(0) );
334 theTarget.SetTotalEnergy ( theTarget.GetMass() );
335 }
336 }
337 else
338 {
339 G4cout <<"Value of frameFlag (1=LAB, 2=CMS): "<<frameFlag;
340 throw G4HadronicException(__FILE__, __LINE__, "G4NeutronHPElasticFS::ApplyYourSelf frameflag incorrect");
341 }
342 // now all in Lab
343// nun den recoil generieren...und energy change, momentum change angeben.
344 theResult.SetEnergyChange(theNeutron.GetKineticEnergy());
345 theResult.SetMomentumChange(theNeutron.GetMomentum().unit());
346 G4DynamicParticle* theRecoil = new G4DynamicParticle;
347 if(targetMass<4.5)
348 {
349 if(targetMass<1)
350 {
351 // proton
352 theRecoil->SetDefinition(G4Proton::Proton());
353 }
354 else if(targetMass<2 )
355 {
356 // deuteron
357 theRecoil->SetDefinition(G4Deuteron::Deuteron());
358 }
359 else if(targetMass<2.999 )
360 {
361 // 3He
362 theRecoil->SetDefinition(G4He3::He3());
363 }
364 else if(targetMass<3 )
365 {
366 // Triton
367 theRecoil->SetDefinition(G4Triton::Triton());
368 }
369 else
370 {
371 // alpha
372 theRecoil->SetDefinition(G4Alpha::Alpha());
373 }
374 }
375 else
376 {
377 theRecoil->SetDefinition(G4ParticleTable::GetParticleTable()
378 ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ)));
379 }
380 theRecoil->SetMomentum(theTarget.GetMomentum());
381 theResult.AddSecondary(theRecoil);
382// G4cout << "G4NeutronHPElasticFS::ApplyYourself 10+"<<G4endl;
383 // postpone the tracking of the primary neutron
384 theResult.SetStatusChange(suspend);
385 return &theResult;
386 }
Note: See TracBrowser for help on using the repository browser.