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