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 example is provided by the Geant4-DNA 26 // This example is provided by the Geant4-DNA collaboration 27 // Any report or published results obtained us 27 // Any report or published results obtained using the Geant4-DNA software 28 // shall cite the following Geant4-DNA collabo 28 // shall cite the following Geant4-DNA collaboration publication: 29 // Med. Phys. 37 (2010) 4692-4708 29 // Med. Phys. 37 (2010) 4692-4708 30 // and papers 30 // and papers 31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507 32 // O. Belov et al. Physica Medica 32 (2016) 15 32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520 33 // The Geant4-DNA web site is available at htt 33 // The Geant4-DNA web site is available at http://geant4-dna.org 34 // << 34 // 35 // ------------------------------------------- 35 // ------------------------------------------------------------------- 36 // November 2016 36 // November 2016 37 // ------------------------------------------- 37 // ------------------------------------------------------------------- 38 // 38 // 39 /// \file PrimaryGeneratorAction.cc << 39 /// \file PrimaryGeneratorAction.cc 40 /// \brief Implementation of the PrimaryGenera 40 /// \brief Implementation of the PrimaryGeneratorAction class 41 41 42 #include "PrimaryGeneratorAction.hh" 42 #include "PrimaryGeneratorAction.hh" 43 << 44 #include "G4SystemOfUnits.hh" 43 #include "G4SystemOfUnits.hh" 45 // 44 // 46 #include "CommandLineParser.hh" 45 #include "CommandLineParser.hh" 47 << 46 #include "G4ParticleTable.hh" 48 #include "G4Box.hh" << 47 #include "G4ParticleGun.hh" >> 48 #include "G4ParticleDefinition.hh" >> 49 #include "G4PhysicalConstants.hh" >> 50 #include "Randomize.hh" 49 #include "G4Event.hh" 51 #include "G4Event.hh" >> 52 #include "G4Box.hh" >> 53 #include "G4Orb.hh" 50 #include "G4LogicalVolume.hh" 54 #include "G4LogicalVolume.hh" 51 #include "G4LogicalVolumeStore.hh" 55 #include "G4LogicalVolumeStore.hh" 52 #include "G4Orb.hh" << 53 #include "G4ParticleDefinition.hh" << 54 #include "G4ParticleGun.hh" << 55 #include "G4PhysicalConstants.hh" << 56 #include "G4PhysicalVolumeStore.hh" << 57 #include "G4Proton.hh" << 58 #include "G4VPhysicalVolume.hh" 56 #include "G4VPhysicalVolume.hh" 59 #include "Randomize.hh" << 57 #include "G4PhysicalVolumeStore.hh" 60 58 61 using namespace G4DNAPARSER; << 59 using namespace G4DNAPARSER ; 62 60 63 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 62 65 PrimaryGeneratorAction::PrimaryGeneratorAction << 63 PrimaryGeneratorAction::PrimaryGeneratorAction(): >> 64 G4VUserPrimaryGeneratorAction() 66 { 65 { 67 G4int n_particle = 1; 66 G4int n_particle = 1; 68 fpParticleGun = new G4ParticleGun(n_particle << 67 fpParticleGun = new G4ParticleGun(n_particle); 69 // 68 // 70 G4ParticleDefinition* particle = G4Proton::P << 69 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); >> 70 G4ParticleDefinition* particle = particleTable->FindParticle("proton"); 71 fpParticleGun->SetParticleDefinition(particl 71 fpParticleGun->SetParticleDefinition(particle); 72 // default gun parameters 72 // default gun parameters 73 fpParticleGun->SetParticleEnergy(10. * MeV); << 73 fpParticleGun->SetParticleEnergy(10.*MeV); 74 fpParticleGun->SetParticleMomentumDirection( << 74 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); 75 fpParticleGun->SetParticlePosition(G4ThreeVe << 75 fpParticleGun->SetParticlePosition(G4ThreeVector(0.,0.,0.*um)); 76 } 76 } 77 77 78 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 79 79 80 PrimaryGeneratorAction::~PrimaryGeneratorActio 80 PrimaryGeneratorAction::~PrimaryGeneratorAction() 81 { 81 { 82 delete fpParticleGun; 82 delete fpParticleGun; 83 } 83 } 84 84 85 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 86 87 void PrimaryGeneratorAction::GeneratePrimaries 87 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 88 { 88 { 89 // Initial kinetic energy of particles beam << 89 90 << 90 // Initial kinetic energy of particles beam as Gaussion distribution! >> 91 91 // G4double Ekin = 10.*MeV; 92 // G4double Ekin = 10.*MeV; 92 // G4double deltaEkin = 9.0955e-5*MeV; << 93 // G4double deltaEkin = 9.0955e-5*MeV; 93 // fpParticleGun->SetParticleEnergy(G4RandGa << 94 // fpParticleGun->SetParticleEnergy(G4RandGauss::shoot(Ekin,deltaEkin)); 94 << 95 95 // In order to avoid dependence of PrimaryGe 96 // In order to avoid dependence of PrimaryGeneratorAction 96 // on DetectorConstruction class we get worl 97 // on DetectorConstruction class we get world volume 97 // from G4LogicalVolumeStore 98 // from G4LogicalVolumeStore 98 99 99 // We included three options for particles d 100 // We included three options for particles direction: 100 << 101 101 G4double mediumRadius = 0.; << 102 G4double mediumRadius = 0.; 102 G4double boundingXHalfLength = 0.; 103 G4double boundingXHalfLength = 0.; 103 G4double boundingYHalfLength = 0.; 104 G4double boundingYHalfLength = 0.; 104 G4double boundingZHalfLength = 0.; << 105 G4double boundingZHalfLength = 0.; 105 G4LogicalVolume* mediumLV = G4LogicalVolumeS << 106 G4LogicalVolume* mediumLV 106 G4LogicalVolume* boundingLV = G4LogicalVolum << 107 = G4LogicalVolumeStore::GetInstance()->GetVolume("Medium"); 107 G4Orb* mediumSphere = 0; << 108 G4LogicalVolume* boundingLV 108 G4Box* boundingSlice = 0; << 109 = G4LogicalVolumeStore::GetInstance()->GetVolume("BoundingSlice"); 109 if (mediumLV && boundingLV) { << 110 G4Orb* mediumSphere = 0; 110 mediumSphere = dynamic_cast<G4Orb*>(medium << 111 G4Box* boundingSlice = 0; 111 boundingSlice = dynamic_cast<G4Box*>(bound << 112 if ( mediumLV && boundingLV) 112 } << 113 { 113 if (mediumSphere && boundingSlice) { << 114 mediumSphere = dynamic_cast< G4Orb*>(mediumLV->GetSolid()); 114 mediumRadius = mediumSphere->GetRadius(); << 115 boundingSlice = dynamic_cast< G4Box*>(boundingLV->GetSolid()); 115 boundingXHalfLength = boundingSlice->GetXH << 116 } 116 boundingYHalfLength = boundingSlice->GetYH << 117 if ( mediumSphere && boundingSlice) >> 118 { >> 119 mediumRadius = mediumSphere->GetRadius(); >> 120 boundingXHalfLength = boundingSlice->GetXHalfLength(); >> 121 boundingYHalfLength = boundingSlice->GetYHalfLength(); 117 boundingZHalfLength = boundingSlice->GetZH 122 boundingZHalfLength = boundingSlice->GetZHalfLength(); 118 << 123 119 /// a) Partilces directed to "square" on t << 124 /// a) Partilces directed to "square" on the XY plane of bounding slice 120 /// (or YZ, XZ) << 125 /// (or YZ, XZ) 121 << 126 122 if (CommandLineParser::GetParser()->GetCom << 127 if ( CommandLineParser::GetParser()->GetCommandIfActive("-sXY")) 123 // INITIAL BEAM DIVERGENCE << 128 { 124 fpParticleGun->SetParticleMomentumDirect << 129 //G4cerr << " Initial beam position uniformly spread on a square! " 125 // // INITIAL BEAM POSITION << 130 // << G4endl; 126 fpParticleGun->SetParticlePosition(G4Thr << 131 // INITIAL BEAM DIVERGENCE 127 CLHEP::RandFlat::shoot(-boundingXHalfL << 132 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); 128 CLHEP::RandFlat::shoot(-boundingYHalfL << 133 // // INITIAL BEAM POSITION 129 } << 134 fpParticleGun->SetParticlePosition(G4ThreeVector( 130 << 135 CLHEP::RandFlat::shoot(-boundingXHalfLength,boundingXHalfLength), 131 /// b) Partilces directed to "disk" on the << 136 CLHEP::RandFlat::shoot(-boundingYHalfLength,boundingYHalfLength), 132 /// bounding slice (or YZ, XZ) << 137 -mediumRadius)); 133 << 138 // Surface area of a square 134 else if (CommandLineParser::GetParser()->G << 139 // fGunArea = 4*boundingXHalfLength*boundingYHalfLength ; 135 fpParticleGun->SetParticleMomentumDirect << 140 //G4cerr << " Particle Fluence Area on a Square (um2) = " 136 G4double x0, y0, z0; << 141 // <<fGunArea / (um*um)<< G4endl; 137 x0 = 100. * mm; << 142 } 138 y0 = 100. * mm; << 143 139 z0 = -mediumRadius; // mediumRadius; << 144 /// b) Partilces directed to "disk" on the XY plane of 140 while (!(std::sqrt(x0 * x0 + y0 * y0) <= << 145 /// bounding slice (or YZ, XZ) 141 x0 = CLHEP::RandFlat::shoot(-mediumRad << 146 142 y0 = CLHEP::RandFlat::shoot(-mediumRad << 147 else if ( CommandLineParser::GetParser()->GetCommandIfActive("-dXY")) 143 } << 148 { 144 fpParticleGun->SetParticlePosition(G4Thr << 149 //G4cerr << " Initial beam position uniformly spread on a disk! " 145 } << 150 // << G4endl; 146 << 151 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.)); 147 /// c) Partilces directed towards the boun << 152 G4double x0,y0,z0; 148 // Select a starting position on a sphere << 153 x0 = 100.*mm; 149 // target volume and neuron morphology << 154 y0 = 100.*mm; 150 else { << 155 z0 = -mediumRadius; // mediumRadius; 151 G4double cosTheta = 2. * G4UniformRand() << 156 while (! (std::sqrt(x0*x0+y0*y0)<= mediumRadius) ) 152 G4double sinTheta = std::sqrt(1. - cosTh << 157 { 153 G4double phi = twopi * G4UniformRand(); << 158 x0 = CLHEP::RandFlat::shoot(-mediumRadius,mediumRadius); 154 G4ThreeVector positionStart(mediumRadius << 159 y0 = CLHEP::RandFlat::shoot(-mediumRadius,mediumRadius); 155 mediumRadius << 160 } 156 fpParticleGun->SetParticlePosition(posit << 161 fpParticleGun->SetParticlePosition(G4ThreeVector(x0, y0,z0)); 157 // To compute the direction, select a po << 162 // Surface area of a disk 158 G4ThreeVector positionDir(boundingXHalfL << 163 // fGunArea = pi*mediumRadius*mediumRadius ; 159 boundingYHalfL << 164 // G4cerr << " Particle Fluence Area on a Disk (um2) = " 160 boundingZHalfL << 165 // <<fGunArea / (um*um)<< G4endl; 161 fpParticleGun->SetParticleMomentumDirect << 166 162 // Surface area of sphere << 167 } 163 fGunArea = 4. * pi * mediumRadius * medi << 168 164 } << 169 /// c) Partilces directed towards the bounding slice (default option!) >> 170 // Select a starting position on a sphere including the >> 171 // target volume and neuron morphology >> 172 else >> 173 { >> 174 //G4cerr << " Initial beam position uniformly spread on a sphere! "<< G4endl; >> 175 G4double cosTheta = 2.*G4UniformRand()-1; >> 176 G4double sinTheta = std::sqrt(1.-cosTheta*cosTheta); >> 177 G4double phi = twopi*G4UniformRand(); >> 178 G4ThreeVector positionStart(mediumRadius*sinTheta*std::cos(phi), >> 179 mediumRadius*sinTheta*std::sin(phi), >> 180 mediumRadius*cosTheta); >> 181 fpParticleGun->SetParticlePosition(positionStart); >> 182 // To compute the direction, select a point inside the target volume >> 183 G4ThreeVector positionDir( >> 184 boundingXHalfLength*(2.*G4UniformRand()-1), >> 185 boundingYHalfLength*(2.*G4UniformRand()-1), >> 186 boundingZHalfLength*(2.*G4UniformRand()-1)); >> 187 fpParticleGun->SetParticleMomentumDirection( >> 188 (positionDir-positionStart).unit()); >> 189 // Surface area of sphere >> 190 fGunArea = 4.*pi*mediumRadius*mediumRadius ; >> 191 // G4cerr << " Particle Fluence Area on sphere (um2) = " >> 192 // <<fGunArea / (um*um)<< G4endl; >> 193 } >> 194 165 } 195 } 166 else { << 196 else >> 197 { 167 G4cerr << "Bounding slice volume not found 198 G4cerr << "Bounding slice volume not found!" << G4endl; 168 G4cerr << "Default particle kinematic used 199 G4cerr << "Default particle kinematic used" << G4endl; 169 } << 200 } 170 201 171 fpParticleGun->GeneratePrimaryVertex(anEvent 202 fpParticleGun->GeneratePrimaryVertex(anEvent); 172 } 203 } 173 204