Geant4_10
DMXParticleSource.cc
Go to the documentation of this file.
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 
70  NumberOfParticlesToBeGenerated = 1;
71  particle_definition = NULL;
72  G4ThreeVector zero(0., 0., 0.);
73  particle_momentum_direction = G4ParticleMomentum(1., 0., 0.);
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 DMXParticleSourceMessenger(this);
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")
188  particle_position = CentreCoords;
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;
254  theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
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 
298  particle_momentum_direction.setX(px);
299  particle_momentum_direction.setY(py);
300  particle_momentum_direction.setZ(pz);
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 {
320  particle_energy = MonoEnergy;
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")
387  SetParticleMomentumDirection(particle_momentum_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 =
398  new G4PrimaryVertex(particle_position,particle_time);
399 
400  if(verbosityLevel >= 2)
401  G4cout << "Creating primaries and assigning to vertex" << G4endl;
402  // create new primaries and set them to the vertex
403  G4double mass = particle_definition->GetPDGMass();
404  G4double energy = particle_energy + mass;
405  G4double pmom = std::sqrt(energy*energy-mass*mass);
406  G4double px = pmom*particle_momentum_direction.x();
407  G4double py = pmom*particle_momentum_direction.y();
408  G4double pz = pmom*particle_momentum_direction.z();
409 
410  if(verbosityLevel >= 1){
411  G4cout << "Particle name: "
412  << particle_definition->GetParticleName() << G4endl;
413  G4cout << " Energy: "<<particle_energy << G4endl;
414  G4cout << " Position: "<<particle_position<< G4endl;
415  G4cout << " Direction: "<<particle_momentum_direction << G4endl;
416  G4cout << " NumberOfParticlesToBeGenerated: "
417  << NumberOfParticlesToBeGenerated << G4endl;
418  }
419 
420 
421  for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
422  G4PrimaryParticle* particle =
423  new G4PrimaryParticle(particle_definition,px,py,pz);
424  particle->SetMass( mass );
425  particle->SetCharge( particle_charge );
426  particle->SetPolarization(particle_polarization.x(),
427  particle_polarization.y(),
428  particle_polarization.z());
429  vertex->SetPrimary( particle );
430  }
431  evt->AddPrimaryVertex( vertex );
432  if(verbosityLevel > 1)
433  G4cout << " Primary Vetex generated "<< G4endl;
434 }
435 
436 
437 
438 
void SetCentreCoords(G4ThreeVector)
void SetAngDistType(G4String)
double x() const
DOUBLE PRECISION function rndm(RDUMMY)
Definition: dpm25nulib.f:1460
void SetHalfZ(G4double)
void AddPrimaryVertex(G4PrimaryVertex *aPrimaryVertex)
Definition: G4Event.hh:143
G4Navigator * GetNavigatorForTracking() const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
void setY(double)
const G4String & GetParticleName() const
double z() const
void setZ(double)
void setX(double)
void SetParticleMomentumDirection(G4ParticleMomentum)
static G4PhysicalVolumeStore * GetInstance()
Double_t y
Definition: plot.C:279
void GeneratePrimaryVertex(G4Event *evt)
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
void SetPosDisType(G4String)
double energy
Definition: plottest35.C:25
const G4String & GetName() const
bool G4bool
Definition: G4Types.hh:79
void SetMass(G4double mas)
static G4TransportationManager * GetTransportationManager()
void SetRadius(G4double)
G4double GetPDGMass() const
void SetPosDisShape(G4String)
Hep3Vector unit() const
void SetCharge(G4double chg)
double y() const
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
Definition: G4Navigator.cc:118
tuple z
Definition: test.py:28
void SetPrimary(G4PrimaryParticle *pp)
#define G4endl
Definition: G4ios.hh:61
void SetEnergyDisType(G4String)
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
void SetParticleDefinition(G4ParticleDefinition *aParticleDefinition)
void SetPolarization(const G4ThreeVector &pol)
void ConfineSourceToVolume(G4String)
void SetMonoEnergy(G4double)