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