Geant4 9.6.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///////////////////////////////////////////////////////////////////////////////
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