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

Last change on this file since 1197 was 816, checked in by garnier, 17 years ago

import all except CVS
File size: 21.4 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26///////////////////////////////////////////////////////////////////////////////
27//
28// MODULE: G4SPSRandomGenerator.cc
29//
30// Version: 1.0
31// Date: 5/02/04
32// Author: Fan Lei
33// Organisation: QinetiQ ltd.
34// Customer: ESA/ESTEC
35//
36///////////////////////////////////////////////////////////////////////////////
37//
38// CHANGE HISTORY
39// --------------
40//
41//
42// Version 1.0, 05/02/2004, Fan Lei, Created.
43// Based on the G4GeneralParticleSource class in Geant4 v6.0
44//
45///////////////////////////////////////////////////////////////////////////////
46//
47#include "G4PrimaryParticle.hh"
48#include "G4Event.hh"
49#include "Randomize.hh"
50#include <cmath>
51#include "G4TransportationManager.hh"
52#include "G4VPhysicalVolume.hh"
53#include "G4PhysicalVolumeStore.hh"
54#include "G4ParticleTable.hh"
55#include "G4ParticleDefinition.hh"
56#include "G4IonTable.hh"
57#include "G4Ions.hh"
58#include "G4TrackingManager.hh"
59#include "G4Track.hh"
60
61#include "G4SPSRandomGenerator.hh"
62
63//G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
64
65G4SPSRandomGenerator::G4SPSRandomGenerator()
66{
67 // Initialise all variables
68
69 // Bias variables
70 XBias = false;
71 IPDFXBias = false;
72 YBias = false;
73 IPDFYBias = false;
74 ZBias = false;
75 IPDFZBias = false;
76 ThetaBias = false;
77 IPDFThetaBias = false;
78 PhiBias = false;
79 IPDFPhiBias = false;
80 EnergyBias = false;
81 IPDFEnergyBias = false;
82 PosThetaBias = false;
83 IPDFPosThetaBias = false;
84 PosPhiBias = false;
85 IPDFPosPhiBias = false;
86 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4] = bweights[5] = bweights[6] = bweights[7]= bweights[8]= 1. ;
87 verbosityLevel = 0 ;
88}
89
90G4SPSRandomGenerator::~G4SPSRandomGenerator()
91{}
92
93//G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
94//{
95// if (instance == 0) instance = new G4SPSRandomGenerator();
96// return instance;
97//}
98
99// Biasing methods
100
101void G4SPSRandomGenerator::SetXBias(G4ThreeVector input)
102{
103 G4double ehi, val;
104 ehi = input.x();
105 val = input.y();
106 XBiasH.InsertValues(ehi, val);
107 XBias = true;
108}
109
110void G4SPSRandomGenerator::SetYBias(G4ThreeVector input)
111{
112 G4double ehi, val;
113 ehi = input.x();
114 val = input.y();
115 YBiasH.InsertValues(ehi, val);
116 YBias = true;
117}
118
119void G4SPSRandomGenerator::SetZBias(G4ThreeVector input)
120{
121 G4double ehi, val;
122 ehi = input.x();
123 val = input.y();
124 ZBiasH.InsertValues(ehi, val);
125 ZBias = true;
126}
127
128void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input)
129{
130 G4double ehi, val;
131 ehi = input.x();
132 val = input.y();
133 ThetaBiasH.InsertValues(ehi, val);
134 ThetaBias = true;
135}
136
137void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input)
138{
139 G4double ehi, val;
140 ehi = input.x();
141 val = input.y();
142 PhiBiasH.InsertValues(ehi, val);
143 PhiBias = true;
144}
145
146void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input)
147{
148 G4double ehi, val;
149 ehi = input.x();
150 val = input.y();
151 EnergyBiasH.InsertValues(ehi, val);
152 EnergyBias = true;
153}
154
155void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input)
156{
157 G4double ehi, val;
158 ehi = input.x();
159 val = input.y();
160 PosThetaBiasH.InsertValues(ehi, val);
161 PosThetaBias = true;
162}
163
164void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input)
165{
166 G4double ehi, val;
167 ehi = input.x();
168 val = input.y();
169 PosPhiBiasH.InsertValues(ehi, val);
170 PosPhiBias = true;
171}
172
173void G4SPSRandomGenerator::ReSetHist(G4String atype)
174{
175 if ( atype == "biasx") {
176 XBias = false ;
177 IPDFXBias = false;
178 XBiasH = IPDFXBiasH = ZeroPhysVector ;}
179 else if ( atype == "biasy") {
180 YBias = false ;
181 IPDFYBias = false;
182 YBiasH = IPDFYBiasH = ZeroPhysVector ;}
183 else if ( atype == "biasz") {
184 ZBias = false ;
185 IPDFZBias = false;
186 ZBiasH = IPDFZBiasH = ZeroPhysVector ;}
187 else if ( atype == "biast") {
188 ThetaBias = false ;
189 IPDFThetaBias = false;
190 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector ;}
191 else if ( atype == "biasp") {
192 PhiBias = false ;
193 IPDFPhiBias = false;
194 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector ;}
195 else if ( atype == "biase") {
196 EnergyBias = false ;
197 IPDFEnergyBias = false;
198 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector ;}
199 else if ( atype == "biaspt") {
200 PosThetaBias = false ;
201 IPDFPosThetaBias = false;
202 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector ;}
203 else if ( atype == "biaspp") {
204 PosPhiBias = false ;
205 IPDFPosPhiBias = false;
206 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector ;}
207 else {
208 G4cout << "Error, histtype not accepted " << G4endl;
209 }
210}
211
212G4double G4SPSRandomGenerator::GenRandX()
213{
214 if(verbosityLevel >= 1)
215 G4cout << "In GenRandX" << G4endl;
216 if(XBias == false)
217 {
218 // X is not biased
219 G4double rndm = G4UniformRand();
220 return(rndm);
221 }
222 else
223 {
224 // X is biased
225 if(IPDFXBias == false)
226 {
227 // IPDF has not been created, so create it
228 G4double bins[1024],vals[1024], sum;
229 G4int ii;
230 G4int maxbin = G4int(XBiasH.GetVectorLength());
231 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
232 vals[0] = XBiasH(size_t(0));
233 sum = vals[0];
234 for(ii=1;ii<maxbin;ii++)
235 {
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;
266 }
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));
278 }
279}
280
281G4double G4SPSRandomGenerator::GenRandY()
282{
283 if(verbosityLevel >= 1)
284 G4cout << "In GenRandY" << G4endl;
285 if(YBias == false)
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) {
325 if (rndm > IPDFYBiasH(biasn2))
326 biasn1 = biasn2;
327 else
328 biasn3 = biasn2;
329 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
330 }
331 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
332 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
333 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
334 G4double NatProb = xaxisu - xaxisl;
335 bweights[1] = NatProb/bweights[1];
336 if(verbosityLevel >= 1)
337 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
338 return(IPDFYBiasH.GetEnergy(rndm));
339 }
340}
341
342G4double G4SPSRandomGenerator::GenRandZ()
343{
344 if(verbosityLevel >= 1)
345 G4cout << "In GenRandZ" << G4endl;
346 if(ZBias == false)
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) {
387 if (rndm > IPDFZBiasH(biasn2))
388 biasn1 = biasn2;
389 else
390 biasn3 = biasn2;
391 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
392 }
393 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
394 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
395 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
396 G4double NatProb = xaxisu - xaxisl;
397 bweights[2] = NatProb/bweights[2];
398 if(verbosityLevel >= 1)
399 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
400 return(IPDFZBiasH.GetEnergy(rndm));
401 }
402}
403
404G4double G4SPSRandomGenerator::GenRandTheta()
405{
406 if(verbosityLevel >= 1)
407 {
408 G4cout << "In GenRandTheta" << G4endl;
409 G4cout << "Verbosity " << verbosityLevel << G4endl;
410 }
411 if(ThetaBias == false)
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) {
452 if (rndm > IPDFThetaBiasH(biasn2))
453 biasn1 = biasn2;
454 else
455 biasn3 = biasn2;
456 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
457 }
458 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
459 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
460 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
461 G4double NatProb = xaxisu - xaxisl;
462 bweights[3] = NatProb/bweights[3];
463 if(verbosityLevel >= 1)
464 G4cout << "Theta bin weight " << bweights[3] << " " << rndm << G4endl;
465 return(IPDFThetaBiasH.GetEnergy(rndm));
466 }
467}
468
469G4double G4SPSRandomGenerator::GenRandPhi()
470{
471 if(verbosityLevel >= 1)
472 G4cout << "In GenRandPhi" << G4endl;
473 if(PhiBias == false)
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) {
514 if (rndm > IPDFPhiBiasH(biasn2))
515 biasn1 = biasn2;
516 else
517 biasn3 = biasn2;
518 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
519 }
520 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
521 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
522 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
523 G4double NatProb = xaxisu - xaxisl;
524 bweights[4] = NatProb/bweights[4];
525 if(verbosityLevel >= 1)
526 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
527 return(IPDFPhiBiasH.GetEnergy(rndm));
528 }
529}
530
531G4double G4SPSRandomGenerator::GenRandEnergy()
532{
533 if(verbosityLevel >= 1)
534 G4cout << "In GenRandEnergy" << G4endl;
535 if(EnergyBias == false)
536 {
537 // Energy is not biased
538 G4double rndm = G4UniformRand();
539 return(rndm);
540 }
541 else {
542 // ENERGY is biased
543 if(IPDFEnergyBias == false) {
544 // IPDF has not been created, so create it
545 G4double bins[1024],vals[1024], sum;
546 G4int ii;
547 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
548 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
549 vals[0] = EnergyBiasH(size_t(0));
550 sum = vals[0];
551 for(ii=1;ii<maxbin;ii++) {
552 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
553 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii-1];
554 sum = sum + EnergyBiasH(size_t(ii));
555 }
556
557 for(ii=0;ii<maxbin;ii++) {
558 vals[ii] = vals[ii]/sum;
559 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
560 }
561 // Make IPDFEnergyBias = true
562 IPDFEnergyBias = true;
563 }
564 // IPDF has been create so carry on
565 G4double rndm = G4UniformRand();
566 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
567 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
568 G4int biasn1 = 0;
569 G4int biasn2 = numberOfBin/2;
570 G4int biasn3 = numberOfBin - 1;
571 while (biasn1 != biasn3 - 1) {
572 if (rndm > IPDFEnergyBiasH(biasn2))
573 biasn1 = biasn2;
574 else
575 biasn3 = biasn2;
576 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
577 }
578 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
579 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
580 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
581 G4double NatProb = xaxisu - xaxisl;
582 bweights[5] = NatProb/bweights[5];
583 if(verbosityLevel >= 1)
584 G4cout << "Energy bin weight " << bweights[5] << " " << rndm << G4endl;
585 return(IPDFEnergyBiasH.GetEnergy(rndm));
586 }
587}
588
589
590G4double G4SPSRandomGenerator::GenRandPosTheta()
591{
592 if(verbosityLevel >= 1)
593 {
594 G4cout << "In GenRandPosTheta" << G4endl;
595 G4cout << "Verbosity " << verbosityLevel << G4endl;
596 }
597 if(PosThetaBias == false)
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) {
638 if (rndm > IPDFPosThetaBiasH(biasn2))
639 biasn1 = biasn2;
640 else
641 biasn3 = biasn2;
642 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
643 }
644 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
645 G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
646 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
647 G4double NatProb = xaxisu - xaxisl;
648 bweights[6] = NatProb/bweights[6];
649 if(verbosityLevel >= 1)
650 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm << G4endl;
651 return(IPDFPosThetaBiasH.GetEnergy(rndm));
652 }
653}
654
655G4double G4SPSRandomGenerator::GenRandPosPhi()
656{
657 if(verbosityLevel >= 1)
658 G4cout << "In GenRandPosPhi" << G4endl;
659 if(PosPhiBias == false)
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) {
700 if (rndm > IPDFPosPhiBiasH(biasn2))
701 biasn1 = biasn2;
702 else
703 biasn3 = biasn2;
704 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2;
705 }
706 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
707 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1));
708 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
709 G4double NatProb = xaxisu - xaxisl;
710 bweights[7] = NatProb/bweights[7];
711 if(verbosityLevel >= 1)
712 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm << G4endl;
713 return(IPDFPosPhiBiasH.GetEnergy(rndm));
714 }
715}
716
717
718
719
720
Note: See TracBrowser for help on using the repository browser.