Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4SPSRandomGenerator.cc
Go to the documentation of this file.
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
51G4SPSRandomGenerator::bweights_t::bweights_t()
52{
53 for ( G4int i=0 ; i<9 ; ++i ) { w[i] = 1; }
54}
55
56G4double& G4SPSRandomGenerator::bweights_t::operator [](const G4int i)
57{
58 return w[i];
59}
60
61G4SPSRandomGenerator::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
88G4SPSRandomGenerator::~G4SPSRandomGenerator()
89{
90 G4MUTEXDESTROY(mutex);
91}
92
93// Biasing methods
94
95void 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
105void 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
115void 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
125void 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
135void 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
145void 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
155void 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
165void 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
175void G4SPSRandomGenerator::SetIntensityWeight(G4double weight)
176{
177 bweights.Get()[8] = weight;
178}
179
180G4double 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
186void G4SPSRandomGenerator::SetVerbosity(G4int a)
187{
188 G4AutoLock l(&mutex);
189 verbosityLevel = a;
190}
191
192namespace
193{
194 G4PhysicsOrderedFreeVector ZeroPhysVector; // for re-set only
195}
196
197void 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
245G4double G4SPSRandomGenerator::GenRandX()
246{
247 if (verbosityLevel >= 1)
248 G4cout << "In GenRandX" << G4endl;
249 if (XBias == false)
250 {
251 // X is not biased
252 G4double rndm = G4UniformRand();
253 return (rndm);
254 }
255 else
256 {
257 // X is biased
258 // This is shared among threads, and we need to initialize
259 // only once. Multiple instances of this class can exists
260 // so we rely on a class-private, thread-private variable
261 // to check if we need an initialiation. We do not use a lock here
262 // because the Boolean on which we check is thread private
263 //
264 if ( local_IPDFXBias.Get().val == false )
265 {
266 // For time that this thread arrived, here
267 // Now two cases are possible: it is the first time
268 // ANY thread has ever initialized this.
269 // Now we need a lock. In any case, the thread local
270 // variable can now be set to true
271 //
272 local_IPDFXBias.Get().val = true;
273 G4AutoLock l(&mutex);
274 if (IPDFXBias == false)
275 {
276 // IPDF has not been created, so create it
277 //
278 G4double bins[1024], vals[1024], sum;
279 G4int ii;
280 G4int maxbin = G4int(XBiasH.GetVectorLength());
281 bins[0] = XBiasH.GetLowEdgeEnergy(std::size_t(0));
282 vals[0] = XBiasH(std::size_t(0));
283 sum = vals[0];
284 for (ii=1; ii<maxbin; ++ii)
285 {
286 bins[ii] = XBiasH.GetLowEdgeEnergy(std::size_t(ii));
287 vals[ii] = XBiasH(std::size_t(ii)) + vals[ii - 1];
288 sum = sum + XBiasH(std::size_t(ii));
289 }
290
291 for (ii=0; ii<maxbin; ++ii)
292 {
293 vals[ii] = vals[ii] / sum;
294 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
295 }
296 IPDFXBias = true;
297 }
298 }
299
300 // IPDF has been create so carry on
301
302 G4double rndm = G4UniformRand();
303
304 // Calculate the weighting: Find the bin that the determined
305 // rndm is in and the weigthing will be the difference in the
306 // natural probability (from the x-axis) divided by the
307 // difference in the biased probability (the area)
308 //
309 std::size_t numberOfBin = IPDFXBiasH.GetVectorLength();
310 G4int biasn1 = 0;
311 G4int biasn2 = numberOfBin / 2;
312 G4int biasn3 = numberOfBin - 1;
313 while (biasn1 != biasn3 - 1)
314 {
315 if (rndm > IPDFXBiasH(biasn2))
316 { biasn1 = biasn2; }
317 else
318 { biasn3 = biasn2; }
319 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
320 }
321
322 // Retrieve the areas and then the x-axis values
323 //
324 bweights_t& w = bweights.Get();
325 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
326 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
327 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
328 G4double NatProb = xaxisu - xaxisl;
329 w[0] = NatProb / w[0];
330 if (verbosityLevel >= 1)
331 {
332 G4cout << "X bin weight " << w[0] << " " << rndm << G4endl;
333 }
334 return (IPDFXBiasH.GetEnergy(rndm));
335 }
336}
337
338G4double G4SPSRandomGenerator::GenRandY()
339{
340 if (verbosityLevel >= 1)
341 {
342 G4cout << "In GenRandY" << G4endl;
343 }
344
345 if (YBias == false) // Y is not biased
346 {
347 G4double rndm = G4UniformRand();
348 return (rndm);
349 }
350 else // Y is biased
351 {
352 if ( local_IPDFYBias.Get().val == false )
353 {
354 local_IPDFYBias.Get().val = true;
355 G4AutoLock l(&mutex);
356 if (IPDFYBias == false)
357 {
358 // IPDF has not been created, so create it
359 //
360 G4double bins[1024], vals[1024], sum;
361 G4int ii;
362 G4int maxbin = G4int(YBiasH.GetVectorLength());
363 bins[0] = YBiasH.GetLowEdgeEnergy(std::size_t(0));
364 vals[0] = YBiasH(std::size_t(0));
365 sum = vals[0];
366 for (ii=1; ii<maxbin; ++ii)
367 {
368 bins[ii] = YBiasH.GetLowEdgeEnergy(std::size_t(ii));
369 vals[ii] = YBiasH(std::size_t(ii)) + vals[ii - 1];
370 sum = sum + YBiasH(std::size_t(ii));
371 }
372
373 for (ii=0; ii<maxbin; ++ii)
374 {
375 vals[ii] = vals[ii] / sum;
376 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
377 }
378 IPDFYBias = true;
379 }
380 }
381
382 // IPDF has been created so carry on
383
384 G4double rndm = G4UniformRand();
385 std::size_t numberOfBin = IPDFYBiasH.GetVectorLength();
386 G4int biasn1 = 0;
387 G4int biasn2 = numberOfBin / 2;
388 G4int biasn3 = numberOfBin - 1;
389 while (biasn1 != biasn3 - 1)
390 {
391 if (rndm > IPDFYBiasH(biasn2))
392 { biasn1 = biasn2; }
393 else
394 { biasn3 = biasn2; }
395 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
396 }
397 bweights_t& w = bweights.Get();
398 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
399 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
400 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
401 G4double NatProb = xaxisu - xaxisl;
402 w[1] = NatProb / w[1];
403 if (verbosityLevel >= 1)
404 {
405 G4cout << "Y bin weight " << w[1] << " " << rndm << G4endl;
406 }
407 return (IPDFYBiasH.GetEnergy(rndm));
408 }
409}
410
411G4double G4SPSRandomGenerator::GenRandZ()
412{
413 if (verbosityLevel >= 1)
414 {
415 G4cout << "In GenRandZ" << G4endl;
416 }
417
418 if (ZBias == false) // Z is not biased
419 {
420 G4double rndm = G4UniformRand();
421 return (rndm);
422 }
423 else // Z is biased
424 {
425 if ( local_IPDFZBias.Get().val == false )
426 {
427 local_IPDFZBias.Get().val = true;
428 G4AutoLock l(&mutex);
429 if (IPDFZBias == false)
430 {
431 // IPDF has not been created, so create it
432 //
433 G4double bins[1024], vals[1024], sum;
434 G4int ii;
435 G4int maxbin = G4int(ZBiasH.GetVectorLength());
436 bins[0] = ZBiasH.GetLowEdgeEnergy(std::size_t(0));
437 vals[0] = ZBiasH(std::size_t(0));
438 sum = vals[0];
439 for (ii=1; ii<maxbin; ++ii)
440 {
441 bins[ii] = ZBiasH.GetLowEdgeEnergy(std::size_t(ii));
442 vals[ii] = ZBiasH(std::size_t(ii)) + vals[ii - 1];
443 sum = sum + ZBiasH(std::size_t(ii));
444 }
445
446 for (ii=0; ii<maxbin; ++ii)
447 {
448 vals[ii] = vals[ii] / sum;
449 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
450 }
451 IPDFZBias = true;
452 }
453 }
454
455 // IPDF has been create so carry on
456
457 G4double rndm = G4UniformRand();
458 std::size_t numberOfBin = IPDFZBiasH.GetVectorLength();
459 G4int biasn1 = 0;
460 G4int biasn2 = numberOfBin / 2;
461 G4int biasn3 = numberOfBin - 1;
462 while (biasn1 != biasn3 - 1)
463 {
464 if (rndm > IPDFZBiasH(biasn2))
465 { biasn1 = biasn2; }
466 else
467 { biasn3 = biasn2; }
468 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
469 }
470 bweights_t& w = bweights.Get();
471 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
472 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
473 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
474 G4double NatProb = xaxisu - xaxisl;
475 w[2] = NatProb / w[2];
476 if (verbosityLevel >= 1)
477 {
478 G4cout << "Z bin weight " << w[2] << " " << rndm << G4endl;
479 }
480 return (IPDFZBiasH.GetEnergy(rndm));
481 }
482}
483
484G4double G4SPSRandomGenerator::GenRandTheta()
485{
486 if (verbosityLevel >= 1)
487 {
488 G4cout << "In GenRandTheta" << G4endl;
489 G4cout << "Verbosity " << verbosityLevel << G4endl;
490 }
491
492 if (ThetaBias == false) // Theta is not biased
493 {
494 G4double rndm = G4UniformRand();
495 return (rndm);
496 }
497 else // Theta is biased
498 {
499 if ( local_IPDFThetaBias.Get().val == false )
500 {
501 local_IPDFThetaBias.Get().val = true;
502 G4AutoLock l(&mutex);
503 if (IPDFThetaBias == false)
504 {
505 // IPDF has not been created, so create it
506 //
507 G4double bins[1024], vals[1024], sum;
508 G4int ii;
509 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
510 bins[0] = ThetaBiasH.GetLowEdgeEnergy(std::size_t(0));
511 vals[0] = ThetaBiasH(std::size_t(0));
512 sum = vals[0];
513 for (ii=1; ii<maxbin; ++ii)
514 {
515 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(std::size_t(ii));
516 vals[ii] = ThetaBiasH(std::size_t(ii)) + vals[ii - 1];
517 sum = sum + ThetaBiasH(std::size_t(ii));
518 }
519
520 for (ii=0; ii<maxbin; ++ii)
521 {
522 vals[ii] = vals[ii] / sum;
523 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
524 }
525 IPDFThetaBias = true;
526 }
527 }
528
529 // IPDF has been create so carry on
530
531 G4double rndm = G4UniformRand();
532 std::size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
533 G4int biasn1 = 0;
534 G4int biasn2 = numberOfBin / 2;
535 G4int biasn3 = numberOfBin - 1;
536 while (biasn1 != biasn3 - 1)
537 {
538 if (rndm > IPDFThetaBiasH(biasn2))
539 { biasn1 = biasn2; }
540 else
541 { biasn3 = biasn2; }
542 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
543 }
544 bweights_t& w = bweights.Get();
545 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
546 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
547 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
548 G4double NatProb = xaxisu - xaxisl;
549 w[3] = NatProb / w[3];
550 if (verbosityLevel >= 1)
551 {
552 G4cout << "Theta bin weight " << w[3] << " " << rndm << G4endl;
553 }
554 return (IPDFThetaBiasH.GetEnergy(rndm));
555 }
556}
557
558G4double G4SPSRandomGenerator::GenRandPhi()
559{
560 if (verbosityLevel >= 1)
561 {
562 G4cout << "In GenRandPhi" << G4endl;
563 }
564
565 if (PhiBias == false) // Phi is not biased
566 {
567 G4double rndm = G4UniformRand();
568 return (rndm);
569 }
570 else // Phi is biased
571 {
572 if ( local_IPDFPhiBias.Get().val == false )
573 {
574 local_IPDFPhiBias.Get().val = true;
575 G4AutoLock l(&mutex);
576 if (IPDFPhiBias == false)
577 {
578 // IPDF has not been created, so create it
579 //
580 G4double bins[1024], vals[1024], sum;
581 G4int ii;
582 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
583 bins[0] = PhiBiasH.GetLowEdgeEnergy(std::size_t(0));
584 vals[0] = PhiBiasH(std::size_t(0));
585 sum = vals[0];
586 for (ii=1; ii<maxbin; ++ii)
587 {
588 bins[ii] = PhiBiasH.GetLowEdgeEnergy(std::size_t(ii));
589 vals[ii] = PhiBiasH(std::size_t(ii)) + vals[ii - 1];
590 sum = sum + PhiBiasH(std::size_t(ii));
591 }
592
593 for (ii=0; ii<maxbin; ++ii)
594 {
595 vals[ii] = vals[ii] / sum;
596 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
597 }
598 IPDFPhiBias = true;
599 }
600 }
601
602 // IPDF has been create so carry on
603
604 G4double rndm = G4UniformRand();
605 std::size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
606 G4int biasn1 = 0;
607 G4int biasn2 = numberOfBin / 2;
608 G4int biasn3 = numberOfBin - 1;
609 while (biasn1 != biasn3 - 1)
610 {
611 if (rndm > IPDFPhiBiasH(biasn2))
612 { biasn1 = biasn2; }
613 else
614 { biasn3 = biasn2; }
615 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
616 }
617 bweights_t& w = bweights.Get();
618 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
619 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
620 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
621 G4double NatProb = xaxisu - xaxisl;
622 w[4] = NatProb / w[4];
623 if (verbosityLevel >= 1)
624 {
625 G4cout << "Phi bin weight " << w[4] << " " << rndm << G4endl;
626 }
627 return (IPDFPhiBiasH.GetEnergy(rndm));
628 }
629}
630
631G4double G4SPSRandomGenerator::GenRandEnergy()
632{
633 if (verbosityLevel >= 1)
634 {
635 G4cout << "In GenRandEnergy" << G4endl;
636 }
637
638 if (EnergyBias == false) // Energy is not biased
639 {
640 G4double rndm = G4UniformRand();
641 return (rndm);
642 }
643 else // Energy is biased
644 {
645 if ( local_IPDFEnergyBias.Get().val == false )
646 {
647 local_IPDFEnergyBias.Get().val = true;
648 G4AutoLock l(&mutex);
649 if (IPDFEnergyBias == false)
650 {
651 // IPDF has not been created, so create it
652 //
653 G4double bins[1024], vals[1024], sum;
654 G4int ii;
655 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
656 bins[0] = EnergyBiasH.GetLowEdgeEnergy(std::size_t(0));
657 vals[0] = EnergyBiasH(std::size_t(0));
658 sum = vals[0];
659 for (ii=1; ii<maxbin; ++ii)
660 {
661 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(std::size_t(ii));
662 vals[ii] = EnergyBiasH(std::size_t(ii)) + vals[ii - 1];
663 sum = sum + EnergyBiasH(std::size_t(ii));
664 }
665 IPDFEnergyBiasH = ZeroPhysVector;
666 for (ii=0; ii<maxbin; ++ii)
667 {
668 vals[ii] = vals[ii] / sum;
669 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
670 }
671 IPDFEnergyBias = true;
672 }
673 }
674
675 // IPDF has been create so carry on
676
677 G4double rndm = G4UniformRand();
678 std::size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
679 G4int biasn1 = 0;
680 G4int biasn2 = numberOfBin / 2;
681 G4int biasn3 = numberOfBin - 1;
682 while (biasn1 != biasn3 - 1)
683 {
684 if (rndm > IPDFEnergyBiasH(biasn2))
685 { biasn1 = biasn2; }
686 else
687 { biasn3 = biasn2; }
688 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
689 }
690 bweights_t& w = bweights.Get();
691 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
692 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
693 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
694 G4double NatProb = xaxisu - xaxisl;
695 w[5] = NatProb / w[5];
696 if (verbosityLevel >= 1)
697 {
698 G4cout << "Energy bin weight " << w[5] << " " << rndm << G4endl;
699 }
700 return (IPDFEnergyBiasH.GetEnergy(rndm));
701 }
702}
703
704G4double G4SPSRandomGenerator::GenRandPosTheta()
705{
706 if (verbosityLevel >= 1)
707 {
708 G4cout << "In GenRandPosTheta" << G4endl;
709 G4cout << "Verbosity " << verbosityLevel << G4endl;
710 }
711
712 if (PosThetaBias == false) // Theta is not biased
713 {
714 G4double rndm = G4UniformRand();
715 return (rndm);
716 }
717 else // Theta is biased
718 {
719 if ( local_IPDFPosThetaBias.Get().val == false )
720 {
721 local_IPDFPosThetaBias.Get().val = true;
722 G4AutoLock l(&mutex);
723 if (IPDFPosThetaBias == false)
724 {
725 // IPDF has not been created, so create it
726 //
727 G4double bins[1024], vals[1024], sum;
728 G4int ii;
729 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
730 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(std::size_t(0));
731 vals[0] = PosThetaBiasH(std::size_t(0));
732 sum = vals[0];
733 for (ii=1; ii<maxbin; ++ii)
734 {
735 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(std::size_t(ii));
736 vals[ii] = PosThetaBiasH(std::size_t(ii)) + vals[ii - 1];
737 sum = sum + PosThetaBiasH(std::size_t(ii));
738 }
739
740 for (ii=0; ii<maxbin; ++ii)
741 {
742 vals[ii] = vals[ii] / sum;
743 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
744 }
745 IPDFPosThetaBias = true;
746 }
747 }
748
749 // IPDF has been create so carry on
750 //
751 G4double rndm = G4UniformRand();
752 std::size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
753 G4int biasn1 = 0;
754 G4int biasn2 = numberOfBin / 2;
755 G4int biasn3 = numberOfBin - 1;
756 while (biasn1 != biasn3 - 1)
757 {
758 if (rndm > IPDFPosThetaBiasH(biasn2))
759 { biasn1 = biasn2; }
760 else
761 { biasn3 = biasn2; }
762 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
763 }
764 bweights_t& w = bweights.Get();
765 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
766 G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2-1));
767 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
768 G4double NatProb = xaxisu - xaxisl;
769 w[6] = NatProb / w[6];
770 if (verbosityLevel >= 1)
771 {
772 G4cout << "PosTheta bin weight " << w[6] << " " << rndm << G4endl;
773 }
774 return (IPDFPosThetaBiasH.GetEnergy(rndm));
775 }
776}
777
778G4double G4SPSRandomGenerator::GenRandPosPhi()
779{
780 if (verbosityLevel >= 1)
781 {
782 G4cout << "In GenRandPosPhi" << G4endl;
783 }
784
785 if (PosPhiBias == false) // PosPhi is not biased
786 {
787 G4double rndm = G4UniformRand();
788 return (rndm);
789 }
790 else // PosPhi is biased
791 {
792 if (local_IPDFPosPhiBias.Get().val == false )
793 {
794 local_IPDFPosPhiBias.Get().val = true;
795 G4AutoLock l(&mutex);
796 if (IPDFPosPhiBias == false)
797 {
798 // IPDF has not been created, so create it
799 //
800 G4double bins[1024], vals[1024], sum;
801 G4int ii;
802 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
803 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(std::size_t(0));
804 vals[0] = PosPhiBiasH(std::size_t(0));
805 sum = vals[0];
806 for (ii=1; ii<maxbin; ++ii)
807 {
808 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(std::size_t(ii));
809 vals[ii] = PosPhiBiasH(std::size_t(ii)) + vals[ii - 1];
810 sum = sum + PosPhiBiasH(std::size_t(ii));
811 }
812
813 for (ii=0; ii<maxbin; ++ii)
814 {
815 vals[ii] = vals[ii] / sum;
816 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
817 }
818 IPDFPosPhiBias = true;
819 }
820 }
821
822 // IPDF has been create so carry on
823
824 G4double rndm = G4UniformRand();
825 std::size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
826 G4int biasn1 = 0;
827 G4int biasn2 = numberOfBin / 2;
828 G4int biasn3 = numberOfBin - 1;
829 while (biasn1 != biasn3 - 1)
830 {
831 if (rndm > IPDFPosPhiBiasH(biasn2))
832 { biasn1 = biasn2; }
833 else
834 { biasn3 = biasn2; }
835 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
836 }
837 bweights_t& w = bweights.Get();
838 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
839 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2 - 1));
840 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(std::size_t(biasn2));
841 G4double NatProb = xaxisu - xaxisl;
842 w[7] = NatProb / w[7];
843 if (verbosityLevel >= 1)
844 {
845 G4cout << "PosPhi bin weight " << w[7] << " " << rndm << G4endl;
846 }
847 return (IPDFPosPhiBiasH.GetEnergy(rndm));
848 }
849}
G4MUTEXDESTROY
#define G4MUTEXDESTROY(mutex)
Definition: G4Threading.hh:90
G4MUTEXINIT
#define G4MUTEXINIT(mutex)
Definition: G4Threading.hh:87
G4double
double G4double
Definition: G4Types.hh:83
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4Cache::Get
value_type & Get() const
Definition: G4Cache.hh:315
G4PhysicsOrderedFreeVector::InsertValues
void InsertValues(G4double energy, G4double value)
Definition: G4PhysicsOrderedFreeVector.cc:88
G4PhysicsOrderedFreeVector::GetEnergy
G4double GetEnergy(G4double aValue)
Definition: G4PhysicsOrderedFreeVector.cc:105
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(std::size_t binNumber) const
Definition: G4PhysicsVector.cc:136
G4PhysicsVector::GetVectorLength
std::size_t GetVectorLength() const
G4SPSRandomGenerator::SetIntensityWeight
void SetIntensityWeight(G4double weight)
Definition: G4SPSRandomGenerator.cc:175
G4SPSRandomGenerator::SetXBias
void SetXBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:95
G4SPSRandomGenerator::SetEnergyBias
void SetEnergyBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:145
G4SPSRandomGenerator::SetPosPhiBias
void SetPosPhiBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:165
G4SPSRandomGenerator::SetThetaBias
void SetThetaBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:125
G4SPSRandomGenerator::GetBiasWeight
G4double GetBiasWeight() const
Definition: G4SPSRandomGenerator.cc:180
G4SPSRandomGenerator::SetYBias
void SetYBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:105
G4SPSRandomGenerator::SetPosThetaBias
void SetPosThetaBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:155
G4SPSRandomGenerator::SetPhiBias
void SetPhiBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:135
G4SPSRandomGenerator::SetZBias
void SetZBias(const G4ThreeVector &)
Definition: G4SPSRandomGenerator.cc:115
G4SPSRandomGenerator::ReSetHist
void ReSetHist(const G4String &)
Definition: G4SPSRandomGenerator.cc:197