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