| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4SPSRandomGenerator class implementation << 26 ///////////////////////////////////////////////////////////////////////////////
>> 27 //
>> 28 // MODULE: G4SPSRandomGenerator.cc
>> 29 //
>> 30 // Version: 1.0
>> 31 // Date: 5/02/04
>> 32 // Author: Fan Lei
>> 33 // Organisation: QinetiQ ltd.
>> 34 // Customer: ESA/ESTEC
>> 35 //
>> 36 ///////////////////////////////////////////////////////////////////////////////
>> 37 //
>> 38 // CHANGE HISTORY
>> 39 // --------------
>> 40 //
>> 41 //
>> 42 // Version 1.0, 05/02/2004, Fan Lei, Created.
>> 43 // Based on the G4GeneralParticleSource class in Geant4 v6.0
>> 44 //
>> 45 ///////////////////////////////////////////////////////////////////////////////
27 // 46 //
28 // Author: Fan Lei, QinetiQ ltd. - 05/02/2004 <<
29 // Customer: ESA/ESTEC <<
30 // Revision: Andrea Dotti, SLAC <<
31 // ------------------------------------------- <<
32 <<
33 #include <cmath> <<
34 <<
35 #include "G4PrimaryParticle.hh" 47 #include "G4PrimaryParticle.hh"
36 #include "G4Event.hh" 48 #include "G4Event.hh"
37 #include "Randomize.hh" 49 #include "Randomize.hh"
>> 50 #include <cmath>
38 #include "G4TransportationManager.hh" 51 #include "G4TransportationManager.hh"
39 #include "G4VPhysicalVolume.hh" 52 #include "G4VPhysicalVolume.hh"
40 #include "G4PhysicalVolumeStore.hh" 53 #include "G4PhysicalVolumeStore.hh"
41 #include "G4ParticleTable.hh" 54 #include "G4ParticleTable.hh"
42 #include "G4ParticleDefinition.hh" 55 #include "G4ParticleDefinition.hh"
43 #include "G4IonTable.hh" 56 #include "G4IonTable.hh"
44 #include "G4Ions.hh" 57 #include "G4Ions.hh"
45 #include "G4TrackingManager.hh" 58 #include "G4TrackingManager.hh"
46 #include "G4Track.hh" 59 #include "G4Track.hh"
47 #include "G4AutoLock.hh" <<
48 60
49 #include "G4SPSRandomGenerator.hh" 61 #include "G4SPSRandomGenerator.hh"
50 62
51 G4SPSRandomGenerator::bweights_t::bweights_t() << 63 //G4SPSRandomGenerator* G4SPSRandomGenerator::instance = 0;
52 { <<
53 for (double & i : w) { i = 1; } <<
54 } <<
55 <<
56 G4double& G4SPSRandomGenerator::bweights_t::op <<
57 { <<
58 return w[i]; <<
59 } <<
60 64
61 G4SPSRandomGenerator::G4SPSRandomGenerator() 65 G4SPSRandomGenerator::G4SPSRandomGenerator()
>> 66 : alpha(0.)
62 { 67 {
63 // Initialise all variables << 68 // Initialise all variables
64 69
65 // Bias variables << 70 // Bias variables
66 // << 71 XBias = false;
67 XBias = false; << 72 IPDFXBias = false;
68 IPDFXBias = false; << 73 YBias = false;
69 YBias = false; << 74 IPDFYBias = false;
70 IPDFYBias = false; << 75 ZBias = false;
71 ZBias = false; << 76 IPDFZBias = false;
72 IPDFZBias = false; << 77 ThetaBias = false;
73 ThetaBias = false; << 78 IPDFThetaBias = false;
74 IPDFThetaBias = false; << 79 PhiBias = false;
75 PhiBias = false; << 80 IPDFPhiBias = false;
76 IPDFPhiBias = false; << 81 EnergyBias = false;
77 EnergyBias = false; << 82 IPDFEnergyBias = false;
78 IPDFEnergyBias = false; << 83 PosThetaBias = false;
79 PosThetaBias = false; << 84 IPDFPosThetaBias = false;
80 IPDFPosThetaBias = false; << 85 PosPhiBias = false;
81 PosPhiBias = false; << 86 IPDFPosPhiBias = false;
82 IPDFPosPhiBias = false; << 87 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4]
83 << 88 = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
84 verbosityLevel = 0; << 89 verbosityLevel = 0;
85 G4MUTEXINIT(mutex); << 90 }
86 } << 91
87 << 92 G4SPSRandomGenerator::~G4SPSRandomGenerator() {
88 G4SPSRandomGenerator::~G4SPSRandomGenerator() << 93 }
89 { << 94
90 G4MUTEXDESTROY(mutex); << 95 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance ()
91 } << 96 //{
>> 97 // if (instance == 0) instance = new G4SPSRandomGenerator();
>> 98 // return instance;
>> 99 //}
92 100
93 // Biasing methods 101 // Biasing methods
94 102
95 void G4SPSRandomGenerator::SetXBias(const G4Th << 103 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
96 { << 104 G4double ehi, val;
97 G4AutoLock l(&mutex); << 105 ehi = input.x();
98 G4double ehi, val; << 106 val = input.y();
99 ehi = input.x(); << 107 XBiasH.InsertValues(ehi, val);
100 val = input.y(); << 108 XBias = true;
101 XBiasH.InsertValues(ehi, val); << 109 }
102 XBias = true; << 110
103 } << 111 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
104 << 112 G4double ehi, val;
105 void G4SPSRandomGenerator::SetYBias(const G4Th << 113 ehi = input.x();
106 { << 114 val = input.y();
107 G4AutoLock l(&mutex); << 115 YBiasH.InsertValues(ehi, val);
108 G4double ehi, val; << 116 YBias = true;
109 ehi = input.x(); << 117 }
110 val = input.y(); << 118
111 YBiasH.InsertValues(ehi, val); << 119 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
112 YBias = true; << 120 G4double ehi, val;
113 } << 121 ehi = input.x();
114 << 122 val = input.y();
115 void G4SPSRandomGenerator::SetZBias(const G4Th << 123 ZBiasH.InsertValues(ehi, val);
116 { << 124 ZBias = true;
117 G4AutoLock l(&mutex); << 125 }
118 G4double ehi, val; << 126
119 ehi = input.x(); << 127 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
120 val = input.y(); << 128 G4double ehi, val;
121 ZBiasH.InsertValues(ehi, val); << 129 ehi = input.x();
122 ZBias = true; << 130 val = input.y();
123 } << 131 ThetaBiasH.InsertValues(ehi, val);
124 << 132 ThetaBias = true;
125 void G4SPSRandomGenerator::SetThetaBias(const << 133 }
126 { << 134
127 G4AutoLock l(&mutex); << 135 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
128 G4double ehi, val; << 136 G4double ehi, val;
129 ehi = input.x(); << 137 ehi = input.x();
130 val = input.y(); << 138 val = input.y();
131 ThetaBiasH.InsertValues(ehi, val); << 139 PhiBiasH.InsertValues(ehi, val);
132 ThetaBias = true; << 140 PhiBias = true;
>> 141 }
>> 142
>> 143 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
>> 144 G4double ehi, val;
>> 145 ehi = input.x();
>> 146 val = input.y();
>> 147 EnergyBiasH.InsertValues(ehi, val);
>> 148 EnergyBias = true;
>> 149 }
>> 150
>> 151 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
>> 152 G4double ehi, val;
>> 153 ehi = input.x();
>> 154 val = input.y();
>> 155 PosThetaBiasH.InsertValues(ehi, val);
>> 156 PosThetaBias = true;
>> 157 }
>> 158
>> 159 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
>> 160 G4double ehi, val;
>> 161 ehi = input.x();
>> 162 val = input.y();
>> 163 PosPhiBiasH.InsertValues(ehi, val);
>> 164 PosPhiBias = true;
>> 165 }
>> 166
>> 167 void G4SPSRandomGenerator::ReSetHist(G4String atype) {
>> 168 if (atype == "biasx") {
>> 169 XBias = false;
>> 170 IPDFXBias = false;
>> 171 XBiasH = IPDFXBiasH = ZeroPhysVector;
>> 172 } else if (atype == "biasy") {
>> 173 YBias = false;
>> 174 IPDFYBias = false;
>> 175 YBiasH = IPDFYBiasH = ZeroPhysVector;
>> 176 } else if (atype == "biasz") {
>> 177 ZBias = false;
>> 178 IPDFZBias = false;
>> 179 ZBiasH = IPDFZBiasH = ZeroPhysVector;
>> 180 } else if (atype == "biast") {
>> 181 ThetaBias = false;
>> 182 IPDFThetaBias = false;
>> 183 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
>> 184 } else if (atype == "biasp") {
>> 185 PhiBias = false;
>> 186 IPDFPhiBias = false;
>> 187 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
>> 188 } else if (atype == "biase") {
>> 189 EnergyBias = false;
>> 190 IPDFEnergyBias = false;
>> 191 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
>> 192 } else if (atype == "biaspt") {
>> 193 PosThetaBias = false;
>> 194 IPDFPosThetaBias = false;
>> 195 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
>> 196 } else if (atype == "biaspp") {
>> 197 PosPhiBias = false;
>> 198 IPDFPosPhiBias = false;
>> 199 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
>> 200 } else {
>> 201 G4cout << "Error, histtype not accepted " << G4endl;
>> 202 }
>> 203 }
>> 204
>> 205 G4double G4SPSRandomGenerator::GenRandX() {
>> 206 if (verbosityLevel >= 1)
>> 207 G4cout << "In GenRandX" << G4endl;
>> 208 if (XBias == false) {
>> 209 // X is not biased
>> 210 G4double rndm = G4UniformRand();
>> 211 return (rndm);
>> 212 } else {
>> 213 // X is biased
>> 214 if (IPDFXBias == false) {
>> 215 // IPDF has not been created, so create it
>> 216 G4double bins[1024], vals[1024], sum;
>> 217 G4int ii;
>> 218 G4int maxbin = G4int(XBiasH.GetVectorLength());
>> 219 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
>> 220 vals[0] = XBiasH(size_t(0));
>> 221 sum = vals[0];
>> 222 for (ii = 1; ii < maxbin; ii++) {
>> 223 bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
>> 224 vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
>> 225 sum = sum + XBiasH(size_t(ii));
>> 226 }
>> 227
>> 228 for (ii = 0; ii < maxbin; ii++) {
>> 229 vals[ii] = vals[ii] / sum;
>> 230 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
>> 231 }
>> 232 // Make IPDFXBias = true
>> 233 IPDFXBias = true;
>> 234 }
>> 235 // IPDF has been create so carry on
>> 236 G4double rndm = G4UniformRand();
>> 237
>> 238 // Calculate the weighting: Find the bin that the determined
>> 239 // rndm is in and the weigthing will be the difference in the
>> 240 // natural probability (from the x-axis) divided by the
>> 241 // difference in the biased probability (the area).
>> 242 size_t numberOfBin = IPDFXBiasH.GetVectorLength();
>> 243 G4int biasn1 = 0;
>> 244 G4int biasn2 = numberOfBin / 2;
>> 245 G4int biasn3 = numberOfBin - 1;
>> 246 while (biasn1 != biasn3 - 1) {
>> 247 if (rndm > IPDFXBiasH(biasn2))
>> 248 biasn1 = biasn2;
>> 249 else
>> 250 biasn3 = biasn2;
>> 251 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 252 }
>> 253 // retrieve the areas and then the x-axis values
>> 254 bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
>> 255 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 256 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 257 G4double NatProb = xaxisu - xaxisl;
>> 258 //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl;
>> 259 //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl;
>> 260 bweights[0] = NatProb / bweights[0];
>> 261 if (verbosityLevel >= 1)
>> 262 G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
>> 263 return (IPDFXBiasH.GetEnergy(rndm));
>> 264 }
>> 265 }
>> 266
>> 267 G4double G4SPSRandomGenerator::GenRandY() {
>> 268 if (verbosityLevel >= 1)
>> 269 G4cout << "In GenRandY" << G4endl;
>> 270 if (YBias == false) {
>> 271 // Y is not biased
>> 272 G4double rndm = G4UniformRand();
>> 273 return (rndm);
>> 274 } else {
>> 275 // Y is biased
>> 276 if (IPDFYBias == false) {
>> 277 // IPDF has not been created, so create it
>> 278 G4double bins[1024], vals[1024], sum;
>> 279 G4int ii;
>> 280 G4int maxbin = G4int(YBiasH.GetVectorLength());
>> 281 bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
>> 282 vals[0] = YBiasH(size_t(0));
>> 283 sum = vals[0];
>> 284 for (ii = 1; ii < maxbin; ii++) {
>> 285 bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
>> 286 vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
>> 287 sum = sum + YBiasH(size_t(ii));
>> 288 }
>> 289
>> 290 for (ii = 0; ii < maxbin; ii++) {
>> 291 vals[ii] = vals[ii] / sum;
>> 292 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
>> 293 }
>> 294 // Make IPDFYBias = true
>> 295 IPDFYBias = true;
>> 296 }
>> 297 // IPDF has been create so carry on
>> 298 G4double rndm = G4UniformRand();
>> 299 size_t numberOfBin = IPDFYBiasH.GetVectorLength();
>> 300 G4int biasn1 = 0;
>> 301 G4int biasn2 = numberOfBin / 2;
>> 302 G4int biasn3 = numberOfBin - 1;
>> 303 while (biasn1 != biasn3 - 1) {
>> 304 if (rndm > IPDFYBiasH(biasn2))
>> 305 biasn1 = biasn2;
>> 306 else
>> 307 biasn3 = biasn2;
>> 308 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 309 }
>> 310 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
>> 311 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 312 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 313 G4double NatProb = xaxisu - xaxisl;
>> 314 bweights[1] = NatProb / bweights[1];
>> 315 if (verbosityLevel >= 1)
>> 316 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
>> 317 return (IPDFYBiasH.GetEnergy(rndm));
>> 318 }
>> 319 }
>> 320
>> 321 G4double G4SPSRandomGenerator::GenRandZ() {
>> 322 if (verbosityLevel >= 1)
>> 323 G4cout << "In GenRandZ" << G4endl;
>> 324 if (ZBias == false) {
>> 325 // Z is not biased
>> 326 G4double rndm = G4UniformRand();
>> 327 return (rndm);
>> 328 } else {
>> 329 // Z is biased
>> 330 if (IPDFZBias == false) {
>> 331 // IPDF has not been created, so create it
>> 332 G4double bins[1024], vals[1024], sum;
>> 333 G4int ii;
>> 334 G4int maxbin = G4int(ZBiasH.GetVectorLength());
>> 335 bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
>> 336 vals[0] = ZBiasH(size_t(0));
>> 337 sum = vals[0];
>> 338 for (ii = 1; ii < maxbin; ii++) {
>> 339 bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
>> 340 vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
>> 341 sum = sum + ZBiasH(size_t(ii));
>> 342 }
>> 343
>> 344 for (ii = 0; ii < maxbin; ii++) {
>> 345 vals[ii] = vals[ii] / sum;
>> 346 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
>> 347 }
>> 348 // Make IPDFZBias = true
>> 349 IPDFZBias = true;
>> 350 }
>> 351 // IPDF has been create so carry on
>> 352 G4double rndm = G4UniformRand();
>> 353 // size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm);
>> 354 size_t numberOfBin = IPDFZBiasH.GetVectorLength();
>> 355 G4int biasn1 = 0;
>> 356 G4int biasn2 = numberOfBin / 2;
>> 357 G4int biasn3 = numberOfBin - 1;
>> 358 while (biasn1 != biasn3 - 1) {
>> 359 if (rndm > IPDFZBiasH(biasn2))
>> 360 biasn1 = biasn2;
>> 361 else
>> 362 biasn3 = biasn2;
>> 363 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 364 }
>> 365 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
>> 366 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 367 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 368 G4double NatProb = xaxisu - xaxisl;
>> 369 bweights[2] = NatProb / bweights[2];
>> 370 if (verbosityLevel >= 1)
>> 371 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
>> 372 return (IPDFZBiasH.GetEnergy(rndm));
>> 373 }
>> 374 }
>> 375
>> 376 G4double G4SPSRandomGenerator::GenRandTheta() {
>> 377 if (verbosityLevel >= 1) {
>> 378 G4cout << "In GenRandTheta" << G4endl;
>> 379 G4cout << "Verbosity " << verbosityLevel << G4endl;
>> 380 }
>> 381 if (ThetaBias == false) {
>> 382 // Theta is not biased
>> 383 G4double rndm = G4UniformRand();
>> 384 return (rndm);
>> 385 } else {
>> 386 // Theta is biased
>> 387 if (IPDFThetaBias == false) {
>> 388 // IPDF has not been created, so create it
>> 389 G4double bins[1024], vals[1024], sum;
>> 390 G4int ii;
>> 391 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
>> 392 bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
>> 393 vals[0] = ThetaBiasH(size_t(0));
>> 394 sum = vals[0];
>> 395 for (ii = 1; ii < maxbin; ii++) {
>> 396 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
>> 397 vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
>> 398 sum = sum + ThetaBiasH(size_t(ii));
>> 399 }
>> 400
>> 401 for (ii = 0; ii < maxbin; ii++) {
>> 402 vals[ii] = vals[ii] / sum;
>> 403 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
>> 404 }
>> 405 // Make IPDFThetaBias = true
>> 406 IPDFThetaBias = true;
>> 407 }
>> 408 // IPDF has been create so carry on
>> 409 G4double rndm = G4UniformRand();
>> 410 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
>> 411 size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
>> 412 G4int biasn1 = 0;
>> 413 G4int biasn2 = numberOfBin / 2;
>> 414 G4int biasn3 = numberOfBin - 1;
>> 415 while (biasn1 != biasn3 - 1) {
>> 416 if (rndm > IPDFThetaBiasH(biasn2))
>> 417 biasn1 = biasn2;
>> 418 else
>> 419 biasn3 = biasn2;
>> 420 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 421 }
>> 422 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
>> 423 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 424 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 425 G4double NatProb = xaxisu - xaxisl;
>> 426 bweights[3] = NatProb / bweights[3];
>> 427 if (verbosityLevel >= 1)
>> 428 G4cout << "Theta bin weight " << bweights[3] << " " << rndm
>> 429 << G4endl;
>> 430 return (IPDFThetaBiasH.GetEnergy(rndm));
>> 431 }
>> 432 }
>> 433
>> 434 G4double G4SPSRandomGenerator::GenRandPhi() {
>> 435 if (verbosityLevel >= 1)
>> 436 G4cout << "In GenRandPhi" << G4endl;
>> 437 if (PhiBias == false) {
>> 438 // Phi is not biased
>> 439 G4double rndm = G4UniformRand();
>> 440 return (rndm);
>> 441 } else {
>> 442 // Phi is biased
>> 443 if (IPDFPhiBias == false) {
>> 444 // IPDF has not been created, so create it
>> 445 G4double bins[1024], vals[1024], sum;
>> 446 G4int ii;
>> 447 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
>> 448 bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
>> 449 vals[0] = PhiBiasH(size_t(0));
>> 450 sum = vals[0];
>> 451 for (ii = 1; ii < maxbin; ii++) {
>> 452 bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
>> 453 vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
>> 454 sum = sum + PhiBiasH(size_t(ii));
>> 455 }
>> 456
>> 457 for (ii = 0; ii < maxbin; ii++) {
>> 458 vals[ii] = vals[ii] / sum;
>> 459 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
>> 460 }
>> 461 // Make IPDFPhiBias = true
>> 462 IPDFPhiBias = true;
>> 463 }
>> 464 // IPDF has been create so carry on
>> 465 G4double rndm = G4UniformRand();
>> 466 // size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm);
>> 467 size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
>> 468 G4int biasn1 = 0;
>> 469 G4int biasn2 = numberOfBin / 2;
>> 470 G4int biasn3 = numberOfBin - 1;
>> 471 while (biasn1 != biasn3 - 1) {
>> 472 if (rndm > IPDFPhiBiasH(biasn2))
>> 473 biasn1 = biasn2;
>> 474 else
>> 475 biasn3 = biasn2;
>> 476 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 477 }
>> 478 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
>> 479 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 480 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 481 G4double NatProb = xaxisu - xaxisl;
>> 482 bweights[4] = NatProb / bweights[4];
>> 483 if (verbosityLevel >= 1)
>> 484 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
>> 485 return (IPDFPhiBiasH.GetEnergy(rndm));
>> 486 }
>> 487 }
>> 488
>> 489 G4double G4SPSRandomGenerator::GenRandEnergy() {
>> 490 if (verbosityLevel >= 1)
>> 491 G4cout << "In GenRandEnergy" << G4endl;
>> 492 if (EnergyBias == false) {
>> 493 // Energy is not biased
>> 494 G4double rndm = G4UniformRand();
>> 495 return (rndm);
>> 496 } else {
>> 497 // ENERGY is biased
>> 498 if (IPDFEnergyBias == false) {
>> 499 // IPDF has not been created, so create it
>> 500 G4double bins[1024], vals[1024], sum;
>> 501 G4int ii;
>> 502 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
>> 503 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
>> 504 vals[0] = EnergyBiasH(size_t(0));
>> 505 sum = vals[0];
>> 506 for (ii = 1; ii < maxbin; ii++) {
>> 507 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
>> 508 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
>> 509 sum = sum + EnergyBiasH(size_t(ii));
>> 510 }
>> 511 IPDFEnergyBiasH = ZeroPhysVector;
>> 512 for (ii = 0; ii < maxbin; ii++) {
>> 513 vals[ii] = vals[ii] / sum;
>> 514 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
>> 515 }
>> 516 // Make IPDFEnergyBias = true
>> 517 IPDFEnergyBias = true;
>> 518 }
>> 519 // IPDF has been create so carry on
>> 520 G4double rndm = G4UniformRand();
>> 521 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm);
>> 522 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
>> 523 G4int biasn1 = 0;
>> 524 G4int biasn2 = numberOfBin / 2;
>> 525 G4int biasn3 = numberOfBin - 1;
>> 526 while (biasn1 != biasn3 - 1) {
>> 527 if (rndm > IPDFEnergyBiasH(biasn2))
>> 528 biasn1 = biasn2;
>> 529 else
>> 530 biasn3 = biasn2;
>> 531 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 532 }
>> 533 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
>> 534 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 535 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 536 G4double NatProb = xaxisu - xaxisl;
>> 537 bweights[5] = NatProb / bweights[5];
>> 538 if (verbosityLevel >= 1)
>> 539 G4cout << "Energy bin weight " << bweights[5] << " " << rndm
>> 540 << G4endl;
>> 541 return (IPDFEnergyBiasH.GetEnergy(rndm));
>> 542 }
>> 543 }
>> 544
>> 545 G4double G4SPSRandomGenerator::GenRandPosTheta() {
>> 546 if (verbosityLevel >= 1) {
>> 547 G4cout << "In GenRandPosTheta" << G4endl;
>> 548 G4cout << "Verbosity " << verbosityLevel << G4endl;
>> 549 }
>> 550 if (PosThetaBias == false) {
>> 551 // Theta is not biased
>> 552 G4double rndm = G4UniformRand();
>> 553 return (rndm);
>> 554 } else {
>> 555 // Theta is biased
>> 556 if (IPDFPosThetaBias == false) {
>> 557 // IPDF has not been created, so create it
>> 558 G4double bins[1024], vals[1024], sum;
>> 559 G4int ii;
>> 560 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
>> 561 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
>> 562 vals[0] = PosThetaBiasH(size_t(0));
>> 563 sum = vals[0];
>> 564 for (ii = 1; ii < maxbin; ii++) {
>> 565 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
>> 566 vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
>> 567 sum = sum + PosThetaBiasH(size_t(ii));
>> 568 }
>> 569
>> 570 for (ii = 0; ii < maxbin; ii++) {
>> 571 vals[ii] = vals[ii] / sum;
>> 572 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
>> 573 }
>> 574 // Make IPDFThetaBias = true
>> 575 IPDFPosThetaBias = true;
>> 576 }
>> 577 // IPDF has been create so carry on
>> 578 G4double rndm = G4UniformRand();
>> 579 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm);
>> 580 size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
>> 581 G4int biasn1 = 0;
>> 582 G4int biasn2 = numberOfBin / 2;
>> 583 G4int biasn3 = numberOfBin - 1;
>> 584 while (biasn1 != biasn3 - 1) {
>> 585 if (rndm > IPDFPosThetaBiasH(biasn2))
>> 586 biasn1 = biasn2;
>> 587 else
>> 588 biasn3 = biasn2;
>> 589 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 590 }
>> 591 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
>> 592 G4double xaxisl =
>> 593 IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 594 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 595 G4double NatProb = xaxisu - xaxisl;
>> 596 bweights[6] = NatProb / bweights[6];
>> 597 if (verbosityLevel >= 1)
>> 598 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
>> 599 << G4endl;
>> 600 return (IPDFPosThetaBiasH.GetEnergy(rndm));
>> 601 }
>> 602 }
>> 603
>> 604 G4double G4SPSRandomGenerator::GenRandPosPhi() {
>> 605 if (verbosityLevel >= 1)
>> 606 G4cout << "In GenRandPosPhi" << G4endl;
>> 607 if (PosPhiBias == false) {
>> 608 // PosPhi is not biased
>> 609 G4double rndm = G4UniformRand();
>> 610 return (rndm);
>> 611 } else {
>> 612 // PosPhi is biased
>> 613 if (IPDFPosPhiBias == false) {
>> 614 // IPDF has not been created, so create it
>> 615 G4double bins[1024], vals[1024], sum;
>> 616 G4int ii;
>> 617 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
>> 618 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
>> 619 vals[0] = PosPhiBiasH(size_t(0));
>> 620 sum = vals[0];
>> 621 for (ii = 1; ii < maxbin; ii++) {
>> 622 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
>> 623 vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
>> 624 sum = sum + PosPhiBiasH(size_t(ii));
>> 625 }
>> 626
>> 627 for (ii = 0; ii < maxbin; ii++) {
>> 628 vals[ii] = vals[ii] / sum;
>> 629 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
>> 630 }
>> 631 // Make IPDFPosPhiBias = true
>> 632 IPDFPosPhiBias = true;
>> 633 }
>> 634 // IPDF has been create so carry on
>> 635 G4double rndm = G4UniformRand();
>> 636 // size_t weight_bin_no = IPDFPosPhiBiasH.FindValueBinLocation(rndm);
>> 637 size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
>> 638 G4int biasn1 = 0;
>> 639 G4int biasn2 = numberOfBin / 2;
>> 640 G4int biasn3 = numberOfBin - 1;
>> 641 while (biasn1 != biasn3 - 1) {
>> 642 if (rndm > IPDFPosPhiBiasH(biasn2))
>> 643 biasn1 = biasn2;
>> 644 else
>> 645 biasn3 = biasn2;
>> 646 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
>> 647 }
>> 648 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
>> 649 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
>> 650 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
>> 651 G4double NatProb = xaxisu - xaxisl;
>> 652 bweights[7] = NatProb / bweights[7];
>> 653 if (verbosityLevel >= 1)
>> 654 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
>> 655 << G4endl;
>> 656 return (IPDFPosPhiBiasH.GetEnergy(rndm));
>> 657 }
133 } 658 }
134 659
135 void G4SPSRandomGenerator::SetPhiBias(const G4 <<
136 { <<
137 G4AutoLock l(&mutex); <<
138 G4double ehi, val; <<
139 ehi = input.x(); <<
140 val = input.y(); <<
141 PhiBiasH.InsertValues(ehi, val); <<
142 PhiBias = true; <<
143 } <<
144 <<
145 void G4SPSRandomGenerator::SetEnergyBias(const <<
146 { <<
147 G4AutoLock l(&mutex); <<
148 G4double ehi, val; <<
149 ehi = input.x(); <<
150 val = input.y(); <<
151 EnergyBiasH.InsertValues(ehi, val); <<
152 EnergyBias = true; <<
153 } <<
154 <<
155 void G4SPSRandomGenerator::SetPosThetaBias(con <<
156 { <<
157 G4AutoLock l(&mutex); <<
158 G4double ehi, val; <<
159 ehi = input.x(); <<
160 val = input.y(); <<
161 PosThetaBiasH.InsertValues(ehi, val); <<
162 PosThetaBias = true; <<
163 } <<
164 <<
165 void G4SPSRandomGenerator::SetPosPhiBias(const <<
166 { <<
167 G4AutoLock l(&mutex); <<
168 G4double ehi, val; <<
169 ehi = input.x(); <<
170 val = input.y(); <<
171 PosPhiBiasH.InsertValues(ehi, val); <<
172 PosPhiBias = true; <<
173 } <<
174 <<
175 void G4SPSRandomGenerator::SetIntensityWeight( <<
176 { <<
177 bweights.Get()[8] = weight; <<
178 } <<
179 <<
180 G4double G4SPSRandomGenerator::GetBiasWeight() <<
181 { <<
182 bweights_t& w = bweights.Get(); <<
183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[ <<
184 } <<
185 <<
186 void G4SPSRandomGenerator::SetVerbosity(G4int <<
187 { <<
188 G4AutoLock l(&mutex); <<
189 verbosityLevel = a; <<
190 } <<
191 <<
192 namespace <<
193 { <<
194 G4PhysicsFreeVector ZeroPhysVector; // for r <<
195 } <<
196 <<
197 void G4SPSRandomGenerator::ReSetHist(const G4S <<
198 { <<
199 G4AutoLock l(&mutex); <<
200 if (atype == "biasx") { <<
201 XBias = false; <<
202 IPDFXBias = false; <<
203 local_IPDFXBias.Get().val = fa <<
204 XBiasH = IPDFXBiasH = ZeroPhys <<
205 } else if (atype == "biasy") { <<
206 YBias = false; <<
207 IPDFYBias = false; <<
208 local_IPDFYBias.Get().val = fa <<
209 YBiasH = IPDFYBiasH = ZeroPhys <<
210 } else if (atype == "biasz") { <<
211 ZBias = false; <<
212 IPDFZBias = false; <<
213 local_IPDFZBias.Get().val = fa <<
214 ZBiasH = IPDFZBiasH = ZeroPhys <<
215 } else if (atype == "biast") { <<
216 ThetaBias = false; <<
217 IPDFThetaBias = false; <<
218 local_IPDFThetaBias.Get().val <<
219 ThetaBiasH = IPDFThetaBiasH = <<
220 } else if (atype == "biasp") { <<
221 PhiBias = false; <<
222 IPDFPhiBias = false; <<
223 local_IPDFPhiBias.Get().val = <<
224 PhiBiasH = IPDFPhiBiasH = Zero <<
225 } else if (atype == "biase") { <<
226 EnergyBias = false; <<
227 IPDFEnergyBias = false; <<
228 local_IPDFEnergyBias.Get().val <<
229 EnergyBiasH = IPDFEnergyBiasH <<
230 } else if (atype == "biaspt") { <<
231 PosThetaBias = false; <<
232 IPDFPosThetaBias = false; <<
233 local_IPDFPosThetaBias.Get().v <<
234 PosThetaBiasH = IPDFPosThetaBi <<
235 } else if (atype == "biaspp") { <<
236 PosPhiBias = false; <<
237 IPDFPosPhiBias = false; <<
238 local_IPDFPosPhiBias.Get().val <<
239 PosPhiBiasH = IPDFPosPhiBiasH <<
240 } else { <<
241 G4cout << "Error, histtype not <<
242 } <<
243 } <<
244 <<
245 G4double G4SPSRandomGenerator::GenRandX() <<
246 { <<
247 if (verbosityLevel >= 1) <<
248 G4cout << "In GenRandX" << G4endl; <<
249 if (!XBias) <<
250 { <<
251 // X is not biased <<
252 G4double rndm = G4UniformRand(); <<
253 return (rndm); <<
254 } <<
255 <<
256 // X is biased <<
257 // This is shared among threads, and we need <<
258 // only once. Multiple instances of this cla <<
259 // so we rely on a class-private, thread-pri <<
260 // to check if we need an initialiation. We <<
261 // because the Boolean on which we check is <<
262 // <<
263 if ( !local_IPDFXBias.Get().val ) <<
264 { <<
265 // For time that this thread arrived, here <<
266 // Now two cases are possible: it is the f <<
267 // ANY thread has ever initialized this. <<
268 // Now we need a lock. In any case, the th <<
269 // variable can now be set to true <<
270 // <<
271 local_IPDFXBias.Get().val = true; <<
272 G4AutoLock l(&mutex); <<
273 if (!IPDFXBias) <<
274 { <<
275 // IPDF has not been created, so create <<
276 // <<
277 G4double bins[1024], vals[1024], sum; <<
278 std::size_t ii; <<
279 std::size_t maxbin = XBiasH.GetVectorLen <<
280 bins[0] = XBiasH.GetLowEdgeEnergy(0); <<
281 vals[0] = XBiasH(0); <<
282 sum = vals[0]; <<
283 for (ii=1; ii<maxbin; ++ii) <<
284 { <<
285 bins[ii] = XBiasH.GetLowEdgeEnergy(ii) <<
286 vals[ii] = XBiasH(ii) + vals[ii - 1]; <<
287 sum = sum + XBiasH(ii); <<
288 } <<
289 <<
290 for (ii=0; ii<maxbin; ++ii) <<
291 { <<
292 vals[ii] = vals[ii] / sum; <<
293 IPDFXBiasH.InsertValues(bins[ii], vals <<
294 } <<
295 IPDFXBias = true; <<
296 } <<
297 } <<
298 <<
299 // IPDF has been create so carry on <<
300 <<
301 G4double rndm = G4UniformRand(); <<
302 <<
303 // Calculate the weighting: Find the bin tha <<
304 // rndm is in and the weigthing will be the <<
305 // natural probability (from the x-axis) div <<
306 // difference in the biased probability (the <<
307 // <<
308 std::size_t numberOfBin = IPDFXBiasH.GetVect <<
309 std::size_t biasn1 = 0; <<
310 std::size_t biasn2 = numberOfBin / 2; <<
311 std::size_t biasn3 = numberOfBin - 1; <<
312 while (biasn1 != biasn3 - 1) <<
313 { <<
314 if (rndm > IPDFXBiasH(biasn2)) <<
315 { biasn1 = biasn2; } <<
316 else <<
317 { biasn3 = biasn2; } <<
318 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / <<
319 } <<
320 <<
321 // Retrieve the areas and then the x-axis va <<
322 // <<
323 bweights_t& w = bweights.Get(); <<
324 w[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn <<
325 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnerg <<
326 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnerg <<
327 G4double NatProb = xaxisu - xaxisl; <<
328 w[0] = NatProb / w[0]; <<
329 if (verbosityLevel >= 1) <<
330 { <<
331 G4cout << "X bin weight " << w[0] << " " < <<
332 } <<
333 return (IPDFXBiasH.GetEnergy(rndm)); <<
334 <<
335 } <<
336 <<
337 G4double G4SPSRandomGenerator::GenRandY() <<
338 { <<
339 if (verbosityLevel >= 1) <<
340 { <<
341 G4cout << "In GenRandY" << G4endl; <<
342 } <<
343 <<
344 if (!YBias) // Y is not biased <<
345 { <<
346 G4double rndm = G4UniformRand(); <<
347 return (rndm); <<
348 } <<
349 // Y is biased <<
350 if ( !local_IPDFYBias.Get().val ) <<
351 { <<
352 local_IPDFYBias.Get().val = true; <<
353 G4AutoLock l(&mutex); <<
354 if (!IPDFYBias) <<
355 { <<
356 // IPDF has not been created, so creat <<
357 // <<
358 G4double bins[1024], vals[1024], sum; <<
359 std::size_t ii; <<
360 std::size_t maxbin = YBiasH.GetVectorL <<
361 bins[0] = YBiasH.GetLowEdgeEnergy(0); <<
362 vals[0] = YBiasH(0); <<
363 sum = vals[0]; <<
364 for (ii=1; ii<maxbin; ++ii) <<
365 { <<
366 bins[ii] = YBiasH.GetLowEdgeEnergy(i <<
367 vals[ii] = YBiasH(ii) + vals[ii - 1] <<
368 sum = sum + YBiasH(ii); <<
369 } <<
370 <<
371 for (ii=0; ii<maxbin; ++ii) <<
372 { <<
373 vals[ii] = vals[ii] / sum; <<
374 IPDFYBiasH.InsertValues(bins[ii], va <<
375 } <<
376 IPDFYBias = true; <<
377 } <<
378 } <<
379 <<
380 // IPDF has been created so carry on <<
381 <<
382 G4double rndm = G4UniformRand(); <<
383 std::size_t numberOfBin = IPDFYBiasH.GetVe <<
384 std::size_t biasn1 = 0; <<
385 std::size_t biasn2 = numberOfBin / 2; <<
386 std::size_t biasn3 = numberOfBin - 1; <<
387 while (biasn1 != biasn3 - 1) <<
388 { <<
389 if (rndm > IPDFYBiasH(biasn2)) <<
390 { biasn1 = biasn2; } <<
391 else <<
392 { biasn3 = biasn2; } <<
393 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
394 } <<
395 bweights_t& w = bweights.Get(); <<
396 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(bia <<
397 G4double xaxisl = IPDFYBiasH.GetLowEdgeEne <<
398 G4double xaxisu = IPDFYBiasH.GetLowEdgeEne <<
399 G4double NatProb = xaxisu - xaxisl; <<
400 w[1] = NatProb / w[1]; <<
401 if (verbosityLevel >= 1) <<
402 { <<
403 G4cout << "Y bin weight " << w[1] << " " <<
404 } <<
405 return (IPDFYBiasH.GetEnergy(rndm)); <<
406 <<
407 } <<
408 <<
409 G4double G4SPSRandomGenerator::GenRandZ() <<
410 { <<
411 if (verbosityLevel >= 1) <<
412 { <<
413 G4cout << "In GenRandZ" << G4endl; <<
414 } <<
415 <<
416 if (!ZBias) // Z is not biased <<
417 { <<
418 G4double rndm = G4UniformRand(); <<
419 return (rndm); <<
420 } <<
421 // Z is biased <<
422 if ( !local_IPDFZBias.Get().val ) <<
423 { <<
424 local_IPDFZBias.Get().val = true; <<
425 G4AutoLock l(&mutex); <<
426 if (!IPDFZBias) <<
427 { <<
428 // IPDF has not been created, so creat <<
429 // <<
430 G4double bins[1024], vals[1024], sum; <<
431 std::size_t ii; <<
432 std::size_t maxbin = ZBiasH.GetVectorL <<
433 bins[0] = ZBiasH.GetLowEdgeEnergy(0); <<
434 vals[0] = ZBiasH(0); <<
435 sum = vals[0]; <<
436 for (ii=1; ii<maxbin; ++ii) <<
437 { <<
438 bins[ii] = ZBiasH.GetLowEdgeEnergy(i <<
439 vals[ii] = ZBiasH(ii) + vals[ii - 1] <<
440 sum = sum + ZBiasH(ii); <<
441 } <<
442 <<
443 for (ii=0; ii<maxbin; ++ii) <<
444 { <<
445 vals[ii] = vals[ii] / sum; <<
446 IPDFZBiasH.InsertValues(bins[ii], va <<
447 } <<
448 IPDFZBias = true; <<
449 } <<
450 } <<
451 <<
452 // IPDF has been create so carry on <<
453 <<
454 G4double rndm = G4UniformRand(); <<
455 std::size_t numberOfBin = IPDFZBiasH.GetVe <<
456 std::size_t biasn1 = 0; <<
457 std::size_t biasn2 = numberOfBin / 2; <<
458 std::size_t biasn3 = numberOfBin - 1; <<
459 while (biasn1 != biasn3 - 1) <<
460 { <<
461 if (rndm > IPDFZBiasH(biasn2)) <<
462 { biasn1 = biasn2; } <<
463 else <<
464 { biasn3 = biasn2; } <<
465 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
466 } <<
467 bweights_t& w = bweights.Get(); <<
468 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(bia <<
469 G4double xaxisl = IPDFZBiasH.GetLowEdgeEne <<
470 G4double xaxisu = IPDFZBiasH.GetLowEdgeEne <<
471 G4double NatProb = xaxisu - xaxisl; <<
472 w[2] = NatProb / w[2]; <<
473 if (verbosityLevel >= 1) <<
474 { <<
475 G4cout << "Z bin weight " << w[2] << " " <<
476 } <<
477 return (IPDFZBiasH.GetEnergy(rndm)); <<
478 <<
479 } <<
480 <<
481 G4double G4SPSRandomGenerator::GenRandTheta() <<
482 { <<
483 if (verbosityLevel >= 1) <<
484 { <<
485 G4cout << "In GenRandTheta" << G4endl; <<
486 G4cout << "Verbosity " << verbosityLevel < <<
487 } <<
488 <<
489 if (!ThetaBias) // Theta is not biased <<
490 { <<
491 G4double rndm = G4UniformRand(); <<
492 return (rndm); <<
493 } <<
494 // Theta is biased <<
495 if ( !local_IPDFThetaBias.Get().val ) <<
496 { <<
497 local_IPDFThetaBias.Get().val = true; <<
498 G4AutoLock l(&mutex); <<
499 if (!IPDFThetaBias) <<
500 { <<
501 // IPDF has not been created, so creat <<
502 // <<
503 G4double bins[1024], vals[1024], sum; <<
504 std::size_t ii; <<
505 std::size_t maxbin = ThetaBiasH.GetVec <<
506 bins[0] = ThetaBiasH.GetLowEdgeEnergy( <<
507 vals[0] = ThetaBiasH(0); <<
508 sum = vals[0]; <<
509 for (ii=1; ii<maxbin; ++ii) <<
510 { <<
511 bins[ii] = ThetaBiasH.GetLowEdgeEner <<
512 vals[ii] = ThetaBiasH(ii) + vals[ii <<
513 sum = sum + ThetaBiasH(ii); <<
514 } <<
515 <<
516 for (ii=0; ii<maxbin; ++ii) <<
517 { <<
518 vals[ii] = vals[ii] / sum; <<
519 IPDFThetaBiasH.InsertValues(bins[ii] <<
520 } <<
521 IPDFThetaBias = true; <<
522 } <<
523 } <<
524 <<
525 // IPDF has been create so carry on <<
526 <<
527 G4double rndm = G4UniformRand(); <<
528 std::size_t numberOfBin = IPDFThetaBiasH.G <<
529 std::size_t biasn1 = 0; <<
530 std::size_t biasn2 = numberOfBin / 2; <<
531 std::size_t biasn3 = numberOfBin - 1; <<
532 while (biasn1 != biasn3 - 1) <<
533 { <<
534 if (rndm > IPDFThetaBiasH(biasn2)) <<
535 { biasn1 = biasn2; } <<
536 else <<
537 { biasn3 = biasn2; } <<
538 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
539 } <<
540 bweights_t& w = bweights.Get(); <<
541 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaB <<
542 G4double xaxisl = IPDFThetaBiasH.GetLowEdg <<
543 G4double xaxisu = IPDFThetaBiasH.GetLowEdg <<
544 G4double NatProb = xaxisu - xaxisl; <<
545 w[3] = NatProb / w[3]; <<
546 if (verbosityLevel >= 1) <<
547 { <<
548 G4cout << "Theta bin weight " << w[3] << <<
549 } <<
550 return (IPDFThetaBiasH.GetEnergy(rndm)); <<
551 <<
552 } <<
553 <<
554 G4double G4SPSRandomGenerator::GenRandPhi() <<
555 { <<
556 if (verbosityLevel >= 1) <<
557 { <<
558 G4cout << "In GenRandPhi" << G4endl; <<
559 } <<
560 <<
561 if (!PhiBias) // Phi is not biased <<
562 { <<
563 G4double rndm = G4UniformRand(); <<
564 return (rndm); <<
565 } <<
566 // Phi is biased <<
567 if ( !local_IPDFPhiBias.Get().val ) <<
568 { <<
569 local_IPDFPhiBias.Get().val = true; <<
570 G4AutoLock l(&mutex); <<
571 if (!IPDFPhiBias) <<
572 { <<
573 // IPDF has not been created, so creat <<
574 // <<
575 G4double bins[1024], vals[1024], sum; <<
576 std::size_t ii; <<
577 std::size_t maxbin = PhiBiasH.GetVecto <<
578 bins[0] = PhiBiasH.GetLowEdgeEnergy(0) <<
579 vals[0] = PhiBiasH(0); <<
580 sum = vals[0]; <<
581 for (ii=1; ii<maxbin; ++ii) <<
582 { <<
583 bins[ii] = PhiBiasH.GetLowEdgeEnergy <<
584 vals[ii] = PhiBiasH(ii) + vals[ii - <<
585 sum = sum + PhiBiasH(ii); <<
586 } <<
587 <<
588 for (ii=0; ii<maxbin; ++ii) <<
589 { <<
590 vals[ii] = vals[ii] / sum; <<
591 IPDFPhiBiasH.InsertValues(bins[ii], <<
592 } <<
593 IPDFPhiBias = true; <<
594 } <<
595 } <<
596 <<
597 // IPDF has been create so carry on <<
598 <<
599 G4double rndm = G4UniformRand(); <<
600 std::size_t numberOfBin = IPDFPhiBiasH.Get <<
601 std::size_t biasn1 = 0; <<
602 std::size_t biasn2 = numberOfBin / 2; <<
603 std::size_t biasn3 = numberOfBin - 1; <<
604 while (biasn1 != biasn3 - 1) <<
605 { <<
606 if (rndm > IPDFPhiBiasH(biasn2)) <<
607 { biasn1 = biasn2; } <<
608 else <<
609 { biasn3 = biasn2; } <<
610 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
611 } <<
612 bweights_t& w = bweights.Get(); <<
613 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH <<
614 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeE <<
615 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeE <<
616 G4double NatProb = xaxisu - xaxisl; <<
617 w[4] = NatProb / w[4]; <<
618 if (verbosityLevel >= 1) <<
619 { <<
620 G4cout << "Phi bin weight " << w[4] << " <<
621 } <<
622 return (IPDFPhiBiasH.GetEnergy(rndm)); <<
623 <<
624 } <<
625 <<
626 G4double G4SPSRandomGenerator::GenRandEnergy() <<
627 { <<
628 if (verbosityLevel >= 1) <<
629 { <<
630 G4cout << "In GenRandEnergy" << G4endl; <<
631 } <<
632 <<
633 if (!EnergyBias) // Energy is not biased <<
634 { <<
635 G4double rndm = G4UniformRand(); <<
636 return (rndm); <<
637 } <<
638 // Energy is biased <<
639 if ( !local_IPDFEnergyBias.Get().val ) <<
640 { <<
641 local_IPDFEnergyBias.Get().val = true; <<
642 G4AutoLock l(&mutex); <<
643 if (!IPDFEnergyBias) <<
644 { <<
645 // IPDF has not been created, so creat <<
646 // <<
647 G4double bins[1024], vals[1024], sum; <<
648 std::size_t ii; <<
649 std::size_t maxbin = EnergyBiasH.GetVe <<
650 bins[0] = EnergyBiasH.GetLowEdgeEnergy <<
651 vals[0] = EnergyBiasH(0); <<
652 sum = vals[0]; <<
653 for (ii=1; ii<maxbin; ++ii) <<
654 { <<
655 bins[ii] = EnergyBiasH.GetLowEdgeEne <<
656 vals[ii] = EnergyBiasH(ii) + vals[ii <<
657 sum = sum + EnergyBiasH(ii); <<
658 } <<
659 IPDFEnergyBiasH = ZeroPhysVector; <<
660 for (ii=0; ii<maxbin; ++ii) <<
661 { <<
662 vals[ii] = vals[ii] / sum; <<
663 IPDFEnergyBiasH.InsertValues(bins[ii <<
664 } <<
665 IPDFEnergyBias = true; <<
666 } <<
667 } <<
668 <<
669 // IPDF has been create so carry on <<
670 <<
671 G4double rndm = G4UniformRand(); <<
672 std::size_t numberOfBin = IPDFEnergyBiasH. <<
673 std::size_t biasn1 = 0; <<
674 std::size_t biasn2 = numberOfBin / 2; <<
675 std::size_t biasn3 = numberOfBin - 1; <<
676 while (biasn1 != biasn3 - 1) <<
677 { <<
678 if (rndm > IPDFEnergyBiasH(biasn2)) <<
679 { biasn1 = biasn2; } <<
680 else <<
681 { biasn3 = biasn2; } <<
682 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
683 } <<
684 bweights_t& w = bweights.Get(); <<
685 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnerg <<
686 G4double xaxisl = IPDFEnergyBiasH.GetLowEd <<
687 G4double xaxisu = IPDFEnergyBiasH.GetLowEd <<
688 G4double NatProb = xaxisu - xaxisl; <<
689 w[5] = NatProb / w[5]; <<
690 if (verbosityLevel >= 1) <<
691 { <<
692 G4cout << "Energy bin weight " << w[5] < <<
693 } <<
694 return (IPDFEnergyBiasH.GetEnergy(rndm)); <<
695 <<
696 } <<
697 <<
698 G4double G4SPSRandomGenerator::GenRandPosTheta <<
699 { <<
700 if (verbosityLevel >= 1) <<
701 { <<
702 G4cout << "In GenRandPosTheta" << G4endl; <<
703 G4cout << "Verbosity " << verbosityLevel < <<
704 } <<
705 <<
706 if (!PosThetaBias) // Theta is not biased <<
707 { <<
708 G4double rndm = G4UniformRand(); <<
709 return (rndm); <<
710 } <<
711 // Theta is biased <<
712 if ( !local_IPDFPosThetaBias.Get().val ) <<
713 { <<
714 local_IPDFPosThetaBias.Get().val = true; <<
715 G4AutoLock l(&mutex); <<
716 if (!IPDFPosThetaBias) <<
717 { <<
718 // IPDF has not been created, so creat <<
719 // <<
720 G4double bins[1024], vals[1024], sum; <<
721 std::size_t ii; <<
722 std::size_t maxbin = PosThetaBiasH.Get <<
723 bins[0] = PosThetaBiasH.GetLowEdgeEner <<
724 vals[0] = PosThetaBiasH(0); <<
725 sum = vals[0]; <<
726 for (ii=1; ii<maxbin; ++ii) <<
727 { <<
728 bins[ii] = PosThetaBiasH.GetLowEdgeE <<
729 vals[ii] = PosThetaBiasH(ii) + vals[ <<
730 sum = sum + PosThetaBiasH(ii); <<
731 } <<
732 <<
733 for (ii=0; ii<maxbin; ++ii) <<
734 { <<
735 vals[ii] = vals[ii] / sum; <<
736 IPDFPosThetaBiasH.InsertValues(bins[ <<
737 } <<
738 IPDFPosThetaBias = true; <<
739 } <<
740 } <<
741 <<
742 // IPDF has been create so carry on <<
743 // <<
744 G4double rndm = G4UniformRand(); <<
745 std::size_t numberOfBin = IPDFPosThetaBias <<
746 std::size_t biasn1 = 0; <<
747 std::size_t biasn2 = numberOfBin / 2; <<
748 std::size_t biasn3 = numberOfBin - 1; <<
749 while (biasn1 != biasn3 - 1) <<
750 { <<
751 if (rndm > IPDFPosThetaBiasH(biasn2)) <<
752 { biasn1 = biasn2; } <<
753 else <<
754 { biasn3 = biasn2; } <<
755 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
756 } <<
757 bweights_t& w = bweights.Get(); <<
758 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPos <<
759 G4double xaxisl = IPDFPosThetaBiasH.GetLow <<
760 G4double xaxisu = IPDFPosThetaBiasH.GetLow <<
761 G4double NatProb = xaxisu - xaxisl; <<
762 w[6] = NatProb / w[6]; <<
763 if (verbosityLevel >= 1) <<
764 { <<
765 G4cout << "PosTheta bin weight " << w[6] <<
766 } <<
767 return (IPDFPosThetaBiasH.GetEnergy(rndm)) <<
768 <<
769 } <<
770 <<
771 G4double G4SPSRandomGenerator::GenRandPosPhi() <<
772 { <<
773 if (verbosityLevel >= 1) <<
774 { <<
775 G4cout << "In GenRandPosPhi" << G4endl; <<
776 } <<
777 <<
778 if (!PosPhiBias) // PosPhi is not biased <<
779 { <<
780 G4double rndm = G4UniformRand(); <<
781 return (rndm); <<
782 } <<
783 // PosPhi is biased <<
784 if (!local_IPDFPosPhiBias.Get().val ) <<
785 { <<
786 local_IPDFPosPhiBias.Get().val = true; <<
787 G4AutoLock l(&mutex); <<
788 if (!IPDFPosPhiBias) <<
789 { <<
790 // IPDF has not been created, so creat <<
791 // <<
792 G4double bins[1024], vals[1024], sum; <<
793 std::size_t ii; <<
794 std::size_t maxbin = PosPhiBiasH.GetVe <<
795 bins[0] = PosPhiBiasH.GetLowEdgeEnergy <<
796 vals[0] = PosPhiBiasH(0); <<
797 sum = vals[0]; <<
798 for (ii=1; ii<maxbin; ++ii) <<
799 { <<
800 bins[ii] = PosPhiBiasH.GetLowEdgeEne <<
801 vals[ii] = PosPhiBiasH(ii) + vals[ii <<
802 sum = sum + PosPhiBiasH(ii); <<
803 } <<
804 <<
805 for (ii=0; ii<maxbin; ++ii) <<
806 { <<
807 vals[ii] = vals[ii] / sum; <<
808 IPDFPosPhiBiasH.InsertValues(bins[ii <<
809 } <<
810 IPDFPosPhiBias = true; <<
811 } <<
812 } <<
813 <<
814 // IPDF has been create so carry on <<
815 <<
816 G4double rndm = G4UniformRand(); <<
817 std::size_t numberOfBin = IPDFPosPhiBiasH. <<
818 std::size_t biasn1 = 0; <<
819 std::size_t biasn2 = numberOfBin / 2; <<
820 std::size_t biasn3 = numberOfBin - 1; <<
821 while (biasn1 != biasn3 - 1) <<
822 { <<
823 if (rndm > IPDFPosPhiBiasH(biasn2)) <<
824 { biasn1 = biasn2; } <<
825 else <<
826 { biasn3 = biasn2; } <<
827 biasn2 = biasn1 + (biasn3 - biasn1 + 1) <<
828 } <<
829 bweights_t& w = bweights.Get(); <<
830 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPh <<
831 G4double xaxisl = IPDFPosPhiBiasH.GetLowEd <<
832 G4double xaxisu = IPDFPosPhiBiasH.GetLowEd <<
833 G4double NatProb = xaxisu - xaxisl; <<
834 w[7] = NatProb / w[7]; <<
835 if (verbosityLevel >= 1) <<
836 { <<
837 G4cout << "PosPhi bin weight " << w[7] < <<
838 } <<
839 return (IPDFPosPhiBiasH.GetEnergy(rndm)); <<
840 <<
841 } <<
842 660