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 // Write these values into the analysis if needed. Have to be written separately on change. >> 93 IORTAnalysisManager::GetInstance()-> >> 94 setBeamMetaData(meanKineticEnergy, sigmaEnergy); >> 95 >> 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 // Increment the event counter >> 135 IORTAnalysisManager::GetInstance()->startNewEvent(); >> 136 125 // **************************************** 137 // **************************************** 126 // Set the beam angular apread 138 // Set the beam angular apread 127 // and spot size 139 // and spot size 128 // beam spot size 140 // beam spot size 129 // **************************************** 141 // **************************************** 130 142 131 // Set the position of the primary particles 143 // Set the position of the primary particles 132 G4double x = X0; 144 G4double x = X0; 133 G4double y = Y0; 145 G4double y = Y0; 134 G4double z = Z0; 146 G4double z = Z0; 135 147 136 if ( sigmaY > 0.0 ) 148 if ( sigmaY > 0.0 ) 137 { 149 { 138 y += G4RandGauss::shoot( Y0, sigmaY ); 150 y += G4RandGauss::shoot( Y0, sigmaY ); 139 } 151 } 140 152 141 if ( sigmaZ > 0.0 ) 153 if ( sigmaZ > 0.0 ) 142 { 154 { 143 z += G4RandGauss::shoot( Z0, sigmaZ ); 155 z += G4RandGauss::shoot( Z0, sigmaZ ); 144 } 156 } 145 157 146 particleGun -> SetParticlePosition(G4ThreeVe 158 particleGun -> SetParticlePosition(G4ThreeVector( x , y , z ) ); 147 159 148 // ***************************************** 160 // ******************************************** 149 // Set the beam energy and energy spread 161 // Set the beam energy and energy spread 150 // ***************************************** 162 // ******************************************** 151 163 152 G4double kineticEnergy = G4RandGauss::shoot( 164 G4double kineticEnergy = G4RandGauss::shoot( meanKineticEnergy, sigmaEnergy ); 153 particleGun -> SetParticleEnergy ( kineticEn 165 particleGun -> SetParticleEnergy ( kineticEnergy ); 154 166 155 // Set the direction of the primary particle 167 // Set the direction of the primary particles 156 168 157 169 158 /* 170 /* 159 G4double momentumX = 1.0; 171 G4double momentumX = 1.0; 160 G4double momentumY = 0.0; 172 G4double momentumY = 0.0; 161 G4double momentumZ = 0.0; 173 G4double momentumZ = 0.0; 162 174 163 if ( sigmaMomentumY > 0.0 ) 175 if ( sigmaMomentumY > 0.0 ) 164 { 176 { 165 momentumY += G4RandGauss::shoot( 0., sig 177 momentumY += G4RandGauss::shoot( 0., sigmaMomentumY ); 166 } 178 } 167 if ( sigmaMomentumZ > 0.0 ) 179 if ( sigmaMomentumZ > 0.0 ) 168 { 180 { 169 momentumZ += G4RandGauss::shoot( 0., sig 181 momentumZ += G4RandGauss::shoot( 0., sigmaMomentumZ ); 170 } 182 } 171 183 172 particleGun -> SetParticleMomentumDirection( 184 particleGun -> SetParticleMomentumDirection( G4ThreeVector(momentumX,momentumY,momentumZ) ); 173 185 174 */ 186 */ 175 187 176 188 177 G4double Mx; 189 G4double Mx; 178 G4double My; 190 G4double My; 179 G4double Mz; 191 G4double Mz; 180 G4double condizione; 192 G4double condizione; 181 193 182 while (true) { 194 while (true) { 183 195 184 //Mx = CLHEP::RandFlat::shoot(0.9,1); 196 //Mx = CLHEP::RandFlat::shoot(0.9,1); 185 //My = CLHEP::RandFlat::shoot(-0.1,0.1); 197 //My = CLHEP::RandFlat::shoot(-0.1,0.1); 186 //Mz = CLHEP::RandFlat::shoot(-0.1,0.1); 198 //Mz = CLHEP::RandFlat::shoot(-0.1,0.1); 187 199 188 Mx = CLHEP::RandFlat::shoot(0.7,1); 200 Mx = CLHEP::RandFlat::shoot(0.7,1); 189 My = CLHEP::RandFlat::shoot(-0.3,0.3); // r 201 My = CLHEP::RandFlat::shoot(-0.3,0.3); // ranges good for 0<Theta<20 190 Mz = CLHEP::RandFlat::shoot(-0.3,0.3); 202 Mz = CLHEP::RandFlat::shoot(-0.3,0.3); 191 203 192 condizione = std::sqrt(Mx*Mx + My*My + Mz*Mz 204 condizione = std::sqrt(Mx*Mx + My*My + Mz*Mz); 193 205 194 206 195 if (condizione < 1) { 207 if (condizione < 1) { 196 Mx = Mx/condizione; 208 Mx = Mx/condizione; 197 My = My/condizione; 209 My = My/condizione; 198 Mz = Mz/condizione; 210 Mz = Mz/condizione; 199 211 200 212 201 if (Mx > std::cos(Theta)) { 213 if (Mx > std::cos(Theta)) { 202 break; 214 break; 203 } 215 } 204 } 216 } 205 } 217 } 206 218 207 219 208 particleGun -> SetParticleMomentumDirection( 220 particleGun -> SetParticleMomentumDirection( G4ThreeVector(Mx,My,Mz) ); 209 221 210 222 211 // Generate a primary particle 223 // Generate a primary particle 212 particleGun -> GeneratePrimaryVertex( anEven 224 particleGun -> GeneratePrimaryVertex( anEvent ); 213 } 225 } 214 226 215 void IORTPrimaryGeneratorAction::SetmeanKineti 227 void IORTPrimaryGeneratorAction::SetmeanKineticEnergy (G4double val ) 216 { 228 { 217 meanKineticEnergy = val; << 229 meanKineticEnergy = val; 218 G4cout << "The mean Kinetic energy of the in << 230 219 << meanKineticEnergy/MeV << G4endl; << 231 // Update the beam-data in the analysis manager >> 232 IORTAnalysisManager::GetInstance()-> >> 233 setBeamMetaData(meanKineticEnergy, sigmaEnergy); >> 234 >> 235 220 } 236 } 221 237 222 void IORTPrimaryGeneratorAction::SetsigmaEnerg 238 void IORTPrimaryGeneratorAction::SetsigmaEnergy (G4double val ) 223 { 239 { 224 sigmaEnergy = val; << 240 sigmaEnergy = val; 225 G4cout << "The sigma of the kinetic energy of << 241 // Update the sigmaenergy in the metadata. 226 << sigmaEnergy/MeV << G4endl; << 242 IORTAnalysisManager::GetInstance()-> >> 243 setBeamMetaData(meanKineticEnergy, sigmaEnergy); 227 } 244 } 228 245 229 void IORTPrimaryGeneratorAction::SetXposition 246 void IORTPrimaryGeneratorAction::SetXposition (G4double val ) 230 { X0 = val;} 247 { X0 = val;} 231 248 232 void IORTPrimaryGeneratorAction::SetYposition 249 void IORTPrimaryGeneratorAction::SetYposition (G4double val ) 233 { Y0 = val;} 250 { Y0 = val;} 234 251 235 void IORTPrimaryGeneratorAction::SetZposition 252 void IORTPrimaryGeneratorAction::SetZposition (G4double val ) 236 { Z0 = val;} 253 { Z0 = val;} 237 254 238 void IORTPrimaryGeneratorAction::SetsigmaY (G4 255 void IORTPrimaryGeneratorAction::SetsigmaY (G4double val ) 239 { sigmaY = val;} 256 { sigmaY = val;} 240 257 241 void IORTPrimaryGeneratorAction::SetsigmaZ (G4 258 void IORTPrimaryGeneratorAction::SetsigmaZ (G4double val ) 242 { sigmaZ = val;} 259 { sigmaZ = val;} 243 260 244 /* 261 /* 245 void IORTPrimaryGeneratorAction::SetsigmaMomen 262 void IORTPrimaryGeneratorAction::SetsigmaMomentumY (G4double val ) 246 { sigmaMomentumY = val;} 263 { sigmaMomentumY = val;} 247 264 248 void IORTPrimaryGeneratorAction::SetsigmaMomen 265 void IORTPrimaryGeneratorAction::SetsigmaMomentumZ (G4double val ) 249 { sigmaMomentumZ = val;} 266 { sigmaMomentumZ = val;} 250 */ 267 */ 251 268 252 void IORTPrimaryGeneratorAction::SetTheta (G4d 269 void IORTPrimaryGeneratorAction::SetTheta (G4double val ) 253 { Theta = val;} 270 { Theta = val;} 254 271 255 272 256 G4double IORTPrimaryGeneratorAction::GetmeanKi 273 G4double IORTPrimaryGeneratorAction::GetmeanKineticEnergy(void) 257 { return meanKineticEnergy;} 274 { return meanKineticEnergy;} 258 275 259 276