| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4SPSRandomGenerator class implementation
27 //
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"
36 #include "G4Event.hh"
37 #include "Randomize.hh"
38 #include "G4TransportationManager.hh"
39 #include "G4VPhysicalVolume.hh"
40 #include "G4PhysicalVolumeStore.hh"
41 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh"
44 #include "G4Ions.hh"
45 #include "G4TrackingManager.hh"
46 #include "G4Track.hh"
47 #include "G4AutoLock.hh"
48
49 #include "G4SPSRandomGenerator.hh"
50
51 G4SPSRandomGenerator::bweights_t::bweights_t()
52 {
53 for (double & i : w) { i = 1; }
54 }
55
56 G4double& G4SPSRandomGenerator::bweights_t::operator [](const G4int i)
57 {
58 return w[i];
59 }
60
61 G4SPSRandomGenerator::G4SPSRandomGenerator()
62 {
63 // 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
84 verbosityLevel = 0;
85 G4MUTEXINIT(mutex);
86 }
87
88 G4SPSRandomGenerator::~G4SPSRandomGenerator()
89 {
90 G4MUTEXDESTROY(mutex);
91 }
92
93 // Biasing methods
94
95 void G4SPSRandomGenerator::SetXBias(const G4ThreeVector& input)
96 {
97 G4AutoLock l(&mutex);
98 G4double ehi, val;
99 ehi = input.x();
100 val = input.y();
101 XBiasH.InsertValues(ehi, val);
102 XBias = true;
103 }
104
105 void G4SPSRandomGenerator::SetYBias(const G4ThreeVector& input)
106 {
107 G4AutoLock l(&mutex);
108 G4double ehi, val;
109 ehi = input.x();
110 val = input.y();
111 YBiasH.InsertValues(ehi, val);
112 YBias = true;
113 }
114
115 void G4SPSRandomGenerator::SetZBias(const G4ThreeVector& input)
116 {
117 G4AutoLock l(&mutex);
118 G4double ehi, val;
119 ehi = input.x();
120 val = input.y();
121 ZBiasH.InsertValues(ehi, val);
122 ZBias = true;
123 }
124
125 void G4SPSRandomGenerator::SetThetaBias(const G4ThreeVector& input)
126 {
127 G4AutoLock l(&mutex);
128 G4double ehi, val;
129 ehi = input.x();
130 val = input.y();
131 ThetaBiasH.InsertValues(ehi, val);
132 ThetaBias = true;
133 }
134
135 void G4SPSRandomGenerator::SetPhiBias(const G4ThreeVector& input)
136 {
137 G4AutoLock l(&mutex);
138 G4double ehi, val;
139 ehi = input.x();
140 val = input.y();
141 PhiBiasH.InsertValues(ehi, val);
142 PhiBias = true;
143 }
144
145 void G4SPSRandomGenerator::SetEnergyBias(const G4ThreeVector& input)
146 {
147 G4AutoLock l(&mutex);
148 G4double ehi, val;
149 ehi = input.x();
150 val = input.y();
151 EnergyBiasH.InsertValues(ehi, val);
152 EnergyBias = true;
153 }
154
155 void G4SPSRandomGenerator::SetPosThetaBias(const G4ThreeVector& input)
156 {
157 G4AutoLock l(&mutex);
158 G4double ehi, val;
159 ehi = input.x();
160 val = input.y();
161 PosThetaBiasH.InsertValues(ehi, val);
162 PosThetaBias = true;
163 }
164
165 void G4SPSRandomGenerator::SetPosPhiBias(const G4ThreeVector& input)
166 {
167 G4AutoLock l(&mutex);
168 G4double ehi, val;
169 ehi = input.x();
170 val = input.y();
171 PosPhiBiasH.InsertValues(ehi, val);
172 PosPhiBias = true;
173 }
174
175 void G4SPSRandomGenerator::SetIntensityWeight(G4double weight)
176 {
177 bweights.Get()[8] = weight;
178 }
179
180 G4double G4SPSRandomGenerator::GetBiasWeight() const
181 {
182 bweights_t& w = bweights.Get();
183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[5] * w[6] * w[7] * w[8];
184 }
185
186 void G4SPSRandomGenerator::SetVerbosity(G4int a)
187 {
188 G4AutoLock l(&mutex);
189 verbosityLevel = a;
190 }
191
192 namespace
193 {
194 G4PhysicsFreeVector ZeroPhysVector; // for re-set only
195 }
196
197 void G4SPSRandomGenerator::ReSetHist(const G4String& atype)
198 {
199 G4AutoLock l(&mutex);
200 if (atype == "biasx") {
201 XBias = false;
202 IPDFXBias = false;
203 local_IPDFXBias.Get().val = false;
204 XBiasH = IPDFXBiasH = ZeroPhysVector;
205 } else if (atype == "biasy") {
206 YBias = false;
207 IPDFYBias = false;
208 local_IPDFYBias.Get().val = false;
209 YBiasH = IPDFYBiasH = ZeroPhysVector;
210 } else if (atype == "biasz") {
211 ZBias = false;
212 IPDFZBias = false;
213 local_IPDFZBias.Get().val = false;
214 ZBiasH = IPDFZBiasH = ZeroPhysVector;
215 } else if (atype == "biast") {
216 ThetaBias = false;
217 IPDFThetaBias = false;
218 local_IPDFThetaBias.Get().val = false;
219 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
220 } else if (atype == "biasp") {
221 PhiBias = false;
222 IPDFPhiBias = false;
223 local_IPDFPhiBias.Get().val = false;
224 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
225 } else if (atype == "biase") {
226 EnergyBias = false;
227 IPDFEnergyBias = false;
228 local_IPDFEnergyBias.Get().val = false;
229 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
230 } else if (atype == "biaspt") {
231 PosThetaBias = false;
232 IPDFPosThetaBias = false;
233 local_IPDFPosThetaBias.Get().val = false;
234 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
235 } else if (atype == "biaspp") {
236 PosPhiBias = false;
237 IPDFPosPhiBias = false;
238 local_IPDFPosPhiBias.Get().val = false;
239 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
240 } else {
241 G4cout << "Error, histtype not accepted " << G4endl;
242 }
243 }
244
245 G4double G4SPSRandomGenerator::GenRandX()
246 {
247 if (verbosityLevel >= 1)
248 G4cout << "In GenRandX" << G4endl;
249 if (!XBias)
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 to initialize
258 // only once. Multiple instances of this class can exists
259 // so we rely on a class-private, thread-private variable
260 // to check if we need an initialiation. We do not use a lock here
261 // because the Boolean on which we check is thread private
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 first time
267 // ANY thread has ever initialized this.
268 // Now we need a lock. In any case, the thread local
269 // variable can now be set to true
270 //
271 local_IPDFXBias.Get().val = true;
272 G4AutoLock l(&mutex);
273 if (!IPDFXBias)
274 {
275 // IPDF has not been created, so create it
276 //
277 G4double bins[1024], vals[1024], sum;
278 std::size_t ii;
279 std::size_t maxbin = XBiasH.GetVectorLength();
280 bins[0] = XBiasH.GetLowEdgeEnergy(0);
281 vals[0] = XBiasH(0);
282 sum = vals[0];
283 for (ii=1; ii<maxbin; ++ii)
284 {
285 bins[ii] = XBiasH.GetLowEdgeEnergy(ii);
286 vals[ii] = XBiasH(ii) + vals[ii - 1];
287 sum = sum + XBiasH(ii);
288 }
289
290 for (ii=0; ii<maxbin; ++ii)
291 {
292 vals[ii] = vals[ii] / sum;
293 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
294 }
295 IPDFXBias = true;
296 }
297 }
298
299 // IPDF has been create so carry on
300
301 G4double rndm = G4UniformRand();
302
303 // Calculate the weighting: Find the bin that the determined
304 // rndm is in and the weigthing will be the difference in the
305 // natural probability (from the x-axis) divided by the
306 // difference in the biased probability (the area)
307 //
308 std::size_t numberOfBin = IPDFXBiasH.GetVectorLength();
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) / 2;
319 }
320
321 // Retrieve the areas and then the x-axis values
322 //
323 bweights_t& w = bweights.Get();
324 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
325 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(biasn2 - 1);
326 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(biasn2);
327 G4double NatProb = xaxisu - xaxisl;
328 w[0] = NatProb / w[0];
329 if (verbosityLevel >= 1)
330 {
331 G4cout << "X bin weight " << w[0] << " " << rndm << G4endl;
332 }
333 return (IPDFXBiasH.GetEnergy(rndm));
334
335 }
336
337 G4double G4SPSRandomGenerator::GenRandY()
338 {
339 if (verbosityLevel >= 1)
340 {
341 G4cout << "In GenRandY" << G4endl;
342 }
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 create it
357 //
358 G4double bins[1024], vals[1024], sum;
359 std::size_t ii;
360 std::size_t maxbin = YBiasH.GetVectorLength();
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(ii);
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], vals[ii]);
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.GetVectorLength();
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 {
389 if (rndm > IPDFYBiasH(biasn2))
390 { biasn1 = biasn2; }
391 else
392 { biasn3 = biasn2; }
393 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
394 }
395 bweights_t& w = bweights.Get();
396 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
397 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(biasn2 - 1);
398 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(biasn2);
399 G4double NatProb = xaxisu - xaxisl;
400 w[1] = NatProb / w[1];
401 if (verbosityLevel >= 1)
402 {
403 G4cout << "Y bin weight " << w[1] << " " << rndm << G4endl;
404 }
405 return (IPDFYBiasH.GetEnergy(rndm));
406
407 }
408
409 G4double G4SPSRandomGenerator::GenRandZ()
410 {
411 if (verbosityLevel >= 1)
412 {
413 G4cout << "In GenRandZ" << G4endl;
414 }
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 create it
429 //
430 G4double bins[1024], vals[1024], sum;
431 std::size_t ii;
432 std::size_t maxbin = ZBiasH.GetVectorLength();
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(ii);
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], vals[ii]);
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.GetVectorLength();
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 {
461 if (rndm > IPDFZBiasH(biasn2))
462 { biasn1 = biasn2; }
463 else
464 { biasn3 = biasn2; }
465 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
466 }
467 bweights_t& w = bweights.Get();
468 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
469 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(biasn2 - 1);
470 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(biasn2);
471 G4double NatProb = xaxisu - xaxisl;
472 w[2] = NatProb / w[2];
473 if (verbosityLevel >= 1)
474 {
475 G4cout << "Z bin weight " << w[2] << " " << rndm << G4endl;
476 }
477 return (IPDFZBiasH.GetEnergy(rndm));
478
479 }
480
481 G4double G4SPSRandomGenerator::GenRandTheta()
482 {
483 if (verbosityLevel >= 1)
484 {
485 G4cout << "In GenRandTheta" << G4endl;
486 G4cout << "Verbosity " << verbosityLevel << G4endl;
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 {
497 local_IPDFThetaBias.Get().val = true;
498 G4AutoLock l(&mutex);
499 if (!IPDFThetaBias)
500 {
501 // IPDF has not been created, so create it
502 //
503 G4double bins[1024], vals[1024], sum;
504 std::size_t ii;
505 std::size_t maxbin = ThetaBiasH.GetVectorLength();
506 bins[0] = ThetaBiasH.GetLowEdgeEnergy(0);
507 vals[0] = ThetaBiasH(0);
508 sum = vals[0];
509 for (ii=1; ii<maxbin; ++ii)
510 {
511 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(ii);
512 vals[ii] = ThetaBiasH(ii) + vals[ii - 1];
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], vals[ii]);
520 }
521 IPDFThetaBias = true;
522 }
523 }
524
525 // IPDF has been create so carry on
526
527 G4double rndm = G4UniformRand();
528 std::size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
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 {
534 if (rndm > IPDFThetaBiasH(biasn2))
535 { biasn1 = biasn2; }
536 else
537 { biasn3 = biasn2; }
538 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
539 }
540 bweights_t& w = bweights.Get();
541 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
542 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(biasn2 - 1);
543 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(biasn2);
544 G4double NatProb = xaxisu - xaxisl;
545 w[3] = NatProb / w[3];
546 if (verbosityLevel >= 1)
547 {
548 G4cout << "Theta bin weight " << w[3] << " " << rndm << G4endl;
549 }
550 return (IPDFThetaBiasH.GetEnergy(rndm));
551
552 }
553
554 G4double G4SPSRandomGenerator::GenRandPhi()
555 {
556 if (verbosityLevel >= 1)
557 {
558 G4cout << "In GenRandPhi" << G4endl;
559 }
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 create it
574 //
575 G4double bins[1024], vals[1024], sum;
576 std::size_t ii;
577 std::size_t maxbin = PhiBiasH.GetVectorLength();
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(ii);
584 vals[ii] = PhiBiasH(ii) + vals[ii - 1];
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], vals[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.GetVectorLength();
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 {
606 if (rndm > IPDFPhiBiasH(biasn2))
607 { biasn1 = biasn2; }
608 else
609 { biasn3 = biasn2; }
610 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
611 }
612 bweights_t& w = bweights.Get();
613 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
614 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(biasn2 - 1);
615 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(biasn2);
616 G4double NatProb = xaxisu - xaxisl;
617 w[4] = NatProb / w[4];
618 if (verbosityLevel >= 1)
619 {
620 G4cout << "Phi bin weight " << w[4] << " " << rndm << G4endl;
621 }
622 return (IPDFPhiBiasH.GetEnergy(rndm));
623
624 }
625
626 G4double G4SPSRandomGenerator::GenRandEnergy()
627 {
628 if (verbosityLevel >= 1)
629 {
630 G4cout << "In GenRandEnergy" << G4endl;
631 }
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 {
641 local_IPDFEnergyBias.Get().val = true;
642 G4AutoLock l(&mutex);
643 if (!IPDFEnergyBias)
644 {
645 // IPDF has not been created, so create it
646 //
647 G4double bins[1024], vals[1024], sum;
648 std::size_t ii;
649 std::size_t maxbin = EnergyBiasH.GetVectorLength();
650 bins[0] = EnergyBiasH.GetLowEdgeEnergy(0);
651 vals[0] = EnergyBiasH(0);
652 sum = vals[0];
653 for (ii=1; ii<maxbin; ++ii)
654 {
655 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(ii);
656 vals[ii] = EnergyBiasH(ii) + vals[ii - 1];
657 sum = sum + EnergyBiasH(ii);
658 }
659 IPDFEnergyBiasH = ZeroPhysVector;
660 for (ii=0; ii<maxbin; ++ii)
661 {
662 vals[ii] = vals[ii] / sum;
663 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
664 }
665 IPDFEnergyBias = true;
666 }
667 }
668
669 // IPDF has been create so carry on
670
671 G4double rndm = G4UniformRand();
672 std::size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
673 std::size_t biasn1 = 0;
674 std::size_t biasn2 = numberOfBin / 2;
675 std::size_t biasn3 = numberOfBin - 1;
676 while (biasn1 != biasn3 - 1)
677 {
678 if (rndm > IPDFEnergyBiasH(biasn2))
679 { biasn1 = biasn2; }
680 else
681 { biasn3 = biasn2; }
682 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
683 }
684 bweights_t& w = bweights.Get();
685 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
686 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(biasn2 - 1);
687 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(biasn2);
688 G4double NatProb = xaxisu - xaxisl;
689 w[5] = NatProb / w[5];
690 if (verbosityLevel >= 1)
691 {
692 G4cout << "Energy bin weight " << w[5] << " " << rndm << G4endl;
693 }
694 return (IPDFEnergyBiasH.GetEnergy(rndm));
695
696 }
697
698 G4double G4SPSRandomGenerator::GenRandPosTheta()
699 {
700 if (verbosityLevel >= 1)
701 {
702 G4cout << "In GenRandPosTheta" << G4endl;
703 G4cout << "Verbosity " << verbosityLevel << G4endl;
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 {
714 local_IPDFPosThetaBias.Get().val = true;
715 G4AutoLock l(&mutex);
716 if (!IPDFPosThetaBias)
717 {
718 // IPDF has not been created, so create it
719 //
720 G4double bins[1024], vals[1024], sum;
721 std::size_t ii;
722 std::size_t maxbin = PosThetaBiasH.GetVectorLength();
723 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(0);
724 vals[0] = PosThetaBiasH(0);
725 sum = vals[0];
726 for (ii=1; ii<maxbin; ++ii)
727 {
728 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(ii);
729 vals[ii] = PosThetaBiasH(ii) + vals[ii - 1];
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[ii], vals[ii]);
737 }
738 IPDFPosThetaBias = true;
739 }
740 }
741
742 // IPDF has been create so carry on
743 //
744 G4double rndm = G4UniformRand();
745 std::size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
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 {
751 if (rndm > IPDFPosThetaBiasH(biasn2))
752 { biasn1 = biasn2; }
753 else
754 { biasn3 = biasn2; }
755 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
756 }
757 bweights_t& w = bweights.Get();
758 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
759 G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(biasn2-1);
760 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(biasn2);
761 G4double NatProb = xaxisu - xaxisl;
762 w[6] = NatProb / w[6];
763 if (verbosityLevel >= 1)
764 {
765 G4cout << "PosTheta bin weight " << w[6] << " " << rndm << G4endl;
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 create it
791 //
792 G4double bins[1024], vals[1024], sum;
793 std::size_t ii;
794 std::size_t maxbin = PosPhiBiasH.GetVectorLength();
795 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(0);
796 vals[0] = PosPhiBiasH(0);
797 sum = vals[0];
798 for (ii=1; ii<maxbin; ++ii)
799 {
800 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(ii);
801 vals[ii] = PosPhiBiasH(ii) + vals[ii - 1];
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], vals[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.GetVectorLength();
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) / 2;
828 }
829 bweights_t& w = bweights.Get();
830 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
831 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(biasn2 - 1);
832 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(biasn2);
833 G4double NatProb = xaxisu - xaxisl;
834 w[7] = NatProb / w[7];
835 if (verbosityLevel >= 1)
836 {
837 G4cout << "PosPhi bin weight " << w[7] << " " << rndm << G4endl;
838 }
839 return (IPDFPosPhiBiasH.GetEnergy(rndm));
840
841 }
842