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