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

#include <G4NeutronHPDiscreteTwoBody.hh>

+ Inheritance diagram for G4NeutronHPDiscreteTwoBody:

Public Member Functions

 G4NeutronHPDiscreteTwoBody ()
 
 ~G4NeutronHPDiscreteTwoBody ()
 
void Init (std::ifstream &aDataFile)
 
G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass)
 
G4double MeanEnergyOfThisInteraction ()
 
- Public Member Functions inherited from G4VNeutronHPEnergyAngular
 G4VNeutronHPEnergyAngular ()
 
virtual ~G4VNeutronHPEnergyAngular ()
 
virtual void Init (std::ifstream &aDataFile)=0
 
virtual G4ReactionProduct * Sample (G4double anEnergy, G4double massCode, G4double mass)=0
 
virtual G4double MeanEnergyOfThisInteraction ()=0
 
void SetNeutron (G4ReactionProduct *aNeutron)
 
void SetTarget (G4ReactionProduct *aTarget)
 
G4ReactionProduct * GetTarget ()
 
G4ReactionProduct * GetNeutron ()
 
G4ReactionProduct * GetCMS ()
 
void SetQValue (G4double aValue)
 
virtual void ClearHistories ()
 

Additional Inherited Members

- Protected Member Functions inherited from G4VNeutronHPEnergyAngular
G4double GetQValue ()
 

Detailed Description

Definition at line 42 of file G4NeutronHPDiscreteTwoBody.hh.

Constructor & Destructor Documentation

◆ G4NeutronHPDiscreteTwoBody()

G4NeutronHPDiscreteTwoBody::G4NeutronHPDiscreteTwoBody ( )
inline

Definition at line 46 of file G4NeutronHPDiscreteTwoBody.hh.

47 {
48 theCoeff = 0;
49 }

◆ ~G4NeutronHPDiscreteTwoBody()

G4NeutronHPDiscreteTwoBody::~G4NeutronHPDiscreteTwoBody ( )
inline

Definition at line 50 of file G4NeutronHPDiscreteTwoBody.hh.

51 {
52 if(theCoeff!=0) delete [] theCoeff;
53 }

Member Function Documentation

◆ Init()

void G4NeutronHPDiscreteTwoBody::Init ( std::ifstream &  aDataFile)
inlinevirtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 55 of file G4NeutronHPDiscreteTwoBody.hh.

56 {
57 aDataFile >> nEnergy;
58 theManager.Init(aDataFile);
59 theCoeff = new G4NeutronHPLegendreTable[nEnergy];
60 for(G4int i=0; i<nEnergy; i++)
61 {
62 G4double energy;
63 G4int aRep, nCoeff;
64 aDataFile >> energy >> aRep >> nCoeff;
65 energy*=CLHEP::eV;
66 G4int nPoints=nCoeff;
67 if(aRep>0) nPoints*=2;
68 //theCoeff[i].Init(energy, nPoints);
69
70 theCoeff[i].Init(energy, nPoints-1);
71 theCoeff[i].SetRepresentation(aRep);
72 for(G4int ii=0; ii<nPoints; ii++)
73 {
74 G4double y;
75 aDataFile >> y;
76 theCoeff[i].SetCoeff(ii, y);
77 }
78 }
79 }
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4InterpolationManager::Init
void Init(G4int aScheme, G4int aRange)
Definition: G4InterpolationManager.hh:85
G4NeutronHPLegendreTable::SetCoeff
void SetCoeff(G4int l, G4double coeff)
Definition: G4NeutronHPLegendreTable.hh:92
G4NeutronHPLegendreTable::Init
void Init(std::ifstream &aDataFile)
Definition: G4NeutronHPLegendreTable.hh:67

◆ MeanEnergyOfThisInteraction()

G4double G4NeutronHPDiscreteTwoBody::MeanEnergyOfThisInteraction ( )
inlinevirtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 82 of file G4NeutronHPDiscreteTwoBody.hh.

82{ return -1; }

◆ Sample()

G4ReactionProduct * G4NeutronHPDiscreteTwoBody::Sample ( G4double  anEnergy,
G4double  massCode,
G4double  mass 
)
virtual

Implements G4VNeutronHPEnergyAngular.

Definition at line 48 of file G4NeutronHPDiscreteTwoBody.cc.

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