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