Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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///////////////////////////////////////////////////////////////////////////////
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///////////////////////////////////////////////////////////////////////////////
46//
47#include "G4PrimaryParticle.hh"
48#include "G4Event.hh"
49#include "Randomize.hh"
50#include <cmath>
51#include "G4TransportationManager.hh"
52#include "G4VPhysicalVolume.hh"
53#include "G4PhysicalVolumeStore.hh"
54#include "G4ParticleTable.hh"
55#include "G4ParticleDefinition.hh"
56#include "G4IonTable.hh"
57#include "G4Ions.hh"
58#include "G4TrackingManager.hh"
59#include "G4Track.hh"
60
61#include "G4SPSRandomGenerator.hh"
62
63//G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
64
65G4SPSRandomGenerator::G4SPSRandomGenerator()
66 : alpha(0.)
67{
68 // Initialise all variables
69
70 // Bias variables
71 XBias = false;
72 IPDFXBias = false;
73 YBias = false;
74 IPDFYBias = false;
75 ZBias = false;
76 IPDFZBias = false;
77 ThetaBias = false;
78 IPDFThetaBias = false;
79 PhiBias = false;
80 IPDFPhiBias = false;
81 EnergyBias = false;
82 IPDFEnergyBias = false;
83 PosThetaBias = false;
84 IPDFPosThetaBias = false;
85 PosPhiBias = false;
86 IPDFPosPhiBias = false;
87 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4]
88 = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
89 verbosityLevel = 0;
90}
91
92G4SPSRandomGenerator::~G4SPSRandomGenerator() {
93}
94
95//G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
96//{
97// if (instance == 0) instance = new G4SPSRandomGenerator();
98// return instance;
99//}
100
101// Biasing methods
102
103void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
104 G4double ehi, val;
105 ehi = input.x();
106 val = input.y();
107 XBiasH.InsertValues(ehi, val);
108 XBias = true;
109}
110
111void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
112 G4double ehi, val;
113 ehi = input.x();
114 val = input.y();
115 YBiasH.InsertValues(ehi, val);
116 YBias = true;
117}
118
119void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
120 G4double ehi, val;
121 ehi = input.x();
122 val = input.y();
123 ZBiasH.InsertValues(ehi, val);
124 ZBias = true;
125}
126
127void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
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(G4ThreeVector input) {
136 G4double ehi, val;
137 ehi = input.x();
138 val = input.y();
139 PhiBiasH.InsertValues(ehi, val);
140 PhiBias = true;
141}
142
143void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
144 G4double ehi, val;
145 ehi = input.x();
146 val = input.y();
147 EnergyBiasH.InsertValues(ehi, val);
148 EnergyBias = true;
149}
150
151void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
152 G4double ehi, val;
153 ehi = input.x();
154 val = input.y();
155 PosThetaBiasH.InsertValues(ehi, val);
156 PosThetaBias = true;
157}
158
159void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
160 G4double ehi, val;
161 ehi = input.x();
162 val = input.y();
163 PosPhiBiasH.InsertValues(ehi, val);
164 PosPhiBias = true;
165}
166
167void G4SPSRandomGenerator::ReSetHist(G4String atype) {
168 if (atype == "biasx") {
169 XBias = false;
170 IPDFXBias = false;
171 XBiasH = IPDFXBiasH = ZeroPhysVector;
172 } else if (atype == "biasy") {
173 YBias = false;
174 IPDFYBias = false;
175 YBiasH = IPDFYBiasH = ZeroPhysVector;
176 } else if (atype == "biasz") {
177 ZBias = false;
178 IPDFZBias = false;
179 ZBiasH = IPDFZBiasH = ZeroPhysVector;
180 } else if (atype == "biast") {
181 ThetaBias = false;
182 IPDFThetaBias = false;
183 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
184 } else if (atype == "biasp") {
185 PhiBias = false;
186 IPDFPhiBias = false;
187 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
188 } else if (atype == "biase") {
189 EnergyBias = false;
190 IPDFEnergyBias = false;
191 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
192 } else if (atype == "biaspt") {
193 PosThetaBias = false;
194 IPDFPosThetaBias = false;
195 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
196 } else if (atype == "biaspp") {
197 PosPhiBias = false;
198 IPDFPosPhiBias = false;
199 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
200 } else {
201 G4cout << "Error, histtype not accepted " << G4endl;
202 }
203}
204
205G4double G4SPSRandomGenerator::GenRandX() {
206 if (verbosityLevel >= 1)
207 G4cout << "In GenRandX" << G4endl;
208 if (XBias == false) {
209 // X is not biased
210 G4double rndm = G4UniformRand();
211 return (rndm);
212 } else {
213 // X is biased
214 if (IPDFXBias == false) {
215 // IPDF has not been created, so create it
216 G4double bins[1024], vals[1024], sum;
217 G4int ii;
218 G4int maxbin = G4int(XBiasH.GetVectorLength());
219 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
220 vals[0] = XBiasH(size_t(0));
221 sum = vals[0];
222 for (ii = 1; ii < maxbin; ii++) {
223 bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
224 vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
225 sum = sum + XBiasH(size_t(ii));
226 }
227
228 for (ii = 0; ii < maxbin; ii++) {
229 vals[ii] = vals[ii] / sum;
230 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
231 }
232 // Make IPDFXBias = true
233 IPDFXBias = true;
234 }
235 // IPDF has been create so carry on
236 G4double rndm = G4UniformRand();
237
238 // Calculate the weighting: Find the bin that the determined
239 // rndm is in and the weigthing will be the difference in the
240 // natural probability (from the x-axis) divided by the
241 // difference in the biased probability (the area).
242 size_t numberOfBin = IPDFXBiasH.GetVectorLength();
243 G4int biasn1 = 0;
244 G4int biasn2 = numberOfBin / 2;
245 G4int biasn3 = numberOfBin - 1;
246 while (biasn1 != biasn3 - 1) {
247 if (rndm > IPDFXBiasH(biasn2))
248 biasn1 = biasn2;
249 else
250 biasn3 = biasn2;
251 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
252 }
253 // retrieve the areas and then the x-axis values
254 bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
255 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
256 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
257 G4double NatProb = xaxisu - xaxisl;
258 //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
259 //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
260 bweights[0] = NatProb / bweights[0];
261 if (verbosityLevel >= 1)
262 G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
263 return (IPDFXBiasH.GetEnergy(rndm));
264 }
265}
266
267G4double G4SPSRandomGenerator::GenRandY() {
268 if (verbosityLevel >= 1)
269 G4cout << "In GenRandY" << G4endl;
270 if (YBias == false) {
271 // Y is not biased
272 G4double rndm = G4UniformRand();
273 return (rndm);
274 } else {
275 // Y is biased
276 if (IPDFYBias == false) {
277 // IPDF has not been created, so create it
278 G4double bins[1024], vals[1024], sum;
279 G4int ii;
280 G4int maxbin = G4int(YBiasH.GetVectorLength());
281 bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
282 vals[0] = YBiasH(size_t(0));
283 sum = vals[0];
284 for (ii = 1; ii < maxbin; ii++) {
285 bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
286 vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
287 sum = sum + YBiasH(size_t(ii));
288 }
289
290 for (ii = 0; ii < maxbin; ii++) {
291 vals[ii] = vals[ii] / sum;
292 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
293 }
294 // Make IPDFYBias = true
295 IPDFYBias = true;
296 }
297 // IPDF has been create so carry on
298 G4double rndm = G4UniformRand();
299 size_t numberOfBin = IPDFYBiasH.GetVectorLength();
300 G4int biasn1 = 0;
301 G4int biasn2 = numberOfBin / 2;
302 G4int biasn3 = numberOfBin - 1;
303 while (biasn1 != biasn3 - 1) {
304 if (rndm > IPDFYBiasH(biasn2))
305 biasn1 = biasn2;
306 else
307 biasn3 = biasn2;
308 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
309 }
310 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
311 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
312 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
313 G4double NatProb = xaxisu - xaxisl;
314 bweights[1] = NatProb / bweights[1];
315 if (verbosityLevel >= 1)
316 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
317 return (IPDFYBiasH.GetEnergy(rndm));
318 }
319}
320
321G4double G4SPSRandomGenerator::GenRandZ() {
322 if (verbosityLevel >= 1)
323 G4cout << "In GenRandZ" << G4endl;
324 if (ZBias == false) {
325 // Z is not biased
326 G4double rndm = G4UniformRand();
327 return (rndm);
328 } else {
329 // Z is biased
330 if (IPDFZBias == false) {
331 // IPDF has not been created, so create it
332 G4double bins[1024], vals[1024], sum;
333 G4int ii;
334 G4int maxbin = G4int(ZBiasH.GetVectorLength());
335 bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
336 vals[0] = ZBiasH(size_t(0));
337 sum = vals[0];
338 for (ii = 1; ii < maxbin; ii++) {
339 bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
340 vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
341 sum = sum + ZBiasH(size_t(ii));
342 }
343
344 for (ii = 0; ii < maxbin; ii++) {
345 vals[ii] = vals[ii] / sum;
346 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
347 }
348 // Make IPDFZBias = true
349 IPDFZBias = true;
350 }
351 // IPDF has been create so carry on
352 G4double rndm = G4UniformRand();
353 // size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
354 size_t numberOfBin = IPDFZBiasH.GetVectorLength();
355 G4int biasn1 = 0;
356 G4int biasn2 = numberOfBin / 2;
357 G4int biasn3 = numberOfBin - 1;
358 while (biasn1 != biasn3 - 1) {
359 if (rndm > IPDFZBiasH(biasn2))
360 biasn1 = biasn2;
361 else
362 biasn3 = biasn2;
363 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
364 }
365 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
366 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
367 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
368 G4double NatProb = xaxisu - xaxisl;
369 bweights[2] = NatProb / bweights[2];
370 if (verbosityLevel >= 1)
371 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
372 return (IPDFZBiasH.GetEnergy(rndm));
373 }
374}
375
376G4double G4SPSRandomGenerator::GenRandTheta() {
377 if (verbosityLevel >= 1) {
378 G4cout << "In GenRandTheta" << G4endl;
379 G4cout << "Verbosity " << verbosityLevel << G4endl;
380 }
381 if (ThetaBias == false) {
382 // Theta is not biased
383 G4double rndm = G4UniformRand();
384 return (rndm);
385 } else {
386 // Theta is biased
387 if (IPDFThetaBias == false) {
388 // IPDF has not been created, so create it
389 G4double bins[1024], vals[1024], sum;
390 G4int ii;
391 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
392 bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
393 vals[0] = ThetaBiasH(size_t(0));
394 sum = vals[0];
395 for (ii = 1; ii < maxbin; ii++) {
396 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
397 vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
398 sum = sum + ThetaBiasH(size_t(ii));
399 }
400
401 for (ii = 0; ii < maxbin; ii++) {
402 vals[ii] = vals[ii] / sum;
403 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
404 }
405 // Make IPDFThetaBias = true
406 IPDFThetaBias = true;
407 }
408 // IPDF has been create so carry on
409 G4double rndm = G4UniformRand();
410 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
411 size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
412 G4int biasn1 = 0;
413 G4int biasn2 = numberOfBin / 2;
414 G4int biasn3 = numberOfBin - 1;
415 while (biasn1 != biasn3 - 1) {
416 if (rndm > IPDFThetaBiasH(biasn2))
417 biasn1 = biasn2;
418 else
419 biasn3 = biasn2;
420 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
421 }
422 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
423 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
424 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
425 G4double NatProb = xaxisu - xaxisl;
426 bweights[3] = NatProb / bweights[3];
427 if (verbosityLevel >= 1)
428 G4cout << "Theta bin weight " << bweights[3] << " " << rndm
429 << G4endl;
430 return (IPDFThetaBiasH.GetEnergy(rndm));
431 }
432}
433
434G4double G4SPSRandomGenerator::GenRandPhi() {
435 if (verbosityLevel >= 1)
436 G4cout << "In GenRandPhi" << G4endl;
437 if (PhiBias == false) {
438 // Phi is not biased
439 G4double rndm = G4UniformRand();
440 return (rndm);
441 } else {
442 // Phi is biased
443 if (IPDFPhiBias == false) {
444 // IPDF has not been created, so create it
445 G4double bins[1024], vals[1024], sum;
446 G4int ii;
447 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
448 bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
449 vals[0] = PhiBiasH(size_t(0));
450 sum = vals[0];
451 for (ii = 1; ii < maxbin; ii++) {
452 bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
453 vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
454 sum = sum + PhiBiasH(size_t(ii));
455 }
456
457 for (ii = 0; ii < maxbin; ii++) {
458 vals[ii] = vals[ii] / sum;
459 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
460 }
461 // Make IPDFPhiBias = true
462 IPDFPhiBias = true;
463 }
464 // IPDF has been create so carry on
465 G4double rndm = G4UniformRand();
466 // size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
467 size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
468 G4int biasn1 = 0;
469 G4int biasn2 = numberOfBin / 2;
470 G4int biasn3 = numberOfBin - 1;
471 while (biasn1 != biasn3 - 1) {
472 if (rndm > IPDFPhiBiasH(biasn2))
473 biasn1 = biasn2;
474 else
475 biasn3 = biasn2;
476 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
477 }
478 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
479 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
480 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
481 G4double NatProb = xaxisu - xaxisl;
482 bweights[4] = NatProb / bweights[4];
483 if (verbosityLevel >= 1)
484 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
485 return (IPDFPhiBiasH.GetEnergy(rndm));
486 }
487}
488
489G4double G4SPSRandomGenerator::GenRandEnergy() {
490 if (verbosityLevel >= 1)
491 G4cout << "In GenRandEnergy" << G4endl;
492 if (EnergyBias == false) {
493 // Energy is not biased
494 G4double rndm = G4UniformRand();
495 return (rndm);
496 } else {
497 // ENERGY is biased
498 if (IPDFEnergyBias == false) {
499 // IPDF has not been created, so create it
500 G4double bins[1024], vals[1024], sum;
501 G4int ii;
502 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
503 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
504 vals[0] = EnergyBiasH(size_t(0));
505 sum = vals[0];
506 for (ii = 1; ii < maxbin; ii++) {
507 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
508 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
509 sum = sum + EnergyBiasH(size_t(ii));
510 }
511 IPDFEnergyBiasH = ZeroPhysVector;
512 for (ii = 0; ii < maxbin; ii++) {
513 vals[ii] = vals[ii] / sum;
514 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
515 }
516 // Make IPDFEnergyBias = true
517 IPDFEnergyBias = true;
518 }
519 // IPDF has been create so carry on
520 G4double rndm = G4UniformRand();
521 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
522 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
523 G4int biasn1 = 0;
524 G4int biasn2 = numberOfBin / 2;
525 G4int biasn3 = numberOfBin - 1;
526 while (biasn1 != biasn3 - 1) {
527 if (rndm > IPDFEnergyBiasH(biasn2))
528 biasn1 = biasn2;
529 else
530 biasn3 = biasn2;
531 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
532 }
533 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
534 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
535 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
536 G4double NatProb = xaxisu - xaxisl;
537 bweights[5] = NatProb / bweights[5];
538 if (verbosityLevel >= 1)
539 G4cout << "Energy bin weight " << bweights[5] << " " << rndm
540 << G4endl;
541 return (IPDFEnergyBiasH.GetEnergy(rndm));
542 }
543}
544
545G4double G4SPSRandomGenerator::GenRandPosTheta() {
546 if (verbosityLevel >= 1) {
547 G4cout << "In GenRandPosTheta" << G4endl;
548 G4cout << "Verbosity " << verbosityLevel << G4endl;
549 }
550 if (PosThetaBias == false) {
551 // Theta is not biased
552 G4double rndm = G4UniformRand();
553 return (rndm);
554 } else {
555 // Theta is biased
556 if (IPDFPosThetaBias == false) {
557 // IPDF has not been created, so create it
558 G4double bins[1024], vals[1024], sum;
559 G4int ii;
560 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
561 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
562 vals[0] = PosThetaBiasH(size_t(0));
563 sum = vals[0];
564 for (ii = 1; ii < maxbin; ii++) {
565 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
566 vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
567 sum = sum + PosThetaBiasH(size_t(ii));
568 }
569
570 for (ii = 0; ii < maxbin; ii++) {
571 vals[ii] = vals[ii] / sum;
572 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
573 }
574 // Make IPDFThetaBias = true
575 IPDFPosThetaBias = true;
576 }
577 // IPDF has been create so carry on
578 G4double rndm = G4UniformRand();
579 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
580 size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
581 G4int biasn1 = 0;
582 G4int biasn2 = numberOfBin / 2;
583 G4int biasn3 = numberOfBin - 1;
584 while (biasn1 != biasn3 - 1) {
585 if (rndm > IPDFPosThetaBiasH(biasn2))
586 biasn1 = biasn2;
587 else
588 biasn3 = biasn2;
589 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
590 }
591 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
592 G4double xaxisl =
593 IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
594 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
595 G4double NatProb = xaxisu - xaxisl;
596 bweights[6] = NatProb / bweights[6];
597 if (verbosityLevel >= 1)
598 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
599 << G4endl;
600 return (IPDFPosThetaBiasH.GetEnergy(rndm));
601 }
602}
603
604G4double G4SPSRandomGenerator::GenRandPosPhi() {
605 if (verbosityLevel >= 1)
606 G4cout << "In GenRandPosPhi" << G4endl;
607 if (PosPhiBias == false) {
608 // PosPhi is not biased
609 G4double rndm = G4UniformRand();
610 return (rndm);
611 } else {
612 // PosPhi is biased
613 if (IPDFPosPhiBias == false) {
614 // IPDF has not been created, so create it
615 G4double bins[1024], vals[1024], sum;
616 G4int ii;
617 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
618 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
619 vals[0] = PosPhiBiasH(size_t(0));
620 sum = vals[0];
621 for (ii = 1; ii < maxbin; ii++) {
622 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
623 vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
624 sum = sum + PosPhiBiasH(size_t(ii));
625 }
626
627 for (ii = 0; ii < maxbin; ii++) {
628 vals[ii] = vals[ii] / sum;
629 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
630 }
631 // Make IPDFPosPhiBias = true
632 IPDFPosPhiBias = true;
633 }
634 // IPDF has been create so carry on
635 G4double rndm = G4UniformRand();
636 // size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
637 size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
638 G4int biasn1 = 0;
639 G4int biasn2 = numberOfBin / 2;
640 G4int biasn3 = numberOfBin - 1;
641 while (biasn1 != biasn3 - 1) {
642 if (rndm > IPDFPosPhiBiasH(biasn2))
643 biasn1 = biasn2;
644 else
645 biasn3 = biasn2;
646 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
647 }
648 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
649 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
650 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
651 G4double NatProb = xaxisu - xaxisl;
652 bweights[7] = NatProb / bweights[7];
653 if (verbosityLevel >= 1)
654 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
655 << G4endl;
656 return (IPDFPosPhiBiasH.GetEnergy(rndm));
657 }
658}
659
G4double
double G4double
Definition: G4Types.hh:64
G4int
int G4int
Definition: G4Types.hh:66
G4endl
#define G4endl
Definition: G4ios.hh:52
G4cout
G4DLLIMPORT std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:53
CLHEP::Hep3Vector::x
double x() const
CLHEP::Hep3Vector::y
double y() const
G4PhysicsOrderedFreeVector::InsertValues
void InsertValues(G4double energy, G4double value)
Definition: G4PhysicsOrderedFreeVector.cc:90
G4PhysicsOrderedFreeVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsOrderedFreeVector.cc:107
G4PhysicsOrderedFreeVector::GetEnergy
G4double GetEnergy(G4double aValue)
Definition: G4PhysicsOrderedFreeVector.cc:112
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4SPSRandomGenerator::SetThetaBias
void SetThetaBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:127
G4SPSRandomGenerator::SetYBias
void SetYBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:111
G4SPSRandomGenerator::SetEnergyBias
void SetEnergyBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:143
G4SPSRandomGenerator::SetPhiBias
void SetPhiBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:135
G4SPSRandomGenerator::SetPosThetaBias
void SetPosThetaBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:151
G4SPSRandomGenerator::SetZBias
void SetZBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:119
G4SPSRandomGenerator::SetPosPhiBias
void SetPosPhiBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:159
G4SPSRandomGenerator::SetXBias
void SetXBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:103