Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4ParticleHPElasticFS.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// 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// P. Arce, June-2014 Conversion neutron_hp to particle_hp
35//
38
40#include "G4SystemOfUnits.hh"
41#include "G4ReactionProduct.hh"
42#include "G4Nucleus.hh"
43#include "G4Proton.hh"
44#include "G4Deuteron.hh"
45#include "G4Triton.hh"
46#include "G4Alpha.hh"
47#include "G4ThreeVector.hh"
48#include "G4LorentzVector.hh"
49#include "G4IonTable.hh"
51#include "G4Pow.hh"
52#include "zlib.h"
54
55
57{
58 secID = G4PhysicsModelCatalog::GetModelID( "model_NeutronHPElastic" );
59
60 hasXsec = false;
61 theCoefficients = 0;
62 theProbArray = 0;
63
64 repFlag = 0;
65 tE_of_repFlag3 = 0.0;
66 targetMass = 0.0;
67 frameFlag = 0;
68}
69
71 G4String& dirName, G4String&,
73{
74 G4String tString = "/FS";
75 G4bool dbool;
77 theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
78 G4String filename = aFile.GetName();
79 SetAZMs( A, Z, M, aFile );
80 //theBaseA = aFile.GetA();
81 //theBaseZ = aFile.GetZ();
82 if (!dbool) {
83 hasAnyData = false;
84 hasFSData = false;
85 hasXsec = false;
86 return;
87 }
88
89 //130205 For compressed data files
90 std::istringstream theData(std::ios::in);
92 //130205 END
93 theData >> repFlag >> targetMass >> frameFlag;
94
95 if (repFlag == 1) {
96 G4int nEnergy;
97 theData >> nEnergy;
98 theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
99 theCoefficients->InitInterpolation(theData);
100 G4double temp, energy;
101 G4int tempdep, nLegendre;
102 G4int i, ii;
103 for (i=0; i < nEnergy; i++) {
104 theData >> temp >> energy >> tempdep >> nLegendre;
105 energy *=eV;
106 theCoefficients->Init(i, energy, nLegendre);
107 theCoefficients->SetTemperature(i, temp);
108 G4double coeff = 0;
109 for (ii = 0; ii < nLegendre; ii++) {
110 // load legendre coefficients.
111 theData >> coeff;
112 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
113 }
114 }
115
116 } else if (repFlag == 2) {
117 G4int nEnergy;
118 theData >> nEnergy;
119 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
120 theProbArray->InitInterpolation(theData);
121 G4double temp, energy;
122 G4int tempdep, nPoints;
123 for (G4int i = 0; i < nEnergy; i++) {
124 theData >> temp >> energy >> tempdep >> nPoints;
125 energy *= eV;
126 theProbArray->InitInterpolation(i, theData);
127 theProbArray->SetT(i, temp);
128 theProbArray->SetX(i, energy);
129 G4double prob, costh;
130 for (G4int ii = 0; ii < nPoints; ii++) {
131 // fill probability arrays.
132 theData >> costh >> prob;
133 theProbArray->SetX(i, ii, costh);
134 theProbArray->SetY(i, ii, prob);
135 }
136 theProbArray->DoneSetXY( i );
137 }
138
139 } else if (repFlag == 3) {
140 G4int nEnergy_Legendre;
141 theData >> nEnergy_Legendre;
142 if (nEnergy_Legendre <= 0 ) {
143 std::stringstream iss;
144 iss << "G4ParticleHPElasticFS::Init Data Error repFlag is 3 but nEnergy_Legendre <= 0";
145 iss << "Z, A and M of problematic file is " << theNDLDataZ << ", "
146 << theNDLDataA << " and " << theNDLDataM << " respectively.";
147 throw G4HadronicException(__FILE__, __LINE__, iss.str() );
148 }
149 theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre );
150 theCoefficients->InitInterpolation( theData );
151 G4double temp, energy;
152 G4int tempdep, nLegendre;
153
154 for (G4int i = 0; i < nEnergy_Legendre; i++) {
155 theData >> temp >> energy >> tempdep >> nLegendre;
156 energy *=eV;
157 theCoefficients->Init( i , energy , nLegendre );
158 theCoefficients->SetTemperature( i , temp );
159 G4double coeff = 0;
160 for (G4int ii = 0; ii < nLegendre; ii++) {
161 // load legendre coefficients.
162 theData >> coeff;
163 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
164 }
165 }
166
167 tE_of_repFlag3 = energy;
168
169 G4int nEnergy_Prob;
170 theData >> nEnergy_Prob;
171 theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob );
172 theProbArray->InitInterpolation( theData );
173 G4int nPoints;
174 for (G4int i = 0; i < nEnergy_Prob; i++) {
175 theData >> temp >> energy >> tempdep >> nPoints;
176 energy *= eV;
177
178 // consistency check
179 if (i == 0)
180 //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
181 if (std::abs(energy - tE_of_repFlag3) / tE_of_repFlag3 > 1.0e-15)
182 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
183
184 theProbArray->InitInterpolation( i , theData );
185 theProbArray->SetT( i , temp );
186 theProbArray->SetX( i , energy );
187 G4double prob, costh;
188 for (G4int ii = 0; ii < nPoints; ii++) {
189 // fill probability arrays.
190 theData >> costh >> prob;
191 theProbArray->SetX( i , ii , costh );
192 theProbArray->SetY( i , ii , prob );
193 }
194 theProbArray->DoneSetXY( i );
195 }
196
197 } else if (repFlag==0) {
198 theData >> frameFlag;
199
200 } else {
201 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
202 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
203 }
204 //130205 For compressed data files(theData changed from ifstream to istringstream)
205 //theData.close();
206}
207
208
211{
212 if (theResult.Get() == NULL) theResult.Put(new G4HadFinalState);
213 theResult.Get()->Clear();
214 G4double eKinetic = theTrack.GetKineticEnergy();
215 const G4HadProjectile *incidentParticle = &theTrack;
216 G4ReactionProduct theNeutron(const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition() ));
217 theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect() );
218 theNeutron.SetKineticEnergy(eKinetic);
219
220 G4ReactionProduct theTarget;
221 G4Nucleus aNucleus;
222 G4ThreeVector neuVelo =
223 (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
224 theTarget =
225 aNucleus.GetBiasedThermalNucleus(targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
226
227 // Neutron and target defined as G4ReactionProducts
228 // Prepare Lorentz transformation to lab
229
230 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
231 G4double nEnergy = theNeutron.GetTotalEnergy();
232 G4ThreeVector the3Target = theTarget.GetMomentum();
233 G4double tEnergy = theTarget.GetTotalEnergy();
234 G4ReactionProduct theCMS;
235 G4double totE = nEnergy+tEnergy;
236 G4ThreeVector the3CMS = the3Target+the3Neutron;
237 theCMS.SetMomentum(the3CMS);
238 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
239 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
240 theCMS.SetMass(sqrts);
241 theCMS.SetTotalEnergy(totE);
242
243 // Data come as function of n-energy in nuclear rest frame
244 G4ReactionProduct boosted;
245 boosted.Lorentz(theNeutron, theTarget);
246 eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
247 G4double cosTh = -2;
248
249 if (repFlag == 1) {
250 cosTh = theCoefficients->SampleElastic(eKinetic);
251
252 } else if (repFlag == 2) {
253 cosTh = theProbArray->Sample(eKinetic);
254
255 } else if (repFlag == 3) {
256 if (eKinetic <= tE_of_repFlag3) {
257 cosTh = theCoefficients->SampleElastic(eKinetic);
258 } else {
259 cosTh = theProbArray->Sample(eKinetic);
260 }
261
262 } else if (repFlag == 0) {
263 cosTh = 2.*G4UniformRand() - 1.;
264
265 } else {
266 G4cout << "Unusable number for repFlag: repFlag=" << repFlag << G4endl;
267 throw G4HadronicException(__FILE__, __LINE__,
268 "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
269 }
270
271 if (cosTh < -1.1) { return 0; }
272
273 G4double phi = twopi*G4UniformRand();
274 G4double cosPhi = std::cos(phi);
275 G4double sinPhi = std::sin(phi);
276 G4double theta = std::acos(cosTh);
277 G4double sinth = std::sin(theta);
278
279 if (frameFlag == 1) {
280 // Projectile scattering values cosTh are in target rest frame
281 // In this frame, do relativistic calculation of scattered projectile and
282 // target 4-momenta
283
284 theNeutron.Lorentz(theNeutron, theTarget);
285 G4double mN = theNeutron.GetMass();
286 G4double Pinit = theNeutron.GetTotalMomentum(); // Incident momentum
287 G4double Einit = theNeutron.GetTotalEnergy(); // Incident energy
288 G4double mT = theTarget.GetMass();
289
290 G4double ratio = mT/mN;
291 G4double sqt = std::sqrt(ratio*ratio - 1.0 + cosTh*cosTh);
292 G4double beta = Pinit/(mT + Einit); // CMS beta
293 G4double denom = 1. - beta*beta*cosTh*cosTh;
294 G4double term1 = cosTh*(Einit*ratio + mN)/(mN*ratio + Einit);
295 G4double pN = beta*mN*(term1 + sqt)/denom;
296
297 // Get the scattered momentum and rotate it in theta and phi
298 G4ThreeVector pDir = theNeutron.GetMomentum()/Pinit;
299 G4double px = pN*pDir.x();
300 G4double py = pN*pDir.y();
301 G4double pz = pN*pDir.z();
302
303 G4ThreeVector pcmRot;
304 pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
305 pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
306 pcmRot.setZ(-px*sinth + pz*cosTh);
307 theNeutron.SetMomentum(pcmRot);
308 G4double eN = std::sqrt(pN*pN + mN*mN); // Scattered neutron energy
309 theNeutron.SetTotalEnergy(eN);
310
311 // Get the scattered target momentum
312 G4ReactionProduct toLab(-1.*theTarget);
313 theTarget.SetMomentum(pDir*Pinit - pcmRot);
314 G4double eT = Einit - eN + mT;
315 theTarget.SetTotalEnergy(eT);
316
317 // Now back to lab frame
318 theNeutron.Lorentz(theNeutron, toLab);
319 theTarget.Lorentz(theTarget, toLab);
320
321 //111005 Protection for not producing 0 kinetic energy target
322 if (theNeutron.GetKineticEnergy() <= 0)
323 theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
324 if (theTarget.GetKineticEnergy() <= 0)
325 theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
326
327 } else if (frameFlag == 2) {
328 // Projectile scattering values cosTh taken from center of mass tabulation
329
330 G4LorentzVector proj(nEnergy, the3Neutron);
331 G4LorentzVector targ(tEnergy, the3Target);
332 G4ThreeVector boostToCM = proj.findBoostToCM(targ);
333 proj.boost(boostToCM);
334 targ.boost(boostToCM);
335
336 // Rotate projectile and target momenta by CM scattering angle
337 // Note: at this point collision axis is not along z axis, due to
338 // momentum given target nucleus by thermal process
339 G4double px = proj.px();
340 G4double py = proj.py();
341 G4double pz = proj.pz();
342
343 G4ThreeVector pcmRot;
344 pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
345 pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
346 pcmRot.setZ(-px*sinth + pz*cosTh);
347 proj.setVect(pcmRot);
348 targ.setVect(-pcmRot);
349
350 // Back to lab frame
351 proj.boost(-boostToCM);
352 targ.boost(-boostToCM);
353
354 theNeutron.SetMomentum(proj.vect() );
355 theNeutron.SetTotalEnergy(proj.e() );
356
357 theTarget.SetMomentum(targ.vect() );
358 theTarget.SetTotalEnergy(targ.e() );
359
360 //080904 Add Protection for very low energy (1e-6eV) scattering
361 if (theNeutron.GetKineticEnergy() <= 0) {
362 theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
363 }
364
365 //080904 Add Protection for very low energy (1e-6eV) scattering
366 if (theTarget.GetKineticEnergy() <= 0) {
367 theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
368 }
369
370 } else {
371 G4cout << "Value of frameFlag (1=LAB, 2=CMS): " << frameFlag;
372 throw G4HadronicException(__FILE__, __LINE__,
373 "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
374 }
375
376 // Everything is now in the lab frame
377 // Set energy change and momentum change
380
381 // Make recoil a G4DynamicParticle
382 G4DynamicParticle* theRecoil = new G4DynamicParticle;
383 theRecoil->SetDefinition(G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ),
384 static_cast<G4int>(theBaseA), 0) );
385 theRecoil->SetMomentum(theTarget.GetMomentum());
386 theResult.Get()->AddSecondary(theRecoil, secID);
387
388 // Postpone the tracking of the primary neutron
390 return theResult.Get();
391}
392
@ suspend
#define M(row, col)
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 G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
Hep3Vector unit() const
double x() const
void setY(double)
void setZ(double)
void setX(double)
HepLorentzVector & boost(double, double, double)
Hep3Vector vect() const
void setVect(const Hep3Vector &)
Hep3Vector findBoostToCM() 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)
void SetMomentum(const G4ThreeVector &momentum)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void SetEnergyChange(G4double anEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
const G4Material * GetMaterial() const
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
static G4IonTable * GetIonTable()
Definition: G4IonTable.cc:170
G4double GetTemperature() const
Definition: G4Material.hh:177
G4ReactionProduct GetBiasedThermalNucleus(G4double aMass, G4ThreeVector aVelocity, G4double temp=-1) const
Definition: G4Nucleus.cc:118
G4HadFinalState * ApplyYourself(const G4HadProjectile &theTrack)
void Init(G4double A, G4double Z, G4int M, G4String &dirName, G4String &aFSType, G4ParticleDefinition *)
void SetAZMs(G4double anA, G4double aZ, G4int aM, G4ParticleHPDataUsed used)
G4Cache< G4HadFinalState * > theResult
void SetCoeff(G4int i, G4int l, G4double coeff)
void Init(G4int i, G4double e, G4int n)
void SetTemperature(G4int i, G4double temp)
void InitInterpolation(std::istream &aDataFile)
G4double SampleElastic(G4double anEnergy)
static G4ParticleHPManager * GetInstance()
void GetDataStream(G4String, std::istringstream &iss)
G4ParticleHPDataUsed GetName(G4int A, G4int Z, G4String base, G4String rest, G4bool &active)
G4double Sample(G4double x)
void SetT(G4int i, G4double x)
void SetX(G4int i, G4double x)
void SetY(G4int i, G4int j, G4double y)
void InitInterpolation(G4int i, std::istream &aDataFile)
static G4int GetModelID(const G4int modelIndex)
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
void SetMomentum(const G4double x, const G4double y, const G4double z)
void SetTotalEnergy(const G4double en)
G4double GetTotalMomentum() const
G4double GetKineticEnergy() const
G4double GetTotalEnergy() const
G4ThreeVector GetMomentum() const
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
void SetKineticEnergy(const G4double en)
G4double GetMass() const
void SetMass(const G4double mas)