| 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 // neutron_hp -- source file 26 // neutron_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 // 080808 Bug fix in serching mu bin and index 30 // 080808 Bug fix in serching mu bin and index for theBuff2b by T. Koi
31 // 31 //
32 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
33 // 33 //
34 // E. Mendoza, Nov. 2020 - bug fix <<
35 // <<
36 #include "G4ParticleHPLabAngularEnergy.hh" 34 #include "G4ParticleHPLabAngularEnergy.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 "G4PhysicalConstants.hh" 35 #include "G4PhysicalConstants.hh"
>> 36 #include "G4SystemOfUnits.hh"
>> 37 #include "Randomize.hh"
>> 38 #include "G4Gamma.hh"
>> 39 #include "G4Electron.hh"
45 #include "G4Positron.hh" 40 #include "G4Positron.hh"
>> 41 #include "G4Neutron.hh"
46 #include "G4Proton.hh" 42 #include "G4Proton.hh"
47 #include "G4SystemOfUnits.hh" << 43 #include "G4Deuteron.hh"
48 #include "G4Triton.hh" 44 #include "G4Triton.hh"
49 #include "Randomize.hh" << 45 #include "G4He3.hh"
>> 46 #include "G4Alpha.hh"
50 47
51 void G4ParticleHPLabAngularEnergy::Init(std::i << 48 void G4ParticleHPLabAngularEnergy::Init(std::istream & aDataFile)
52 { 49 {
53 aDataFile >> nEnergies; 50 aDataFile >> nEnergies;
54 theManager.Init(aDataFile); 51 theManager.Init(aDataFile);
55 const std::size_t esize = nEnergies > 0 ? nE << 52 theEnergies = new G4double[nEnergies];
56 theEnergies = new G4double[esize]; << 53 nCosTh = new G4int[nEnergies];
57 nCosTh = new G4int[esize]; << 54 theData = new G4ParticleHPVector * [nEnergies];
58 theData = new G4ParticleHPVector*[esize]; << 55 theSecondManager = new G4InterpolationManager [nEnergies];
59 theSecondManager = new G4InterpolationManage << 56 for(G4int i=0; i<nEnergies; i++)
60 for (G4int i = 0; i < nEnergies; ++i) { << 57 {
61 aDataFile >> theEnergies[i]; 58 aDataFile >> theEnergies[i];
62 theEnergies[i] *= eV; << 59 theEnergies[i]*=eV;
63 aDataFile >> nCosTh[i]; 60 aDataFile >> nCosTh[i];
64 theSecondManager[i].Init(aDataFile); << 61 theSecondManager[i].Init(aDataFile);
65 const std::size_t dsize = nCosTh[i] > 0 ? << 62 theData[i] = new G4ParticleHPVector[nCosTh[i]];
66 theData[i] = new G4ParticleHPVector[dsize] <<
67 G4double label; 63 G4double label;
68 for (std::size_t ii = 0; ii < dsize; ++ii) << 64 for(G4int ii=0; ii<nCosTh[i]; ii++)
>> 65 {
69 aDataFile >> label; 66 aDataFile >> label;
70 theData[i][ii].SetLabel(label); 67 theData[i][ii].SetLabel(label);
71 theData[i][ii].Init(aDataFile, eV); 68 theData[i][ii].Init(aDataFile, eV);
72 } 69 }
73 } 70 }
74 } 71 }
75 72
76 G4ReactionProduct* G4ParticleHPLabAngularEnerg << 73 G4ReactionProduct * G4ParticleHPLabAngularEnergy::Sample(G4double anEnergy, G4double massCode, G4double )
77 <<
78 { 74 {
79 auto result = new G4ReactionProduct; << 75 G4ReactionProduct * result = new G4ReactionProduct;
80 auto Z = static_cast<G4int>(massCode / 1000) << 76 G4int Z = static_cast<G4int>(massCode/1000);
81 auto A = static_cast<G4int>(massCode - 1000 << 77 G4int A = static_cast<G4int>(massCode-1000*Z);
82 << 78
83 if (massCode == 0) { << 79 if(massCode==0)
84 result->SetDefinition(G4Gamma::Gamma()); << 80 {
85 } << 81 result->SetDefinition(G4Gamma::Gamma());
86 else if (A == 0) { << 82 }
87 result->SetDefinition(G4Electron::Electron << 83 else if(A==0)
88 if (Z == 1) result->SetDefinition(G4Positr << 84 {
89 } << 85 result->SetDefinition(G4Electron::Electron());
90 else if (A == 1) { << 86 if(Z==1) result->SetDefinition(G4Positron::Positron());
91 result->SetDefinition(G4Neutron::Neutron() << 87 }
92 if (Z == 1) result->SetDefinition(G4Proton << 88 else if(A==1)
93 } << 89 {
94 else if (A == 2) { << 90 result->SetDefinition(G4Neutron::Neutron());
95 result->SetDefinition(G4Deuteron::Deuteron << 91 if(Z==1) result->SetDefinition(G4Proton::Proton());
96 } << 92 }
97 else if (A == 3) { << 93 else if(A==2)
98 result->SetDefinition(G4Triton::Triton()); << 94 {
99 if (Z == 2) result->SetDefinition(G4He3::H << 95 result->SetDefinition(G4Deuteron::Deuteron());
100 } << 96 }
101 else if (A == 4) { << 97 else if(A==3)
102 result->SetDefinition(G4Alpha::Alpha()); << 98 {
103 if (Z != 2) throw G4HadronicException(__FI << 99 result->SetDefinition(G4Triton::Triton());
104 } << 100 if(Z==2) result->SetDefinition(G4He3::He3());
105 else { << 101 }
106 throw G4HadronicException(__FILE__, __LINE << 102 else if(A==4)
107 "G4ParticleHPLab << 103 {
108 } << 104 result->SetDefinition(G4Alpha::Alpha());
109 << 105 if(Z!=2) throw G4HadronicException(__FILE__, __LINE__, "Unknown ion case 1");
110 // get theta, E << 106 }
111 G4double cosTh, secEnergy; << 107 else
112 G4int i, it(0); << 108 {
113 // find the energy bin << 109 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPLabAngularEnergy: Unknown ion case 2");
114 for (i = 0; i < nEnergies; i++) { << 110 }
115 it = i; << 111
116 if (anEnergy < theEnergies[i]) break; << 112 // get theta, E
117 } << 113 G4double cosTh, secEnergy;
118 << 114 G4int i, it(0);
119 if (it == 0) // it marks the energy bin --> << 115 // find the energy bin
120 // to high energies (??) << 116 for(i=0; i<nEnergies; i++)
121 { << 117 {
122 // Do not sample between incident neutron << 118 it = i;
123 //---------------------------------------- << 119 if ( anEnergy < theEnergies[i] ) break;
124 // CosTheta: << 120 }
125 G4ParticleHPVector theThVec; << 121 //080808
126 theThVec.SetInterpolationManager(theSecond << 122 //if ( it == 0 || it == nEnergies-1 ) // it marks the energy bin
127 for (i = 0; i < nCosTh[it]; i++) { << 123 if ( it == 0 ) // it marks the energy bin
128 theThVec.SetX(i, theData[it][i].GetLabel << 124 {
129 theThVec.SetY(i, theData[it][i].GetInteg << 125 G4cout << "080808 Something unexpected is happen in G4ParticleHPLabAngularEnergy " << G4endl;
130 } << 126 // integrate the prob for each costh, and select theta.
131 cosTh = theThVec.Sample(); << 127 G4double * running = new G4double [nCosTh[it]];
132 //---------------------------------------- << 128 running[0]=0;
133 << 129 for(i=0;i<nCosTh[it]; i++)
134 //---------------------------------------- << 130 {
135 // Now the secondary energy: << 131 if(i!=0) running[i] = running[i-1];
136 G4double x, x1, x2, y1, y2, y, E; << 132 running[i]+=theData[it][i].GetIntegral(); // Does interpolated integral.
137 G4int i1(0); << 133 }
138 for (i = 1; i < nCosTh[it]; i++) { << 134 G4double random = running[nCosTh[it]-1]*G4UniformRand();
139 i1 = i; << 135 G4int ith(0);
140 if (cosTh < theData[it][i].GetLabel()) b << 136 for(i=0;i<nCosTh[it]; i++)
141 } << 137 {
142 // now get the prob at this energy for the << 138 ith = i;
143 x = cosTh; << 139 if(random<running[i]) break;
144 x1 = theData[it][i1 - 1].GetLabel(); << 140 }
145 x2 = theData[it][i1].GetLabel(); << 141 //080807
146 G4ParticleHPVector theBuff1a; << 142 //if ( ith == 0 || ith == nCosTh[it]-1 ) //ith marks the angluar bin
147 theBuff1a.SetInterpolationManager(theData[ << 143 if ( ith == 0 ) //ith marks the angluar bin
148 for (i = 0; i < theData[it][i1 - 1].GetVec << 144 {
149 E = theData[it][i1 - 1].GetX(i); << 145 cosTh = theData[it][ith].GetLabel();
150 y1 = theData[it][i1 - 1].GetY(i); << 146 secEnergy = theData[it][ith].Sample();
151 y2 = theData[it][i1].GetY(E); << 147 currentMeanEnergy = theData[it][ith].GetMeanX();
152 y = theInt.Interpolate(theSecondManager[ << 148 }
153 theBuff1a.SetData(i, E, y); << 149 else
154 } << 150 {
155 G4ParticleHPVector theBuff2a; << 151 //080808
156 theBuff2a.SetInterpolationManager(theData[ << 152 //G4double x1 = theData[it][ith-1].GetIntegral();
157 for (i = 0; i < theData[it][i1].GetVectorL << 153 //G4double x2 = theData[it][ith].GetIntegral();
158 E = theData[it][i1].GetX(i); << 154 G4double x1 = running [ ith-1 ];
159 y1 = theData[it][i1 - 1].GetY(E); << 155 G4double x2 = running [ ith ];
160 y2 = theData[it][i1].GetY(i); << 156 G4double x = random;
161 y = theInt.Interpolate(theSecondManager[ << 157 G4double y1 = theData[it][ith-1].GetLabel();
162 theBuff2a.SetData(i, E, y); // wrong E, << 158 G4double y2 = theData[it][ith].GetLabel();
163 } << 159 cosTh = theInt.Interpolate(theSecondManager[it].GetInverseScheme(ith),
164 G4ParticleHPVector theStore; << 160 x, x1, x2, y1, y2);
165 theStore.Merge(&theBuff1a, &theBuff2a); << 161 G4ParticleHPVector theBuff1;
166 secEnergy = theStore.Sample(); << 162 theBuff1.SetInterpolationManager(theData[it][ith-1].GetInterpolationManager());
167 currentMeanEnergy = theStore.GetMeanX(); << 163 G4ParticleHPVector theBuff2;
168 //---------------------------------------- << 164 theBuff2.SetInterpolationManager(theData[it][ith].GetInterpolationManager());
169 } << 165 x1=y1;
170 else // this is the small big else. << 166 x2=y2;
171 { << 167 G4double y, mu;
172 G4double x, x1, x2, y1, y2, y, tmp, E; << 168 for(i=0;i<theData[it][ith-1].GetVectorLength(); i++)
173 // integrate the prob for each costh, and << 169 {
174 G4ParticleHPVector run1; << 170 mu = theData[it][ith-1].GetX(i);
175 for (i = 0; i < nCosTh[it - 1]; i++) { << 171 y1 = theData[it][ith-1].GetY(i);
176 run1.SetX(i, theData[it - 1][i].GetLabel << 172 y2 = theData[it][ith].GetY(mu);
177 run1.SetY(i, theData[it - 1][i].GetInteg << 173
178 } << 174 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
179 G4ParticleHPVector run2; << 175 cosTh, x1,x2,y1,y2);
180 for (i = 0; i < nCosTh[it]; i++) { << 176 theBuff1.SetData(i, mu, y);
181 run2.SetX(i, theData[it][i].GetLabel()); << 177 }
182 run2.SetY(i, theData[it][i].GetIntegral( << 178 for(i=0;i<theData[it][ith].GetVectorLength(); i++)
183 } << 179 {
184 // get the distributions for the correct n << 180 mu = theData[it][ith].GetX(i);
185 x = anEnergy; << 181 y1 = theData[it][ith-1].GetY(mu);
186 x1 = theEnergies[it - 1]; << 182 y2 = theData[it][ith].GetY(i);
187 x2 = theEnergies[it]; << 183 y = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
188 G4ParticleHPVector thBuff1; // to be inte << 184 cosTh, x1,x2,y1,y2);
189 thBuff1.SetInterpolationManager(theSecondM << 185 theBuff2.SetData(i, mu, y);
190 for (i = 0; i < run1.GetVectorLength(); i+ << 186 }
191 tmp = run1.GetX(i); // theta << 187 G4ParticleHPVector theStore;
192 y1 = run1.GetY(i); // integral << 188 theStore.Merge(&theBuff1, &theBuff2);
193 y2 = run2.GetY(tmp); << 189 secEnergy = theStore.Sample();
194 y = theInt.Interpolate(theManager.GetSch << 190 currentMeanEnergy = theStore.GetMeanX();
195 thBuff1.SetData(i, tmp, y); << 191 }
196 } << 192 delete [] running;
197 G4ParticleHPVector thBuff2; << 193 }
198 thBuff2.SetInterpolationManager(theSecondM << 194 else // this is the small big else.
199 for (i = 0; i < run2.GetVectorLength(); i+ << 195 {
200 tmp = run2.GetX(i); // theta << 196 G4double x, x1, x2, y1, y2, y, tmp, E;
201 y1 = run1.GetY(tmp); // integral << 197 // integrate the prob for each costh, and select theta.
202 y2 = run2.GetY(i); << 198 G4ParticleHPVector run1;
203 y = theInt.Interpolate(theManager.GetSch << 199 run1.SetY(0, 0.);
204 thBuff2.SetData(i, tmp, y); << 200 for(i=0;i<nCosTh[it-1]; i++)
205 } << 201 {
206 G4ParticleHPVector theThVec; << 202 if(i!=0) run1.SetY(i, run1.GetY(i-1));
207 theThVec.Merge(&thBuff1, &thBuff2); // ta << 203 run1.SetX(i, theData[it-1][i].GetLabel());
208 cosTh = theThVec.Sample(); << 204 run1.SetY(i, run1.GetY(i)+theData[it-1][i].GetIntegral());
209 /* << 205 }
210 if(massCode>0.99 && massCode<1.01){theThVe << 206 G4ParticleHPVector run2;
211 G4double random = (theThVec.GetY(theThVec. << 207 run2.SetY(0, 0.);
212 -theThVec.GetY(0)) *G << 208 for(i=0;i<nCosTh[it]; i++)
213 G4cout<<" -- "<<theThVec.GetY(theThVec.Get << 209 {
214 "<<random<<G4endl; G4int ith(0); for(i=1;i << 210 if(i!=0) run2.SetY(i, run2.GetY(i-1));
215 { << 211 run2.SetX(i, theData[it][i].GetLabel());
216 ith = i; << 212 run2.SetY(i, run2.GetY(i)+theData[it][i].GetIntegral());
217 if(random<theThVec.GetY(i)-theThVec.GetY << 213 }
218 } << 214 // get the distributions for the correct neutron energy
219 { << 215 x = anEnergy;
220 // calculate theta << 216 x1 = theEnergies[it-1];
221 G4double xx, xx1, xx2, yy1, yy2; << 217 x2 = theEnergies[it];
222 xx = random; << 218 G4ParticleHPVector thBuff1; // to be interpolated as run1.
223 xx1 = theThVec.GetY(ith-1)-theThVec.GetY << 219 thBuff1.SetInterpolationManager(theSecondManager[it-1]);
224 xx2 = theThVec.GetY(ith)-theThVec.GetY(0 << 220 for(i=0; i<run1.GetVectorLength(); i++)
225 yy1 = theThVec.GetX(ith-1); // std::cos( << 221 {
226 yy2 = theThVec.GetX(ith); << 222 tmp = run1.GetX(i); //theta
227 cosTh = theInt.Interpolate(theSecondMana << 223 y1 = run1.GetY(i); // integral
228 xx, xx1,xx2,y << 224 y2 = run2.GetY(tmp);
229 } << 225 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
230 */ << 226 thBuff1.SetData(i, tmp, y);
231 G4int i1(0), i2(0); << 227 }
232 // get the indixes of the vectors close to << 228 G4ParticleHPVector thBuff2;
233 // first it-1 !!!! i.e. low in energy << 229 thBuff2.SetInterpolationManager(theSecondManager[it]);
234 for (i = 1; i < nCosTh[it - 1]; i++) { << 230 for(i=0; i<run2.GetVectorLength(); i++)
235 i1 = i; << 231 {
236 if (cosTh < theData[it - 1][i].GetLabel( << 232 tmp = run2.GetX(i); //theta
237 } << 233 y1 = run1.GetY(tmp); // integral
238 // now get the prob at this energy for the << 234 y2 = run2.GetY(i);
239 x = cosTh; << 235 y = theInt.Lin(x, x1,x2,y1,y2);
240 x1 = theData[it - 1][i1 - 1].GetLabel(); << 236 thBuff2.SetData(i, tmp, y);
241 x2 = theData[it - 1][i1].GetLabel(); << 237 }
242 G4ParticleHPVector theBuff1a; << 238 G4ParticleHPVector theThVec;
243 theBuff1a.SetInterpolationManager(theData[ << 239 theThVec.Merge(&thBuff1 ,&thBuff2); // takes care of interpolation
244 for (i = 0; i < theData[it - 1][i1 - 1].Ge << 240 G4double random = (theThVec.GetY(theThVec.GetVectorLength()-1)
245 E = theData[it - 1][i1 - 1].GetX(i); << 241 -theThVec.GetY(0)) *G4UniformRand();
246 y1 = theData[it - 1][i1 - 1].GetY(i); << 242 G4int ith(0);
247 y2 = theData[it - 1][i1].GetY(E); << 243 for(i=1;i<theThVec.GetVectorLength(); i++)
248 y = theInt.Interpolate(theSecondManager[ << 244 {
249 theBuff1a.SetData(i, E, y); // wrong E, << 245 ith = i;
250 } << 246 if(random<theThVec.GetY(i)-theThVec.GetY(0)) break;
251 G4ParticleHPVector theBuff2a; << 247 }
252 theBuff2a.SetInterpolationManager(theData[ << 248 {
253 for (i = 0; i < theData[it - 1][i1].GetVec << 249 // calculate theta
254 E = theData[it - 1][i1].GetX(i); << 250 G4double xx, xx1, xx2, yy1, yy2;
255 y1 = theData[it - 1][i1 - 1].GetY(E); << 251 xx = random;
256 y2 = theData[it - 1][i1].GetY(i); << 252 xx1 = theThVec.GetY(ith-1)-theThVec.GetY(0); // integrals
257 y = theInt.Interpolate(theSecondManager[ << 253 xx2 = theThVec.GetY(ith)-theThVec.GetY(0);
258 theBuff2a.SetData(i, E, y); // wrong E, << 254 yy1 = theThVec.GetX(ith-1); // std::cos(theta)
259 } << 255 yy2 = theThVec.GetX(ith);
260 G4ParticleHPVector theStore1; << 256 cosTh = theInt.Interpolate(theSecondManager[it].GetScheme(ith),
261 theStore1.Merge(&theBuff1a, &theBuff2a); << 257 xx, xx1,xx2,yy1,yy2);
262 << 258 }
263 // get the indixes of the vectors close to << 259 G4int i1(0), i2(0);
264 // then it !!!! i.e. high in energy << 260 // get the indixes of the vectors close to theta for low energy
265 for (i = 1; i < nCosTh[it]; i++) { << 261 // first it-1 !!!! i.e. low in energy
266 i2 = i; << 262 for(i=0; i<nCosTh[it-1]; i++)
267 if (cosTh < theData[it][i2].GetLabel()) << 263 {
268 } // sonderfaelle mit i1 oder i2 head on << 264 i1 = i;
269 x1 = theData[it][i2 - 1].GetLabel(); << 265 if(cosTh<theData[it-1][i].GetLabel()) break;
270 x2 = theData[it][i2].GetLabel(); << 266 }
271 G4ParticleHPVector theBuff1b; << 267 // now get the prob at this energy for the right theta value
272 theBuff1b.SetInterpolationManager(theData[ << 268 x = cosTh;
273 for (i = 0; i < theData[it][i2 - 1].GetVec << 269 x1 = theData[it-1][i1-1].GetLabel();
274 E = theData[it][i2 - 1].GetX(i); << 270 x2 = theData[it-1][i1].GetLabel();
275 y1 = theData[it][i2 - 1].GetY(i); << 271 G4ParticleHPVector theBuff1a;
276 y2 = theData[it][i2].GetY(E); << 272 theBuff1a.SetInterpolationManager(theData[it-1][i1-1].GetInterpolationManager());
277 y = theInt.Interpolate(theSecondManager[ << 273 for(i=0;i<theData[it-1][i1-1].GetVectorLength(); i++)
278 theBuff1b.SetData(i, E, y); // wrong E, << 274 {
279 } << 275 E = theData[it-1][i1-1].GetX(i);
280 G4ParticleHPVector theBuff2b; << 276 y1 = theData[it-1][i1-1].GetY(i);
281 theBuff2b.SetInterpolationManager(theData[ << 277 y2 = theData[it-1][i1].GetY(E);
282 // 080808 i1 -> i2 << 278 y = theInt.Lin(x, x1,x2,y1,y2);
283 // for(i=0;i<theData[it][i1].GetVectorLeng << 279 theBuff1a.SetData(i, E, y); // wrong E, right theta.
284 for (i = 0; i < theData[it][i2].GetVectorL << 280 }
285 // E = theData[it][i1].GetX(i); << 281 G4ParticleHPVector theBuff2a;
286 // y1 = theData[it][i1-1].GetY(E); << 282 theBuff2a.SetInterpolationManager(theData[it-1][i1].GetInterpolationManager());
287 // y2 = theData[it][i1].GetY(i); << 283 for(i=0;i<theData[it-1][i1].GetVectorLength(); i++)
288 E = theData[it][i2].GetX(i); << 284 {
289 y1 = theData[it][i2 - 1].GetY(E); << 285 E = theData[it-1][i1].GetX(i);
290 y2 = theData[it][i2].GetY(i); << 286 y1 = theData[it-1][i1-1].GetY(E);
291 y = theInt.Interpolate(theSecondManager[ << 287 y2 = theData[it-1][i1].GetY(i);
292 theBuff2b.SetData(i, E, y); // wrong E, << 288 y = theInt.Lin(x, x1,x2,y1,y2);
293 } << 289 theBuff2a.SetData(i, E, y); // wrong E, right theta.
294 G4ParticleHPVector theStore2; << 290 }
295 theStore2.Merge(&theBuff1b, &theBuff2b); << 291 G4ParticleHPVector theStore1;
296 // now get to the right energy. << 292 theStore1.Merge(&theBuff1a, &theBuff2a); // wrong E, right theta, complete binning
297 << 293
298 x = anEnergy; << 294 // get the indixes of the vectors close to theta for high energy
299 x1 = theEnergies[it - 1]; << 295 // then it !!!! i.e. high in energy
300 x2 = theEnergies[it]; << 296 for(i=0; i<nCosTh[it]; i++)
301 G4ParticleHPVector theOne1; << 297 {
302 theOne1.SetInterpolationManager(theStore1. << 298 i2 = i;
303 for (i = 0; i < theStore1.GetVectorLength( << 299 if(cosTh<theData[it][i2].GetLabel()) break;
304 E = theStore1.GetX(i); << 300 } // sonderfaelle mit i1 oder i2 head on fehlen. @@@@@
305 y1 = theStore1.GetY(i); << 301 x1 = theData[it][i2-1].GetLabel();
306 y2 = theStore2.GetY(E); << 302 x2 = theData[it][i2].GetLabel();
307 y = theInt.Interpolate(theManager.GetSch << 303 G4ParticleHPVector theBuff1b;
308 theOne1.SetData(i, E, y); // both corre << 304 theBuff1b.SetInterpolationManager(theData[it][i2-1].GetInterpolationManager());
309 } << 305 for(i=0;i<theData[it][i2-1].GetVectorLength(); i++)
310 G4ParticleHPVector theOne2; << 306 {
311 theOne2.SetInterpolationManager(theStore2. << 307 E = theData[it][i2-1].GetX(i);
312 for (i = 0; i < theStore2.GetVectorLength( << 308 y1 = theData[it][i2-1].GetY(i);
313 E = theStore2.GetX(i); << 309 y2 = theData[it][i2].GetY(E);
314 y1 = theStore1.GetY(E); << 310 y = theInt.Lin(x, x1,x2,y1,y2);
315 y2 = theStore2.GetY(i); << 311 theBuff1b.SetData(i, E, y); // wrong E, right theta.
316 y = theInt.Interpolate(theManager.GetSch << 312 }
317 theOne2.SetData(i, E, y); // both corre << 313 G4ParticleHPVector theBuff2b;
318 } << 314 theBuff2b.SetInterpolationManager(theData[it][i2].GetInterpolationManager());
319 G4ParticleHPVector theOne; << 315 //080808 i1 -> i2
320 theOne.Merge(&theOne1, &theOne2); // both << 316 //for(i=0;i<theData[it][i1].GetVectorLength(); i++)
321 << 317 for(i=0;i<theData[it][i2].GetVectorLength(); i++)
322 secEnergy = theOne.Sample(); << 318 {
323 currentMeanEnergy = theOne.GetMeanX(); << 319 //E = theData[it][i1].GetX(i);
324 } << 320 //y1 = theData[it][i1-1].GetY(E);
325 << 321 //y2 = theData[it][i1].GetY(i);
326 // now do random direction in phi, and fill << 322 E = theData[it][i2].GetX(i);
327 << 323 y1 = theData[it][i2-1].GetY(E);
328 result->SetKineticEnergy(secEnergy); << 324 y2 = theData[it][i2].GetY(i);
329 << 325 y = theInt.Lin(x, x1,x2,y1,y2);
330 G4double phi = twopi * G4UniformRand(); << 326 theBuff2b.SetData(i, E, y); // wrong E, right theta.
331 G4double theta = std::acos(cosTh); << 327 }
332 G4double sinth = std::sin(theta); << 328 G4ParticleHPVector theStore2;
333 G4double mtot = result->GetTotalMomentum(); << 329 theStore2.Merge(&theBuff1b, &theBuff2b); // wrong E, right theta, complete binning
334 G4ThreeVector tempVector(mtot * sinth * std: << 330 // now get to the right energy.
335 mtot * std::cos(the << 331
336 result->SetMomentum(tempVector); << 332 x = anEnergy;
337 << 333 x1 = theEnergies[it-1];
338 return result; << 334 x2 = theEnergies[it];
>> 335 G4ParticleHPVector theOne1;
>> 336 theOne1.SetInterpolationManager(theStore1.GetInterpolationManager());
>> 337 for(i=0; i<theStore1.GetVectorLength(); i++)
>> 338 {
>> 339 E = theStore1.GetX(i);
>> 340 y1 = theStore1.GetY(i);
>> 341 y2 = theStore2.GetY(E);
>> 342 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
>> 343 theOne1.SetData(i, E, y); // both correct
>> 344 }
>> 345 G4ParticleHPVector theOne2;
>> 346 theOne2.SetInterpolationManager(theStore2.GetInterpolationManager());
>> 347 for(i=0; i<theStore2.GetVectorLength(); i++)
>> 348 {
>> 349 E = theStore2.GetX(i);
>> 350 y1 = theStore1.GetY(E);
>> 351 y2 = theStore2.GetY(i);
>> 352 y = theInt.Interpolate(theManager.GetScheme(it), x, x1,x2,y1,y2);
>> 353 theOne2.SetData(i, E, y); // both correct
>> 354 }
>> 355 G4ParticleHPVector theOne;
>> 356 theOne.Merge(&theOne1, &theOne2); // both correct, complete binning
>> 357
>> 358 secEnergy = theOne.Sample();
>> 359 currentMeanEnergy = theOne.GetMeanX();
>> 360 }
>> 361
>> 362 // now do random direction in phi, and fill the result.
>> 363
>> 364 result->SetKineticEnergy(secEnergy);
>> 365
>> 366 G4double phi = twopi*G4UniformRand();
>> 367 G4double theta = std::acos(cosTh);
>> 368 G4double sinth = std::sin(theta);
>> 369 G4double mtot = result->GetTotalMomentum();
>> 370 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
>> 371 result->SetMomentum(tempVector);
>> 372
>> 373 return result;
339 } 374 }
340 375