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