| 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 // Hadrontherapy advanced example for Geant4 << 26 // $Id: HadrontherapyDetectorSD.cc;
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // Last modified: G.A.P.Cirrone March 2008;
>> 28 //
>> 29 // See more at: http://geant4infn.wikispaces.com/HadrontherapyExample
>> 30 //
>> 31 // ----------------------------------------------------------------------------
>> 32 // GEANT 4 - Hadrontherapy example
>> 33 // ----------------------------------------------------------------------------
>> 34 // Code developed by:
>> 35 //
>> 36 // G.A.P. Cirrone(a)*, F. Di Rosa(a), S. Guatelli(b), G. Russo(a)
>> 37 //
>> 38 // (a) Laboratori Nazionali del Sud
>> 39 // of the National Institute for Nuclear Physics, Catania, Italy
>> 40 // (b) National Institute for Nuclear Physics Section of Genova, genova, Italy
>> 41 //
>> 42 // * cirrone@lns.infn.it
>> 43 // ----------------------------------------------------------------------------
28 44
29 #include "HadrontherapyDetectorSD.hh" 45 #include "HadrontherapyDetectorSD.hh"
>> 46 #include "HadrontherapyAnalysisManager.hh"
30 #include "HadrontherapyDetectorHit.hh" 47 #include "HadrontherapyDetectorHit.hh"
31 #include "G4Step.hh" 48 #include "G4Step.hh"
32 #include "G4VTouchable.hh" 49 #include "G4VTouchable.hh"
33 #include "G4TouchableHistory.hh" 50 #include "G4TouchableHistory.hh"
34 #include "G4SDManager.hh" 51 #include "G4SDManager.hh"
35 #include "HadrontherapyMatrix.hh" <<
36 #include "HadrontherapyLet.hh" <<
37 #include "G4Track.hh" <<
38 #include "G4SystemOfUnits.hh" <<
39 #include "HadrontherapyMatrix.hh" <<
40 52
41 << 53 HadrontherapyDetectorSD::HadrontherapyDetectorSD(G4String name):G4VSensitiveDetector(name)
42 #include "G4SteppingManager.hh" << 54 {
43 #include "G4TrackVector.hh" << 55 G4String HCname;
44 #include "HadrontherapySteppingAction.hh" << 56 collectionName.insert(HCname="HadrontherapyDetectorHitsCollection");
45 #include "G4ios.hh" << 57
46 #include "G4SteppingManager.hh" << 58 HitsCollection = NULL;
47 #include "G4Track.hh" << 59 G4String sensitiveDetectorName = name;
48 #include "G4Step.hh" <<
49 #include "G4StepPoint.hh" <<
50 #include "G4TrackStatus.hh" <<
51 #include "G4TrackVector.hh" <<
52 #include "G4ParticleDefinition.hh" <<
53 #include "G4ParticleTypes.hh" <<
54 #include "G4UserEventAction.hh" <<
55 #include "G4TransportationManager.hh" <<
56 #include "G4VSensitiveDetector.hh" <<
57 #include "HadrontherapyRunAction.hh" <<
58 #include "G4SystemOfUnits.hh" <<
59 #include "HadrontherapyRBE.hh" <<
60 #include <G4AccumulableManager.hh> <<
61 <<
62 <<
63 ////////////////////////////////////////////// <<
64 HadrontherapyDetectorSD::HadrontherapyDetector <<
65 G4VSensitiveDetector(name) <<
66 { <<
67 G4String HCname; <<
68 collectionName.insert(HCname="Hadrontherap <<
69 HitsCollection = NULL; <<
70 sensitiveDetectorName = name; <<
71 <<
72 } 60 }
73 61
74 ////////////////////////////////////////////// <<
75 HadrontherapyDetectorSD::~HadrontherapyDetecto 62 HadrontherapyDetectorSD::~HadrontherapyDetectorSD()
76 {} << 63 {
>> 64 }
77 65
78 ////////////////////////////////////////////// <<
79 void HadrontherapyDetectorSD::Initialize(G4HCo 66 void HadrontherapyDetectorSD::Initialize(G4HCofThisEvent*)
80 { << 67 {
81 << 68 HitsCollection = new HadrontherapyDetectorHitsCollection(sensitiveDetectorName,
82 HitsCollection = new HadrontherapyDetector << 69 collectionName[0]);
83 <<
84 } 70 }
85 71
86 ////////////////////////////////////////////// << 72
87 G4bool HadrontherapyDetectorSD::ProcessHits(G4 << 73 G4bool HadrontherapyDetectorSD::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist)
88 { 74 {
89 << 75 if(!ROhist)
90 << 76 return false;
91 if (aStep -> GetPreStepPoint() -> GetPhysi << 77
92 << 78 if(aStep -> GetPreStepPoint() -> GetPhysicalVolume() -> GetName() != "DetectorPhys")
93 << 79 return false;
94 // Get kinetic energy << 80
95 G4Track * theTrack = aStep -> GetTrack() << 81 G4double energyDeposit = aStep -> GetTotalEnergyDeposit();
96 G4double kineticEnergy = theTrack->GetKine << 82
97 G4ParticleDefinition *particleDef = theTra << 83 if(energyDeposit == 0.) return false;
98 //Get particle name << 84
99 G4String particleName = particleDef -> Ge << 85 // Read voxel indexes: i is the x index, k is the z index
100 << 86
101 // Get particle PDG code << 87 G4int k = ROhist -> GetReplicaNumber(0);
102 G4int pdg = particleDef ->GetPDGEncoding() << 88 G4int i = ROhist -> GetReplicaNumber(2);
103 << 89 G4int j = ROhist -> GetReplicaNumber(1);
104 // Get unique track_id (in an event) << 90
105 G4int trackID = theTrack -> GetTrackID(); << 91 G4String particleName = aStep -> GetTrack() -> GetDynamicParticle() ->
106 << 92 GetDefinition() -> GetParticleName();
107 G4double energyDeposit = aStep -> GetTotal << 93
108 << 94 if(energyDeposit != 0)
109 G4double DX = aStep -> GetStepLength(); << 95 {
110 G4int Z = particleDef-> GetAtomicNumber(); << 96 // Create a hit with the information of position is in the detector
111 G4int A = particleDef-> GetAtomicMass(); << 97 HadrontherapyDetectorHit* detectorHit = new HadrontherapyDetectorHit();
112 G4StepPoint* PreStep = aStep->GetPreStepPo << 98 detectorHit -> SetEdepAndPosition(i, j, k, energyDeposit);
113 << 99 HitsCollection -> insert(detectorHit);
114 // Read voxel indexes: i is the x index, k <<
115 const G4VTouchable* touchable = aStep->Get <<
116 G4int k = touchable->GetReplicaNumber(0); <<
117 G4int i = touchable->GetReplicaNumber(2); <<
118 G4int j = touchable->GetReplicaNumber(1); <<
119 <<
120 G4TouchableHandle touchPreStep = PreStep-> <<
121 G4VPhysicalVolume* volumePre = touchPreSte <<
122 G4String namePre = volumePre->GetName(); <<
123 <<
124 <<
125 <<
126 HadrontherapyMatrix* matrix = Hadrontherap <<
127 HadrontherapyLet* let = HadrontherapyLet:: <<
128 <<
129 G4int* hitTrack = matrix -> GetHitTrack(i, <<
130 <<
131 <<
132 // ******************** let ************* <<
133 if (let) <<
134 { <<
135 if ( !(Z==0 && A==1) ) // All but not <<
136 { <<
137 if( energyDeposit>0. && DX >0. )// <<
138 { <<
139 if (pdg !=22 && pdg !=11) // n <<
140 { <<
141 <<
142 // Get the pre-step kineti <<
143 G4double eKinPre = aStep - <<
144 // Get the post-step kinet <<
145 G4double eKinPost = aStep <<
146 // Get the step average ki <<
147 G4double eKinMean = (eKinP <<
148 <<
149 // get the material <<
150 const G4Material * materia <<
151 <<
152 // get the secondary patic <<
153 G4Step fstep = *theTrack - <<
154 // store all the secondary <<
155 const std::vector<const G4 <<
156 <<
157 size_t SecondarySize = (*s <<
158 G4double EnergySecondary = <<
159 <<
160 // get secondary electrons <<
161 if (SecondarySize) // calc <<
162 { <<
163 for (size_t numsec = 0 <<
164 { <<
165 //Get the PDG code <<
166 G4int PDGSecondary <<
167 <<
168 if(PDGSecondary == <<
169 { <<
170 // calculate t <<
171 EnergySecondar <<
172 } <<
173 } <<
174 <<
175 } <<
176 <<
177 // call the let filldatas <<
178 let -> FillEnergySpectrum( <<
179 <<
180 } <<
181 } <<
182 } <<
183 } <<
184 <<
185 <<
186 if (matrix) <<
187 { <<
188 <<
189 // Increment Fluences & accumulate ene <<
190 // Hit voxels are marked with track_id <<
191 //G4int* hitTrack = matrix -> GetHitTr <<
192 if ( *hitTrack != trackID ) <<
193 { <<
194 *hitTrack = trackID; <<
195 /* <<
196 * Fill FLUENCE data for every sin <<
197 * Exclude e-, neutrons, gamma, .. <<
198 */ <<
199 if ( Z >= 1) matrix -> Fill(trackI <<
200 <<
201 } <<
202 <<
203 if(energyDeposit != 0) <<
204 { <<
205 /* <<
206 * This method will fill a dose m <<
207 * A method of the HadrontherapyM <<
208 * is called automatically at the <<
209 * It permits to store all dose/f <<
210 */ <<
211 // if (A==1 && Z==1) // primary an <<
212 if ( Z>=1 ) // exclude e-, neu <<
213 matrix -> Fill(trackID, partic <<
214 /* <<
215 * Create a hit with the informati <<
216 */ <<
217 HadrontherapyDetectorHit* detector <<
218 detectorHit -> SetEdepAndPosition( <<
219 HitsCollection -> insert(detectorH <<
220 } <<
221 } 100 }
222 101
223 auto rbe = HadrontherapyRBE::GetInstance() << 102 // Energy deposit of secondary particles along X (integrated on Y and Z)
224 if (rbe->IsCalculationEnabled()) << 103 #ifdef G4ANALYSIS_USE
225 { <<
226 if (!fRBEAccumulable) <<
227 { <<
228 fRBEAccumulable = dynamic_cast<Had <<
229 if (!fRBEAccumulable) <<
230 { <<
231 G4Exception("HadrontherapyDete <<
232 } <<
233 } <<
234 if (A>0) //protect against gammas, e- , etc <<
235 { <<
236 fRBEAccumulable->Accumulate(kineticEnerg <<
237 } <<
238 } <<
239 104
>> 105 HadrontherapyAnalysisManager* analysis =
>> 106 HadrontherapyAnalysisManager::getInstance();
>> 107
>> 108 if(energyDeposit != 0)
>> 109 {
>> 110
>> 111 if(aStep -> GetTrack() -> GetTrackID()!= 1)
>> 112 {
>> 113 if (particleName == "proton")
>> 114 analysis -> SecondaryProtonEnergyDeposit(i, energyDeposit/MeV);
>> 115
>> 116 if (particleName == "neutron")
>> 117 analysis -> SecondaryNeutronEnergyDeposit(i, energyDeposit/MeV);
>> 118
>> 119 if (particleName == "alpha")
>> 120 analysis -> SecondaryAlphaEnergyDeposit(i, energyDeposit/MeV);
>> 121
>> 122 if (particleName == "gamma")
>> 123 analysis -> SecondaryGammaEnergyDeposit(i, energyDeposit/MeV);
>> 124
>> 125 if (particleName == "e-")
>> 126 analysis -> SecondaryElectronEnergyDeposit(i, energyDeposit/MeV);
>> 127
>> 128 if (particleName == "triton")
>> 129 analysis -> SecondaryTritonEnergyDeposit(i, energyDeposit/MeV);
>> 130
>> 131 if (particleName == "deuteron")
>> 132 analysis -> SecondaryDeuteronEnergyDeposit(i, energyDeposit/MeV);
>> 133
>> 134 if (particleName == "pi+" || particleName == "pi-" || particleName == "pi0")
>> 135 analysis -> SecondaryPionEnergyDeposit(i, energyDeposit/MeV);
>> 136 }
>> 137 }
>> 138 #endif
240 139
241 return true; << 140 return true;
242 } 141 }
243 142
244 ////////////////////////////////////////////// <<
245 void HadrontherapyDetectorSD::EndOfEvent(G4HCo 143 void HadrontherapyDetectorSD::EndOfEvent(G4HCofThisEvent* HCE)
246 { 144 {
247 << 145 static G4int HCID = -1;
248 static G4int HCID = -1; << 146 if(HCID < 0)
249 if(HCID < 0) << 147 {
250 { << 148 HCID = GetCollectionID(0);
251 HCID = GetCollectionID(0); <<
252 } 149 }
253 << 150 HCE -> AddHitsCollection(HCID,HitsCollection);
254 HCE -> AddHitsCollection(HCID,HitsCollecti <<
255 } 151 }
256 152
257 153