| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: GammaRayTelPrimaryGeneratorAction.cc,v 1.10 2006/06/29 15:56:58 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01 $
27 // ------------------------------------------- 29 // ------------------------------------------------------------
28 // GEANT 4 class implementation file 30 // GEANT 4 class implementation file
29 // CERN Geneva Switzerland 31 // CERN Geneva Switzerland
30 // 32 //
31 // 33 //
32 // ------------ GammaRayTelPrimaryGenerat 34 // ------------ GammaRayTelPrimaryGeneratorAction ------
33 // by G.Santin, F.Longo & R.Giannit 35 // by G.Santin, F.Longo & R.Giannitrapani (13 nov 2000)
34 // 36 //
35 // ******************************************* 37 // ************************************************************
36 38
>> 39 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 41
37 #include "G4RunManager.hh" 42 #include "G4RunManager.hh"
38 #include "GammaRayTelPrimaryGeneratorAction.hh 43 #include "GammaRayTelPrimaryGeneratorAction.hh"
39 44
40 #include "GammaRayTelDetectorConstruction.hh" 45 #include "GammaRayTelDetectorConstruction.hh"
41 #include "GammaRayTelPrimaryGeneratorMessenger 46 #include "GammaRayTelPrimaryGeneratorMessenger.hh"
42 47
43 #include "G4PhysicalConstants.hh" <<
44 #include "G4SystemOfUnits.hh" <<
45 #include "G4Event.hh" 48 #include "G4Event.hh"
46 #include "G4ParticleGun.hh" 49 #include "G4ParticleGun.hh"
47 #include "G4GeneralParticleSource.hh" 50 #include "G4GeneralParticleSource.hh"
48 #include "G4ParticleTable.hh" 51 #include "G4ParticleTable.hh"
49 #include "G4ParticleDefinition.hh" 52 #include "G4ParticleDefinition.hh"
50 #include "Randomize.hh" 53 #include "Randomize.hh"
51 54
52 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 56
54 GammaRayTelPrimaryGeneratorAction::GammaRayTel << 57 GammaRayTelPrimaryGeneratorAction::GammaRayTelPrimaryGeneratorAction()
55 detector = static_cast<const GammaRayTelDe << 58 :rndmFlag("off"),nSourceType(0),nSpectrumType(0)
56 << 59 {
57 // create a messenger for this class << 60 G4RunManager* runManager = G4RunManager::GetRunManager();
58 << 61 GammaRayTelDetector =
59 gunMessenger = new GammaRayTelPrimaryGener << 62 (GammaRayTelDetectorConstruction*)(runManager->GetUserDetectorConstruction());
60 << 63
61 constexpr auto NUMBER_OF_PARTICLES{1}; << 64
62 particleGun = new G4ParticleGun(NUMBER_OF_ << 65 //create a messenger for this class
63 << 66
64 // default particle kinematic << 67 gunMessenger = new GammaRayTelPrimaryGeneratorMessenger(this);
65 << 68
66 auto *particleTable = G4ParticleTable::Get << 69 G4int n_particle = 1;
67 auto *particle = particleTable->FindPartic << 70
68 particleGun->SetParticleDefinition(particl << 71 if (sourceGun)
69 particleGun->SetParticleMomentumDirection( << 72 {
70 << 73 particleGun = new G4ParticleGun(n_particle);
71 constexpr auto PARTICLE_ENERGY{30. * MeV}; << 74 // default particle kinematic
72 particleGun->SetParticleEnergy(PARTICLE_EN << 75
73 << 76 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
74 auto position = 0.5 * (detector->GetWorldS << 77 G4String particleName;
75 particleGun->SetParticlePosition(G4ThreeVe << 78 G4ParticleDefinition* particle
76 particleSource = new G4GeneralParticleSour << 79 = particleTable->FindParticle(particleName="e-");
>> 80 particleGun->SetParticleDefinition(particle);
>> 81 particleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,-1.));
>> 82 particleGun->SetParticleEnergy(30.*MeV);
>> 83 G4double position = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
>> 84 particleGun->SetParticlePosition(G4ThreeVector(0.*cm,0.*cm,position));
>> 85 }
>> 86 else
>> 87 {
>> 88 particleSource = new G4GeneralParticleSource();
>> 89 }
>> 90
77 } 91 }
78 92
79 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 94
81 GammaRayTelPrimaryGeneratorAction::~GammaRayTe << 95 GammaRayTelPrimaryGeneratorAction::~GammaRayTelPrimaryGeneratorAction()
>> 96 {
>> 97 if (sourceGun)
82 delete particleGun; 98 delete particleGun;
>> 99 else
83 delete particleSource; 100 delete particleSource;
84 delete gunMessenger; << 101
>> 102 delete gunMessenger;
85 } 103 }
86 104
87 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 106
89 void GammaRayTelPrimaryGeneratorAction::Genera << 107 void GammaRayTelPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
90 if (sourceGun) { << 108 {
91 G4cout << "Using G4ParticleGun... " << << 109 if (sourceGun)
92 << 110 {
93 // this function is called at the begi << 111
94 // << 112 //this function is called at the begining of event
95 G4double x0 = 0. * cm; << 113 //
96 G4double y0 = 0. * cm; << 114 G4double z0 = 0.5*(GammaRayTelDetector->GetWorldSizeZ());
97 G4double z0 = 0.5 * (detector->GetWorl << 115 G4double x0 = 0.*cm, y0 = 0.*cm;
98 << 116
99 G4ThreeVector pos0; << 117 G4ThreeVector pos0;
100 auto vertex0 = G4ThreeVector(x0, y0, z << 118 G4ThreeVector dir0;
101 auto momentumDirection0 = G4ThreeVecto << 119 G4ThreeVector vertex0 = G4ThreeVector(x0,y0,z0);
102 << 120
103 G4double theta; << 121 dir0 = G4ThreeVector(0.,0.,-1.);
104 G4double phi; << 122
105 G4double y = 0.; << 123 G4double theta, phi, y, f;
106 G4double f = 0.; << 124 G4double theta0=0.;
107 G4double theta0 = 0.; << 125 G4double phi0=0.;
108 G4double phi0 = 0.; << 126
109 << 127 switch(nSourceType) {
110 switch (sourceType) { << 128 case 0:
111 case 0: << 129 particleGun->SetParticlePosition(vertex0);
112 particleGun->SetParticlePosition(v << 130 particleGun->SetParticleMomentumDirection(dir0);
113 particleGun->SetParticleMomentumDi << 131 break;
114 break; << 132 case 1:
115 case 1: << 133 // GS: Generate random position on the 4PIsphere to create a unif. distrib.
116 // GS: Generate random position on << 134 // GS: on the sphere
117 // GS: on the sphere << 135 phi = G4UniformRand() * twopi;
118 phi = G4UniformRand() * twopi; << 136 do {
119 do { << 137 y = G4UniformRand()*1.0;
120 y = G4UniformRand() * 1.0; << 138 theta = G4UniformRand() * pi;
121 theta = G4UniformRand() * pi; << 139 f = std::sin(theta);
122 f = std::sin(theta); << 140 } while (y > f);
123 } while (y > f); << 141 vertex0 = G4ThreeVector(1.,0.,0.);
124 vertex0 = G4ThreeVector(1., 0., 0. << 142 vertex0.setMag(dVertexRadius);
125 vertex0.setMag(vertexRadius); << 143 vertex0.setTheta(theta);
126 vertex0.setTheta(theta); << 144 vertex0.setPhi(phi);
127 vertex0.setPhi(phi); << 145 particleGun->SetParticlePosition(vertex0);
128 particleGun->SetParticlePosition(v << 146
129 << 147 dir0 = G4ThreeVector(1.,0.,0.);
130 momentumDirection0 = G4ThreeVector << 148 do {
131 << 149 phi = G4UniformRand() * twopi;
132 do { << 150 do {
133 phi = G4UniformRand() * twopi; << 151 y = G4UniformRand()*1.0;
134 do { << 152 theta = G4UniformRand() * pi;
135 y = G4UniformRand() * 1.0; << 153 f = std::sin(theta);
136 theta = G4UniformRand() * << 154 } while (y > f);
137 f = std::sin(theta); << 155 dir0.setPhi(phi);
138 } while (y > f); << 156 dir0.setTheta(theta);
139 momentumDirection0.setPhi(phi) << 157 } while (vertex0.dot(dir0) >= -0.7 * vertex0.mag());
140 momentumDirection0.setTheta(th << 158 particleGun->SetParticleMomentumDirection((G4ParticleMomentum)dir0);
141 } while (vertex0.dot(momentumDirec << 159
142 << 160 break;
143 particleGun->SetParticleMomentumDi << 161 case 2:
144 << 162 // GS: Generate random position on the upper semi-sphere z>0 to create a unif. distrib.
145 break; << 163 // GS: on a plane
146 case 2: << 164 phi = G4UniformRand() * twopi;
147 // GS: Generate random position on << 165 do {
148 // GS: on a plane << 166 y = G4UniformRand()*1.0;
149 phi = G4UniformRand() * twopi; << 167 theta = G4UniformRand() * halfpi;
150 << 168 f = std::sin(theta) * std::cos(theta);
151 do { << 169 } while (y > f);
152 y = G4UniformRand() * 1.0; << 170 vertex0 = G4ThreeVector(1.,0.,0.);
153 theta = G4UniformRand() * half << 171
154 f = std::sin(theta) * std::cos << 172 G4double xy = GammaRayTelDetector->GetWorldSizeXY();
155 } while (y > f); << 173 G4double z = GammaRayTelDetector->GetWorldSizeZ();
156 << 174
157 vertex0 = G4ThreeVector(1., 0., 0. << 175 if (dVertexRadius > xy*0.5)
158 << 176 {
159 auto xy = detector->GetWorldSizeXY << 177 G4cout << "vertexRadius too big " << G4endl;
160 auto z = detector->GetWorldSizeZ() << 178 G4cout << "vertexRadius setted to " << xy*0.45 << G4endl;
161 << 179 dVertexRadius = xy*0.45;
162 if (vertexRadius > xy * 0.5) { << 180 }
163 G4cout << "vertexRadius too bi << 181
164 G4cout << "vertexRadius set to << 182 if (dVertexRadius > z*0.5)
165 vertexRadius = xy * 0.45; << 183 {
166 } << 184 G4cout << "vertexRadius too high " << G4endl;
167 << 185 G4cout << "vertexRadius setted to " << z*0.45 << G4endl;
168 if (vertexRadius > z * 0.5) { << 186 dVertexRadius = z*0.45;
169 G4cout << "vertexRadius too hi << 187 }
170 G4cout << "vertexRadius set to << 188
171 vertexRadius = z * 0.45; << 189
172 } << 190 vertex0.setMag(dVertexRadius);
173 << 191 vertex0.setTheta(theta);
174 vertex0.setMag(vertexRadius); << 192 vertex0.setPhi(phi);
175 vertex0.setTheta(theta); << 193
176 vertex0.setPhi(phi); << 194 // GS: Get the user defined direction for the primaries and
177 << 195 // GS: Rotate the random position according to the user defined direction for the particle
178 // GS: Get the user defined direct << 196
179 // GS: Rotate the random position << 197 dir0 = particleGun->GetParticleMomentumDirection();
180 << 198 if (dir0.mag() > 0.001)
181 momentumDirection0 = particleGun-> << 199 {
182 if (momentumDirection0.mag() > 0.0 << 200 theta0 = dir0.theta();
183 theta0 = momentumDirection0.th << 201 phi0 = dir0.phi();
184 phi0 = momentumDirection0.phi( << 202 }
185 } << 203
186 << 204 if (theta0!=0.)
187 if (theta0 != 0.) { << 205 {
188 G4ThreeVector rotationAxis(1., << 206 G4ThreeVector rotationAxis(1.,0.,0.);
189 rotationAxis.setPhi(phi0 + hal << 207 rotationAxis.setPhi(phi0+halfpi);
190 vertex0.rotate(theta0 + pi, ro << 208 vertex0.rotate(theta0+pi,rotationAxis);
191 } << 209 }
192 particleGun->SetParticlePosition(v << 210 particleGun->SetParticlePosition(vertex0);
193 break; << 211 break;
194 } << 212 }
195 << 213
196 constexpr auto INITIAL_PARTICLE_ENERGY << 214
197 G4double particleEnergy = INITIAL_PART << 215 G4double pEnergy;
198 << 216
199 switch (spectrumType) { << 217 switch(nSpectrumType) {
200 case 0: // Uniform energy (1 GeV - 10 << 218 case 0:
201 y = G4UniformRand(); << 219 break;
202 particleEnergy = y * 9.0 * GeV + 1 << 220 case 1:
203 G4cout << "Particle energy: " << p << 221 break;
204 break; << 222 case 2:
205 case 1: // Logarithmic energy << 223 do {
206 y = G4UniformRand(); << 224 y = G4UniformRand()*100000.0;
207 particleEnergy = std::pow(10, y) * << 225 pEnergy = G4UniformRand() * 10. * GeV;
208 G4cout << "Particle energy: " << p << 226 f = std::pow(pEnergy * (1/GeV), -4.);
209 break; << 227 } while (y > f);
210 case 2: // Power law (-4) << 228
211 do { << 229 particleGun->SetParticleEnergy(pEnergy);
212 y = G4UniformRand() * 100000.0 << 230
213 particleEnergy = G4UniformRand << 231 break;
214 f = std::pow(particleEnergy * << 232 case 3:
215 } while (y > f); << 233 break;
216 // particleGun->SetParticleEnergy( << 234 }
217 break; << 235
218 case 3: // Monochromatic << 236 particleGun->GeneratePrimaryVertex(anEvent);
219 particleEnergy = particleGun->GetP <<
220 // 100 MeV; <<
221 G4cout << "Particle energy: " << p <<
222 break; <<
223 } <<
224 particleGun->SetParticleEnergy(particl <<
225 G4cout << "Particle: " << particleGun- <<
226 particleGun->GeneratePrimaryVertex(eve <<
227 } else { <<
228 particleSource->GeneratePrimaryVertex( <<
229 } 237 }
>> 238 else
>> 239 {
>> 240 particleSource->GeneratePrimaryVertex(anEvent);
>> 241 }
>> 242
230 } 243 }
>> 244
>> 245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 246
>> 247
>> 248
>> 249
>> 250
>> 251
>> 252
>> 253
231 254