| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 080612 Bug fix contribution from Benoit Pir
31 // 080709 Bug fix Sampling Legendre expansion
32 // 101110 Bug fix in MF=6, LAW=2 case; contrib
33 //
34 // P. Arce, June-2014 Conversion neutron_hp to
35 //
36 #include "G4ParticleHPDiscreteTwoBody.hh"
37
38 #include "G4Alpha.hh"
39 #include "G4Deuteron.hh"
40 #include "G4Electron.hh"
41 #include "G4Gamma.hh"
42 #include "G4He3.hh"
43 #include "G4Neutron.hh"
44 #include "G4ParticleHPLegendreStore.hh"
45 #include "G4ParticleHPVector.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4Positron.hh"
48 #include "G4Proton.hh"
49 #include "G4Triton.hh"
50
51 G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscr
52 {
53 if (G4ParticleHPManager::GetInstance()->GetP
54 bCheckDiffCoeffRepr = false;
55 }
56
57 G4ParticleHPDiscreteTwoBody::~G4ParticleHPDisc
58 {
59 delete[] theCoeff;
60 }
61
62 void G4ParticleHPDiscreteTwoBody::Init(std::is
63 {
64 aDataFile >> nEnergy;
65 theManager.Init(aDataFile);
66 const std::size_t tsize = nEnergy > 0 ? nEne
67 theCoeff = new G4ParticleHPLegendreTable[tsi
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 << " G4Parti
73 //<< " " << aRep << " " << nCoeff << G4end
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 }
86
87 G4ReactionProduct* G4ParticleHPDiscreteTwoBody
88
89 { // Interpolation still only for the most us
90 auto result = new G4ReactionProduct;
91 auto Z = static_cast<G4int>(massCode / 1000)
92 auto A = static_cast<G4int>(massCode - 1000
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(G4Positr
100 }
101 else if (A == 1) {
102 result->SetDefinition(G4Neutron::Neutron()
103 if (Z == 1) result->SetDefinition(G4Proton
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::H
111 }
112 else if (A == 4) {
113 result->SetDefinition(G4Alpha::Alpha());
114 if (Z != 2) throw G4HadronicException(__FI
115 }
116 else {
117 throw G4HadronicException(__FILE__, __LINE
118 "G4ParticleHPDis
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) br
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 Pir
136 cosTh = theStore.SampleDiscreteTwoBody(a
137 }
138 else if (theCoeff[it].GetRepresentation()
139 {
140 G4ParticleHPVector theStore;
141 G4InterpolationManager aManager;
142 aManager.Init(LINLIN, theCoeff[it].GetNu
143 theStore.SetInterpolationManager(aManage
144 for (i = 0; i < theCoeff[it].GetNumberOf
145 // 101110
146 // theStore.SetX(i, theCoeff[it].GetCo
147 // theStore.SetY(i, theCoeff[it].GetCo
148 theStore.SetX(i / 2, theCoeff[it].GetC
149 theStore.SetY(i / 2, theCoeff[it].GetC
150 }
151 cosTh = theStore.Sample();
152 }
153 else if (theCoeff[it].GetRepresentation()
154 {
155 G4ParticleHPVector theStore;
156 G4InterpolationManager aManager;
157 aManager.Init(LOGLIN, theCoeff[it].GetNu
158 theStore.SetInterpolationManager(aManage
159 for (i = 0; i < theCoeff[it].GetNumberOf
160 // 101110
161 // theStore.SetX(i, theCoeff[it].GetCo
162 // theStore.SetY(i, theCoeff[it].GetCo
163 theStore.SetX(i / 2, theCoeff[it].GetC
164 theStore.SetY(i / 2, theCoeff[it].GetC
165 }
166 cosTh = theStore.Sample();
167 }
168 else {
169 throw G4HadronicException(__FILE__, __LI
170 "unknown repre
171 }
172 }
173 else {
174 if (!bCheckDiffCoeffRepr
175 || theCoeff[it].GetRepresentation() ==
176 {
177 if (theCoeff[it].GetRepresentation() ==
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)
184 theStore.SetManager(aManager);
185 // 080709 TKDB
186 cosTh = theStore.SampleDiscreteTwoBody
187 }
188 else if (theCoeff[it].GetRepresentation(
189 {
190 G4ParticleHPVector theBuff1;
191 G4InterpolationManager aManager1;
192 aManager1.Init(LINLIN, theCoeff[it - 1
193 theBuff1.SetInterpolationManager(aMana
194 for (i = 0; i < theCoeff[it - 1].GetNu
195 // 101110
196 theBuff1.SetX(i / 2, theCoeff[it - 1
197 theBuff1.SetY(i / 2, theCoeff[it - 1
198 }
199 G4ParticleHPVector theBuff2;
200 G4InterpolationManager aManager2;
201 aManager2.Init(LINLIN, theCoeff[it].Ge
202 theBuff2.SetInterpolationManager(aMana
203 for (i = 0; i < theCoeff[it].GetNumber
204 theBuff2.SetX(i / 2, theCoeff[it].Ge
205 theBuff2.SetY(i / 2, theCoeff[it].Ge
206 }
207
208 G4double x1 = theCoeff[it - 1].GetEner
209 G4double x2 = theCoeff[it].GetEnergy()
210 G4double x = anEnergy;
211 G4double y1, y2, y, mu;
212
213 G4ParticleHPVector theStore1;
214 theStore1.SetInterpolationManager(aMan
215 G4ParticleHPVector theStore2;
216 theStore2.SetInterpolationManager(aMan
217 G4ParticleHPVector theStore;
218
219 // for fixed mu get p1, p2 and interpo
220 for (i = 0; i < theBuff1.GetVectorLeng
221 mu = theBuff1.GetX(i);
222 y1 = theBuff1.GetY(i);
223 y2 = theBuff2.GetY(mu);
224 y = theInt.Interpolate(theManager.Ge
225 theStore1.SetData(i, mu, y);
226 }
227 for (i = 0; i < theBuff2.GetVectorLeng
228 mu = theBuff2.GetX(i);
229 y1 = theBuff2.GetY(i);
230 y2 = theBuff1.GetY(mu);
231 y = theInt.Interpolate(theManager.Ge
232 theStore2.SetData(i, mu, y);
233 }
234 theStore.Merge(&theStore1, &theStore2)
235 cosTh = theStore.Sample();
236 }
237 else if (theCoeff[it].GetRepresentation(
238 {
239 G4ParticleHPVector theBuff1;
240 G4InterpolationManager aManager1;
241 aManager1.Init(LOGLIN, theCoeff[it - 1
242 theBuff1.SetInterpolationManager(aMana
243 for (i = 0; i < theCoeff[it - 1].GetNu
244 // 101110
245 theBuff1.SetX(i / 2, theCoeff[it - 1
246 theBuff1.SetY(i / 2, theCoeff[it - 1
247 }
248
249 G4ParticleHPVector theBuff2;
250 G4InterpolationManager aManager2;
251 aManager2.Init(LOGLIN, theCoeff[it].Ge
252 theBuff2.SetInterpolationManager(aMana
253 for (i = 0; i < theCoeff[it].GetNumber
254 // 101110
255 theBuff2.SetX(i / 2, theCoeff[it].Ge
256 theBuff2.SetY(i / 2, theCoeff[it].Ge
257 }
258
259 G4double x1 = theCoeff[it - 1].GetEner
260 G4double x2 = theCoeff[it].GetEnergy()
261 G4double x = anEnergy;
262 G4double y1, y2, y, mu;
263
264 G4ParticleHPVector theStore1;
265 theStore1.SetInterpolationManager(aMan
266 G4ParticleHPVector theStore2;
267 theStore2.SetInterpolationManager(aMan
268 G4ParticleHPVector theStore;
269
270 // for fixed mu get p1, p2 and interpo
271 for (i = 0; i < theBuff1.GetVectorLeng
272 mu = theBuff1.GetX(i);
273 y1 = theBuff1.GetY(i);
274 y2 = theBuff2.GetY(mu);
275 y = theInt.Interpolate(theManager.Ge
276 theStore1.SetData(i, mu, y);
277 }
278 for (i = 0; i < theBuff2.GetVectorLeng
279 mu = theBuff2.GetX(i);
280 y1 = theBuff2.GetY(i);
281 y2 = theBuff1.GetY(mu);
282 y = theInt.Interpolate(theManager.Ge
283 theStore2.SetData(i, mu, y);
284 }
285 theStore.Merge(&theStore1, &theStore2)
286 cosTh = theStore.Sample();
287 }
288 else {
289 throw G4HadronicException(__FILE__, __
290 "Two neighbo
291 }
292 }
293 else {
294 G4cout << " theCoeff[it].GetRepresent ME
295 << &theCoeff[it - 1] << G4endl;
296 G4cout << " theCoeff[it].GetRepresent "
297 << " != " << theCoeff[it - 1].Get
298
299 throw G4HadronicException(__FILE__, __LI
300 "unknown repre
301 }
302 }
303
304 // now get the energy from kinematics and Q-
305
306 // G4double restEnergy = anEnergy+GetQValue(
307
308 // assumed to be in CMS @@@@@@@@@@@@@@@@@
309
310 // 080612TK contribution from Benoit Pirard
311 // G4double residualMass = GetTarget()->Ge
312 // - result->GetMass
313 // G4double kinE = restEnergy/(1+result->Get
314 G4double A1 = GetTarget()->GetMass() / GetPr
315 G4double A1prim = result->GetMass() / GetPro
316 // G4double E1 = (A1+1)*(A1+1)/A1/A1*an
317 // Bug fix Bugzilla #1815
318 G4double E1 = anEnergy;
319 G4double kinE = (A1 + 1 - A1prim) / (A1 + 1)
320
321 result->SetKineticEnergy(kinE); // non rela
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.
326 G4double mtot = result->GetTotalMomentum();
327 G4ThreeVector tempVector(mtot * sinth * std:
328 mtot * cosTh);
329 result->SetMomentum(tempVector);
330 return result;
331 }
332