| 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 // G4SPSRandomGenerator class implementation 26 // G4SPSRandomGenerator class implementation
27 // 27 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004 28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004
29 // Customer: ESA/ESTEC 29 // Customer: ESA/ESTEC
30 // Revision: Andrea Dotti, SLAC 30 // Revision: Andrea Dotti, SLAC
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include <cmath> 33 #include <cmath>
34 34
35 #include "G4PrimaryParticle.hh" 35 #include "G4PrimaryParticle.hh"
36 #include "G4Event.hh" 36 #include "G4Event.hh"
37 #include "Randomize.hh" 37 #include "Randomize.hh"
38 #include "G4TransportationManager.hh" 38 #include "G4TransportationManager.hh"
39 #include "G4VPhysicalVolume.hh" 39 #include "G4VPhysicalVolume.hh"
40 #include "G4PhysicalVolumeStore.hh" 40 #include "G4PhysicalVolumeStore.hh"
41 #include "G4ParticleTable.hh" 41 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh" 42 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh" 43 #include "G4IonTable.hh"
44 #include "G4Ions.hh" 44 #include "G4Ions.hh"
45 #include "G4TrackingManager.hh" 45 #include "G4TrackingManager.hh"
46 #include "G4Track.hh" 46 #include "G4Track.hh"
47 #include "G4AutoLock.hh" 47 #include "G4AutoLock.hh"
48 48
49 #include "G4SPSRandomGenerator.hh" 49 #include "G4SPSRandomGenerator.hh"
50 50
51 G4SPSRandomGenerator::bweights_t::bweights_t() 51 G4SPSRandomGenerator::bweights_t::bweights_t()
52 { 52 {
53 for (double & i : w) { i = 1; } << 53 for ( G4int i=0 ; i<9 ; ++i ) { w[i] = 1; }
54 } 54 }
55 55
56 G4double& G4SPSRandomGenerator::bweights_t::op 56 G4double& G4SPSRandomGenerator::bweights_t::operator [](const G4int i)
57 { 57 {
58 return w[i]; 58 return w[i];
59 } 59 }
60 60
61 G4SPSRandomGenerator::G4SPSRandomGenerator() 61 G4SPSRandomGenerator::G4SPSRandomGenerator()
62 { 62 {
63 // Initialise all variables 63 // Initialise all variables
64 64
65 // Bias variables 65 // Bias variables
66 // 66 //
67 XBias = false; 67 XBias = false;
68 IPDFXBias = false; 68 IPDFXBias = false;
69 YBias = false; 69 YBias = false;
70 IPDFYBias = false; 70 IPDFYBias = false;
71 ZBias = false; 71 ZBias = false;
72 IPDFZBias = false; 72 IPDFZBias = false;
73 ThetaBias = false; 73 ThetaBias = false;
74 IPDFThetaBias = false; 74 IPDFThetaBias = false;
75 PhiBias = false; 75 PhiBias = false;
76 IPDFPhiBias = false; 76 IPDFPhiBias = false;
77 EnergyBias = false; 77 EnergyBias = false;
78 IPDFEnergyBias = false; 78 IPDFEnergyBias = false;
79 PosThetaBias = false; 79 PosThetaBias = false;
80 IPDFPosThetaBias = false; 80 IPDFPosThetaBias = false;
81 PosPhiBias = false; 81 PosPhiBias = false;
82 IPDFPosPhiBias = false; 82 IPDFPosPhiBias = false;
83 83
84 verbosityLevel = 0; 84 verbosityLevel = 0;
85 G4MUTEXINIT(mutex); 85 G4MUTEXINIT(mutex);
86 } 86 }
87 87
88 G4SPSRandomGenerator::~G4SPSRandomGenerator() 88 G4SPSRandomGenerator::~G4SPSRandomGenerator()
89 { 89 {
90 G4MUTEXDESTROY(mutex); 90 G4MUTEXDESTROY(mutex);
91 } 91 }
92 92
93 // Biasing methods 93 // Biasing methods
94 94
95 void G4SPSRandomGenerator::SetXBias(const G4Th 95 void G4SPSRandomGenerator::SetXBias(const G4ThreeVector& input)
96 { 96 {
97 G4AutoLock l(&mutex); 97 G4AutoLock l(&mutex);
98 G4double ehi, val; 98 G4double ehi, val;
99 ehi = input.x(); 99 ehi = input.x();
100 val = input.y(); 100 val = input.y();
101 XBiasH.InsertValues(ehi, val); 101 XBiasH.InsertValues(ehi, val);
102 XBias = true; 102 XBias = true;
103 } 103 }
104 104
105 void G4SPSRandomGenerator::SetYBias(const G4Th 105 void G4SPSRandomGenerator::SetYBias(const G4ThreeVector& input)
106 { 106 {
107 G4AutoLock l(&mutex); 107 G4AutoLock l(&mutex);
108 G4double ehi, val; 108 G4double ehi, val;
109 ehi = input.x(); 109 ehi = input.x();
110 val = input.y(); 110 val = input.y();
111 YBiasH.InsertValues(ehi, val); 111 YBiasH.InsertValues(ehi, val);
112 YBias = true; 112 YBias = true;
113 } 113 }
114 114
115 void G4SPSRandomGenerator::SetZBias(const G4Th 115 void G4SPSRandomGenerator::SetZBias(const G4ThreeVector& input)
116 { 116 {
117 G4AutoLock l(&mutex); 117 G4AutoLock l(&mutex);
118 G4double ehi, val; 118 G4double ehi, val;
119 ehi = input.x(); 119 ehi = input.x();
120 val = input.y(); 120 val = input.y();
121 ZBiasH.InsertValues(ehi, val); 121 ZBiasH.InsertValues(ehi, val);
122 ZBias = true; 122 ZBias = true;
123 } 123 }
124 124
125 void G4SPSRandomGenerator::SetThetaBias(const 125 void G4SPSRandomGenerator::SetThetaBias(const G4ThreeVector& input)
126 { 126 {
127 G4AutoLock l(&mutex); 127 G4AutoLock l(&mutex);
128 G4double ehi, val; 128 G4double ehi, val;
129 ehi = input.x(); 129 ehi = input.x();
130 val = input.y(); 130 val = input.y();
131 ThetaBiasH.InsertValues(ehi, val); 131 ThetaBiasH.InsertValues(ehi, val);
132 ThetaBias = true; 132 ThetaBias = true;
133 } 133 }
134 134
135 void G4SPSRandomGenerator::SetPhiBias(const G4 135 void G4SPSRandomGenerator::SetPhiBias(const G4ThreeVector& input)
136 { 136 {
137 G4AutoLock l(&mutex); 137 G4AutoLock l(&mutex);
138 G4double ehi, val; 138 G4double ehi, val;
139 ehi = input.x(); 139 ehi = input.x();
140 val = input.y(); 140 val = input.y();
141 PhiBiasH.InsertValues(ehi, val); 141 PhiBiasH.InsertValues(ehi, val);
142 PhiBias = true; 142 PhiBias = true;
143 } 143 }
144 144
145 void G4SPSRandomGenerator::SetEnergyBias(const 145 void G4SPSRandomGenerator::SetEnergyBias(const G4ThreeVector& input)
146 { 146 {
147 G4AutoLock l(&mutex); 147 G4AutoLock l(&mutex);
148 G4double ehi, val; 148 G4double ehi, val;
149 ehi = input.x(); 149 ehi = input.x();
150 val = input.y(); 150 val = input.y();
151 EnergyBiasH.InsertValues(ehi, val); 151 EnergyBiasH.InsertValues(ehi, val);
152 EnergyBias = true; 152 EnergyBias = true;
153 } 153 }
154 154
155 void G4SPSRandomGenerator::SetPosThetaBias(con 155 void G4SPSRandomGenerator::SetPosThetaBias(const G4ThreeVector& input)
156 { 156 {
157 G4AutoLock l(&mutex); 157 G4AutoLock l(&mutex);
158 G4double ehi, val; 158 G4double ehi, val;
159 ehi = input.x(); 159 ehi = input.x();
160 val = input.y(); 160 val = input.y();
161 PosThetaBiasH.InsertValues(ehi, val); 161 PosThetaBiasH.InsertValues(ehi, val);
162 PosThetaBias = true; 162 PosThetaBias = true;
163 } 163 }
164 164
165 void G4SPSRandomGenerator::SetPosPhiBias(const 165 void G4SPSRandomGenerator::SetPosPhiBias(const G4ThreeVector& input)
166 { 166 {
167 G4AutoLock l(&mutex); 167 G4AutoLock l(&mutex);
168 G4double ehi, val; 168 G4double ehi, val;
169 ehi = input.x(); 169 ehi = input.x();
170 val = input.y(); 170 val = input.y();
171 PosPhiBiasH.InsertValues(ehi, val); 171 PosPhiBiasH.InsertValues(ehi, val);
172 PosPhiBias = true; 172 PosPhiBias = true;
173 } 173 }
174 174
175 void G4SPSRandomGenerator::SetIntensityWeight( 175 void G4SPSRandomGenerator::SetIntensityWeight(G4double weight)
176 { 176 {
177 bweights.Get()[8] = weight; 177 bweights.Get()[8] = weight;
178 } 178 }
179 179
180 G4double G4SPSRandomGenerator::GetBiasWeight() 180 G4double G4SPSRandomGenerator::GetBiasWeight() const
181 { 181 {
182 bweights_t& w = bweights.Get(); 182 bweights_t& w = bweights.Get();
183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[ 183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[5] * w[6] * w[7] * w[8];
184 } 184 }
185 185
186 void G4SPSRandomGenerator::SetVerbosity(G4int 186 void G4SPSRandomGenerator::SetVerbosity(G4int a)
187 { 187 {
188 G4AutoLock l(&mutex); 188 G4AutoLock l(&mutex);
189 verbosityLevel = a; 189 verbosityLevel = a;
190 } 190 }
191 191
192 namespace 192 namespace
193 { 193 {
194 G4PhysicsFreeVector ZeroPhysVector; // for r << 194 G4PhysicsOrderedFreeVector ZeroPhysVector; // for re-set only
195 } 195 }
196 196
197 void G4SPSRandomGenerator::ReSetHist(const G4S 197 void G4SPSRandomGenerator::ReSetHist(const G4String& atype)
198 { 198 {
199 G4AutoLock l(&mutex); 199 G4AutoLock l(&mutex);
200 if (atype == "biasx") { 200 if (atype == "biasx") {
201 XBias = false; 201 XBias = false;
202 IPDFXBias = false; 202 IPDFXBias = false;
203 local_IPDFXBias.Get().val = fa 203 local_IPDFXBias.Get().val = false;
204 XBiasH = IPDFXBiasH = ZeroPhys 204 XBiasH = IPDFXBiasH = ZeroPhysVector;
205 } else if (atype == "biasy") { 205 } else if (atype == "biasy") {
206 YBias = false; 206 YBias = false;
207 IPDFYBias = false; 207 IPDFYBias = false;
208 local_IPDFYBias.Get().val = fa 208 local_IPDFYBias.Get().val = false;
209 YBiasH = IPDFYBiasH = ZeroPhys 209 YBiasH = IPDFYBiasH = ZeroPhysVector;
210 } else if (atype == "biasz") { 210 } else if (atype == "biasz") {
211 ZBias = false; 211 ZBias = false;
212 IPDFZBias = false; 212 IPDFZBias = false;
213 local_IPDFZBias.Get().val = fa 213 local_IPDFZBias.Get().val = false;
214 ZBiasH = IPDFZBiasH = ZeroPhys 214 ZBiasH = IPDFZBiasH = ZeroPhysVector;
215 } else if (atype == "biast") { 215 } else if (atype == "biast") {
216 ThetaBias = false; 216 ThetaBias = false;
217 IPDFThetaBias = false; 217 IPDFThetaBias = false;
218 local_IPDFThetaBias.Get().val 218 local_IPDFThetaBias.Get().val = false;
219 ThetaBiasH = IPDFThetaBiasH = 219 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
220 } else if (atype == "biasp") { 220 } else if (atype == "biasp") {
221 PhiBias = false; 221 PhiBias = false;
222 IPDFPhiBias = false; 222 IPDFPhiBias = false;
223 local_IPDFPhiBias.Get().val = 223 local_IPDFPhiBias.Get().val = false;
224 PhiBiasH = IPDFPhiBiasH = Zero 224 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
225 } else if (atype == "biase") { 225 } else if (atype == "biase") {
226 EnergyBias = false; 226 EnergyBias = false;
227 IPDFEnergyBias = false; 227 IPDFEnergyBias = false;
228 local_IPDFEnergyBias.Get().val 228 local_IPDFEnergyBias.Get().val = false;
229 EnergyBiasH = IPDFEnergyBiasH 229 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
230 } else if (atype == "biaspt") { 230 } else if (atype == "biaspt") {
231 PosThetaBias = false; 231 PosThetaBias = false;
232 IPDFPosThetaBias = false; 232 IPDFPosThetaBias = false;
233 local_IPDFPosThetaBias.Get().v 233 local_IPDFPosThetaBias.Get().val = false;
234 PosThetaBiasH = IPDFPosThetaBi 234 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
235 } else if (atype == "biaspp") { 235 } else if (atype == "biaspp") {
236 PosPhiBias = false; 236 PosPhiBias = false;
237 IPDFPosPhiBias = false; 237 IPDFPosPhiBias = false;
238 local_IPDFPosPhiBias.Get().val 238 local_IPDFPosPhiBias.Get().val = false;
239 PosPhiBiasH = IPDFPosPhiBiasH 239 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
240 } else { 240 } else {
241 G4cout << "Error, histtype not 241 G4cout << "Error, histtype not accepted " << G4endl;
242 } 242 }
243 } 243 }
244 244
245 G4double G4SPSRandomGenerator::GenRandX() 245 G4double G4SPSRandomGenerator::GenRandX()
246 { 246 {
247 if (verbosityLevel >= 1) 247 if (verbosityLevel >= 1)
248 G4cout << "In GenRandX" << G4endl; 248 G4cout << "In GenRandX" << G4endl;
249 if (!XBias) << 249 if (XBias == false)
250 { 250 {
251 // X is not biased 251 // X is not biased
252 G4double rndm = G4UniformRand(); 252 G4double rndm = G4UniformRand();
253 return (rndm); 253 return (rndm);
254 } 254 }
255 << 255 else
256 // X is biased <<
257 // This is shared among threads, and we need <<
258 // only once. Multiple instances of this cla <<
259 // so we rely on a class-private, thread-pri <<
260 // to check if we need an initialiation. We <<
261 // because the Boolean on which we check is <<
262 // <<
263 if ( !local_IPDFXBias.Get().val ) <<
264 { 256 {
265 // For time that this thread arrived, here << 257 // X is biased
266 // Now two cases are possible: it is the f << 258 // This is shared among threads, and we need to initialize
267 // ANY thread has ever initialized this. << 259 // only once. Multiple instances of this class can exists
268 // Now we need a lock. In any case, the th << 260 // so we rely on a class-private, thread-private variable
269 // variable can now be set to true << 261 // to check if we need an initialiation. We do not use a lock here
>> 262 // because the Boolean on which we check is thread private
270 // 263 //
271 local_IPDFXBias.Get().val = true; << 264 if ( local_IPDFXBias.Get().val == false )
272 G4AutoLock l(&mutex); <<
273 if (!IPDFXBias) <<
274 { 265 {
275 // IPDF has not been created, so create << 266 // For time that this thread arrived, here
>> 267 // Now two cases are possible: it is the first time
>> 268 // ANY thread has ever initialized this.
>> 269 // Now we need a lock. In any case, the thread local
>> 270 // variable can now be set to true
276 // 271 //
277 G4double bins[1024], vals[1024], sum; << 272 local_IPDFXBias.Get().val = true;
278 std::size_t ii; << 273 G4AutoLock l(&mutex);
279 std::size_t maxbin = XBiasH.GetVectorLen << 274 if (IPDFXBias == false)
280 bins[0] = XBiasH.GetLowEdgeEnergy(0); <<
281 vals[0] = XBiasH(0); <<
282 sum = vals[0]; <<
283 for (ii=1; ii<maxbin; ++ii) <<
284 { 275 {
285 bins[ii] = XBiasH.GetLowEdgeEnergy(ii) << 276 // IPDF has not been created, so create it
286 vals[ii] = XBiasH(ii) + vals[ii - 1]; << 277 //
287 sum = sum + XBiasH(ii); << 278 G4double bins[1024], vals[1024], sum;
288 } << 279 G4int ii;
>> 280 G4int maxbin = G4int(XBiasH.GetVectorLength());
>> 281 bins[0] = XBiasH.GetLowEdgeEnergy(std::size_t(0));
>> 282 vals[0] = XBiasH(std::size_t(0));
>> 283 sum = vals[0];
>> 284 for (ii=1; ii<maxbin; ++ii)
>> 285 {
>> 286 bins[ii] = XBiasH.GetLowEdgeEnergy(std::size_t(ii));
>> 287 vals[ii] = XBiasH(std::size_t(ii)) + vals[ii - 1];
>> 288 sum = sum + XBiasH(std::size_t(ii));
>> 289 }
289 290
290 for (ii=0; ii<maxbin; ++ii) << 291 for (ii=0; ii<maxbin; ++ii)
291 { << 292 {
292 vals[ii] = vals[ii] / sum; << 293 vals[ii] = vals[ii] / sum;
293 IPDFXBiasH.InsertValues(bins[ii], vals << 294 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
>> 295 }
>> 296 IPDFXBias = true;
294 } 297 }
295 IPDFXBias = true; <<
296 } 298 }
297 } << 299
298 << 300 // IPDF has been create so carry on
299 // IPDF has been create so carry on <<
300 301
301 G4double rndm = G4UniformRand(); << 302 G4double rndm = G4UniformRand();
302 303
303 // Calculate the weighting: Find the bin tha << 304 // Calculate the weighting: Find the bin that the determined
304 // rndm is in and the weigthing will be the << 305 // rndm is in and the weigthing will be the difference in the
305 // natural probability (from the x-axis) div << 306 // natural probability (from the x-axis) divided by the
306 // difference in the biased probability (the << 307 // difference in the biased probability (the area)
307 // << 308 //
308 std::size_t numberOfBin = IPDFXBiasH.GetVect << 309 std::size_t numberOfBin = IPDFXBiasH.GetVectorLength();
309 std::size_t biasn1 = 0; << 310 G4int biasn1 = 0;
310 std::size_t biasn2 = numberOfBin / 2; << 311 G4int biasn2 = numberOfBin / 2;
311 std::size_t biasn3 = numberOfBin - 1; << 312 G4int biasn3 = numberOfBin - 1;
312 while (biasn1 != biasn3 - 1) << 313 while (biasn1 != biasn3 - 1)
313 { << 314 {
314 if (rndm > IPDFXBiasH(biasn2)) << 315 if (rndm > IPDFXBiasH(biasn2))
315 { biasn1 = biasn2; } << 316 { biasn1 = biasn2; }
316 else << 317 else
317 { biasn3 = biasn2; } << 318 { biasn3 = biasn2; }
318 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / << 319 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
319 } << 320 }
320 321
321 // Retrieve the areas and then the x-axis va << 322 // Retrieve the areas and then the x-axis values
322 // << 323 //
323 bweights_t& w = bweights.Get(); << 324 bweights_t& w = bweights.Get();
324 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn << 325 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
325 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnerg << 326 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
326 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnerg << 327 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
327 G4double NatProb = xaxisu - xaxisl; << 328 G4double NatProb = xaxisu - xaxisl;
328 w[0] = NatProb / w[0]; << 329 w[0] = NatProb / w[0];
329 if (verbosityLevel >= 1) << 330 if (verbosityLevel >= 1)
330 { << 331 {
331 G4cout << "X bin weight " << w[0] << " " < << 332 G4cout << "X bin weight " << w[0] << " " << rndm << G4endl;
>> 333 }
>> 334 return (IPDFXBiasH.GetEnergy(rndm));
332 } 335 }
333 return (IPDFXBiasH.GetEnergy(rndm)); <<
334 <<
335 } 336 }
336 337
337 G4double G4SPSRandomGenerator::GenRandY() 338 G4double G4SPSRandomGenerator::GenRandY()
338 { 339 {
339 if (verbosityLevel >= 1) 340 if (verbosityLevel >= 1)
340 { 341 {
341 G4cout << "In GenRandY" << G4endl; 342 G4cout << "In GenRandY" << G4endl;
342 } 343 }
343 344
344 if (!YBias) // Y is not biased << 345 if (YBias == false) // Y is not biased
345 { 346 {
346 G4double rndm = G4UniformRand(); 347 G4double rndm = G4UniformRand();
347 return (rndm); 348 return (rndm);
348 } 349 }
349 // Y is biased << 350 else // Y is biased
350 if ( !local_IPDFYBias.Get().val ) << 351 {
>> 352 if ( local_IPDFYBias.Get().val == false )
351 { 353 {
352 local_IPDFYBias.Get().val = true; 354 local_IPDFYBias.Get().val = true;
353 G4AutoLock l(&mutex); 355 G4AutoLock l(&mutex);
354 if (!IPDFYBias) << 356 if (IPDFYBias == false)
355 { 357 {
356 // IPDF has not been created, so creat 358 // IPDF has not been created, so create it
357 // 359 //
358 G4double bins[1024], vals[1024], sum; 360 G4double bins[1024], vals[1024], sum;
359 std::size_t ii; << 361 G4int ii;
360 std::size_t maxbin = YBiasH.GetVectorL << 362 G4int maxbin = G4int(YBiasH.GetVectorLength());
361 bins[0] = YBiasH.GetLowEdgeEnergy(0); << 363 bins[0] = YBiasH.GetLowEdgeEnergy(std::size_t(0));
362 vals[0] = YBiasH(0); << 364 vals[0] = YBiasH(std::size_t(0));
363 sum = vals[0]; 365 sum = vals[0];
364 for (ii=1; ii<maxbin; ++ii) 366 for (ii=1; ii<maxbin; ++ii)
365 { 367 {
366 bins[ii] = YBiasH.GetLowEdgeEnergy(i << 368 bins[ii] = YBiasH.GetLowEdgeEnergy(std::size_t(ii));
367 vals[ii] = YBiasH(ii) + vals[ii - 1] << 369 vals[ii] = YBiasH(std::size_t(ii)) + vals[ii - 1];
368 sum = sum + YBiasH(ii); << 370 sum = sum + YBiasH(std::size_t(ii));
369 } 371 }
370 372
371 for (ii=0; ii<maxbin; ++ii) 373 for (ii=0; ii<maxbin; ++ii)
372 { 374 {
373 vals[ii] = vals[ii] / sum; 375 vals[ii] = vals[ii] / sum;
374 IPDFYBiasH.InsertValues(bins[ii], va 376 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
375 } 377 }
376 IPDFYBias = true; 378 IPDFYBias = true;
377 } 379 }
378 } 380 }
379 381
380 // IPDF has been created so carry on 382 // IPDF has been created so carry on
381 383
382 G4double rndm = G4UniformRand(); 384 G4double rndm = G4UniformRand();
383 std::size_t numberOfBin = IPDFYBiasH.GetVe 385 std::size_t numberOfBin = IPDFYBiasH.GetVectorLength();
384 std::size_t biasn1 = 0; << 386 G4int biasn1 = 0;
385 std::size_t biasn2 = numberOfBin / 2; << 387 G4int biasn2 = numberOfBin / 2;
386 std::size_t biasn3 = numberOfBin - 1; << 388 G4int biasn3 = numberOfBin - 1;
387 while (biasn1 != biasn3 - 1) 389 while (biasn1 != biasn3 - 1)
388 { 390 {
389 if (rndm > IPDFYBiasH(biasn2)) 391 if (rndm > IPDFYBiasH(biasn2))
390 { biasn1 = biasn2; } 392 { biasn1 = biasn2; }
391 else 393 else
392 { biasn3 = biasn2; } 394 { biasn3 = biasn2; }
393 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 395 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
394 } 396 }
395 bweights_t& w = bweights.Get(); 397 bweights_t& w = bweights.Get();
396 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(bia 398 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
397 G4double xaxisl = IPDFYBiasH.GetLowEdgeEne << 399 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
398 G4double xaxisu = IPDFYBiasH.GetLowEdgeEne << 400 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
399 G4double NatProb = xaxisu - xaxisl; 401 G4double NatProb = xaxisu - xaxisl;
400 w[1] = NatProb / w[1]; 402 w[1] = NatProb / w[1];
401 if (verbosityLevel >= 1) 403 if (verbosityLevel >= 1)
402 { 404 {
403 G4cout << "Y bin weight " << w[1] << " " 405 G4cout << "Y bin weight " << w[1] << " " << rndm << G4endl;
404 } 406 }
405 return (IPDFYBiasH.GetEnergy(rndm)); 407 return (IPDFYBiasH.GetEnergy(rndm));
406 << 408 }
407 } 409 }
408 410
409 G4double G4SPSRandomGenerator::GenRandZ() 411 G4double G4SPSRandomGenerator::GenRandZ()
410 { 412 {
411 if (verbosityLevel >= 1) 413 if (verbosityLevel >= 1)
412 { 414 {
413 G4cout << "In GenRandZ" << G4endl; 415 G4cout << "In GenRandZ" << G4endl;
414 } 416 }
415 417
416 if (!ZBias) // Z is not biased << 418 if (ZBias == false) // Z is not biased
417 { 419 {
418 G4double rndm = G4UniformRand(); 420 G4double rndm = G4UniformRand();
419 return (rndm); 421 return (rndm);
420 } 422 }
421 // Z is biased << 423 else // Z is biased
422 if ( !local_IPDFZBias.Get().val ) << 424 {
>> 425 if ( local_IPDFZBias.Get().val == false )
423 { 426 {
424 local_IPDFZBias.Get().val = true; 427 local_IPDFZBias.Get().val = true;
425 G4AutoLock l(&mutex); 428 G4AutoLock l(&mutex);
426 if (!IPDFZBias) << 429 if (IPDFZBias == false)
427 { 430 {
428 // IPDF has not been created, so creat 431 // IPDF has not been created, so create it
429 // 432 //
430 G4double bins[1024], vals[1024], sum; 433 G4double bins[1024], vals[1024], sum;
431 std::size_t ii; << 434 G4int ii;
432 std::size_t maxbin = ZBiasH.GetVectorL << 435 G4int maxbin = G4int(ZBiasH.GetVectorLength());
433 bins[0] = ZBiasH.GetLowEdgeEnergy(0); << 436 bins[0] = ZBiasH.GetLowEdgeEnergy(std::size_t(0));
434 vals[0] = ZBiasH(0); << 437 vals[0] = ZBiasH(std::size_t(0));
435 sum = vals[0]; 438 sum = vals[0];
436 for (ii=1; ii<maxbin; ++ii) 439 for (ii=1; ii<maxbin; ++ii)
437 { 440 {
438 bins[ii] = ZBiasH.GetLowEdgeEnergy(i << 441 bins[ii] = ZBiasH.GetLowEdgeEnergy(std::size_t(ii));
439 vals[ii] = ZBiasH(ii) + vals[ii - 1] << 442 vals[ii] = ZBiasH(std::size_t(ii)) + vals[ii - 1];
440 sum = sum + ZBiasH(ii); << 443 sum = sum + ZBiasH(std::size_t(ii));
441 } 444 }
442 445
443 for (ii=0; ii<maxbin; ++ii) 446 for (ii=0; ii<maxbin; ++ii)
444 { 447 {
445 vals[ii] = vals[ii] / sum; 448 vals[ii] = vals[ii] / sum;
446 IPDFZBiasH.InsertValues(bins[ii], va 449 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
447 } 450 }
448 IPDFZBias = true; 451 IPDFZBias = true;
449 } 452 }
450 } 453 }
451 454
452 // IPDF has been create so carry on 455 // IPDF has been create so carry on
453 456
454 G4double rndm = G4UniformRand(); 457 G4double rndm = G4UniformRand();
455 std::size_t numberOfBin = IPDFZBiasH.GetVe 458 std::size_t numberOfBin = IPDFZBiasH.GetVectorLength();
456 std::size_t biasn1 = 0; << 459 G4int biasn1 = 0;
457 std::size_t biasn2 = numberOfBin / 2; << 460 G4int biasn2 = numberOfBin / 2;
458 std::size_t biasn3 = numberOfBin - 1; << 461 G4int biasn3 = numberOfBin - 1;
459 while (biasn1 != biasn3 - 1) 462 while (biasn1 != biasn3 - 1)
460 { 463 {
461 if (rndm > IPDFZBiasH(biasn2)) 464 if (rndm > IPDFZBiasH(biasn2))
462 { biasn1 = biasn2; } 465 { biasn1 = biasn2; }
463 else 466 else
464 { biasn3 = biasn2; } 467 { biasn3 = biasn2; }
465 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 468 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
466 } 469 }
467 bweights_t& w = bweights.Get(); 470 bweights_t& w = bweights.Get();
468 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(bia 471 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
469 G4double xaxisl = IPDFZBiasH.GetLowEdgeEne << 472 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
470 G4double xaxisu = IPDFZBiasH.GetLowEdgeEne << 473 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
471 G4double NatProb = xaxisu - xaxisl; 474 G4double NatProb = xaxisu - xaxisl;
472 w[2] = NatProb / w[2]; 475 w[2] = NatProb / w[2];
473 if (verbosityLevel >= 1) 476 if (verbosityLevel >= 1)
474 { 477 {
475 G4cout << "Z bin weight " << w[2] << " " 478 G4cout << "Z bin weight " << w[2] << " " << rndm << G4endl;
476 } 479 }
477 return (IPDFZBiasH.GetEnergy(rndm)); 480 return (IPDFZBiasH.GetEnergy(rndm));
478 << 481 }
479 } 482 }
480 483
481 G4double G4SPSRandomGenerator::GenRandTheta() 484 G4double G4SPSRandomGenerator::GenRandTheta()
482 { 485 {
483 if (verbosityLevel >= 1) 486 if (verbosityLevel >= 1)
484 { 487 {
485 G4cout << "In GenRandTheta" << G4endl; 488 G4cout << "In GenRandTheta" << G4endl;
486 G4cout << "Verbosity " << verbosityLevel < 489 G4cout << "Verbosity " << verbosityLevel << G4endl;
487 } 490 }
488 491
489 if (!ThetaBias) // Theta is not biased << 492 if (ThetaBias == false) // Theta is not biased
490 { 493 {
491 G4double rndm = G4UniformRand(); 494 G4double rndm = G4UniformRand();
492 return (rndm); 495 return (rndm);
493 } 496 }
494 // Theta is biased << 497 else // Theta is biased
495 if ( !local_IPDFThetaBias.Get().val ) << 498 {
>> 499 if ( local_IPDFThetaBias.Get().val == false )
496 { 500 {
497 local_IPDFThetaBias.Get().val = true; 501 local_IPDFThetaBias.Get().val = true;
498 G4AutoLock l(&mutex); 502 G4AutoLock l(&mutex);
499 if (!IPDFThetaBias) << 503 if (IPDFThetaBias == false)
500 { 504 {
501 // IPDF has not been created, so creat 505 // IPDF has not been created, so create it
502 // 506 //
503 G4double bins[1024], vals[1024], sum; 507 G4double bins[1024], vals[1024], sum;
504 std::size_t ii; << 508 G4int ii;
505 std::size_t maxbin = ThetaBiasH.GetVec << 509 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
506 bins[0] = ThetaBiasH.GetLowEdgeEnergy( << 510 bins[0] = ThetaBiasH.GetLowEdgeEnergy(std::size_t(0));
507 vals[0] = ThetaBiasH(0); << 511 vals[0] = ThetaBiasH(std::size_t(0));
508 sum = vals[0]; 512 sum = vals[0];
509 for (ii=1; ii<maxbin; ++ii) 513 for (ii=1; ii<maxbin; ++ii)
510 { 514 {
511 bins[ii] = ThetaBiasH.GetLowEdgeEner << 515 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(std::size_t(ii));
512 vals[ii] = ThetaBiasH(ii) + vals[ii << 516 vals[ii] = ThetaBiasH(std::size_t(ii)) + vals[ii - 1];
513 sum = sum + ThetaBiasH(ii); << 517 sum = sum + ThetaBiasH(std::size_t(ii));
514 } 518 }
515 519
516 for (ii=0; ii<maxbin; ++ii) 520 for (ii=0; ii<maxbin; ++ii)
517 { 521 {
518 vals[ii] = vals[ii] / sum; 522 vals[ii] = vals[ii] / sum;
519 IPDFThetaBiasH.InsertValues(bins[ii] 523 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
520 } 524 }
521 IPDFThetaBias = true; 525 IPDFThetaBias = true;
522 } 526 }
523 } 527 }
524 528
525 // IPDF has been create so carry on 529 // IPDF has been create so carry on
526 530
527 G4double rndm = G4UniformRand(); 531 G4double rndm = G4UniformRand();
528 std::size_t numberOfBin = IPDFThetaBiasH.G 532 std::size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
529 std::size_t biasn1 = 0; << 533 G4int biasn1 = 0;
530 std::size_t biasn2 = numberOfBin / 2; << 534 G4int biasn2 = numberOfBin / 2;
531 std::size_t biasn3 = numberOfBin - 1; << 535 G4int biasn3 = numberOfBin - 1;
532 while (biasn1 != biasn3 - 1) 536 while (biasn1 != biasn3 - 1)
533 { 537 {
534 if (rndm > IPDFThetaBiasH(biasn2)) 538 if (rndm > IPDFThetaBiasH(biasn2))
535 { biasn1 = biasn2; } 539 { biasn1 = biasn2; }
536 else 540 else
537 { biasn3 = biasn2; } 541 { biasn3 = biasn2; }
538 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 542 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
539 } 543 }
540 bweights_t& w = bweights.Get(); 544 bweights_t& w = bweights.Get();
541 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaB 545 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
542 G4double xaxisl = IPDFThetaBiasH.GetLowEdg << 546 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
543 G4double xaxisu = IPDFThetaBiasH.GetLowEdg << 547 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
544 G4double NatProb = xaxisu - xaxisl; 548 G4double NatProb = xaxisu - xaxisl;
545 w[3] = NatProb / w[3]; 549 w[3] = NatProb / w[3];
546 if (verbosityLevel >= 1) 550 if (verbosityLevel >= 1)
547 { 551 {
548 G4cout << "Theta bin weight " << w[3] << 552 G4cout << "Theta bin weight " << w[3] << " " << rndm << G4endl;
549 } 553 }
550 return (IPDFThetaBiasH.GetEnergy(rndm)); 554 return (IPDFThetaBiasH.GetEnergy(rndm));
551 << 555 }
552 } 556 }
553 557
554 G4double G4SPSRandomGenerator::GenRandPhi() 558 G4double G4SPSRandomGenerator::GenRandPhi()
555 { 559 {
556 if (verbosityLevel >= 1) 560 if (verbosityLevel >= 1)
557 { 561 {
558 G4cout << "In GenRandPhi" << G4endl; 562 G4cout << "In GenRandPhi" << G4endl;
559 } 563 }
560 564
561 if (!PhiBias) // Phi is not biased << 565 if (PhiBias == false) // Phi is not biased
562 { 566 {
563 G4double rndm = G4UniformRand(); 567 G4double rndm = G4UniformRand();
564 return (rndm); 568 return (rndm);
565 } 569 }
566 // Phi is biased << 570 else // Phi is biased
567 if ( !local_IPDFPhiBias.Get().val ) << 571 {
>> 572 if ( local_IPDFPhiBias.Get().val == false )
568 { 573 {
569 local_IPDFPhiBias.Get().val = true; 574 local_IPDFPhiBias.Get().val = true;
570 G4AutoLock l(&mutex); 575 G4AutoLock l(&mutex);
571 if (!IPDFPhiBias) << 576 if (IPDFPhiBias == false)
572 { 577 {
573 // IPDF has not been created, so creat 578 // IPDF has not been created, so create it
574 // 579 //
575 G4double bins[1024], vals[1024], sum; 580 G4double bins[1024], vals[1024], sum;
576 std::size_t ii; << 581 G4int ii;
577 std::size_t maxbin = PhiBiasH.GetVecto << 582 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
578 bins[0] = PhiBiasH.GetLowEdgeEnergy(0) << 583 bins[0] = PhiBiasH.GetLowEdgeEnergy(std::size_t(0));
579 vals[0] = PhiBiasH(0); << 584 vals[0] = PhiBiasH(std::size_t(0));
580 sum = vals[0]; 585 sum = vals[0];
581 for (ii=1; ii<maxbin; ++ii) 586 for (ii=1; ii<maxbin; ++ii)
582 { 587 {
583 bins[ii] = PhiBiasH.GetLowEdgeEnergy << 588 bins[ii] = PhiBiasH.GetLowEdgeEnergy(std::size_t(ii));
584 vals[ii] = PhiBiasH(ii) + vals[ii - << 589 vals[ii] = PhiBiasH(std::size_t(ii)) + vals[ii - 1];
585 sum = sum + PhiBiasH(ii); << 590 sum = sum + PhiBiasH(std::size_t(ii));
586 } 591 }
587 592
588 for (ii=0; ii<maxbin; ++ii) 593 for (ii=0; ii<maxbin; ++ii)
589 { 594 {
590 vals[ii] = vals[ii] / sum; 595 vals[ii] = vals[ii] / sum;
591 IPDFPhiBiasH.InsertValues(bins[ii], 596 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
592 } 597 }
593 IPDFPhiBias = true; 598 IPDFPhiBias = true;
594 } 599 }
595 } 600 }
596 601
597 // IPDF has been create so carry on 602 // IPDF has been create so carry on
598 603
599 G4double rndm = G4UniformRand(); 604 G4double rndm = G4UniformRand();
600 std::size_t numberOfBin = IPDFPhiBiasH.Get 605 std::size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
601 std::size_t biasn1 = 0; << 606 G4int biasn1 = 0;
602 std::size_t biasn2 = numberOfBin / 2; << 607 G4int biasn2 = numberOfBin / 2;
603 std::size_t biasn3 = numberOfBin - 1; << 608 G4int biasn3 = numberOfBin - 1;
604 while (biasn1 != biasn3 - 1) 609 while (biasn1 != biasn3 - 1)
605 { 610 {
606 if (rndm > IPDFPhiBiasH(biasn2)) 611 if (rndm > IPDFPhiBiasH(biasn2))
607 { biasn1 = biasn2; } 612 { biasn1 = biasn2; }
608 else 613 else
609 { biasn3 = biasn2; } 614 { biasn3 = biasn2; }
610 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 615 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
611 } 616 }
612 bweights_t& w = bweights.Get(); 617 bweights_t& w = bweights.Get();
613 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH 618 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
614 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeE << 619 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
615 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeE << 620 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
616 G4double NatProb = xaxisu - xaxisl; 621 G4double NatProb = xaxisu - xaxisl;
617 w[4] = NatProb / w[4]; 622 w[4] = NatProb / w[4];
618 if (verbosityLevel >= 1) 623 if (verbosityLevel >= 1)
619 { 624 {
620 G4cout << "Phi bin weight " << w[4] << " 625 G4cout << "Phi bin weight " << w[4] << " " << rndm << G4endl;
621 } 626 }
622 return (IPDFPhiBiasH.GetEnergy(rndm)); 627 return (IPDFPhiBiasH.GetEnergy(rndm));
623 << 628 }
624 } 629 }
625 630
626 G4double G4SPSRandomGenerator::GenRandEnergy() 631 G4double G4SPSRandomGenerator::GenRandEnergy()
627 { 632 {
628 if (verbosityLevel >= 1) 633 if (verbosityLevel >= 1)
629 { 634 {
630 G4cout << "In GenRandEnergy" << G4endl; 635 G4cout << "In GenRandEnergy" << G4endl;
631 } 636 }
632 637
633 if (!EnergyBias) // Energy is not biased << 638 if (EnergyBias == false) // Energy is not biased
634 { 639 {
635 G4double rndm = G4UniformRand(); 640 G4double rndm = G4UniformRand();
636 return (rndm); 641 return (rndm);
637 } 642 }
638 // Energy is biased << 643 else // Energy is biased
639 if ( !local_IPDFEnergyBias.Get().val ) << 644 {
>> 645 if ( local_IPDFEnergyBias.Get().val == false )
640 { 646 {
641 local_IPDFEnergyBias.Get().val = true; 647 local_IPDFEnergyBias.Get().val = true;
642 G4AutoLock l(&mutex); 648 G4AutoLock l(&mutex);
643 if (!IPDFEnergyBias) << 649 if (IPDFEnergyBias == false)
644 { 650 {
645 // IPDF has not been created, so creat 651 // IPDF has not been created, so create it
646 // 652 //
647 G4double bins[1024], vals[1024], sum; 653 G4double bins[1024], vals[1024], sum;
648 std::size_t ii; << 654 G4int ii;
649 std::size_t maxbin = EnergyBiasH.GetVe << 655 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
650 bins[0] = EnergyBiasH.GetLowEdgeEnergy << 656 bins[0] = EnergyBiasH.GetLowEdgeEnergy(std::size_t(0));
651 vals[0] = EnergyBiasH(0); << 657 vals[0] = EnergyBiasH(std::size_t(0));
652 sum = vals[0]; 658 sum = vals[0];
653 for (ii=1; ii<maxbin; ++ii) 659 for (ii=1; ii<maxbin; ++ii)
654 { 660 {
655 bins[ii] = EnergyBiasH.GetLowEdgeEne << 661 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(std::size_t(ii));
656 vals[ii] = EnergyBiasH(ii) + vals[ii << 662 vals[ii] = EnergyBiasH(std::size_t(ii)) + vals[ii - 1];
657 sum = sum + EnergyBiasH(ii); << 663 sum = sum + EnergyBiasH(std::size_t(ii));
658 } 664 }
659 IPDFEnergyBiasH = ZeroPhysVector; 665 IPDFEnergyBiasH = ZeroPhysVector;
660 for (ii=0; ii<maxbin; ++ii) 666 for (ii=0; ii<maxbin; ++ii)
661 { 667 {
662 vals[ii] = vals[ii] / sum; 668 vals[ii] = vals[ii] / sum;
663 IPDFEnergyBiasH.InsertValues(bins[ii 669 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
664 } 670 }
665 IPDFEnergyBias = true; 671 IPDFEnergyBias = true;
666 } 672 }
667 } 673 }
668 674
669 // IPDF has been create so carry on 675 // IPDF has been create so carry on
670 676
671 G4double rndm = G4UniformRand(); 677 G4double rndm = G4UniformRand();
672 std::size_t numberOfBin = IPDFEnergyBiasH. 678 std::size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
673 std::size_t biasn1 = 0; << 679 G4int biasn1 = 0;
674 std::size_t biasn2 = numberOfBin / 2; << 680 G4int biasn2 = numberOfBin / 2;
675 std::size_t biasn3 = numberOfBin - 1; << 681 G4int biasn3 = numberOfBin - 1;
676 while (biasn1 != biasn3 - 1) 682 while (biasn1 != biasn3 - 1)
677 { 683 {
678 if (rndm > IPDFEnergyBiasH(biasn2)) 684 if (rndm > IPDFEnergyBiasH(biasn2))
679 { biasn1 = biasn2; } 685 { biasn1 = biasn2; }
680 else 686 else
681 { biasn3 = biasn2; } 687 { biasn3 = biasn2; }
682 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 688 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
683 } 689 }
684 bweights_t& w = bweights.Get(); 690 bweights_t& w = bweights.Get();
685 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnerg 691 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
686 G4double xaxisl = IPDFEnergyBiasH.GetLowEd << 692 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
687 G4double xaxisu = IPDFEnergyBiasH.GetLowEd << 693 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
688 G4double NatProb = xaxisu - xaxisl; 694 G4double NatProb = xaxisu - xaxisl;
689 w[5] = NatProb / w[5]; 695 w[5] = NatProb / w[5];
690 if (verbosityLevel >= 1) 696 if (verbosityLevel >= 1)
691 { 697 {
692 G4cout << "Energy bin weight " << w[5] < 698 G4cout << "Energy bin weight " << w[5] << " " << rndm << G4endl;
693 } 699 }
694 return (IPDFEnergyBiasH.GetEnergy(rndm)); 700 return (IPDFEnergyBiasH.GetEnergy(rndm));
695 << 701 }
696 } 702 }
697 703
698 G4double G4SPSRandomGenerator::GenRandPosTheta 704 G4double G4SPSRandomGenerator::GenRandPosTheta()
699 { 705 {
700 if (verbosityLevel >= 1) 706 if (verbosityLevel >= 1)
701 { 707 {
702 G4cout << "In GenRandPosTheta" << G4endl; 708 G4cout << "In GenRandPosTheta" << G4endl;
703 G4cout << "Verbosity " << verbosityLevel < 709 G4cout << "Verbosity " << verbosityLevel << G4endl;
704 } 710 }
705 711
706 if (!PosThetaBias) // Theta is not biased << 712 if (PosThetaBias == false) // Theta is not biased
707 { 713 {
708 G4double rndm = G4UniformRand(); 714 G4double rndm = G4UniformRand();
709 return (rndm); 715 return (rndm);
710 } 716 }
711 // Theta is biased << 717 else // Theta is biased
712 if ( !local_IPDFPosThetaBias.Get().val ) << 718 {
>> 719 if ( local_IPDFPosThetaBias.Get().val == false )
713 { 720 {
714 local_IPDFPosThetaBias.Get().val = true; 721 local_IPDFPosThetaBias.Get().val = true;
715 G4AutoLock l(&mutex); 722 G4AutoLock l(&mutex);
716 if (!IPDFPosThetaBias) << 723 if (IPDFPosThetaBias == false)
717 { 724 {
718 // IPDF has not been created, so creat 725 // IPDF has not been created, so create it
719 // 726 //
720 G4double bins[1024], vals[1024], sum; 727 G4double bins[1024], vals[1024], sum;
721 std::size_t ii; << 728 G4int ii;
722 std::size_t maxbin = PosThetaBiasH.Get << 729 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
723 bins[0] = PosThetaBiasH.GetLowEdgeEner << 730 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(std::size_t(0));
724 vals[0] = PosThetaBiasH(0); << 731 vals[0] = PosThetaBiasH(std::size_t(0));
725 sum = vals[0]; 732 sum = vals[0];
726 for (ii=1; ii<maxbin; ++ii) 733 for (ii=1; ii<maxbin; ++ii)
727 { 734 {
728 bins[ii] = PosThetaBiasH.GetLowEdgeE << 735 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(std::size_t(ii));
729 vals[ii] = PosThetaBiasH(ii) + vals[ << 736 vals[ii] = PosThetaBiasH(std::size_t(ii)) + vals[ii - 1];
730 sum = sum + PosThetaBiasH(ii); << 737 sum = sum + PosThetaBiasH(std::size_t(ii));
731 } 738 }
732 739
733 for (ii=0; ii<maxbin; ++ii) 740 for (ii=0; ii<maxbin; ++ii)
734 { 741 {
735 vals[ii] = vals[ii] / sum; 742 vals[ii] = vals[ii] / sum;
736 IPDFPosThetaBiasH.InsertValues(bins[ 743 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
737 } 744 }
738 IPDFPosThetaBias = true; 745 IPDFPosThetaBias = true;
739 } 746 }
740 } 747 }
741 748
742 // IPDF has been create so carry on 749 // IPDF has been create so carry on
743 // 750 //
744 G4double rndm = G4UniformRand(); 751 G4double rndm = G4UniformRand();
745 std::size_t numberOfBin = IPDFPosThetaBias 752 std::size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
746 std::size_t biasn1 = 0; << 753 G4int biasn1 = 0;
747 std::size_t biasn2 = numberOfBin / 2; << 754 G4int biasn2 = numberOfBin / 2;
748 std::size_t biasn3 = numberOfBin - 1; << 755 G4int biasn3 = numberOfBin - 1;
749 while (biasn1 != biasn3 - 1) 756 while (biasn1 != biasn3 - 1)
750 { 757 {
751 if (rndm > IPDFPosThetaBiasH(biasn2)) 758 if (rndm > IPDFPosThetaBiasH(biasn2))
752 { biasn1 = biasn2; } 759 { biasn1 = biasn2; }
753 else 760 else
754 { biasn3 = biasn2; } 761 { biasn3 = biasn2; }
755 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 762 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
756 } 763 }
757 bweights_t& w = bweights.Get(); 764 bweights_t& w = bweights.Get();
758 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPos 765 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
759 G4double xaxisl = IPDFPosThetaBiasH.GetLow << 766 G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2-1));
760 G4double xaxisu = IPDFPosThetaBiasH.GetLow << 767 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
761 G4double NatProb = xaxisu - xaxisl; 768 G4double NatProb = xaxisu - xaxisl;
762 w[6] = NatProb / w[6]; 769 w[6] = NatProb / w[6];
763 if (verbosityLevel >= 1) 770 if (verbosityLevel >= 1)
764 { 771 {
765 G4cout << "PosTheta bin weight " << w[6] 772 G4cout << "PosTheta bin weight " << w[6] << " " << rndm << G4endl;
766 } 773 }
767 return (IPDFPosThetaBiasH.GetEnergy(rndm)) 774 return (IPDFPosThetaBiasH.GetEnergy(rndm));
768 << 775 }
769 } 776 }
770 777
771 G4double G4SPSRandomGenerator::GenRandPosPhi() 778 G4double G4SPSRandomGenerator::GenRandPosPhi()
772 { 779 {
773 if (verbosityLevel >= 1) 780 if (verbosityLevel >= 1)
774 { 781 {
775 G4cout << "In GenRandPosPhi" << G4endl; 782 G4cout << "In GenRandPosPhi" << G4endl;
776 } 783 }
777 784
778 if (!PosPhiBias) // PosPhi is not biased << 785 if (PosPhiBias == false) // PosPhi is not biased
779 { 786 {
780 G4double rndm = G4UniformRand(); 787 G4double rndm = G4UniformRand();
781 return (rndm); 788 return (rndm);
782 } 789 }
783 // PosPhi is biased << 790 else // PosPhi is biased
784 if (!local_IPDFPosPhiBias.Get().val ) << 791 {
>> 792 if (local_IPDFPosPhiBias.Get().val == false )
785 { 793 {
786 local_IPDFPosPhiBias.Get().val = true; 794 local_IPDFPosPhiBias.Get().val = true;
787 G4AutoLock l(&mutex); 795 G4AutoLock l(&mutex);
788 if (!IPDFPosPhiBias) << 796 if (IPDFPosPhiBias == false)
789 { 797 {
790 // IPDF has not been created, so creat 798 // IPDF has not been created, so create it
791 // 799 //
792 G4double bins[1024], vals[1024], sum; 800 G4double bins[1024], vals[1024], sum;
793 std::size_t ii; << 801 G4int ii;
794 std::size_t maxbin = PosPhiBiasH.GetVe << 802 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
795 bins[0] = PosPhiBiasH.GetLowEdgeEnergy << 803 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(std::size_t(0));
796 vals[0] = PosPhiBiasH(0); << 804 vals[0] = PosPhiBiasH(std::size_t(0));
797 sum = vals[0]; 805 sum = vals[0];
798 for (ii=1; ii<maxbin; ++ii) 806 for (ii=1; ii<maxbin; ++ii)
799 { 807 {
800 bins[ii] = PosPhiBiasH.GetLowEdgeEne << 808 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(std::size_t(ii));
801 vals[ii] = PosPhiBiasH(ii) + vals[ii << 809 vals[ii] = PosPhiBiasH(std::size_t(ii)) + vals[ii - 1];
802 sum = sum + PosPhiBiasH(ii); << 810 sum = sum + PosPhiBiasH(std::size_t(ii));
803 } 811 }
804 812
805 for (ii=0; ii<maxbin; ++ii) 813 for (ii=0; ii<maxbin; ++ii)
806 { 814 {
807 vals[ii] = vals[ii] / sum; 815 vals[ii] = vals[ii] / sum;
808 IPDFPosPhiBiasH.InsertValues(bins[ii 816 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
809 } 817 }
810 IPDFPosPhiBias = true; 818 IPDFPosPhiBias = true;
811 } 819 }
812 } 820 }
813 821
814 // IPDF has been create so carry on 822 // IPDF has been create so carry on
815 823
816 G4double rndm = G4UniformRand(); 824 G4double rndm = G4UniformRand();
817 std::size_t numberOfBin = IPDFPosPhiBiasH. 825 std::size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
818 std::size_t biasn1 = 0; << 826 G4int biasn1 = 0;
819 std::size_t biasn2 = numberOfBin / 2; << 827 G4int biasn2 = numberOfBin / 2;
820 std::size_t biasn3 = numberOfBin - 1; << 828 G4int biasn3 = numberOfBin - 1;
821 while (biasn1 != biasn3 - 1) 829 while (biasn1 != biasn3 - 1)
822 { 830 {
823 if (rndm > IPDFPosPhiBiasH(biasn2)) 831 if (rndm > IPDFPosPhiBiasH(biasn2))
824 { biasn1 = biasn2; } 832 { biasn1 = biasn2; }
825 else 833 else
826 { biasn3 = biasn2; } 834 { biasn3 = biasn2; }
827 biasn2 = biasn1 + (biasn3 - biasn1 + 1) 835 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
828 } 836 }
829 bweights_t& w = bweights.Get(); 837 bweights_t& w = bweights.Get();
830 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPh 838 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
831 G4double xaxisl = IPDFPosPhiBiasH.GetLowEd << 839 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
832 G4double xaxisu = IPDFPosPhiBiasH.GetLowEd << 840 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
833 G4double NatProb = xaxisu - xaxisl; 841 G4double NatProb = xaxisu - xaxisl;
834 w[7] = NatProb / w[7]; 842 w[7] = NatProb / w[7];
835 if (verbosityLevel >= 1) 843 if (verbosityLevel >= 1)
836 { 844 {
837 G4cout << "PosPhi bin weight " << w[7] < 845 G4cout << "PosPhi bin weight " << w[7] << " " << rndm << G4endl;
838 } 846 }
839 return (IPDFPosPhiBiasH.GetEnergy(rndm)); 847 return (IPDFPosPhiBiasH.GetEnergy(rndm));
840 << 848 }
841 } 849 }
842 850