00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00027
00028
00029
00030
00031
00032
00033
00034
00035
00037
00038
00039
00040
00041
00042
00043
00044
00046
00047 #include "G4PrimaryParticle.hh"
00048 #include "G4Event.hh"
00049 #include "Randomize.hh"
00050 #include <cmath>
00051 #include "G4TransportationManager.hh"
00052 #include "G4VPhysicalVolume.hh"
00053 #include "G4PhysicalVolumeStore.hh"
00054 #include "G4ParticleTable.hh"
00055 #include "G4ParticleDefinition.hh"
00056 #include "G4IonTable.hh"
00057 #include "G4Ions.hh"
00058 #include "G4TrackingManager.hh"
00059 #include "G4Track.hh"
00060
00061 #include "G4SPSRandomGenerator.hh"
00062
00063
00064
00065 G4SPSRandomGenerator::G4SPSRandomGenerator()
00066 {
00067
00068
00069
00070 XBias = false;
00071 IPDFXBias = false;
00072 YBias = false;
00073 IPDFYBias = false;
00074 ZBias = false;
00075 IPDFZBias = false;
00076 ThetaBias = false;
00077 IPDFThetaBias = false;
00078 PhiBias = false;
00079 IPDFPhiBias = false;
00080 EnergyBias = false;
00081 IPDFEnergyBias = false;
00082 PosThetaBias = false;
00083 IPDFPosThetaBias = false;
00084 PosPhiBias = false;
00085 IPDFPosPhiBias = false;
00086 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4]
00087 = bweights[5] = bweights[6] = bweights[7] = bweights[8] = 1.;
00088 verbosityLevel = 0;
00089 }
00090
00091 G4SPSRandomGenerator::~G4SPSRandomGenerator() {
00092 }
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) {
00103 G4double ehi, val;
00104 ehi = input.x();
00105 val = input.y();
00106 XBiasH.InsertValues(ehi, val);
00107 XBias = true;
00108 }
00109
00110 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) {
00111 G4double ehi, val;
00112 ehi = input.x();
00113 val = input.y();
00114 YBiasH.InsertValues(ehi, val);
00115 YBias = true;
00116 }
00117
00118 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) {
00119 G4double ehi, val;
00120 ehi = input.x();
00121 val = input.y();
00122 ZBiasH.InsertValues(ehi, val);
00123 ZBias = true;
00124 }
00125
00126 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) {
00127 G4double ehi, val;
00128 ehi = input.x();
00129 val = input.y();
00130 ThetaBiasH.InsertValues(ehi, val);
00131 ThetaBias = true;
00132 }
00133
00134 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) {
00135 G4double ehi, val;
00136 ehi = input.x();
00137 val = input.y();
00138 PhiBiasH.InsertValues(ehi, val);
00139 PhiBias = true;
00140 }
00141
00142 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) {
00143 G4double ehi, val;
00144 ehi = input.x();
00145 val = input.y();
00146 EnergyBiasH.InsertValues(ehi, val);
00147 EnergyBias = true;
00148 }
00149
00150 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) {
00151 G4double ehi, val;
00152 ehi = input.x();
00153 val = input.y();
00154 PosThetaBiasH.InsertValues(ehi, val);
00155 PosThetaBias = true;
00156 }
00157
00158 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) {
00159 G4double ehi, val;
00160 ehi = input.x();
00161 val = input.y();
00162 PosPhiBiasH.InsertValues(ehi, val);
00163 PosPhiBias = true;
00164 }
00165
00166 void G4SPSRandomGenerator::ReSetHist(G4String atype) {
00167 if (atype == "biasx") {
00168 XBias = false;
00169 IPDFXBias = false;
00170 XBiasH = IPDFXBiasH = ZeroPhysVector;
00171 } else if (atype == "biasy") {
00172 YBias = false;
00173 IPDFYBias = false;
00174 YBiasH = IPDFYBiasH = ZeroPhysVector;
00175 } else if (atype == "biasz") {
00176 ZBias = false;
00177 IPDFZBias = false;
00178 ZBiasH = IPDFZBiasH = ZeroPhysVector;
00179 } else if (atype == "biast") {
00180 ThetaBias = false;
00181 IPDFThetaBias = false;
00182 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector;
00183 } else if (atype == "biasp") {
00184 PhiBias = false;
00185 IPDFPhiBias = false;
00186 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector;
00187 } else if (atype == "biase") {
00188 EnergyBias = false;
00189 IPDFEnergyBias = false;
00190 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector;
00191 } else if (atype == "biaspt") {
00192 PosThetaBias = false;
00193 IPDFPosThetaBias = false;
00194 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector;
00195 } else if (atype == "biaspp") {
00196 PosPhiBias = false;
00197 IPDFPosPhiBias = false;
00198 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector;
00199 } else {
00200 G4cout << "Error, histtype not accepted " << G4endl;
00201 }
00202 }
00203
00204 G4double G4SPSRandomGenerator::GenRandX() {
00205 if (verbosityLevel >= 1)
00206 G4cout << "In GenRandX" << G4endl;
00207 if (XBias == false) {
00208
00209 G4double rndm = G4UniformRand();
00210 return (rndm);
00211 } else {
00212
00213 if (IPDFXBias == false) {
00214
00215 G4double bins[1024], vals[1024], sum;
00216 G4int ii;
00217 G4int maxbin = G4int(XBiasH.GetVectorLength());
00218 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0));
00219 vals[0] = XBiasH(size_t(0));
00220 sum = vals[0];
00221 for (ii = 1; ii < maxbin; ii++) {
00222 bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii));
00223 vals[ii] = XBiasH(size_t(ii)) + vals[ii - 1];
00224 sum = sum + XBiasH(size_t(ii));
00225 }
00226
00227 for (ii = 0; ii < maxbin; ii++) {
00228 vals[ii] = vals[ii] / sum;
00229 IPDFXBiasH.InsertValues(bins[ii], vals[ii]);
00230 }
00231
00232 IPDFXBias = true;
00233 }
00234
00235 G4double rndm = G4UniformRand();
00236
00237
00238
00239
00240
00241 size_t numberOfBin = IPDFXBiasH.GetVectorLength();
00242 G4int biasn1 = 0;
00243 G4int biasn2 = numberOfBin / 2;
00244 G4int biasn3 = numberOfBin - 1;
00245 while (biasn1 != biasn3 - 1) {
00246 if (rndm > IPDFXBiasH(biasn2))
00247 biasn1 = biasn2;
00248 else
00249 biasn3 = biasn2;
00250 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00251 }
00252
00253 bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1);
00254 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00255 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2));
00256 G4double NatProb = xaxisu - xaxisl;
00257
00258
00259 bweights[0] = NatProb / bweights[0];
00260 if (verbosityLevel >= 1)
00261 G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl;
00262 return (IPDFXBiasH.GetEnergy(rndm));
00263 }
00264 }
00265
00266 G4double G4SPSRandomGenerator::GenRandY() {
00267 if (verbosityLevel >= 1)
00268 G4cout << "In GenRandY" << G4endl;
00269 if (YBias == false) {
00270
00271 G4double rndm = G4UniformRand();
00272 return (rndm);
00273 } else {
00274
00275 if (IPDFYBias == false) {
00276
00277 G4double bins[1024], vals[1024], sum;
00278 G4int ii;
00279 G4int maxbin = G4int(YBiasH.GetVectorLength());
00280 bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0));
00281 vals[0] = YBiasH(size_t(0));
00282 sum = vals[0];
00283 for (ii = 1; ii < maxbin; ii++) {
00284 bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii));
00285 vals[ii] = YBiasH(size_t(ii)) + vals[ii - 1];
00286 sum = sum + YBiasH(size_t(ii));
00287 }
00288
00289 for (ii = 0; ii < maxbin; ii++) {
00290 vals[ii] = vals[ii] / sum;
00291 IPDFYBiasH.InsertValues(bins[ii], vals[ii]);
00292 }
00293
00294 IPDFYBias = true;
00295 }
00296
00297 G4double rndm = G4UniformRand();
00298 size_t numberOfBin = IPDFYBiasH.GetVectorLength();
00299 G4int biasn1 = 0;
00300 G4int biasn2 = numberOfBin / 2;
00301 G4int biasn3 = numberOfBin - 1;
00302 while (biasn1 != biasn3 - 1) {
00303 if (rndm > IPDFYBiasH(biasn2))
00304 biasn1 = biasn2;
00305 else
00306 biasn3 = biasn2;
00307 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00308 }
00309 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1);
00310 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00311 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2));
00312 G4double NatProb = xaxisu - xaxisl;
00313 bweights[1] = NatProb / bweights[1];
00314 if (verbosityLevel >= 1)
00315 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl;
00316 return (IPDFYBiasH.GetEnergy(rndm));
00317 }
00318 }
00319
00320 G4double G4SPSRandomGenerator::GenRandZ() {
00321 if (verbosityLevel >= 1)
00322 G4cout << "In GenRandZ" << G4endl;
00323 if (ZBias == false) {
00324
00325 G4double rndm = G4UniformRand();
00326 return (rndm);
00327 } else {
00328
00329 if (IPDFZBias == false) {
00330
00331 G4double bins[1024], vals[1024], sum;
00332 G4int ii;
00333 G4int maxbin = G4int(ZBiasH.GetVectorLength());
00334 bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0));
00335 vals[0] = ZBiasH(size_t(0));
00336 sum = vals[0];
00337 for (ii = 1; ii < maxbin; ii++) {
00338 bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii));
00339 vals[ii] = ZBiasH(size_t(ii)) + vals[ii - 1];
00340 sum = sum + ZBiasH(size_t(ii));
00341 }
00342
00343 for (ii = 0; ii < maxbin; ii++) {
00344 vals[ii] = vals[ii] / sum;
00345 IPDFZBiasH.InsertValues(bins[ii], vals[ii]);
00346 }
00347
00348 IPDFZBias = true;
00349 }
00350
00351 G4double rndm = G4UniformRand();
00352
00353 size_t numberOfBin = IPDFZBiasH.GetVectorLength();
00354 G4int biasn1 = 0;
00355 G4int biasn2 = numberOfBin / 2;
00356 G4int biasn3 = numberOfBin - 1;
00357 while (biasn1 != biasn3 - 1) {
00358 if (rndm > IPDFZBiasH(biasn2))
00359 biasn1 = biasn2;
00360 else
00361 biasn3 = biasn2;
00362 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00363 }
00364 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1);
00365 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00366 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2));
00367 G4double NatProb = xaxisu - xaxisl;
00368 bweights[2] = NatProb / bweights[2];
00369 if (verbosityLevel >= 1)
00370 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl;
00371 return (IPDFZBiasH.GetEnergy(rndm));
00372 }
00373 }
00374
00375 G4double G4SPSRandomGenerator::GenRandTheta() {
00376 if (verbosityLevel >= 1) {
00377 G4cout << "In GenRandTheta" << G4endl;
00378 G4cout << "Verbosity " << verbosityLevel << G4endl;
00379 }
00380 if (ThetaBias == false) {
00381
00382 G4double rndm = G4UniformRand();
00383 return (rndm);
00384 } else {
00385
00386 if (IPDFThetaBias == false) {
00387
00388 G4double bins[1024], vals[1024], sum;
00389 G4int ii;
00390 G4int maxbin = G4int(ThetaBiasH.GetVectorLength());
00391 bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0));
00392 vals[0] = ThetaBiasH(size_t(0));
00393 sum = vals[0];
00394 for (ii = 1; ii < maxbin; ii++) {
00395 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii));
00396 vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii - 1];
00397 sum = sum + ThetaBiasH(size_t(ii));
00398 }
00399
00400 for (ii = 0; ii < maxbin; ii++) {
00401 vals[ii] = vals[ii] / sum;
00402 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]);
00403 }
00404
00405 IPDFThetaBias = true;
00406 }
00407
00408 G4double rndm = G4UniformRand();
00409
00410 size_t numberOfBin = IPDFThetaBiasH.GetVectorLength();
00411 G4int biasn1 = 0;
00412 G4int biasn2 = numberOfBin / 2;
00413 G4int biasn3 = numberOfBin - 1;
00414 while (biasn1 != biasn3 - 1) {
00415 if (rndm > IPDFThetaBiasH(biasn2))
00416 biasn1 = biasn2;
00417 else
00418 biasn3 = biasn2;
00419 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00420 }
00421 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1);
00422 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00423 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
00424 G4double NatProb = xaxisu - xaxisl;
00425 bweights[3] = NatProb / bweights[3];
00426 if (verbosityLevel >= 1)
00427 G4cout << "Theta bin weight " << bweights[3] << " " << rndm
00428 << G4endl;
00429 return (IPDFThetaBiasH.GetEnergy(rndm));
00430 }
00431 }
00432
00433 G4double G4SPSRandomGenerator::GenRandPhi() {
00434 if (verbosityLevel >= 1)
00435 G4cout << "In GenRandPhi" << G4endl;
00436 if (PhiBias == false) {
00437
00438 G4double rndm = G4UniformRand();
00439 return (rndm);
00440 } else {
00441
00442 if (IPDFPhiBias == false) {
00443
00444 G4double bins[1024], vals[1024], sum;
00445 G4int ii;
00446 G4int maxbin = G4int(PhiBiasH.GetVectorLength());
00447 bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0));
00448 vals[0] = PhiBiasH(size_t(0));
00449 sum = vals[0];
00450 for (ii = 1; ii < maxbin; ii++) {
00451 bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii));
00452 vals[ii] = PhiBiasH(size_t(ii)) + vals[ii - 1];
00453 sum = sum + PhiBiasH(size_t(ii));
00454 }
00455
00456 for (ii = 0; ii < maxbin; ii++) {
00457 vals[ii] = vals[ii] / sum;
00458 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]);
00459 }
00460
00461 IPDFPhiBias = true;
00462 }
00463
00464 G4double rndm = G4UniformRand();
00465
00466 size_t numberOfBin = IPDFPhiBiasH.GetVectorLength();
00467 G4int biasn1 = 0;
00468 G4int biasn2 = numberOfBin / 2;
00469 G4int biasn3 = numberOfBin - 1;
00470 while (biasn1 != biasn3 - 1) {
00471 if (rndm > IPDFPhiBiasH(biasn2))
00472 biasn1 = biasn2;
00473 else
00474 biasn3 = biasn2;
00475 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00476 }
00477 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1);
00478 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00479 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
00480 G4double NatProb = xaxisu - xaxisl;
00481 bweights[4] = NatProb / bweights[4];
00482 if (verbosityLevel >= 1)
00483 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl;
00484 return (IPDFPhiBiasH.GetEnergy(rndm));
00485 }
00486 }
00487
00488 G4double G4SPSRandomGenerator::GenRandEnergy() {
00489 if (verbosityLevel >= 1)
00490 G4cout << "In GenRandEnergy" << G4endl;
00491 if (EnergyBias == false) {
00492
00493 G4double rndm = G4UniformRand();
00494 return (rndm);
00495 } else {
00496
00497 if (IPDFEnergyBias == false) {
00498
00499 G4double bins[1024], vals[1024], sum;
00500 G4int ii;
00501 G4int maxbin = G4int(EnergyBiasH.GetVectorLength());
00502 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0));
00503 vals[0] = EnergyBiasH(size_t(0));
00504 sum = vals[0];
00505 for (ii = 1; ii < maxbin; ii++) {
00506 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii));
00507 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii - 1];
00508 sum = sum + EnergyBiasH(size_t(ii));
00509 }
00510 IPDFEnergyBiasH = ZeroPhysVector;
00511 for (ii = 0; ii < maxbin; ii++) {
00512 vals[ii] = vals[ii] / sum;
00513 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]);
00514 }
00515
00516 IPDFEnergyBias = true;
00517 }
00518
00519 G4double rndm = G4UniformRand();
00520
00521 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength();
00522 G4int biasn1 = 0;
00523 G4int biasn2 = numberOfBin / 2;
00524 G4int biasn3 = numberOfBin - 1;
00525 while (biasn1 != biasn3 - 1) {
00526 if (rndm > IPDFEnergyBiasH(biasn2))
00527 biasn1 = biasn2;
00528 else
00529 biasn3 = biasn2;
00530 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00531 }
00532 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1);
00533 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00534 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2));
00535 G4double NatProb = xaxisu - xaxisl;
00536 bweights[5] = NatProb / bweights[5];
00537 if (verbosityLevel >= 1)
00538 G4cout << "Energy bin weight " << bweights[5] << " " << rndm
00539 << G4endl;
00540 return (IPDFEnergyBiasH.GetEnergy(rndm));
00541 }
00542 }
00543
00544 G4double G4SPSRandomGenerator::GenRandPosTheta() {
00545 if (verbosityLevel >= 1) {
00546 G4cout << "In GenRandPosTheta" << G4endl;
00547 G4cout << "Verbosity " << verbosityLevel << G4endl;
00548 }
00549 if (PosThetaBias == false) {
00550
00551 G4double rndm = G4UniformRand();
00552 return (rndm);
00553 } else {
00554
00555 if (IPDFPosThetaBias == false) {
00556
00557 G4double bins[1024], vals[1024], sum;
00558 G4int ii;
00559 G4int maxbin = G4int(PosThetaBiasH.GetVectorLength());
00560 bins[0] = PosThetaBiasH.GetLowEdgeEnergy(size_t(0));
00561 vals[0] = PosThetaBiasH(size_t(0));
00562 sum = vals[0];
00563 for (ii = 1; ii < maxbin; ii++) {
00564 bins[ii] = PosThetaBiasH.GetLowEdgeEnergy(size_t(ii));
00565 vals[ii] = PosThetaBiasH(size_t(ii)) + vals[ii - 1];
00566 sum = sum + PosThetaBiasH(size_t(ii));
00567 }
00568
00569 for (ii = 0; ii < maxbin; ii++) {
00570 vals[ii] = vals[ii] / sum;
00571 IPDFPosThetaBiasH.InsertValues(bins[ii], vals[ii]);
00572 }
00573
00574 IPDFPosThetaBias = true;
00575 }
00576
00577 G4double rndm = G4UniformRand();
00578
00579 size_t numberOfBin = IPDFPosThetaBiasH.GetVectorLength();
00580 G4int biasn1 = 0;
00581 G4int biasn2 = numberOfBin / 2;
00582 G4int biasn3 = numberOfBin - 1;
00583 while (biasn1 != biasn3 - 1) {
00584 if (rndm > IPDFPosThetaBiasH(biasn2))
00585 biasn1 = biasn2;
00586 else
00587 biasn3 = biasn2;
00588 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00589 }
00590 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1);
00591 G4double xaxisl =
00592 IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00593 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2));
00594 G4double NatProb = xaxisu - xaxisl;
00595 bweights[6] = NatProb / bweights[6];
00596 if (verbosityLevel >= 1)
00597 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm
00598 << G4endl;
00599 return (IPDFPosThetaBiasH.GetEnergy(rndm));
00600 }
00601 }
00602
00603 G4double G4SPSRandomGenerator::GenRandPosPhi() {
00604 if (verbosityLevel >= 1)
00605 G4cout << "In GenRandPosPhi" << G4endl;
00606 if (PosPhiBias == false) {
00607
00608 G4double rndm = G4UniformRand();
00609 return (rndm);
00610 } else {
00611
00612 if (IPDFPosPhiBias == false) {
00613
00614 G4double bins[1024], vals[1024], sum;
00615 G4int ii;
00616 G4int maxbin = G4int(PosPhiBiasH.GetVectorLength());
00617 bins[0] = PosPhiBiasH.GetLowEdgeEnergy(size_t(0));
00618 vals[0] = PosPhiBiasH(size_t(0));
00619 sum = vals[0];
00620 for (ii = 1; ii < maxbin; ii++) {
00621 bins[ii] = PosPhiBiasH.GetLowEdgeEnergy(size_t(ii));
00622 vals[ii] = PosPhiBiasH(size_t(ii)) + vals[ii - 1];
00623 sum = sum + PosPhiBiasH(size_t(ii));
00624 }
00625
00626 for (ii = 0; ii < maxbin; ii++) {
00627 vals[ii] = vals[ii] / sum;
00628 IPDFPosPhiBiasH.InsertValues(bins[ii], vals[ii]);
00629 }
00630
00631 IPDFPosPhiBias = true;
00632 }
00633
00634 G4double rndm = G4UniformRand();
00635
00636 size_t numberOfBin = IPDFPosPhiBiasH.GetVectorLength();
00637 G4int biasn1 = 0;
00638 G4int biasn2 = numberOfBin / 2;
00639 G4int biasn3 = numberOfBin - 1;
00640 while (biasn1 != biasn3 - 1) {
00641 if (rndm > IPDFPosPhiBiasH(biasn2))
00642 biasn1 = biasn2;
00643 else
00644 biasn3 = biasn2;
00645 biasn2 = biasn1 + (biasn3 - biasn1 + 1) / 2;
00646 }
00647 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1);
00648 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2 - 1));
00649 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2));
00650 G4double NatProb = xaxisu - xaxisl;
00651 bweights[7] = NatProb / bweights[7];
00652 if (verbosityLevel >= 1)
00653 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm
00654 << G4endl;
00655 return (IPDFPosPhiBiasH.GetEnergy(rndm));
00656 }
00657 }
00658