Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // This example is provided by the Geant4-DNA collaboration 27 // Any report or published results obtained using the Geant4-DNA software 28 // shall cite the following Geant4-DNA collaboration publication: 29 // Med. Phys. 37 (2010) 4692-4708 30 // and papers 31 // M. Batmunkh et al. J Radiat Res Appl Sci 8 (2015) 498-507 32 // O. Belov et al. Physica Medica 32 (2016) 1510-1520 33 // The Geant4-DNA web site is available at http://geant4-dna.org 34 // 35 // ------------------------------------------------------------------- 36 // November 2016 37 // ------------------------------------------------------------------- 38 // 39 /// \file PrimaryGeneratorAction.cc 40 /// \brief Implementation of the PrimaryGeneratorAction class 41 42 #include "PrimaryGeneratorAction.hh" 43 44 #include "G4SystemOfUnits.hh" 45 // 46 #include "CommandLineParser.hh" 47 48 #include "G4Box.hh" 49 #include "G4Event.hh" 50 #include "G4LogicalVolume.hh" 51 #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" 59 #include "Randomize.hh" 60 61 using namespace G4DNAPARSER; 62 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 64 65 PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction() 66 { 67 G4int n_particle = 1; 68 fpParticleGun = new G4ParticleGun(n_particle); 69 // 70 G4ParticleDefinition* particle = G4Proton::Proton(); 71 fpParticleGun->SetParticleDefinition(particle); 72 // default gun parameters 73 fpParticleGun->SetParticleEnergy(10. * MeV); 74 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.)); 75 fpParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 0.)); 76 } 77 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 79 80 PrimaryGeneratorAction::~PrimaryGeneratorAction() 81 { 82 delete fpParticleGun; 83 } 84 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 86 87 void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) 88 { 89 // Initial kinetic energy of particles beam as Gaussion distribution! 90 91 // G4double Ekin = 10.*MeV; 92 // G4double deltaEkin = 9.0955e-5*MeV; 93 // fpParticleGun->SetParticleEnergy(G4RandGauss::shoot(Ekin,deltaEkin)); 94 95 // In order to avoid dependence of PrimaryGeneratorAction 96 // on DetectorConstruction class we get world volume 97 // from G4LogicalVolumeStore 98 99 // We included three options for particles direction: 100 101 G4double mediumRadius = 0.; 102 G4double boundingXHalfLength = 0.; 103 G4double boundingYHalfLength = 0.; 104 G4double boundingZHalfLength = 0.; 105 G4LogicalVolume* mediumLV = G4LogicalVolumeStore::GetInstance()->GetVolume("Medium"); 106 G4LogicalVolume* boundingLV = G4LogicalVolumeStore::GetInstance()->GetVolume("BoundingSlice"); 107 G4Orb* mediumSphere = 0; 108 G4Box* boundingSlice = 0; 109 if (mediumLV && boundingLV) { 110 mediumSphere = dynamic_cast<G4Orb*>(mediumLV->GetSolid()); 111 boundingSlice = dynamic_cast<G4Box*>(boundingLV->GetSolid()); 112 } 113 if (mediumSphere && boundingSlice) { 114 mediumRadius = mediumSphere->GetRadius(); 115 boundingXHalfLength = boundingSlice->GetXHalfLength(); 116 boundingYHalfLength = boundingSlice->GetYHalfLength(); 117 boundingZHalfLength = boundingSlice->GetZHalfLength(); 118 119 /// a) Partilces directed to "square" on the XY plane of bounding slice 120 /// (or YZ, XZ) 121 122 if (CommandLineParser::GetParser()->GetCommandIfActive("-sXY")) { 123 // INITIAL BEAM DIVERGENCE 124 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.)); 125 // // INITIAL BEAM POSITION 126 fpParticleGun->SetParticlePosition(G4ThreeVector( 127 CLHEP::RandFlat::shoot(-boundingXHalfLength, boundingXHalfLength), 128 CLHEP::RandFlat::shoot(-boundingYHalfLength, boundingYHalfLength), -mediumRadius)); 129 } 130 131 /// b) Partilces directed to "disk" on the XY plane of 132 /// bounding slice (or YZ, XZ) 133 134 else if (CommandLineParser::GetParser()->GetCommandIfActive("-dXY")) { 135 fpParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.)); 136 G4double x0, y0, z0; 137 x0 = 100. * mm; 138 y0 = 100. * mm; 139 z0 = -mediumRadius; // mediumRadius; 140 while (!(std::sqrt(x0 * x0 + y0 * y0) <= mediumRadius)) { 141 x0 = CLHEP::RandFlat::shoot(-mediumRadius, mediumRadius); 142 y0 = CLHEP::RandFlat::shoot(-mediumRadius, mediumRadius); 143 } 144 fpParticleGun->SetParticlePosition(G4ThreeVector(x0, y0, z0)); 145 } 146 147 /// c) Partilces directed towards the bounding slice (default option!) 148 // Select a starting position on a sphere including the 149 // target volume and neuron morphology 150 else { 151 G4double cosTheta = 2. * G4UniformRand() - 1; 152 G4double sinTheta = std::sqrt(1. - cosTheta * cosTheta); 153 G4double phi = twopi * G4UniformRand(); 154 G4ThreeVector positionStart(mediumRadius * sinTheta * std::cos(phi), 155 mediumRadius * sinTheta * std::sin(phi), mediumRadius * cosTheta); 156 fpParticleGun->SetParticlePosition(positionStart); 157 // To compute the direction, select a point inside the target volume 158 G4ThreeVector positionDir(boundingXHalfLength * (2. * G4UniformRand() - 1), 159 boundingYHalfLength * (2. * G4UniformRand() - 1), 160 boundingZHalfLength * (2. * G4UniformRand() - 1)); 161 fpParticleGun->SetParticleMomentumDirection((positionDir - positionStart).unit()); 162 // Surface area of sphere 163 fGunArea = 4. * pi * mediumRadius * mediumRadius; 164 } 165 } 166 else { 167 G4cerr << "Bounding slice volume not found!" << G4endl; 168 G4cerr << "Default particle kinematic used" << G4endl; 169 } 170 171 fpParticleGun->GeneratePrimaryVertex(anEvent); 172 } 173