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