| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // particle_hp -- source file 26 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 080612 Bug fix contribution from Benoit Pir 30 // 080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
31 // 080709 Bug fix Sampling Legendre expansion 31 // 080709 Bug fix Sampling Legendre expansion by T. Koi
32 // 101110 Bug fix in MF=6, LAW=2 case; contrib 32 // 101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
33 // 33 //
34 // P. Arce, June-2014 Conversion neutron_hp to 34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
35 // 35 //
36 #include "G4ParticleHPDiscreteTwoBody.hh" 36 #include "G4ParticleHPDiscreteTwoBody.hh"
37 37
38 #include "G4Alpha.hh" 38 #include "G4Alpha.hh"
39 #include "G4Deuteron.hh" 39 #include "G4Deuteron.hh"
40 #include "G4Electron.hh" 40 #include "G4Electron.hh"
41 #include "G4Gamma.hh" 41 #include "G4Gamma.hh"
42 #include "G4He3.hh" 42 #include "G4He3.hh"
43 #include "G4Neutron.hh" 43 #include "G4Neutron.hh"
44 #include "G4ParticleHPLegendreStore.hh" 44 #include "G4ParticleHPLegendreStore.hh"
45 #include "G4ParticleHPVector.hh" 45 #include "G4ParticleHPVector.hh"
46 #include "G4ParticleHPManager.hh" 46 #include "G4ParticleHPManager.hh"
47 #include "G4Positron.hh" 47 #include "G4Positron.hh"
48 #include "G4Proton.hh" 48 #include "G4Proton.hh"
49 #include "G4Triton.hh" 49 #include "G4Triton.hh"
50 50
51 G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscr 51 G4ParticleHPDiscreteTwoBody::G4ParticleHPDiscreteTwoBody()
52 { 52 {
53 if (G4ParticleHPManager::GetInstance()->GetP 53 if (G4ParticleHPManager::GetInstance()->GetPHPCheck())
54 bCheckDiffCoeffRepr = false; 54 bCheckDiffCoeffRepr = false;
55 } 55 }
56 56
57 G4ParticleHPDiscreteTwoBody::~G4ParticleHPDisc 57 G4ParticleHPDiscreteTwoBody::~G4ParticleHPDiscreteTwoBody()
58 { 58 {
59 delete[] theCoeff; 59 delete[] theCoeff;
60 } 60 }
61 61
62 void G4ParticleHPDiscreteTwoBody::Init(std::is 62 void G4ParticleHPDiscreteTwoBody::Init(std::istream& aDataFile)
63 { 63 {
64 aDataFile >> nEnergy; 64 aDataFile >> nEnergy;
65 theManager.Init(aDataFile); 65 theManager.Init(aDataFile);
66 const std::size_t tsize = nEnergy > 0 ? nEne 66 const std::size_t tsize = nEnergy > 0 ? nEnergy : 1;
67 theCoeff = new G4ParticleHPLegendreTable[tsi 67 theCoeff = new G4ParticleHPLegendreTable[tsize];
68 for (std::size_t i = 0; i < tsize; ++i) { 68 for (std::size_t i = 0; i < tsize; ++i) {
69 G4double energy; 69 G4double energy;
70 G4int aRep, nCoeff; 70 G4int aRep, nCoeff;
71 aDataFile >> energy >> aRep >> nCoeff; 71 aDataFile >> energy >> aRep >> nCoeff;
72 // G4cout << this << " " << i << " G4Parti 72 // G4cout << this << " " << i << " G4ParticleHPDiscreteTwoBody READ DATA " << energy
73 //<< " " << aRep << " " << nCoeff << G4end 73 //<< " " << aRep << " " << nCoeff << G4endl;
74 energy *= CLHEP::eV; 74 energy *= CLHEP::eV;
75 G4int nPoints = nCoeff; 75 G4int nPoints = nCoeff;
76 if (aRep > 0) nPoints *= 2; 76 if (aRep > 0) nPoints *= 2;
77 theCoeff[i].Init(energy, nPoints - 1); 77 theCoeff[i].Init(energy, nPoints - 1);
78 theCoeff[i].SetRepresentation(aRep); 78 theCoeff[i].SetRepresentation(aRep);
79 for (G4int ii = 0; ii < nPoints; ++ii) { 79 for (G4int ii = 0; ii < nPoints; ++ii) {
80 G4double y; 80 G4double y;
81 aDataFile >> y; 81 aDataFile >> y;
82 theCoeff[i].SetCoeff(ii, y); 82 theCoeff[i].SetCoeff(ii, y);
83 } 83 }
84 } 84 }
85 } 85 }
86 86
87 G4ReactionProduct* G4ParticleHPDiscreteTwoBody 87 G4ReactionProduct* G4ParticleHPDiscreteTwoBody::Sample(G4double anEnergy, G4double massCode,
88 88 G4double)
89 { // Interpolation still only for the most us 89 { // Interpolation still only for the most used parts; rest to be Done @@@@@
90 auto result = new G4ReactionProduct; 90 auto result = new G4ReactionProduct;
91 auto Z = static_cast<G4int>(massCode / 1000) 91 auto Z = static_cast<G4int>(massCode / 1000);
92 auto A = static_cast<G4int>(massCode - 1000 92 auto A = static_cast<G4int>(massCode - 1000 * Z);
93 93
94 if (massCode == 0) { 94 if (massCode == 0) {
95 result->SetDefinition(G4Gamma::Gamma()); 95 result->SetDefinition(G4Gamma::Gamma());
96 } 96 }
97 else if (A == 0) { 97 else if (A == 0) {
98 result->SetDefinition(G4Electron::Electron 98 result->SetDefinition(G4Electron::Electron());
99 if (Z == 1) result->SetDefinition(G4Positr 99 if (Z == 1) result->SetDefinition(G4Positron::Positron());
100 } 100 }
101 else if (A == 1) { 101 else if (A == 1) {
102 result->SetDefinition(G4Neutron::Neutron() 102 result->SetDefinition(G4Neutron::Neutron());
103 if (Z == 1) result->SetDefinition(G4Proton 103 if (Z == 1) result->SetDefinition(G4Proton::Proton());
104 } 104 }
105 else if (A == 2) { 105 else if (A == 2) {
106 result->SetDefinition(G4Deuteron::Deuteron 106 result->SetDefinition(G4Deuteron::Deuteron());
107 } 107 }
108 else if (A == 3) { 108 else if (A == 3) {
109 result->SetDefinition(G4Triton::Triton()); 109 result->SetDefinition(G4Triton::Triton());
110 if (Z == 2) result->SetDefinition(G4He3::H 110 if (Z == 2) result->SetDefinition(G4He3::He3());
111 } 111 }
112 else if (A == 4) { 112 else if (A == 4) {
113 result->SetDefinition(G4Alpha::Alpha()); 113 result->SetDefinition(G4Alpha::Alpha());
114 if (Z != 2) throw G4HadronicException(__FI 114 if (Z != 2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
115 } 115 }
116 else { 116 else {
117 throw G4HadronicException(__FILE__, __LINE 117 throw G4HadronicException(__FILE__, __LINE__,
118 "G4ParticleHPDis 118 "G4ParticleHPDiscreteTwoBody: Unknown ion case 2");
119 } 119 }
120 120
121 // get cosine(theta) 121 // get cosine(theta)
122 G4int i(0), it(0); 122 G4int i(0), it(0);
123 G4double cosTh(0); 123 G4double cosTh(0);
124 for (i = 0; i < nEnergy; i++) { 124 for (i = 0; i < nEnergy; i++) {
125 it = i; 125 it = i;
126 if (theCoeff[i].GetEnergy() > anEnergy) br 126 if (theCoeff[i].GetEnergy() > anEnergy) break;
127 } 127 }
128 if (it == 0 || it == nEnergy - 1) { 128 if (it == 0 || it == nEnergy - 1) {
129 if (theCoeff[it].GetRepresentation() == 0) 129 if (theCoeff[it].GetRepresentation() == 0) {
130 // TK Legendre expansion 130 // TK Legendre expansion
131 G4ParticleHPLegendreStore theStore(1); 131 G4ParticleHPLegendreStore theStore(1);
132 theStore.SetCoeff(0, theCoeff); 132 theStore.SetCoeff(0, theCoeff);
133 theStore.SetManager(theManager); 133 theStore.SetManager(theManager);
134 // cosTh = theStore.SampleMax(anEnergy); 134 // cosTh = theStore.SampleMax(anEnergy);
135 // 080612TK contribution from Benoit Pir 135 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
136 cosTh = theStore.SampleDiscreteTwoBody(a 136 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
137 } 137 }
138 else if (theCoeff[it].GetRepresentation() 138 else if (theCoeff[it].GetRepresentation() == 12) // means LINLIN
139 { 139 {
140 G4ParticleHPVector theStore; 140 G4ParticleHPVector theStore;
141 G4InterpolationManager aManager; 141 G4InterpolationManager aManager;
142 aManager.Init(LINLIN, theCoeff[it].GetNu 142 aManager.Init(LINLIN, theCoeff[it].GetNumberOfPoly() / 2);
143 theStore.SetInterpolationManager(aManage 143 theStore.SetInterpolationManager(aManager);
144 for (i = 0; i < theCoeff[it].GetNumberOf 144 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
145 // 101110 145 // 101110
146 // theStore.SetX(i, theCoeff[it].GetCo 146 // theStore.SetX(i, theCoeff[it].GetCoeff(i));
147 // theStore.SetY(i, theCoeff[it].GetCo 147 // theStore.SetY(i, theCoeff[it].GetCoeff(i));
148 theStore.SetX(i / 2, theCoeff[it].GetC 148 theStore.SetX(i / 2, theCoeff[it].GetCoeff(i));
149 theStore.SetY(i / 2, theCoeff[it].GetC 149 theStore.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
150 } 150 }
151 cosTh = theStore.Sample(); 151 cosTh = theStore.Sample();
152 } 152 }
153 else if (theCoeff[it].GetRepresentation() 153 else if (theCoeff[it].GetRepresentation() == 14) // this is LOGLIN
154 { 154 {
155 G4ParticleHPVector theStore; 155 G4ParticleHPVector theStore;
156 G4InterpolationManager aManager; 156 G4InterpolationManager aManager;
157 aManager.Init(LOGLIN, theCoeff[it].GetNu 157 aManager.Init(LOGLIN, theCoeff[it].GetNumberOfPoly() / 2);
158 theStore.SetInterpolationManager(aManage 158 theStore.SetInterpolationManager(aManager);
159 for (i = 0; i < theCoeff[it].GetNumberOf 159 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
160 // 101110 160 // 101110
161 // theStore.SetX(i, theCoeff[it].GetCo 161 // theStore.SetX(i, theCoeff[it].GetCoeff(i));
162 // theStore.SetY(i, theCoeff[it].GetCo 162 // theStore.SetY(i, theCoeff[it].GetCoeff(i));
163 theStore.SetX(i / 2, theCoeff[it].GetC 163 theStore.SetX(i / 2, theCoeff[it].GetCoeff(i));
164 theStore.SetY(i / 2, theCoeff[it].GetC 164 theStore.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
165 } 165 }
166 cosTh = theStore.Sample(); 166 cosTh = theStore.Sample();
167 } 167 }
168 else { 168 else {
169 throw G4HadronicException(__FILE__, __LI 169 throw G4HadronicException(__FILE__, __LINE__,
170 "unknown repre 170 "unknown representation type in Two-body scattering");
171 } 171 }
172 } 172 }
173 else { 173 else {
174 if (!bCheckDiffCoeffRepr 174 if (!bCheckDiffCoeffRepr
175 || theCoeff[it].GetRepresentation() == 175 || theCoeff[it].GetRepresentation() == theCoeff[it - 1].GetRepresentation())
176 { 176 {
177 if (theCoeff[it].GetRepresentation() == 177 if (theCoeff[it].GetRepresentation() == 0) {
178 // TK Legendre expansion 178 // TK Legendre expansion
179 G4ParticleHPLegendreStore theStore(2); 179 G4ParticleHPLegendreStore theStore(2);
180 theStore.SetCoeff(0, &(theCoeff[it - 1 180 theStore.SetCoeff(0, &(theCoeff[it - 1]));
181 theStore.SetCoeff(1, &(theCoeff[it])); 181 theStore.SetCoeff(1, &(theCoeff[it]));
182 G4InterpolationManager aManager; 182 G4InterpolationManager aManager;
183 aManager.Init(theManager.GetScheme(it) 183 aManager.Init(theManager.GetScheme(it), 2);
184 theStore.SetManager(aManager); 184 theStore.SetManager(aManager);
185 // 080709 TKDB 185 // 080709 TKDB
186 cosTh = theStore.SampleDiscreteTwoBody 186 cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
187 } 187 }
188 else if (theCoeff[it].GetRepresentation( 188 else if (theCoeff[it].GetRepresentation() == 12) // LINLIN
189 { 189 {
190 G4ParticleHPVector theBuff1; 190 G4ParticleHPVector theBuff1;
191 G4InterpolationManager aManager1; 191 G4InterpolationManager aManager1;
192 aManager1.Init(LINLIN, theCoeff[it - 1 192 aManager1.Init(LINLIN, theCoeff[it - 1].GetNumberOfPoly() / 2);
193 theBuff1.SetInterpolationManager(aMana 193 theBuff1.SetInterpolationManager(aManager1);
194 for (i = 0; i < theCoeff[it - 1].GetNu 194 for (i = 0; i < theCoeff[it - 1].GetNumberOfPoly(); i += 2) {
195 // 101110 195 // 101110
196 theBuff1.SetX(i / 2, theCoeff[it - 1 196 theBuff1.SetX(i / 2, theCoeff[it - 1].GetCoeff(i));
197 theBuff1.SetY(i / 2, theCoeff[it - 1 197 theBuff1.SetY(i / 2, theCoeff[it - 1].GetCoeff(i + 1));
198 } 198 }
199 G4ParticleHPVector theBuff2; 199 G4ParticleHPVector theBuff2;
200 G4InterpolationManager aManager2; 200 G4InterpolationManager aManager2;
201 aManager2.Init(LINLIN, theCoeff[it].Ge 201 aManager2.Init(LINLIN, theCoeff[it].GetNumberOfPoly() / 2);
202 theBuff2.SetInterpolationManager(aMana 202 theBuff2.SetInterpolationManager(aManager2);
203 for (i = 0; i < theCoeff[it].GetNumber 203 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
204 theBuff2.SetX(i / 2, theCoeff[it].Ge 204 theBuff2.SetX(i / 2, theCoeff[it].GetCoeff(i));
205 theBuff2.SetY(i / 2, theCoeff[it].Ge 205 theBuff2.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
206 } 206 }
207 207
208 G4double x1 = theCoeff[it - 1].GetEner 208 G4double x1 = theCoeff[it - 1].GetEnergy();
209 G4double x2 = theCoeff[it].GetEnergy() 209 G4double x2 = theCoeff[it].GetEnergy();
210 G4double x = anEnergy; 210 G4double x = anEnergy;
211 G4double y1, y2, y, mu; 211 G4double y1, y2, y, mu;
212 212
213 G4ParticleHPVector theStore1; 213 G4ParticleHPVector theStore1;
214 theStore1.SetInterpolationManager(aMan 214 theStore1.SetInterpolationManager(aManager1);
215 G4ParticleHPVector theStore2; 215 G4ParticleHPVector theStore2;
216 theStore2.SetInterpolationManager(aMan 216 theStore2.SetInterpolationManager(aManager2);
217 G4ParticleHPVector theStore; 217 G4ParticleHPVector theStore;
218 218
219 // for fixed mu get p1, p2 and interpo 219 // for fixed mu get p1, p2 and interpolate according to x
220 for (i = 0; i < theBuff1.GetVectorLeng 220 for (i = 0; i < theBuff1.GetVectorLength(); ++i) {
221 mu = theBuff1.GetX(i); 221 mu = theBuff1.GetX(i);
222 y1 = theBuff1.GetY(i); 222 y1 = theBuff1.GetY(i);
223 y2 = theBuff2.GetY(mu); 223 y2 = theBuff2.GetY(mu);
224 y = theInt.Interpolate(theManager.Ge 224 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
225 theStore1.SetData(i, mu, y); 225 theStore1.SetData(i, mu, y);
226 } 226 }
227 for (i = 0; i < theBuff2.GetVectorLeng 227 for (i = 0; i < theBuff2.GetVectorLength(); ++i) {
228 mu = theBuff2.GetX(i); 228 mu = theBuff2.GetX(i);
229 y1 = theBuff2.GetY(i); 229 y1 = theBuff2.GetY(i);
230 y2 = theBuff1.GetY(mu); 230 y2 = theBuff1.GetY(mu);
231 y = theInt.Interpolate(theManager.Ge 231 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
232 theStore2.SetData(i, mu, y); 232 theStore2.SetData(i, mu, y);
233 } 233 }
234 theStore.Merge(&theStore1, &theStore2) 234 theStore.Merge(&theStore1, &theStore2); // merge takes care of interpolationschemes
235 cosTh = theStore.Sample(); 235 cosTh = theStore.Sample();
236 } 236 }
237 else if (theCoeff[it].GetRepresentation( 237 else if (theCoeff[it].GetRepresentation() == 14) // TK LOG_LIN
238 { 238 {
239 G4ParticleHPVector theBuff1; 239 G4ParticleHPVector theBuff1;
240 G4InterpolationManager aManager1; 240 G4InterpolationManager aManager1;
241 aManager1.Init(LOGLIN, theCoeff[it - 1 241 aManager1.Init(LOGLIN, theCoeff[it - 1].GetNumberOfPoly() / 2);
242 theBuff1.SetInterpolationManager(aMana 242 theBuff1.SetInterpolationManager(aManager1);
243 for (i = 0; i < theCoeff[it - 1].GetNu 243 for (i = 0; i < theCoeff[it - 1].GetNumberOfPoly(); i += 2) {
244 // 101110 244 // 101110
245 theBuff1.SetX(i / 2, theCoeff[it - 1 245 theBuff1.SetX(i / 2, theCoeff[it - 1].GetCoeff(i));
246 theBuff1.SetY(i / 2, theCoeff[it - 1 246 theBuff1.SetY(i / 2, theCoeff[it - 1].GetCoeff(i + 1));
247 } 247 }
248 248
249 G4ParticleHPVector theBuff2; 249 G4ParticleHPVector theBuff2;
250 G4InterpolationManager aManager2; 250 G4InterpolationManager aManager2;
251 aManager2.Init(LOGLIN, theCoeff[it].Ge 251 aManager2.Init(LOGLIN, theCoeff[it].GetNumberOfPoly() / 2);
252 theBuff2.SetInterpolationManager(aMana 252 theBuff2.SetInterpolationManager(aManager2);
253 for (i = 0; i < theCoeff[it].GetNumber 253 for (i = 0; i < theCoeff[it].GetNumberOfPoly(); i += 2) {
254 // 101110 254 // 101110
255 theBuff2.SetX(i / 2, theCoeff[it].Ge 255 theBuff2.SetX(i / 2, theCoeff[it].GetCoeff(i));
256 theBuff2.SetY(i / 2, theCoeff[it].Ge 256 theBuff2.SetY(i / 2, theCoeff[it].GetCoeff(i + 1));
257 } 257 }
258 258
259 G4double x1 = theCoeff[it - 1].GetEner 259 G4double x1 = theCoeff[it - 1].GetEnergy();
260 G4double x2 = theCoeff[it].GetEnergy() 260 G4double x2 = theCoeff[it].GetEnergy();
261 G4double x = anEnergy; 261 G4double x = anEnergy;
262 G4double y1, y2, y, mu; 262 G4double y1, y2, y, mu;
263 263
264 G4ParticleHPVector theStore1; 264 G4ParticleHPVector theStore1;
265 theStore1.SetInterpolationManager(aMan 265 theStore1.SetInterpolationManager(aManager1);
266 G4ParticleHPVector theStore2; 266 G4ParticleHPVector theStore2;
267 theStore2.SetInterpolationManager(aMan 267 theStore2.SetInterpolationManager(aManager2);
268 G4ParticleHPVector theStore; 268 G4ParticleHPVector theStore;
269 269
270 // for fixed mu get p1, p2 and interpo 270 // for fixed mu get p1, p2 and interpolate according to x
271 for (i = 0; i < theBuff1.GetVectorLeng 271 for (i = 0; i < theBuff1.GetVectorLength(); i++) {
272 mu = theBuff1.GetX(i); 272 mu = theBuff1.GetX(i);
273 y1 = theBuff1.GetY(i); 273 y1 = theBuff1.GetY(i);
274 y2 = theBuff2.GetY(mu); 274 y2 = theBuff2.GetY(mu);
275 y = theInt.Interpolate(theManager.Ge 275 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
276 theStore1.SetData(i, mu, y); 276 theStore1.SetData(i, mu, y);
277 } 277 }
278 for (i = 0; i < theBuff2.GetVectorLeng 278 for (i = 0; i < theBuff2.GetVectorLength(); i++) {
279 mu = theBuff2.GetX(i); 279 mu = theBuff2.GetX(i);
280 y1 = theBuff2.GetY(i); 280 y1 = theBuff2.GetY(i);
281 y2 = theBuff1.GetY(mu); 281 y2 = theBuff1.GetY(mu);
282 y = theInt.Interpolate(theManager.Ge 282 y = theInt.Interpolate(theManager.GetScheme(it), x, x1, x2, y1, y2);
283 theStore2.SetData(i, mu, y); 283 theStore2.SetData(i, mu, y);
284 } 284 }
285 theStore.Merge(&theStore1, &theStore2) 285 theStore.Merge(&theStore1, &theStore2);
286 cosTh = theStore.Sample(); 286 cosTh = theStore.Sample();
287 } 287 }
288 else { 288 else {
289 throw G4HadronicException(__FILE__, __ 289 throw G4HadronicException(__FILE__, __LINE__,
290 "Two neighbo 290 "Two neighbouring distributions with different interpolation");
291 } 291 }
292 } 292 }
293 else { 293 else {
294 G4cout << " theCoeff[it].GetRepresent ME 294 G4cout << " theCoeff[it].GetRepresent MEM " << it << " " << &theCoeff[it] << " "
295 << &theCoeff[it - 1] << G4endl; 295 << &theCoeff[it - 1] << G4endl;
296 G4cout << " theCoeff[it].GetRepresent " 296 G4cout << " theCoeff[it].GetRepresent " << it << " " << theCoeff[it].GetRepresentation()
297 << " != " << theCoeff[it - 1].Get 297 << " != " << theCoeff[it - 1].GetRepresentation() << G4endl;
298 298
299 throw G4HadronicException(__FILE__, __LI 299 throw G4HadronicException(__FILE__, __LINE__,
300 "unknown repre 300 "unknown representation type in Two-body scattering, case 2");
301 } 301 }
302 } 302 }
303 303
304 // now get the energy from kinematics and Q- 304 // now get the energy from kinematics and Q-value.
305 305
306 // G4double restEnergy = anEnergy+GetQValue( 306 // G4double restEnergy = anEnergy+GetQValue();
307 307
308 // assumed to be in CMS @@@@@@@@@@@@@@@@@ 308 // assumed to be in CMS @@@@@@@@@@@@@@@@@
309 309
310 // 080612TK contribution from Benoit Pirard 310 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
311 // G4double residualMass = GetTarget()->Ge 311 // G4double residualMass = GetTarget()->GetMass() + GetNeutron()->GetMass()
312 // - result->GetMass 312 // - result->GetMass() - GetQValue();
313 // G4double kinE = restEnergy/(1+result->Get 313 // G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
314 G4double A1 = GetTarget()->GetMass() / GetPr 314 G4double A1 = GetTarget()->GetMass() / GetProjectileRP()->GetMass();
315 G4double A1prim = result->GetMass() / GetPro 315 G4double A1prim = result->GetMass() / GetProjectileRP()->GetMass();
316 // G4double E1 = (A1+1)*(A1+1)/A1/A1*an 316 // G4double E1 = (A1+1)*(A1+1)/A1/A1*anEnergy;
317 // Bug fix Bugzilla #1815 317 // Bug fix Bugzilla #1815
318 G4double E1 = anEnergy; 318 G4double E1 = anEnergy;
319 G4double kinE = (A1 + 1 - A1prim) / (A1 + 1) 319 G4double kinE = (A1 + 1 - A1prim) / (A1 + 1) / (A1 + 1) * (A1 * E1 + (1 + A1) * GetQValue());
320 320
321 result->SetKineticEnergy(kinE); // non rela 321 result->SetKineticEnergy(kinE); // non relativistic @@
322 if (cosTh > 1.0) { cosTh = 1.0; } 322 if (cosTh > 1.0) { cosTh = 1.0; }
323 else if(cosTh < -1.0) { cosTh = -1.0; } 323 else if(cosTh < -1.0) { cosTh = -1.0; }
324 G4double phi = CLHEP::twopi * G4UniformRand( 324 G4double phi = CLHEP::twopi * G4UniformRand();
325 G4double sinth = std::sqrt((1.0 + cosTh)*(1. 325 G4double sinth = std::sqrt((1.0 + cosTh)*(1.0 - cosTh));
326 G4double mtot = result->GetTotalMomentum(); 326 G4double mtot = result->GetTotalMomentum();
327 G4ThreeVector tempVector(mtot * sinth * std: 327 G4ThreeVector tempVector(mtot * sinth * std::cos(phi), mtot * sinth * std::sin(phi),
328 mtot * cosTh); 328 mtot * cosTh);
329 result->SetMomentum(tempVector); 329 result->SetMomentum(tempVector);
330 return result; 330 return result;
331 } 331 }
332 332