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() 62 { 66 { 63 // Initialise all variables 67 // Initialise all variables 64 68 65 // Bias variables 69 // Bias variables 66 // << 67 XBias = false; 70 XBias = false; 68 IPDFXBias = false; 71 IPDFXBias = false; 69 YBias = false; 72 YBias = false; 70 IPDFYBias = false; 73 IPDFYBias = false; 71 ZBias = false; 74 ZBias = false; 72 IPDFZBias = false; 75 IPDFZBias = false; 73 ThetaBias = false; 76 ThetaBias = false; 74 IPDFThetaBias = false; 77 IPDFThetaBias = false; 75 PhiBias = false; 78 PhiBias = false; 76 IPDFPhiBias = false; 79 IPDFPhiBias = false; 77 EnergyBias = false; 80 EnergyBias = false; 78 IPDFEnergyBias = false; 81 IPDFEnergyBias = false; 79 PosThetaBias = false; 82 PosThetaBias = false; 80 IPDFPosThetaBias = false; 83 IPDFPosThetaBias = false; 81 PosPhiBias = false; 84 PosPhiBias = false; 82 IPDFPosPhiBias = false; 85 IPDFPosPhiBias = false; 83 << 86 bweights[0] = bweights[1] = bweights[2] = bweights[3] = bweights[4] = bweights[5] = bweights[6] = bweights[7]= bweights[8]= 1. ; 84 verbosityLevel = 0; << 87 verbosityLevel = 0 ; 85 G4MUTEXINIT(mutex); << 86 } 88 } 87 89 88 G4SPSRandomGenerator::~G4SPSRandomGenerator() 90 G4SPSRandomGenerator::~G4SPSRandomGenerator() 89 { << 91 {} 90 G4MUTEXDESTROY(mutex); << 92 91 } << 93 //G4SPSRandomGenerator* G4SPSRandomGenerator::getInstance () >> 94 //{ >> 95 // if (instance == 0) instance = new G4SPSRandomGenerator(); >> 96 // return instance; >> 97 //} 92 98 93 // Biasing methods 99 // Biasing methods 94 100 95 void G4SPSRandomGenerator::SetXBias(const G4Th << 101 void G4SPSRandomGenerator::SetXBias(G4ThreeVector input) 96 { << 102 { 97 G4AutoLock l(&mutex); << 98 G4double ehi, val; 103 G4double ehi, val; 99 ehi = input.x(); 104 ehi = input.x(); 100 val = input.y(); 105 val = input.y(); 101 XBiasH.InsertValues(ehi, val); 106 XBiasH.InsertValues(ehi, val); 102 XBias = true; 107 XBias = true; 103 } 108 } 104 109 105 void G4SPSRandomGenerator::SetYBias(const G4Th << 110 void G4SPSRandomGenerator::SetYBias(G4ThreeVector input) 106 { 111 { 107 G4AutoLock l(&mutex); << 108 G4double ehi, val; 112 G4double ehi, val; 109 ehi = input.x(); 113 ehi = input.x(); 110 val = input.y(); 114 val = input.y(); 111 YBiasH.InsertValues(ehi, val); 115 YBiasH.InsertValues(ehi, val); 112 YBias = true; 116 YBias = true; 113 } 117 } 114 118 115 void G4SPSRandomGenerator::SetZBias(const G4Th << 119 void G4SPSRandomGenerator::SetZBias(G4ThreeVector input) 116 { 120 { 117 G4AutoLock l(&mutex); << 118 G4double ehi, val; 121 G4double ehi, val; 119 ehi = input.x(); 122 ehi = input.x(); 120 val = input.y(); 123 val = input.y(); 121 ZBiasH.InsertValues(ehi, val); 124 ZBiasH.InsertValues(ehi, val); 122 ZBias = true; 125 ZBias = true; 123 } 126 } 124 127 125 void G4SPSRandomGenerator::SetThetaBias(const << 128 void G4SPSRandomGenerator::SetThetaBias(G4ThreeVector input) 126 { 129 { 127 G4AutoLock l(&mutex); << 128 G4double ehi, val; 130 G4double ehi, val; 129 ehi = input.x(); 131 ehi = input.x(); 130 val = input.y(); 132 val = input.y(); 131 ThetaBiasH.InsertValues(ehi, val); 133 ThetaBiasH.InsertValues(ehi, val); 132 ThetaBias = true; 134 ThetaBias = true; 133 } 135 } 134 136 135 void G4SPSRandomGenerator::SetPhiBias(const G4 << 137 void G4SPSRandomGenerator::SetPhiBias(G4ThreeVector input) 136 { 138 { 137 G4AutoLock l(&mutex); << 138 G4double ehi, val; 139 G4double ehi, val; 139 ehi = input.x(); 140 ehi = input.x(); 140 val = input.y(); 141 val = input.y(); 141 PhiBiasH.InsertValues(ehi, val); 142 PhiBiasH.InsertValues(ehi, val); 142 PhiBias = true; 143 PhiBias = true; 143 } 144 } 144 145 145 void G4SPSRandomGenerator::SetEnergyBias(const << 146 void G4SPSRandomGenerator::SetEnergyBias(G4ThreeVector input) 146 { 147 { 147 G4AutoLock l(&mutex); << 148 G4double ehi, val; 148 G4double ehi, val; 149 ehi = input.x(); 149 ehi = input.x(); 150 val = input.y(); 150 val = input.y(); 151 EnergyBiasH.InsertValues(ehi, val); 151 EnergyBiasH.InsertValues(ehi, val); 152 EnergyBias = true; 152 EnergyBias = true; 153 } 153 } 154 154 155 void G4SPSRandomGenerator::SetPosThetaBias(con << 155 void G4SPSRandomGenerator::SetPosThetaBias(G4ThreeVector input) 156 { 156 { 157 G4AutoLock l(&mutex); << 158 G4double ehi, val; 157 G4double ehi, val; 159 ehi = input.x(); 158 ehi = input.x(); 160 val = input.y(); 159 val = input.y(); 161 PosThetaBiasH.InsertValues(ehi, val); 160 PosThetaBiasH.InsertValues(ehi, val); 162 PosThetaBias = true; 161 PosThetaBias = true; 163 } 162 } 164 163 165 void G4SPSRandomGenerator::SetPosPhiBias(const << 164 void G4SPSRandomGenerator::SetPosPhiBias(G4ThreeVector input) 166 { 165 { 167 G4AutoLock l(&mutex); << 168 G4double ehi, val; 166 G4double ehi, val; 169 ehi = input.x(); 167 ehi = input.x(); 170 val = input.y(); 168 val = input.y(); 171 PosPhiBiasH.InsertValues(ehi, val); 169 PosPhiBiasH.InsertValues(ehi, val); 172 PosPhiBias = true; 170 PosPhiBias = true; 173 } 171 } 174 172 175 void G4SPSRandomGenerator::SetIntensityWeight( << 173 void G4SPSRandomGenerator::ReSetHist(G4String atype) 176 { 174 { 177 bweights.Get()[8] = weight; << 175 if ( atype == "biasx") { 178 } << 176 XBias = false ; 179 << 177 IPDFXBias = false; 180 G4double G4SPSRandomGenerator::GetBiasWeight() << 178 XBiasH = IPDFXBiasH = ZeroPhysVector ;} 181 { << 179 else if ( atype == "biasy") { 182 bweights_t& w = bweights.Get(); << 180 YBias = false ; 183 return w[0] * w[1] * w[2] * w[3] * w[4] * w[ << 181 IPDFYBias = false; 184 } << 182 YBiasH = IPDFYBiasH = ZeroPhysVector ;} 185 << 183 else if ( atype == "biasz") { 186 void G4SPSRandomGenerator::SetVerbosity(G4int << 184 ZBias = false ; 187 { << 185 IPDFZBias = false; 188 G4AutoLock l(&mutex); << 186 ZBiasH = IPDFZBiasH = ZeroPhysVector ;} 189 verbosityLevel = a; << 187 else if ( atype == "biast") { 190 } << 188 ThetaBias = false ; 191 << 189 IPDFThetaBias = false; 192 namespace << 190 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector ;} 193 { << 191 else if ( atype == "biasp") { 194 G4PhysicsFreeVector ZeroPhysVector; // for r << 192 PhiBias = false ; 195 } << 193 IPDFPhiBias = false; 196 << 194 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector ;} 197 void G4SPSRandomGenerator::ReSetHist(const G4S << 195 else if ( atype == "biase") { 198 { << 196 EnergyBias = false ; 199 G4AutoLock l(&mutex); << 197 IPDFEnergyBias = false; 200 if (atype == "biasx") { << 198 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector ;} 201 XBias = false; << 199 else if ( atype == "biaspt") { 202 IPDFXBias = false; << 200 PosThetaBias = false ; 203 local_IPDFXBias.Get().val = fa << 201 IPDFPosThetaBias = false; 204 XBiasH = IPDFXBiasH = ZeroPhys << 202 PosThetaBiasH = IPDFPosThetaBiasH = ZeroPhysVector ;} 205 } else if (atype == "biasy") { << 203 else if ( atype == "biaspp") { 206 YBias = false; << 204 PosPhiBias = false ; 207 IPDFYBias = false; << 205 IPDFPosPhiBias = false; 208 local_IPDFYBias.Get().val = fa << 206 PosPhiBiasH = IPDFPosPhiBiasH = ZeroPhysVector ;} 209 YBiasH = IPDFYBiasH = ZeroPhys << 207 else { 210 } else if (atype == "biasz") { << 208 G4cout << "Error, histtype not accepted " << G4endl; 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 } 209 } 243 } 210 } 244 211 245 G4double G4SPSRandomGenerator::GenRandX() 212 G4double G4SPSRandomGenerator::GenRandX() 246 { 213 { 247 if (verbosityLevel >= 1) << 214 if(verbosityLevel >= 1) 248 G4cout << "In GenRandX" << G4endl; 215 G4cout << "In GenRandX" << G4endl; 249 if (!XBias) << 216 if(XBias == false) 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 { 217 { 275 // IPDF has not been created, so create << 218 // X is not biased 276 // << 219 G4double rndm = G4UniformRand(); 277 G4double bins[1024], vals[1024], sum; << 220 return(rndm); 278 std::size_t ii; << 221 } 279 std::size_t maxbin = XBiasH.GetVectorLen << 222 else 280 bins[0] = XBiasH.GetLowEdgeEnergy(0); << 223 { 281 vals[0] = XBiasH(0); << 224 // X is biased 282 sum = vals[0]; << 225 if(IPDFXBias == false) 283 for (ii=1; ii<maxbin; ++ii) << 226 { 284 { << 227 // IPDF has not been created, so create it 285 bins[ii] = XBiasH.GetLowEdgeEnergy(ii) << 228 G4double bins[1024],vals[1024], sum; 286 vals[ii] = XBiasH(ii) + vals[ii - 1]; << 229 G4int ii; 287 sum = sum + XBiasH(ii); << 230 G4int maxbin = G4int(XBiasH.GetVectorLength()); 288 } << 231 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0)); 289 << 232 vals[0] = XBiasH(size_t(0)); 290 for (ii=0; ii<maxbin; ++ii) << 233 sum = vals[0]; 291 { << 234 for(ii=1;ii<maxbin;ii++) 292 vals[ii] = vals[ii] / sum; << 235 { 293 IPDFXBiasH.InsertValues(bins[ii], vals << 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; 294 } 266 } 295 IPDFXBias = true; << 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)); 296 } 278 } 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 } 279 } 336 280 337 G4double G4SPSRandomGenerator::GenRandY() 281 G4double G4SPSRandomGenerator::GenRandY() 338 { 282 { 339 if (verbosityLevel >= 1) << 283 if(verbosityLevel >= 1) 340 { << 341 G4cout << "In GenRandY" << G4endl; 284 G4cout << "In GenRandY" << G4endl; 342 } << 285 if(YBias == false) 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 { 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) { 389 if (rndm > IPDFYBiasH(biasn2)) 325 if (rndm > IPDFYBiasH(biasn2)) 390 { biasn1 = biasn2; } << 326 biasn1 = biasn2; 391 else 327 else 392 { biasn3 = biasn2; } << 328 biasn3 = biasn2; 393 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 329 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 394 } << 330 } 395 bweights_t& w = bweights.Get(); << 331 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1); 396 w[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(bia << 332 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 397 G4double xaxisl = IPDFYBiasH.GetLowEdgeEne << 333 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2)); 398 G4double xaxisu = IPDFYBiasH.GetLowEdgeEne << 334 G4double NatProb = xaxisu - xaxisl; 399 G4double NatProb = xaxisu - xaxisl; << 335 bweights[1] = NatProb/bweights[1]; 400 w[1] = NatProb / w[1]; << 336 if(verbosityLevel >= 1) 401 if (verbosityLevel >= 1) << 337 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl; 402 { << 338 return(IPDFYBiasH.GetEnergy(rndm)); 403 G4cout << "Y bin weight " << w[1] << " " << 404 } 339 } 405 return (IPDFYBiasH.GetEnergy(rndm)); << 406 << 407 } 340 } 408 341 409 G4double G4SPSRandomGenerator::GenRandZ() 342 G4double G4SPSRandomGenerator::GenRandZ() 410 { 343 { 411 if (verbosityLevel >= 1) << 344 if(verbosityLevel >= 1) 412 { << 413 G4cout << "In GenRandZ" << G4endl; 345 G4cout << "In GenRandZ" << G4endl; 414 } << 346 if(ZBias == false) 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 { 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) { 461 if (rndm > IPDFZBiasH(biasn2)) 387 if (rndm > IPDFZBiasH(biasn2)) 462 { biasn1 = biasn2; } << 388 biasn1 = biasn2; 463 else 389 else 464 { biasn3 = biasn2; } << 390 biasn3 = biasn2; 465 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 391 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 466 } << 392 } 467 bweights_t& w = bweights.Get(); << 393 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1); 468 w[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(bia << 394 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 469 G4double xaxisl = IPDFZBiasH.GetLowEdgeEne << 395 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2)); 470 G4double xaxisu = IPDFZBiasH.GetLowEdgeEne << 396 G4double NatProb = xaxisu - xaxisl; 471 G4double NatProb = xaxisu - xaxisl; << 397 bweights[2] = NatProb/bweights[2]; 472 w[2] = NatProb / w[2]; << 398 if(verbosityLevel >= 1) 473 if (verbosityLevel >= 1) << 399 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl; 474 { << 400 return(IPDFZBiasH.GetEnergy(rndm)); 475 G4cout << "Z bin weight " << w[2] << " " << 476 } 401 } 477 return (IPDFZBiasH.GetEnergy(rndm)); << 478 << 479 } 402 } 480 403 481 G4double G4SPSRandomGenerator::GenRandTheta() 404 G4double G4SPSRandomGenerator::GenRandTheta() 482 { 405 { 483 if (verbosityLevel >= 1) << 406 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 { 407 { 497 local_IPDFThetaBias.Get().val = true; << 408 G4cout << "In GenRandTheta" << G4endl; 498 G4AutoLock l(&mutex); << 409 G4cout << "Verbosity " << verbosityLevel << G4endl; 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 } 410 } 524 << 411 if(ThetaBias == false) 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 { 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) { 534 if (rndm > IPDFThetaBiasH(biasn2)) 452 if (rndm > IPDFThetaBiasH(biasn2)) 535 { biasn1 = biasn2; } << 453 biasn1 = biasn2; 536 else 454 else 537 { biasn3 = biasn2; } << 455 biasn3 = biasn2; 538 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 456 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 539 } << 457 } 540 bweights_t& w = bweights.Get(); << 458 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1); 541 w[3] = IPDFThetaBiasH(biasn2) - IPDFThetaB << 459 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 542 G4double xaxisl = IPDFThetaBiasH.GetLowEdg << 460 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2)); 543 G4double xaxisu = IPDFThetaBiasH.GetLowEdg << 461 G4double NatProb = xaxisu - xaxisl; 544 G4double NatProb = xaxisu - xaxisl; << 462 bweights[3] = NatProb/bweights[3]; 545 w[3] = NatProb / w[3]; << 463 if(verbosityLevel >= 1) 546 if (verbosityLevel >= 1) << 464 G4cout << "Theta bin weight " << bweights[3] << " " << rndm << G4endl; 547 { << 465 return(IPDFThetaBiasH.GetEnergy(rndm)); 548 G4cout << "Theta bin weight " << w[3] << << 549 } 466 } 550 return (IPDFThetaBiasH.GetEnergy(rndm)); << 551 << 552 } 467 } 553 468 554 G4double G4SPSRandomGenerator::GenRandPhi() 469 G4double G4SPSRandomGenerator::GenRandPhi() 555 { 470 { 556 if (verbosityLevel >= 1) << 471 if(verbosityLevel >= 1) 557 { << 558 G4cout << "In GenRandPhi" << G4endl; 472 G4cout << "In GenRandPhi" << G4endl; 559 } << 473 if(PhiBias == false) 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 { 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) { 606 if (rndm > IPDFPhiBiasH(biasn2)) 514 if (rndm > IPDFPhiBiasH(biasn2)) 607 { biasn1 = biasn2; } << 515 biasn1 = biasn2; 608 else 516 else 609 { biasn3 = biasn2; } << 517 biasn3 = biasn2; 610 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 518 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 611 } << 519 } 612 bweights_t& w = bweights.Get(); << 520 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1); 613 w[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH << 521 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 614 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeE << 522 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2)); 615 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeE << 523 G4double NatProb = xaxisu - xaxisl; 616 G4double NatProb = xaxisu - xaxisl; << 524 bweights[4] = NatProb/bweights[4]; 617 w[4] = NatProb / w[4]; << 525 if(verbosityLevel >= 1) 618 if (verbosityLevel >= 1) << 526 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl; 619 { << 527 return(IPDFPhiBiasH.GetEnergy(rndm)); 620 G4cout << "Phi bin weight " << w[4] << " << 621 } 528 } 622 return (IPDFPhiBiasH.GetEnergy(rndm)); << 623 << 624 } 529 } 625 530 626 G4double G4SPSRandomGenerator::GenRandEnergy() 531 G4double G4SPSRandomGenerator::GenRandEnergy() 627 { 532 { 628 if (verbosityLevel >= 1) << 533 if(verbosityLevel >= 1) 629 { << 630 G4cout << "In GenRandEnergy" << G4endl; 534 G4cout << "In GenRandEnergy" << G4endl; 631 } << 535 if(EnergyBias == false) 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 { 536 { 641 local_IPDFEnergyBias.Get().val = true; << 537 // Energy is not biased 642 G4AutoLock l(&mutex); << 538 G4double rndm = G4UniformRand(); 643 if (!IPDFEnergyBias) << 539 return(rndm); 644 { << 540 } 645 // IPDF has not been created, so creat << 541 else { 646 // << 542 // ENERGY is biased 647 G4double bins[1024], vals[1024], sum; << 543 if(IPDFEnergyBias == false) { 648 std::size_t ii; << 544 // IPDF has not been created, so create it 649 std::size_t maxbin = EnergyBiasH.GetVe << 545 G4double bins[1024],vals[1024], sum; 650 bins[0] = EnergyBiasH.GetLowEdgeEnergy << 546 G4int ii; 651 vals[0] = EnergyBiasH(0); << 547 G4int maxbin = G4int(EnergyBiasH.GetVectorLength()); 652 sum = vals[0]; << 548 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0)); 653 for (ii=1; ii<maxbin; ++ii) << 549 vals[0] = EnergyBiasH(size_t(0)); 654 { << 550 sum = vals[0]; 655 bins[ii] = EnergyBiasH.GetLowEdgeEne << 551 for(ii=1;ii<maxbin;ii++) { 656 vals[ii] = EnergyBiasH(ii) + vals[ii << 552 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii)); 657 sum = sum + EnergyBiasH(ii); << 553 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii-1]; 658 } << 554 sum = sum + EnergyBiasH(size_t(ii)); 659 IPDFEnergyBiasH = ZeroPhysVector; << 555 } 660 for (ii=0; ii<maxbin; ++ii) << 556 661 { << 557 for(ii=0;ii<maxbin;ii++) { 662 vals[ii] = vals[ii] / sum; << 558 vals[ii] = vals[ii]/sum; 663 IPDFEnergyBiasH.InsertValues(bins[ii << 559 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]); 664 } << 665 IPDFEnergyBias = true; << 666 } 560 } >> 561 // Make IPDFEnergyBias = true >> 562 IPDFEnergyBias = true; 667 } 563 } 668 << 669 // IPDF has been create so carry on 564 // IPDF has been create so carry on 670 << 671 G4double rndm = G4UniformRand(); 565 G4double rndm = G4UniformRand(); 672 std::size_t numberOfBin = IPDFEnergyBiasH. << 566 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm); 673 std::size_t biasn1 = 0; << 567 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength(); 674 std::size_t biasn2 = numberOfBin / 2; << 568 G4int biasn1 = 0; 675 std::size_t biasn3 = numberOfBin - 1; << 569 G4int biasn2 = numberOfBin/2; 676 while (biasn1 != biasn3 - 1) << 570 G4int biasn3 = numberOfBin - 1; 677 { << 571 while (biasn1 != biasn3 - 1) { 678 if (rndm > IPDFEnergyBiasH(biasn2)) 572 if (rndm > IPDFEnergyBiasH(biasn2)) 679 { biasn1 = biasn2; } << 573 biasn1 = biasn2; 680 else 574 else 681 { biasn3 = biasn2; } << 575 biasn3 = biasn2; 682 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 576 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 683 } 577 } 684 bweights_t& w = bweights.Get(); << 578 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1); 685 w[5] = IPDFEnergyBiasH(biasn2) - IPDFEnerg << 579 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 686 G4double xaxisl = IPDFEnergyBiasH.GetLowEd << 580 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2)); 687 G4double xaxisu = IPDFEnergyBiasH.GetLowEd << 688 G4double NatProb = xaxisu - xaxisl; 581 G4double NatProb = xaxisu - xaxisl; 689 w[5] = NatProb / w[5]; << 582 bweights[5] = NatProb/bweights[5]; 690 if (verbosityLevel >= 1) << 583 if(verbosityLevel >= 1) 691 { << 584 G4cout << "Energy bin weight " << bweights[5] << " " << rndm << G4endl; 692 G4cout << "Energy bin weight " << w[5] < << 585 return(IPDFEnergyBiasH.GetEnergy(rndm)); 693 } << 586 } 694 return (IPDFEnergyBiasH.GetEnergy(rndm)); << 695 << 696 } 587 } 697 588 >> 589 698 G4double G4SPSRandomGenerator::GenRandPosTheta 590 G4double G4SPSRandomGenerator::GenRandPosTheta() 699 { 591 { 700 if (verbosityLevel >= 1) << 592 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 { 593 { 714 local_IPDFPosThetaBias.Get().val = true; << 594 G4cout << "In GenRandPosTheta" << G4endl; 715 G4AutoLock l(&mutex); << 595 G4cout << "Verbosity " << verbosityLevel << G4endl; 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 } 596 } 741 << 597 if(PosThetaBias == false) 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 { 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) { 751 if (rndm > IPDFPosThetaBiasH(biasn2)) 638 if (rndm > IPDFPosThetaBiasH(biasn2)) 752 { biasn1 = biasn2; } << 639 biasn1 = biasn2; 753 else 640 else 754 { biasn3 = biasn2; } << 641 biasn3 = biasn2; 755 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 642 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 756 } << 643 } 757 bweights_t& w = bweights.Get(); << 644 bweights[6] = IPDFPosThetaBiasH(biasn2) - IPDFPosThetaBiasH(biasn2 - 1); 758 w[6] = IPDFPosThetaBiasH(biasn2) - IPDFPos << 645 G4double xaxisl = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 759 G4double xaxisl = IPDFPosThetaBiasH.GetLow << 646 G4double xaxisu = IPDFPosThetaBiasH.GetLowEdgeEnergy(size_t(biasn2)); 760 G4double xaxisu = IPDFPosThetaBiasH.GetLow << 647 G4double NatProb = xaxisu - xaxisl; 761 G4double NatProb = xaxisu - xaxisl; << 648 bweights[6] = NatProb/bweights[6]; 762 w[6] = NatProb / w[6]; << 649 if(verbosityLevel >= 1) 763 if (verbosityLevel >= 1) << 650 G4cout << "PosTheta bin weight " << bweights[6] << " " << rndm << G4endl; 764 { << 651 return(IPDFPosThetaBiasH.GetEnergy(rndm)); 765 G4cout << "PosTheta bin weight " << w[6] << 766 } 652 } 767 return (IPDFPosThetaBiasH.GetEnergy(rndm)) << 768 << 769 } 653 } 770 654 771 G4double G4SPSRandomGenerator::GenRandPosPhi() 655 G4double G4SPSRandomGenerator::GenRandPosPhi() 772 { 656 { 773 if (verbosityLevel >= 1) << 657 if(verbosityLevel >= 1) 774 { << 775 G4cout << "In GenRandPosPhi" << G4endl; 658 G4cout << "In GenRandPosPhi" << G4endl; 776 } << 659 if(PosPhiBias == false) 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 { 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) { 823 if (rndm > IPDFPosPhiBiasH(biasn2)) 700 if (rndm > IPDFPosPhiBiasH(biasn2)) 824 { biasn1 = biasn2; } << 701 biasn1 = biasn2; 825 else 702 else 826 { biasn3 = biasn2; } << 703 biasn3 = biasn2; 827 biasn2 = biasn1 + (biasn3 - biasn1 + 1) << 704 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 828 } << 705 } 829 bweights_t& w = bweights.Get(); << 706 bweights[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPhiBiasH(biasn2 - 1); 830 w[7] = IPDFPosPhiBiasH(biasn2) - IPDFPosPh << 707 G4double xaxisl = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); 831 G4double xaxisl = IPDFPosPhiBiasH.GetLowEd << 708 G4double xaxisu = IPDFPosPhiBiasH.GetLowEdgeEnergy(size_t(biasn2)); 832 G4double xaxisu = IPDFPosPhiBiasH.GetLowEd << 709 G4double NatProb = xaxisu - xaxisl; 833 G4double NatProb = xaxisu - xaxisl; << 710 bweights[7] = NatProb/bweights[7]; 834 w[7] = NatProb / w[7]; << 711 if(verbosityLevel >= 1) 835 if (verbosityLevel >= 1) << 712 G4cout << "PosPhi bin weight " << bweights[7] << " " << rndm << G4endl; 836 { << 713 return(IPDFPosPhiBiasH.GetEnergy(rndm)); 837 G4cout << "PosPhi bin weight " << w[7] < << 838 } 714 } 839 return (IPDFPosPhiBiasH.GetEnergy(rndm)); << 840 << 841 } 715 } >> 716 >> 717 >> 718 >> 719 >> 720 842 721