Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4NeutronHPAngular Class Reference

#include <G4NeutronHPAngular.hh>

Public Member Functions

 G4NeutronHPAngular ()
 
 ~G4NeutronHPAngular ()
 
void Init (std::ifstream &aDataFile)
 
void SampleAndUpdate (G4ReactionProduct &aNeutron)
 
void SetTarget (const G4ReactionProduct &aTarget)
 
void SetNeutron (const G4ReactionProduct &aNeutron)
 
G4double GetTargetMass ()
 

Detailed Description

Definition at line 42 of file G4NeutronHPAngular.hh.

Constructor & Destructor Documentation

◆ G4NeutronHPAngular()

G4NeutronHPAngular::G4NeutronHPAngular ( )
inline

Definition at line 46 of file G4NeutronHPAngular.hh.

47 {
48//TKDB110511
49 //theAngularDistributionType = 0;
50 //theIsoFlag = false;
51 theIsoFlag = true;
52// TKDB
53 theCoefficients = 0;
54 theProbArray = 0;
55 }

◆ ~G4NeutronHPAngular()

G4NeutronHPAngular::~G4NeutronHPAngular ( )
inline

Definition at line 57 of file G4NeutronHPAngular.hh.

58 {
59// TKDB
60 delete theCoefficients;
61 delete theProbArray;
62 }

Member Function Documentation

◆ GetTargetMass()

G4double G4NeutronHPAngular::GetTargetMass ( )
inline

Definition at line 72 of file G4NeutronHPAngular.hh.

72{ return targetMass; }

Referenced by G4NeutronHPInelasticBaseFS::BaseApply().

◆ Init()

void G4NeutronHPAngular::Init ( std::ifstream &  aDataFile)

Definition at line 39 of file G4NeutronHPAngular.cc.

40{
41// G4cout << "here we are entering the Angular Init"<<G4endl;
42 aDataFile >> theAngularDistributionType >> targetMass;
43 aDataFile >> frameFlag;
44 if(theAngularDistributionType == 0 )
45 {
46 theIsoFlag = true;
47 }
48 else if(theAngularDistributionType==1)
49 {
50 theIsoFlag = false;
51 G4int nEnergy;
52 aDataFile >> nEnergy;
53 theCoefficients = new G4NeutronHPLegendreStore(nEnergy);
54 theCoefficients->InitInterpolation(aDataFile);
55 G4double temp, energy;
56 G4int tempdep, nLegendre;
57 G4int i, ii;
58 for (i=0; i<nEnergy; i++)
59 {
60 aDataFile >> temp >> energy >> tempdep >> nLegendre;
61 energy *=eV;
62 theCoefficients->Init(i, energy, nLegendre);
63 theCoefficients->SetTemperature(i, temp);
64 G4double coeff=0;
65 for(ii=0; ii<nLegendre; ii++)
66 {
67 aDataFile >> coeff;
68 theCoefficients->SetCoeff(i, ii+1, coeff);
69 }
70 }
71 }
72 else if (theAngularDistributionType==2)
73 {
74 theIsoFlag = false;
75 G4int nEnergy;
76 aDataFile >> nEnergy;
77 theProbArray = new G4NeutronHPPartial(nEnergy, nEnergy);
78 theProbArray->InitInterpolation(aDataFile);
79 G4double temp, energy;
80 G4int tempdep;
81 for(G4int i=0; i<nEnergy; i++)
82 {
83 aDataFile >> temp >> energy >> tempdep;
84 energy *= eV;
85 theProbArray->SetT(i, temp);
86 theProbArray->SetX(i, energy);
87 theProbArray->InitData(i, aDataFile);
88 }
89 }
90 else
91 {
92 theIsoFlag = false;
93 G4cout << "unknown distribution found for Angular"<<G4endl;
94 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
95 }
96}
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4NeutronHPLegendreStore::SetCoeff
void SetCoeff(G4int i, G4int l, G4double coeff)
Definition: G4NeutronHPLegendreStore.hh:59
G4NeutronHPLegendreStore::Init
void Init(G4int i, G4double e, G4int n)
Definition: G4NeutronHPLegendreStore.hh:52
G4NeutronHPLegendreStore::InitInterpolation
void InitInterpolation(std::ifstream &aDataFile)
Definition: G4NeutronHPLegendreStore.hh:79
G4NeutronHPLegendreStore::SetTemperature
void SetTemperature(G4int i, G4double temp)
Definition: G4NeutronHPLegendreStore.hh:58
G4NeutronHPPartial::InitData
void InitData(G4int i, std::ifstream &aDataFile, G4double unit=1.)
Definition: G4NeutronHPPartial.hh:84
G4NeutronHPPartial::InitInterpolation
void InitInterpolation(G4int i, std::ifstream &aDataFile)
Definition: G4NeutronHPPartial.hh:61
G4NeutronHPPartial::SetT
void SetT(G4int i, G4double x)
Definition: G4NeutronHPPartial.hh:109
G4NeutronHPPartial::SetX
void SetX(G4int i, G4double x)
Definition: G4NeutronHPPartial.hh:108

