source: trunk/source/event/src/G4SPSRandomGenerator.cc@ 1353

Last change on this file since 1353 was 1347, checked in by garnier, 15 years ago

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