Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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 55 of file G4ParticleHPAngular.hh.

56 {
57 theIsoFlag = true;
58 theCoefficients = nullptr;
59 theProbArray = nullptr;
60 toBeCached val;
61 fCache.Put(val);
62
63 theAngularDistributionType = 0;
64 frameFlag = 0;
65 targetMass = 0.0;
66 }

◆ ~G4ParticleHPAngular()

G4ParticleHPAngular::~G4ParticleHPAngular ( )
inline

Definition at line 68 of file G4ParticleHPAngular.hh.

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

Member Function Documentation

◆ GetTargetMass()

G4double G4ParticleHPAngular::GetTargetMass ( )
inline

Definition at line 85 of file G4ParticleHPAngular.hh.

85{ return targetMass; }

◆ Init()

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

Definition at line 42 of file G4ParticleHPAngular.cc.

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

◆ SampleAndUpdate()

void G4ParticleHPAngular::SampleAndUpdate ( G4ReactionProduct & anIncidentParticle)

Definition at line 94 of file G4ParticleHPAngular.cc.

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

◆ SetProjectileRP()

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

Definition at line 80 of file G4ParticleHPAngular.hh.

81 {
82 fCache.Get().theProjectileRP = &anIncidentParticleRP;
83 }

◆ SetTarget()

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

Definition at line 78 of file G4ParticleHPAngular.hh.

78{ fCache.Get().theTarget = &aTarget; }
The documentation for this class was generated from the following files: