| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4SPSRandomGenerator class implementation << 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 ///////////////////////////////////////////////////////////////////////////////
27 // 46 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Revision: Andrea Dotti, SLAC <<
31 // ------------------------------------------- <<
32 <<
33 #include <cmath> <<
34 <<
35 #include "G4PrimaryParticle.hh" 47 #include "G4PrimaryParticle.hh"
36 #include "G4Event.hh" 48 #include "G4Event.hh"
37 #include "Randomize.hh" 49 #include "Randomize.hh"
>> 50 #include <cmath>
38 #include "G4TransportationManager.hh" 51 #include "G4TransportationManager.hh"
39 #include "G4VPhysicalVolume.hh" 52 #include "G4VPhysicalVolume.hh"
40 #include "G4PhysicalVolumeStore.hh" 53 #include "G4PhysicalVolumeStore.hh"
41 #include "G4ParticleTable.hh" 54 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh" 56 #include "G4IonTable.hh"
44 #include "G4Ions.hh" 57 #include "G4Ions.hh"
45 #include "G4TrackingManager.hh" 58 #include "G4TrackingManager.hh"
46 #include "G4Track.hh" 59 #include "G4Track.hh"
47 #include "G4AutoLock.hh" <<
48 60
49 #include "G4SPSRandomGenerator.hh" 61 #include "G4SPSRandomGenerator.hh"
50 62
51 G4SPSRandomGenerator::bweights_t::bweights_t() << 63 //G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
52 { << 64
53 for (double & i : w) { i = 1; } << 65 G4SPSRandomGenerator::G4SPSRandomGenerator() {
54 } << 66 // Initialise all variables
55 << 67
56 G4double& G4SPSRandomGenerator::bweights_t::op << 68 // Bias variables
57 { << 69 XBias = false;
58 return w[i]; << 70 IPDFXBias = false;
59 } << 71 YBias = false;
60 << 72 IPDFYBias = false;
61 G4SPSRandomGenerator::G4SPSRandomGenerator() << 73 ZBias = false;
62 { << 74 IPDFZBias = false;
63 // Initialise all variables << 75 ThetaBias = false;
64 << 76 IPDFThetaBias = false;
65 // Bias variables << 77 PhiBias = false;
66 // << 78 IPDFPhiBias = false;
67 XBias = false; << 79 EnergyBias = false;
68 IPDFXBias = false; << 80 IPDFEnergyBias = false;
69 YBias = false; << 81 PosThetaBias = false;
70 IPDFYBias = false; << 82 IPDFPosThetaBias = false;
71 ZBias = false; << 83 PosPhiBias = false;
72 IPDFZBias = false; << 84 IPDFPosPhiBias = false;
73 ThetaBias = false; << 85 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4]
74 IPDFThetaBias = false; << 86 = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
75 PhiBias = false; << 87 verbosityLevel = 0;
76 IPDFPhiBias = false; << 88 }
77 EnergyBias = false; << 89
78 IPDFEnergyBias = false; << 90 G4SPSRandomGenerator::~G4SPSRandomGenerator() {
79 PosThetaBias = false; << 91 }
80 IPDFPosThetaBias = false; << 92
81 PosPhiBias = false; << 93 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
82 IPDFPosPhiBias = false; << 94 //{
83 << 95 // if (instance == 0) instance = new G4SPSRandomGenerator();
84 verbosityLevel = 0; << 96 // return instance;
85 G4MUTEXINIT(mutex); << 97 //}
86 } <<
87 <<
88 G4SPSRandomGenerator::~G4SPSRandomGenerator() <<
89 { <<
90 G4MUTEXDESTROY(mutex); <<
91 } <<
92 98
93 // Biasing methods 99 // Biasing methods
94 100
95 void G4SPSRandomGenerator::SetXBias(const G4Th << 101 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
96 { << 102 G4double ehi, val;
97 G4AutoLock l(&mutex); << 103 ehi = input.x();
98 G4double ehi, val; << 104 val = input.y();
99 ehi = input.x(); << 105 XBiasH.InsertValues(ehi, val);
100 val = input.y(); << 106 XBias = true;
101 XBiasH.InsertValues(ehi, val); << 107 }
102 XBias = true; << 108
103 } << 109 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
104 << 110 G4double ehi, val;
105 void G4SPSRandomGenerator::SetYBias(const G4Th << 111 ehi = input.x();
106 { << 112 val = input.y();
107 G4AutoLock l(&mutex); << 113 YBiasH.InsertValues(ehi, val);
108 G4double ehi, val; << 114 YBias = true;
109 ehi = input.x(); << 115 }
110 val = input.y(); << 116
111 YBiasH.InsertValues(ehi, val); << 117 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
112 YBias = true; << 118 G4double ehi, val;
113 } << 119 ehi = input.x();
114 << 120 val = input.y();
115 void G4SPSRandomGenerator::SetZBias(const G4Th << 121 ZBiasH.InsertValues(ehi, val);
116 { << 122 ZBias = true;
117 G4AutoLock l(&mutex); << 123 }
118 G4double ehi, val; << 124
119 ehi = input.x(); << 125 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
120 val = input.y(); << 126 G4double ehi, val;
121 ZBiasH.InsertValues(ehi, val); << 127 ehi = input.x();
122 ZBias = true; << 128 val = input.y();
123 } << 129 ThetaBiasH.InsertValues(ehi, val);
124 << 130 ThetaBias = true;
125 void G4SPSRandomGenerator::SetThetaBias(const << 131 }
126 { << 132
127 G4AutoLock l(&mutex); << 133 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
128 G4double ehi, val; << 134 G4double ehi, val;
129 ehi = input.x(); << 135 ehi = input.x();
130 val = input.y(); << 136 val = input.y();
131 ThetaBiasH.InsertValues(ehi, val); << 137 PhiBiasH.InsertValues(ehi, val);
132 ThetaBias = true; << 138 PhiBias = true;
133 } << 139 }
134 << 140
135 void G4SPSRandomGenerator::SetPhiBias(const G4 << 141 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
136 { << 142 G4double ehi, val;
137 G4AutoLock l(&mutex); << 143 ehi = input.x();
138 G4double ehi, val; << 144 val = input.y();
139 ehi = input.x(); << 145 EnergyBiasH.InsertValues(ehi, val);
140 val = input.y(); << 146 EnergyBias = true;
141 PhiBiasH.InsertValues(ehi, val); << 147 }
142 PhiBias = true; << 148
143 } << 149 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
144 << 150 G4double ehi, val;
145 void G4SPSRandomGenerator::SetEnergyBias(const << 151 ehi = input.x();
146 { << 152 val = input.y();
147 G4AutoLock l(&mutex); << 153 PosThetaBiasH.InsertValues(ehi, val);
148 G4double ehi, val; << 154 PosThetaBias = true;
149 ehi = input.x(); << 155 }
150 val = input.y(); << 156
151 EnergyBiasH.InsertValues(ehi, val); << 157 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
152 EnergyBias = true; << 158 G4double ehi, val;
153 } << 159 ehi = input.x();
154 << 160 val = input.y();
155 void G4SPSRandomGenerator::SetPosThetaBias(con << 161 PosPhiBiasH.InsertValues(ehi, val);
156 { << 162 PosPhiBias = true;
157 G4AutoLock l(&mutex); << 163 }
158 G4double ehi, val; << 164
159 ehi = input.x(); << 165 void G4SPSRandomGenerator::ReSetHist(G4String atype) {
160 val = input.y(); << 166 if (atype == "biasx") {
161 PosThetaBiasH.InsertValues(ehi, val); << 167 XBias = false;
162 PosThetaBias = true; << 168 IPDFXBias = false;
163 } << 169 XBiasH = IPDFXBiasH = ZeroPhysVector;
164 << 170 } else if (atype == "biasy") {
165 void G4SPSRandomGenerator::SetPosPhiBias(const << 171 YBias = false;
166 { << 172 IPDFYBias = false;
167 G4AutoLock l(&mutex); << 173 YBiasH = IPDFYBiasH = ZeroPhysVector;
168 G4double ehi, val; << 174 } else if (atype == "biasz") {
169 ehi = input.x(); << 175 ZBias = false;
170 val = input.y(); << 176 IPDFZBias = false;
171 PosPhiBiasH.InsertValues(ehi, val); << 177 ZBiasH = IPDFZBiasH = ZeroPhysVector;
172 PosPhiBias = true; << 178 } else if (atype == "biast") {
173 } << 179 ThetaBias = false;
174 << 180 IPDFThetaBias = false;
175 void G4SPSRandomGenerator::SetIntensityWeight( << 181 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
176 { << 182 } else if (atype == "biasp") {
177 bweights.Get()[8] = weight; << 183 PhiBias = false;
178 } << 184 IPDFPhiBias = false;
179 << 185 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
180 G4double G4SPSRandomGenerator::GetBiasWeight() << 186 } else if (atype == "biase") {
181 { << 187 EnergyBias = false;
182 bweights_t& w = bweights.Get(); << 188 IPDFEnergyBias = false;
183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[ << 189 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
184 } << 190 } else if (atype == "biaspt") {
185 << 191 PosThetaBias = false;
186 void G4SPSRandomGenerator::SetVerbosity(G4int << 192 IPDFPosThetaBias = false;
187 { << 193 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
188 G4AutoLock l(&mutex); << 194 } else if (atype == "biaspp") {
189 verbosityLevel = a; << 195 PosPhiBias = false;
190 } << 196 IPDFPosPhiBias = false;
191 << 197 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
192 namespace << 198 } else {
193 { << 199 G4cout << "Error, histtype not accepted " << G4endl;
194 G4PhysicsFreeVector ZeroPhysVector; // for r << 200 }
195 } << 201 }
196 << 202
197 void G4SPSRandomGenerator::ReSetHist(const G4S << 203 G4double G4SPSRandomGenerator::GenRandX() {
198 { << 204 if (verbosityLevel >= 1)
199 G4AutoLock l(&mutex); << 205 G4cout << "In GenRandX" << G4endl;
200 if (atype == "biasx") { << 206 if (XBias == false) {
201 XBias = false; << 207 // X is not biased
202 IPDFXBias = false; << 208 G4double rndm = G4UniformRand();
203 local_IPDFXBias.Get().val = fa << 209 return (rndm);
204 XBiasH = IPDFXBiasH = ZeroPhys << 210 } else {
205 } else if (atype == "biasy") { << 211 // X is biased
206 YBias = false; << 212 if (IPDFXBias == false) {
207 IPDFYBias = false; << 213 // IPDF has not been created, so create it
208 local_IPDFYBias.Get().val = fa << 214 G4double bins[1024], vals[1024], sum;
209 YBiasH = IPDFYBiasH = ZeroPhys << 215 G4int ii;
210 } else if (atype == "biasz") { << 216 G4int maxbin = G4int(XBiasH.GetVectorLength());
211 ZBias = false; << 217 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
212 IPDFZBias = false; << 218 vals[0] = XBiasH(size_t(0));
213 local_IPDFZBias.Get().val = fa << 219 sum = vals[0];
214 ZBiasH = IPDFZBiasH = ZeroPhys << 220 for (ii = 1; ii < maxbin; ii++) {
215 } else if (atype == "biast") { << 221 bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
216 ThetaBias = false; << 222 vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
217 IPDFThetaBias = false; << 223 sum = sum + XBiasH(size_t(ii));
218 local_IPDFThetaBias.Get().val << 224 }
219 ThetaBiasH = IPDFThetaBiasH = << 225
220 } else if (atype == "biasp") { << 226 for (ii = 0; ii < maxbin; ii++) {
221 PhiBias = false; << 227 vals[ii] = vals[ii] / sum;
222 IPDFPhiBias = false; << 228 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
223 local_IPDFPhiBias.Get().val = << 229 }
224 PhiBiasH = IPDFPhiBiasH = Zero << 230 // Make IPDFXBias = true
225 } else if (atype == "biase") { << 231 IPDFXBias = true;
226 EnergyBias = false; << 232 }
227 IPDFEnergyBias = false; << 233 // IPDF has been create so carry on
228 local_IPDFEnergyBias.Get().val << 234 G4double rndm = G4UniformRand();
229 EnergyBiasH = IPDFEnergyBiasH << 235
230 } else if (atype == "biaspt") { << 236 // Calculate the weighting: Find the bin that the determined
231 PosThetaBias = false; << 237 // rndm is in and the weigthing will be the difference in the
232 IPDFPosThetaBias = false; << 238 // natural probability (from the x-axis) divided by the
233 local_IPDFPosThetaBias.Get().v << 239 // difference in the biased probability (the area).
234 PosThetaBiasH = IPDFPosThetaBi << 240 size_t numberOfBin = IPDFXBiasH.GetVectorLength();
235 } else if (atype == "biaspp") { << 241 G4int biasn1 = 0;
236 PosPhiBias = false; << 242 G4int biasn2 = numberOfBin / 2;
237 IPDFPosPhiBias = false; << 243 G4int biasn3 = numberOfBin - 1;
238 local_IPDFPosPhiBias.Get().val << 244 while (biasn1 != biasn3 - 1) {
239 PosPhiBiasH = IPDFPosPhiBiasH << 245 if (rndm > IPDFXBiasH(biasn2))
240 } else { << 246 biasn1 = biasn2;
241 G4cout << "Error, histtype not << 247 else
242 } << 248 biasn3 = biasn2;
243 } << 249 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
244 << 250 }
245 G4double G4SPSRandomGenerator::GenRandX() << 251 // retrieve the areas and then the x-axis values
246 { << 252 bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
247 if (verbosityLevel >= 1) << 253 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
248 G4cout << "In GenRandX" << G4endl; << 254 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
249 if (!XBias) << 255 G4double NatProb = xaxisu - xaxisl;
250 { << 256 //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
251 // X is not biased << 257 //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
252 G4double rndm = G4UniformRand(); << 258 bweights[0] = NatProb / bweights[0];
253 return (rndm); << 259 if (verbosityLevel >= 1)
254 } << 260 G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
255 << 261 return (IPDFXBiasH.GetEnergy(rndm));
256 // X is biased << 262 }
257 // This is shared among threads, and we need << 263 }
258 // only once. Multiple instances of this cla << 264
259 // so we rely on a class-private, thread-pri << 265 G4double G4SPSRandomGenerator::GenRandY() {
260 // to check if we need an initialiation. We << 266 if (verbosityLevel >= 1)
261 // because the Boolean on which we check is << 267 G4cout << "In GenRandY" << G4endl;
262 // << 268 if (YBias == false) {
263 if ( !local_IPDFXBias.Get().val ) << 269 // Y is not biased
264 { << 270 G4double rndm = G4UniformRand();
265 // For time that this thread arrived, here << 271 return (rndm);
266 // Now two cases are possible: it is the f << 272 } else {
267 // ANY thread has ever initialized this. << 273 // Y is biased
268 // Now we need a lock. In any case, the th << 274 if (IPDFYBias == false) {
269 // variable can now be set to true << 275 // IPDF has not been created, so create it
270 // << 276 G4double bins[1024], vals[1024], sum;
271 local_IPDFXBias.Get().val = true; << 277 G4int ii;
272 G4AutoLock l(&mutex); << 278 G4int maxbin = G4int(YBiasH.GetVectorLength());
273 if (!IPDFXBias) << 279 bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
274 { << 280 vals[0] = YBiasH(size_t(0));
275 // IPDF has not been created, so create << 281 sum = vals[0];
276 // << 282 for (ii = 1; ii < maxbin; ii++) {
277 G4double bins[1024], vals[1024], sum; << 283 bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
278 std::size_t ii; << 284 vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
279 std::size_t maxbin = XBiasH.GetVectorLen << 285 sum = sum + YBiasH(size_t(ii));
280 bins[0] = XBiasH.GetLowEdgeEnergy(0); << 286 }
281 vals[0] = XBiasH(0); << 287
282 sum = vals[0]; << 288 for (ii = 0; ii < maxbin; ii++) {
283 for (ii=1; ii<maxbin; ++ii) << 289 vals[ii] = vals[ii] / sum;
284 { << 290 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
285 bins[ii] = XBiasH.GetLowEdgeEnergy(ii) << 291 }
286 vals[ii] = XBiasH(ii) + vals[ii - 1]; << 292 // Make IPDFYBias = true
287 sum = sum + XBiasH(ii); << 293 IPDFYBias = true;
288 } << 294 }
289 << 295 // IPDF has been create so carry on
290 for (ii=0; ii<maxbin; ++ii) << 296 G4double rndm = G4UniformRand();
291 { << 297 size_t numberOfBin = IPDFYBiasH.GetVectorLength();
292 vals[ii] = vals[ii] / sum; << 298 G4int biasn1 = 0;
293 IPDFXBiasH.InsertValues(bins[ii], vals << 299 G4int biasn2 = numberOfBin / 2;
294 } << 300 G4int biasn3 = numberOfBin - 1;
295 IPDFXBias = true; << 301 while (biasn1 != biasn3 - 1) {
296 } << 302 if (rndm > IPDFYBiasH(biasn2))
297 } << 303 biasn1 = biasn2;
298 << 304 else
299 // IPDF has been create so carry on << 305 biasn3 = biasn2;
300 << 306 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
301 G4double rndm = G4UniformRand(); << 307 }
302 << 308 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
303 // Calculate the weighting: Find the bin tha << 309 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
304 // rndm is in and the weigthing will be the << 310 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
305 // natural probability (from the x-axis) div << 311 G4double NatProb = xaxisu - xaxisl;
306 // difference in the biased probability (the << 312 bweights[1] = NatProb / bweights[1];
307 // << 313 if (verbosityLevel >= 1)
308 std::size_t numberOfBin = IPDFXBiasH.GetVect << 314 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
309 std::size_t biasn1 = 0; << 315 return (IPDFYBiasH.GetEnergy(rndm));
310 std::size_t biasn2 = numberOfBin / 2; << 316 }
311 std::size_t biasn3 = numberOfBin - 1; << 317 }
312 while (biasn1 != biasn3 - 1) << 318
313 { << 319 G4double G4SPSRandomGenerator::GenRandZ() {
314 if (rndm > IPDFXBiasH(biasn2)) << 320 if (verbosityLevel >= 1)
315 { biasn1 = biasn2; } << 321 G4cout << "In GenRandZ" << G4endl;
316 else << 322 if (ZBias == false) {
317 { biasn3 = biasn2; } << 323 // Z is not biased
318 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / << 324 G4double rndm = G4UniformRand();
319 } << 325 return (rndm);
320 << 326 } else {
321 // Retrieve the areas and then the x-axis va << 327 // Z is biased
322 // << 328 if (IPDFZBias == false) {
323 bweights_t& w = bweights.Get(); << 329 // IPDF has not been created, so create it
324 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn << 330 G4double bins[1024], vals[1024], sum;
325 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnerg << 331 G4int ii;
326 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnerg << 332 G4int maxbin = G4int(ZBiasH.GetVectorLength());
327 G4double NatProb = xaxisu - xaxisl; << 333 bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
328 w[0] = NatProb / w[0]; << 334 vals[0] = ZBiasH(size_t(0));
329 if (verbosityLevel >= 1) << 335 sum = vals[0];
330 { << 336 for (ii = 1; ii < maxbin; ii++) {
331 G4cout << "X bin weight " << w[0] << " " < << 337 bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
332 } << 338 vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
333 return (IPDFXBiasH.GetEnergy(rndm)); << 339 sum = sum + ZBiasH(size_t(ii));
334 << 340 }
335 } << 341
336 << 342 for (ii = 0; ii < maxbin; ii++) {
337 G4double G4SPSRandomGenerator::GenRandY() << 343 vals[ii] = vals[ii] / sum;
338 { << 344 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
339 if (verbosityLevel >= 1) << 345 }
340 { << 346 // Make IPDFZBias = true
341 G4cout << "In GenRandY" << G4endl; << 347 IPDFZBias = true;
342 } << 348 }
343 << 349 // IPDF has been create so carry on
344 if (!YBias) // Y is not biased << 350 G4double rndm = G4UniformRand();
345 { << 351 // size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
346 G4double rndm = G4UniformRand(); << 352 size_t numberOfBin = IPDFZBiasH.GetVectorLength();
347 return (rndm); << 353 G4int biasn1 = 0;
348 } << 354 G4int biasn2 = numberOfBin / 2;
349 // Y is biased << 355 G4int biasn3 = numberOfBin - 1;
350 if ( !local_IPDFYBias.Get().val ) << 356 while (biasn1 != biasn3 - 1) {
351 { << 357 if (rndm > IPDFZBiasH(biasn2))
352 local_IPDFYBias.Get().val = true; << 358 biasn1 = biasn2;
353 G4AutoLock l(&mutex); << 359 else
354 if (!IPDFYBias) << 360 biasn3 = biasn2;
355 { << 361 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
356 // IPDF has not been created, so creat << 362 }
357 // << 363 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
358 G4double bins[1024], vals[1024], sum; << 364 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
359 std::size_t ii; << 365 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
360 std::size_t maxbin = YBiasH.GetVectorL << 366 G4double NatProb = xaxisu - xaxisl;
361 bins[0] = YBiasH.GetLowEdgeEnergy(0); << 367 bweights[2] = NatProb / bweights[2];
362 vals[0] = YBiasH(0); << 368 if (verbosityLevel >= 1)
363 sum = vals[0]; << 369 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
364 for (ii=1; ii<maxbin; ++ii) << 370 return (IPDFZBiasH.GetEnergy(rndm));
365 { << 371 }
366 bins[ii] = YBiasH.GetLowEdgeEnergy(i << 372 }
367 vals[ii] = YBiasH(ii) + vals[ii - 1] << 373
368 sum = sum + YBiasH(ii); << 374 G4double G4SPSRandomGenerator::GenRandTheta() {
369 } << 375 if (verbosityLevel >= 1) {
370 << 376 G4cout << "In GenRandTheta" << G4endl;
371 for (ii=0; ii<maxbin; ++ii) << 377 G4cout << "Verbosity " << verbosityLevel << G4endl;
372 { << 378 }
373 vals[ii] = vals[ii] / sum; << 379 if (ThetaBias == false) {
374 IPDFYBiasH.InsertValues(bins[ii], va << 380 // Theta is not biased
375 } << 381 G4double rndm = G4UniformRand();
376 IPDFYBias = true; << 382 return (rndm);
377 } << 383 } else {
378 } << 384 // Theta is biased
379 << 385 if (IPDFThetaBias == false) {
380 // IPDF has been created so carry on << 386 // IPDF has not been created, so create it
381 << 387 G4double bins[1024], vals[1024], sum;
382 G4double rndm = G4UniformRand(); << 388 G4int ii;
383 std::size_t numberOfBin = IPDFYBiasH.GetVe << 389 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
384 std::size_t biasn1 = 0; << 390 bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
385 std::size_t biasn2 = numberOfBin / 2; << 391 vals[0] = ThetaBiasH(size_t(0));
386 std::size_t biasn3 = numberOfBin - 1; << 392 sum = vals[0];
387 while (biasn1 != biasn3 - 1) << 393 for (ii = 1; ii < maxbin; ii++) {
388 { << 394 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
389 if (rndm > IPDFYBiasH(biasn2)) << 395 vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
390 { biasn1 = biasn2; } << 396 sum = sum + ThetaBiasH(size_t(ii));
391 else << 397 }
392 { biasn3 = biasn2; } << 398
393 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 399 for (ii = 0; ii < maxbin; ii++) {
394 } << 400 vals[ii] = vals[ii] / sum;
395 bweights_t& w = bweights.Get(); << 401 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
396 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(bia << 402 }
397 G4double xaxisl = IPDFYBiasH.GetLowEdgeEne << 403 // Make IPDFThetaBias = true
398 G4double xaxisu = IPDFYBiasH.GetLowEdgeEne << 404 IPDFThetaBias = true;
399 G4double NatProb = xaxisu - xaxisl; << 405 }
400 w[1] = NatProb / w[1]; << 406 // IPDF has been create so carry on
401 if (verbosityLevel >= 1) << 407 G4double rndm = G4UniformRand();
402 { << 408 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
403 G4cout << "Y bin weight " << w[1] << " " << 409 size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
404 } << 410 G4int biasn1 = 0;
405 return (IPDFYBiasH.GetEnergy(rndm)); << 411 G4int biasn2 = numberOfBin / 2;
406 << 412 G4int biasn3 = numberOfBin - 1;
407 } << 413 while (biasn1 != biasn3 - 1) {
408 << 414 if (rndm > IPDFThetaBiasH(biasn2))
409 G4double G4SPSRandomGenerator::GenRandZ() << 415 biasn1 = biasn2;
410 { << 416 else
411 if (verbosityLevel >= 1) << 417 biasn3 = biasn2;
412 { << 418 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
413 G4cout << "In GenRandZ" << G4endl; << 419 }
414 } << 420 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
415 << 421 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
416 if (!ZBias) // Z is not biased << 422 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
417 { << 423 G4double NatProb = xaxisu - xaxisl;
418 G4double rndm = G4UniformRand(); << 424 bweights[3] = NatProb / bweights[3];
419 return (rndm); << 425 if (verbosityLevel >= 1)
420 } << 426 G4cout << "Theta bin weight " << bweights[3] << " " << rndm
421 // Z is biased << 427 << G4endl;
422 if ( !local_IPDFZBias.Get().val ) << 428 return (IPDFThetaBiasH.GetEnergy(rndm));
423 { << 429 }
424 local_IPDFZBias.Get().val = true; << 430 }
425 G4AutoLock l(&mutex); << 431
426 if (!IPDFZBias) << 432 G4double G4SPSRandomGenerator::GenRandPhi() {
427 { << 433 if (verbosityLevel >= 1)
428 // IPDF has not been created, so creat << 434 G4cout << "In GenRandPhi" << G4endl;
429 // << 435 if (PhiBias == false) {
430 G4double bins[1024], vals[1024], sum; << 436 // Phi is not biased
431 std::size_t ii; << 437 G4double rndm = G4UniformRand();
432 std::size_t maxbin = ZBiasH.GetVectorL << 438 return (rndm);
433 bins[0] = ZBiasH.GetLowEdgeEnergy(0); << 439 } else {
434 vals[0] = ZBiasH(0); << 440 // Phi is biased
435 sum = vals[0]; << 441 if (IPDFPhiBias == false) {
436 for (ii=1; ii<maxbin; ++ii) << 442 // IPDF has not been created, so create it
437 { << 443 G4double bins[1024], vals[1024], sum;
438 bins[ii] = ZBiasH.GetLowEdgeEnergy(i << 444 G4int ii;
439 vals[ii] = ZBiasH(ii) + vals[ii - 1] << 445 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
440 sum = sum + ZBiasH(ii); << 446 bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
441 } << 447 vals[0] = PhiBiasH(size_t(0));
442 << 448 sum = vals[0];
443 for (ii=0; ii<maxbin; ++ii) << 449 for (ii = 1; ii < maxbin; ii++) {
444 { << 450 bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
445 vals[ii] = vals[ii] / sum; << 451 vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
446 IPDFZBiasH.InsertValues(bins[ii], va << 452 sum = sum + PhiBiasH(size_t(ii));
447 } << 453 }
448 IPDFZBias = true; << 454
449 } << 455 for (ii = 0; ii < maxbin; ii++) {
450 } << 456 vals[ii] = vals[ii] / sum;
451 << 457 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
452 // IPDF has been create so carry on << 458 }
453 << 459 // Make IPDFPhiBias = true
454 G4double rndm = G4UniformRand(); << 460 IPDFPhiBias = true;
455 std::size_t numberOfBin = IPDFZBiasH.GetVe << 461 }
456 std::size_t biasn1 = 0; << 462 // IPDF has been create so carry on
457 std::size_t biasn2 = numberOfBin / 2; << 463 G4double rndm = G4UniformRand();
458 std::size_t biasn3 = numberOfBin - 1; << 464 // size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
459 while (biasn1 != biasn3 - 1) << 465 size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
460 { << 466 G4int biasn1 = 0;
461 if (rndm > IPDFZBiasH(biasn2)) << 467 G4int biasn2 = numberOfBin / 2;
462 { biasn1 = biasn2; } << 468 G4int biasn3 = numberOfBin - 1;
463 else << 469 while (biasn1 != biasn3 - 1) {
464 { biasn3 = biasn2; } << 470 if (rndm > IPDFPhiBiasH(biasn2))
465 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 471 biasn1 = biasn2;
466 } << 472 else
467 bweights_t& w = bweights.Get(); << 473 biasn3 = biasn2;
468 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(bia << 474 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
469 G4double xaxisl = IPDFZBiasH.GetLowEdgeEne << 475 }
470 G4double xaxisu = IPDFZBiasH.GetLowEdgeEne << 476 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
471 G4double NatProb = xaxisu - xaxisl; << 477 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
472 w[2] = NatProb / w[2]; << 478 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
473 if (verbosityLevel >= 1) << 479 G4double NatProb = xaxisu - xaxisl;
474 { << 480 bweights[4] = NatProb / bweights[4];
475 G4cout << "Z bin weight " << w[2] << " " << 481 if (verbosityLevel >= 1)
476 } << 482 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
477 return (IPDFZBiasH.GetEnergy(rndm)); << 483 return (IPDFPhiBiasH.GetEnergy(rndm));
478 << 484 }
479 } << 485 }
480 << 486
481 G4double G4SPSRandomGenerator::GenRandTheta() << 487 G4double G4SPSRandomGenerator::GenRandEnergy() {
482 { << 488 if (verbosityLevel >= 1)
483 if (verbosityLevel >= 1) << 489 G4cout << "In GenRandEnergy" << G4endl;
484 { << 490 if (EnergyBias == false) {
485 G4cout << "In GenRandTheta" << G4endl; << 491 // Energy is not biased
486 G4cout << "Verbosity " << verbosityLevel < << 492 G4double rndm = G4UniformRand();
487 } << 493 return (rndm);
488 << 494 } else {
489 if (!ThetaBias) // Theta is not biased << 495 // ENERGY is biased
490 { << 496 if (IPDFEnergyBias == false) {
491 G4double rndm = G4UniformRand(); << 497 // IPDF has not been created, so create it
492 return (rndm); << 498 G4double bins[1024], vals[1024], sum;
493 } << 499 G4int ii;
494 // Theta is biased << 500 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
495 if ( !local_IPDFThetaBias.Get().val ) << 501 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
496 { << 502 vals[0] = EnergyBiasH(size_t(0));
497 local_IPDFThetaBias.Get().val = true; << 503 sum = vals[0];
498 G4AutoLock l(&mutex); << 504 for (ii = 1; ii < maxbin; ii++) {
499 if (!IPDFThetaBias) << 505 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
500 { << 506 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
501 // IPDF has not been created, so creat << 507 sum = sum + EnergyBiasH(size_t(ii));
502 // << 508 }
503 G4double bins[1024], vals[1024], sum; << 509 IPDFEnergyBiasH = ZeroPhysVector;
504 std::size_t ii; << 510 for (ii = 0; ii < maxbin; ii++) {
505 std::size_t maxbin = ThetaBiasH.GetVec << 511 vals[ii] = vals[ii] / sum;
506 bins[0] = ThetaBiasH.GetLowEdgeEnergy( << 512 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
507 vals[0] = ThetaBiasH(0); << 513 }
508 sum = vals[0]; << 514 // Make IPDFEnergyBias = true
509 for (ii=1; ii<maxbin; ++ii) << 515 IPDFEnergyBias = true;
510 { << 516 }
511 bins[ii] = ThetaBiasH.GetLowEdgeEner << 517 // IPDF has been create so carry on
512 vals[ii] = ThetaBiasH(ii) + vals[ii << 518 G4double rndm = G4UniformRand();
513 sum = sum + ThetaBiasH(ii); << 519 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
514 } << 520 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
515 << 521 G4int biasn1 = 0;
516 for (ii=0; ii<maxbin; ++ii) << 522 G4int biasn2 = numberOfBin / 2;
517 { << 523 G4int biasn3 = numberOfBin - 1;
518 vals[ii] = vals[ii] / sum; << 524 while (biasn1 != biasn3 - 1) {
519 IPDFThetaBiasH.InsertValues(bins[ii] << 525 if (rndm > IPDFEnergyBiasH(biasn2))
520 } << 526 biasn1 = biasn2;
521 IPDFThetaBias = true; << 527 else
522 } << 528 biasn3 = biasn2;
523 } << 529 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
524 << 530 }
525 // IPDF has been create so carry on << 531 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
526 << 532 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
527 G4double rndm = G4UniformRand(); << 533 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
528 std::size_t numberOfBin = IPDFThetaBiasH.G << 534 G4double NatProb = xaxisu - xaxisl;
529 std::size_t biasn1 = 0; << 535 bweights[5] = NatProb / bweights[5];
530 std::size_t biasn2 = numberOfBin / 2; << 536 if (verbosityLevel >= 1)
531 std::size_t biasn3 = numberOfBin - 1; << 537 G4cout << "Energy bin weight " << bweights[5] << " " << rndm
532 while (biasn1 != biasn3 - 1) << 538 << G4endl;
533 { << 539 return (IPDFEnergyBiasH.GetEnergy(rndm));
534 if (rndm > IPDFThetaBiasH(biasn2)) << 540 }
535 { biasn1 = biasn2; } << 541 }
536 else << 542
537 { biasn3 = biasn2; } << 543 G4double G4SPSRandomGenerator::GenRandPosTheta() {
538 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 544 if (verbosityLevel >= 1) {
539 } << 545 G4cout << "In GenRandPosTheta" << G4endl;
540 bweights_t& w = bweights.Get(); << 546 G4cout << "Verbosity " << verbosityLevel << G4endl;
541 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaB << 547 }
542 G4double xaxisl = IPDFThetaBiasH.GetLowEdg << 548 if (PosThetaBias == false) {
543 G4double xaxisu = IPDFThetaBiasH.GetLowEdg << 549 // Theta is not biased
544 G4double NatProb = xaxisu - xaxisl; << 550 G4double rndm = G4UniformRand();
545 w[3] = NatProb / w[3]; << 551 return (rndm);
546 if (verbosityLevel >= 1) << 552 } else {
547 { << 553 // Theta is biased
548 G4cout << "Theta bin weight " << w[3] << << 554 if (IPDFPosThetaBias == false) {
549 } << 555 // IPDF has not been created, so create it
550 return (IPDFThetaBiasH.GetEnergy(rndm)); << 556 G4double bins[1024], vals[1024], sum;
551 << 557 G4int ii;
552 } << 558 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
553 << 559 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
554 G4double G4SPSRandomGenerator::GenRandPhi() << 560 vals[0] = PosThetaBiasH(size_t(0));
555 { << 561 sum = vals[0];
556 if (verbosityLevel >= 1) << 562 for (ii = 1; ii < maxbin; ii++) {
557 { << 563 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
558 G4cout << "In GenRandPhi" << G4endl; << 564 vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
559 } << 565 sum = sum + PosThetaBiasH(size_t(ii));
560 << 566 }
561 if (!PhiBias) // Phi is not biased << 567
562 { << 568 for (ii = 0; ii < maxbin; ii++) {
563 G4double rndm = G4UniformRand(); << 569 vals[ii] = vals[ii] / sum;
564 return (rndm); << 570 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
565 } << 571 }
566 // Phi is biased << 572 // Make IPDFThetaBias = true
567 if ( !local_IPDFPhiBias.Get().val ) << 573 IPDFPosThetaBias = true;
568 { << 574 }
569 local_IPDFPhiBias.Get().val = true; << 575 // IPDF has been create so carry on
570 G4AutoLock l(&mutex); << 576 G4double rndm = G4UniformRand();
571 if (!IPDFPhiBias) << 577 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
572 { << 578 size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
573 // IPDF has not been created, so creat << 579 G4int biasn1 = 0;
574 // << 580 G4int biasn2 = numberOfBin / 2;
575 G4double bins[1024], vals[1024], sum; << 581 G4int biasn3 = numberOfBin - 1;
576 std::size_t ii; << 582 while (biasn1 != biasn3 - 1) {
577 std::size_t maxbin = PhiBiasH.GetVecto << 583 if (rndm > IPDFPosThetaBiasH(biasn2))
578 bins[0] = PhiBiasH.GetLowEdgeEnergy(0) << 584 biasn1 = biasn2;
579 vals[0] = PhiBiasH(0); << 585 else
580 sum = vals[0]; << 586 biasn3 = biasn2;
581 for (ii=1; ii<maxbin; ++ii) << 587 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
582 { << 588 }
583 bins[ii] = PhiBiasH.GetLowEdgeEnergy << 589 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
584 vals[ii] = PhiBiasH(ii) + vals[ii - << 590 G4double xaxisl =
585 sum = sum + PhiBiasH(ii); << 591 IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
586 } << 592 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
587 << 593 G4double NatProb = xaxisu - xaxisl;
588 for (ii=0; ii<maxbin; ++ii) << 594 bweights[6] = NatProb / bweights[6];
589 { << 595 if (verbosityLevel >= 1)
590 vals[ii] = vals[ii] / sum; << 596 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
591 IPDFPhiBiasH.InsertValues(bins[ii], << 597 << G4endl;
592 } << 598 return (IPDFPosThetaBiasH.GetEnergy(rndm));
593 IPDFPhiBias = true; << 599 }
594 } << 600 }
595 } << 601
596 << 602 G4double G4SPSRandomGenerator::GenRandPosPhi() {
597 // IPDF has been create so carry on << 603 if (verbosityLevel >= 1)
598 << 604 G4cout << "In GenRandPosPhi" << G4endl;
599 G4double rndm = G4UniformRand(); << 605 if (PosPhiBias == false) {
600 std::size_t numberOfBin = IPDFPhiBiasH.Get << 606 // PosPhi is not biased
601 std::size_t biasn1 = 0; << 607 G4double rndm = G4UniformRand();
602 std::size_t biasn2 = numberOfBin / 2; << 608 return (rndm);
603 std::size_t biasn3 = numberOfBin - 1; << 609 } else {
604 while (biasn1 != biasn3 - 1) << 610 // PosPhi is biased
605 { << 611 if (IPDFPosPhiBias == false) {
606 if (rndm > IPDFPhiBiasH(biasn2)) << 612 // IPDF has not been created, so create it
607 { biasn1 = biasn2; } << 613 G4double bins[1024], vals[1024], sum;
608 else << 614 G4int ii;
609 { biasn3 = biasn2; } << 615 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
610 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 616 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
611 } << 617 vals[0] = PosPhiBiasH(size_t(0));
612 bweights_t& w = bweights.Get(); << 618 sum = vals[0];
613 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH << 619 for (ii = 1; ii < maxbin; ii++) {
614 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeE << 620 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
615 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeE << 621 vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
616 G4double NatProb = xaxisu - xaxisl; << 622 sum = sum + PosPhiBiasH(size_t(ii));
617 w[4] = NatProb / w[4]; << 623 }
618 if (verbosityLevel >= 1) << 624
619 { << 625 for (ii = 0; ii < maxbin; ii++) {
620 G4cout << "Phi bin weight " << w[4] << " << 626 vals[ii] = vals[ii] / sum;
621 } << 627 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
622 return (IPDFPhiBiasH.GetEnergy(rndm)); << 628 }
623 << 629 // Make IPDFPosPhiBias = true
624 } << 630 IPDFPosPhiBias = true;
625 << 631 }
626 G4double G4SPSRandomGenerator::GenRandEnergy() << 632 // IPDF has been create so carry on
627 { << 633 G4double rndm = G4UniformRand();
628 if (verbosityLevel >= 1) << 634 // size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
629 { << 635 size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
630 G4cout << "In GenRandEnergy" << G4endl; << 636 G4int biasn1 = 0;
631 } << 637 G4int biasn2 = numberOfBin / 2;
632 << 638 G4int biasn3 = numberOfBin - 1;
633 if (!EnergyBias) // Energy is not biased << 639 while (biasn1 != biasn3 - 1) {
634 { << 640 if (rndm > IPDFPosPhiBiasH(biasn2))
635 G4double rndm = G4UniformRand(); << 641 biasn1 = biasn2;
636 return (rndm); << 642 else
637 } << 643 biasn3 = biasn2;
638 // Energy is biased << 644 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
639 if ( !local_IPDFEnergyBias.Get().val ) << 645 }
640 { << 646 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
641 local_IPDFEnergyBias.Get().val = true; << 647 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
642 G4AutoLock l(&mutex); << 648 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
643 if (!IPDFEnergyBias) << 649 G4double NatProb = xaxisu - xaxisl;
644 { << 650 bweights[7] = NatProb / bweights[7];
645 // IPDF has not been created, so creat << 651 if (verbosityLevel >= 1)
646 // << 652 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
647 G4double bins[1024], vals[1024], sum; << 653 << G4endl;
648 std::size_t ii; << 654 return (IPDFPosPhiBiasH.GetEnergy(rndm));
649 std::size_t maxbin = EnergyBiasH.GetVe << 655 }
650 bins[0] = EnergyBiasH.GetLowEdgeEnergy <<
651 vals[0] = EnergyBiasH(0); <<
652 sum = vals[0]; <<
653 for (ii=1; ii<maxbin; ++ii) <<
654 { <<
655 bins[ii] = EnergyBiasH.GetLowEdgeEne <<
656 vals[ii] = EnergyBiasH(ii) + vals[ii <<
657 sum = sum + EnergyBiasH(ii); <<
658 } <<
659 IPDFEnergyBiasH = ZeroPhysVector; <<
660 for (ii=0; ii<maxbin; ++ii) <<
661 { <<
662 vals[ii] = vals[ii] / sum; <<
663 IPDFEnergyBiasH.InsertValues(bins[ii <<
664 } <<
665 IPDFEnergyBias = true; <<
666 } <<
667 } <<
668 <<
669 // IPDF has been create so carry on <<
670 <<
671 G4double rndm = G4UniformRand(); <<
672 std::size_t numberOfBin = IPDFEnergyBiasH. <<
673 std::size_t biasn1 = 0; <<
674 std::size_t biasn2 = numberOfBin / 2; <<
675 std::size_t biasn3 = numberOfBin - 1; <<
676 while (biasn1 != biasn3 - 1) <<
677 { <<
678 if (rndm > IPDFEnergyBiasH(biasn2)) <<
679 { biasn1 = biasn2; } <<
680 else <<
681 { biasn3 = biasn2; } <<
682 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
683 } <<
684 bweights_t& w = bweights.Get(); <<
685 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnerg <<
686 G4double xaxisl = IPDFEnergyBiasH.GetLowEd <<
687 G4double xaxisu = IPDFEnergyBiasH.GetLowEd <<
688 G4double NatProb = xaxisu - xaxisl; <<
689 w[5] = NatProb / w[5]; <<
690 if (verbosityLevel >= 1) <<
691 { <<
692 G4cout << "Energy bin weight " << w[5] < <<
693 } <<
694 return (IPDFEnergyBiasH.GetEnergy(rndm)); <<
695 <<
696 } <<
697 <<
698 G4double G4SPSRandomGenerator::GenRandPosTheta <<
699 { <<
700 if (verbosityLevel >= 1) <<
701 { <<
702 G4cout << "In GenRandPosTheta" << G4endl; <<
703 G4cout << "Verbosity " << verbosityLevel < <<
704 } <<
705 <<
706 if (!PosThetaBias) // Theta is not biased <<
707 { <<
708 G4double rndm = G4UniformRand(); <<
709 return (rndm); <<
710 } <<
711 // Theta is biased <<
712 if ( !local_IPDFPosThetaBias.Get().val ) <<
713 { <<
714 local_IPDFPosThetaBias.Get().val = true; <<
715 G4AutoLock l(&mutex); <<
716 if (!IPDFPosThetaBias) <<
717 { <<
718 // IPDF has not been created, so creat <<
719 // <<
720 G4double bins[1024], vals[1024], sum; <<
721 std::size_t ii; <<
722 std::size_t maxbin = PosThetaBiasH.Get <<
723 bins[0] = PosThetaBiasH.GetLowEdgeEner <<
724 vals[0] = PosThetaBiasH(0); <<
725 sum = vals[0]; <<
726 for (ii=1; ii<maxbin; ++ii) <<
727 { <<
728 bins[ii] = PosThetaBiasH.GetLowEdgeE <<
729 vals[ii] = PosThetaBiasH(ii) + vals[ <<
730 sum = sum + PosThetaBiasH(ii); <<
731 } <<
732 <<
733 for (ii=0; ii<maxbin; ++ii) <<
734 { <<
735 vals[ii] = vals[ii] / sum; <<
736 IPDFPosThetaBiasH.InsertValues(bins[ <<
737 } <<
738 IPDFPosThetaBias = true; <<
739 } <<
740 } <<
741 <<
742 // IPDF has been create so carry on <<
743 // <<
744 G4double rndm = G4UniformRand(); <<
745 std::size_t numberOfBin = IPDFPosThetaBias <<
746 std::size_t biasn1 = 0; <<
747 std::size_t biasn2 = numberOfBin / 2; <<
748 std::size_t biasn3 = numberOfBin - 1; <<
749 while (biasn1 != biasn3 - 1) <<
750 { <<
751 if (rndm > IPDFPosThetaBiasH(biasn2)) <<
752 { biasn1 = biasn2; } <<
753 else <<
754 { biasn3 = biasn2; } <<
755 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
756 } <<
757 bweights_t& w = bweights.Get(); <<
758 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPos <<
759 G4double xaxisl = IPDFPosThetaBiasH.GetLow <<
760 G4double xaxisu = IPDFPosThetaBiasH.GetLow <<
761 G4double NatProb = xaxisu - xaxisl; <<
762 w[6] = NatProb / w[6]; <<
763 if (verbosityLevel >= 1) <<
764 { <<
765 G4cout << "PosTheta bin weight " << w[6] <<
766 } <<
767 return (IPDFPosThetaBiasH.GetEnergy(rndm)) <<
768 <<
769 } <<
770 <<
771 G4double G4SPSRandomGenerator::GenRandPosPhi() <<
772 { <<
773 if (verbosityLevel >= 1) <<
774 { <<
775 G4cout << "In GenRandPosPhi" << G4endl; <<
776 } <<
777 <<
778 if (!PosPhiBias) // PosPhi is not biased <<
779 { <<
780 G4double rndm = G4UniformRand(); <<
781 return (rndm); <<
782 } <<
783 // PosPhi is biased <<
784 if (!local_IPDFPosPhiBias.Get().val ) <<
785 { <<
786 local_IPDFPosPhiBias.Get().val = true; <<
787 G4AutoLock l(&mutex); <<
788 if (!IPDFPosPhiBias) <<
789 { <<
790 // IPDF has not been created, so creat <<
791 // <<
792 G4double bins[1024], vals[1024], sum; <<
793 std::size_t ii; <<
794 std::size_t maxbin = PosPhiBiasH.GetVe <<
795 bins[0] = PosPhiBiasH.GetLowEdgeEnergy <<
796 vals[0] = PosPhiBiasH(0); <<
797 sum = vals[0]; <<
798 for (ii=1; ii<maxbin; ++ii) <<
799 { <<
800 bins[ii] = PosPhiBiasH.GetLowEdgeEne <<
801 vals[ii] = PosPhiBiasH(ii) + vals[ii <<
802 sum = sum + PosPhiBiasH(ii); <<
803 } <<
804 <<
805 for (ii=0; ii<maxbin; ++ii) <<
806 { <<
807 vals[ii] = vals[ii] / sum; <<
808 IPDFPosPhiBiasH.InsertValues(bins[ii <<
809 } <<
810 IPDFPosPhiBias = true; <<
811 } <<
812 } <<
813 <<
814 // IPDF has been create so carry on <<
815 <<
816 G4double rndm = G4UniformRand(); <<
817 std::size_t numberOfBin = IPDFPosPhiBiasH. <<
818 std::size_t biasn1 = 0; <<
819 std::size_t biasn2 = numberOfBin / 2; <<
820 std::size_t biasn3 = numberOfBin - 1; <<
821 while (biasn1 != biasn3 - 1) <<
822 { <<
823 if (rndm > IPDFPosPhiBiasH(biasn2)) <<
824 { biasn1 = biasn2; } <<
825 else <<
826 { biasn3 = biasn2; } <<
827 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
828 } <<
829 bweights_t& w = bweights.Get(); <<
830 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPh <<
831 G4double xaxisl = IPDFPosPhiBiasH.GetLowEd <<
832 G4double xaxisu = IPDFPosPhiBiasH.GetLowEd <<
833 G4double NatProb = xaxisu - xaxisl; <<
834 w[7] = NatProb / w[7]; <<
835 if (verbosityLevel >= 1) <<
836 { <<
837 G4cout << "PosPhi bin weight " << w[7] < <<
838 } <<
839 return (IPDFPosPhiBiasH.GetEnergy(rndm)); <<
840 <<
841 } 656 }
>> 657
842 658