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