Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer << 3 // * DISCLAIMER * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers * 6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which * 7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: * 8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license * 9 // * include a list of copyright holders. << 10 // * 9 // * * 11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file << 14 // * use. * 16 // * for the full disclaimer and the limitatio << 17 // * 15 // * * 18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the * 19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. * 20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work * 21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this * 22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. * 23 // * acceptance of all terms of the Geant4 Sof << 24 // ******************************************* 21 // ******************************************************************** 25 // 22 // 26 // G4GeneralParticleSource class implementatio << 23 /////////////////////////////////////////////////////////////////////////////// 27 // 24 // 28 // Author: Fan Lei, QinetiQ ltd - 05/02/2004 << 25 // MODULE: G4GeneralParticleSource.cc 29 // Customer: ESA/ESTEC << 26 // 30 // Version: 2.0 << 27 // Version: 1.1 31 // ------------------------------------------- << 28 // Date: 18/10/00 32 << 29 // Author: C Ferguson, F Lei, P Truscott >> 30 // Organisation: University of Southampton / DERA >> 31 // Customer: ESA/ESTEC >> 32 // >> 33 // Documentation avaialable at http://www.space.dera.gov.uk/space_env/gspm.html >> 34 // These include: >> 35 // User Requirement Document (URD) >> 36 // Software Specification Documents (SSD) >> 37 // Software User Manual (SUM) >> 38 // Technical Note (TN) on the physics and algorithms >> 39 // >> 40 /////////////////////////////////////////////////////////////////////////////// >> 41 // $Id: G4GeneralParticleSource.cc,v 1.14 2001/10/19 16:48:28 flei Exp $ >> 42 // GEANT4 tag $Name: geant4-05-00 $ >> 43 /////////////////////////////////////////////////////////////////////////////// >> 44 // >> 45 // CHANGE HISTORY >> 46 // -------------- >> 47 // Further changes is recorded in the CVS HIstory file >> 48 // >> 49 // 9-May-2001 F. Lei >> 50 // added all the g4pariclegun commands >> 51 // >> 52 // 26-01-2001, F. Lei >> 53 // replace: >> 54 // posphi = acos(tx/sin(posthe)); >> 55 // >> 56 // with: >> 57 // if (posthe != 0. && posthe != pi) >> 58 // posphi = acos(tx/sin(posthe)); >> 59 // else >> 60 // posphi = 0.0; >> 61 // >> 62 // >> 63 // 10-Nov-2000, F. Lei >> 64 // some bug fixing: >> 65 // i) dclared ' G4int count' in all ****Interpolation functions >> 66 // ii) added ' return (0.) ' to GenerateUserDefTheta and GenerateUserDefPhi >> 67 // as default. >> 68 // >> 69 // 09-Nov-2000, F Lei >> 70 // Changed to a fast implementation for generating iso and cos angular directions >> 71 // >> 72 // Version 1.1, 18 October 2000, Modified to inherit from G4VPrimaryGenerator. >> 73 // New name at the request of M. Asai. >> 74 // >> 75 // >> 76 // Version 1.0, 28 February 2000, C Ferguson, Created. >> 77 // >> 78 /////////////////////////////////////////////////////////////////////////////// >> 79 // >> 80 #include "G4PrimaryParticle.hh" 33 #include "G4Event.hh" 81 #include "G4Event.hh" 34 #include "Randomize.hh" 82 #include "Randomize.hh" >> 83 #include <math.h> >> 84 #include "G4TransportationManager.hh" >> 85 #include "G4VPhysicalVolume.hh" >> 86 #include "G4PhysicalVolumeStore.hh" >> 87 #include "G4ParticleTable.hh" >> 88 #include "G4ParticleDefinition.hh" >> 89 #include "G4IonTable.hh" >> 90 #include "G4Ions.hh" >> 91 #include "G4TrackingManager.hh" >> 92 #include "G4Track.hh" 35 #include "G4GeneralParticleSource.hh" 93 #include "G4GeneralParticleSource.hh" 36 #include "G4SingleParticleSource.hh" << 37 #include "G4UnitsTable.hh" << 38 94 39 #include "G4GeneralParticleSourceData.hh" << 95 G4GeneralParticleSource::G4GeneralParticleSource() >> 96 { >> 97 // Initialise all variables >> 98 // Position distribution Variables >> 99 >> 100 NumberOfParticlesToBeGenerated = 1; >> 101 particle_definition = NULL; >> 102 G4ThreeVector zero; >> 103 particle_momentum_direction = G4ParticleMomentum(1,0,0); >> 104 particle_energy = 1.0*MeV; >> 105 particle_position = zero; >> 106 particle_time = 0.0; >> 107 particle_polarization = zero; >> 108 particle_charge = 0.0; >> 109 >> 110 SourcePosType = "Point"; >> 111 Shape = "NULL"; >> 112 halfx = 0.; >> 113 halfy = 0.; >> 114 halfz = 0.; >> 115 Radius = 0.; >> 116 Radius0 = 0.; >> 117 SR = 0.; >> 118 SX = 0.; >> 119 SY = 0.; >> 120 ParAlpha = 0.; >> 121 ParTheta = 0.; >> 122 ParPhi = 0.; >> 123 CentreCoords = G4ThreeVector(0., 0., 0.); >> 124 Rotx = HepXHat; >> 125 Roty = HepYHat; >> 126 Rotz = HepZHat; >> 127 Confine = false; //If true confines source distribution to VolName >> 128 VolName = "NULL"; >> 129 SideRefVec1 = HepXHat; // x-axis >> 130 SideRefVec2 = HepYHat; // y-axis >> 131 SideRefVec3 = HepZHat; // z-axis >> 132 >> 133 // Angular distribution variables. >> 134 AngDistType = "planar"; >> 135 AngRef1 = HepXHat; >> 136 AngRef2 = HepYHat; >> 137 AngRef3 = HepZHat; >> 138 MinTheta = 0.; >> 139 MaxTheta = pi; >> 140 MinPhi = 0.; >> 141 MaxPhi = twopi; >> 142 DR = 0.; >> 143 DX = 0.; >> 144 DY = 0.; >> 145 UserDistType = "NULL"; >> 146 IPDFThetaExist = false; >> 147 IPDFPhiExist = false; >> 148 >> 149 // Energy Distribution variables >> 150 EnergyDisType = "Mono"; >> 151 MonoEnergy = 1*MeV; >> 152 Emin = 0.; >> 153 Emax = 0.; >> 154 alpha = 0.; >> 155 Ezero = 0.; >> 156 SE = 0.; >> 157 Temp = 0.; >> 158 grad = 0.; >> 159 cept = 0.; >> 160 EnergySpec = true; // true - energy spectra, false - momentum spectra >> 161 DiffSpec = true; // true - differential spec, false integral spec >> 162 IntType = "NULL"; // Interpolation type >> 163 UserWRTSurface = false; // Any user-defined distribution is wrt co-ordinate >> 164 UserAngRef = false; // angular distribution is to mother or surface-normal co-or. >> 165 IPDFEnergyExist = false; >> 166 IPDFArbExist = false; >> 167 >> 168 // Bias variables >> 169 XBias = false; >> 170 IPDFXBias = false; >> 171 YBias = false; >> 172 IPDFYBias = false; >> 173 ZBias = false; >> 174 IPDFZBias = false; >> 175 ThetaBias = false; >> 176 IPDFThetaBias = false; >> 177 PhiBias = false; >> 178 IPDFPhiBias = false; >> 179 EnergyBias = false; >> 180 IPDFEnergyBias = false; >> 181 >> 182 ArbEmin = 0.; >> 183 ArbEmax = 0.; >> 184 >> 185 // verbosity >> 186 verbosityLevel = 0; >> 187 >> 188 theMessenger = new G4GeneralParticleSourceMessenger(this); >> 189 gNavigator = G4TransportationManager::GetTransportationManager() >> 190 ->GetNavigatorForTracking(); >> 191 } 40 192 41 #include "G4Threading.hh" << 193 G4GeneralParticleSource::~G4GeneralParticleSource() 42 #include "G4AutoLock.hh" << 194 { >> 195 delete theMessenger; >> 196 } 43 197 44 namespace << 198 void G4GeneralParticleSource::SetPosDisType(G4String PosType) 45 { 199 { 46 G4Mutex messangerInit = G4MUTEX_INITIALIZER; << 200 SourcePosType = PosType; 47 } 201 } 48 202 49 G4GeneralParticleSource::G4GeneralParticleSour << 203 void G4GeneralParticleSource::SetPosDisShape(G4String shapeType) 50 { 204 { 51 GPSData = G4GeneralParticleSourceData::Insta << 205 Shape = shapeType; 52 // currentSource = GPSData->GetCurrentSource << 206 } 53 // currentSourceIdx = G4int(GPSData->GetSour << 54 207 55 // Messenger is special, only a worker shoul << 208 void G4GeneralParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre) 56 // Singleton pattern << 209 { 57 // << 210 CentreCoords = coordsOfCentre; 58 theMessenger = G4GeneralParticleSourceMessen << 211 } 59 212 60 // Some initialization should be done only o << 213 void G4GeneralParticleSource::SetPosRot1(G4ThreeVector posrot1) 61 // << 214 { 62 G4AutoLock l(&messangerInit); << 215 // This should be x' 63 static G4bool onlyOnce = false; << 216 Rotx = posrot1; 64 if ( !onlyOnce ) << 217 if(verbosityLevel == 2) 65 { << 218 { 66 theMessenger->SetParticleGun(GPSData->GetC << 219 G4cout << "Vector x' " << Rotx << G4endl; 67 IntensityNormalization(); << 220 } 68 onlyOnce = true; << 221 GenerateRotationMatrices(); >> 222 } >> 223 >> 224 void G4GeneralParticleSource::SetPosRot2(G4ThreeVector posrot2) >> 225 { >> 226 // This is a vector in the plane x'y' but need not >> 227 // be y' >> 228 Roty = posrot2; >> 229 if(verbosityLevel == 2) >> 230 { >> 231 G4cout << "The vector in the x'-y' plane " << Roty << G4endl; >> 232 } >> 233 GenerateRotationMatrices(); >> 234 } >> 235 >> 236 void G4GeneralParticleSource::SetHalfX(G4double xhalf) >> 237 { >> 238 halfx = xhalf; >> 239 } >> 240 >> 241 void G4GeneralParticleSource::SetHalfY(G4double yhalf) >> 242 { >> 243 halfy = yhalf; >> 244 } >> 245 >> 246 void G4GeneralParticleSource::SetHalfZ(G4double zhalf) >> 247 { >> 248 halfz = zhalf; >> 249 } >> 250 >> 251 void G4GeneralParticleSource::SetRadius(G4double rad) >> 252 { >> 253 Radius = rad; >> 254 } >> 255 >> 256 void G4GeneralParticleSource::SetRadius0(G4double rad) >> 257 { >> 258 Radius0 = rad; >> 259 } >> 260 >> 261 void G4GeneralParticleSource::SetBeamSigmaInR(G4double r) >> 262 { >> 263 SR = r; >> 264 } >> 265 >> 266 void G4GeneralParticleSource::SetBeamSigmaInX(G4double r) >> 267 { >> 268 SX = r; >> 269 } >> 270 >> 271 void G4GeneralParticleSource::SetBeamSigmaInY(G4double r) >> 272 { >> 273 SY = r; >> 274 } >> 275 >> 276 void G4GeneralParticleSource::SetParAlpha(G4double paralp) >> 277 { >> 278 ParAlpha = paralp; >> 279 } >> 280 >> 281 void G4GeneralParticleSource::SetParTheta(G4double parthe) >> 282 { >> 283 ParTheta = parthe; >> 284 } >> 285 >> 286 void G4GeneralParticleSource::SetParPhi(G4double parphi) >> 287 { >> 288 ParPhi = parphi; >> 289 } >> 290 >> 291 void G4GeneralParticleSource::GenerateRotationMatrices() >> 292 { >> 293 // This takes in 2 vectors, x' and one in the plane x'-y', >> 294 // and from these takes a cross product to calculate z'. >> 295 // Then a cross product is taken between x' and z' to give >> 296 // y'. >> 297 Rotx = Rotx.unit(); // x' >> 298 Roty = Roty.unit(); // vector in x'y' plane >> 299 Rotz = Rotx.cross(Roty); // z' >> 300 Roty = Rotz.cross(Rotx); // y' >> 301 if(verbosityLevel == 2) >> 302 { >> 303 G4cout << "The new axes, x', y', z' " << Rotx << " " << Roty << " " << Rotz << G4endl; >> 304 } >> 305 } >> 306 >> 307 void G4GeneralParticleSource::ConfineSourceToVolume(G4String Vname) >> 308 { >> 309 VolName = Vname; >> 310 if(verbosityLevel == 2) >> 311 G4cout << VolName << G4endl; >> 312 G4VPhysicalVolume *tempPV = NULL; >> 313 G4PhysicalVolumeStore *PVStore = 0; >> 314 G4String theRequiredVolumeName = VolName; >> 315 PVStore = G4PhysicalVolumeStore::GetInstance(); >> 316 G4int i = 0; >> 317 G4bool found = false; >> 318 if(verbosityLevel == 2) >> 319 G4cout << PVStore->size() << G4endl; >> 320 while (!found && i<G4int(PVStore->size())) { >> 321 tempPV = (*PVStore)[i]; >> 322 found = tempPV->GetName() == theRequiredVolumeName; >> 323 if(verbosityLevel == 2) >> 324 G4cout << i << " " << " " << tempPV->GetName() << " " << theRequiredVolumeName << " " << found << G4endl; >> 325 if (!found) >> 326 {i++;} 69 } 327 } >> 328 // found = true then the volume exists else it doesnt. >> 329 if(found == true) >> 330 { >> 331 if(verbosityLevel >= 1) >> 332 G4cout << "Volume " << VolName << " exists" << G4endl; >> 333 Confine = true; >> 334 } >> 335 else >> 336 { >> 337 G4cout << " **** Error: Volume does not exist **** " << G4endl; >> 338 G4cout << " Ignoring confine condition" << G4endl; >> 339 Confine = false; >> 340 VolName = "NULL"; >> 341 } >> 342 70 } 343 } 71 344 72 G4GeneralParticleSource::~G4GeneralParticleSou << 345 void G4GeneralParticleSource::GeneratePointSource() 73 { 346 { 74 theMessenger->Destroy(); << 347 // Generates Points given the point source. >> 348 if(SourcePosType == "Point") >> 349 particle_position = CentreCoords; >> 350 else >> 351 if(verbosityLevel >= 1) >> 352 G4cout << "Error SourcePosType is not set to Point" << G4endl; 75 } 353 } 76 354 77 void G4GeneralParticleSource::AddaSource(G4dou << 355 void G4GeneralParticleSource::GeneratePointsInBeam() 78 { 356 { 79 GPSData->Lock(); << 357 G4double x, y, z, r, phi; 80 358 81 GPSData->AddASource(aV); << 359 G4ThreeVector RandPos; 82 theMessenger->SetParticleGun(GPSData->GetCur << 360 G4double tempx, tempy, tempz; >> 361 z = 0.; >> 362 // Private Method to create points in a plane >> 363 if(Shape == "Circle") >> 364 { >> 365 r = G4RandGauss::shoot(0.0,SR) ; >> 366 phi = twopi * G4UniformRand(); >> 367 x = r*cos(phi) ; >> 368 y = r*sin(phi) ; >> 369 } >> 370 else >> 371 { >> 372 // all other cases default to Rectangle case >> 373 x = G4RandGauss::shoot(0.0,SX); >> 374 y = G4RandGauss::shoot(0.0,SY); >> 375 } >> 376 // Apply Rotation Matrix >> 377 // x * Rotx, y * Roty and z * Rotz >> 378 if(verbosityLevel == 2) >> 379 { >> 380 G4cout << "Raw position " << x << "," << y << "," << z << G4endl; >> 381 } >> 382 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); >> 383 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); >> 384 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); 83 385 84 // TODO: But do we really normalize here aft << 386 RandPos.setX(tempx); 85 IntensityNormalization(); << 387 RandPos.setY(tempy); >> 388 RandPos.setZ(tempz); 86 389 87 GPSData->Unlock(); << 390 // Translate >> 391 particle_position = CentreCoords + RandPos; >> 392 if(verbosityLevel >= 1) >> 393 { >> 394 if(verbosityLevel == 2) >> 395 { >> 396 G4cout << "Rotated Position " << RandPos << G4endl; >> 397 } >> 398 G4cout << "Rotated and Translated position " << particle_position << G4endl; >> 399 } 88 } 400 } 89 401 90 void G4GeneralParticleSource::IntensityNormali << 402 void G4GeneralParticleSource::GeneratePointsInPlane() 91 { 403 { 92 GPSData->IntensityNormalise(); << 404 G4double x, y, z; 93 normalised=GPSData->Normalised(); << 405 G4double expression; >> 406 G4ThreeVector RandPos; >> 407 G4double tempx, tempy, tempz; >> 408 x = y = z = 0.; >> 409 >> 410 if(SourcePosType != "Plane" && verbosityLevel >= 1) >> 411 G4cout << "Error: SourcePosType is not Plane" << G4endl; >> 412 >> 413 // Private Method to create points in a plane >> 414 if(Shape == "Circle") >> 415 { >> 416 x = Radius + 100.; >> 417 y = Radius + 100.; >> 418 while(sqrt((x*x) + (y*y)) > Radius) >> 419 { >> 420 x = GenRandX(); >> 421 y = GenRandY(); >> 422 >> 423 x = (x*2.*Radius) - Radius; >> 424 y = (y*2.*Radius) - Radius; >> 425 } >> 426 } >> 427 else if(Shape == "Annulus") >> 428 { >> 429 x = Radius + 100.; >> 430 y = Radius + 100.; >> 431 while(sqrt((x*x) + (y*y)) > Radius || sqrt((x*x) + (y*y)) < Radius0 ) >> 432 { >> 433 x = GenRandX(); >> 434 y = GenRandY(); >> 435 >> 436 x = (x*2.*Radius) - Radius; >> 437 y = (y*2.*Radius) - Radius; >> 438 } >> 439 } >> 440 else if(Shape == "Ellipse") >> 441 { >> 442 expression = 20.; >> 443 while(expression > 1.) >> 444 { >> 445 x = GenRandX(); >> 446 y = GenRandY(); >> 447 >> 448 x = (x*2.*halfx) - halfx; >> 449 y = (y*2.*halfy) - halfy; >> 450 >> 451 expression = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy)); >> 452 } >> 453 } >> 454 else if(Shape == "Square") >> 455 { >> 456 x = GenRandX(); >> 457 y = GenRandY(); >> 458 x = (x*2.*halfx) - halfx; >> 459 y = (y*2.*halfy) - halfy; >> 460 } >> 461 else if(Shape == "Rectangle") >> 462 { >> 463 x = GenRandX(); >> 464 y = GenRandY(); >> 465 x = (x*2.*halfx) - halfx; >> 466 y = (y*2.*halfy) - halfy; >> 467 } >> 468 else >> 469 G4cout << "Shape not one of the plane types" << G4endl; >> 470 >> 471 // Apply Rotation Matrix >> 472 // x * Rotx, y * Roty and z * Rotz >> 473 if(verbosityLevel == 2) >> 474 { >> 475 G4cout << "Raw position " << x << "," << y << "," << z << G4endl; >> 476 } >> 477 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); >> 478 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); >> 479 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); >> 480 >> 481 RandPos.setX(tempx); >> 482 RandPos.setY(tempy); >> 483 RandPos.setZ(tempz); >> 484 >> 485 // Translate >> 486 particle_position = CentreCoords + RandPos; >> 487 if(verbosityLevel >= 1) >> 488 { >> 489 if(verbosityLevel == 2) >> 490 { >> 491 G4cout << "Rotated Position " << RandPos << G4endl; >> 492 } >> 493 G4cout << "Rotated and Translated position " << particle_position << G4endl; >> 494 } >> 495 >> 496 // For Cosine-Law make SideRefVecs = to Rotation matrix vectors >> 497 SideRefVec1 = Rotx; >> 498 SideRefVec2 = Roty; >> 499 SideRefVec3 = Rotz; >> 500 // If rotation matrix z' point to origin then invert the matrix >> 501 // So that SideRefVecs point away. >> 502 if((CentreCoords.x() > 0. && Rotz.x() < 0.) >> 503 || (CentreCoords.x() < 0. && Rotz.x() > 0.) >> 504 || (CentreCoords.y() > 0. && Rotz.y() < 0.) >> 505 || (CentreCoords.y() < 0. && Rotz.y() > 0.) >> 506 || (CentreCoords.z() > 0. && Rotz.z() < 0.) >> 507 || (CentreCoords.z() < 0. && Rotz.z() > 0.)) >> 508 { >> 509 // Invert y and z. >> 510 SideRefVec2 = -SideRefVec2; >> 511 SideRefVec3 = -SideRefVec3; >> 512 } >> 513 if(verbosityLevel == 2) >> 514 { >> 515 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; >> 516 } 94 } 517 } 95 518 96 void G4GeneralParticleSource::ListSource() << 519 void G4GeneralParticleSource::GeneratePointsOnSurface() 97 { 520 { 98 G4cout << "The number of particle sources is << 521 //Private method to create points on a surface 99 << GPSData->GetIntensityVectorSize() << 522 G4double theta, phi; 100 G4cout << " Multiple Vertex sources: " << GP << 523 G4double x, y, z; 101 G4cout << " Flat Sampling flag: " << GPSData << 524 x = y = z = 0.; 102 const G4int currentIdx = GPSData->GetCurrent << 525 G4ThreeVector RandPos; 103 for(G4int i=0; i<GPSData->GetIntensityVector << 526 // G4double tempx, tempy, tempz; 104 { << 527 105 G4cout << "\tsource " << i << " with inten << 528 if(SourcePosType != "Surface" && verbosityLevel >= 1) 106 << GPSData->GetIntensity(i) << G4en << 529 G4cout << "Error SourcePosType not Surface" << G4endl; 107 const G4SingleParticleSource* thisSrc = GP << 530 108 G4cout << " \t\tNum Particles: "<<thisSrc- << 531 if(Shape == "Sphere") 109 << "; Particle type: " << 532 { 110 << thisSrc->GetParticleDefinition() << 533 G4double tantheta; 111 G4cout << " \t\tEnergy: " << 534 theta = GenRandTheta(); 112 << G4BestUnit(thisSrc->GetEneDist()->GetM << 535 phi = GenRandPhi(); 113 G4cout << " \t\tDirection: " << 536 114 << thisSrc->GetAngDist()->GetDirect << 537 theta = acos(1. - 2.*theta); // theta isotropic 115 G4cout << G4BestUnit(thisSrc->GetPosDist() << 538 phi = phi * 2. * pi; 116 << G4endl; << 539 tantheta = tan(theta); 117 G4cout << " \t\tAngular Distribution: " << 540 118 << thisSrc->GetAngDist()->GetDistTy << 541 x = Radius * sin(theta) * cos(phi); 119 G4cout << " \t\tEnergy Distribution: " << 542 y = Radius * sin(theta) * sin(phi); 120 << thisSrc->GetEneDist()->GetEnergy << 543 z = Radius * cos(theta); 121 G4cout << " \t\tPosition Distribution Type << 544 122 << thisSrc->GetPosDist()->GetPosDis << 545 RandPos.setX(x); 123 G4cout << "; Position Shape: " << 546 RandPos.setY(y); 124 << thisSrc->GetPosDist()->GetPosDis << 547 RandPos.setZ(z); >> 548 >> 549 // Cosine-law (not a good idea to use this here) >> 550 G4ThreeVector zdash(x,y,z); >> 551 zdash = zdash.unit(); >> 552 G4ThreeVector xdash = Rotz.cross(zdash); >> 553 G4ThreeVector ydash = xdash.cross(zdash); >> 554 SideRefVec1 = xdash; >> 555 SideRefVec2 = ydash; >> 556 SideRefVec3 = zdash; >> 557 } >> 558 else if(Shape == "Ellipsoid") >> 559 { >> 560 G4double theta, phi, minphi, maxphi, middlephi; >> 561 G4double answer, constant; >> 562 >> 563 constant = pi/(halfx*halfx) + pi/(halfy*halfy) + >> 564 twopi/(halfz*halfz); >> 565 >> 566 // simplified approach >> 567 theta = GenRandTheta(); >> 568 phi = GenRandPhi(); >> 569 >> 570 theta = acos(1. - 2.*theta); >> 571 minphi = 0.; >> 572 maxphi = twopi; >> 573 while(maxphi-minphi > 0.) >> 574 { >> 575 middlephi = (maxphi+minphi)/2.; >> 576 answer = (1./(halfx*halfx))*(middlephi/2. + sin(2*middlephi)/4.) >> 577 + (1./(halfy*halfy))*(middlephi/2. - sin(2*middlephi)/4.) >> 578 + middlephi/(halfz*halfz); >> 579 answer = answer/constant; >> 580 if(answer > phi) maxphi = middlephi; >> 581 if(answer < phi) minphi = middlephi; >> 582 if(fabs(answer-phi) <= 0.00001) >> 583 { >> 584 minphi = maxphi +1; >> 585 phi = middlephi; >> 586 } >> 587 } >> 588 >> 589 x = sin(theta)*cos(phi); >> 590 y = sin(theta)*sin(phi); >> 591 z = cos(theta); >> 592 // x,y and z form a unit vector. Put this onto the ellipse. >> 593 G4double lhs; >> 594 // solve for x >> 595 G4double numYinX = y/x; >> 596 G4double numZinX = z/x; >> 597 G4double tempxvar; >> 598 tempxvar= 1./(halfx*halfx)+(numYinX*numYinX)/(halfy*halfy) >> 599 + (numZinX*numZinX)/(halfz*halfz); >> 600 >> 601 tempxvar = 1./tempxvar; >> 602 G4double coordx = sqrt(tempxvar); >> 603 >> 604 //solve for y >> 605 G4double numXinY = x/y; >> 606 G4double numZinY = z/y; >> 607 G4double tempyvar; >> 608 tempyvar=(numXinY*numXinY)/(halfx*halfx)+1./(halfy*halfy) >> 609 +(numZinY*numZinY)/(halfz*halfz); >> 610 tempyvar = 1./tempyvar; >> 611 G4double coordy = sqrt(tempyvar); >> 612 >> 613 //solve for z >> 614 G4double numXinZ = x/z; >> 615 G4double numYinZ = y/z; >> 616 G4double tempzvar; >> 617 tempzvar=(numXinZ*numXinZ)/(halfx*halfx) >> 618 +(numYinZ*numYinZ)/(halfy*halfy)+1./(halfz*halfz); >> 619 tempzvar = 1./tempzvar; >> 620 G4double coordz = sqrt(tempzvar); >> 621 >> 622 lhs = sqrt((coordx*coordx)/(halfx*halfx) + >> 623 (coordy*coordy)/(halfy*halfy) + >> 624 (coordz*coordz)/(halfz*halfz)); >> 625 >> 626 if(fabs(lhs-1.) > 0.001 && verbosityLevel >= 1) >> 627 G4cout << "Error: theta, phi not really on ellipsoid" << G4endl; >> 628 >> 629 // coordx, coordy and coordz are all positive >> 630 G4double TestRandVar = G4UniformRand(); >> 631 if(TestRandVar > 0.5) >> 632 { >> 633 coordx = -coordx; >> 634 } >> 635 TestRandVar = G4UniformRand(); >> 636 if(TestRandVar > 0.5) >> 637 { >> 638 coordy = -coordy; >> 639 } >> 640 TestRandVar = G4UniformRand(); >> 641 if(TestRandVar > 0.5) >> 642 { >> 643 coordz = -coordz; >> 644 } >> 645 >> 646 RandPos.setX(coordx); >> 647 RandPos.setY(coordy); >> 648 RandPos.setZ(coordz); >> 649 >> 650 // Cosine-law (not a good idea to use this here) >> 651 G4ThreeVector zdash(coordx,coordy,coordz); >> 652 zdash = zdash.unit(); >> 653 G4ThreeVector xdash = Rotz.cross(zdash); >> 654 G4ThreeVector ydash = xdash.cross(zdash); >> 655 SideRefVec1 = xdash; >> 656 SideRefVec2 = ydash; >> 657 SideRefVec3 = zdash; >> 658 } >> 659 else if(Shape == "Cylinder") >> 660 { >> 661 G4double AreaTop, AreaBot, AreaLat; >> 662 G4double AreaTotal, prob1, prob2, prob3; >> 663 G4double testrand; >> 664 >> 665 // User giver Radius and z-half length >> 666 // Calculate surface areas, maybe move this to >> 667 // a different routine. >> 668 >> 669 AreaTop = pi * Radius * Radius; >> 670 AreaBot = AreaTop; >> 671 AreaLat = 2. * pi * Radius * 2. * halfz; >> 672 AreaTotal = AreaTop + AreaBot + AreaLat; >> 673 >> 674 prob1 = AreaTop / AreaTotal; >> 675 prob2 = AreaBot / AreaTotal; >> 676 prob3 = 1.00 - prob1 - prob2; >> 677 if(fabs(prob3 - (AreaLat/AreaTotal)) >= 0.001) >> 678 { >> 679 if(verbosityLevel >= 1) >> 680 G4cout << AreaLat/AreaTotal << " " << prob3<<G4endl; >> 681 G4cout << "Error in prob3" << G4endl; >> 682 } >> 683 >> 684 // Decide surface to calculate point on. >> 685 >> 686 testrand = G4UniformRand(); >> 687 if(testrand <= prob1) >> 688 { >> 689 //Point on Top surface >> 690 z = halfz; >> 691 x = Radius + 100.; >> 692 y = Radius + 100.; >> 693 while(((x*x)+(y*y)) > (Radius*Radius)) >> 694 { >> 695 x = GenRandX(); >> 696 y = GenRandY(); >> 697 >> 698 x = x * 2. * Radius; >> 699 y = y * 2. * Radius; >> 700 x = x - Radius; >> 701 y = y - Radius; >> 702 } >> 703 // Cosine law >> 704 SideRefVec1 = Rotx; >> 705 SideRefVec2 = Roty; >> 706 SideRefVec3 = Rotz; >> 707 } >> 708 else if((testrand > prob1) && (testrand <= (prob1 + prob2))) >> 709 { >> 710 //Point on Bottom surface >> 711 z = -halfz; >> 712 x = Radius + 100.; >> 713 y = Radius + 100.; >> 714 while(((x*x)+(y*y)) > (Radius*Radius)) >> 715 { >> 716 x = GenRandX(); >> 717 y = GenRandY(); >> 718 >> 719 x = x * 2. * Radius; >> 720 y = y * 2. * Radius; >> 721 x = x - Radius; >> 722 y = y - Radius; >> 723 } >> 724 // Cosine law >> 725 SideRefVec1 = Rotx; >> 726 SideRefVec2 = -Roty; >> 727 SideRefVec3 = -Rotz; >> 728 } >> 729 else if(testrand > (prob1+prob2)) >> 730 { >> 731 G4double rand; >> 732 //Point on Lateral Surface >> 733 >> 734 rand = G4UniformRand(); >> 735 rand = rand * 2. * pi; >> 736 >> 737 x = Radius * cos(rand); >> 738 y = Radius * sin(rand); >> 739 >> 740 z = GenRandZ(); >> 741 >> 742 z = z * 2. * halfz; >> 743 z = z - halfz; >> 744 >> 745 // Cosine law >> 746 G4ThreeVector zdash(x,y,0.); >> 747 zdash = zdash.unit(); >> 748 G4ThreeVector xdash = Rotz.cross(zdash); >> 749 G4ThreeVector ydash = xdash.cross(zdash); >> 750 SideRefVec1 = xdash; >> 751 SideRefVec2 = ydash; >> 752 SideRefVec3 = zdash; >> 753 } >> 754 else >> 755 G4cout << "Error: testrand " << testrand << G4endl; >> 756 >> 757 RandPos.setX(x); >> 758 RandPos.setY(y); >> 759 RandPos.setZ(z); >> 760 >> 761 } >> 762 else if(Shape == "Para") >> 763 { >> 764 G4double testrand; >> 765 //Right Parallelepiped. >> 766 // User gives x,y,z half lengths and ParAlpha >> 767 // ParTheta and ParPhi >> 768 // +x = <1, -x >1 & <2, +y >2 & <3, -y >3 &<4 >> 769 // +z >4 & < 5, -z >5 &<6. >> 770 testrand = G4UniformRand(); >> 771 G4double AreaX = halfy * halfz * 4.; >> 772 G4double AreaY = halfx * halfz * 4.; >> 773 G4double AreaZ = halfx * halfy * 4.; >> 774 G4double AreaTotal = 2*(AreaX + AreaY + AreaZ); >> 775 G4double Probs[6]; >> 776 Probs[0] = AreaX/AreaTotal; >> 777 Probs[1] = Probs[0] + AreaX/AreaTotal; >> 778 Probs[2] = Probs[1] + AreaY/AreaTotal; >> 779 Probs[3] = Probs[2] + AreaY/AreaTotal; >> 780 Probs[4] = Probs[3] + AreaZ/AreaTotal; >> 781 Probs[5] = Probs[4] + AreaZ/AreaTotal; >> 782 >> 783 x = GenRandX(); >> 784 y = GenRandY(); >> 785 z = GenRandZ(); >> 786 >> 787 x = x * halfx * 2.; >> 788 x = x - halfx; >> 789 y = y * halfy * 2.; >> 790 y = y - halfy; >> 791 z = z * halfz * 2.; >> 792 z = z - halfz; >> 793 // Pick a side first >> 794 if(testrand < Probs[0]) >> 795 { >> 796 // side is +x >> 797 x = halfx + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha); >> 798 y = y + z*tan(ParTheta)*sin(ParPhi); >> 799 z = z; >> 800 // Cosine-law >> 801 G4ThreeVector xdash(halfz*tan(ParTheta)*cos(ParPhi), >> 802 halfz*tan(ParTheta)*sin(ParPhi), >> 803 halfz/cos(ParPhi)); >> 804 G4ThreeVector ydash(halfy*tan(ParAlpha), halfy, 0.0); >> 805 xdash = xdash.unit(); >> 806 ydash = ydash.unit(); >> 807 G4ThreeVector zdash = xdash.cross(ydash); >> 808 SideRefVec1 = xdash; >> 809 SideRefVec2 = ydash; >> 810 SideRefVec3 = zdash; >> 811 } >> 812 else if(testrand >= Probs[0] && testrand < Probs[1]) >> 813 { >> 814 // side is -x >> 815 x = -halfx + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha); >> 816 y = y + z*tan(ParTheta)*sin(ParPhi); >> 817 z = z; >> 818 // Cosine-law >> 819 G4ThreeVector xdash(halfz*tan(ParTheta)*cos(ParPhi), >> 820 halfz*tan(ParTheta)*sin(ParPhi), >> 821 halfz/cos(ParPhi)); >> 822 G4ThreeVector ydash(halfy*tan(ParAlpha), halfy, 0.0); >> 823 xdash = xdash.unit(); >> 824 ydash = ydash.unit(); >> 825 G4ThreeVector zdash = xdash.cross(ydash); >> 826 SideRefVec1 = xdash; >> 827 SideRefVec2 = -ydash; >> 828 SideRefVec3 = -zdash; >> 829 } >> 830 else if(testrand >= Probs[1] && testrand < Probs[2]) >> 831 { >> 832 // side is +y >> 833 x = x + z*tan(ParTheta)*cos(ParPhi) + halfy*tan(ParAlpha); >> 834 y = halfy + z*tan(ParTheta)*sin(ParPhi); >> 835 z = z; >> 836 // Cosine-law >> 837 G4ThreeVector ydash(halfz*tan(ParTheta)*cos(ParPhi), >> 838 halfz*tan(ParTheta)*sin(ParPhi), >> 839 halfz/cos(ParPhi)); >> 840 ydash = ydash.unit(); >> 841 G4ThreeVector xdash = Roty.cross(ydash); >> 842 G4ThreeVector zdash = xdash.cross(ydash); >> 843 SideRefVec1 = xdash; >> 844 SideRefVec2 = ydash; >> 845 SideRefVec3 = zdash; >> 846 } >> 847 else if(testrand >= Probs[2] && testrand < Probs[3]) >> 848 { >> 849 // side is -y >> 850 x = x + z*tan(ParTheta)*cos(ParPhi) - halfy*tan(ParAlpha); >> 851 y = -halfy + z*tan(ParTheta)*sin(ParPhi); >> 852 z = z; >> 853 // Cosine-law >> 854 G4ThreeVector ydash(halfz*tan(ParTheta)*cos(ParPhi), >> 855 halfz*tan(ParTheta)*sin(ParPhi), >> 856 halfz/cos(ParPhi)); >> 857 ydash = ydash.unit(); >> 858 G4ThreeVector xdash = Roty.cross(ydash); >> 859 G4ThreeVector zdash = xdash.cross(ydash); >> 860 SideRefVec1 = xdash; >> 861 SideRefVec2 = -ydash; >> 862 SideRefVec3 = -zdash; >> 863 } >> 864 else if(testrand >= Probs[3] && testrand < Probs[4]) >> 865 { >> 866 // side is +z >> 867 z = halfz; >> 868 y = y + halfz*sin(ParPhi)*tan(ParTheta); >> 869 x = x + halfz*cos(ParPhi)*tan(ParTheta) + y*tan(ParAlpha); >> 870 // Cosine-law >> 871 SideRefVec1 = Rotx; >> 872 SideRefVec2 = Roty; >> 873 SideRefVec3 = Rotz; >> 874 } >> 875 else if(testrand >= Probs[4] && testrand < Probs[5]) >> 876 { >> 877 // side is -z >> 878 z = -halfz; >> 879 y = y - halfz*sin(ParPhi)*tan(ParTheta); >> 880 x = x - halfz*cos(ParPhi)*tan(ParTheta) + y*tan(ParAlpha); >> 881 // Cosine-law >> 882 SideRefVec1 = Rotx; >> 883 SideRefVec2 = -Roty; >> 884 SideRefVec3 = -Rotz; >> 885 } >> 886 else >> 887 { >> 888 G4cout << "Error: testrand out of range" << G4endl; >> 889 if(verbosityLevel >= 1) >> 890 G4cout << "testrand=" << testrand << " Probs[5]=" << Probs[5] <<G4endl; >> 891 } >> 892 >> 893 RandPos.setX(x); >> 894 RandPos.setY(y); >> 895 RandPos.setZ(z); >> 896 } >> 897 >> 898 // Apply Rotation Matrix >> 899 // x * Rotx, y * Roty and z * Rotz >> 900 if(verbosityLevel == 2) >> 901 G4cout << "Raw position " << RandPos << G4endl; >> 902 >> 903 x=(RandPos.x()*Rotx.x())+(RandPos.y()*Roty.x())+(RandPos.z()*Rotz.x()); >> 904 y=(RandPos.x()*Rotx.y())+(RandPos.y()*Roty.y())+(RandPos.z()*Rotz.y()); >> 905 z=(RandPos.x()*Rotx.z())+(RandPos.y()*Roty.z())+(RandPos.z()*Rotz.z()); >> 906 >> 907 RandPos.setX(x); >> 908 RandPos.setY(y); >> 909 RandPos.setZ(z); >> 910 >> 911 // Translate >> 912 particle_position = CentreCoords + RandPos; >> 913 >> 914 if(verbosityLevel >= 1) >> 915 { >> 916 if(verbosityLevel == 2) >> 917 G4cout << "Rotated position " << RandPos << G4endl; >> 918 G4cout << "Rotated and translated position " << particle_position << G4endl; >> 919 } >> 920 if(verbosityLevel == 2) >> 921 { >> 922 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; >> 923 } >> 924 } >> 925 >> 926 void G4GeneralParticleSource::GeneratePointsInVolume() >> 927 { >> 928 G4ThreeVector RandPos; >> 929 G4double tempx, tempy, tempz; >> 930 G4double x, y, z; >> 931 x = y = z = 0.; >> 932 if(SourcePosType != "Volume" && verbosityLevel >= 1) >> 933 G4cout << "Error SourcePosType not Volume" << G4endl; >> 934 //Private method to create points in a volume >> 935 if(Shape == "Sphere") >> 936 { >> 937 x = Radius*2.; >> 938 y = Radius*2.; >> 939 z = Radius*2.; >> 940 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) >> 941 { >> 942 x = GenRandX(); >> 943 y = GenRandY(); >> 944 z = GenRandZ(); >> 945 >> 946 x = (x*2.*Radius) - Radius; >> 947 y = (y*2.*Radius) - Radius; >> 948 z = (z*2.*Radius) - Radius; >> 949 } >> 950 } >> 951 else if(Shape == "Ellipsoid") >> 952 { >> 953 G4double temp; >> 954 temp = 100.; >> 955 while(temp > 1.) >> 956 { >> 957 x = GenRandX(); >> 958 y = GenRandY(); >> 959 z = GenRandZ(); >> 960 >> 961 x = (x*2.*halfx) - halfx; >> 962 y = (y*2.*halfy) - halfy; >> 963 z = (z*2.*halfz) - halfz; >> 964 >> 965 temp = ((x*x)/(halfx*halfx)) + ((y*y)/(halfy*halfy)) >> 966 + ((z*z)/(halfz*halfz)); >> 967 } >> 968 } >> 969 else if(Shape == "Cylinder") >> 970 { >> 971 x = Radius*2.; >> 972 y = Radius*2.; >> 973 while(((x*x)+(y*y)) > (Radius*Radius)) >> 974 { >> 975 x = GenRandX(); >> 976 y = GenRandY(); >> 977 z = GenRandZ(); >> 978 >> 979 x = (x*2.*Radius) - Radius; >> 980 y = (y*2.*Radius) - Radius; >> 981 z = (z*2.*halfz) - halfz; >> 982 } >> 983 } >> 984 else if(Shape == "Para") >> 985 { >> 986 x = GenRandX(); >> 987 y = GenRandY(); >> 988 z = GenRandZ(); >> 989 x = (x*2.*halfx) - halfx; >> 990 y = (y*2.*halfy) - halfy; >> 991 z = (z*2.*halfz) - halfz; >> 992 x = x + z*tan(ParTheta)*cos(ParPhi) + y*tan(ParAlpha); >> 993 y = y + z*tan(ParTheta)*sin(ParPhi); >> 994 z = z; >> 995 } >> 996 else >> 997 G4cout << "Error: Volume Shape Doesnt Exist" << G4endl; >> 998 >> 999 RandPos.setX(x); >> 1000 RandPos.setY(y); >> 1001 RandPos.setZ(z); >> 1002 >> 1003 // Apply Rotation Matrix >> 1004 // x * Rotx, y * Roty and z * Rotz >> 1005 tempx = (x * Rotx.x()) + (y * Roty.x()) + (z * Rotz.x()); >> 1006 tempy = (x * Rotx.y()) + (y * Roty.y()) + (z * Rotz.y()); >> 1007 tempz = (x * Rotx.z()) + (y * Roty.z()) + (z * Rotz.z()); >> 1008 >> 1009 RandPos.setX(tempx); >> 1010 RandPos.setY(tempy); >> 1011 RandPos.setZ(tempz); >> 1012 >> 1013 // Translate >> 1014 particle_position = CentreCoords + RandPos; >> 1015 >> 1016 if(verbosityLevel == 2) >> 1017 { >> 1018 G4cout << "Raw position " << x << "," << y << "," << z << G4endl; >> 1019 G4cout << "Rotated position " << RandPos << G4endl; >> 1020 } >> 1021 if(verbosityLevel >= 1) >> 1022 G4cout << "Rotated and translated position " << particle_position << G4endl; >> 1023 >> 1024 // Cosine-law (not a good idea to use this here) >> 1025 G4ThreeVector zdash(tempx,tempy,tempz); >> 1026 zdash = zdash.unit(); >> 1027 G4ThreeVector xdash = Rotz.cross(zdash); >> 1028 G4ThreeVector ydash = xdash.cross(zdash); >> 1029 SideRefVec1 = xdash; >> 1030 SideRefVec2 = ydash; >> 1031 SideRefVec3 = zdash; >> 1032 >> 1033 if(verbosityLevel == 2) >> 1034 { >> 1035 G4cout << "Reference vectors for cosine-law " << SideRefVec1 << " " << SideRefVec2 << " " << SideRefVec3 << G4endl; >> 1036 } >> 1037 } >> 1038 >> 1039 G4bool G4GeneralParticleSource::IsSourceConfined() >> 1040 { >> 1041 // Method to check point is within the volume specified >> 1042 if(Confine == false) >> 1043 G4cout << "Error: Confine is false" << G4endl; >> 1044 G4ThreeVector null(0.,0.,0.); >> 1045 G4ThreeVector *ptr; >> 1046 ptr = &null; >> 1047 >> 1048 // Check particle_position is within VolName, if so true, >> 1049 // else false >> 1050 G4VPhysicalVolume *theVolume; >> 1051 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true); >> 1052 G4String theVolName = theVolume->GetName(); >> 1053 if(theVolName == VolName) >> 1054 { >> 1055 if(verbosityLevel >= 1) >> 1056 G4cout << "Particle is in volume " << VolName << G4endl; >> 1057 return(true); >> 1058 } >> 1059 else >> 1060 return(false); >> 1061 } >> 1062 >> 1063 // Now follows the angular distribution methods. >> 1064 >> 1065 void G4GeneralParticleSource::SetAngDistType(G4String atype) >> 1066 { >> 1067 if(atype != "iso" && atype != "cos" && atype != "user" && atype != "planar" >> 1068 && atype != "beam1d" && atype != "beam2d") >> 1069 G4cout << "Error, distribution must be iso, cos, planar, beam1d, beam2d or user" << G4endl; >> 1070 else >> 1071 AngDistType = atype; >> 1072 if (AngDistType == "cos") MaxTheta = pi/2. ; >> 1073 if (AngDistType == "user") { >> 1074 UDefThetaH = IPDFThetaH = ZeroPhysVector ; >> 1075 IPDFThetaExist = false ; >> 1076 UDefPhiH = IPDFPhiH = ZeroPhysVector ; >> 1077 IPDFPhiExist = false ; 125 } 1078 } >> 1079 } >> 1080 >> 1081 void G4GeneralParticleSource::DefineAngRefAxes(G4String refname, G4ThreeVector ref) >> 1082 { >> 1083 >> 1084 if(refname == "angref1") >> 1085 AngRef1 = ref.unit(); // x' >> 1086 else if(refname == "angref2") >> 1087 AngRef2 = ref.unit(); // vector in x'y' plane >> 1088 >> 1089 // User defines x' (AngRef1) and a vector in the x'y' >> 1090 // plane (AngRef2). Then, AngRef1 x AngRef2 = AngRef3 >> 1091 // the z' vector. Then, AngRef3 x AngRef1 = AngRef2 >> 1092 // which will now be y'. >> 1093 >> 1094 AngRef3 = AngRef1.cross(AngRef2); // z' >> 1095 AngRef2 = AngRef3.cross(AngRef1); // y' >> 1096 UserAngRef = true ; >> 1097 if(verbosityLevel == 2) >> 1098 { >> 1099 G4cout << "Angular distribution rotation axes " << AngRef1 << " " << AngRef2 << " " << AngRef3 << G4endl; >> 1100 } >> 1101 } 126 1102 127 // Set back previous source << 1103 void G4GeneralParticleSource::SetMinTheta(G4double mint) 128 GPSData->GetCurrentSource(currentIdx); << 1104 { >> 1105 MinTheta = mint; 129 } 1106 } 130 1107 131 void G4GeneralParticleSource::SetCurrentSource << 1108 void G4GeneralParticleSource::SetMinPhi(G4double minp) 132 { 1109 { 133 G4int id = aV; << 1110 MinPhi = minp; 134 if ( id < GPSData->GetIntensityVectorSize() << 1111 } 135 { << 1112 136 // currentSourceIdx = aV; << 1113 void G4GeneralParticleSource::SetMaxTheta(G4double maxt) 137 // currentSource = GPSData->GetCurrentSour << 1114 { 138 theMessenger->SetParticleGun(GPSData->GetC << 1115 MaxTheta = maxt; >> 1116 } >> 1117 >> 1118 void G4GeneralParticleSource::SetMaxPhi(G4double maxp) >> 1119 { >> 1120 MaxPhi = maxp; >> 1121 } >> 1122 >> 1123 void G4GeneralParticleSource::SetBeamSigmaInAngR(G4double r) >> 1124 { >> 1125 DR = r; >> 1126 } >> 1127 >> 1128 void G4GeneralParticleSource::SetBeamSigmaInAngX(G4double r) >> 1129 { >> 1130 DX = r; >> 1131 } >> 1132 >> 1133 void G4GeneralParticleSource::SetBeamSigmaInAngY(G4double r) >> 1134 { >> 1135 DY = r; >> 1136 } >> 1137 >> 1138 void G4GeneralParticleSource::UserDefAngTheta(G4ThreeVector input) >> 1139 { >> 1140 if(UserDistType == "NULL") UserDistType = "theta"; >> 1141 if(UserDistType == "phi") UserDistType = "both"; >> 1142 G4double thi, val; >> 1143 thi = input.x(); >> 1144 val = input.y(); >> 1145 if(verbosityLevel >= 1) >> 1146 G4cout << "In UserDefAngTheta" << G4endl; >> 1147 UDefThetaH.InsertValues(thi, val); >> 1148 } >> 1149 >> 1150 void G4GeneralParticleSource::UserDefAngPhi(G4ThreeVector input) >> 1151 { >> 1152 if(UserDistType == "NULL") UserDistType = "phi"; >> 1153 if(UserDistType == "theta") UserDistType = "both"; >> 1154 G4double phhi, val; >> 1155 phhi = input.x(); >> 1156 val = input.y(); >> 1157 if(verbosityLevel >= 1) >> 1158 G4cout << "In UserDefAngPhi" << G4endl; >> 1159 UDefPhiH.InsertValues(phhi, val); >> 1160 } >> 1161 >> 1162 void G4GeneralParticleSource::SetUserWRTSurface(G4bool wrtSurf) >> 1163 { >> 1164 // if UserWRTSurface = true then the user wants momenta with respect >> 1165 // to the surface normals. >> 1166 // When doing this theta has to be 0-90 only otherwise there will be >> 1167 // errors, which currently are flagged anywhere. >> 1168 UserWRTSurface = wrtSurf; >> 1169 } >> 1170 >> 1171 void G4GeneralParticleSource::SetUseUserAngAxis(G4bool userang) >> 1172 { >> 1173 // if UserAngRef = true the angular distribution is defined wrt >> 1174 // the user defined co-ordinates >> 1175 UserAngRef = userang; >> 1176 } >> 1177 >> 1178 void G4GeneralParticleSource::GenerateBeamFlux() >> 1179 { >> 1180 // generates isotropic flux. >> 1181 // No vectors are needed. >> 1182 G4double theta, phi; >> 1183 G4double px, py, pz; >> 1184 if (AngDistType == "beam1d") >> 1185 { >> 1186 theta = G4RandGauss::shoot(0.0,DR); >> 1187 phi = twopi * G4UniformRand(); >> 1188 } >> 1189 else >> 1190 { >> 1191 px = G4RandGauss::shoot(0.0,DX); >> 1192 py = G4RandGauss::shoot(0.0,DY); >> 1193 theta = sqrt (px*px + py*py); >> 1194 if (theta != 0.) { >> 1195 phi = acos(px/theta); >> 1196 if ( py < 0.) phi = -phi; >> 1197 } >> 1198 else >> 1199 { >> 1200 phi = 0.0; >> 1201 } >> 1202 } >> 1203 px = -sin(theta) * cos(phi); >> 1204 py = -sin(theta) * sin(phi); >> 1205 pz = -cos(theta); >> 1206 // Apply Angular Rotation Matrix >> 1207 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1208 G4double finx, finy, finz; >> 1209 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1210 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1211 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1212 G4double ResMag = sqrt((finx*finx) + (finy*finy) + (finz*finz)); >> 1213 finx = finx/ResMag; >> 1214 finy = finy/ResMag; >> 1215 finz = finz/ResMag; >> 1216 >> 1217 particle_momentum_direction.setX(finx); >> 1218 particle_momentum_direction.setY(finy); >> 1219 particle_momentum_direction.setZ(finz); >> 1220 >> 1221 // particle_momentum_direction now holds unit momentum vector. >> 1222 if(verbosityLevel >= 1) >> 1223 G4cout << "Generating beam vector: " << particle_momentum_direction << G4endl; >> 1224 } >> 1225 >> 1226 void G4GeneralParticleSource::GenerateIsotropicFlux() >> 1227 { >> 1228 // generates isotropic flux. >> 1229 // No vectors are needed. >> 1230 G4double rndm, rndm2; >> 1231 G4double px, py, pz; >> 1232 >> 1233 // generate rand nos. but make sure in theta/phi limits >> 1234 /* Colin's old stuff >> 1235 Theta = 10.; // This is well above the pi upper limit. >> 1236 Phi = 10.; // Still well above 2pi upper limit. >> 1237 >> 1238 while(Theta > MaxTheta || Theta < MinTheta) >> 1239 { >> 1240 rndm = GenRandTheta(); >> 1241 Theta = acos(1. - 2.*rndm); >> 1242 } >> 1243 >> 1244 while(Phi > MaxPhi || Phi < MinPhi) >> 1245 { >> 1246 rndm2 = GenRandPhi(); >> 1247 Phi = twopi * rndm2; >> 1248 } >> 1249 >> 1250 px = sin(Theta) * cos(Phi); >> 1251 py = sin(Theta) * sin(Phi); >> 1252 pz = cos(Theta); >> 1253 >> 1254 */ >> 1255 // More efficient implementation by F. Lei >> 1256 // >> 1257 G4double sintheta, sinphi,costheta,cosphi; >> 1258 rndm = GenRandTheta(); >> 1259 costheta = cos(MinTheta) - rndm * (cos(MinTheta) - cos(MaxTheta)); >> 1260 sintheta = sqrt(1. - costheta*costheta); >> 1261 >> 1262 rndm2 = GenRandPhi(); >> 1263 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; >> 1264 sinphi = sin(Phi); >> 1265 cosphi = cos(Phi); >> 1266 >> 1267 px = -sintheta * cosphi; >> 1268 py = -sintheta * sinphi; >> 1269 pz = -costheta; >> 1270 >> 1271 // end of F.Lei implementation >> 1272 >> 1273 // px = -px; >> 1274 //py = -py; >> 1275 //pz = -pz; >> 1276 >> 1277 // pmag = sqrt((px*px) + (py*py) + (pz*pz)); >> 1278 >> 1279 //px = px/pmag; >> 1280 //py = py/pmag; >> 1281 //pz = pz/pmag; >> 1282 >> 1283 /* this implmentation is commented out by f.lei >> 1284 >> 1285 //Need to rotate the particle_momentum_direction round such that the >> 1286 //start position is at the north pole, so that all particles go into >> 1287 // the centre of the geomtries. >> 1288 // particle_position stores the particle's position. >> 1289 // use particle position to calculate theta and phi - posthe and posphi >> 1290 >> 1291 G4double posthe, posphi; >> 1292 // G4cout << "particle_position " << particle_position << G4endl; >> 1293 G4double tx, ty, tz, tt; >> 1294 tx = particle_position.x(); >> 1295 ty = particle_position.y(); >> 1296 tz = particle_position.z(); >> 1297 tt = sqrt(tx*tx + ty*ty + tz*tz); >> 1298 tx = tx/tt; >> 1299 ty = ty/tt; >> 1300 tz = tz/tt; >> 1301 // G4cout << "unit position " << tx << " " << ty << " " << tz << G4endl; >> 1302 if (tt != 0.) { >> 1303 posthe = acos(tz); 139 } 1304 } 140 else << 1305 else { 141 { << 1306 posthe = 0.; 142 G4ExceptionDescription msg; << 1307 } 143 msg << "Trying to set source to index " << << 1308 if (posthe != 0. && posthe != pi){ 144 << GPSData->GetIntensityVectorSize() < << 1309 posphi = acos(tx/sin(posthe)); 145 G4Exception("G4GeneralParticleSoruce::SetC << 1310 } 146 FatalException, msg); << 1311 else{ >> 1312 posphi = 0.0; >> 1313 } >> 1314 // G4cout << "Posthe and posphi " << posthe << " " << posphi << G4endl; >> 1315 G4double finx, finy, finz; >> 1316 finx = (px*cos(posthe)*cos(posphi)) - (py*sin(posphi)) + (pz*sin(posthe)*cos(posphi)); >> 1317 finy = (px*cos(posthe)*sin(posphi)) + (py*cos(posphi)) + (pz*sin(posthe)*sin(posphi)); >> 1318 finz = (-px*sin(posthe)) + (pz*cos(posthe)); >> 1319 // G4cout << "finx... " << finx << " " << finy << " " << finz << G4endl; >> 1320 */ >> 1321 >> 1322 // New implementation >> 1323 // for volume and ponit source use mother or user defined co-ordinates >> 1324 // for plane and surface source user surface-normal or userdefined co-ordinates >> 1325 // >> 1326 G4double finx, finy, finz; >> 1327 if (SourcePosType == "Point" || SourcePosType == "Volume") { >> 1328 if (UserAngRef){ >> 1329 // Apply Rotation Matrix >> 1330 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1331 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1332 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1333 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1334 } else { >> 1335 finx = px; >> 1336 finy = py; >> 1337 finz = pz; >> 1338 } >> 1339 } else { // for plane and surface source >> 1340 if (UserAngRef){ >> 1341 // Apply Rotation Matrix >> 1342 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1343 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1344 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1345 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1346 } else { >> 1347 finx = (px*SideRefVec1.x()) + (py*SideRefVec2.x()) + (pz*SideRefVec3.x()); >> 1348 finy = (px*SideRefVec1.y()) + (py*SideRefVec2.y()) + (pz*SideRefVec3.y()); >> 1349 finz = (px*SideRefVec1.z()) + (py*SideRefVec2.z()) + (pz*SideRefVec3.z()); >> 1350 } 147 } 1351 } >> 1352 G4double ResMag = sqrt((finx*finx) + (finy*finy) + (finz*finz)); >> 1353 finx = finx/ResMag; >> 1354 finy = finy/ResMag; >> 1355 finz = finz/ResMag; >> 1356 >> 1357 particle_momentum_direction.setX(finx); >> 1358 particle_momentum_direction.setY(finy); >> 1359 particle_momentum_direction.setZ(finz); >> 1360 >> 1361 // particle_momentum_direction now holds unit momentum vector. >> 1362 if(verbosityLevel >= 1) >> 1363 G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl; 148 } 1364 } 149 1365 150 void G4GeneralParticleSource::SetCurrentSource << 1366 void G4GeneralParticleSource::GenerateCosineLawFlux() 151 { 1367 { 152 GPSData->Lock(); << 1368 // Method to generate flux distributed with a cosine law 153 GPSData->SetCurrentSourceIntensity(aV); << 1369 // such as is in TIMM. 154 GPSData->Unlock(); << 1370 G4double px, py, pz; 155 normalised = GPSData->Normalised(); << 1371 G4double rndm, rndm2; >> 1372 // G4double resultx, resulty, resultz; >> 1373 >> 1374 /* Colin's old implementation which is very slow >> 1375 >> 1376 Theta = 10.; // Well above MaxTheta >> 1377 Phi = 10.; // Well above MaxPhi >> 1378 >> 1379 while(Theta > MaxTheta || Theta < MinTheta) >> 1380 { >> 1381 rndm = GenRandTheta(); >> 1382 Theta = asin(sqrt(MaxTheta*rndm)); >> 1383 } >> 1384 >> 1385 while(Phi > MaxPhi || Phi < MinPhi) >> 1386 { >> 1387 rndm2 = GenRandPhi(); >> 1388 Phi = twopi * rndm2; >> 1389 } >> 1390 >> 1391 px = sin(Theta) * cos(Phi); >> 1392 py = sin(Theta) * sin(Phi); >> 1393 pz = cos(Theta); >> 1394 */ >> 1395 // More efficient implementation by F. Lei >> 1396 // >> 1397 G4double sintheta, sinphi,costheta,cosphi; >> 1398 rndm = GenRandTheta(); >> 1399 sintheta = sqrt( rndm * (sin(MaxTheta)*sin(MaxTheta) - sin(MinTheta)*sin(MinTheta) ) >> 1400 +sin(MinTheta)*sin(MinTheta) ); >> 1401 costheta = sqrt(1. -sintheta*sintheta); >> 1402 >> 1403 rndm2 = GenRandPhi(); >> 1404 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; >> 1405 sinphi = sin(Phi); >> 1406 cosphi = cos(Phi); >> 1407 >> 1408 px = -sintheta * cosphi; >> 1409 py = -sintheta * sinphi; >> 1410 pz = -costheta; >> 1411 >> 1412 // end of F.Lei implementation >> 1413 >> 1414 // px = -px; >> 1415 //py = -py; >> 1416 //pz = -pz; >> 1417 //pmag = sqrt((px*px) + (py*py) + (pz*pz)); >> 1418 //G4double pxh = px/pmag; >> 1419 //G4double pyh = py/pmag; >> 1420 //G4double pzh = pz/pmag; >> 1421 >> 1422 /* commented out by F. Lei >> 1423 if(verbosityLevel == 2) >> 1424 { >> 1425 G4cout <<"SideRefVecs " <<SideRefVec1<<SideRefVec2<<SideRefVec3<<G4endl; >> 1426 G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl; >> 1427 } >> 1428 resultx = (pxh*SideRefVec1.x()) + (pyh*SideRefVec2.x()) + >> 1429 (pzh*SideRefVec3.x()); >> 1430 >> 1431 resulty = (pxh*SideRefVec1.y()) + (pyh*SideRefVec2.y()) + >> 1432 (pzh*SideRefVec3.y()); >> 1433 >> 1434 resultz = (pxh*SideRefVec1.z()) + (pyh*SideRefVec2.z()) + >> 1435 (pzh*SideRefVec3.z()); >> 1436 >> 1437 G4double ResMag = sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz)); >> 1438 resultx = resultx/ResMag; >> 1439 resulty = resulty/ResMag; >> 1440 resultz = resultz/ResMag; >> 1441 >> 1442 particle_momentum_direction.setX(resultx); >> 1443 particle_momentum_direction.setY(resulty); >> 1444 particle_momentum_direction.setZ(resultz); >> 1445 >> 1446 */ >> 1447 // New implementation >> 1448 // for volume and ponit source use mother or user defined co-ordinates >> 1449 // for plane and surface source user surface-normal or userdefined co-ordinates >> 1450 // >> 1451 G4double finx, finy, finz; >> 1452 if (SourcePosType == "Point" || SourcePosType == "Volume") { >> 1453 if (UserAngRef){ >> 1454 // Apply Rotation Matrix >> 1455 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1456 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1457 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1458 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1459 } else { >> 1460 finx = px; >> 1461 finy = py; >> 1462 finz = pz; >> 1463 } >> 1464 } else { // for plane and surface source >> 1465 if (UserAngRef){ >> 1466 // Apply Rotation Matrix >> 1467 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1468 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1469 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1470 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1471 } else { >> 1472 finx = (px*SideRefVec1.x()) + (py*SideRefVec2.x()) + (pz*SideRefVec3.x()); >> 1473 finy = (px*SideRefVec1.y()) + (py*SideRefVec2.y()) + (pz*SideRefVec3.y()); >> 1474 finz = (px*SideRefVec1.z()) + (py*SideRefVec2.z()) + (pz*SideRefVec3.z()); >> 1475 } >> 1476 } >> 1477 G4double ResMag = sqrt((finx*finx) + (finy*finy) + (finz*finz)); >> 1478 finx = finx/ResMag; >> 1479 finy = finy/ResMag; >> 1480 finz = finz/ResMag; >> 1481 >> 1482 particle_momentum_direction.setX(finx); >> 1483 particle_momentum_direction.setY(finy); >> 1484 particle_momentum_direction.setZ(finz); >> 1485 >> 1486 // particle_momentum_direction now contains unit momentum vector. >> 1487 if(verbosityLevel >= 1) >> 1488 { >> 1489 G4cout << "Resultant cosine-law unit momentum vector " << particle_momentum_direction << G4endl; >> 1490 } 156 } 1491 } 157 1492 158 void G4GeneralParticleSource::ClearAll() << 1493 void G4GeneralParticleSource::GeneratePlanarFlux() 159 { 1494 { 160 GPSData->ClearSources(); << 1495 // particle_momentum_direction now contains unit momentum vector. 161 normalised=GPSData->Normalised(); << 1496 // nothing need be done here as the m-directions have been set directly >> 1497 // under this option >> 1498 if(verbosityLevel >= 1) >> 1499 { >> 1500 G4cout << "Resultant Planar wave momentum vector " << particle_momentum_direction << G4endl; >> 1501 } 162 } 1502 } 163 1503 164 void G4GeneralParticleSource::DeleteaSource(G4 << 1504 void G4GeneralParticleSource::GenerateUserDefFlux() 165 { 1505 { 166 G4int id = aV; << 1506 G4double rndm, px, py, pz, pmag; 167 if ( id <= GPSData->GetIntensityVectorSize() << 1507 168 { << 1508 if(UserDistType == "NULL") 169 GPSData->DeleteASource(aV); << 1509 G4cout << "Error: UserDistType undefined" << G4endl; 170 normalised=GPSData->Normalised(); << 1510 else if(UserDistType == "theta") >> 1511 { >> 1512 Theta = GenerateUserDefTheta(); >> 1513 Phi = 10.; >> 1514 while(Phi > MaxPhi || Phi < MinPhi) >> 1515 { >> 1516 rndm = GenRandPhi(); >> 1517 Phi = twopi * rndm; >> 1518 } >> 1519 } >> 1520 else if(UserDistType == "phi") >> 1521 { >> 1522 Theta = 10.; >> 1523 while(Theta > MaxTheta || Theta < MinTheta) >> 1524 { >> 1525 rndm = GenRandTheta(); >> 1526 Theta = acos(1. - (2. * rndm)); >> 1527 } >> 1528 Phi = GenerateUserDefPhi(); >> 1529 } >> 1530 else if(UserDistType == "both") >> 1531 { >> 1532 Theta = GenerateUserDefTheta(); >> 1533 Phi = GenerateUserDefPhi(); >> 1534 } >> 1535 >> 1536 px = -sin(Theta) * cos(Phi); >> 1537 py = -sin(Theta) * sin(Phi); >> 1538 pz = -cos(Theta); >> 1539 >> 1540 pmag = sqrt((px*px) + (py*py) + (pz*pz)); >> 1541 >> 1542 if(UserWRTSurface == false) >> 1543 { >> 1544 G4double finx, finy, finz; >> 1545 if (UserAngRef) { >> 1546 // Apply Rotation Matrix >> 1547 // x * AngRef1, y * AngRef2 and z * AngRef3 >> 1548 finx = (px * AngRef1.x()) + (py * AngRef2.x()) + (pz * AngRef3.x()); >> 1549 finy = (px * AngRef1.y()) + (py * AngRef2.y()) + (pz * AngRef3.y()); >> 1550 finz = (px * AngRef1.z()) + (py * AngRef2.z()) + (pz * AngRef3.z()); >> 1551 } else { // use mother co-ordinates >> 1552 finx = px; >> 1553 finy = py; >> 1554 finz = pz; >> 1555 } >> 1556 G4double ResMag = sqrt((finx*finx) + (finy*finy) + (finz*finz)); >> 1557 finx = finx/ResMag; >> 1558 finy = finy/ResMag; >> 1559 finz = finz/ResMag; >> 1560 >> 1561 particle_momentum_direction.setX(finx); >> 1562 particle_momentum_direction.setY(finy); >> 1563 particle_momentum_direction.setZ(finz); >> 1564 } >> 1565 else if(UserWRTSurface == true) >> 1566 { >> 1567 G4double pxh = px/pmag; >> 1568 G4double pyh = py/pmag; >> 1569 G4double pzh = pz/pmag; >> 1570 >> 1571 if(verbosityLevel == 2) >> 1572 { >> 1573 G4cout <<"SideRefVecs " <<SideRefVec1<<SideRefVec2<<SideRefVec3<<G4endl; >> 1574 G4cout <<"Raw Unit vector "<<pxh<<","<<pyh<<","<<pzh<<G4endl; >> 1575 } >> 1576 G4double resultx = (pxh*SideRefVec1.x()) + (pyh*SideRefVec2.x()) + >> 1577 (pzh*SideRefVec3.x()); >> 1578 >> 1579 G4double resulty = (pxh*SideRefVec1.y()) + (pyh*SideRefVec2.y()) + >> 1580 (pzh*SideRefVec3.y()); >> 1581 >> 1582 G4double resultz = (pxh*SideRefVec1.z()) + (pyh*SideRefVec2.z()) + >> 1583 (pzh*SideRefVec3.z()); >> 1584 >> 1585 G4double ResMag = sqrt((resultx*resultx) + (resulty*resulty) + (resultz*resultz)); >> 1586 resultx = resultx/ResMag; >> 1587 resulty = resulty/ResMag; >> 1588 resultz = resultz/ResMag; >> 1589 >> 1590 particle_momentum_direction.setX(resultx); >> 1591 particle_momentum_direction.setY(resulty); >> 1592 particle_momentum_direction.setZ(resultz); >> 1593 } >> 1594 >> 1595 // particle_momentum_direction now contains unit momentum vector. >> 1596 if(verbosityLevel >= 1) >> 1597 { >> 1598 G4cout << "Final User Defined momentum vector " << particle_momentum_direction << G4endl; >> 1599 } >> 1600 } >> 1601 >> 1602 G4double G4GeneralParticleSource::GenerateUserDefTheta() >> 1603 { >> 1604 // Create cumulative histogram if not already done so. Then use RandFlat >> 1605 //::shoot to generate the output Theta value. >> 1606 if(UserDistType == "NULL" || UserDistType == "phi") >> 1607 { >> 1608 // No user defined theta distribution >> 1609 G4cout << "Error ***********************" << G4endl; >> 1610 G4cout << "UserDistType = " << UserDistType << G4endl; >> 1611 return (0.); >> 1612 } >> 1613 else >> 1614 { >> 1615 // UserDistType = theta or both and so a theta distribution >> 1616 // is defined. This should be integrated if not already done. >> 1617 if(IPDFThetaExist == false) >> 1618 { >> 1619 // IPDF has not been created, so create it >> 1620 G4double bins[256],vals[256], sum; >> 1621 G4int ii; >> 1622 G4int maxbin = G4int(UDefThetaH.GetVectorLength()); >> 1623 bins[0] = UDefThetaH.GetLowEdgeEnergy(size_t(0)); >> 1624 vals[0] = UDefThetaH(size_t(0)); >> 1625 sum = vals[0]; >> 1626 for(ii=1;ii<maxbin;ii++) >> 1627 { >> 1628 bins[ii] = UDefThetaH.GetLowEdgeEnergy(size_t(ii)); >> 1629 vals[ii] = UDefThetaH(size_t(ii)) + vals[ii-1]; >> 1630 sum = sum + UDefThetaH(size_t(ii)); >> 1631 } >> 1632 >> 1633 for(ii=0;ii<maxbin;ii++) >> 1634 { >> 1635 vals[ii] = vals[ii]/sum; >> 1636 IPDFThetaH.InsertValues(bins[ii], vals[ii]); >> 1637 } >> 1638 // Make IPDFThetaExist = true >> 1639 IPDFThetaExist = true; >> 1640 } >> 1641 // IPDF has been create so carry on >> 1642 G4double rndm = G4UniformRand(); >> 1643 return(IPDFThetaH.GetEnergy(rndm)); >> 1644 } >> 1645 } >> 1646 >> 1647 G4double G4GeneralParticleSource::GenerateUserDefPhi() >> 1648 { >> 1649 // Create cumulative histogram if not already done so. Then use RandFlat >> 1650 //::shoot to generate the output Theta value. >> 1651 >> 1652 if(UserDistType == "NULL" || UserDistType == "theta") >> 1653 { >> 1654 // No user defined phi distribution >> 1655 G4cout << "Error ***********************" << G4endl; >> 1656 G4cout << "UserDistType = " << UserDistType << G4endl; >> 1657 return(0.); >> 1658 } >> 1659 else >> 1660 { >> 1661 // UserDistType = phi or both and so a phi distribution >> 1662 // is defined. This should be integrated if not already done. >> 1663 if(IPDFPhiExist == false) >> 1664 { >> 1665 // IPDF has not been created, so create it >> 1666 G4double bins[256],vals[256], sum; >> 1667 G4int ii; >> 1668 G4int maxbin = G4int(UDefPhiH.GetVectorLength()); >> 1669 bins[0] = UDefPhiH.GetLowEdgeEnergy(size_t(0)); >> 1670 vals[0] = UDefPhiH(size_t(0)); >> 1671 sum = vals[0]; >> 1672 for(ii=1;ii<maxbin;ii++) >> 1673 { >> 1674 bins[ii] = UDefPhiH.GetLowEdgeEnergy(size_t(ii)); >> 1675 vals[ii] = UDefPhiH(size_t(ii)) + vals[ii-1]; >> 1676 sum = sum + UDefPhiH(size_t(ii)); >> 1677 } >> 1678 >> 1679 for(ii=0;ii<maxbin;ii++) >> 1680 { >> 1681 vals[ii] = vals[ii]/sum; >> 1682 IPDFPhiH.InsertValues(bins[ii], vals[ii]); >> 1683 } >> 1684 // Make IPDFPhiExist = true >> 1685 IPDFPhiExist = true; >> 1686 } >> 1687 // IPDF has been create so carry on >> 1688 G4double rndm = G4UniformRand(); >> 1689 return(IPDFPhiH.GetEnergy(rndm)); >> 1690 } >> 1691 } >> 1692 >> 1693 // Now follows Energy distribution Methods >> 1694 >> 1695 void G4GeneralParticleSource::SetEnergyDisType(G4String DisType) >> 1696 { >> 1697 EnergyDisType = DisType; >> 1698 if (EnergyDisType == "User"){ >> 1699 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; >> 1700 IPDFEnergyExist = false ; >> 1701 } else if ( EnergyDisType == "Arb"){ >> 1702 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ; >> 1703 IPDFArbExist = false; >> 1704 } else if (EnergyDisType == "Epn"){ >> 1705 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; >> 1706 IPDFEnergyExist = false ; >> 1707 EpnEnergyH = ZeroPhysVector ; 171 } 1708 } >> 1709 } >> 1710 >> 1711 void G4GeneralParticleSource::SetEmin(G4double emi) >> 1712 { >> 1713 Emin = emi; >> 1714 } >> 1715 >> 1716 void G4GeneralParticleSource::SetEmax(G4double ema) >> 1717 { >> 1718 Emax = ema; >> 1719 } >> 1720 >> 1721 void G4GeneralParticleSource::SetMonoEnergy(G4double menergy) >> 1722 { >> 1723 MonoEnergy = menergy; >> 1724 Emin = menergy; >> 1725 Emax = menergy; >> 1726 } >> 1727 >> 1728 void G4GeneralParticleSource::SetBeamSigmaInE(G4double e) >> 1729 { >> 1730 SE = e; >> 1731 } >> 1732 void G4GeneralParticleSource::SetAlpha(G4double alp) >> 1733 { >> 1734 alpha = alp; >> 1735 } >> 1736 >> 1737 void G4GeneralParticleSource::SetTemp(G4double tem) >> 1738 { >> 1739 Temp = tem; >> 1740 } >> 1741 >> 1742 void G4GeneralParticleSource::SetEzero(G4double eze) >> 1743 { >> 1744 Ezero = eze; >> 1745 } >> 1746 >> 1747 void G4GeneralParticleSource::SetGradient(G4double gr) >> 1748 { >> 1749 grad = gr; >> 1750 } >> 1751 >> 1752 void G4GeneralParticleSource::SetInterCept(G4double c) >> 1753 { >> 1754 cept = c; >> 1755 } >> 1756 >> 1757 void G4GeneralParticleSource::UserEnergyHisto(G4ThreeVector input) >> 1758 { >> 1759 G4double ehi, val; >> 1760 ehi = input.x(); >> 1761 val = input.y(); >> 1762 if(verbosityLevel == 2) >> 1763 G4cout << " " << ehi << " " << val << G4endl; >> 1764 if(verbosityLevel == 2) >> 1765 G4cout << "In UserEnergyHisto" << G4endl; >> 1766 UDefEnergyH.InsertValues(ehi, val); >> 1767 Emax = ehi; >> 1768 } >> 1769 >> 1770 void G4GeneralParticleSource::ArbEnergyHisto(G4ThreeVector input) >> 1771 { >> 1772 G4double ehi, val; >> 1773 ehi = input.x(); >> 1774 val = input.y(); >> 1775 if(verbosityLevel == 2) >> 1776 G4cout << "In ArbEnergyHisto" << G4endl; >> 1777 ArbEnergyH.InsertValues(ehi, val); >> 1778 } >> 1779 >> 1780 void G4GeneralParticleSource::EpnEnergyHisto(G4ThreeVector input) >> 1781 { >> 1782 G4double ehi, val; >> 1783 ehi = input.x(); >> 1784 val = input.y(); >> 1785 if(verbosityLevel == 2) >> 1786 G4cout << "In EpnEnergyHisto" << G4endl; >> 1787 EpnEnergyH.InsertValues(ehi, val); >> 1788 Emax = ehi; >> 1789 Epnflag = true; >> 1790 } >> 1791 >> 1792 void G4GeneralParticleSource::Calculate() >> 1793 { >> 1794 if(EnergyDisType == "Cdg") >> 1795 CalculateCdgSpectrum(); >> 1796 else if(EnergyDisType == "Bbody") >> 1797 CalculateBbodySpectrum(); >> 1798 } >> 1799 >> 1800 void G4GeneralParticleSource::CalculateCdgSpectrum() >> 1801 { >> 1802 // This uses the spectrum from The INTEGRAL Mass Model (TIMM) >> 1803 // to generate a Cosmic Diffuse X/gamma ray spectrum. >> 1804 G4double pfact[2] = {8.5, 112}; >> 1805 G4double spind[2] = {1.4, 2.3}; >> 1806 G4double ene_line[3] = {1.*keV, 18.*keV, 1E6*keV}; >> 1807 G4int n_par; >> 1808 >> 1809 ene_line[0] = Emin; >> 1810 if(Emin < 18*keV) >> 1811 { >> 1812 n_par = 2; >> 1813 ene_line[2] = Emax; >> 1814 if(Emax < 18*keV) >> 1815 { >> 1816 n_par = 1; >> 1817 ene_line[1] = Emax; >> 1818 } >> 1819 } 172 else 1820 else 173 { << 1821 { 174 G4cout << " source index is invalid " << G << 1822 n_par = 1; 175 G4cout << " it shall be <= " << 1823 pfact[0] = 112.; 176 << GPSData->GetIntensityVectorSize( << 1824 spind[0] = 2.3; >> 1825 ene_line[1] = Emax; >> 1826 } >> 1827 >> 1828 // Create a cumulative histogram. >> 1829 CDGhist[0] = 0.; >> 1830 G4double omalpha; >> 1831 G4int i = 0; >> 1832 >> 1833 while(i < n_par) >> 1834 { >> 1835 omalpha = 1. - spind[i]; >> 1836 CDGhist[i+1] = CDGhist[i] + (pfact[i]/omalpha)* >> 1837 (pow(ene_line[i+1],omalpha)-pow(ene_line[i],omalpha)); >> 1838 i++; >> 1839 } >> 1840 >> 1841 // Normalise histo and divide by 1000 to make MeV. >> 1842 i = 0; >> 1843 while(i < n_par) >> 1844 { >> 1845 CDGhist[i+1] = CDGhist[i+1]/CDGhist[n_par]; >> 1846 // G4cout << CDGhist[i] << CDGhist[n_par] << G4endl; >> 1847 i++; >> 1848 } >> 1849 } >> 1850 >> 1851 void G4GeneralParticleSource::CalculateBbodySpectrum() >> 1852 { >> 1853 // create bbody spectrum >> 1854 // Proved very hard to integrate indefinitely, so different >> 1855 // method. User inputs emin, emax and T. These are used to >> 1856 // create a 10,000 bin histogram. >> 1857 // Use photon density spectrum = 2 nu**2/c**2 * (exp(h nu/kT)-1) >> 1858 // = 2 E**2/h**2c**2 times the exponential >> 1859 G4double erange = Emax - Emin; >> 1860 G4double steps = erange/10000.; >> 1861 G4double Bbody_y[10000]; >> 1862 G4double k = 8.6181e-11; //Boltzmann const in MeV/K >> 1863 G4double h = 4.1362e-21; // Plancks const in MeV s >> 1864 G4double c = 3e8; // Speed of light >> 1865 G4double h2 = h*h; >> 1866 G4double c2 = c*c; >> 1867 G4int count = 0; >> 1868 G4double sum = 0.; >> 1869 BBHist[0] = 0.; >> 1870 while(count < 10000) >> 1871 { >> 1872 Bbody_x[count] = Emin + G4double(count*steps); >> 1873 Bbody_y[count] = (2.*pow(Bbody_x[count],2.))/ >> 1874 (h2*c2*(exp(Bbody_x[count]/(k*Temp)) - 1.)); >> 1875 sum = sum + Bbody_y[count]; >> 1876 BBHist[count+1] = BBHist[count] + Bbody_y[count]; >> 1877 count++; >> 1878 } >> 1879 >> 1880 Bbody_x[10000] = Emax; >> 1881 // Normalise cumulative histo. >> 1882 count = 0; >> 1883 while(count<10001) >> 1884 { >> 1885 BBHist[count] = BBHist[count]/sum; >> 1886 count++; >> 1887 } >> 1888 } >> 1889 >> 1890 void G4GeneralParticleSource::InputEnergySpectra(G4bool value) >> 1891 { >> 1892 // Allows user to specifiy spectrum is momentum >> 1893 EnergySpec = value; // false if momentum >> 1894 if(verbosityLevel == 2) >> 1895 G4cout << "EnergySpec has value " << EnergySpec << G4endl; >> 1896 } >> 1897 >> 1898 void G4GeneralParticleSource::InputDifferentialSpectra(G4bool value) >> 1899 { >> 1900 // Allows user to specify integral or differential spectra >> 1901 DiffSpec = value; // true = differential, false = integral >> 1902 if(verbosityLevel == 2) >> 1903 G4cout << "Diffspec has value " << DiffSpec << G4endl; >> 1904 } >> 1905 >> 1906 void G4GeneralParticleSource::ArbInterpolate(G4String IType) >> 1907 { >> 1908 if(EnergyDisType != "Arb") >> 1909 G4cout << "Error: this is for arbitrary distributions" << G4endl; >> 1910 IntType = IType; >> 1911 >> 1912 >> 1913 // Calcuate Emin and Emax, mainly for use in debugging >> 1914 G4int len = G4int(ArbEnergyH.GetVectorLength()); >> 1915 ArbEmax = ArbEnergyH.GetLowEdgeEnergy(size_t(len-1)); >> 1916 ArbEmin = ArbEnergyH.GetLowEdgeEnergy(size_t(0)); >> 1917 G4cout << "ArbEmin, ArbEmax " << ArbEmin << " " << ArbEmax << G4endl; >> 1918 >> 1919 // Now interpolate points >> 1920 if(IntType == "Lin") >> 1921 LinearInterpolation(); >> 1922 if(IntType == "Log") >> 1923 LogInterpolation(); >> 1924 if(IntType == "Exp") >> 1925 ExpInterpolation(); >> 1926 if(IntType == "Spline") >> 1927 SplineInterpolation(); >> 1928 } >> 1929 >> 1930 void G4GeneralParticleSource::LinearInterpolation() >> 1931 { >> 1932 // Method to do linear interpolation on the Arb points >> 1933 // Calculate equation of each line segment, max 256. >> 1934 // Calculate Area under each segment >> 1935 // Create a cumulative array which is then normalised Arb_Cum_Area >> 1936 G4double Area_seg[256]; // Stores area under each segment >> 1937 G4double sum = 0., Arb_x[256], Arb_y[256], Arb_Cum_Area[256]; >> 1938 G4int i, count; >> 1939 G4int maxi = ArbEnergyH.GetVectorLength(); >> 1940 for(i=0;i<maxi;i++) >> 1941 { >> 1942 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); >> 1943 Arb_y[i] = ArbEnergyH(size_t(i)); >> 1944 } >> 1945 // Points are now in x,y arrays. If the spectrum is integral it has to be >> 1946 // made differential and if momentum it has to be made energy. >> 1947 if(EnergySpec == false) >> 1948 { >> 1949 // change currently stored values (emin etc) which are actually momenta >> 1950 // to energies. >> 1951 if(particle_definition == NULL) >> 1952 G4cout << "Error: particle not defined" << G4endl; >> 1953 else >> 1954 { >> 1955 // Apply Energy**2 = p**2c**2 + m0**2c**4 >> 1956 // p should be entered as E/c i.e. without the division by c >> 1957 // being done - energy equivalent. >> 1958 G4double mass = particle_definition->GetPDGMass(); >> 1959 >> 1960 // multiply the function (Arb_y) up by the bin width >> 1961 // to make the function counts/s (i.e. get rid of momentum >> 1962 // dependence). >> 1963 for(count=0;count<maxi;count++) >> 1964 { >> 1965 Arb_y[count] = Arb_y[count] * (Arb_x[count+1] - Arb_x[count]); >> 1966 } >> 1967 // Change Arb_x to energy, plus divide by energy bin width >> 1968 // to make Arb_y counts/s/energy >> 1969 for(count=0;count<maxi;count++) >> 1970 { >> 1971 Arb_x[count] = sqrt((Arb_x[count]*Arb_x[count]) >> 1972 + (mass*mass)); >> 1973 } >> 1974 for(count=0;count<maxi;count++) >> 1975 { >> 1976 Arb_y[count] = Arb_y[count]/(Arb_x[count+1] - Arb_x[count]); >> 1977 } >> 1978 } >> 1979 } >> 1980 if(DiffSpec == false) >> 1981 { >> 1982 // Converts integral point-wise spectra to Differential >> 1983 for( count=1;count<=maxi;count++) >> 1984 { >> 1985 Arb_y[count] = Arb_y[count] - Arb_y[count-1]; >> 1986 } >> 1987 } >> 1988 >> 1989 i=1; >> 1990 Arb_grad[0] = 0.; >> 1991 Arb_cept[0] = 0.; >> 1992 Area_seg[0] = 0.; >> 1993 Arb_Cum_Area[0] = 0.; >> 1994 while(i < maxi) >> 1995 { >> 1996 // calc gradient and intercept for each segment >> 1997 Arb_grad[i] = (Arb_y[i] - Arb_y[i-1]) / (Arb_x[i] - Arb_x[i-1]); >> 1998 if(verbosityLevel == 2) >> 1999 G4cout << Arb_grad[i] << G4endl; >> 2000 if(Arb_grad[i] > 0.) >> 2001 { >> 2002 if(verbosityLevel == 2) >> 2003 G4cout << "Arb_grad is positive" << G4endl; >> 2004 Arb_cept[i] = Arb_y[i] - (Arb_grad[i] * Arb_x[i]); >> 2005 } >> 2006 else if(Arb_grad[i] < 0.) >> 2007 { >> 2008 if(verbosityLevel == 2) >> 2009 G4cout << "Arb_grad is negative" << G4endl; >> 2010 Arb_cept[i] = Arb_y[i] + (-Arb_grad[i] * Arb_x[i]); >> 2011 } >> 2012 else >> 2013 { >> 2014 if(verbosityLevel == 2) >> 2015 G4cout << "Arb_grad is 0." << G4endl; >> 2016 Arb_cept[i] = Arb_y[i]; >> 2017 } >> 2018 >> 2019 Area_seg[i] = ((Arb_grad[i]/2)*(Arb_x[i]*Arb_x[i] - Arb_x[i-1]*Arb_x[i-1]) + Arb_cept[i]*(Arb_x[i] - Arb_x[i-1])); >> 2020 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 2021 sum = sum + Area_seg[i]; >> 2022 if(verbosityLevel == 2) >> 2023 G4cout << Arb_x[i] << Arb_y[i] << Area_seg[i] << sum << Arb_grad[i] << G4endl; >> 2024 i++; >> 2025 } >> 2026 >> 2027 i=0; >> 2028 while(i < maxi) >> 2029 { >> 2030 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; // normalisation >> 2031 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); >> 2032 i++; >> 2033 } >> 2034 >> 2035 if(verbosityLevel >= 1) >> 2036 { >> 2037 G4cout << "Leaving LinearInterpolation" << G4endl; >> 2038 ArbEnergyH.DumpValues(); >> 2039 IPDFArbEnergyH.DumpValues(); >> 2040 } >> 2041 } >> 2042 >> 2043 void G4GeneralParticleSource::LogInterpolation() >> 2044 { >> 2045 // Interpolation based on Logarithmic equations >> 2046 // Generate equations of line segments >> 2047 // y = Ax**alpha => log y = alpha*logx + logA >> 2048 // Find area under line segments >> 2049 // create normalised, cumulative array Arb_Cum_Area >> 2050 G4double Area_seg[256], Arb_x[256], Arb_y[256], Arb_Cum_Area[256]; >> 2051 G4double alp, sum=0.; >> 2052 Arb_Cum_Area[0] = 0.; >> 2053 if(verbosityLevel == 2) >> 2054 ArbEnergyH.DumpValues(); >> 2055 >> 2056 G4int i, count; >> 2057 G4int maxi = ArbEnergyH.GetVectorLength(); >> 2058 for(i=0;i<maxi;i++) >> 2059 { >> 2060 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); >> 2061 Arb_y[i] = ArbEnergyH(size_t(i)); >> 2062 } >> 2063 // Points are now in x,y arrays. If the spectrum is integral it has to be >> 2064 // made differential and if momentum it has to be made energy. >> 2065 if(EnergySpec == false) >> 2066 { >> 2067 // change currently stored values (emin etc) which are actually momenta >> 2068 // to energies. >> 2069 if(particle_definition == NULL) >> 2070 G4cout << "Error: particle not defined" << G4endl; >> 2071 else >> 2072 { >> 2073 // Apply Energy**2 = p**2c**2 + m0**2c**4 >> 2074 // p should be entered as E/c i.e. without the division by c >> 2075 // being done - energy equivalent. >> 2076 G4double mass = particle_definition->GetPDGMass(); >> 2077 >> 2078 // multiply the function (Arb_y) up by the bin width >> 2079 // to make the function counts/s (i.e. get rid of momentum >> 2080 // dependence). >> 2081 for(count=0;count<maxi;count++) >> 2082 { >> 2083 Arb_y[count] = Arb_y[count] * (Arb_x[count+1] - Arb_x[count]); >> 2084 } >> 2085 // Change Arb_x to energy, plus divide by energy bin width >> 2086 // to make Arb_y counts/s/energy >> 2087 for(count=0;count<maxi;count++) >> 2088 { >> 2089 Arb_x[count] = sqrt((Arb_x[count]*Arb_x[count]) >> 2090 + (mass*mass)); >> 2091 } >> 2092 for(count=0;count<maxi;count++) >> 2093 { >> 2094 Arb_y[count] = Arb_y[count]/(Arb_x[count+1] - Arb_x[count]); >> 2095 } >> 2096 } >> 2097 } >> 2098 if(DiffSpec == false) >> 2099 { >> 2100 // Converts integral point-wise spectra to Differential >> 2101 for(count=1;count<=maxi;count++) >> 2102 { >> 2103 Arb_y[count] = Arb_y[count] - Arb_y[count-1]; >> 2104 } >> 2105 } >> 2106 >> 2107 i=1; >> 2108 Arb_alpha[0] = 0.; >> 2109 Arb_Const[0] = 0.; >> 2110 Area_seg[0] = 0.; >> 2111 if(Arb_x[0] <= 0. || Arb_y[0] <= 0.) >> 2112 { >> 2113 G4cout << "You should not use log interpolation with points <= 0." << G4endl; >> 2114 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl; >> 2115 if(Arb_x[0] <= 0.) >> 2116 Arb_x[0] = 1e-20; >> 2117 if(Arb_y[0] <= 0.) >> 2118 Arb_y[0] = 1e-20; >> 2119 } >> 2120 >> 2121 while(i <maxi) >> 2122 { >> 2123 // Incase points are negative or zero >> 2124 if(Arb_x[i] <= 0. || Arb_y[i] <= 0.) >> 2125 { >> 2126 G4cout << "You should not use log interpolation with points <= 0." << G4endl; >> 2127 G4cout << "These will be changed to 1e-20, which may cause problems" << G4endl; >> 2128 if(Arb_x[i] <= 0.) >> 2129 Arb_x[i] = 1e-20; >> 2130 if(Arb_y[i] <= 0.) >> 2131 Arb_y[i] = 1e-20; >> 2132 } >> 2133 >> 2134 Arb_alpha[i] = (log10(Arb_y[i])-log10(Arb_y[i-1]))/(log10(Arb_x[i])-log10(Arb_x[i-1])); >> 2135 Arb_Const[i] = Arb_y[i]/(pow(Arb_x[i],Arb_alpha[i])); >> 2136 alp = Arb_alpha[i] + 1; >> 2137 Area_seg[i] = (Arb_Const[i]/alp) * (pow(Arb_x[i],alp) - pow(Arb_x[i-1],alp)); >> 2138 sum = sum + Area_seg[i]; >> 2139 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 2140 if(verbosityLevel == 2) >> 2141 G4cout << Arb_alpha[i] << Arb_Const[i] << Area_seg[i] << G4endl; >> 2142 i++; >> 2143 } >> 2144 >> 2145 i=0; >> 2146 while(i<maxi) >> 2147 { >> 2148 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; >> 2149 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); >> 2150 i++; >> 2151 } >> 2152 if(verbosityLevel >= 1) >> 2153 G4cout << "Leaving LogInterpolation " << G4endl; >> 2154 } >> 2155 >> 2156 void G4GeneralParticleSource::ExpInterpolation() >> 2157 { >> 2158 // Interpolation based on Exponential equations >> 2159 // Generate equations of line segments >> 2160 // y = Ae**-(x/e0) => ln y = -x/e0 + lnA >> 2161 // Find area under line segments >> 2162 // create normalised, cumulative array Arb_Cum_Area >> 2163 G4double Area_seg[256], Arb_x[256], Arb_y[256], Arb_Cum_Area[256]; >> 2164 G4double sum=0.; >> 2165 Arb_Cum_Area[0] = 0.; >> 2166 if(verbosityLevel == 2) >> 2167 ArbEnergyH.DumpValues(); >> 2168 >> 2169 G4int i, count; >> 2170 G4int maxi = ArbEnergyH.GetVectorLength(); >> 2171 for(i=0;i<maxi;i++) >> 2172 { >> 2173 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); >> 2174 Arb_y[i] = ArbEnergyH(size_t(i)); >> 2175 } >> 2176 // Points are now in x,y arrays. If the spectrum is integral it has to be >> 2177 // made differential and if momentum it has to be made energy. >> 2178 if(EnergySpec == false) >> 2179 { >> 2180 // change currently stored values (emin etc) which are actually momenta >> 2181 // to energies. >> 2182 if(particle_definition == NULL) >> 2183 G4cout << "Error: particle not defined" << G4endl; >> 2184 else >> 2185 { >> 2186 // Apply Energy**2 = p**2c**2 + m0**2c**4 >> 2187 // p should be entered as E/c i.e. without the division by c >> 2188 // being done - energy equivalent. >> 2189 G4double mass = particle_definition->GetPDGMass(); >> 2190 >> 2191 // multiply the function (Arb_y) up by the bin width >> 2192 // to make the function counts/s (i.e. get rid of momentum >> 2193 // dependence). >> 2194 for( count=0;count<maxi;count++) >> 2195 { >> 2196 Arb_y[count] = Arb_y[count] * (Arb_x[count+1] - Arb_x[count]); >> 2197 } >> 2198 // Change Arb_x to energy, plus divide by energy bin width >> 2199 // to make Arb_y counts/s/energy >> 2200 for(count=0;count<maxi;count++) >> 2201 { >> 2202 Arb_x[count] = sqrt((Arb_x[count]*Arb_x[count]) >> 2203 + (mass*mass)); >> 2204 } >> 2205 for(count=0;count<maxi;count++) >> 2206 { >> 2207 Arb_y[count] = Arb_y[count]/(Arb_x[count+1] - Arb_x[count]); >> 2208 } >> 2209 } >> 2210 } >> 2211 if(DiffSpec == false) >> 2212 { >> 2213 // Converts integral point-wise spectra to Differential >> 2214 for(count=1;count<=maxi;count++) >> 2215 { >> 2216 Arb_y[count] = Arb_y[count] - Arb_y[count-1]; >> 2217 } >> 2218 } >> 2219 >> 2220 i=1; >> 2221 Arb_ezero[0] = 0.; >> 2222 Arb_Const[0] = 0.; >> 2223 Area_seg[0] = 0.; >> 2224 while(i < maxi) >> 2225 { >> 2226 G4double test = log(Arb_y[i]) - log(Arb_y[i-1]); >> 2227 if(test > 0. || test < 0.) >> 2228 { >> 2229 Arb_ezero[i] = -(Arb_x[i] - Arb_x[i-1])/(log(Arb_y[i]) - log(Arb_y[i-1])); >> 2230 Arb_Const[i] = Arb_y[i]/(exp(-Arb_x[i]/Arb_ezero[i])); >> 2231 Area_seg[i]=-(Arb_Const[i]*Arb_ezero[i])*(exp(-Arb_x[i]/Arb_ezero[i]) >> 2232 -exp(-Arb_x[i-1]/Arb_ezero[i])); >> 2233 } >> 2234 else >> 2235 { >> 2236 G4cout << "Flat line segment: problem" << G4endl; >> 2237 Arb_ezero[i] = 0.; >> 2238 Arb_Const[i] = 0.; >> 2239 Area_seg[i] = 0.; >> 2240 } >> 2241 sum = sum + Area_seg[i]; >> 2242 Arb_Cum_Area[i] = Arb_Cum_Area[i-1] + Area_seg[i]; >> 2243 if(verbosityLevel == 2) >> 2244 G4cout << Arb_ezero[i] << Arb_Const[i] << Area_seg[i] << G4endl; >> 2245 i++; >> 2246 } >> 2247 >> 2248 i=0; >> 2249 while(i<maxi) >> 2250 { >> 2251 Arb_Cum_Area[i] = Arb_Cum_Area[i]/sum; >> 2252 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_Cum_Area[i]); >> 2253 i++; >> 2254 } >> 2255 if(verbosityLevel >= 1) >> 2256 G4cout << "Leaving ExpInterpolation " << G4endl; >> 2257 } >> 2258 >> 2259 void G4GeneralParticleSource::SplineInterpolation() >> 2260 { >> 2261 // Interpolation using Splines. >> 2262 // Create Normalised arrays, make x 0->1 and y hold >> 2263 // the function (Energy) >> 2264 G4double Arb_x[256], Arb_y[256]; >> 2265 G4double sum = 0.; >> 2266 G4int i, count; >> 2267 if(verbosityLevel == 2) >> 2268 ArbEnergyH.DumpValues(); >> 2269 >> 2270 G4int maxi = ArbEnergyH.GetVectorLength(); >> 2271 for(i=0;i<maxi;i++) >> 2272 { >> 2273 Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); >> 2274 // if(i == 0) >> 2275 Arb_y[i] = ArbEnergyH(size_t(i)); >> 2276 // else >> 2277 // Arb_y[i] = Arb_y[i-1] + ArbEnergyH(size_t(i)); >> 2278 // sum = sum + ArbEnergyH(size_t(i)); >> 2279 } >> 2280 // Points are now in x,y arrays. If the spectrum is integral it has to be >> 2281 // made differential and if momentum it has to be made energy. >> 2282 if(EnergySpec == false) >> 2283 { >> 2284 // change currently stored values (emin etc) which are actually momenta >> 2285 // to energies. >> 2286 if(particle_definition == NULL) >> 2287 G4cout << "Error: particle not defined" << G4endl; >> 2288 else >> 2289 { >> 2290 // Apply Energy**2 = p**2c**2 + m0**2c**4 >> 2291 // p should be entered as E/c i.e. without the division by c >> 2292 // being done - energy equivalent. >> 2293 G4double mass = particle_definition->GetPDGMass(); >> 2294 >> 2295 // multiply the function (Arb_y) up by the bin width >> 2296 // to make the function counts/s (i.e. get rid of momentum >> 2297 // dependence). >> 2298 for(count=0;count<maxi;count++) >> 2299 { >> 2300 Arb_y[count] = Arb_y[count] * (Arb_x[count+1] - Arb_x[count]); >> 2301 } >> 2302 // Change Arb_x to energy, plus divide by energy bin width >> 2303 // to make Arb_y counts/s/energy >> 2304 for(count=0;count<maxi;count++) >> 2305 { >> 2306 Arb_x[count] = sqrt((Arb_x[count]*Arb_x[count]) >> 2307 + (mass*mass)); >> 2308 } >> 2309 for(count=0;count<maxi;count++) >> 2310 { >> 2311 Arb_y[count] = Arb_y[count]/(Arb_x[count+1] - Arb_x[count]); >> 2312 } >> 2313 } >> 2314 } >> 2315 if(DiffSpec == false) >> 2316 { >> 2317 // Converts integral point-wise spectra to Differential >> 2318 for(count=1;count<=maxi;count++) >> 2319 { >> 2320 Arb_y[count] = Arb_y[count] - Arb_y[count-1]; >> 2321 } >> 2322 } >> 2323 >> 2324 for(i=1;i<maxi;i++) >> 2325 { >> 2326 // Arb_x[i] = ArbEnergyH.GetLowEdgeEnergy(size_t(i)); >> 2327 // if(i == 0) >> 2328 // Arb_y[i] = ArbEnergyH(size_t(i)); >> 2329 // else >> 2330 Arb_y[i] = Arb_y[i-1] + ArbEnergyH(size_t(i)); >> 2331 sum = sum + ArbEnergyH(size_t(i)); >> 2332 } >> 2333 >> 2334 for(i=0;i<maxi;i++) >> 2335 { >> 2336 Arb_y[i] = Arb_y[i]/sum; >> 2337 } >> 2338 >> 2339 if((Arb_x[0] != 0.) || (Arb_y[0] != 0.)) >> 2340 { >> 2341 for(i=maxi;i>=0;i--) >> 2342 { >> 2343 if(i == 0) >> 2344 { >> 2345 Arb_x[0] = 0.; >> 2346 Arb_y[0] = 0.; >> 2347 } >> 2348 else >> 2349 { >> 2350 Arb_x[i] = Arb_x[i-1]; >> 2351 Arb_y[i] = Arb_y[i-1]; >> 2352 } >> 2353 } >> 2354 } >> 2355 >> 2356 for(i=0;i<=maxi;i++) >> 2357 { >> 2358 if(verbosityLevel == 2) >> 2359 G4cout << i <<" "<< Arb_x[i] << " " << Arb_y[i] << G4endl; >> 2360 IPDFArbEnergyH.InsertValues(Arb_x[i], Arb_y[i]); >> 2361 } >> 2362 >> 2363 // Should now have normalised cumulative probabilities in Arb_y >> 2364 // and energy values in Arb_x. >> 2365 maxi = maxi + 1; >> 2366 // Put y into x and x into y. The spline interpolation will then >> 2367 // go through x-axis to find where to interpolate (cum probability) >> 2368 // then generate a y (which will now be energy). >> 2369 SplineInt = new G4DataInterpolation(Arb_y,Arb_x,maxi,1e30,1e30); >> 2370 if(verbosityLevel == 2) >> 2371 { >> 2372 G4cout << SplineInt << G4endl; >> 2373 G4cout << SplineInt->LocateArgument(1.0) << G4endl; >> 2374 } >> 2375 if(verbosityLevel >= 1) >> 2376 G4cout << "Leaving SplineInterpolation " << G4endl; >> 2377 } >> 2378 >> 2379 void G4GeneralParticleSource::GenerateMonoEnergetic() >> 2380 { >> 2381 // Method to generate MonoEnergetic particles. >> 2382 particle_energy = MonoEnergy; >> 2383 } >> 2384 >> 2385 void G4GeneralParticleSource::GenerateGaussEnergies() >> 2386 { >> 2387 // Method to generate Gaussian particles. >> 2388 particle_energy = G4RandGauss::shoot(MonoEnergy,SE); >> 2389 if (particle_energy < 0) particle_energy = 0.; >> 2390 } >> 2391 >> 2392 void G4GeneralParticleSource::GenerateLinearEnergies() >> 2393 { >> 2394 G4double rndm; >> 2395 G4double emaxsq = pow(Emax,2.); //Emax squared >> 2396 G4double eminsq = pow(Emin,2.); //Emin squared >> 2397 G4double intersq = pow(cept,2.); //cept squared >> 2398 >> 2399 rndm = GenRandEnergy(); >> 2400 >> 2401 G4double bracket = ((grad/2.)*(emaxsq - eminsq) + cept*(Emax-Emin)); >> 2402 bracket = bracket * rndm; >> 2403 bracket = bracket + (grad/2.)*eminsq + cept*Emin; >> 2404 // Now have a quad of form m/2 E**2 + cE - bracket = 0 >> 2405 bracket = -bracket; >> 2406 // G4cout << "BRACKET" << bracket << G4endl; >> 2407 if(grad != 0.) >> 2408 { >> 2409 G4double sqbrack = (intersq - 4*(grad/2.)*(bracket)); >> 2410 // G4cout << "SQBRACK" << sqbrack << G4endl; >> 2411 sqbrack = sqrt(sqbrack); >> 2412 G4double root1 = -cept + sqbrack; >> 2413 root1 = root1/(2.*(grad/2.)); >> 2414 >> 2415 G4double root2 = -cept - sqbrack; >> 2416 root2 = root2/(2.*(grad/2.)); >> 2417 >> 2418 // G4cout << root1 << " roots " << root2 << G4endl; >> 2419 >> 2420 if(root1 > Emin && root1 < Emax) >> 2421 particle_energy = root1; >> 2422 if(root2 > Emin && root2 < Emax) >> 2423 particle_energy = root2; >> 2424 } >> 2425 else if(grad == 0.) >> 2426 // have equation of form cE - bracket =0 >> 2427 particle_energy = bracket/cept; >> 2428 >> 2429 if(particle_energy < 0.) >> 2430 particle_energy = -particle_energy; >> 2431 >> 2432 if(verbosityLevel >= 1) >> 2433 G4cout << "Energy is " << particle_energy << G4endl; >> 2434 } >> 2435 >> 2436 void G4GeneralParticleSource::GeneratePowEnergies() >> 2437 { >> 2438 // Method to generate particle energies distributed as >> 2439 // a powerlaw >> 2440 >> 2441 G4double rndm; >> 2442 G4double emina, emaxa; >> 2443 >> 2444 emina = pow(Emin,alpha+1); >> 2445 emaxa = pow(Emax,alpha+1); >> 2446 >> 2447 rndm = GenRandEnergy(); >> 2448 >> 2449 if(alpha != -1.) >> 2450 { >> 2451 particle_energy = ((rndm*(emaxa - emina)) + emina); >> 2452 particle_energy = pow(particle_energy,(1./(alpha+1.))); >> 2453 } >> 2454 else if(alpha == -1.) >> 2455 { >> 2456 particle_energy = (log(Emin) + rndm*(log(Emax) - log(Emin))); >> 2457 particle_energy = exp(particle_energy); >> 2458 } >> 2459 if(verbosityLevel >= 1) >> 2460 G4cout << "Energy is " << particle_energy << G4endl; >> 2461 } >> 2462 >> 2463 void G4GeneralParticleSource::GenerateExpEnergies() >> 2464 { >> 2465 // Method to generate particle energies distributed according >> 2466 // to an exponential curve. >> 2467 G4double rndm; >> 2468 rndm = GenRandEnergy(); >> 2469 >> 2470 particle_energy = -Ezero*(log(rndm*(exp(-Emax/Ezero) - exp(-Emin/Ezero)) + >> 2471 exp(-Emin/Ezero))); >> 2472 if(verbosityLevel >= 1) >> 2473 G4cout << "Energy is " << particle_energy << G4endl; >> 2474 } >> 2475 >> 2476 void G4GeneralParticleSource::GenerateBremEnergies() >> 2477 { >> 2478 // Method to generate particle energies distributed according >> 2479 // to a Bremstrahlung equation of >> 2480 // form I = const*((kT)**1/2)*E*(e**(-E/kT)) >> 2481 >> 2482 G4double rndm; >> 2483 rndm = GenRandEnergy(); >> 2484 G4double expmax, expmin, k; >> 2485 >> 2486 k = 8.6181e-11; // Boltzmann's const in MeV/K >> 2487 G4double ksq = pow(k,2.); // k squared >> 2488 G4double Tsq = pow(Temp,2.); // Temp squared >> 2489 >> 2490 expmax = exp(-Emax/(k*Temp)); >> 2491 expmin = exp(-Emin/(k*Temp)); >> 2492 >> 2493 // If either expmax or expmin are zero then this will cause problems >> 2494 // Most probably this will be because T is too low or E is too high >> 2495 >> 2496 if(expmax == 0.) >> 2497 G4cout << "*****EXPMAX=0. Choose different E's or Temp" << G4endl; >> 2498 if(expmin == 0.) >> 2499 G4cout << "*****EXPMIN=0. Choose different E's or Temp" << G4endl; >> 2500 >> 2501 G4double tempvar = rndm *((-k)*Temp*(Emax*expmax - Emin*expmin) - >> 2502 (ksq*Tsq*(expmax-expmin))); >> 2503 >> 2504 G4double bigc = (tempvar - k*Temp*Emin*expmin - ksq*Tsq*expmin)/(-k*Temp); >> 2505 >> 2506 // This gives an equation of form: Ee(-E/kT) + kTe(-E/kT) - C =0 >> 2507 // Solve this iteratively, step from Emin to Emax in 1000 steps >> 2508 // and take the best solution. >> 2509 >> 2510 G4double erange = Emax - Emin; >> 2511 G4double steps = erange/1000.; >> 2512 G4int i; >> 2513 G4double etest, diff, err; >> 2514 >> 2515 err = 100000.; >> 2516 >> 2517 for(i=1; i<1000; i++) >> 2518 { >> 2519 etest = Emin + (i-1)*steps; >> 2520 >> 2521 diff = etest*(exp(-etest/(k*Temp))) + k*Temp*(exp(-etest/(k*Temp))) - bigc; >> 2522 >> 2523 if(diff < 0.) >> 2524 diff = -diff; >> 2525 >> 2526 if(diff < err) >> 2527 { >> 2528 err = diff; >> 2529 particle_energy = etest; >> 2530 } >> 2531 } >> 2532 if(verbosityLevel >= 1) >> 2533 G4cout << "Energy is " << particle_energy << G4endl; >> 2534 } >> 2535 >> 2536 void G4GeneralParticleSource::GenerateBbodyEnergies() >> 2537 { >> 2538 // BBody_x holds Energies, and BBHist holds the cumulative histo. >> 2539 // binary search to find correct bin then lin interpolation. >> 2540 // Use the earlier defined histogram + RandGeneral method to generate >> 2541 // random numbers following the histos distribution. >> 2542 G4double rndm; >> 2543 G4int nabove, nbelow = 0, middle; >> 2544 nabove = 10001; >> 2545 rndm = GenRandEnergy(); >> 2546 >> 2547 // Binary search to find bin that rndm is in >> 2548 while(nabove-nbelow > 1) >> 2549 { >> 2550 middle = (nabove + nbelow)/2; >> 2551 if(rndm == BBHist[middle]) break; >> 2552 if(rndm < BBHist[middle]) nabove = middle; >> 2553 else nbelow = middle; >> 2554 } >> 2555 >> 2556 // Now interpolate in that bin to find the correct output value. >> 2557 G4double x1, x2, y1, y2, m, q; >> 2558 x1 = Bbody_x[nbelow]; >> 2559 x2 = Bbody_x[nbelow+1]; >> 2560 y1 = BBHist[nbelow]; >> 2561 y2 = BBHist[nbelow+1]; >> 2562 m = (y2-y1)/(x2-x1); >> 2563 q = y1 - m*x1; >> 2564 >> 2565 particle_energy = (rndm - q)/m; >> 2566 >> 2567 if(verbosityLevel >= 1) >> 2568 { >> 2569 G4cout << "Energy is " << particle_energy << G4endl; >> 2570 } >> 2571 } >> 2572 >> 2573 void G4GeneralParticleSource::GenerateCdgEnergies() >> 2574 { >> 2575 // Gen random numbers, compare with values in cumhist >> 2576 // to find appropriate part of spectrum and then >> 2577 // generate energy in the usual inversion way. >> 2578 // G4double pfact[2] = {8.5, 112}; >> 2579 // G4double spind[2] = {1.4, 2.3}; >> 2580 // G4double ene_line[3] = {1., 18., 1E6}; >> 2581 G4double rndm, rndm2; >> 2582 G4double ene_line[3]; >> 2583 G4double omalpha[2]; >> 2584 if(Emin < 18*keV && Emax < 18*keV) >> 2585 { >> 2586 omalpha[0] = 1. - 1.4; >> 2587 ene_line[0] = Emin; >> 2588 ene_line[1] = Emax; >> 2589 } >> 2590 if(Emin < 18*keV && Emax > 18*keV) >> 2591 { >> 2592 omalpha[0] = 1. - 1.4; >> 2593 omalpha[1] = 1. - 2.3; >> 2594 ene_line[0] = Emin; >> 2595 ene_line[1] = 18.; >> 2596 ene_line[2] = Emax; >> 2597 } >> 2598 if(Emin > 18*keV) >> 2599 { >> 2600 omalpha[0] = 1. - 2.3; >> 2601 ene_line[0] = Emin; >> 2602 ene_line[1] = Emax; >> 2603 } >> 2604 rndm = GenRandEnergy(); >> 2605 rndm2 = GenRandEnergy(); >> 2606 >> 2607 G4int i = 0; >> 2608 while( rndm >= CDGhist[i]) >> 2609 { >> 2610 i++; >> 2611 } >> 2612 // Generate final energy. >> 2613 particle_energy = (pow(ene_line[i-1],omalpha[i-1]) + (pow(ene_line[i],omalpha[i-1]) >> 2614 - pow(ene_line[i-1],omalpha[i-1]))*rndm2); >> 2615 particle_energy = pow(particle_energy,(1./omalpha[i-1])); >> 2616 >> 2617 if(verbosityLevel >= 1) >> 2618 G4cout << "Energy is " << particle_energy << G4endl; >> 2619 } >> 2620 >> 2621 void G4GeneralParticleSource::GenUserHistEnergies() >> 2622 { >> 2623 // Histograms are DIFFERENTIAL. >> 2624 // G4cout << "In GenUserHistEnergies " << G4endl; >> 2625 if(IPDFEnergyExist == false) >> 2626 { >> 2627 G4int ii; >> 2628 G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); >> 2629 G4double bins[256], vals[256], sum; >> 2630 sum=0.; >> 2631 // UDefEnergyH.DumpValues(); >> 2632 G4double mass = particle_definition->GetPDGMass(); >> 2633 // G4cout << mass << G4endl; >> 2634 // G4cout << EnergySpec << " " << DiffSpec << G4endl; >> 2635 if((EnergySpec == false) && (particle_definition == NULL)) >> 2636 G4cout << "Error: particle definition is NULL" << G4endl; >> 2637 >> 2638 if(maxbin > 256) >> 2639 { >> 2640 G4cout << "Maxbin > 256" << G4endl; >> 2641 G4cout << "Setting maxbin to 256, other bins are lost" << G4endl; >> 2642 } >> 2643 >> 2644 if(DiffSpec == false) >> 2645 G4cout << "Histograms are Differential!!! " << G4endl; >> 2646 else >> 2647 { >> 2648 // G4cout << "Here 2" << G4endl; >> 2649 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); >> 2650 vals[0] = UDefEnergyH(size_t(0)); >> 2651 sum = vals[0]; >> 2652 for(ii=1;ii<maxbin;ii++) >> 2653 { >> 2654 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); >> 2655 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1]; >> 2656 sum = sum + UDefEnergyH(size_t(ii)); >> 2657 } >> 2658 } >> 2659 >> 2660 if(EnergySpec == false) >> 2661 { >> 2662 // multiply the function (vals) up by the bin width >> 2663 // to make the function counts/s (i.e. get rid of momentum >> 2664 // dependence). >> 2665 for(ii=1;ii<maxbin;ii++) >> 2666 { >> 2667 // G4cout << vals[ii] << " " << bins[ii] << " " << bins[ii-1] << G4endl; >> 2668 vals[ii] = vals[ii] * (bins[ii] - bins[ii-1]); >> 2669 } >> 2670 // Put energy bins into new histo, plus divide by energy bin width >> 2671 // to make evals counts/s/energy >> 2672 for(ii=0;ii<maxbin;ii++) >> 2673 { >> 2674 bins[ii] = sqrt((bins[ii]*bins[ii]) + (mass*mass)); >> 2675 // G4cout << bins[ii] << " " << mass << G4endl; >> 2676 } >> 2677 for(ii=1;ii<maxbin;ii++) >> 2678 { >> 2679 // G4cout << vals[ii] << " " << bins[ii] << " " << bins[ii-1] << G4endl; >> 2680 vals[ii] = vals[ii]/(bins[ii] - bins[ii-1]); >> 2681 } >> 2682 sum = vals[maxbin-1]; >> 2683 } >> 2684 >> 2685 for(ii=0;ii<maxbin;ii++) >> 2686 { >> 2687 vals[ii] = vals[ii]/sum; >> 2688 IPDFEnergyH.InsertValues(bins[ii], vals[ii]); >> 2689 } >> 2690 >> 2691 // Make IPDFEnergyExist = true >> 2692 IPDFEnergyExist = true; >> 2693 if(verbosityLevel == 2) >> 2694 IPDFEnergyH.DumpValues(); >> 2695 >> 2696 } >> 2697 >> 2698 // IPDF has been create so carry on >> 2699 G4double rndm = GenRandEnergy(); >> 2700 particle_energy = IPDFEnergyH.GetEnergy(rndm); >> 2701 >> 2702 if(verbosityLevel >= 1) >> 2703 G4cout << "Energy is " << particle_energy << G4endl; >> 2704 } >> 2705 >> 2706 void G4GeneralParticleSource::GenArbPointEnergies() >> 2707 { >> 2708 if(verbosityLevel >= 1) >> 2709 G4cout << "In GenArbPointEnergies" << G4endl; >> 2710 G4double rndm; >> 2711 rndm = GenRandEnergy(); >> 2712 >> 2713 if(IntType != "Spline") >> 2714 { >> 2715 // IPDFArbEnergyH.DumpValues(); >> 2716 // Find the Bin >> 2717 // have x, y, no of points, and cumulative area distribution >> 2718 G4int nabove, nbelow = 0, middle; >> 2719 nabove = IPDFArbEnergyH.GetVectorLength(); >> 2720 // G4cout << nabove << G4endl; >> 2721 // Binary search to find bin that rndm is in >> 2722 while(nabove-nbelow > 1) >> 2723 { >> 2724 middle = (nabove + nbelow)/2; >> 2725 if(rndm == IPDFArbEnergyH(size_t(middle))) break; >> 2726 if(rndm < IPDFArbEnergyH(size_t(middle))) nabove = middle; >> 2727 else nbelow = middle; >> 2728 } >> 2729 if(IntType == "Lin") >> 2730 { >> 2731 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 2732 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 2733 grad = Arb_grad[nbelow+1]; >> 2734 cept = Arb_cept[nbelow+1]; >> 2735 // G4cout << rndm << " " << Emax << " " << Emin << " " << grad << " " << cept << G4endl; >> 2736 GenerateLinearEnergies(); >> 2737 } >> 2738 else if(IntType == "Log") >> 2739 { >> 2740 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 2741 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 2742 alpha = Arb_alpha[nbelow+1]; >> 2743 // G4cout << rndm << " " << Emax << " " << Emin << " " << alpha << G4endl; >> 2744 GeneratePowEnergies(); >> 2745 } >> 2746 else if(IntType == "Exp") >> 2747 { >> 2748 Emax = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow+1)); >> 2749 Emin = IPDFArbEnergyH.GetLowEdgeEnergy(size_t(nbelow)); >> 2750 Ezero = Arb_ezero[nbelow+1]; >> 2751 // G4cout << rndm << " " << Emax << " " << Emin << " " << Ezero << G4endl; >> 2752 GenerateExpEnergies(); >> 2753 } >> 2754 } >> 2755 else if(IntType == "Spline") >> 2756 { >> 2757 if(verbosityLevel == 2) >> 2758 G4cout << "IntType = Spline " << rndm << G4endl; >> 2759 // in SplineInterpolation created SplineInt >> 2760 // Now generate a random number put it into CubicSplineInterpolation >> 2761 // and you should get out an energy!?! >> 2762 particle_energy = SplineInt->CubicSplineInterpolation(rndm); >> 2763 if(verbosityLevel >= 1) >> 2764 G4cout << "Energy is " << particle_energy << G4endl; >> 2765 } >> 2766 else >> 2767 G4cout << "Error: IntType unknown type" << G4endl; >> 2768 } >> 2769 >> 2770 void G4GeneralParticleSource::GenEpnHistEnergies() >> 2771 { >> 2772 // G4cout << "In GenEpnHistEnergies " << Epnflag << G4endl; >> 2773 >> 2774 // Firstly convert to energy if not already done. >> 2775 if(Epnflag == true) >> 2776 // epnflag = true means spectrum is epn, false means e. >> 2777 { >> 2778 // convert to energy by multiplying by A number >> 2779 ConvertEPNToEnergy(); >> 2780 // EpnEnergyH will be replace by UDefEnergyH. >> 2781 // UDefEnergyH.DumpValues(); >> 2782 } >> 2783 >> 2784 // G4cout << "Creating IPDFEnergy if not already done so" << G4endl; >> 2785 if(IPDFEnergyExist == false) >> 2786 { >> 2787 // IPDF has not been created, so create it >> 2788 G4double bins[256],vals[256], sum; >> 2789 G4int ii; >> 2790 G4int maxbin = G4int(UDefEnergyH.GetVectorLength()); >> 2791 bins[0] = UDefEnergyH.GetLowEdgeEnergy(size_t(0)); >> 2792 vals[0] = UDefEnergyH(size_t(0)); >> 2793 sum = vals[0]; >> 2794 for(ii=1;ii<maxbin;ii++) >> 2795 { >> 2796 bins[ii] = UDefEnergyH.GetLowEdgeEnergy(size_t(ii)); >> 2797 vals[ii] = UDefEnergyH(size_t(ii)) + vals[ii-1]; >> 2798 sum = sum + UDefEnergyH(size_t(ii)); >> 2799 } >> 2800 >> 2801 for(ii=0;ii<maxbin;ii++) >> 2802 { >> 2803 vals[ii] = vals[ii]/sum; >> 2804 IPDFEnergyH.InsertValues(bins[ii], vals[ii]); >> 2805 } >> 2806 // Make IPDFEpnExist = true >> 2807 IPDFEnergyExist = true; >> 2808 } >> 2809 // IPDFEnergyH.DumpValues(); >> 2810 // IPDF has been create so carry on >> 2811 G4double rndm = GenRandEnergy(); >> 2812 particle_energy = IPDFEnergyH.GetEnergy(rndm); >> 2813 >> 2814 if(verbosityLevel >= 1) >> 2815 G4cout << "Energy is " << particle_energy << G4endl; >> 2816 } >> 2817 >> 2818 void G4GeneralParticleSource::ConvertEPNToEnergy() >> 2819 { >> 2820 // Use this before particle generation to convert the >> 2821 // currently stored histogram from energy/nucleon >> 2822 // to energy. >> 2823 // G4cout << "In ConvertEpntoEnergy " << G4endl; >> 2824 if(particle_definition==NULL) >> 2825 G4cout << "Error: particle not defined" << G4endl; >> 2826 else >> 2827 { >> 2828 // Need to multiply histogram by the number of nucleons. >> 2829 // Baryon Number looks to hold the no. of nucleons. >> 2830 G4int Bary = particle_definition->GetBaryonNumber(); >> 2831 // G4cout << "Baryon No. " << Bary << G4endl; >> 2832 // Change values in histogram, Read it out, delete it, re-create it >> 2833 G4int count, maxcount; >> 2834 maxcount = G4int(EpnEnergyH.GetVectorLength()); >> 2835 // G4cout << maxcount << G4endl; >> 2836 G4double ebins[256],evals[256]; >> 2837 if(maxcount > 256) >> 2838 { >> 2839 G4cout << "Histogram contains more than 256 bins!" << G4endl; >> 2840 G4cout << "Those above 256 will be ignored" << G4endl; >> 2841 maxcount = 256; >> 2842 } >> 2843 for(count=0;count<maxcount;count++) >> 2844 { >> 2845 // Read out >> 2846 ebins[count] = EpnEnergyH.GetLowEdgeEnergy(size_t(count)); >> 2847 evals[count] = EpnEnergyH(size_t(count)); >> 2848 } >> 2849 >> 2850 // Multiply the channels by the nucleon number to give energies >> 2851 for(count=0;count<maxcount;count++) >> 2852 { >> 2853 ebins[count] = ebins[count] * Bary; >> 2854 } >> 2855 >> 2856 // Set Emin and Emax >> 2857 Emin = ebins[0]; >> 2858 Emax = ebins[maxcount-1]; >> 2859 >> 2860 // Put energy bins into new histogram - UDefEnergyH. >> 2861 for(count=0;count<maxcount;count++) >> 2862 { >> 2863 UDefEnergyH.InsertValues(ebins[count], evals[count]); >> 2864 } >> 2865 Epnflag = false; //so that you dont repeat this method. >> 2866 } >> 2867 } >> 2868 >> 2869 // Biasing methods >> 2870 >> 2871 void G4GeneralParticleSource::SetXBias(G4ThreeVector input) >> 2872 { >> 2873 G4double ehi, val; >> 2874 ehi = input.x(); >> 2875 val = input.y(); >> 2876 XBiasH.InsertValues(ehi, val); >> 2877 XBias = true; >> 2878 } >> 2879 >> 2880 void G4GeneralParticleSource::SetYBias(G4ThreeVector input) >> 2881 { >> 2882 G4double ehi, val; >> 2883 ehi = input.x(); >> 2884 val = input.y(); >> 2885 YBiasH.InsertValues(ehi, val); >> 2886 YBias = true; >> 2887 } >> 2888 >> 2889 void G4GeneralParticleSource::SetZBias(G4ThreeVector input) >> 2890 { >> 2891 G4double ehi, val; >> 2892 ehi = input.x(); >> 2893 val = input.y(); >> 2894 ZBiasH.InsertValues(ehi, val); >> 2895 ZBias = true; >> 2896 } >> 2897 >> 2898 void G4GeneralParticleSource::SetThetaBias(G4ThreeVector input) >> 2899 { >> 2900 G4double ehi, val; >> 2901 ehi = input.x(); >> 2902 val = input.y(); >> 2903 ThetaBiasH.InsertValues(ehi, val); >> 2904 ThetaBias = true; >> 2905 } >> 2906 >> 2907 void G4GeneralParticleSource::SetPhiBias(G4ThreeVector input) >> 2908 { >> 2909 G4double ehi, val; >> 2910 ehi = input.x(); >> 2911 val = input.y(); >> 2912 PhiBiasH.InsertValues(ehi, val); >> 2913 PhiBias = true; >> 2914 } >> 2915 >> 2916 void G4GeneralParticleSource::SetEnergyBias(G4ThreeVector input) >> 2917 { >> 2918 G4double ehi, val; >> 2919 ehi = input.x(); >> 2920 val = input.y(); >> 2921 EnergyBiasH.InsertValues(ehi, val); >> 2922 EnergyBias = true; >> 2923 } >> 2924 >> 2925 >> 2926 void G4GeneralParticleSource::ReSetHist(G4String atype) >> 2927 { >> 2928 if (atype == "theta") { >> 2929 UDefThetaH = IPDFThetaH = ZeroPhysVector ; >> 2930 IPDFThetaExist = false ;} >> 2931 else if (atype == "phi"){ >> 2932 UDefPhiH = IPDFPhiH = ZeroPhysVector ; >> 2933 IPDFPhiExist = false ;} >> 2934 else if (atype == "energy"){ >> 2935 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; >> 2936 IPDFEnergyExist = false ;} >> 2937 else if ( atype == "arb"){ >> 2938 ArbEnergyH =IPDFArbEnergyH = ZeroPhysVector ; >> 2939 IPDFArbExist = false;} >> 2940 else if ( atype == "epn"){ >> 2941 UDefEnergyH = IPDFEnergyH = ZeroPhysVector ; >> 2942 IPDFEnergyExist = false ; >> 2943 EpnEnergyH = ZeroPhysVector ;} >> 2944 else if ( atype == "biasx") { >> 2945 XBias = false ; >> 2946 IPDFXBias = false; >> 2947 XBiasH = IPDFXBiasH = ZeroPhysVector ;} >> 2948 else if ( atype == "biasy") { >> 2949 YBias = false ; >> 2950 IPDFYBias = false; >> 2951 YBiasH = IPDFYBiasH = ZeroPhysVector ;} >> 2952 else if ( atype == "biasz") { >> 2953 ZBias = false ; >> 2954 IPDFZBias = false; >> 2955 ZBiasH = IPDFZBiasH = ZeroPhysVector ;} >> 2956 else if ( atype == "biast") { >> 2957 ThetaBias = false ; >> 2958 IPDFThetaBias = false; >> 2959 ThetaBiasH = IPDFThetaBiasH = ZeroPhysVector ;} >> 2960 else if ( atype == "biasp") { >> 2961 PhiBias = false ; >> 2962 IPDFPhiBias = false; >> 2963 PhiBiasH = IPDFPhiBiasH = ZeroPhysVector ;} >> 2964 else if ( atype == "biase") { >> 2965 EnergyBias = false ; >> 2966 IPDFEnergyBias = false; >> 2967 EnergyBiasH = IPDFEnergyBiasH = ZeroPhysVector ;} >> 2968 else { >> 2969 G4cout << "Error, histtype not accepted " << G4endl; 177 } 2970 } 178 } 2971 } 179 2972 180 void G4GeneralParticleSource::GeneratePrimaryV << 2973 >> 2974 >> 2975 G4double G4GeneralParticleSource::GenRandX() 181 { 2976 { 182 if (!GPSData->GetMultipleVertex()) << 2977 if(verbosityLevel >= 1) 183 { << 2978 G4cout << "In GenRandX" << G4endl; 184 G4SingleParticleSource* currentSource = GP << 2979 if(XBias == false) 185 if (GPSData->GetIntensityVectorSize() > 1) << 186 { 2980 { 187 // Try to minimize locks << 2981 // X is not biased 188 if (! normalised ) << 2982 G4double rndm = G4UniformRand(); 189 { << 2983 return(rndm); 190 // According to local variable, normal << 2984 } 191 // Check with underlying shared resour << 2985 else 192 // thread could have already normalize << 2986 { 193 GPSData->Lock(); << 2987 // X is biased 194 G4bool norm = GPSData->Normalised(); << 2988 if(IPDFXBias == false) 195 if (!norm) << 2989 { 196 { << 2990 // IPDF has not been created, so create it 197 IntensityNormalization(); << 2991 G4double bins[256],vals[256], sum; 198 } << 2992 G4int ii; 199 // This takes care of the case in whic << 2993 G4int maxbin = G4int(XBiasH.GetVectorLength()); 200 // is False and the underlying resourc << 2994 bins[0] = XBiasH.GetLowEdgeEnergy(size_t(0)); 201 normalised = GPSData->Normalised(); << 2995 vals[0] = XBiasH(size_t(0)); 202 GPSData->Unlock(); << 2996 sum = vals[0]; >> 2997 for(ii=1;ii<maxbin;ii++) >> 2998 { >> 2999 bins[ii] = XBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3000 vals[ii] = XBiasH(size_t(ii)) + vals[ii-1]; >> 3001 sum = sum + XBiasH(size_t(ii)); >> 3002 } >> 3003 >> 3004 for(ii=0;ii<maxbin;ii++) >> 3005 { >> 3006 vals[ii] = vals[ii]/sum; >> 3007 IPDFXBiasH.InsertValues(bins[ii], vals[ii]); >> 3008 } >> 3009 // Make IPDFXBias = true >> 3010 IPDFXBias = true; >> 3011 } >> 3012 // IPDF has been create so carry on >> 3013 G4double rndm = G4UniformRand(); >> 3014 >> 3015 // Calculate the weighting: Find the bin that the determined >> 3016 // rndm is in and the weigthing will be the difference in the >> 3017 // natural probability (from the x-axis) divided by the >> 3018 // difference in the biased probability (the area). >> 3019 size_t numberOfBin = IPDFXBiasH.GetVectorLength(); >> 3020 G4int biasn1 = 0; >> 3021 G4int biasn2 = numberOfBin/2; >> 3022 G4int biasn3 = numberOfBin - 1; >> 3023 while (biasn1 != biasn3 - 1) { >> 3024 if (rndm > IPDFXBiasH(biasn2)) >> 3025 biasn1 = biasn2; >> 3026 else >> 3027 biasn3 = biasn2; >> 3028 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 203 } 3029 } >> 3030 // retrieve the areas and then the x-axis values >> 3031 bweights[0] = IPDFXBiasH(biasn2) - IPDFXBiasH(biasn2 - 1); >> 3032 G4double xaxisl = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3033 G4double xaxisu = IPDFXBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3034 G4double NatProb = xaxisu - xaxisl; >> 3035 //G4cout << "X Bin weight " << bweights[0] << " " << rndm << G4endl; >> 3036 //G4cout << "lower and upper xaxis vals "<<xaxisl<<" "<<xaxisu<<G4endl; >> 3037 bweights[0] = NatProb/bweights[0]; >> 3038 if(verbosityLevel >= 1) >> 3039 G4cout << "X bin weight " << bweights[0] << " " << rndm << G4endl; >> 3040 return(IPDFXBiasH.GetEnergy(rndm)); >> 3041 } >> 3042 } >> 3043 >> 3044 G4double G4GeneralParticleSource::GenRandY() >> 3045 { >> 3046 if(verbosityLevel >= 1) >> 3047 G4cout << "In GenRandY" << G4endl; >> 3048 if(YBias == false) >> 3049 { >> 3050 // Y is not biased 204 G4double rndm = G4UniformRand(); 3051 G4double rndm = G4UniformRand(); 205 G4int i = 0 ; << 3052 return(rndm); 206 if (! GPSData->GetFlatSampling() ) << 3053 } 207 { << 3054 else 208 while ( rndm > GPSData->GetSourceProba << 3055 { 209 currentSource = GPSData->GetCurrentS << 3056 // Y is biased >> 3057 if(IPDFYBias == false) >> 3058 { >> 3059 // IPDF has not been created, so create it >> 3060 G4double bins[256],vals[256], sum; >> 3061 G4int ii; >> 3062 G4int maxbin = G4int(YBiasH.GetVectorLength()); >> 3063 bins[0] = YBiasH.GetLowEdgeEnergy(size_t(0)); >> 3064 vals[0] = YBiasH(size_t(0)); >> 3065 sum = vals[0]; >> 3066 for(ii=1;ii<maxbin;ii++) >> 3067 { >> 3068 bins[ii] = YBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3069 vals[ii] = YBiasH(size_t(ii)) + vals[ii-1]; >> 3070 sum = sum + YBiasH(size_t(ii)); >> 3071 } >> 3072 >> 3073 for(ii=0;ii<maxbin;ii++) >> 3074 { >> 3075 vals[ii] = vals[ii]/sum; >> 3076 IPDFYBiasH.InsertValues(bins[ii], vals[ii]); >> 3077 } >> 3078 // Make IPDFYBias = true >> 3079 IPDFYBias = true; >> 3080 } >> 3081 // IPDF has been create so carry on >> 3082 G4double rndm = G4UniformRand(); >> 3083 size_t numberOfBin = IPDFYBiasH.GetVectorLength(); >> 3084 G4int biasn1 = 0; >> 3085 G4int biasn2 = numberOfBin/2; >> 3086 G4int biasn3 = numberOfBin - 1; >> 3087 while (biasn1 != biasn3 - 1) { >> 3088 if (rndm > IPDFYBiasH(biasn2)) >> 3089 biasn1 = biasn2; >> 3090 else >> 3091 biasn3 = biasn2; >> 3092 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 210 } 3093 } >> 3094 bweights[1] = IPDFYBiasH(biasn2) - IPDFYBiasH(biasn2 - 1); >> 3095 G4double xaxisl = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3096 G4double xaxisu = IPDFYBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3097 G4double NatProb = xaxisu - xaxisl; >> 3098 bweights[1] = NatProb/bweights[1]; >> 3099 if(verbosityLevel >= 1) >> 3100 G4cout << "Y bin weight " << bweights[1] << " " << rndm << G4endl; >> 3101 return(IPDFYBiasH.GetEnergy(rndm)); >> 3102 } >> 3103 } >> 3104 >> 3105 G4double G4GeneralParticleSource::GenRandZ() >> 3106 { >> 3107 if(verbosityLevel >= 1) >> 3108 G4cout << "In GenRandZ" << G4endl; >> 3109 if(ZBias == false) >> 3110 { >> 3111 // Z is not biased >> 3112 G4double rndm = G4UniformRand(); >> 3113 return(rndm); >> 3114 } >> 3115 else >> 3116 { >> 3117 // Z is biased >> 3118 if(IPDFZBias == false) >> 3119 { >> 3120 // IPDF has not been created, so create it >> 3121 G4double bins[256],vals[256], sum; >> 3122 G4int ii; >> 3123 G4int maxbin = G4int(ZBiasH.GetVectorLength()); >> 3124 bins[0] = ZBiasH.GetLowEdgeEnergy(size_t(0)); >> 3125 vals[0] = ZBiasH(size_t(0)); >> 3126 sum = vals[0]; >> 3127 for(ii=1;ii<maxbin;ii++) >> 3128 { >> 3129 bins[ii] = ZBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3130 vals[ii] = ZBiasH(size_t(ii)) + vals[ii-1]; >> 3131 sum = sum + ZBiasH(size_t(ii)); >> 3132 } >> 3133 >> 3134 for(ii=0;ii<maxbin;ii++) >> 3135 { >> 3136 vals[ii] = vals[ii]/sum; >> 3137 IPDFZBiasH.InsertValues(bins[ii], vals[ii]); >> 3138 } >> 3139 // Make IPDFZBias = true >> 3140 IPDFZBias = true; >> 3141 } >> 3142 // IPDF has been create so carry on >> 3143 G4double rndm = G4UniformRand(); >> 3144 // size_t weight_bin_no = IPDFZBiasH.FindValueBinLocation(rndm); >> 3145 size_t numberOfBin = IPDFZBiasH.GetVectorLength(); >> 3146 G4int biasn1 = 0; >> 3147 G4int biasn2 = numberOfBin/2; >> 3148 G4int biasn3 = numberOfBin - 1; >> 3149 while (biasn1 != biasn3 - 1) { >> 3150 if (rndm > IPDFZBiasH(biasn2)) >> 3151 biasn1 = biasn2; 211 else 3152 else 212 { << 3153 biasn3 = biasn2; 213 i = G4int (GPSData->GetIntensityVector << 3154 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 214 currentSource = GPSData->GetCurrentSou << 3155 } >> 3156 bweights[2] = IPDFZBiasH(biasn2) - IPDFZBiasH(biasn2 - 1); >> 3157 G4double xaxisl = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3158 G4double xaxisu = IPDFZBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3159 G4double NatProb = xaxisu - xaxisl; >> 3160 bweights[2] = NatProb/bweights[2]; >> 3161 if(verbosityLevel >= 1) >> 3162 G4cout << "Z bin weight " << bweights[2] << " " << rndm << G4endl; >> 3163 return(IPDFZBiasH.GetEnergy(rndm)); >> 3164 } >> 3165 } >> 3166 >> 3167 G4double G4GeneralParticleSource::GenRandTheta() >> 3168 { >> 3169 if(verbosityLevel >= 1) >> 3170 { >> 3171 G4cout << "In GenRandTheta" << G4endl; >> 3172 G4cout << "Verbosity " << verbosityLevel << G4endl; >> 3173 } >> 3174 if(ThetaBias == false) >> 3175 { >> 3176 // Theta is not biased >> 3177 G4double rndm = G4UniformRand(); >> 3178 return(rndm); >> 3179 } >> 3180 else >> 3181 { >> 3182 // Theta is biased >> 3183 if(IPDFThetaBias == false) >> 3184 { >> 3185 // IPDF has not been created, so create it >> 3186 G4double bins[256],vals[256], sum; >> 3187 G4int ii; >> 3188 G4int maxbin = G4int(ThetaBiasH.GetVectorLength()); >> 3189 bins[0] = ThetaBiasH.GetLowEdgeEnergy(size_t(0)); >> 3190 vals[0] = ThetaBiasH(size_t(0)); >> 3191 sum = vals[0]; >> 3192 for(ii=1;ii<maxbin;ii++) >> 3193 { >> 3194 bins[ii] = ThetaBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3195 vals[ii] = ThetaBiasH(size_t(ii)) + vals[ii-1]; >> 3196 sum = sum + ThetaBiasH(size_t(ii)); >> 3197 } >> 3198 >> 3199 for(ii=0;ii<maxbin;ii++) >> 3200 { >> 3201 vals[ii] = vals[ii]/sum; >> 3202 IPDFThetaBiasH.InsertValues(bins[ii], vals[ii]); >> 3203 } >> 3204 // Make IPDFThetaBias = true >> 3205 IPDFThetaBias = true; >> 3206 } >> 3207 // IPDF has been create so carry on >> 3208 G4double rndm = G4UniformRand(); >> 3209 // size_t weight_bin_no = IPDFThetaBiasH.FindValueBinLocation(rndm); >> 3210 size_t numberOfBin = IPDFThetaBiasH.GetVectorLength(); >> 3211 G4int biasn1 = 0; >> 3212 G4int biasn2 = numberOfBin/2; >> 3213 G4int biasn3 = numberOfBin - 1; >> 3214 while (biasn1 != biasn3 - 1) { >> 3215 if (rndm > IPDFThetaBiasH(biasn2)) >> 3216 biasn1 = biasn2; >> 3217 else >> 3218 biasn3 = biasn2; >> 3219 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; 215 } 3220 } >> 3221 bweights[3] = IPDFThetaBiasH(biasn2) - IPDFThetaBiasH(biasn2 - 1); >> 3222 G4double xaxisl = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3223 G4double xaxisu = IPDFThetaBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3224 G4double NatProb = xaxisu - xaxisl; >> 3225 bweights[3] = NatProb/bweights[3]; >> 3226 if(verbosityLevel >= 1) >> 3227 G4cout << "Theta bin weight " << bweights[3] << " " << rndm << G4endl; >> 3228 return(IPDFThetaBiasH.GetEnergy(rndm)); >> 3229 } >> 3230 } >> 3231 >> 3232 G4double G4GeneralParticleSource::GenRandPhi() >> 3233 { >> 3234 if(verbosityLevel >= 1) >> 3235 G4cout << "In GenRandPhi" << G4endl; >> 3236 if(PhiBias == false) >> 3237 { >> 3238 // Phi is not biased >> 3239 G4double rndm = G4UniformRand(); >> 3240 return(rndm); 216 } 3241 } 217 currentSource->GeneratePrimaryVertex(evt); << 218 } << 219 else 3242 else 220 { << 221 for (G4int i = 0; i < GPSData->GetIntensi << 222 { 3243 { 223 GPSData->GetCurrentSource(i)->GeneratePr << 3244 // Phi is biased >> 3245 if(IPDFPhiBias == false) >> 3246 { >> 3247 // IPDF has not been created, so create it >> 3248 G4double bins[256],vals[256], sum; >> 3249 G4int ii; >> 3250 G4int maxbin = G4int(PhiBiasH.GetVectorLength()); >> 3251 bins[0] = PhiBiasH.GetLowEdgeEnergy(size_t(0)); >> 3252 vals[0] = PhiBiasH(size_t(0)); >> 3253 sum = vals[0]; >> 3254 for(ii=1;ii<maxbin;ii++) >> 3255 { >> 3256 bins[ii] = PhiBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3257 vals[ii] = PhiBiasH(size_t(ii)) + vals[ii-1]; >> 3258 sum = sum + PhiBiasH(size_t(ii)); >> 3259 } >> 3260 >> 3261 for(ii=0;ii<maxbin;ii++) >> 3262 { >> 3263 vals[ii] = vals[ii]/sum; >> 3264 IPDFPhiBiasH.InsertValues(bins[ii], vals[ii]); >> 3265 } >> 3266 // Make IPDFPhiBias = true >> 3267 IPDFPhiBias = true; >> 3268 } >> 3269 // IPDF has been create so carry on >> 3270 G4double rndm = G4UniformRand(); >> 3271 // size_t weight_bin_no = IPDFPhiBiasH.FindValueBinLocation(rndm); >> 3272 size_t numberOfBin = IPDFPhiBiasH.GetVectorLength(); >> 3273 G4int biasn1 = 0; >> 3274 G4int biasn2 = numberOfBin/2; >> 3275 G4int biasn3 = numberOfBin - 1; >> 3276 while (biasn1 != biasn3 - 1) { >> 3277 if (rndm > IPDFPhiBiasH(biasn2)) >> 3278 biasn1 = biasn2; >> 3279 else >> 3280 biasn3 = biasn2; >> 3281 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; >> 3282 } >> 3283 bweights[4] = IPDFPhiBiasH(biasn2) - IPDFPhiBiasH(biasn2 - 1); >> 3284 G4double xaxisl = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3285 G4double xaxisu = IPDFPhiBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3286 G4double NatProb = xaxisu - xaxisl; >> 3287 bweights[4] = NatProb/bweights[4]; >> 3288 if(verbosityLevel >= 1) >> 3289 G4cout << "Phi bin weight " << bweights[4] << " " << rndm << G4endl; >> 3290 return(IPDFPhiBiasH.GetEnergy(rndm)); 224 } 3291 } >> 3292 } >> 3293 >> 3294 G4double G4GeneralParticleSource::GenRandEnergy() >> 3295 { >> 3296 if(verbosityLevel >= 1) >> 3297 G4cout << "In GenRandEnergy" << G4endl; >> 3298 if(EnergyBias == false) >> 3299 { >> 3300 // Energy is not biased >> 3301 G4double rndm = G4UniformRand(); >> 3302 return(rndm); >> 3303 } >> 3304 else >> 3305 { >> 3306 // ENERGY is biased >> 3307 if(IPDFEnergyBias == false) >> 3308 { >> 3309 // IPDF has not been created, so create it >> 3310 G4double bins[256],vals[256], sum; >> 3311 G4int ii; >> 3312 G4int maxbin = G4int(EnergyBiasH.GetVectorLength()); >> 3313 bins[0] = EnergyBiasH.GetLowEdgeEnergy(size_t(0)); >> 3314 vals[0] = EnergyBiasH(size_t(0)); >> 3315 sum = vals[0]; >> 3316 for(ii=1;ii<maxbin;ii++) >> 3317 { >> 3318 bins[ii] = EnergyBiasH.GetLowEdgeEnergy(size_t(ii)); >> 3319 vals[ii] = EnergyBiasH(size_t(ii)) + vals[ii-1]; >> 3320 sum = sum + EnergyBiasH(size_t(ii)); >> 3321 } >> 3322 >> 3323 for(ii=0;ii<maxbin;ii++) >> 3324 { >> 3325 vals[ii] = vals[ii]/sum; >> 3326 IPDFEnergyBiasH.InsertValues(bins[ii], vals[ii]); >> 3327 } >> 3328 // Make IPDFEnergyBias = true >> 3329 IPDFEnergyBias = true; >> 3330 } >> 3331 // IPDF has been create so carry on >> 3332 G4double rndm = G4UniformRand(); >> 3333 // size_t weight_bin_no = IPDFEnergyBiasH.FindValueBinLocation(rndm); >> 3334 size_t numberOfBin = IPDFEnergyBiasH.GetVectorLength(); >> 3335 G4int biasn1 = 0; >> 3336 G4int biasn2 = numberOfBin/2; >> 3337 G4int biasn3 = numberOfBin - 1; >> 3338 while (biasn1 != biasn3 - 1) { >> 3339 if (rndm > IPDFEnergyBiasH(biasn2)) >> 3340 biasn1 = biasn2; >> 3341 else >> 3342 biasn3 = biasn2; >> 3343 biasn2 = biasn1 + (biasn3 - biasn1 + 1)/2; >> 3344 } >> 3345 bweights[5] = IPDFEnergyBiasH(biasn2) - IPDFEnergyBiasH(biasn2 - 1); >> 3346 G4double xaxisl = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2-1)); >> 3347 G4double xaxisu = IPDFEnergyBiasH.GetLowEdgeEnergy(size_t(biasn2)); >> 3348 G4double NatProb = xaxisu - xaxisl; >> 3349 bweights[5] = NatProb/bweights[5]; >> 3350 if(verbosityLevel >= 1) >> 3351 G4cout << "Energy bin weight " << bweights[5] << " " << rndm << G4endl; >> 3352 return(IPDFEnergyBiasH.GetEnergy(rndm)); >> 3353 } >> 3354 } >> 3355 >> 3356 // Verbosity >> 3357 void G4GeneralParticleSource::SetVerbosity(int vL) >> 3358 { >> 3359 verbosityLevel = vL; >> 3360 G4cout << "Verbosity Set to: " << verbosityLevel << G4endl; >> 3361 } >> 3362 >> 3363 void G4GeneralParticleSource::SetParticleDefinition >> 3364 (G4ParticleDefinition* aParticleDefinition) >> 3365 { >> 3366 particle_definition = aParticleDefinition; >> 3367 particle_charge = particle_definition->GetPDGCharge(); >> 3368 } >> 3369 >> 3370 // SR1.3 >> 3371 //void G4GeneralParticleSource::SetNucleus (Nucleus theIon1) >> 3372 //{ >> 3373 // theIon = theIon1; >> 3374 >> 3375 // G4IonTable *theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable()); >> 3376 // G4ParticleDefinition *aIon = NULL; >> 3377 >> 3378 // G4int A = theIon.GetA(); >> 3379 // G4int Z = theIon.GetZ(); >> 3380 // G4double E = theIon.GetE(); >> 3381 >> 3382 // aIon = theIonTable->GetIon (Z, A, E); >> 3383 >> 3384 // SetParticleDefinition(aIon); >> 3385 //} >> 3386 >> 3387 >> 3388 void G4GeneralParticleSource::GeneratePrimaryVertex(G4Event *evt) >> 3389 { >> 3390 if(particle_definition==NULL) return; >> 3391 >> 3392 // Position stuff >> 3393 G4bool srcconf = false; >> 3394 G4int LoopCount = 0; >> 3395 while(srcconf == false) >> 3396 { >> 3397 if(SourcePosType == "Point") >> 3398 GeneratePointSource(); >> 3399 else if(SourcePosType == "Beam") >> 3400 GeneratePointsInBeam(); >> 3401 else if(SourcePosType == "Plane") >> 3402 GeneratePointsInPlane(); >> 3403 else if(SourcePosType == "Surface") >> 3404 GeneratePointsOnSurface(); >> 3405 else if(SourcePosType == "Volume") >> 3406 GeneratePointsInVolume(); >> 3407 else >> 3408 { >> 3409 G4cout << "Error: SourcePosType undefined" << G4endl; >> 3410 G4cout << "Generating point source" << G4endl; >> 3411 GeneratePointSource(); >> 3412 } >> 3413 if(Confine == true) >> 3414 { >> 3415 srcconf = IsSourceConfined(); >> 3416 // if source in confined srcconf = true terminating the loop >> 3417 // if source isnt confined srcconf = false and loop continues >> 3418 } >> 3419 else if(Confine == false) >> 3420 srcconf = true; // terminate loop >> 3421 >> 3422 LoopCount++; >> 3423 if(LoopCount == 100000) >> 3424 { >> 3425 G4cout << "*************************************" << G4endl; >> 3426 G4cout << "LoopCount = 100000" << G4endl; >> 3427 G4cout << "Either the source distribution >> confinement" << G4endl; >> 3428 G4cout << "or any confining volume may not overlap with" << G4endl; >> 3429 G4cout << "the source distribution or any confining volumes" << G4endl; >> 3430 G4cout << "may not exist"<< G4endl; >> 3431 G4cout << "If you have set confine then this will be ignored" <<G4endl; >> 3432 G4cout << "for this event." << G4endl; >> 3433 G4cout << "*************************************" << G4endl; >> 3434 srcconf = true; //Avoids an infinite loop >> 3435 } >> 3436 } >> 3437 // Angular stuff >> 3438 if(AngDistType == "iso") >> 3439 GenerateIsotropicFlux(); >> 3440 else if(AngDistType == "cos") >> 3441 GenerateCosineLawFlux(); >> 3442 else if(AngDistType == "planar") >> 3443 GeneratePlanarFlux(); >> 3444 else if(AngDistType == "beam1d" || AngDistType == "beam2d" ) >> 3445 GenerateBeamFlux(); >> 3446 else if(AngDistType == "user") >> 3447 GenerateUserDefFlux(); >> 3448 else >> 3449 G4cout << "Error: AngDistType has unusual value" << G4endl; >> 3450 // Energy stuff >> 3451 if(EnergyDisType == "Mono") >> 3452 GenerateMonoEnergetic(); >> 3453 else if(EnergyDisType == "Lin") >> 3454 GenerateLinearEnergies(); >> 3455 else if(EnergyDisType == "Pow") >> 3456 GeneratePowEnergies(); >> 3457 else if(EnergyDisType == "Exp") >> 3458 GenerateExpEnergies(); >> 3459 else if(EnergyDisType == "Gauss") >> 3460 GenerateGaussEnergies(); >> 3461 else if(EnergyDisType == "Brem") >> 3462 GenerateBremEnergies(); >> 3463 else if(EnergyDisType == "Bbody") >> 3464 GenerateBbodyEnergies(); >> 3465 else if(EnergyDisType == "Cdg") >> 3466 GenerateCdgEnergies(); >> 3467 else if(EnergyDisType == "User") >> 3468 GenUserHistEnergies(); >> 3469 else if(EnergyDisType == "Arb") >> 3470 GenArbPointEnergies(); >> 3471 else if(EnergyDisType == "Epn") >> 3472 GenEpnHistEnergies(); >> 3473 else >> 3474 G4cout << "Error: EnergyDisType has unusual value" << G4endl; >> 3475 >> 3476 // create a new vertex >> 3477 G4PrimaryVertex* vertex = >> 3478 new G4PrimaryVertex(particle_position,particle_time); >> 3479 >> 3480 if(verbosityLevel == 2) >> 3481 G4cout << "Creating primaries and assigning to vertex" << G4endl; >> 3482 // create new primaries and set them to the vertex >> 3483 G4double mass = particle_definition->GetPDGMass(); >> 3484 G4double energy = particle_energy + mass; >> 3485 G4double pmom = sqrt(energy*energy-mass*mass); >> 3486 G4double px = pmom*particle_momentum_direction.x(); >> 3487 G4double py = pmom*particle_momentum_direction.y(); >> 3488 G4double pz = pmom*particle_momentum_direction.z(); >> 3489 >> 3490 if(verbosityLevel > 1){ >> 3491 G4cout << "Particle name: "<<particle_definition->GetParticleName() << G4endl; >> 3492 G4cout << " Energy: "<<particle_energy << G4endl; >> 3493 G4cout << " Position: "<<particle_position<< G4endl; >> 3494 G4cout << " Direction: "<<particle_momentum_direction << G4endl; >> 3495 G4cout << " NumberOfParticlesToBeGenerated: "<<NumberOfParticlesToBeGenerated << G4endl; 225 } 3496 } >> 3497 for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) >> 3498 { >> 3499 G4PrimaryParticle* particle = >> 3500 new G4PrimaryParticle(particle_definition,px,py,pz); >> 3501 particle->SetMass( mass ); >> 3502 particle->SetCharge( particle_charge ); >> 3503 particle->SetPolarization(particle_polarization.x(), >> 3504 particle_polarization.y(), >> 3505 particle_polarization.z()); >> 3506 vertex->SetPrimary( particle ); >> 3507 >> 3508 // Set bweight equal to the multiple of all non-zero weights >> 3509 bweight = 1.; >> 3510 for(int bib=0; bib<6; bib++) >> 3511 { >> 3512 if(bweights[bib] > 0.) bweight = bweight * bweights[bib]; >> 3513 } >> 3514 // bweight will now contain the events final weighting. >> 3515 // now pass it to the primary vertex >> 3516 vertex->SetWeight(bweight); >> 3517 } >> 3518 evt->AddPrimaryVertex( vertex ); >> 3519 if(verbosityLevel > 1) >> 3520 G4cout << " Primary Vetex generated !"<< G4endl; 226 } 3521 } >> 3522 >> 3523 >> 3524 >> 3525 >> 3526 >> 3527 >> 3528 >> 3529 >> 3530 227 3531