Referenced by G4FissionLibrary::Init(), G4NeutronHPFSFissionFS::Init(), G4NeutronHPInelasticCompFS::Init(), G4NeutronHPFissionBaseFS::Init(), and G4NeutronHPInelasticBaseFS::Init().

◆ SampleAndUpdate()

void G4NeutronHPAngular::SampleAndUpdate ( G4ReactionProduct &  aNeutron)

Definition at line 98 of file G4NeutronHPAngular.cc.

99{
100
101 //********************************************************************
102 //EMendoza -> sampling can be isotropic in LAB or in CMS
103 /*
104 if(theIsoFlag)
105 {
106// G4cout << "Angular result "<<aHadron.GetTotalMomentum()<<" ";
107// @@@ add code for isotropic emission in CMS.
108 G4double costheta = 2.*G4UniformRand()-1;
109 G4double theta = std::acos(costheta);
110 G4double phi = twopi*G4UniformRand();
111 G4double sinth = std::sin(theta);
112 G4double en = aHadron.GetTotalMomentum();
113 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
114 aHadron.SetMomentum( temp );
115 aHadron.Lorentz(aHadron, -1.*theTarget);
116 }
117 else
118 {
119 */
120 //********************************************************************
121 if(frameFlag == 1) // LAB
122 {
123 G4double en = aHadron.GetTotalMomentum();
124 G4ReactionProduct boosted;
125 boosted.Lorentz(theNeutron, theTarget);
126 G4double kineticEnergy = boosted.GetKineticEnergy();
127 G4double cosTh = 0.0;
128 //********************************************************************
129 //EMendoza --> sampling can be also isotropic
130 /*
131 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
132 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
133 */
134 //********************************************************************
135 if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
136 else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
137 else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
138 else{
139 G4cout << "unknown distribution found for Angular"<<G4endl;
140 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
141 }
142 //********************************************************************
143 G4double theta = std::acos(cosTh);
144 G4double phi = twopi*G4UniformRand();
145 G4double sinth = std::sin(theta);
146 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
147 aHadron.SetMomentum( temp );
148 }
149 else if(frameFlag == 2) // costh in CMS
150 {
151 G4ReactionProduct boostedN;
152 boostedN.Lorentz(theNeutron, theTarget);
153 G4double kineticEnergy = boostedN.GetKineticEnergy();
154
155 G4double cosTh = 0.0;
156 //********************************************************************
157 //EMendoza --> sampling can be also isotropic
158 /*
159 if(theAngularDistributionType == 1) cosTh = theCoefficients->SampleMax(kineticEnergy);
160 if(theAngularDistributionType == 2) cosTh = theProbArray->Sample(kineticEnergy);
161 */
162 //********************************************************************
163 if(theIsoFlag){cosTh =2.*G4UniformRand()-1;}
164 else if(theAngularDistributionType == 1) {cosTh = theCoefficients->SampleMax(kineticEnergy);}
165 else if(theAngularDistributionType == 2) {cosTh = theProbArray->Sample(kineticEnergy);}
166 else{
167 G4cout << "unknown distribution found for Angular"<<G4endl;
168 throw G4HadronicException(__FILE__, __LINE__, "unknown distribution needs implementation!!!");
169 }
170 //********************************************************************
171//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern)
172/*
173 if(theAngularDistributionType == 1) // LAB
174 {
175 G4double en = aHadron.GetTotalMomentum();
176 G4ReactionProduct boosted;
177 boosted.Lorentz(theNeutron, theTarget);
178 G4double kineticEnergy = boosted.GetKineticEnergy();
179 G4double cosTh = theCoefficients->SampleMax(kineticEnergy);
180 G4double theta = std::acos(cosTh);
181 G4double phi = twopi*G4UniformRand();
182 G4double sinth = std::sin(theta);
183 G4ThreeVector temp(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
184 aHadron.SetMomentum( temp );
185 }
186 else if(theAngularDistributionType == 2) // costh in CMS {
187 }
188*/
189
190// G4ReactionProduct boostedN;
191// boostedN.Lorentz(theNeutron, theTarget);
192// G4double kineticEnergy = boostedN.GetKineticEnergy();
193// G4double cosTh = theProbArray->Sample(kineticEnergy);
194
195 G4double theta = std::acos(cosTh);
196 G4double phi = twopi*G4UniformRand();
197 G4double sinth = std::sin(theta);
198
199 G4ThreeVector temp(sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) ); //CMS
200
201//080612TK bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #5
202/*
203 G4double en = aHadron.GetTotalEnergy(); // Target rest
204
205 // get trafo from Target rest frame to CMS
206 G4ReactionProduct boostedT;
207 boostedT.Lorentz(theTarget, theTarget);
208
209 G4ThreeVector the3Neutron = boostedN.GetMomentum();
210 G4double nEnergy = boostedN.GetTotalEnergy();
211 G4ThreeVector the3Target = boostedT.GetMomentum();
212 G4double tEnergy = boostedT.GetTotalEnergy();
213 G4double totE = nEnergy+tEnergy;
214 G4ThreeVector the3trafo = -the3Target-the3Neutron;
215 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
216 trafo.SetMomentum(the3trafo);
217 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
218 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
219 trafo.SetMass(sqrts);
220 trafo.SetTotalEnergy(totE);
221
222 G4double gamma = trafo.GetTotalEnergy()/trafo.GetMass();
223 G4double cosalpha = temp*trafo.GetMomentum()/trafo.GetTotalMomentum()/temp.mag();
224 G4double fac = cosalpha*trafo.GetTotalMomentum()/trafo.GetMass();
225 fac*=gamma;
226
227 G4double mom;
228// For G4FPE_DEBUG ON
229 G4double mom2 = ( en*fac*en*fac -
230 (fac*fac - gamma*gamma)*
231 (en*en - gamma*gamma*aHadron.GetMass()*aHadron.GetMass())
232 );
233 if ( mom2 > 0.0 )
234 mom = std::sqrt( mom2 );
235 else
236 mom = 0.0;
237
238 mom = -en*fac - mom;
239 mom /= (fac*fac-gamma*gamma);
240 temp = mom*temp;
241
242 aHadron.SetMomentum( temp ); // now all in CMS
243 aHadron.SetTotalEnergy( std::sqrt( mom*mom + aHadron.GetMass()*aHadron.GetMass() ) );
244 aHadron.Lorentz(aHadron, trafo); // now in target rest frame
245*/
246 // Determination of the hadron kinetic energy in CMS
247 // aHadron.GetKineticEnergy() is actually the residual kinetic energy in CMS (or target frame)
248 // kineticEnergy is incident neutron kinetic energy in CMS (or target frame)
249 G4double QValue = aHadron.GetKineticEnergy() - kineticEnergy;
250 G4double A1 = theTarget.GetMass()/boostedN.GetMass();
251 G4double A1prim = aHadron.GetMass()/ boostedN.GetMass();
252 G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*kineticEnergy+(1+A1)*QValue);
253 G4double totalE = kinE + aHadron.GetMass();
254 G4double mom2 = totalE*totalE - aHadron.GetMass()*aHadron.GetMass();
255 G4double mom;
256 if ( mom2 > 0.0 ) mom = std::sqrt( mom2 );
257 else mom = 0.0;
258
259 aHadron.SetMomentum( mom*temp ); // Set momentum in CMS
260 aHadron.SetKineticEnergy(kinE); // Set kinetic energy in CMS
261
262 // get trafo from Target rest frame to CMS
263 G4ReactionProduct boostedT;
264 boostedT.Lorentz(theTarget, theTarget);
265
266 G4ThreeVector the3Neutron = boostedN.GetMomentum();
267 G4double nEnergy = boostedN.GetTotalEnergy();
268 G4ThreeVector the3Target = boostedT.GetMomentum();
269 G4double tEnergy = boostedT.GetTotalEnergy();
270 G4double totE = nEnergy+tEnergy;
271 G4ThreeVector the3trafo = -the3Target-the3Neutron;
272 G4ReactionProduct trafo; // for transformation from CMS to target rest frame
273 trafo.SetMomentum(the3trafo);
274 G4double cmsMom = std::sqrt(the3trafo*the3trafo);
275 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
276 trafo.SetMass(sqrts);
277 trafo.SetTotalEnergy(totE);
278
279 aHadron.Lorentz(aHadron, trafo);
280
281 }
282 else
283 {
284 throw G4HadronicException(__FILE__, __LINE__, "Tried to sample non isotropic neutron angular");
285 }
286 aHadron.Lorentz(aHadron, -1.*theTarget);
287// G4cout << aHadron.GetMomentum()<<" ";
288// G4cout << aHadron.GetTotalMomentum()<<G4endl;
289}
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
G4NeutronHPLegendreStore::SampleMax
G4double SampleMax(G4double energy)
Definition: G4NeutronHPLegendreStore.cc:117
G4NeutronHPPartial::Sample
G4double Sample(G4double x)
Definition: G4NeutronHPPartial.cc:103
G4ReactionProduct::SetMomentum
void SetMomentum(const G4double x, const G4double y, const G4double z)
Definition: G4ReactionProduct.cc:166
G4ReactionProduct::SetTotalEnergy
void SetTotalEnergy(const G4double en)
Definition: G4ReactionProduct.hh:138
G4ReactionProduct::GetKineticEnergy
G4double GetKineticEnergy() const
Definition: G4ReactionProduct.hh:135
G4ReactionProduct::GetTotalEnergy
G4double GetTotalEnergy() const
Definition: G4ReactionProduct.hh:126
G4ReactionProduct::GetMomentum
G4ThreeVector GetMomentum() const
Definition: G4ReactionProduct.hh:120
G4ReactionProduct::Lorentz
void Lorentz(const G4ReactionProduct &p1, const G4ReactionProduct &p2)
Definition: G4ReactionProduct.cc:200
G4ReactionProduct::GetMass
G4double GetMass() const
Definition: G4ReactionProduct.hh:147
G4ReactionProduct::SetMass
void SetMass(const G4double mas)
Definition: G4ReactionProduct.hh:144

Referenced by G4NeutronHPFissionBaseFS::ApplyYourself(), G4NeutronHPFSFissionFS::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), and G4NeutronHPInelasticCompFS::CompositeApply().

◆ SetNeutron()

void G4NeutronHPAngular::SetNeutron ( const G4ReactionProduct &  aNeutron)
inline

Definition at line 70 of file G4NeutronHPAngular.hh.

70{ theNeutron = aNeutron; }

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFissionBaseFS::SetNeutron(), and G4NeutronHPFSFissionFS::SetNeutron().

◆ SetTarget()

void G4NeutronHPAngular::SetTarget ( const G4ReactionProduct &  aTarget)
inline

Definition at line 68 of file G4NeutronHPAngular.hh.

68{ theTarget = aTarget; }

Referenced by G4FissionLibrary::ApplyYourself(), G4NeutronHPInelasticBaseFS::BaseApply(), G4NeutronHPInelasticCompFS::InitDistributionInitialState(), G4NeutronHPFissionBaseFS::SetTarget(), and G4NeutronHPFSFissionFS::SetTarget().

The documentation for this class was generated from the following files: