Context Navigation

G4SPSRandomGenerator.cc @ 1347

Last change on this file since 1347 was 1347, checked in by garnier, 13 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
65	G4SPSRandomGenerator::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
90	G4SPSRandomGenerator::~G4SPSRandomGenerator() {
91	}
92
93	//G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
94	//{
95	// if (instance == 0) instance = new G4SPSRandomGenerator();
96	// return instance;
97	//}
98
99	// Biasing methods
100
101	void 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
109	void 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
117	void 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
125	void 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
133	void 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
141	void 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
149	void 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
157	void 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
165	void 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
203	G4double 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
265	G4double 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
319	G4double 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
374	G4double 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
432	G4double 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
487	G4double 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
543	G4double 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
602	G4double 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.

Context Navigation

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

Download in other formats: