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
G4ParticleHPDiscreteTwoBody Class Reference

#include <G4ParticleHPDiscreteTwoBody.hh>

+ Inheritance diagram for G4ParticleHPDiscreteTwoBody:

Public Member Functions

 G4ParticleHPDiscreteTwoBody ()
 
 ~G4ParticleHPDiscreteTwoBody () override
 
void Init (std::istream &aDataFile) override
 
G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass) override
 
G4double MeanEnergyOfThisInteraction () override
 
- Public Member Functions inherited from G4VParticleHPEnergyAngular
 G4VParticleHPEnergyAngular ()
 
virtual ~G4VParticleHPEnergyAngular ()=default
 
void SetProjectileRP (G4ReactionProduct *aIncidentParticleRP)
 
void SetTarget (G4ReactionProduct *aTarget)
 
G4ReactionProduct * GetTarget () const
 
G4ReactionProduct * GetProjectileRP () const
 
G4ReactionProduct * GetCMS ()
 
void SetQValue (G4double aValue)
 
virtual void ClearHistories ()
 
 G4VParticleHPEnergyAngular (G4VParticleHPEnergyAngular &)=delete
 
G4VParticleHPEnergyAngular & operator= (const G4VParticleHPEnergyAngular &right)=delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4VParticleHPEnergyAngular
G4double GetQValue () const
 

Detailed Description

Definition at line 43 of file G4ParticleHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

◆ G4ParticleHPDiscreteTwoBody()

G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscreteTwoBody ( )

Definition at line 51 of file G4ParticleHPDiscreteTwoBody.cc.

52{
53 if (G4ParticleHPManager::GetInstance()->GetPHPCheck())
54 bCheckDiffCoeffRepr = false;
55}
G4ParticleHPManager::GetInstance
static G4ParticleHPManager * GetInstance()
Definition G4ParticleHPManager.cc:97

◆ ~G4ParticleHPDiscreteTwoBody()

G4ParticleHPDiscreteTwoBody::~G4ParticleHPDiscreteTwoBody ( )
override

Definition at line 57 of file G4ParticleHPDiscreteTwoBody.cc.

58{
59 delete[] theCoeff;
60}

Member Function Documentation

◆ Init()

void G4ParticleHPDiscreteTwoBody::Init ( std::istream & aDataFile)
overridevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 62 of file G4ParticleHPDiscreteTwoBody.cc.

63{
64 aDataFile >> nEnergy;
65 theManager.Init(aDataFile);
66 const std::size_t tsize = nEnergy > 0 ? nEnergy : 1;
67 theCoeff = new G4ParticleHPLegendreTable[tsize];
68 for (std::size_t i = 0; i < tsize; ++i) {
69 G4double energy;
70 G4int aRep, nCoeff;
71 aDataFile >> energy >> aRep >> nCoeff;
72 // G4cout << this << " " << i << " G4ParticleHPDiscreteTwoBody READ DATA " << energy
73 //<< " " << aRep << " " << nCoeff << G4endl;
74 energy *= CLHEP::eV;
75 G4int nPoints = nCoeff;
76 if (aRep > 0) nPoints *= 2;
77 theCoeff[i].Init(energy, nPoints - 1);
78 theCoeff[i].SetRepresentation(aRep);
79 for (G4int ii = 0; ii < nPoints; ++ii) {
80 G4double y;
81 aDataFile >> y;
82 theCoeff[i].SetCoeff(ii, y);
83 }
84 }
85}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4INCL::KinematicsUtils::energy
G4double energy(const ThreeVector &p, const G4double m)
Definition G4INCLKinematicsUtils.cc:185

◆ MeanEnergyOfThisInteraction()

G4double G4ParticleHPDiscreteTwoBody::MeanEnergyOfThisInteraction ( )
inlineoverridevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 52 of file G4ParticleHPDiscreteTwoBody.hh.

52{ return -1.0; }

◆ Sample()

G4ReactionProduct * G4ParticleHPDiscreteTwoBody::Sample ( G4double anEnergy,
G4double massCode,
G4double mass )
overridevirtual

Implements G4VParticleHPEnergyAngular.

Definition at line 87 of file G4ParticleHPDiscreteTwoBody.cc.

