| Geant4 Cross Reference |
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 class implementation file
29 // CERN Geneva Switzerland
30 //
31 //
32 // ------------ GammaRayTelPrimaryGeneratorAction ------
33 // by G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
34 //
35 // ************************************************************
36
37 #include "G4RunManager.hh"
38 #include "GammaRayTelPrimaryGeneratorAction.hh"
39
40 #include "GammaRayTelDetectorConstruction.hh"
41 #include "GammaRayTelPrimaryGeneratorMessenger.hh"
42
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4Event.hh"
46 #include "G4ParticleGun.hh"
47 #include "G4GeneralParticleSource.hh"
48 #include "G4ParticleTable.hh"
49 #include "G4ParticleDefinition.hh"
50 #include "Randomize.hh"
51
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction() {
55 detector = static_cast<const GammaRayTelDetectorConstruction*>(G4RunManager::GetRunManager()->GetUserDetectorConstruction());
56
57 // create a messenger for this class
58
59 gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
60
61 constexpr auto NUMBER_OF_PARTICLES{1};
62 particleGun = new G4ParticleGun(NUMBER_OF_PARTICLES);
63
64 // default particle kinematic
65
66 auto *particleTable = G4ParticleTable::GetParticleTable();
67 auto *particle = particleTable->FindParticle("e-");
68 particleGun->SetParticleDefinition(particle);
69 particleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., -1.));
70
71 constexpr auto PARTICLE_ENERGY{30. * MeV};
72 particleGun->SetParticleEnergy(PARTICLE_ENERGY);
73
74 auto position = 0.5 * (detector->GetWorldSizeZ());
75 particleGun->SetParticlePosition(G4ThreeVector(0. * cm, 0. * cm, position));
76 particleSource = new G4GeneralParticleSource();
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction() {
82 delete particleGun;
83 delete particleSource;
84 delete gunMessenger;
85 }
86
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88
89 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event *event) {
90 if (sourceGun) {
91 G4cout << "Using G4ParticleGun... " << G4endl;
92
93 // this function is called at the begining of event
94 //
95 G4double x0 = 0. * cm;
96 G4double y0 = 0. * cm;
97 G4double z0 = 0.5 * (detector->GetWorldSizeZ());
98
99 G4ThreeVector pos0;
100 auto vertex0 = G4ThreeVector(x0, y0, z0);
101 auto momentumDirection0 = G4ThreeVector(0., 0., -1.);
102
103 G4double theta;
104 G4double phi;
105 G4double y = 0.;
106 G4double f = 0.;
107 G4double theta0 = 0.;
108 G4double phi0 = 0.;
109
110 switch (sourceType) {
111 case 0:
112 particleGun->SetParticlePosition(vertex0);
113 particleGun->SetParticleMomentumDirection(momentumDirection0);
114 break;
115 case 1:
116 // GS: Generate random position on the 4PIsphere to create a uniform distribution
117 // GS: on the sphere
118 phi = G4UniformRand() * twopi;
119 do {
120 y = G4UniformRand() * 1.0;
121 theta = G4UniformRand() * pi;
122 f = std::sin(theta);
123 } while (y > f);
124 vertex0 = G4ThreeVector(1., 0., 0.);
125 vertex0.setMag(vertexRadius);
126 vertex0.setTheta(theta);
127 vertex0.setPhi(phi);
128 particleGun->SetParticlePosition(vertex0);
129
130 momentumDirection0 = G4ThreeVector(1., 0., 0.);
131
132 do {
133 phi = G4UniformRand() * twopi;
134 do {
135 y = G4UniformRand() * 1.0;
136 theta = G4UniformRand() * pi;
137 f = std::sin(theta);
138 } while (y > f);
139 momentumDirection0.setPhi(phi);
140 momentumDirection0.setTheta(theta);
141 } while (vertex0.dot(momentumDirection0) >= -0.7 * vertex0.mag());
142
143 particleGun->SetParticleMomentumDirection((G4ParticleMomentum) momentumDirection0);
144
145 break;
146 case 2:
147 // GS: Generate random position on the upper semi-sphere z > 0 to create a uniform distribution
148 // GS: on a plane
149 phi = G4UniformRand() * twopi;
150
151 do {
152 y = G4UniformRand() * 1.0;
153 theta = G4UniformRand() * halfpi;
154 f = std::sin(theta) * std::cos(theta);
155 } while (y > f);
156
157 vertex0 = G4ThreeVector(1., 0., 0.);
158
159 auto xy = detector->GetWorldSizeXY();
160 auto z = detector->GetWorldSizeZ();
161
162 if (vertexRadius > xy * 0.5) {
163 G4cout << "vertexRadius too big " << G4endl;
164 G4cout << "vertexRadius set to " << xy * 0.45 << G4endl;
165 vertexRadius = xy * 0.45;
166 }
167
168 if (vertexRadius > z * 0.5) {
169 G4cout << "vertexRadius too high " << G4endl;
170 G4cout << "vertexRadius set to " << z * 0.45 << G4endl;
171 vertexRadius = z * 0.45;
172 }
173
174 vertex0.setMag(vertexRadius);
175 vertex0.setTheta(theta);
176 vertex0.setPhi(phi);
177
178 // GS: Get the user defined direction for the primaries and
179 // GS: Rotate the random position according to the user defined direction for the particle
180
181 momentumDirection0 = particleGun->GetParticleMomentumDirection();
182 if (momentumDirection0.mag() > 0.001) {
183 theta0 = momentumDirection0.theta();
184 phi0 = momentumDirection0.phi();
185 }
186
187 if (theta0 != 0.) {
188 G4ThreeVector rotationAxis(1., 0., 0.);
189 rotationAxis.setPhi(phi0 + halfpi);
190 vertex0.rotate(theta0 + pi, rotationAxis);
191 }
192 particleGun->SetParticlePosition(vertex0);
193 break;
194 }
195
196 constexpr auto INITIAL_PARTICLE_ENERGY{100 * MeV};
197 G4double particleEnergy = INITIAL_PARTICLE_ENERGY;
198
199 switch (spectrumType) {
200 case 0: // Uniform energy (1 GeV - 10 GeV)
201 y = G4UniformRand();
202 particleEnergy = y * 9.0 * GeV + 1.0 * GeV;
203 G4cout << "Particle energy: " << particleEnergy / GeV << " LIN" << G4endl;
204 break;
205 case 1: // Logarithmic energy
206 y = G4UniformRand();
207 particleEnergy = std::pow(10, y) * GeV;
208 G4cout << "Particle energy: " << particleEnergy / GeV << " LOG" << G4endl;
209 break;
210 case 2: // Power law (-4)
211 do {
212 y = G4UniformRand() * 100000.0;
213 particleEnergy = G4UniformRand() * 10. * GeV;
214 f = std::pow(particleEnergy * (1 / GeV), -4.);
215 } while (y > f);
216 // particleGun->SetParticleEnergy(particleEnergy);
217 break;
218 case 3: // Monochromatic
219 particleEnergy = particleGun->GetParticleEnergy();
220 // 100 MeV;
221 G4cout << "Particle energy: " << particleEnergy << " MONOCHROMATIC" << G4endl;
222 break;
223 }
224 particleGun->SetParticleEnergy(particleEnergy);
225 G4cout << "Particle: " << particleGun->GetParticleDefinition()->GetParticleName() << G4endl;
226 particleGun->GeneratePrimaryVertex(event);
227 } else {
228 particleSource->GeneratePrimaryVertex(event);
229 }
230 }
231