Context Navigation

G4SPSRandomGenerator.cc@ 1355

Last change on this file since 1355 was 1347, checked in by garnier, 15 years ago
geant4 tag 9.4
File size: 20.4 KB

Rev	Line
[816]	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
[1347]	65	G4SPSRandomGenerator::G4SPSRandomGenerator() {
	66	// Initialise all variables
[816]	67
[1347]	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;
[816]	88	}
	89
[1347]	90	G4SPSRandomGenerator::~G4SPSRandomGenerator() {
	91	}
[816]	92
	93	//G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
	94	//{
	95	// if (instance == 0) instance = new G4SPSRandomGenerator();
	96	// return instance;
	97	//}
	98
	99	// Biasing methods
	100
[1347]	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;
[816]	107	}
	108
[1347]	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;
[816]	115	}
	116
[1347]	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;
[816]	123	}
	124
[1347]	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;
[816]	131	}
	132
[1347]	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;
[816]	139	}
	140
[1347]	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;
[816]	147	}
	148
[1347]	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;
[816]	155	}
	156
[1347]	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;
[816]	163	}
	164
[1347]	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	}
[816]	201	}
	202
[1347]	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	}
[816]	225
[1347]	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));
[816]	262	}
	263	}
	264
[1347]	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	}
[816]	287
[1347]	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));
[816]	316	}
	317	}
	318
[1347]	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	}
[816]	341
[1347]	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));
[816]	371	}
	372	}
	373
[1347]	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	}
[816]	398
[1347]	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));
[816]	429	}
	430	}
	431
[1347]	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	}
[816]	454
[1347]	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));
[816]	484	}
	485	}
	486
[1347]	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	}
[816]	541	}
	542
[1347]	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	}
[816]	567
[1347]	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));
[816]	599	}
	600	}
	601
[1347]	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	}
[816]	624
[1347]	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));
[816]	655	}
	656	}
	657

Note: See TracBrowser for help on using the repository browser.

Context Navigation

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

Download in other formats: