| 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-04-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