| 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