Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/underground_physics/src/DMXParticleSource.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/underground_physics/src/DMXParticleSource.cc (Version 11.3.0) and /examples/advanced/underground_physics/src/DMXParticleSource.cc (Version 11.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 // -------------------------------------------     27 // --------------------------------------------------------------
 28 //   GEANT 4 - Underground Dark Matter Detecto     28 //   GEANT 4 - Underground Dark Matter Detector Advanced Example
 29 //                                                 29 //
 30 //      For information related to this code c     30 //      For information related to this code contact: Alex Howard
 31 //      e-mail: alexander.howard@cern.ch           31 //      e-mail: alexander.howard@cern.ch
 32 // -------------------------------------------     32 // --------------------------------------------------------------
 33 // Comments                                        33 // Comments
 34 //                                                 34 //
 35 //                  Underground Advanced           35 //                  Underground Advanced
 36 //               by A. Howard and H. Araujo        36 //               by A. Howard and H. Araujo 
 37 //                    (27th November 2001)         37 //                    (27th November 2001)
 38 //                                                 38 //
 39 //                                                 39 //
 40 // ParticleSource program                          40 // ParticleSource program
 41 // -------------------------------------------     41 // --------------------------------------------------------------
 42 //////////////////////////////////////////////     42 //////////////////////////////////////////////////////////////////////////////
 43 // This particle source is a shortened version     43 // This particle source is a shortened version of G4GeneralParticleSource by
 44 // C Ferguson, F Lei & P Truscott (University      44 // C Ferguson, F Lei & P Truscott (University of Southampton / DERA), with
 45 // some minor modifications.                       45 // some minor modifications.
 46 //////////////////////////////////////////////     46 //////////////////////////////////////////////////////////////////////////////
 47                                                    47 
 48 #include <cmath>                                   48 #include <cmath>
 49                                                    49 
 50 #include "DMXParticleSource.hh"                    50 #include "DMXParticleSource.hh"
 51                                                    51 
 52 #include "G4PhysicalConstants.hh"                  52 #include "G4PhysicalConstants.hh"
 53 #include "G4SystemOfUnits.hh"                      53 #include "G4SystemOfUnits.hh"
 54 #include "G4PrimaryParticle.hh"                    54 #include "G4PrimaryParticle.hh"
 55 #include "G4Event.hh"                              55 #include "G4Event.hh"
 56 #include "Randomize.hh"                            56 #include "Randomize.hh"
 57 #include "G4TransportationManager.hh"              57 #include "G4TransportationManager.hh"
 58 #include "G4VPhysicalVolume.hh"                    58 #include "G4VPhysicalVolume.hh"
 59 #include "G4PhysicalVolumeStore.hh"                59 #include "G4PhysicalVolumeStore.hh"
 60 #include "G4ParticleTable.hh"                      60 #include "G4ParticleTable.hh"
 61 #include "G4ParticleDefinition.hh"                 61 #include "G4ParticleDefinition.hh"
 62 #include "G4IonTable.hh"                           62 #include "G4IonTable.hh"
 63 #include "G4Ions.hh"                               63 #include "G4Ions.hh"
 64 #include "G4TrackingManager.hh"                    64 #include "G4TrackingManager.hh"
 65 #include "G4Track.hh"                              65 #include "G4Track.hh"
 66                                                    66 
 67                                                    67 
 68 DMXParticleSource::DMXParticleSource() {           68 DMXParticleSource::DMXParticleSource() {
 69                                                    69 
 70   NumberOfParticlesToBeGenerated = 1;              70   NumberOfParticlesToBeGenerated = 1;
 71   particle_definition = nullptr;                   71   particle_definition = nullptr;
 72   G4ThreeVector zero(0., 0., 0.);                  72   G4ThreeVector zero(0., 0., 0.);
 73   particle_momentum_direction = G4ParticleMome     73   particle_momentum_direction = G4ParticleMomentum(1., 0., 0.);
 74   particle_energy = 1.0*MeV;                       74   particle_energy = 1.0*MeV;
 75   particle_position = zero;                        75   particle_position = zero;
 76   particle_time = 0.0;                             76   particle_time = 0.0;
 77   particle_polarization = zero;                    77   particle_polarization = zero;
 78   particle_charge = 0.0;                           78   particle_charge = 0.0;
 79                                                    79 
 80   SourcePosType = "Volume";                        80   SourcePosType = "Volume";
 81   Shape = "NULL";                                  81   Shape = "NULL";
 82   halfz = 0.;                                      82   halfz = 0.;
 83   Radius = 0.;                                     83   Radius = 0.;
 84   CentreCoords = zero;                             84   CentreCoords = zero;
 85   Confine = false;                                 85   Confine = false;
 86   VolName = "NULL";                                86   VolName = "NULL";
 87                                                    87 
 88   AngDistType = "iso";                             88   AngDistType = "iso"; 
 89   MinTheta = 0.;                                   89   MinTheta = 0.;
 90   MaxTheta = pi;                                   90   MaxTheta = pi;
 91   MinPhi = 0.;                                     91   MinPhi = 0.;
 92   MaxPhi = twopi;                                  92   MaxPhi = twopi;
 93                                                    93 
 94   EnergyDisType = "Mono";                          94   EnergyDisType = "Mono";
 95   MonoEnergy = 1*MeV;                              95   MonoEnergy = 1*MeV;
 96                                                    96 
 97   verbosityLevel = 0;                              97   verbosityLevel = 0;
 98                                                    98 
 99   theMessenger = new DMXParticleSourceMessenge     99   theMessenger = new DMXParticleSourceMessenger(this);
100   gNavigator = G4TransportationManager::GetTra    100   gNavigator = G4TransportationManager::GetTransportationManager()
101     ->GetNavigatorForTracking();                  101     ->GetNavigatorForTracking();
102 }                                                 102 }
103                                                   103 
104 DMXParticleSource::~DMXParticleSource()           104 DMXParticleSource::~DMXParticleSource()
105 {                                                 105 {
106   delete theMessenger;                            106   delete theMessenger;
107 }                                                 107 }
108                                                   108 
109 void DMXParticleSource::SetPosDisType(G4String    109 void DMXParticleSource::SetPosDisType(G4String PosType) 
110 {                                                 110 {
111   SourcePosType = PosType;                        111   SourcePosType = PosType;
112 }                                                 112 }
113                                                   113 
114 void DMXParticleSource::SetPosDisShape(G4Strin    114 void DMXParticleSource::SetPosDisShape(G4String shapeType)
115 {                                                 115 {
116   Shape = shapeType;                              116   Shape = shapeType;
117 }                                                 117 }
118                                                   118 
119 void DMXParticleSource::SetCentreCoords(G4Thre    119 void DMXParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre)
120 {                                                 120 {
121   CentreCoords = coordsOfCentre;                  121   CentreCoords = coordsOfCentre;
122 }                                                 122 }
123                                                   123 
124 void DMXParticleSource::SetHalfZ(G4double zhal    124 void DMXParticleSource::SetHalfZ(G4double zhalf)
125 {                                                 125 {
126   halfz = zhalf;                                  126   halfz = zhalf;
127 }                                                 127 }
128                                                   128 
129 void DMXParticleSource::SetRadius(G4double rad    129 void DMXParticleSource::SetRadius(G4double radius)
130 {                                                 130 {
131   Radius = radius;                                131   Radius = radius;
132 }                                                 132 }
133                                                   133 
134 void DMXParticleSource::ConfineSourceToVolume(    134 void DMXParticleSource::ConfineSourceToVolume(G4String Vname)
135 {                                                 135 {
136   VolName = Vname;                                136   VolName = Vname;
137   if(verbosityLevel == 2) G4cout << VolName <<    137   if(verbosityLevel == 2) G4cout << VolName << G4endl;
138                                                   138 
139   // checks if selected volume exists             139   // checks if selected volume exists
140   G4VPhysicalVolume *tempPV      = nullptr;       140   G4VPhysicalVolume *tempPV      = nullptr;
141   G4PhysicalVolumeStore *PVStore = nullptr;       141   G4PhysicalVolumeStore *PVStore = nullptr;
142   G4String theRequiredVolumeName = VolName;       142   G4String theRequiredVolumeName = VolName;
143   PVStore = G4PhysicalVolumeStore::GetInstance    143   PVStore = G4PhysicalVolumeStore::GetInstance();
144   G4bool found = false;                           144   G4bool found = false;
145   if(verbosityLevel == 2) G4cout << PVStore->s    145   if(verbosityLevel == 2) G4cout << PVStore->size() << G4endl;
146                                                   146 
147   tempPV = PVStore->GetVolume(theRequiredVolum    147   tempPV = PVStore->GetVolume(theRequiredVolumeName, false);
148   if (tempPV != nullptr) { found = true; }        148   if (tempPV != nullptr) { found = true; }
149                                                   149 
150   // found = true then the volume exists else     150   // found = true then the volume exists else it doesnt.
151   if(found == true) {                             151   if(found == true) {
152     if(verbosityLevel >= 1)                       152     if(verbosityLevel >= 1)
153       G4cout << "Volume " << VolName << " exis    153       G4cout << "Volume " << VolName << " exists" << G4endl;
154     Confine = true;                               154     Confine = true;
155   }                                               155   }
156   else if(VolName=="NULL")                        156   else if(VolName=="NULL")
157     Confine = false;                              157     Confine = false;
158   else {                                          158   else {
159     G4cout << " **** Error: Volume does not ex    159     G4cout << " **** Error: Volume does not exist **** " << G4endl;
160     G4cout << " Ignoring confine condition" <<    160     G4cout << " Ignoring confine condition" << G4endl;
161     VolName = "NULL";                             161     VolName = "NULL";
162     Confine = false;                              162     Confine = false;
163   }                                               163   }
164 }                                                 164 }
165                                                   165 
166 void DMXParticleSource::SetAngDistType(G4Strin    166 void DMXParticleSource::SetAngDistType(G4String atype)
167 {                                                 167 {
168   AngDistType = atype;                            168   AngDistType = atype;
169 }                                                 169 }
170                                                   170 
171 void DMXParticleSource::GeneratePointSource()     171 void DMXParticleSource::GeneratePointSource()
172 {                                                 172 {
173   // Generates Points given the point source.     173   // Generates Points given the point source.
174   if(SourcePosType == "Point")                    174   if(SourcePosType == "Point")
175     particle_position = CentreCoords;             175     particle_position = CentreCoords;
176   else                                            176   else
177     if(verbosityLevel >= 1)                       177     if(verbosityLevel >= 1)
178       G4cout << "Error SourcePosType is not se    178       G4cout << "Error SourcePosType is not set to Point" << G4endl;
179 }                                                 179 }
180                                                   180 
181 void DMXParticleSource::GeneratePointsInVolume    181 void DMXParticleSource::GeneratePointsInVolume()
182 {                                                 182 {
183   G4ThreeVector RandPos;                          183   G4ThreeVector RandPos;
184   G4double x=0., y=0., z=0.;                      184   G4double x=0., y=0., z=0.;
185                                                   185   
186   if(SourcePosType != "Volume" && verbosityLev    186   if(SourcePosType != "Volume" && verbosityLevel >= 1)
187     G4cout << "Error SourcePosType not Volume"    187     G4cout << "Error SourcePosType not Volume" << G4endl;
188                                                   188   
189   if(Shape == "Sphere") {                         189   if(Shape == "Sphere") {
190     x = Radius*2.;                                190     x = Radius*2.;
191     y = Radius*2.;                                191     y = Radius*2.;
192     z = Radius*2.;                                192     z = Radius*2.;
193     while(((x*x)+(y*y)+(z*z)) > (Radius*Radius    193     while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) {
194       x = G4UniformRand();                        194       x = G4UniformRand();
195       y = G4UniformRand();                        195       y = G4UniformRand();
196       z = G4UniformRand();                        196       z = G4UniformRand();
197                                                   197       
198       x = (x*2.*Radius) - Radius;                 198       x = (x*2.*Radius) - Radius;
199       y = (y*2.*Radius) - Radius;                 199       y = (y*2.*Radius) - Radius;
200       z = (z*2.*Radius) - Radius;                 200       z = (z*2.*Radius) - Radius;
201     }                                             201     }
202   }                                               202   }
203                                                   203 
204   else if(Shape == "Cylinder") {                  204   else if(Shape == "Cylinder") {
205     x = Radius*2.;                                205     x = Radius*2.;
206     y = Radius*2.;                                206     y = Radius*2.;
207     while(((x*x)+(y*y)) > (Radius*Radius)) {      207     while(((x*x)+(y*y)) > (Radius*Radius)) {
208       x = G4UniformRand();                        208       x = G4UniformRand();
209       y = G4UniformRand();                        209       y = G4UniformRand();
210       z = G4UniformRand();                        210       z = G4UniformRand();
211       x = (x*2.*Radius) - Radius;                 211       x = (x*2.*Radius) - Radius;
212       y = (y*2.*Radius) - Radius;                 212       y = (y*2.*Radius) - Radius;
213       z = (z*2.*halfz) - halfz;                   213       z = (z*2.*halfz) - halfz;
214     }                                             214     }
215   }                                               215   }
216                                                   216   
217   else                                            217   else
218     G4cout << "Error: Volume Shape Does Not Ex    218     G4cout << "Error: Volume Shape Does Not Exist" << G4endl;
219                                                   219 
220   RandPos.setX(x);                                220   RandPos.setX(x);
221   RandPos.setY(y);                                221   RandPos.setY(y);
222   RandPos.setZ(z);                                222   RandPos.setZ(z);
223   particle_position = CentreCoords + RandPos;     223   particle_position = CentreCoords + RandPos;
224                                                   224 
225 }                                                 225 }
226                                                   226 
227 G4bool DMXParticleSource::IsSourceConfined()      227 G4bool DMXParticleSource::IsSourceConfined()
228 {                                                 228 {
229                                                   229 
230   // Method to check point is within the volum    230   // Method to check point is within the volume specified
231   if(Confine == false)                            231   if(Confine == false)
232     G4cout << "Error: Confine is false" << G4e    232     G4cout << "Error: Confine is false" << G4endl;
233   G4ThreeVector null_vec(0.,0.,0.);               233   G4ThreeVector null_vec(0.,0.,0.);
234   G4ThreeVector *ptr = &null_vec;                 234   G4ThreeVector *ptr = &null_vec;
235                                                   235 
236   // Check particle_position is within VolName    236   // Check particle_position is within VolName
237   G4VPhysicalVolume *theVolume;                   237   G4VPhysicalVolume *theVolume;
238   theVolume=gNavigator->LocateGlobalPointAndSe    238   theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
239   G4String theVolName = theVolume->GetName();     239   G4String theVolName = theVolume->GetName();
240   if(theVolName == VolName) {                     240   if(theVolName == VolName) {
241     if(verbosityLevel >= 1)                       241     if(verbosityLevel >= 1)
242       G4cout << "Particle is in volume " << Vo    242       G4cout << "Particle is in volume " << VolName << G4endl;
243     return(true);                                 243     return(true);
244   }                                               244   }
245   else                                            245   else
246     return(false);                                246     return(false);
247 }                                                 247 }
248                                                   248 
249 void DMXParticleSource::SetParticleMomentumDir    249 void DMXParticleSource::SetParticleMomentumDirection
250    (G4ParticleMomentum aDirection) {              250    (G4ParticleMomentum aDirection) {
251                                                   251 
252   particle_momentum_direction =  aDirection.un    252   particle_momentum_direction =  aDirection.unit();
253 }                                                 253 }
254                                                   254 
255 void DMXParticleSource::GenerateIsotropicFlux(    255 void DMXParticleSource::GenerateIsotropicFlux()
256 {                                                 256 {
257                                                   257 
258   G4double rndm, rndm2;                           258   G4double rndm, rndm2;
259   G4double px, py, pz;                            259   G4double px, py, pz;
260                                                   260 
261   G4double sintheta, sinphi, costheta, cosphi;    261   G4double sintheta, sinphi, costheta, cosphi;
262   rndm = G4UniformRand();                         262   rndm = G4UniformRand();
263   costheta = std::cos(MinTheta) - rndm * (std:    263   costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta)
264                                 - std::cos(Max    264                                 - std::cos(MaxTheta));
265   sintheta = std::sqrt(1. - costheta*costheta)    265   sintheta = std::sqrt(1. - costheta*costheta);
266                                                   266   
267   rndm2 = G4UniformRand();                        267   rndm2 = G4UniformRand();
268   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;       268   Phi = MinPhi + (MaxPhi - MinPhi) * rndm2; 
269   sinphi = std::sin(Phi);                         269   sinphi = std::sin(Phi);
270   cosphi = std::cos(Phi);                         270   cosphi = std::cos(Phi);
271                                                   271 
272   px = -sintheta * cosphi;                        272   px = -sintheta * cosphi;
273   py = -sintheta * sinphi;                        273   py = -sintheta * sinphi;
274   pz = -costheta;                                 274   pz = -costheta;
275                                                   275 
276   G4double ResMag = std::sqrt((px*px) + (py*py    276   G4double ResMag = std::sqrt((px*px) + (py*py) + (pz*pz));
277   px = px/ResMag;                                 277   px = px/ResMag;
278   py = py/ResMag;                                 278   py = py/ResMag;
279   pz = pz/ResMag;                                 279   pz = pz/ResMag;
280                                                   280 
281   particle_momentum_direction.setX(px);           281   particle_momentum_direction.setX(px);
282   particle_momentum_direction.setY(py);           282   particle_momentum_direction.setY(py);
283   particle_momentum_direction.setZ(pz);           283   particle_momentum_direction.setZ(pz);
284                                                   284 
285   // particle_momentum_direction now holds uni    285   // particle_momentum_direction now holds unit momentum vector.
286   if(verbosityLevel >= 2)                         286   if(verbosityLevel >= 2)
287     G4cout << "Generating isotropic vector: "     287     G4cout << "Generating isotropic vector: "
288            << particle_momentum_direction << G    288            << particle_momentum_direction << G4endl;
289 }                                                 289 }
290                                                   290 
291 void DMXParticleSource::SetEnergyDisType(G4Str    291 void DMXParticleSource::SetEnergyDisType(G4String DisType)
292 {                                                 292 {
293   EnergyDisType = DisType;                        293   EnergyDisType = DisType;
294 }                                                 294 }
295                                                   295 
296 void DMXParticleSource::SetMonoEnergy(G4double    296 void DMXParticleSource::SetMonoEnergy(G4double menergy)
297 {                                                 297 {
298   MonoEnergy = menergy;                           298   MonoEnergy = menergy;
299 }                                                 299 }
300                                                   300 
301 void DMXParticleSource::GenerateMonoEnergetic(    301 void DMXParticleSource::GenerateMonoEnergetic()
302 {                                                 302 {
303   particle_energy = MonoEnergy;                   303   particle_energy = MonoEnergy;
304 }                                                 304 }
305                                                   305 
306 void DMXParticleSource::SetVerbosity(int vL)      306 void DMXParticleSource::SetVerbosity(int vL)
307 {                                                 307 {
308   verbosityLevel = vL;                            308   verbosityLevel = vL;
309   G4cout << "Verbosity Set to: " << verbosityL    309   G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
310 }                                                 310 }
311                                                   311 
312 void DMXParticleSource::SetParticleDefinition     312 void DMXParticleSource::SetParticleDefinition
313   (G4ParticleDefinition* aParticleDefinition)     313   (G4ParticleDefinition* aParticleDefinition)
314 {                                                 314 {
315   particle_definition = aParticleDefinition;      315   particle_definition = aParticleDefinition;
316   particle_charge = particle_definition->GetPD    316   particle_charge = particle_definition->GetPDGCharge();
317 }                                                 317 }
318                                                   318 
319 void DMXParticleSource::GeneratePrimaryVertex(    319 void DMXParticleSource::GeneratePrimaryVertex(G4Event *evt)
320 {                                                 320 {
321                                                   321 
322   if(particle_definition==nullptr) {              322   if(particle_definition==nullptr) {
323     G4cout << "No particle has been defined!"     323     G4cout << "No particle has been defined!" << G4endl;
324     return;                                       324     return;
325   }                                               325   }
326                                                   326   
327   // Position                                     327   // Position
328   G4bool srcconf = false;                         328   G4bool srcconf = false;
329   G4int LoopCount = 0;                            329   G4int LoopCount = 0;
330                                                   330   
331   while(srcconf == false)  {                      331   while(srcconf == false)  {
332     if(SourcePosType == "Point")                  332     if(SourcePosType == "Point")
333       GeneratePointSource();                      333       GeneratePointSource();
334     else if(SourcePosType == "Volume")            334     else if(SourcePosType == "Volume")
335       GeneratePointsInVolume();                   335       GeneratePointsInVolume();
336     else {                                        336     else {
337       G4cout << "Error: SourcePosType undefine    337       G4cout << "Error: SourcePosType undefined" << G4endl;
338       G4cout << "Generating point source" << G    338       G4cout << "Generating point source" << G4endl;
339       GeneratePointSource();                      339       GeneratePointSource();
340     }                                             340     }
341     if(Confine == true) {                         341     if(Confine == true) {
342       srcconf = IsSourceConfined();               342       srcconf = IsSourceConfined();
343       // if source in confined srcconf = true     343       // if source in confined srcconf = true terminating the loop
344       // if source isnt confined srcconf = fal    344       // if source isnt confined srcconf = false and loop continues
345     }                                             345     }
346     else if(Confine == false)                     346     else if(Confine == false)
347       srcconf = true; // terminate loop           347       srcconf = true; // terminate loop
348                                                   348     
349     ++LoopCount;                                  349     ++LoopCount;
350     if(LoopCount == 100000) {                     350     if(LoopCount == 100000) {
351       G4cout << "*****************************    351       G4cout << "*************************************" << G4endl;
352         G4cout << "LoopCount = 100000" << G4en    352         G4cout << "LoopCount = 100000" << G4endl;
353         G4cout << "Either the source distribut    353         G4cout << "Either the source distribution >> confinement" << G4endl;
354         G4cout << "or any confining volume may    354         G4cout << "or any confining volume may not overlap with" << G4endl;
355         G4cout << "the source distribution or     355         G4cout << "the source distribution or any confining volumes" << G4endl;
356         G4cout << "may not exist"<< G4endl;       356         G4cout << "may not exist"<< G4endl;
357         G4cout << "If you have set confine the    357         G4cout << "If you have set confine then this will be ignored" <<G4endl;
358         G4cout << "for this event." << G4endl;    358         G4cout << "for this event." << G4endl;
359         G4cout << "***************************    359         G4cout << "*************************************" << G4endl;
360         srcconf = true; //Avoids an infinite l    360         srcconf = true; //Avoids an infinite loop
361       }                                           361       }
362   }                                               362   }
363                                                   363 
364   // Angular stuff                                364   // Angular stuff
365   if(AngDistType == "iso")                        365   if(AngDistType == "iso")
366     GenerateIsotropicFlux();                      366     GenerateIsotropicFlux();
367   else if(AngDistType == "direction")             367   else if(AngDistType == "direction")
368     SetParticleMomentumDirection(particle_mome    368     SetParticleMomentumDirection(particle_momentum_direction);
369   else                                            369   else
370     G4cout << "Error: AngDistType has unusual     370     G4cout << "Error: AngDistType has unusual value" << G4endl;
371   // Energy stuff                                 371   // Energy stuff
372   if(EnergyDisType == "Mono")                     372   if(EnergyDisType == "Mono")
373     GenerateMonoEnergetic();                      373     GenerateMonoEnergetic();
374   else                                            374   else
375     G4cout << "Error: EnergyDisType has unusua    375     G4cout << "Error: EnergyDisType has unusual value" << G4endl;
376                                                   376   
377   // create a new vertex                          377   // create a new vertex
378   G4PrimaryVertex* vertex =                       378   G4PrimaryVertex* vertex = 
379     new G4PrimaryVertex(particle_position,part    379     new G4PrimaryVertex(particle_position,particle_time);
380                                                   380 
381   if(verbosityLevel >= 2)                         381   if(verbosityLevel >= 2)
382     G4cout << "Creating primaries and assignin    382     G4cout << "Creating primaries and assigning to vertex" << G4endl;
383   // create new primaries and set them to the     383   // create new primaries and set them to the vertex
384   G4double mass =  particle_definition->GetPDG    384   G4double mass =  particle_definition->GetPDGMass();
385   G4double energy = particle_energy + mass;       385   G4double energy = particle_energy + mass;
386   G4double pmom = std::sqrt(energy*energy-mass    386   G4double pmom = std::sqrt(energy*energy-mass*mass);
387   G4double px = pmom*particle_momentum_directi    387   G4double px = pmom*particle_momentum_direction.x();
388   G4double py = pmom*particle_momentum_directi    388   G4double py = pmom*particle_momentum_direction.y();
389   G4double pz = pmom*particle_momentum_directi    389   G4double pz = pmom*particle_momentum_direction.z();
390                                                   390   
391   if(verbosityLevel >= 1){                        391   if(verbosityLevel >= 1){
392     G4cout << "Particle name: "                   392     G4cout << "Particle name: " 
393            << particle_definition->GetParticle    393            << particle_definition->GetParticleName() << G4endl; 
394     G4cout << "       Energy: "<<particle_ener    394     G4cout << "       Energy: "<<particle_energy << G4endl;
395     G4cout << "     Position: "<<particle_posi    395     G4cout << "     Position: "<<particle_position<< G4endl; 
396     G4cout << "    Direction: "<<particle_mome    396     G4cout << "    Direction: "<<particle_momentum_direction << G4endl;
397     G4cout << " NumberOfParticlesToBeGenerated    397     G4cout << " NumberOfParticlesToBeGenerated: "
398            << NumberOfParticlesToBeGenerated <    398            << NumberOfParticlesToBeGenerated << G4endl;
399   }                                               399   }
400                                                   400 
401   for( G4int i=0; i<NumberOfParticlesToBeGener    401   for( G4int i=0; i<NumberOfParticlesToBeGenerated; ++i ) {
402     G4PrimaryParticle* particle =                 402     G4PrimaryParticle* particle =
403       new G4PrimaryParticle(particle_definitio    403       new G4PrimaryParticle(particle_definition,px,py,pz);
404     particle->SetMass( mass );                    404     particle->SetMass( mass );
405     particle->SetCharge( particle_charge );       405     particle->SetCharge( particle_charge );
406     particle->SetPolarization(particle_polariz    406     particle->SetPolarization(particle_polarization.x(),
407                               particle_polariz    407                               particle_polarization.y(),
408                               particle_polariz    408                               particle_polarization.z());
409     vertex->SetPrimary( particle );               409     vertex->SetPrimary( particle );
410   }                                               410   }
411   evt->AddPrimaryVertex( vertex );                411   evt->AddPrimaryVertex( vertex );
412   if(verbosityLevel > 1)                          412   if(verbosityLevel > 1)
413     G4cout << " Primary Vetex generated "<< G4    413     G4cout << " Primary Vetex generated "<< G4endl;   
414 }                                                 414 }
415                                                   415