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