ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DMXParticleSource.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DMXParticleSource.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // --------------------------------------------------------------
28 // GEANT 4 - Underground Dark Matter Detector Advanced Example
29 //
30 // For information related to this code contact: Alex Howard
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 // --------------------------------------------------------------
43 // This particle source is a shortened version of G4GeneralParticleSource by
44 // C Ferguson, F Lei & P Truscott (University of Southampton / DERA), with
45 // some minor modifications.
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"
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 
69 
71  particle_definition = NULL;
72  G4ThreeVector zero(0., 0., 0.);
74  particle_energy = 1.0*MeV;
76  particle_time = 0.0;
78  particle_charge = 0.0;
79 
80  SourcePosType = "Volume";
81  Shape = "NULL";
82  halfz = 0.;
83  Radius = 0.;
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 
102 }
103 
105 {
106  delete theMessenger;
107 }
108 
110 {
111  SourcePosType = PosType;
112 }
113 
115 {
116  Shape = shapeType;
117 }
118 
120 {
121  CentreCoords = coordsOfCentre;
122 }
123 
125 {
126  halfz = zhalf;
127 }
128 
130 {
131  Radius = radius;
132 }
133 
135 {
136  VolName = Vname;
137  if(verbosityLevel == 2) G4cout << VolName << G4endl;
138 
139  // checks if selected volume exists
140  G4VPhysicalVolume *tempPV = NULL;
141  G4PhysicalVolumeStore *PVStore = 0;
142  G4String theRequiredVolumeName = VolName;
144  G4int i = 0;
145  G4bool found = false;
146  if(verbosityLevel == 2) G4cout << PVStore->size() << G4endl;
147 
148  // recasting required since PVStore->size() is actually a signed int...
149  while (!found && i<(G4int)PVStore->size())
150  {
151  tempPV = (*PVStore)[i];
152  found = tempPV->GetName() == theRequiredVolumeName;
153  if(verbosityLevel == 2)
154  G4cout << i << " " << " " << tempPV->GetName()
155  << " " << theRequiredVolumeName << " " << found << G4endl;
156  if (!found)
157  {i++;}
158  }
159 
160  // found = true then the volume exists else it doesnt.
161  if(found == true) {
162  if(verbosityLevel >= 1)
163  G4cout << "Volume " << VolName << " exists" << G4endl;
164  Confine = true;
165  }
166  else if(VolName=="NULL")
167  Confine = false;
168  else {
169  G4cout << " **** Error: Volume does not exist **** " << G4endl;
170  G4cout << " Ignoring confine condition" << G4endl;
171  VolName = "NULL";
172  Confine = false;
173  }
174 
175 }
176 
177 
179 {
180  AngDistType = atype;
181 }
182 
183 
185 {
186  // Generates Points given the point source.
187  if(SourcePosType == "Point")
189  else
190  if(verbosityLevel >= 1)
191  G4cout << "Error SourcePosType is not set to Point" << G4endl;
192 }
193 
194 
196 {
197  G4ThreeVector RandPos;
198  G4double x=0., y=0., z=0.;
199 
200  if(SourcePosType != "Volume" && verbosityLevel >= 1)
201  G4cout << "Error SourcePosType not Volume" << G4endl;
202 
203  if(Shape == "Sphere") {
204  x = Radius*2.;
205  y = Radius*2.;
206  z = Radius*2.;
207  while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) {
208  x = G4UniformRand();
209  y = G4UniformRand();
210  z = G4UniformRand();
211 
212  x = (x*2.*Radius) - Radius;
213  y = (y*2.*Radius) - Radius;
214  z = (z*2.*Radius) - Radius;
215  }
216  }
217 
218  else if(Shape == "Cylinder") {
219  x = Radius*2.;
220  y = Radius*2.;
221  while(((x*x)+(y*y)) > (Radius*Radius)) {
222  x = G4UniformRand();
223  y = G4UniformRand();
224  z = G4UniformRand();
225  x = (x*2.*Radius) - Radius;
226  y = (y*2.*Radius) - Radius;
227  z = (z*2.*halfz) - halfz;
228  }
229  }
230 
231  else
232  G4cout << "Error: Volume Shape Does Not Exist" << G4endl;
233 
234  RandPos.setX(x);
235  RandPos.setY(y);
236  RandPos.setZ(z);
237  particle_position = CentreCoords + RandPos;
238 
239 }
240 
241 
243 {
244 
245  // Method to check point is within the volume specified
246  if(Confine == false)
247  G4cout << "Error: Confine is false" << G4endl;
248  G4ThreeVector null(0.,0.,0.);
249  G4ThreeVector *ptr;
250  ptr = &null;
251 
252  // Check particle_position is within VolName
253  G4VPhysicalVolume *theVolume;
255  G4String theVolName = theVolume->GetName();
256  if(theVolName == VolName) {
257  if(verbosityLevel >= 1)
258  G4cout << "Particle is in volume " << VolName << G4endl;
259  return(true);
260  }
261  else
262  return(false);
263 }
264 
265 
267  (G4ParticleMomentum aDirection) {
268 
269  particle_momentum_direction = aDirection.unit();
270 }
271 
272 
274 {
275 
276  G4double rndm, rndm2;
277  G4double px, py, pz;
278 
279  G4double sintheta, sinphi, costheta, cosphi;
280  rndm = G4UniformRand();
281  costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
282  sintheta = std::sqrt(1. - costheta*costheta);
283 
284  rndm2 = G4UniformRand();
285  Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
286  sinphi = std::sin(Phi);
287  cosphi = std::cos(Phi);
288 
289  px = -sintheta * cosphi;
290  py = -sintheta * sinphi;
291  pz = -costheta;
292 
293  G4double ResMag = std::sqrt((px*px) + (py*py) + (pz*pz));
294  px = px/ResMag;
295  py = py/ResMag;
296  pz = pz/ResMag;
297 
301 
302  // particle_momentum_direction now holds unit momentum vector.
303  if(verbosityLevel >= 2)
304  G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
305 }
306 
307 
309 {
310  EnergyDisType = DisType;
311 }
312 
314 {
315  MonoEnergy = menergy;
316 }
317 
319 {
321 }
322 
323 
325 {
326  verbosityLevel = vL;
327  G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
328 }
329 
331  (G4ParticleDefinition* aParticleDefinition)
332 {
333  particle_definition = aParticleDefinition;
334  particle_charge = particle_definition->GetPDGCharge();
335 }
336 
337 
339 {
340 
341  if(particle_definition==NULL) {
342  G4cout << "No particle has been defined!" << G4endl;
343  return;
344  }
345 
346  // Position
347  G4bool srcconf = false;
348  G4int LoopCount = 0;
349 
350  while(srcconf == false) {
351  if(SourcePosType == "Point")
353  else if(SourcePosType == "Volume")
355  else {
356  G4cout << "Error: SourcePosType undefined" << G4endl;
357  G4cout << "Generating point source" << G4endl;
359  }
360  if(Confine == true) {
361  srcconf = IsSourceConfined();
362  // if source in confined srcconf = true terminating the loop
363  // if source isnt confined srcconf = false and loop continues
364  }
365  else if(Confine == false)
366  srcconf = true; // terminate loop
367 
368  LoopCount++;
369  if(LoopCount == 100000) {
370  G4cout << "*************************************" << G4endl;
371  G4cout << "LoopCount = 100000" << G4endl;
372  G4cout << "Either the source distribution >> confinement" << G4endl;
373  G4cout << "or any confining volume may not overlap with" << G4endl;
374  G4cout << "the source distribution or any confining volumes" << G4endl;
375  G4cout << "may not exist"<< G4endl;
376  G4cout << "If you have set confine then this will be ignored" <<G4endl;
377  G4cout << "for this event." << G4endl;
378  G4cout << "*************************************" << G4endl;
379  srcconf = true; //Avoids an infinite loop
380  }
381  }
382 
383  // Angular stuff
384  if(AngDistType == "iso")
386  else if(AngDistType == "direction")
388  else
389  G4cout << "Error: AngDistType has unusual value" << G4endl;
390  // Energy stuff
391  if(EnergyDisType == "Mono")
393  else
394  G4cout << "Error: EnergyDisType has unusual value" << G4endl;
395 
396  // create a new vertex
397  G4PrimaryVertex* vertex =
399 
400  if(verbosityLevel >= 2)
401  G4cout << "Creating primaries and assigning to vertex" << G4endl;
402  // create new primaries and set them to the vertex
405  G4double pmom = std::sqrt(energy*energy-mass*mass);
409 
410  if(verbosityLevel >= 1){
411  G4cout << "Particle name: "
413  G4cout << " Energy: "<<particle_energy << G4endl;
414  G4cout << " Position: "<<particle_position<< G4endl;
415  G4cout << " Direction: "<<particle_momentum_direction << G4endl;
416  G4cout << " NumberOfParticlesToBeGenerated: "
418  }
419 
420 
421  for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
424  particle->SetMass( mass );
425  particle->SetCharge( particle_charge );
429  vertex->SetPrimary( particle );
430  }
431  evt->AddPrimaryVertex( vertex );
432  if(verbosityLevel > 1)
433  G4cout << " Primary Vetex generated "<< G4endl;
434 }
435 
436 
437 
438