source: trunk/examples/advanced/underground_physics/src/DMXParticleSource.cc@ 1022

Last change on this file since 1022 was 807, checked in by garnier, 17 years ago

update
File size: 12.4 KB
RevLine 
[807]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 "DMXParticleSource.hh"
49
50#include "G4PrimaryParticle.hh"
51#include "G4Event.hh"
52#include "Randomize.hh"
53#include <cmath>
54#include "G4TransportationManager.hh"
55#include "G4VPhysicalVolume.hh"
56#include "G4PhysicalVolumeStore.hh"
57#include "G4ParticleTable.hh"
58#include "G4ParticleDefinition.hh"
59#include "G4IonTable.hh"
60#include "G4Ions.hh"
61#include "G4TrackingManager.hh"
62#include "G4Track.hh"
63
64
65DMXParticleSource::DMXParticleSource() {
66
67 NumberOfParticlesToBeGenerated = 1;
68 particle_definition = NULL;
69 G4ThreeVector zero(0., 0., 0.);
70 particle_momentum_direction = G4ParticleMomentum(1., 0., 0.);
71 particle_energy = 1.0*MeV;
72 particle_position = zero;
73 particle_time = 0.0;
74 particle_polarization = zero;
75 particle_charge = 0.0;
76
77 SourcePosType = "Volume";
78 Shape = "NULL";
79 halfz = 0.;
80 Radius = 0.;
81 CentreCoords = zero;
82 Confine = false;
83 VolName = "NULL";
84
85 AngDistType = "iso";
86 MinTheta = 0.;
87 MaxTheta = pi;
88 MinPhi = 0.;
89 MaxPhi = twopi;
90
91 EnergyDisType = "Mono";
92 MonoEnergy = 1*MeV;
93
94 verbosityLevel = 0;
95
96 theMessenger = new DMXParticleSourceMessenger(this);
97 gNavigator = G4TransportationManager::GetTransportationManager()
98 ->GetNavigatorForTracking();
99}
100
101DMXParticleSource::~DMXParticleSource()
102{
103 delete theMessenger;
104}
105
106void DMXParticleSource::SetPosDisType(G4String PosType)
107{
108 SourcePosType = PosType;
109}
110
111void DMXParticleSource::SetPosDisShape(G4String shapeType)
112{
113 Shape = shapeType;
114}
115
116void DMXParticleSource::SetCentreCoords(G4ThreeVector coordsOfCentre)
117{
118 CentreCoords = coordsOfCentre;
119}
120
121void DMXParticleSource::SetHalfZ(G4double zhalf)
122{
123 halfz = zhalf;
124}
125
126void DMXParticleSource::SetRadius(G4double rad)
127{
128 Radius = rad;
129}
130
131void DMXParticleSource::ConfineSourceToVolume(G4String Vname)
132{
133 VolName = Vname;
134 if(verbosityLevel == 2) G4cout << VolName << G4endl;
135
136 // checks if selected volume exists
137 G4VPhysicalVolume *tempPV = NULL;
138 G4PhysicalVolumeStore *PVStore = 0;
139 G4String theRequiredVolumeName = VolName;
140 PVStore = G4PhysicalVolumeStore::GetInstance();
141 G4int i = 0;
142 G4bool found = false;
143 if(verbosityLevel == 2) G4cout << PVStore->size() << G4endl;
144
145 // recasting required since PVStore->size() is actually a signed int...
146 while (!found && i<(G4int)PVStore->size())
147 {
148 tempPV = (*PVStore)[i];
149 found = tempPV->GetName() == theRequiredVolumeName;
150 if(verbosityLevel == 2)
151 G4cout << i << " " << " " << tempPV->GetName()
152 << " " << theRequiredVolumeName << " " << found << G4endl;
153 if (!found)
154 {i++;}
155 }
156
157 // found = true then the volume exists else it doesnt.
158 if(found == true) {
159 if(verbosityLevel >= 1)
160 G4cout << "Volume " << VolName << " exists" << G4endl;
161 Confine = true;
162 }
163 else if(VolName=="NULL")
164 Confine = false;
165 else {
166 G4cout << " **** Error: Volume does not exist **** " << G4endl;
167 G4cout << " Ignoring confine condition" << G4endl;
168 VolName = "NULL";
169 Confine = false;
170 }
171
172}
173
174
175void DMXParticleSource::SetAngDistType(G4String atype)
176{
177 AngDistType = atype;
178}
179
180
181void DMXParticleSource::GeneratePointSource()
182{
183 // Generates Points given the point source.
184 if(SourcePosType == "Point")
185 particle_position = CentreCoords;
186 else
187 if(verbosityLevel >= 1)
188 G4cout << "Error SourcePosType is not set to Point" << G4endl;
189}
190
191
192void DMXParticleSource::GeneratePointsInVolume()
193{
194 G4ThreeVector RandPos;
195 G4double x=0., y=0., z=0.;
196
197 if(SourcePosType != "Volume" && verbosityLevel >= 1)
198 G4cout << "Error SourcePosType not Volume" << G4endl;
199
200 if(Shape == "Sphere") {
201 x = Radius*2.;
202 y = Radius*2.;
203 z = Radius*2.;
204 while(((x*x)+(y*y)+(z*z)) > (Radius*Radius)) {
205 x = G4UniformRand();
206 y = G4UniformRand();
207 z = G4UniformRand();
208
209 x = (x*2.*Radius) - Radius;
210 y = (y*2.*Radius) - Radius;
211 z = (z*2.*Radius) - Radius;
212 }
213 }
214
215 else if(Shape == "Cylinder") {
216 x = Radius*2.;
217 y = Radius*2.;
218 while(((x*x)+(y*y)) > (Radius*Radius)) {
219 x = G4UniformRand();
220 y = G4UniformRand();
221 z = G4UniformRand();
222 x = (x*2.*Radius) - Radius;
223 y = (y*2.*Radius) - Radius;
224 z = (z*2.*halfz) - halfz;
225 }
226 }
227
228 else
229 G4cout << "Error: Volume Shape Does Not Exist" << G4endl;
230
231 RandPos.setX(x);
232 RandPos.setY(y);
233 RandPos.setZ(z);
234 particle_position = CentreCoords + RandPos;
235
236}
237
238
239G4bool DMXParticleSource::IsSourceConfined()
240{
241
242 // Method to check point is within the volume specified
243 if(Confine == false)
244 G4cout << "Error: Confine is false" << G4endl;
245 G4ThreeVector null(0.,0.,0.);
246 G4ThreeVector *ptr;
247 ptr = &null;
248
249 // Check particle_position is within VolName
250 G4VPhysicalVolume *theVolume;
251 theVolume=gNavigator->LocateGlobalPointAndSetup(particle_position,ptr,true);
252 G4String theVolName = theVolume->GetName();
253 if(theVolName == VolName) {
254 if(verbosityLevel >= 1)
255 G4cout << "Particle is in volume " << VolName << G4endl;
256 return(true);
257 }
258 else
259 return(false);
260}
261
262
263void DMXParticleSource::SetParticleMomentumDirection
264 (G4ParticleMomentum aDirection) {
265
266 particle_momentum_direction = aDirection.unit();
267}
268
269
270void DMXParticleSource::GenerateIsotropicFlux()
271{
272
273 G4double rndm, rndm2;
274 G4double px, py, pz;
275
276 G4double sintheta, sinphi, costheta, cosphi;
277 rndm = G4UniformRand();
278 costheta = std::cos(MinTheta) - rndm * (std::cos(MinTheta) - std::cos(MaxTheta));
279 sintheta = std::sqrt(1. - costheta*costheta);
280
281 rndm2 = G4UniformRand();
282 Phi = MinPhi + (MaxPhi - MinPhi) * rndm2;
283 sinphi = std::sin(Phi);
284 cosphi = std::cos(Phi);
285
286 px = -sintheta * cosphi;
287 py = -sintheta * sinphi;
288 pz = -costheta;
289
290 G4double ResMag = std::sqrt((px*px) + (py*py) + (pz*pz));
291 px = px/ResMag;
292 py = py/ResMag;
293 pz = pz/ResMag;
294
295 particle_momentum_direction.setX(px);
296 particle_momentum_direction.setY(py);
297 particle_momentum_direction.setZ(pz);
298
299 // particle_momentum_direction now holds unit momentum vector.
300 if(verbosityLevel >= 2)
301 G4cout << "Generating isotropic vector: " << particle_momentum_direction << G4endl;
302}
303
304
305void DMXParticleSource::SetEnergyDisType(G4String DisType)
306{
307 EnergyDisType = DisType;
308}
309
310void DMXParticleSource::SetMonoEnergy(G4double menergy)
311{
312 MonoEnergy = menergy;
313}
314
315void DMXParticleSource::GenerateMonoEnergetic()
316{
317 particle_energy = MonoEnergy;
318}
319
320
321void DMXParticleSource::SetVerbosity(int vL)
322{
323 verbosityLevel = vL;
324 G4cout << "Verbosity Set to: " << verbosityLevel << G4endl;
325}
326
327void DMXParticleSource::SetParticleDefinition
328 (G4ParticleDefinition* aParticleDefinition)
329{
330 particle_definition = aParticleDefinition;
331 particle_charge = particle_definition->GetPDGCharge();
332}
333
334
335void DMXParticleSource::GeneratePrimaryVertex(G4Event *evt)
336{
337
338 if(particle_definition==NULL) {
339 G4cout << "No particle has been defined!" << G4endl;
340 return;
341 }
342
343 // Position
344 G4bool srcconf = false;
345 G4int LoopCount = 0;
346
347 while(srcconf == false) {
348 if(SourcePosType == "Point")
349 GeneratePointSource();
350 else if(SourcePosType == "Volume")
351 GeneratePointsInVolume();
352 else {
353 G4cout << "Error: SourcePosType undefined" << G4endl;
354 G4cout << "Generating point source" << G4endl;
355 GeneratePointSource();
356 }
357 if(Confine == true) {
358 srcconf = IsSourceConfined();
359 // if source in confined srcconf = true terminating the loop
360 // if source isnt confined srcconf = false and loop continues
361 }
362 else if(Confine == false)
363 srcconf = true; // terminate loop
364
365 LoopCount++;
366 if(LoopCount == 100000) {
367 G4cout << "*************************************" << G4endl;
368 G4cout << "LoopCount = 100000" << G4endl;
369 G4cout << "Either the source distribution >> confinement" << G4endl;
370 G4cout << "or any confining volume may not overlap with" << G4endl;
371 G4cout << "the source distribution or any confining volumes" << G4endl;
372 G4cout << "may not exist"<< G4endl;
373 G4cout << "If you have set confine then this will be ignored" <<G4endl;
374 G4cout << "for this event." << G4endl;
375 G4cout << "*************************************" << G4endl;
376 srcconf = true; //Avoids an infinite loop
377 }
378 }
379
380 // Angular stuff
381 if(AngDistType == "iso")
382 GenerateIsotropicFlux();
383 else if(AngDistType == "direction")
384 SetParticleMomentumDirection(particle_momentum_direction);
385 else
386 G4cout << "Error: AngDistType has unusual value" << G4endl;
387 // Energy stuff
388 if(EnergyDisType == "Mono")
389 GenerateMonoEnergetic();
390 else
391 G4cout << "Error: EnergyDisType has unusual value" << G4endl;
392
393 // create a new vertex
394 G4PrimaryVertex* vertex =
395 new G4PrimaryVertex(particle_position,particle_time);
396
397 if(verbosityLevel >= 2)
398 G4cout << "Creating primaries and assigning to vertex" << G4endl;
399 // create new primaries and set them to the vertex
400 G4double mass = particle_definition->GetPDGMass();
401 G4double energy = particle_energy + mass;
402 G4double pmom = std::sqrt(energy*energy-mass*mass);
403 G4double px = pmom*particle_momentum_direction.x();
404 G4double py = pmom*particle_momentum_direction.y();
405 G4double pz = pmom*particle_momentum_direction.z();
406
407 if(verbosityLevel >= 1){
408 G4cout << "Particle name: "
409 << particle_definition->GetParticleName() << G4endl;
410 G4cout << " Energy: "<<particle_energy << G4endl;
411 G4cout << " Position: "<<particle_position<< G4endl;
412 G4cout << " Direction: "<<particle_momentum_direction << G4endl;
413 G4cout << " NumberOfParticlesToBeGenerated: "
414 << NumberOfParticlesToBeGenerated << G4endl;
415 }
416
417
418 for( G4int i=0; i<NumberOfParticlesToBeGenerated; i++ ) {
419 G4PrimaryParticle* particle =
420 new G4PrimaryParticle(particle_definition,px,py,pz);
421 particle->SetMass( mass );
422 particle->SetCharge( particle_charge );
423 particle->SetPolarization(particle_polarization.x(),
424 particle_polarization.y(),
425 particle_polarization.z());
426 vertex->SetPrimary( particle );
427 }
428 evt->AddPrimaryVertex( vertex );
429 if(verbosityLevel > 1)
430 G4cout << " Primary Vetex generated "<< G4endl;
431}
432
433
434
435
Note: See TracBrowser for help on using the repository browser.