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

#include <G4ParticleHPAngular.hh>

Public Member Functions

 G4ParticleHPAngular ()
 
 ~G4ParticleHPAngular ()
 
void Init (std::istream &aDataFile)
 
void SampleAndUpdate (G4ReactionProduct &anIncidentParticle)
 
void SetTarget (const G4ReactionProduct &aTarget)
 
void SetProjectileRP (const G4ReactionProduct &anIncidentParticleRP)
 
G4double GetTargetMass ()
 

Detailed Description

Definition at line 45 of file G4ParticleHPAngular.hh.

Constructor & Destructor Documentation

◆ G4ParticleHPAngular()

G4ParticleHPAngular::G4ParticleHPAngular ( )
inline

Definition at line 56 of file G4ParticleHPAngular.hh.

57 {
58 theIsoFlag = true;
59 theCoefficients = 0;
60 theProbArray = 0;
61 toBeCached val;
62 fCache.Put( val );
63
64 theAngularDistributionType = 0;
65 frameFlag = 0;
66 targetMass = 0.0;
67 }
G4Cache::Put
void Put(const value_type &val) const
Definition: G4Cache.hh:321

◆ ~G4ParticleHPAngular()

G4ParticleHPAngular::~G4ParticleHPAngular ( )
inline

Definition at line 69 of file G4ParticleHPAngular.hh.

70 {
71 delete theCoefficients;
72 delete theProbArray;
73 }

Member Function Documentation

◆ GetTargetMass()

G4double G4ParticleHPAngular::GetTargetMass ( )
inline

Definition at line 89 of file G4ParticleHPAngular.hh.

90 {
91 return targetMass;
92 }

Referenced by G4ParticleHPInelasticBaseFS::BaseApply().

◆ Init()

void G4ParticleHPAngular::Init ( std::istream &  aDataFile)

Definition at line 41 of file G4ParticleHPAngular.cc.

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

Referenced by G4FissionLibrary::Init(), G4ParticleHPFSFissionFS::Init(), G4ParticleHPInelasticCompFS::Init(), G4ParticleHPFissionBaseFS::Init(), and G4ParticleHPInelasticBaseFS::Init().

◆ SampleAndUpdate()

void G4ParticleHPAngular::SampleAndUpdate ( G4ReactionProduct &  anIncidentParticle)

Definition at line 100 of file G4ParticleHPAngular.cc.

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

Referenced by G4ParticleHPFissionBaseFS::ApplyYourself(), G4ParticleHPFSFissionFS::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), and G4ParticleHPInelasticCompFS::CompositeApply().

◆ SetProjectileRP()

void G4ParticleHPAngular::SetProjectileRP ( const G4ReactionProduct &  anIncidentParticleRP)
inline

Definition at line 84 of file G4ParticleHPAngular.hh.

85 {
86 fCache.Get().theProjectileRP = &anIncidentParticleRP;
87 }

Referenced by G4FissionLibrary::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), G4ParticleHPInelasticCompFS::InitDistributionInitialState(), G4ParticleHPFissionBaseFS::SetNeutronRP(), and G4ParticleHPFSFissionFS::SetNeutronRP().

◆ SetTarget()

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

Definition at line 79 of file G4ParticleHPAngular.hh.

80 {
81 fCache.Get().theTarget = &aTarget;
82 }

Referenced by G4FissionLibrary::ApplyYourself(), G4ParticleHPInelasticBaseFS::BaseApply(), G4ParticleHPInelasticCompFS::InitDistributionInitialState(), G4ParticleHPFissionBaseFS::SetTarget(), and G4ParticleHPFSFissionFS::SetTarget().

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