89{ // Interpolation still only for the most used parts; rest to be Done @@@@@
90 auto result = new G4ReactionProduct;
91 auto Z = static_cast<G4int>(massCode / 1000);
92 auto A = static_cast<G4int>(massCode - 1000 * Z);
93
94 if (massCode == 0) {
95 result->SetDefinition(G4Gamma::Gamma());
96 }
97 else if (A == 0) {
98 result->SetDefinition(G4Electron::Electron());
99 if (Z == 1) result->SetDefinition(G4Positron::Positron());
100 }
101 else if (A == 1) {
102 result->SetDefinition(G4Neutron::Neutron());
103 if (Z == 1) result->SetDefinition(G4Proton::Proton());
104 }
105 else if (A == 2) {
106 result->SetDefinition(G4Deuteron::Deuteron());
107 }
108 else if (A == 3) {
109 result->SetDefinition(G4Triton::Triton());
110 if (Z == 2) result->SetDefinition(G4He3::He3());
111 }
112 else if (A == 4) {
113 result->SetDefinition(G4Alpha::Alpha());
114 if (Z != 2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
115 }
116 else {
117 throw G4HadronicException(__FILE__, __LINE__,
118 "G4ParticleHPDiscreteTwoBody: Unknown ion case 2");
119 }
120
121 // get cosine(theta)
122 G4int i(0), it(0);
123 G4double cosTh(0);
124 for (i = 0; i < nEnergy; i++) {
125 it = i;
126 if (theCoeff[i].GetEnergy() > anEnergy) break;
127 }
128 if (it == 0 || it == nEnergy - 1) {
129 if (theCoeff[it].GetRepresentation() == 0) {
130 // TK Legendre expansion
131 G4ParticleHPLegendreStore theStore(1);
132 theStore.SetCoeff(0, theCoeff);
133 theStore.SetManager(theManager);
134 // cosTh = theStore.SampleMax(anEnergy);
135 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
136 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
137 }
138 else if (theCoeff[it].GetRepresentation() == 12) // means LINLIN
139 {
140 G4ParticleHPVector theStore;
141 G4InterpolationManager aManager;
142 aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly() / 2);
143 theStore.SetInterpolationManager(aManager);
144 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
145 // 101110
146 // theStore.SetX(i, theCoeff[it].GetCoeff(i));
147 // theStore.SetY(i, theCoeff[it].GetCoeff(i));
148 theStore.SetX(i / 2, theCoeff[it].GetCoeff(i));
149 theStore.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
150 }
151 cosTh = theStore.Sample();
152 }
153 else if (theCoeff[it].GetRepresentation() == 14) // this is LOGLIN
154 {
155 G4ParticleHPVector theStore;
156 G4InterpolationManager aManager;
157 aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly() / 2);
158 theStore.SetInterpolationManager(aManager);
159 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
160 // 101110
161 // theStore.SetX(i, theCoeff[it].GetCoeff(i));
162 // theStore.SetY(i, theCoeff[it].GetCoeff(i));
163 theStore.SetX(i / 2, theCoeff[it].GetCoeff(i));
164 theStore.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
165 }
166 cosTh = theStore.Sample();
167 }
168 else {
169 throw G4HadronicException(__FILE__, __LINE__,
170 "unknown representation type in Two-body scattering");
171 }
172 }
173 else {
174 if (!bCheckDiffCoeffRepr
175 || theCoeff[it].GetRepresentation() == theCoeff[it - 1].GetRepresentation())
176 {
177 if (theCoeff[it].GetRepresentation() == 0) {
178 // TK Legendre expansion
179 G4ParticleHPLegendreStore theStore(2);
180 theStore.SetCoeff(0, &(theCoeff[it - 1]));
181 theStore.SetCoeff(1, &(theCoeff[it]));
182 G4InterpolationManager aManager;
183 aManager.Init(theManager.GetScheme(it), 2);
184 theStore.SetManager(aManager);
185 // 080709 TKDB
186 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
187 }
188 else if (theCoeff[it].GetRepresentation() == 12) // LINLIN
189 {
190 G4ParticleHPVector theBuff1;
191 G4InterpolationManager aManager1;
192 aManager1.Init(LINLIN, theCoeff[it - 1].GetNumberOfPoly() / 2);
193 theBuff1.SetInterpolationManager(aManager1);
194 for (i = 0; i < theCoeff[it - 1].GetNumberOfPoly(); i += 2) {
195 // 101110
196 theBuff1.SetX(i / 2, theCoeff[it - 1].GetCoeff(i));
197 theBuff1.SetY(i / 2, theCoeff[it - 1].GetCoeff(i + 1));
198 }
199 G4ParticleHPVector theBuff2;
200 G4InterpolationManager aManager2;
201 aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly() / 2);
202 theBuff2.SetInterpolationManager(aManager2);
203 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
204 theBuff2.SetX(i / 2, theCoeff[it].GetCoeff(i));
205 theBuff2.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
206 }
207
208 G4double x1 = theCoeff[it - 1].GetEnergy();
209 G4double x2 = theCoeff[it].GetEnergy();
210 G4double x = anEnergy;
211 G4double y1, y2, y, mu;
212
213 G4ParticleHPVector theStore1;
214 theStore1.SetInterpolationManager(aManager1);
215 G4ParticleHPVector theStore2;
216 theStore2.SetInterpolationManager(aManager2);
217 G4ParticleHPVector theStore;
218
219 // for fixed mu get p1, p2 and interpolate according to x
220 for (i = 0; i < theBuff1.GetVectorLength(); ++i) {
221 mu = theBuff1.GetX(i);
222 y1 = theBuff1.GetY(i);
223 y2 = theBuff2.GetY(mu);
224 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
225 theStore1.SetData(i, mu, y);
226 }
227 for (i = 0; i < theBuff2.GetVectorLength(); ++i) {
228 mu = theBuff2.GetX(i);
229 y1 = theBuff2.GetY(i);
230 y2 = theBuff1.GetY(mu);
231 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
232 theStore2.SetData(i, mu, y);
233 }
234 theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
235 cosTh = theStore.Sample();
236 }
237 else if (theCoeff[it].GetRepresentation() == 14) // TK LOG_LIN
238 {
239 G4ParticleHPVector theBuff1;
240 G4InterpolationManager aManager1;
241 aManager1.Init(LOGLIN, theCoeff[it - 1].GetNumberOfPoly() / 2);
242 theBuff1.SetInterpolationManager(aManager1);
243 for (i = 0; i < theCoeff[it - 1].GetNumberOfPoly(); i += 2) {
244 // 101110
245 theBuff1.SetX(i / 2, theCoeff[it - 1].GetCoeff(i));
246 theBuff1.SetY(i / 2, theCoeff[it - 1].GetCoeff(i + 1));
247 }
248
249 G4ParticleHPVector theBuff2;
250 G4InterpolationManager aManager2;
251 aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly() / 2);
252 theBuff2.SetInterpolationManager(aManager2);
253 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
254 // 101110
255 theBuff2.SetX(i / 2, theCoeff[it].GetCoeff(i));
256 theBuff2.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
257 }
258
259 G4double x1 = theCoeff[it - 1].GetEnergy();
260 G4double x2 = theCoeff[it].GetEnergy();
261 G4double x = anEnergy;
262 G4double y1, y2, y, mu;
263
264 G4ParticleHPVector theStore1;
265 theStore1.SetInterpolationManager(aManager1);
266 G4ParticleHPVector theStore2;
267 theStore2.SetInterpolationManager(aManager2);
268 G4ParticleHPVector theStore;
269
270 // for fixed mu get p1, p2 and interpolate according to x
271 for (i = 0; i < theBuff1.GetVectorLength(); i++) {
272 mu = theBuff1.GetX(i);
273 y1 = theBuff1.GetY(i);
274 y2 = theBuff2.GetY(mu);
275 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
276 theStore1.SetData(i, mu, y);
277 }
278 for (i = 0; i < theBuff2.GetVectorLength(); i++) {
279 mu = theBuff2.GetX(i);
280 y1 = theBuff2.GetY(i);
281 y2 = theBuff1.GetY(mu);
282 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
283 theStore2.SetData(i, mu, y);
284 }
285 theStore.Merge(&theStore1, &theStore2);
286 cosTh = theStore.Sample();
287 }
288 else {
289 throw G4HadronicException(__FILE__, __LINE__,
290 "Two neighbouring distributions with different interpolation");
291 }
292 }
293 else {
294 G4cout << " theCoeff[it].GetRepresent MEM " << it << " " << &theCoeff[it] << " "
295 << &theCoeff[it - 1] << G4endl;
296 G4cout << " theCoeff[it].GetRepresent " << it << " " << theCoeff[it].GetRepresentation()
297 << " != " << theCoeff[it - 1].GetRepresentation() << G4endl;
298
299 throw G4HadronicException(__FILE__, __LINE__,
300 "unknown representation type in Two-body scattering, case 2");
301 }
302 }
303
304 // now get the energy from kinematics and Q-value.
305
306 // G4double restEnergy = anEnergy+GetQValue();
307
308 // assumed to be in CMS @@@@@@@@@@@@@@@@@
309
310 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
311 // G4double residualMass = GetTarget()->GetMass() + GetNeutron()->GetMass()
312 // - result->GetMass() - GetQValue();
313 // G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
314 G4double A1 = GetTarget()->GetMass() / GetProjectileRP()->GetMass();
315 G4double A1prim = result->GetMass() / GetProjectileRP()->GetMass();
316 // G4double E1 = (A1+1)*(A1+1)/A1/A1*anEnergy;
317 // Bug fix Bugzilla #1815
318 G4double E1 = anEnergy;
319 G4double kinE = (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1) * (A1 * E1 + (1 + A1) * GetQValue());
320
321 result->SetKineticEnergy(kinE); // non relativistic @@
322 if (cosTh > 1.0) { cosTh = 1.0; }
323 else if(cosTh < -1.0) { cosTh = -1.0; }
324 G4double phi = CLHEP::twopi * G4UniformRand();
325 G4double sinth = std::sqrt((1.0 + cosTh)*(1.0 - cosTh));
326 G4double mtot = result->GetTotalMomentum();
327 G4ThreeVector tempVector(mtot * sinth * std::cos(phi), mtot * sinth * std::sin(phi),
328 mtot * cosTh);
329 result->SetMomentum(tempVector);
330 return result;
331}
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
A
const G4double A[17]
Definition G4WaterStopping.cc:53
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
G4Alpha::Alpha
static G4Alpha * Alpha()
Definition G4Alpha.cc:83
G4Deuteron::Deuteron
static G4Deuteron * Deuteron()
Definition G4Deuteron.cc:90
G4Electron::Electron
static G4Electron * Electron()
Definition G4Electron.cc:91
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition G4Gamma.cc:81
G4He3::He3
static G4He3 * He3()
Definition G4He3.cc:90
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition G4InterpolationManager.hh:82
G4Neutron::Neutron
static G4Neutron * Neutron()
Definition G4Neutron.cc:101
G4ParticleHPVector::SetY
void SetY(G4int i, G4double x)
Definition G4ParticleHPVector.hh:106
G4ParticleHPVector::SetX
void SetX(G4int i, G4double e)
Definition G4ParticleHPVector.hh:96
G4ParticleHPVector::SetData
void SetData(G4int i, G4double x, G4double y)
Definition G4ParticleHPVector.hh:89
G4ParticleHPVector::SetInterpolationManager
void SetInterpolationManager(const G4InterpolationManager &aManager)
Definition G4ParticleHPVector.hh:449
G4ParticleHPVector::GetY
G4double GetY(G4double x)
Definition G4ParticleHPVector.hh:184
G4ParticleHPVector::GetX
G4double GetX(G4int i) const
Definition G4ParticleHPVector.hh:120
G4ParticleHPVector::Merge
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)
Definition G4ParticleHPVector.hh:248
G4ParticleHPVector::GetVectorLength
G4int GetVectorLength() const
Definition G4ParticleHPVector.hh:185
G4Positron::Positron
static G4Positron * Positron()
Definition G4Positron.cc:90
G4Proton::Proton
static G4Proton * Proton()
Definition G4Proton.cc:90
G4ReactionProduct::GetMass
G4double GetMass() const
Definition G4ReactionProduct.hh:150
G4Triton::Triton
static G4Triton * Triton()
Definition G4Triton.cc:90
G4VParticleHPEnergyAngular::GetTarget
G4ReactionProduct * GetTarget() const
Definition G4VParticleHPEnergyAngular.hh:72
G4VParticleHPEnergyAngular::GetProjectileRP
G4ReactionProduct * GetProjectileRP() const
Definition G4VParticleHPEnergyAngular.hh:74
The documentation for this class was generated from the following files: