Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros 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 #include "G4AutoLock.hh"
61 
62 #include "G4SPSRandomGenerator.hh"
63 
64 //G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
65 
66 G4SPSRandomGenerator::bweights_t::bweights_t() {
67  for ( int i = 0 ; i < 9 ; ++i ) w[i] = 1;
68 }
69 
70 G4double&
71 G4SPSRandomGenerator::bweights_t::operator [](const int i) { return w[i]; }
72 
73 G4SPSRandomGenerator::G4SPSRandomGenerator()
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;
95  G4MUTEXINIT(mutex);
96 }
97 
98 G4SPSRandomGenerator::~G4SPSRandomGenerator() {
99  G4MUTEXDESTROY(mutex);
100 }
101 
102 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
103 //{
104 // if (instance == 0) instance = new G4SPSRandomGenerator();
105 // return instance;
106 //}
107 
108 // Biasing methods
109 
110 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
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 
119 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
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 
128 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
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 
137 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
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 
146 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
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 
155 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
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 
164 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
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 
173 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
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 
182 void G4SPSRandomGenerator::SetIntensityWeight(G4double weight) {
183  bweights.Get()[8] = weight;
184 }
185 
186 G4double G4SPSRandomGenerator::GetBiasWeight() {
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 
193 void G4SPSRandomGenerator::SetVerbosity(G4int a) {
194  G4AutoLock l(&mutex);
195  verbosityLevel = a;
196 }
197 
198 
199 namespace {
200  G4PhysicsOrderedFreeVector ZeroPhysVector; // for re-set only
201 }
202 
203 void G4SPSRandomGenerator::ReSetHist(G4String atype) {
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 
250 G4double G4SPSRandomGenerator::GenRandX() {
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 
327 G4double G4SPSRandomGenerator::GenRandY() {
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 
385 G4double G4SPSRandomGenerator::GenRandZ() {
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 
445 G4double G4SPSRandomGenerator::GenRandTheta() {
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;
463  G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
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 
508 G4double G4SPSRandomGenerator::GenRandPhi() {
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 
568 G4double G4SPSRandomGenerator::GenRandEnergy() {
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;
584  G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
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 
629 G4double G4SPSRandomGenerator::GenRandPosTheta() {
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;
647  G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
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 
693 G4double G4SPSRandomGenerator::GenRandPosPhi() {
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;
709  G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
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 }
CLHEP::Hep3Vector::x
double x() const
G4Track.hh
G4Cache::Get
value_type & Get() const
Definition: G4Cache.hh:282
G4MUTEXINIT
#define G4MUTEXINIT(mutex)
Definition: G4Threading.hh:177
G4PhysicsOrderedFreeVector::InsertValues
void InsertValues(G4double energy, G4double value)
Definition: G4PhysicsOrderedFreeVector.cc:79
G4SPSRandomGenerator::SetEnergyBias
void SetEnergyBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:155
G4Ions.hh
G4SPSRandomGenerator::SetYBias
void SetYBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:119
G4PhysicsVector::GetVectorLength
size_t GetVectorLength() const
G4PhysicsVector::GetLowEdgeEnergy
G4double GetLowEdgeEnergy(size_t binNumber) const
Definition: G4PhysicsVector.cc:151
G4int
int G4int
Definition: G4Types.hh:78
G4SPSRandomGenerator::SetPosThetaBias
void SetPosThetaBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:164
G4SPSRandomGenerator::SetThetaBias
void SetThetaBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:137
G4IonTable.hh
G4SPSRandomGenerator::SetZBias
void SetZBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:128
G4SPSRandomGenerator::SetXBias
void SetXBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:110
G4PrimaryParticle.hh
G4SPSRandomGenerator.hh
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4cout
G4GLOB_DLL std::ostream G4cout
G4ParticleTable.hh
G4ParticleDefinition.hh
Randomize.hh
G4TransportationManager.hh
G4PhysicsOrderedFreeVector::GetEnergy
G4double GetEnergy(G4double aValue)
Definition: G4PhysicsOrderedFreeVector.cc:96
G4SPSRandomGenerator::SetPhiBias
void SetPhiBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:146
G4PhysicalVolumeStore.hh
G4SPSRandomGenerator::SetIntensityWeight
void SetIntensityWeight(G4double weight)
Definition: G4SPSRandomGenerator.cc:182
CLHEP::Hep3Vector::y
double y() const
G4endl
#define G4endl
Definition: G4ios.hh:61
G4VPhysicalVolume.hh
G4double
double G4double
Definition: G4Types.hh:76
G4TrackingManager.hh
G4Event.hh
G4MUTEXDESTROY
#define G4MUTEXDESTROY(mutex)
Definition: G4Threading.hh:178
G4AutoLock.hh
G4SPSRandomGenerator::SetPosPhiBias
void SetPosPhiBias(G4ThreeVector)
Definition: G4SPSRandomGenerator.cc:173