| 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 // This is the *BASIC* version of IORT, a Gean 26 // This is the *BASIC* version of IORT, a Geant4-based application
27 // 27 //
28 // Main Authors: G.Russo(a,b), C.Casarino*(c), 28 // Main Authors: G.Russo(a,b), C.Casarino*(c), G.C. Candiano(c), G.A.P. Cirrone(d), F.Romano(d)
29 // Contributor Authors: S.Guatelli(e) 29 // Contributor Authors: S.Guatelli(e)
30 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d 30 // Past Authors: G.Arnetta(c), S.E.Mazzaglia(d)
31 // 31 //
32 // (a) Fondazione Istituto San Raffaele G.Gi 32 // (a) Fondazione Istituto San Raffaele G.Giglio, Cefalù, Italy
33 // (b) IBFM-CNR , Segrate (Milano), Italy 33 // (b) IBFM-CNR , Segrate (Milano), Italy
34 // (c) LATO (Laboratorio di Tecnologie Oncol 34 // (c) LATO (Laboratorio di Tecnologie Oncologiche), Cefalù, Italy
35 // (d) Laboratori Nazionali del Sud of the I 35 // (d) Laboratori Nazionali del Sud of the INFN, Catania, Italy
36 // (e) University of Wollongong, Australia << 36 // (e) University of Wallongong, Australia
37 // 37 //
38 // *Corresponding author, email to carlo.cas 38 // *Corresponding author, email to carlo.casarino@polooncologicocefalu.it
39 ////////////////////////////////////////////// 39 //////////////////////////////////////////////////////////////////////////////////////////////
40 40
41 #include "G4SystemOfUnits.hh" << 41 #include <CLHEP/Units/SystemOfUnits.h>
>> 42
42 #include "IORTPrimaryGeneratorAction.hh" 43 #include "IORTPrimaryGeneratorAction.hh"
43 #include "IORTPrimaryGeneratorMessenger.hh" 44 #include "IORTPrimaryGeneratorMessenger.hh"
44 45
45 #include "globals.hh" 46 #include "globals.hh"
46 #include "G4Event.hh" 47 #include "G4Event.hh"
47 #include "G4ParticleGun.hh" 48 #include "G4ParticleGun.hh"
48 #include "G4ParticleTable.hh" 49 #include "G4ParticleTable.hh"
49 #include "G4ParticleDefinition.hh" 50 #include "G4ParticleDefinition.hh"
50 #include "Randomize.hh" 51 #include "Randomize.hh"
51 << 52
>> 53 #include "IORTAnalysisManager.hh"
>> 54
52 IORTPrimaryGeneratorAction::IORTPrimaryGenerat 55 IORTPrimaryGeneratorAction::IORTPrimaryGeneratorAction()
53 { 56 {
54 // Define the messenger 57 // Define the messenger
55 gunMessenger = new IORTPrimaryGeneratorMesse 58 gunMessenger = new IORTPrimaryGeneratorMessenger(this);
56 59
57 particleGun = new G4ParticleGun(); 60 particleGun = new G4ParticleGun();
58 61
59 SetDefaultPrimaryParticle(); 62 SetDefaultPrimaryParticle();
60 } 63 }
61 64
62 IORTPrimaryGeneratorAction::~IORTPrimaryGenera 65 IORTPrimaryGeneratorAction::~IORTPrimaryGeneratorAction()
63 { 66 {
64 delete particleGun; 67 delete particleGun;
65 68
66 delete gunMessenger; 69 delete gunMessenger;
67 } 70 }
68 71
69 void IORTPrimaryGeneratorAction::SetDefaultPri 72 void IORTPrimaryGeneratorAction::SetDefaultPrimaryParticle()
70 { 73 {
71 // **************************** 74 // ****************************
72 // Default primary particle 75 // Default primary particle
73 // **************************** 76 // ****************************
74 77
75 // Define primary particles: electrons // pr 78 // Define primary particles: electrons // protons
76 G4ParticleTable* particleTable = G4ParticleT 79 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
77 G4ParticleDefinition* particle = particleTab 80 G4ParticleDefinition* particle = particleTable -> FindParticle("e-"); // ("proton")
78 particleGun -> SetParticleDefinition(particl 81 particleGun -> SetParticleDefinition(particle);
79 82
80 // Define the energy of primary particles: 83 // Define the energy of primary particles:
81 // gaussian distribution with mean energy = 84 // gaussian distribution with mean energy = 10.0 *MeV
82 // and sigma = 400.0 *keV 85 // and sigma = 400.0 *keV
83 G4double defaultMeanKineticEnergy = 10.0 *CL 86 G4double defaultMeanKineticEnergy = 10.0 *CLHEP::MeV;
84 meanKineticEnergy = defaultMeanKineticEnergy 87 meanKineticEnergy = defaultMeanKineticEnergy;
85 88
86 G4double defaultsigmaEnergy = 100.0 *CLHEP:: 89 G4double defaultsigmaEnergy = 100.0 *CLHEP::keV;
87 sigmaEnergy = defaultsigmaEnergy; 90 sigmaEnergy = defaultsigmaEnergy;
88 << 91
>> 92 #ifdef G4ANALYSIS_USE_ROOT
>> 93 // Write these values into the analysis if needed. Have to be written separately on change.
>> 94 IORTAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
>> 95 #endif
>> 96
89 // Define the parameters of the initial posi 97 // Define the parameters of the initial position:
90 // the y, z coordinates have a gaussian dist << 98 // the y, z coordinates have a gaussian distribution
>> 99
91 100
92 G4double defaultX0 = -862.817 *CLHEP::mm; 101 G4double defaultX0 = -862.817 *CLHEP::mm;
93 X0 = defaultX0; 102 X0 = defaultX0;
94 103
95 G4double defaultY0 = 0.0 *CLHEP::mm; 104 G4double defaultY0 = 0.0 *CLHEP::mm;
96 Y0 = defaultY0; 105 Y0 = defaultY0;
97 106
98 G4double defaultZ0 = 0.0 *CLHEP::mm; 107 G4double defaultZ0 = 0.0 *CLHEP::mm;
99 Z0 = defaultZ0; 108 Z0 = defaultZ0;
100 109
101 G4double defaultsigmaY = 1. *CLHEP::mm; 110 G4double defaultsigmaY = 1. *CLHEP::mm;
102 sigmaY = defaultsigmaY; 111 sigmaY = defaultsigmaY;
103 112
104 G4double defaultsigmaZ = 1. *CLHEP::mm; 113 G4double defaultsigmaZ = 1. *CLHEP::mm;
105 sigmaZ = defaultsigmaZ; 114 sigmaZ = defaultsigmaZ;
106 115
107 // Define the parameters of the momentum of 116 // Define the parameters of the momentum of primary particles:
108 // The momentum along the y and z axis has a 117 // The momentum along the y and z axis has a gaussian distribution
109 118
110 /* 119 /*
111 G4double defaultsigmaMomentumY = 0.0; 120 G4double defaultsigmaMomentumY = 0.0;
112 sigmaMomentumY = defaultsigmaMomentumY; 121 sigmaMomentumY = defaultsigmaMomentumY;
113 122
114 G4double defaultsigmaMomentumZ = 0.0; 123 G4double defaultsigmaMomentumZ = 0.0;
115 sigmaMomentumZ = defaultsigmaMomentumZ; 124 sigmaMomentumZ = defaultsigmaMomentumZ;
116 */ 125 */
117 126
118 G4double defaultTheta = 6.0 *CLHEP::deg; 127 G4double defaultTheta = 6.0 *CLHEP::deg;
119 Theta = defaultTheta; 128 Theta = defaultTheta;
120 129
121 } 130 }
122 131
123 void IORTPrimaryGeneratorAction::GeneratePrima 132 void IORTPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
124 { 133 {
>> 134 #ifdef G4ANALYSIS_USE_ROOT
>> 135 // Increment the event counter
>> 136 IORTAnalysisManager::GetInstance()->startNewEvent();
>> 137 #endif
>> 138
125 // **************************************** 139 // ****************************************
126 // Set the beam angular apread 140 // Set the beam angular apread
127 // and spot size 141 // and spot size
128 // beam spot size 142 // beam spot size
129 // **************************************** 143 // ****************************************
130 144
131 // Set the position of the primary particles 145 // Set the position of the primary particles
132 G4double x = X0; 146 G4double x = X0;
133 G4double y = Y0; 147 G4double y = Y0;
134 G4double z = Z0; 148 G4double z = Z0;
135 149
136 if ( sigmaY > 0.0 ) 150 if ( sigmaY > 0.0 )
137 { 151 {
138 y += G4RandGauss::shoot( Y0, sigmaY ); 152 y += G4RandGauss::shoot( Y0, sigmaY );
139 } 153 }
140 154
141 if ( sigmaZ > 0.0 ) 155 if ( sigmaZ > 0.0 )
142 { 156 {
143 z += G4RandGauss::shoot( Z0, sigmaZ ); 157 z += G4RandGauss::shoot( Z0, sigmaZ );
144 } 158 }
145 159
146 particleGun -> SetParticlePosition(G4ThreeVe 160 particleGun -> SetParticlePosition(G4ThreeVector( x , y , z ) );
147 161
148 // ***************************************** 162 // ********************************************
149 // Set the beam energy and energy spread 163 // Set the beam energy and energy spread
150 // ***************************************** 164 // ********************************************
151 165
152 G4double kineticEnergy = G4RandGauss::shoot( 166 G4double kineticEnergy = G4RandGauss::shoot( meanKineticEnergy, sigmaEnergy );
153 particleGun -> SetParticleEnergy ( kineticEn 167 particleGun -> SetParticleEnergy ( kineticEnergy );
154 168
155 // Set the direction of the primary particle 169 // Set the direction of the primary particles
156 170
157 171
158 /* 172 /*
159 G4double momentumX = 1.0; 173 G4double momentumX = 1.0;
160 G4double momentumY = 0.0; 174 G4double momentumY = 0.0;
161 G4double momentumZ = 0.0; 175 G4double momentumZ = 0.0;
162 176
163 if ( sigmaMomentumY > 0.0 ) 177 if ( sigmaMomentumY > 0.0 )
164 { 178 {
165 momentumY += G4RandGauss::shoot( 0., sig 179 momentumY += G4RandGauss::shoot( 0., sigmaMomentumY );
166 } 180 }
167 if ( sigmaMomentumZ > 0.0 ) 181 if ( sigmaMomentumZ > 0.0 )
168 { 182 {
169 momentumZ += G4RandGauss::shoot( 0., sig 183 momentumZ += G4RandGauss::shoot( 0., sigmaMomentumZ );
170 } 184 }
171 185
172 particleGun -> SetParticleMomentumDirection( 186 particleGun -> SetParticleMomentumDirection( G4ThreeVector(momentumX,momentumY,momentumZ) );
173 187
174 */ 188 */
175 189
176 190
177 G4double Mx; 191 G4double Mx;
178 G4double My; 192 G4double My;
179 G4double Mz; 193 G4double Mz;
180 G4double condizione; 194 G4double condizione;
181 195
182 while (true) { 196 while (true) {
183 197
184 //Mx = CLHEP::RandFlat::shoot(0.9,1); 198 //Mx = CLHEP::RandFlat::shoot(0.9,1);
185 //My = CLHEP::RandFlat::shoot(-0.1,0.1); 199 //My = CLHEP::RandFlat::shoot(-0.1,0.1);
186 //Mz = CLHEP::RandFlat::shoot(-0.1,0.1); 200 //Mz = CLHEP::RandFlat::shoot(-0.1,0.1);
187 201
188 Mx = CLHEP::RandFlat::shoot(0.7,1); 202 Mx = CLHEP::RandFlat::shoot(0.7,1);
189 My = CLHEP::RandFlat::shoot(-0.3,0.3); // r 203 My = CLHEP::RandFlat::shoot(-0.3,0.3); // ranges good for 0<Theta<20
190 Mz = CLHEP::RandFlat::shoot(-0.3,0.3); 204 Mz = CLHEP::RandFlat::shoot(-0.3,0.3);
191 205
192 condizione = std::sqrt(Mx*Mx + My*My + Mz*Mz 206 condizione = std::sqrt(Mx*Mx + My*My + Mz*Mz);
193 207
194 208
195 if (condizione < 1) { 209 if (condizione < 1) {
196 Mx = Mx/condizione; 210 Mx = Mx/condizione;
197 My = My/condizione; 211 My = My/condizione;
198 Mz = Mz/condizione; 212 Mz = Mz/condizione;
199 213
200 214
201 if (Mx > std::cos(Theta)) { 215 if (Mx > std::cos(Theta)) {
202 break; 216 break;
203 } 217 }
204 } 218 }
205 } 219 }
206 220
207 221
208 particleGun -> SetParticleMomentumDirection( 222 particleGun -> SetParticleMomentumDirection( G4ThreeVector(Mx,My,Mz) );
209 223
210 224
211 // Generate a primary particle 225 // Generate a primary particle
212 particleGun -> GeneratePrimaryVertex( anEven 226 particleGun -> GeneratePrimaryVertex( anEvent );
213 } 227 }
214 228
215 void IORTPrimaryGeneratorAction::SetmeanKineti 229 void IORTPrimaryGeneratorAction::SetmeanKineticEnergy (G4double val )
216 { 230 {
217 meanKineticEnergy = val; << 231 meanKineticEnergy = val;
218 G4cout << "The mean Kinetic energy of the in << 232 #ifdef G4ANALYSIS_USE_ROOT
219 << meanKineticEnergy/MeV << G4endl; << 233 // Update the beam-data in the analysis manager
>> 234 IORTAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
>> 235 #endif
>> 236
220 } 237 }
221 238
222 void IORTPrimaryGeneratorAction::SetsigmaEnerg 239 void IORTPrimaryGeneratorAction::SetsigmaEnergy (G4double val )
223 { 240 {
224 sigmaEnergy = val; << 241 sigmaEnergy = val;
225 G4cout << "The sigma of the kinetic energy of << 242 #ifdef G4ANALYSIS_USE_ROOT
226 << sigmaEnergy/MeV << G4endl; << 243 // Update the sigmaenergy in the metadata.
>> 244 IORTAnalysisManager::GetInstance()->setBeamMetaData(meanKineticEnergy, sigmaEnergy);
>> 245 #endif
227 } 246 }
228 247
229 void IORTPrimaryGeneratorAction::SetXposition 248 void IORTPrimaryGeneratorAction::SetXposition (G4double val )
230 { X0 = val;} 249 { X0 = val;}
231 250
232 void IORTPrimaryGeneratorAction::SetYposition 251 void IORTPrimaryGeneratorAction::SetYposition (G4double val )
233 { Y0 = val;} 252 { Y0 = val;}
234 253
235 void IORTPrimaryGeneratorAction::SetZposition 254 void IORTPrimaryGeneratorAction::SetZposition (G4double val )
236 { Z0 = val;} 255 { Z0 = val;}
237 256
238 void IORTPrimaryGeneratorAction::SetsigmaY (G4 257 void IORTPrimaryGeneratorAction::SetsigmaY (G4double val )
239 { sigmaY = val;} 258 { sigmaY = val;}
240 259
241 void IORTPrimaryGeneratorAction::SetsigmaZ (G4 260 void IORTPrimaryGeneratorAction::SetsigmaZ (G4double val )
242 { sigmaZ = val;} 261 { sigmaZ = val;}
243 262
244 /* 263 /*
245 void IORTPrimaryGeneratorAction::SetsigmaMomen 264 void IORTPrimaryGeneratorAction::SetsigmaMomentumY (G4double val )
246 { sigmaMomentumY = val;} 265 { sigmaMomentumY = val;}
247 266
248 void IORTPrimaryGeneratorAction::SetsigmaMomen 267 void IORTPrimaryGeneratorAction::SetsigmaMomentumZ (G4double val )
249 { sigmaMomentumZ = val;} 268 { sigmaMomentumZ = val;}
250 */ 269 */
251 270
252 void IORTPrimaryGeneratorAction::SetTheta (G4d 271 void IORTPrimaryGeneratorAction::SetTheta (G4double val )
253 { Theta = val;} 272 { Theta = val;}
254 273
255 274
256 G4double IORTPrimaryGeneratorAction::GetmeanKi 275 G4double IORTPrimaryGeneratorAction::GetmeanKineticEnergy(void)
257 { return meanKineticEnergy;} 276 { return meanKineticEnergy;}
258 277
259 278