Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SPSRandomGenerator.cc
Go to the documentation of this file.
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 //
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 //
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 //
46 //
47 #include "G4PrimaryParticle.hh"
48 #include "G4Event.hh"
49 #include "Randomize.hh"
50 #include <cmath>
51 #include "G4TransportationManager.hh"
52 #include "G4VPhysicalVolume.hh"
53 #include "G4PhysicalVolumeStore.hh"
54 #include "G4ParticleTable.hh"
55 #include "G4ParticleDefinition.hh"
56 #include "G4IonTable.hh"
57 #include "G4Ions.hh"
58 #include "G4TrackingManager.hh"
59 #include "G4Track.hh"
60 
61 #include "G4SPSRandomGenerator.hh"
62 
63 //G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
64 
65 G4SPSRandomGenerator::G4SPSRandomGenerator()
66 {
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]
87  = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
88  verbosityLevel = 0;
89 }
90 
91 G4SPSRandomGenerator::~G4SPSRandomGenerator() {
92 }
93 
94 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
95 //{
96 // if (instance == 0) instance = new G4SPSRandomGenerator();
97 // return instance;
98 //}
99 
100 // Biasing methods
101 
102 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
103  G4double ehi, val;
104  ehi = input.x();
105  val = input.y();
106  XBiasH.InsertValues(ehi, val);
107  XBias = true;
108 }
109 
110 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
111  G4double ehi, val;
112  ehi = input.x();
113  val = input.y();
114  YBiasH.InsertValues(ehi, val);
115  YBias = true;
116 }
117 
118 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
119  G4double ehi, val;
120  ehi = input.x();
121  val = input.y();
122  ZBiasH.InsertValues(ehi, val);
123  ZBias = true;
124 }
125 
126 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
127  G4double ehi, val;
128  ehi = input.x();
129  val = input.y();
130  ThetaBiasH.InsertValues(ehi, val);
131  ThetaBias = true;
132 }
133 
134 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
135  G4double ehi, val;
136  ehi = input.x();
137  val = input.y();
138  PhiBiasH.InsertValues(ehi, val);
139  PhiBias = true;
140 }
141 
142 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
143  G4double ehi, val;
144  ehi = input.x();
145  val = input.y();
146  EnergyBiasH.InsertValues(ehi, val);
147  EnergyBias = true;
148 }
149 
150 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
151  G4double ehi, val;
152  ehi = input.x();
153  val = input.y();
154  PosThetaBiasH.InsertValues(ehi, val);
155  PosThetaBias = true;
156 }
157 
158 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
159  G4double ehi, val;
160  ehi = input.x();
161  val = input.y();
162  PosPhiBiasH.InsertValues(ehi, val);
163  PosPhiBias = true;
164 }
165 
166 void G4SPSRandomGenerator::ReSetHist(G4String atype) {
167  if (atype == "biasx") {
168  XBias = false;
169  IPDFXBias = false;
170  XBiasH = IPDFXBiasH = ZeroPhysVector;
171  } else if (atype == "biasy") {
172  YBias = false;
173  IPDFYBias = false;
174  YBiasH = IPDFYBiasH = ZeroPhysVector;
175  } else if (atype == "biasz") {
176  ZBias = false;
177  IPDFZBias = false;
178  ZBiasH = IPDFZBiasH = ZeroPhysVector;
179  } else if (atype == "biast") {
180  ThetaBias = false;
181  IPDFThetaBias = false;
182  ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
183  } else if (atype == "biasp") {
184  PhiBias = false;
185  IPDFPhiBias = false;
186  PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
187  } else if (atype == "biase") {
188  EnergyBias = false;
189  IPDFEnergyBias = false;
190  EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
191  } else if (atype == "biaspt") {
192  PosThetaBias = false;
193  IPDFPosThetaBias = false;
194  PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
195  } else if (atype == "biaspp") {
196  PosPhiBias = false;
197  IPDFPosPhiBias = false;
198  PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
199  } else {
200  G4cout << "Error, histtype not accepted " << G4endl;
201  }
202 }
203 
204 G4double G4SPSRandomGenerator::GenRandX() {
205  if (verbosityLevel >= 1)
206  G4cout << "In GenRandX" << G4endl;
207  if (XBias == false) {
208  // X is not biased
209  G4double rndm = G4UniformRand();
210  return (rndm);
211  } else {
212  // X is biased
213  if (IPDFXBias == false) {
214  // IPDF has not been created, so create it
215  G4double bins[1024], vals[1024], sum;
216  G4int ii;
217  G4int maxbin = G4int(XBiasH.GetVectorLength());
218  bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
219  vals[0] = XBiasH(size_t(0));
220  sum = vals[0];
221  for (ii = 1; ii < maxbin; ii++) {
222  bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
223  vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
224  sum = sum + XBiasH(size_t(ii));
225  }
226 
227  for (ii = 0; ii < maxbin; ii++) {
228  vals[ii] = vals[ii] / sum;
229  IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
230  }
231  // Make IPDFXBias = true
232  IPDFXBias = true;
233  }
234  // IPDF has been create so carry on
235  G4double rndm = G4UniformRand();
236 
237  // Calculate the weighting: Find the bin that the determined
238  // rndm is in and the weigthing will be the difference in the
239  // natural probability (from the x-axis) divided by the
240  // difference in the biased probability (the area).
241  size_t numberOfBin = IPDFXBiasH.GetVectorLength();
242  G4int biasn1 = 0;
243  G4int biasn2 = numberOfBin / 2;
244  G4int biasn3 = numberOfBin - 1;
245  while (biasn1 != biasn3 - 1) {
246  if (rndm > IPDFXBiasH(biasn2))
247  biasn1 = biasn2;
248  else
249  biasn3 = biasn2;
250  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
251  }
252  // retrieve the areas and then the x-axis values
253  bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
254  G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
255  G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
256  G4double NatProb = xaxisu - xaxisl;
257  //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
258  //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
259  bweights[0] = NatProb / bweights[0];
260  if (verbosityLevel >= 1)
261  G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
262  return (IPDFXBiasH.GetEnergy(rndm));
263  }
264 }
265 
266 G4double G4SPSRandomGenerator::GenRandY() {
267  if (verbosityLevel >= 1)
268  G4cout << "In GenRandY" << G4endl;
269  if (YBias == false) {
270  // Y is not biased
271  G4double rndm = G4UniformRand();
272  return (rndm);
273  } else {
274  // Y is biased
275  if (IPDFYBias == false) {
276  // IPDF has not been created, so create it
277  G4double bins[1024], vals[1024], sum;
278  G4int ii;
279  G4int maxbin = G4int(YBiasH.GetVectorLength());
280  bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
281  vals[0] = YBiasH(size_t(0));
282  sum = vals[0];
283  for (ii = 1; ii < maxbin; ii++) {
284  bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
285  vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
286  sum = sum + YBiasH(size_t(ii));
287  }
288 
289  for (ii = 0; ii < maxbin; ii++) {
290  vals[ii] = vals[ii] / sum;
291  IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
292  }
293  // Make IPDFYBias = true
294  IPDFYBias = true;
295  }
296  // IPDF has been create so carry on
297  G4double rndm = G4UniformRand();
298  size_t numberOfBin = IPDFYBiasH.GetVectorLength();
299  G4int biasn1 = 0;
300  G4int biasn2 = numberOfBin / 2;
301  G4int biasn3 = numberOfBin - 1;
302  while (biasn1 != biasn3 - 1) {
303  if (rndm > IPDFYBiasH(biasn2))
304  biasn1 = biasn2;
305  else
306  biasn3 = biasn2;
307  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
308  }
309  bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
310  G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
311  G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
312  G4double NatProb = xaxisu - xaxisl;
313  bweights[1] = NatProb / bweights[1];
314  if (verbosityLevel >= 1)
315  G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
316  return (IPDFYBiasH.GetEnergy(rndm));
317  }
318 }
319 
320 G4double G4SPSRandomGenerator::GenRandZ() {
321  if (verbosityLevel >= 1)
322  G4cout << "In GenRandZ" << G4endl;
323  if (ZBias == false) {
324  // Z is not biased
325  G4double rndm = G4UniformRand();
326  return (rndm);
327  } else {
328  // Z is biased
329  if (IPDFZBias == false) {
330  // IPDF has not been created, so create it
331  G4double bins[1024], vals[1024], sum;
332  G4int ii;
333  G4int maxbin = G4int(ZBiasH.GetVectorLength());
334  bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
335  vals[0] = ZBiasH(size_t(0));
336  sum = vals[0];
337  for (ii = 1; ii < maxbin; ii++) {
338  bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
339  vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
340  sum = sum + ZBiasH(size_t(ii));
341  }
342 
343  for (ii = 0; ii < maxbin; ii++) {
344  vals[ii] = vals[ii] / sum;
345  IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
346  }
347  // Make IPDFZBias = true
348  IPDFZBias = true;
349  }
350  // IPDF has been create so carry on
351  G4double rndm = G4UniformRand();
352  // size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
353  size_t numberOfBin = IPDFZBiasH.GetVectorLength();
354  G4int biasn1 = 0;
355  G4int biasn2 = numberOfBin / 2;
356  G4int biasn3 = numberOfBin - 1;
357  while (biasn1 != biasn3 - 1) {
358  if (rndm > IPDFZBiasH(biasn2))
359  biasn1 = biasn2;
360  else
361  biasn3 = biasn2;
362  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
363  }
364  bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
365  G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
366  G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
367  G4double NatProb = xaxisu - xaxisl;
368  bweights[2] = NatProb / bweights[2];
369  if (verbosityLevel >= 1)
370  G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
371  return (IPDFZBiasH.GetEnergy(rndm));
372  }
373 }
374 
375 G4double G4SPSRandomGenerator::GenRandTheta() {
376  if (verbosityLevel >= 1) {
377  G4cout << "In GenRandTheta" << G4endl;
378  G4cout << "Verbosity " << verbosityLevel << G4endl;
379  }
380  if (ThetaBias == false) {
381  // Theta is not biased
382  G4double rndm = G4UniformRand();
383  return (rndm);
384  } else {
385  // Theta is biased
386  if (IPDFThetaBias == false) {
387  // IPDF has not been created, so create it
388  G4double bins[1024], vals[1024], sum;
389  G4int ii;
390  G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
391  bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
392  vals[0] = ThetaBiasH(size_t(0));
393  sum = vals[0];
394  for (ii = 1; ii < maxbin; ii++) {
395  bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
396  vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
397  sum = sum + ThetaBiasH(size_t(ii));
398  }
399 
400  for (ii = 0; ii < maxbin; ii++) {
401  vals[ii] = vals[ii] / sum;
402  IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
403  }
404  // Make IPDFThetaBias = true
405  IPDFThetaBias = true;
406  }
407  // IPDF has been create so carry on
408  G4double rndm = G4UniformRand();
409  // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
410  size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
411  G4int biasn1 = 0;
412  G4int biasn2 = numberOfBin / 2;
413  G4int biasn3 = numberOfBin - 1;
414  while (biasn1 != biasn3 - 1) {
415  if (rndm > IPDFThetaBiasH(biasn2))
416  biasn1 = biasn2;
417  else
418  biasn3 = biasn2;
419  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
420  }
421  bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
422  G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
423  G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
424  G4double NatProb = xaxisu - xaxisl;
425  bweights[3] = NatProb / bweights[3];
426  if (verbosityLevel >= 1)
427  G4cout << "Theta bin weight " << bweights[3] << " " << rndm
428  << G4endl;
429  return (IPDFThetaBiasH.GetEnergy(rndm));
430  }
431 }
432 
433 G4double G4SPSRandomGenerator::GenRandPhi() {
434  if (verbosityLevel >= 1)
435  G4cout << "In GenRandPhi" << G4endl;
436  if (PhiBias == false) {
437  // Phi is not biased
438  G4double rndm = G4UniformRand();
439  return (rndm);
440  } else {
441  // Phi is biased
442  if (IPDFPhiBias == false) {
443  // IPDF has not been created, so create it
444  G4double bins[1024], vals[1024], sum;
445  G4int ii;
446  G4int maxbin = G4int(PhiBiasH.GetVectorLength());
447  bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
448  vals[0] = PhiBiasH(size_t(0));
449  sum = vals[0];
450  for (ii = 1; ii < maxbin; ii++) {
451  bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
452  vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
453  sum = sum + PhiBiasH(size_t(ii));
454  }
455 
456  for (ii = 0; ii < maxbin; ii++) {
457  vals[ii] = vals[ii] / sum;
458  IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
459  }
460  // Make IPDFPhiBias = true
461  IPDFPhiBias = true;
462  }
463  // IPDF has been create so carry on
464  G4double rndm = G4UniformRand();
465  // size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
466  size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
467  G4int biasn1 = 0;
468  G4int biasn2 = numberOfBin / 2;
469  G4int biasn3 = numberOfBin - 1;
470  while (biasn1 != biasn3 - 1) {
471  if (rndm > IPDFPhiBiasH(biasn2))
472  biasn1 = biasn2;
473  else
474  biasn3 = biasn2;
475  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
476  }
477  bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
478  G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
479  G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
480  G4double NatProb = xaxisu - xaxisl;
481  bweights[4] = NatProb / bweights[4];
482  if (verbosityLevel >= 1)
483  G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
484  return (IPDFPhiBiasH.GetEnergy(rndm));
485  }
486 }
487 
488 G4double G4SPSRandomGenerator::GenRandEnergy() {
489  if (verbosityLevel >= 1)
490  G4cout << "In GenRandEnergy" << G4endl;
491  if (EnergyBias == false) {
492  // Energy is not biased
493  G4double rndm = G4UniformRand();
494  return (rndm);
495  } else {
496  // ENERGY is biased
497  if (IPDFEnergyBias == false) {
498  // IPDF has not been created, so create it
499  G4double bins[1024], vals[1024], sum;
500  G4int ii;
501  G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
502  bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
503  vals[0] = EnergyBiasH(size_t(0));
504  sum = vals[0];
505  for (ii = 1; ii < maxbin; ii++) {
506  bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
507  vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
508  sum = sum + EnergyBiasH(size_t(ii));
509  }
510  IPDFEnergyBiasH = ZeroPhysVector;
511  for (ii = 0; ii < maxbin; ii++) {
512  vals[ii] = vals[ii] / sum;
513  IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
514  }
515  // Make IPDFEnergyBias = true
516  IPDFEnergyBias = true;
517  }
518  // IPDF has been create so carry on
519  G4double rndm = G4UniformRand();
520  // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
521  size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
522  G4int biasn1 = 0;
523  G4int biasn2 = numberOfBin / 2;
524  G4int biasn3 = numberOfBin - 1;
525  while (biasn1 != biasn3 - 1) {
526  if (rndm > IPDFEnergyBiasH(biasn2))
527  biasn1 = biasn2;
528  else
529  biasn3 = biasn2;
530  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
531  }
532  bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
533  G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
534  G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
535  G4double NatProb = xaxisu - xaxisl;
536  bweights[5] = NatProb / bweights[5];
537  if (verbosityLevel >= 1)
538  G4cout << "Energy bin weight " << bweights[5] << " " << rndm
539  << G4endl;
540  return (IPDFEnergyBiasH.GetEnergy(rndm));
541  }
542 }
543 
544 G4double G4SPSRandomGenerator::GenRandPosTheta() {
545  if (verbosityLevel >= 1) {
546  G4cout << "In GenRandPosTheta" << G4endl;
547  G4cout << "Verbosity " << verbosityLevel << G4endl;
548  }
549  if (PosThetaBias == false) {
550  // Theta is not biased
551  G4double rndm = G4UniformRand();
552  return (rndm);
553  } else {
554  // Theta is biased
555  if (IPDFPosThetaBias == false) {
556  // IPDF has not been created, so create it
557  G4double bins[1024], vals[1024], sum;
558  G4int ii;
559  G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
560  bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
561  vals[0] = PosThetaBiasH(size_t(0));
562  sum = vals[0];
563  for (ii = 1; ii < maxbin; ii++) {
564  bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
565  vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
566  sum = sum + PosThetaBiasH(size_t(ii));
567  }
568 
569  for (ii = 0; ii < maxbin; ii++) {
570  vals[ii] = vals[ii] / sum;
571  IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
572  }
573  // Make IPDFThetaBias = true
574  IPDFPosThetaBias = true;
575  }
576  // IPDF has been create so carry on
577  G4double rndm = G4UniformRand();
578  // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
579  size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
580  G4int biasn1 = 0;
581  G4int biasn2 = numberOfBin / 2;
582  G4int biasn3 = numberOfBin - 1;
583  while (biasn1 != biasn3 - 1) {
584  if (rndm > IPDFPosThetaBiasH(biasn2))
585  biasn1 = biasn2;
586  else
587  biasn3 = biasn2;
588  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
589  }
590  bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
591  G4double xaxisl =
592  IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
593  G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
594  G4double NatProb = xaxisu - xaxisl;
595  bweights[6] = NatProb / bweights[6];
596  if (verbosityLevel >= 1)
597  G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
598  << G4endl;
599  return (IPDFPosThetaBiasH.GetEnergy(rndm));
600  }
601 }
602 
603 G4double G4SPSRandomGenerator::GenRandPosPhi() {
604  if (verbosityLevel >= 1)
605  G4cout << "In GenRandPosPhi" << G4endl;
606  if (PosPhiBias == false) {
607  // PosPhi is not biased
608  G4double rndm = G4UniformRand();
609  return (rndm);
610  } else {
611  // PosPhi is biased
612  if (IPDFPosPhiBias == false) {
613  // IPDF has not been created, so create it
614  G4double bins[1024], vals[1024], sum;
615  G4int ii;
616  G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
617  bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
618  vals[0] = PosPhiBiasH(size_t(0));
619  sum = vals[0];
620  for (ii = 1; ii < maxbin; ii++) {
621  bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
622  vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
623  sum = sum + PosPhiBiasH(size_t(ii));
624  }
625 
626  for (ii = 0; ii < maxbin; ii++) {
627  vals[ii] = vals[ii] / sum;
628  IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
629  }
630  // Make IPDFPosPhiBias = true
631  IPDFPosPhiBias = true;
632  }
633  // IPDF has been create so carry on
634  G4double rndm = G4UniformRand();
635  // size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
636  size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
637  G4int biasn1 = 0;
638  G4int biasn2 = numberOfBin / 2;
639  G4int biasn3 = numberOfBin - 1;
640  while (biasn1 != biasn3 - 1) {
641  if (rndm > IPDFPosPhiBiasH(biasn2))
642  biasn1 = biasn2;
643  else
644  biasn3 = biasn2;
645  biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
646  }
647  bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
648  G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
649  G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
650  G4double NatProb = xaxisu - xaxisl;
651  bweights[7] = NatProb / bweights[7];
652  if (verbosityLevel >= 1)
653  G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
654  << G4endl;
655  return (IPDFPosPhiBiasH.GetEnergy(rndm));
656  }
657 }
658 
Generated on Sat May 25 2013 14:33:05 for Geant4 by   doxygen 1.8.4