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 ]

  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 // --------------------------------------------------------------
 42 //////////////////////////////////////////////////////////////////////////////
 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.
 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 = 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);
100   gNavigator = G4TransportationManager::GetTransportationManager()
101     ->GetNavigatorForTracking();
102 }
103 
104 DMXParticleSource::~DMXParticleSource()
105 {
106   delete theMessenger;
107 }
108 
109 void DMXParticleSource::SetPosDisType(G4String PosType) 
110 {
111   SourcePosType = PosType;
112 }
113 
114 void DMXParticleSource::SetPosDisShape(G4String shapeType)
115 {
116   Shape = shapeType;
117 }
118 
119 void DMXParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre)
120 {
121   CentreCoords = coordsOfCentre;
122 }
123 
124 void DMXParticleSource::SetHalfZ(G4double zhalf)
125 {
126   halfz = zhalf;
127 }
128 
129 void DMXParticleSource::SetRadius(G4double radius)
130 {
131   Radius = radius;
132 }
133 
134 void DMXParticleSource::ConfineSourceToVolume(G4String Vname)
135 {
136   VolName = Vname;
137   if(verbosityLevel == 2) G4cout << VolName << G4endl;
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->size() << G4endl;
146 
147   tempPV = PVStore->GetVolume(theRequiredVolumeName, false);
148   if (tempPV != nullptr) { found = true; }
149 
150   // found = true then the volume exists else it doesnt.
151   if(found == true) {
152     if(verbosityLevel >= 1)
153       G4cout << "Volume " << VolName << " exists" << G4endl;
154     Confine = true;
155   }
156   else if(VolName=="NULL")
157     Confine = false;
158   else {
159     G4cout << " **** Error: Volume does not exist **** " << G4endl;
160     G4cout << " Ignoring confine condition" << G4endl;
161     VolName = "NULL";
162     Confine = false;
163   }
164 }
165 
166 void DMXParticleSource::SetAngDistType(G4String atype)
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 set to Point" << G4endl;
179 }
180 
181 void DMXParticleSource::GeneratePointsInVolume()
182 {
183   G4ThreeVector RandPos;
184   G4double x=0., y=0., z=0.;
185   
186   if(SourcePosType != "Volume" && verbosityLevel >= 1)
187     G4cout << "Error SourcePosType not Volume" << G4endl;
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 Exist" << G4endl;
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 volume specified
231   if(Confine == false)
232     G4cout << "Error: Confine is false" << G4endl;
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->LocateGlobalPointAndSetup(particle_position,ptr,true);
239   G4String theVolName = theVolume->GetName();
240   if(theVolName == VolName) {
241     if(verbosityLevel >= 1)
242       G4cout << "Particle is in volume " << VolName << G4endl;
243     return(true);
244   }
245   else
246     return(false);
247 }
248 
249 void DMXParticleSource::SetParticleMomentumDirection
250    (G4ParticleMomentum aDirection) {
251 
252   particle_momentum_direction =  aDirection.unit();
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::cos(MinTheta)
264                                 - std::cos(MaxTheta));
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) + (pz*pz));
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 unit momentum vector.
286   if(verbosityLevel >= 2)
287     G4cout << "Generating isotropic vector: "
288            << particle_momentum_direction << G4endl;
289 }
290 
291 void DMXParticleSource::SetEnergyDisType(G4String DisType)
292 {
293   EnergyDisType = DisType;
294 }
295 
296 void DMXParticleSource::SetMonoEnergy(G4double menergy)
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: " << verbosityLevel << G4endl;
310 }
311 
312 void DMXParticleSource::SetParticleDefinition
313   (G4ParticleDefinition* aParticleDefinition)
314 {
315   particle_definition = aParticleDefinition;
316   particle_charge = particle_definition->GetPDGCharge();
317 }
318 
319 void DMXParticleSource::GeneratePrimaryVertex(G4Event *evt)
320 {
321 
322   if(particle_definition==nullptr) {
323     G4cout << "No particle has been defined!" << G4endl;
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 undefined" << G4endl;
338       G4cout << "Generating point source" << G4endl;
339       GeneratePointSource();
340     }
341     if(Confine == true) {
342       srcconf = IsSourceConfined();
343       // if source in confined srcconf = true terminating the loop
344       // if source isnt confined srcconf = false and loop continues
345     }
346     else if(Confine == false)
347       srcconf = true; // terminate loop
348     
349     ++LoopCount;
350     if(LoopCount == 100000) {
351       G4cout << "*************************************" << G4endl;
352         G4cout << "LoopCount = 100000" << G4endl;
353         G4cout << "Either the source distribution >> confinement" << G4endl;
354         G4cout << "or any confining volume may not overlap with" << G4endl;
355         G4cout << "the source distribution or any confining volumes" << G4endl;
356         G4cout << "may not exist"<< G4endl;
357         G4cout << "If you have set confine then this will be ignored" <<G4endl;
358         G4cout << "for this event." << G4endl;
359         G4cout << "*************************************" << G4endl;
360         srcconf = true; //Avoids an infinite loop
361       }
362   }
363 
364   // Angular stuff
365   if(AngDistType == "iso")
366     GenerateIsotropicFlux();
367   else if(AngDistType == "direction")
368     SetParticleMomentumDirection(particle_momentum_direction);
369   else
370     G4cout << "Error: AngDistType has unusual value" << G4endl;
371   // Energy stuff
372   if(EnergyDisType == "Mono")
373     GenerateMonoEnergetic();
374   else
375     G4cout << "Error: EnergyDisType has unusual value" << G4endl;
376   
377   // create a new vertex
378   G4PrimaryVertex* vertex = 
379     new G4PrimaryVertex(particle_position,particle_time);
380 
381   if(verbosityLevel >= 2)
382     G4cout << "Creating primaries and assigning to vertex" << G4endl;
383   // create new primaries and set them to the vertex
384   G4double mass =  particle_definition->GetPDGMass();
385   G4double energy = particle_energy + mass;
386   G4double pmom = std::sqrt(energy*energy-mass*mass);
387   G4double px = pmom*particle_momentum_direction.x();
388   G4double py = pmom*particle_momentum_direction.y();
389   G4double pz = pmom*particle_momentum_direction.z();
390   
391   if(verbosityLevel >= 1){
392     G4cout << "Particle name: " 
393            << particle_definition->GetParticleName() << G4endl; 
394     G4cout << "       Energy: "<<particle_energy << G4endl;
395     G4cout << "     Position: "<<particle_position<< G4endl; 
396     G4cout << "    Direction: "<<particle_momentum_direction << G4endl;
397     G4cout << " NumberOfParticlesToBeGenerated: "
398            << NumberOfParticlesToBeGenerated << G4endl;
399   }
400 
401   for( G4int i=0; i<NumberOfParticlesToBeGenerated; ++i ) {
402     G4PrimaryParticle* particle =
403       new G4PrimaryParticle(particle_definition,px,py,pz);
404     particle->SetMass( mass );
405     particle->SetCharge( particle_charge );
406     particle->SetPolarization(particle_polarization.x(),
407                               particle_polarization.y(),
408                               particle_polarization.z());
409     vertex->SetPrimary( particle );
410   }
411   evt->AddPrimaryVertex( vertex );
412   if(verbosityLevel > 1)
413     G4cout << " Primary Vetex generated "<< G4endl;   
414 }
415