| 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 // This is the *BASIC* version of Hadrontherapy, a Geant4-based application
27 // See more at: https://twiki.cern.ch/twiki/bi << 27 // See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
>> 28 //
>> 29 // Visit the Hadrontherapy web site (http://www.lns.infn.it/link/Hadrontherapy) to request
>> 30 // the *COMPLETE* version of this program, together with its documentation;
>> 31 // Hadrontherapy (both basic and full version) are supported by the Italian INFN
>> 32 // Institute in the framework of the MC-INFN Group
>> 33 //
>> 34
28 35
29 #include "G4UnitsTable.hh" 36 #include "G4UnitsTable.hh"
30 #include "G4SDManager.hh" 37 #include "G4SDManager.hh"
31 #include "G4RunManager.hh" 38 #include "G4RunManager.hh"
32 #include "G4GeometryManager.hh" 39 #include "G4GeometryManager.hh"
33 #include "G4SolidStore.hh" 40 #include "G4SolidStore.hh"
34 #include "G4PhysicalVolumeStore.hh" 41 #include "G4PhysicalVolumeStore.hh"
35 #include "G4LogicalVolumeStore.hh" 42 #include "G4LogicalVolumeStore.hh"
36 #include "G4Box.hh" 43 #include "G4Box.hh"
37 #include "G4LogicalVolume.hh" 44 #include "G4LogicalVolume.hh"
38 #include "G4ThreeVector.hh" 45 #include "G4ThreeVector.hh"
39 #include "G4PVPlacement.hh" 46 #include "G4PVPlacement.hh"
40 #include "globals.hh" 47 #include "globals.hh"
41 #include "G4Transform3D.hh" 48 #include "G4Transform3D.hh"
42 #include "G4RotationMatrix.hh" 49 #include "G4RotationMatrix.hh"
43 #include "G4Colour.hh" 50 #include "G4Colour.hh"
44 #include "G4UserLimits.hh" 51 #include "G4UserLimits.hh"
45 #include "G4UnitsTable.hh" 52 #include "G4UnitsTable.hh"
46 #include "G4VisAttributes.hh" 53 #include "G4VisAttributes.hh"
47 #include "G4NistManager.hh" 54 #include "G4NistManager.hh"
>> 55
48 #include "HadrontherapyDetectorConstruction.hh 56 #include "HadrontherapyDetectorConstruction.hh"
49 #include "HadrontherapyDetectorROGeometry.hh" 57 #include "HadrontherapyDetectorROGeometry.hh"
50 #include "HadrontherapyDetectorMessenger.hh" 58 #include "HadrontherapyDetectorMessenger.hh"
51 #include "HadrontherapyDetectorSD.hh" 59 #include "HadrontherapyDetectorSD.hh"
52 #include "HadrontherapyMatrix.hh" 60 #include "HadrontherapyMatrix.hh"
53 #include "HadrontherapyLet.hh" 61 #include "HadrontherapyLet.hh"
54 #include "PassiveProtonBeamLine.hh" 62 #include "PassiveProtonBeamLine.hh"
55 #include "BESTPassiveProtonBeamLine.hh" << 63 #include "HadrontherapyAnalysisManager.hh"
56 #include "HadrontherapyMatrix.hh" <<
57 <<
58 #include "HadrontherapyRBE.hh" <<
59 #include "G4SystemOfUnits.hh" 64 #include "G4SystemOfUnits.hh"
60 65
61 #include <cmath> 66 #include <cmath>
62 67
63 <<
64 <<
65 HadrontherapyDetectorConstruction* Hadronthera 68 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::instance = 0;
66 ////////////////////////////////////////////// 69 /////////////////////////////////////////////////////////////////////////////
67 HadrontherapyDetectorConstruction::Hadronthera 70 HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
68 : motherPhys(physicalTreatmentRoom), // pointe << 71 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
69 detectorSD(0), detectorROGeometry(0), matrix(0 << 72 detectorSD(0), detectorROGeometry(0), matrix(0),
70 phantom(0), detector(0), << 73 phantom(0), detector(0),
71 phantomLogicalVolume(0), detectorLogicalVolume << 74 phantomLogicalVolume(0), detectorLogicalVolume(0),
72 phantomPhysicalVolume(0), detectorPhysicalVolu << 75 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
73 aRegion(0) << 76 aRegion(0)
74 { << 77 {
75 << 78 HadrontherapyAnalysisManager::GetInstance();
76 << 79
77 /* NOTE! that the HadrontherapyDetectorCon << 80 /* NOTE! that the HadrontherapyDetectorConstruction class
78 * does NOT inherit from G4VUserDetectorCo << 81 * does NOT inherit from G4VUserDetectorConstruction G4 class
79 * So the Construct() mandatory virtual me << 82 * So the Construct() mandatory virtual method is inside another geometric class
80 * like the passiveProtonBeamLIne, ... << 83 * like the passiveProtonBeamLIne, ...
81 */ << 84 */
82 << 85
83 // Messenger to change parameters of the p << 86 // Messenger to change parameters of the phantom/detector geometry
84 detectorMessenger = new HadrontherapyDetec << 87 detectorMessenger = new HadrontherapyDetectorMessenger(this);
85 << 88
86 // Default detector voxels size << 89 // Default detector voxels size
87 // 200 slabs along the beam direction (X) << 90 // 200 slabs along the beam direction (X)
88 sizeOfVoxelAlongX = 200 *um; << 91 sizeOfVoxelAlongX = 200 *um;
89 sizeOfVoxelAlongY = 4 *cm; << 92 sizeOfVoxelAlongY = 4 *cm;
90 sizeOfVoxelAlongZ = 4 *cm; << 93 sizeOfVoxelAlongZ = 4 *cm;
91 << 94
92 // Define here the material of the water p << 95 // Define here the material of the water phantom and of the detector
93 SetPhantomMaterial("G4_WATER"); << 96 SetPhantomMaterial("G4_WATER");
94 // Construct geometry (messenger commands) << 97 // Construct geometry (messenger commands)
95 // SetDetectorSize(4.*cm, 4.*cm, 4.*cm); << 98 SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
96 SetDetectorSize(4. *cm, 4. *cm, 4. *cm); << 99 SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
97 SetPhantomSize(40. *cm, 40. *cm, 40. *cm); << 100 SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
98 << 101 SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 18. *cm, 18. *cm));
99 SetPhantomPosition(G4ThreeVector(20. *cm, << 102 SetDetectorPosition();
100 SetDetectorToPhantomPosition(G4ThreeVector << 103 //GetDetectorToWorldPosition();
101 SetDetectorPosition(); << 104
102 //GetDetectorToWorldPosition(); << 105
103 << 106 // Write virtual parameters to the real ones and check for consistency
104 // Write virtual parameters to the real on << 107 UpdateGeometry();
105 UpdateGeometry(); <<
106 <<
107 <<
108 <<
109 } 108 }
110 109
111 ////////////////////////////////////////////// 110 /////////////////////////////////////////////////////////////////////////////
112 HadrontherapyDetectorConstruction::~Hadronther 111 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
113 { << 112 {
114 delete detectorROGeometry; << 113 delete detectorROGeometry;
115 delete matrix; << 114 delete matrix;
116 delete detectorMessenger; 115 delete detectorMessenger;
117 } 116 }
118 117
119 ////////////////////////////////////////////// << 118
120 HadrontherapyDetectorConstruction* Hadronthera 119 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::GetInstance()
121 { 120 {
122 return instance; << 121
>> 122 return instance;
123 } 123 }
124 124
>> 125
125 ////////////////////////////////////////////// 126 /////////////////////////////////////////////////////////////////////////////
126 // ConstructPhantom() is the method that const << 127 // ConstructPhantom() is the method that reconstuct a water box (called phantom
127 // (or water phantom) in the usual Medical phy << 128 // (or water phantom) in the usual Medical physicists slang).
128 // A water phantom can be considered a good ap << 129 // A water phantom can be considered a good
>> 130 // approximation of a an human body.
129 void HadrontherapyDetectorConstruction::Constr 131 void HadrontherapyDetectorConstruction::ConstructPhantom()
130 { 132 {
131 // Definition of the solid volume of the P 133 // Definition of the solid volume of the Phantom
132 phantom = new G4Box("Phantom", << 134 phantom = new G4Box("Phantom",
133 phantomSizeX/2, << 135 phantomSizeX/2,
134 phantomSizeY/2, << 136 phantomSizeY/2,
135 phantomSizeZ/2); << 137 phantomSizeZ/2);
136 << 138
137 // Definition of the logical volume of the << 139 // Definition of the logical volume of the Phantom
138 phantomLogicalVolume = new G4LogicalVolume << 140 phantomLogicalVolume = new G4LogicalVolume(phantom,
139 << 141 phantomMaterial,
140 << 142 "phantomLog", 0, 0, 0);
141 << 143
142 // Definition of the physics volume of the 144 // Definition of the physics volume of the Phantom
143 phantomPhysicalVolume = new G4PVPlacement( 145 phantomPhysicalVolume = new G4PVPlacement(0,
144 << 146 phantomPosition,
145 << 147 "phantomPhys",
146 << 148 phantomLogicalVolume,
147 << 149 motherPhys,
148 << 150 false,
149 << 151 0);
150 << 152
151 // Visualisation attributes of the phantom << 153 // Visualisation attributes of the phantom
152 red = new G4VisAttributes(G4Colour(255/255 154 red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
153 red -> SetVisibility(true); 155 red -> SetVisibility(true);
154 red -> SetForceSolid(true); 156 red -> SetForceSolid(true);
155 red -> SetForceWireframe(true); 157 red -> SetForceWireframe(true);
156 phantomLogicalVolume -> SetVisAttributes(r 158 phantomLogicalVolume -> SetVisAttributes(red);
157 } 159 }
158 160
159 ////////////////////////////////////////////// 161 /////////////////////////////////////////////////////////////////////////////
160 // ConstructDetector() is the method the recon << 162 // ConstructDetector() it the method the reconstruct a detector region
161 // inside the water phantom. It is a volume, l << 163 // inside the water phantom. It is a volume, located inside the water phantom
>> 164 // and with two coincident faces:
162 // 165 //
163 // ************************** 166 // **************************
164 // * water phantom * << 167 // * water phantom *
165 // * * 168 // * *
166 // * * 169 // * *
167 // *--------------- * 170 // *--------------- *
168 // Beam * - * 171 // Beam * - *
169 // -----> * detector - * 172 // -----> * detector - *
170 // * - * 173 // * - *
171 // *--------------- * 174 // *--------------- *
172 // * * 175 // * *
173 // * * 176 // * *
174 // * * 177 // * *
175 // ************************** 178 // **************************
176 // 179 //
177 // The detector can be dived in slices or voxe << 180 // The detector is the volume that can be dived in slices or voxelized
178 // and inside it different quantities (dose di << 181 // and in it we can collect a number of usefull information:
179 // can be stored. << 182 // dose distribution, fluence distribution, LET and so on
180 void HadrontherapyDetectorConstruction::Constr 183 void HadrontherapyDetectorConstruction::ConstructDetector()
181 <<
182 { 184 {
>> 185
183 // Definition of the solid volume of the D 186 // Definition of the solid volume of the Detector
184 detector = new G4Box("Detector", << 187 detector = new G4Box("Detector",
185 << 188 detectorSizeX/2,
186 phantomSizeX/2, << 189 detectorSizeY/2,
187 << 190 detectorSizeZ/2);
188 phantomSizeY/2, <<
189 <<
190 phantomSizeZ/2); <<
191 191
192 // Definition of the logic volume of the P 192 // Definition of the logic volume of the Phantom
193 detectorLogicalVolume = new G4LogicalVolum 193 detectorLogicalVolume = new G4LogicalVolume(detector,
194 << 194 detectorMaterial,
195 << 195 "DetectorLog",
196 << 196 0,0,0);
197 // Definition of the physical volume of th << 197 // Definition of the physical volume of the Phantom
198 detectorPhysicalVolume = new G4PVPlacement << 198 detectorPhysicalVolume = new G4PVPlacement(0,
199 << 199 detectorPosition, // Setted by displacement
200 << 200 "DetectorPhys",
201 << 201 detectorLogicalVolume,
202 << 202 phantomPhysicalVolume,
203 << 203 false,0);
204 << 204
205 // Visualisation attributes of the detecto << 205 // Visualisation attributes of the detector
206 skyBlue = new G4VisAttributes( G4Colour(13 206 skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. , 235/255. ));
207 skyBlue -> SetVisibility(true); 207 skyBlue -> SetVisibility(true);
208 skyBlue -> SetForceSolid(true); 208 skyBlue -> SetForceSolid(true);
209 //skyBlue -> SetForceWireframe(true); 209 //skyBlue -> SetForceWireframe(true);
210 detectorLogicalVolume -> SetVisAttributes( 210 detectorLogicalVolume -> SetVisAttributes(skyBlue);
211 << 211
212 // ************** << 212 // **************
213 // ************** << 213 // Cut per Region
214 // Cut per Region << 214 // **************
215 // ************** << 215
216 // << 216 // A smaller cut is fixed in the phantom to calculate the energy deposit with the
217 // A smaller cut is fixed in the phantom t << 217 // required accuracy
218 // required accuracy <<
219 if (!aRegion) 218 if (!aRegion)
220 { 219 {
221 aRegion = new G4Region("DetectorLog"); << 220 aRegion = new G4Region("DetectorLog");
222 detectorLogicalVolume -> SetRegion(aRe << 221 detectorLogicalVolume -> SetRegion(aRegion);
223 aRegion->AddRootLogicalVolume( detecto << 222 aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
224 } 223 }
>> 224
225 } 225 }
226 226
>> 227
227 ////////////////////////////////////////////// 228 ///////////////////////////////////////////////////////////////////////
>> 229
228 void HadrontherapyDetectorConstruction::Initia 230 void HadrontherapyDetectorConstruction::InitializeDetectorROGeometry(
229 << 231 HadrontherapyDetectorROGeometry* RO,
230 << 232 G4ThreeVector
231 << 233 detectorToWorldPosition)
232 { << 234 {
233 RO->Initialize(detectorToWorldPosition, << 235 RO->Initialize(detectorToWorldPosition,
234 detectorSizeX/2, << 236 detectorSizeX/2,
235 detectorSizeY/2, << 237 detectorSizeY/2,
236 detectorSizeZ/2, << 238 detectorSizeZ/2,
237 numberOfVoxelsAlongX, << 239 numberOfVoxelsAlongX,
238 numberOfVoxelsAlongY, << 240 numberOfVoxelsAlongY,
239 numberOfVoxelsAlongZ); << 241 numberOfVoxelsAlongZ);
240 } <<
241 void HadrontherapyDetectorConstruction::Virtua <<
242 { <<
243 <<
244 //Virtual plane <<
245 VirtualLayerPosition = G4ThreeVector(0*cm, <<
246 NewSource= Varbool; <<
247 if(NewSource == true) <<
248 { <<
249 // std::cout<<"trr"<<std::endl; <<
250 G4Material* airNist = G4NistManager:: <<
251 <<
252 solidVirtualLayer = new G4Box("Virtual <<
253 1.*um, <<
254 20.*cm, <<
255 40.*cm); <<
256 <<
257 logicVirtualLayer = new G4LogicalVolum <<
258 <<
259 <<
260 <<
261 <<
262 physVirtualLayer= new G4PVPlacement(0, <<
263 "V <<
264 lo <<
265 mo <<
266 fa <<
267 0) <<
268 <<
269 logicVirtualLayer -> SetVisAttributes( <<
270 } <<
271 <<
272 <<
273 <<
274 <<
275 } 242 }
276 243
277 244
278 ////////////////////////////////////////////// << 245
>> 246
>> 247
>> 248 /////////////////////////////////////////////////////////
>> 249
279 void HadrontherapyDetectorConstruction::Param 250 void HadrontherapyDetectorConstruction::ParametersCheck()
280 { 251 {
281 // Check phantom/detector sizes & relative 252 // Check phantom/detector sizes & relative position
282 if (!IsInside(detectorSizeX, << 253 if (!IsInside(detectorSizeX,
283 detectorSizeY, << 254 detectorSizeY,
284 detectorSizeZ, << 255 detectorSizeZ,
285 phantomSizeX, << 256 phantomSizeX,
286 phantomSizeY, << 257 phantomSizeY,
287 phantomSizeZ, << 258 phantomSizeZ,
288 detectorToPhantomPosition << 259 detectorToPhantomPosition
289 )) << 260 ))
290 G4Exception("HadrontherapyDetectorCons << 261 G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
291 << 262
292 // Check Detector sizes respect to the vox 263 // Check Detector sizes respect to the voxel ones
293 << 264
294 if ( detectorSizeX < sizeOfVoxelAlongX) { 265 if ( detectorSizeX < sizeOfVoxelAlongX) {
295 G4Exception("HadrontherapyDetectorCons << 266 G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
296 } 267 }
297 if ( detectorSizeY < sizeOfVoxelAlongY) { 268 if ( detectorSizeY < sizeOfVoxelAlongY) {
298 G4Exception(" HadrontherapyDetectorCon << 269 G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
299 } 270 }
300 if ( detectorSizeZ < sizeOfVoxelAlongZ) { 271 if ( detectorSizeZ < sizeOfVoxelAlongZ) {
301 G4Exception(" HadrontherapyDetectorCon << 272 G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
302 } 273 }
>> 274
303 } 275 }
>> 276 /////////////////
>> 277 // MESSENGERS //
>> 278 ////////////////
304 279
305 ////////////////////////////////////////////// <<
306 G4bool HadrontherapyDetectorConstruction::SetP 280 G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
307 { 281 {
308 << 282
309 if (G4Material* pMat = G4NistManager::Inst 283 if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
310 { 284 {
311 phantomMaterial = pMat; << 285 phantomMaterial = pMat;
312 detectorMaterial = pMat; << 286 detectorMaterial = pMat;
313 if (detectorLogicalVolume && phantomLo << 287 if (detectorLogicalVolume && phantomLogicalVolume)
314 { << 288 {
315 detectorLogicalVolume -> SetMateri << 289 detectorLogicalVolume -> SetMaterial(pMat);
316 phantomLogicalVolume -> SetMateri << 290 phantomLogicalVolume -> SetMaterial(pMat);
317 << 291
318 G4RunManager::GetRunManager() -> P << 292 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
319 G4RunManager::GetRunManager() -> G << 293 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
320 G4cout << "The material of Phantom << 294 G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
321 } << 295 }
322 } 296 }
323 else 297 else
324 { 298 {
325 G4cout << "WARNING: material \"" << ma << 299 G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
326 " table [located in $G4INSTALL/source/ << 300 " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
327 G4cout << "Use command \"/parameter/ni << 301 G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
328 return false; << 302 return false;
329 } 303 }
330 << 304
331 return true; 305 return true;
332 } 306 }
333 ////////////////////////////////////////////// 307 /////////////////////////////////////////////////////////////////////////////
334 void HadrontherapyDetectorConstruction::SetPha 308 void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
335 { 309 {
336 if (sizeX > 0.) phantomSizeX = sizeX; 310 if (sizeX > 0.) phantomSizeX = sizeX;
337 if (sizeY > 0.) phantomSizeY = sizeY; 311 if (sizeY > 0.) phantomSizeY = sizeY;
338 if (sizeZ > 0.) phantomSizeZ = sizeZ; 312 if (sizeZ > 0.) phantomSizeZ = sizeZ;
339 } 313 }
340 << 314 /////////////////////////////////////////////////////////////////////////////
341 ////////////////////////////////////////////// 315 /////////////////////////////////////////////////////////////////////////////
342 void HadrontherapyDetectorConstruction::SetDet 316 void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
343 { 317 {
344 if (sizeX > 0.) {detectorSizeX = sizeX;} 318 if (sizeX > 0.) {detectorSizeX = sizeX;}
345 if (sizeY > 0.) {detectorSizeY = sizeY;} 319 if (sizeY > 0.) {detectorSizeY = sizeY;}
346 if (sizeZ > 0.) {detectorSizeZ = sizeZ;} 320 if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
347 SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxe 321 SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
348 } 322 }
349 <<
350 ////////////////////////////////////////////// 323 /////////////////////////////////////////////////////////////////////////////
>> 324
351 void HadrontherapyDetectorConstruction::SetVox 325 void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
352 { 326 {
353 if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX 327 if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
354 if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY 328 if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
355 if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ 329 if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
356 } 330 }
357 <<
358 ////////////////////////////////////////////// <<
359 void HadrontherapyDetectorConstruction::SetPha 331 void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
360 { 332 {
361 phantomPosition = pos; 333 phantomPosition = pos;
362 } 334 }
363 335
364 ////////////////////////////////////////////// 336 /////////////////////////////////////////////////////////////////////////////
365 void HadrontherapyDetectorConstruction::SetDet 337 void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
366 { 338 {
367 detectorToPhantomPosition = displ; 339 detectorToPhantomPosition = displ;
368 } 340 }
369 <<
370 void HadrontherapyDetectorConstruction::SetVir <<
371 { <<
372 <<
373 VirtualLayerPosition = position; <<
374 physVirtualLayer->SetTranslation(VirtualLa <<
375 <<
376 } <<
377 ////////////////////////////////////////////// 341 /////////////////////////////////////////////////////////////////////////////
378 void HadrontherapyDetectorConstruction::Update 342 void HadrontherapyDetectorConstruction::UpdateGeometry()
379 { 343 {
380 /* << 344 /*
381 * Check parameters consistency << 345 * Check parameters consistency
382 */ << 346 */
383 ParametersCheck(); << 347 ParametersCheck();
384 << 348
385 G4GeometryManager::GetInstance() -> OpenGe << 349 G4GeometryManager::GetInstance() -> OpenGeometry();
386 if (phantom) << 350 if (phantom)
387 { 351 {
388 phantom -> SetXHalfLength(phantomSizeX << 352 phantom -> SetXHalfLength(phantomSizeX/2);
389 phantom -> SetYHalfLength(phantomSizeY << 353 phantom -> SetYHalfLength(phantomSizeY/2);
390 phantom -> SetZHalfLength(phantomSizeZ << 354 phantom -> SetZHalfLength(phantomSizeZ/2);
391 << 355 phantomPhysicalVolume -> SetTranslation(phantomPosition);
392 phantomPhysicalVolume -> SetTranslatio <<
393 } 356 }
394 else ConstructPhantom(); << 357 else ConstructPhantom();
395 << 358
396 << 359 // Get the center of the detector
397 // Get the center of the detector << 360 SetDetectorPosition();
398 SetDetectorPosition(); << 361 if (detector)
399 if (detector) <<
400 { <<
401 <<
402 detector -> SetXHalfLength(detectorSiz <<
403 detector -> SetYHalfLength(detectorSiz <<
404 detector -> SetZHalfLength(detectorSiz <<
405 <<
406 detectorPhysicalVolume -> SetTranslati <<
407 } <<
408 else ConstructDetector(); <<
409 <<
410 //std::cout<<NewSource<<std::endl; <<
411 /*if(NewSource) <<
412 { <<
413 std::cout<<"via"<<std::endl; <<
414 }*/ <<
415 <<
416 <<
417 // std::cout<<"i"<<std::endl; <<
418 // std::cout<<VirtualLayerPosition<<std::e <<
419 // physVirtualLayer->SetTranslation(Virtua <<
420 <<
421 <<
422 <<
423 <<
424 <<
425 // Round to nearest integer number of voxe <<
426 <<
427 numberOfVoxelsAlongX = G4lrint(detectorSiz <<
428 sizeOfVoxelAlongX = ( detectorSizeX / numb <<
429 numberOfVoxelsAlongY = G4lrint(detectorSiz <<
430 sizeOfVoxelAlongY = ( detectorSizeY / numb <<
431 numberOfVoxelsAlongZ = G4lrint(detectorSiz <<
432 sizeOfVoxelAlongZ = ( detectorSizeZ / numb <<
433 PassiveProtonBeamLine *ppbl= (PassiveProto <<
434 <<
435 G4RunManager::GetRunManager()->GetUserDete <<
436 <<
437 HadrontherapyDetectorROGeometry* RO = (Had <<
438 <<
439 //Set parameters, either for the Construct <<
440 RO->Initialize(GetDetectorToWorldPosition( <<
441 detectorSizeX/2, <<
442 detectorSizeY/2, <<
443 detectorSizeZ/2, <<
444 numberOfVoxelsAlongX, <<
445 numberOfVoxelsAlongY, <<
446 numberOfVoxelsAlongZ); <<
447 <<
448 //This method below has an effect only if <<
449 RO->UpdateROGeometry(); <<
450 <<
451 <<
452 <<
453 volumeOfVoxel = sizeOfVoxelAlongX * sizeOf <<
454 massOfVoxel = detectorMaterial -> GetDensi <<
455 // This will clear the existing matrix (t <<
456 matrix = HadrontherapyMatrix::GetInstance( <<
457 <<
458 <<
459 <<
460 <<
461 <<
462 // Initialize RBE <<
463 HadrontherapyRBE::CreateInstance(numberOfV <<
464 <<
465 // Comment out the line below if let calcu <<
466 // Initialize LET with energy of primaries <<
467 if ( (let = HadrontherapyLet::GetInstance( <<
468 { 362 {
469 HadrontherapyLet::GetInstance() -> Ini << 363 detector -> SetXHalfLength(detectorSizeX/2);
>> 364 detector -> SetYHalfLength(detectorSizeY/2);
>> 365 detector -> SetZHalfLength(detectorSizeZ/2);
>> 366 detectorPhysicalVolume -> SetTranslation(detectorPosition);
470 } 367 }
471 << 368 else ConstructDetector();
472 << 369
473 // Initialize analysis << 370 // Round to nearest integer number of voxel
474 // Inform the kernel about the new geometr << 371
475 G4RunManager::GetRunManager() -> GeometryH << 372 numberOfVoxelsAlongX = G4lrint(detectorSizeX / sizeOfVoxelAlongX);
476 G4RunManager::GetRunManager() -> PhysicsHa << 373
477 << 374 sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
478 PrintParameters(); << 375
479 << 376 numberOfVoxelsAlongY = G4lrint(detectorSizeY / sizeOfVoxelAlongY);
480 // CheckOverlaps(); << 377
481 } << 378 sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
482 << 379
483 ////////////////////////////////////////////// << 380 numberOfVoxelsAlongZ = G4lrint(detectorSizeZ / sizeOfVoxelAlongZ);
484 //Check of the geometry << 381
485 ////////////////////////////////////////////// << 382 sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
486 void HadrontherapyDetectorConstruction::CheckO << 383
487 { << 384 PassiveProtonBeamLine *ppbl= (PassiveProtonBeamLine*)
488 G4PhysicalVolumeStore* thePVStore = G4Phys << 385
489 G4cout << thePVStore->size() << " physical << 386 G4RunManager::GetRunManager()->GetUserDetectorConstruction();
490 G4bool overlapFlag = false; << 387
491 G4int res=1000; << 388
492 G4double tol=0.; //tolerance << 389 HadrontherapyDetectorROGeometry* RO = (HadrontherapyDetectorROGeometry*) ppbl->GetParallelWorld(0);
493 for (size_t i=0;i<thePVStore->size();i++) << 390
>> 391 //Set parameters, either for the Construct() or for the UpdateROGeometry()
>> 392 RO->Initialize(GetDetectorToWorldPosition(),
>> 393 detectorSizeX/2,
>> 394 detectorSizeY/2,
>> 395 detectorSizeZ/2,
>> 396 numberOfVoxelsAlongX,
>> 397 numberOfVoxelsAlongY,
>> 398 numberOfVoxelsAlongZ);
>> 399
>> 400 //This method below has an effect only if the RO geometry is already built.
>> 401 RO->UpdateROGeometry();
>> 402
>> 403 volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
>> 404 massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
>> 405 // This will clear the existing matrix (together with all data inside it)!
>> 406 matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
>> 407 numberOfVoxelsAlongY,
>> 408 numberOfVoxelsAlongZ,
>> 409 massOfVoxel);
>> 410
>> 411
>> 412 // Comment out the line below if let calculation is not needed!
>> 413 // Initialize LET with energy of primaries and clear data inside
>> 414 if ( (let = HadrontherapyLet::GetInstance(this)) )
494 { 415 {
495 //overlapFlag = (*thePVStore)[i]->Chec << 416 HadrontherapyLet::GetInstance() -> Initialize();
496 overlapFlag = (*thePVStore)[i]->CheckO << 417 }
497 if (overlapFlag) << 418
498 G4cout << "Check: there are overlappin << 419
>> 420 // Initialize analysis
>> 421 #ifdef G4ANALYSIS_USE_ROOT
>> 422 HadrontherapyAnalysisManager* analysis = HadrontherapyAnalysisManager::GetInstance();
>> 423 analysis -> flush(); // Finalize the root file
>> 424 analysis -> book();
>> 425 #endif
>> 426 // Inform the kernel about the new geometry
>> 427 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
>> 428 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
>> 429
>> 430 PrintParameters();
499 } 431 }
500 432
501 ////////////////////////////////////////////// <<
502 void HadrontherapyDetectorConstruction::PrintP 433 void HadrontherapyDetectorConstruction::PrintParameters()
503 { 434 {
504 <<
505 G4cout << "The (X,Y,Z) dimensions of the p <<
506 G4BestUnit( phantom -> GetXHalfLength()*2. <<
507 G4BestUnit( phantom -> GetYHalfLength()*2. <<
508 G4BestUnit( phantom -> GetZHalfLength()*2. <<
509 <<
510 G4cout << "The (X,Y,Z) dimensions of the d <<
511 G4BestUnit( detector -> GetXHalfLength()*2 <<
512 G4BestUnit( detector -> GetYHalfLength()*2 <<
513 G4BestUnit( detector -> GetZHalfLength()*2 <<
514 <<
515 G4cout << "Displacement between Phantom an <<
516 G4cout << "DX= "<< G4BestUnit(phantomPosit <<
517 "DY= "<< G4BestUnit(phantomPosition.getY() <<
518 "DZ= "<< G4BestUnit(phantomPosition.getZ() <<
519 <<
520 G4cout << "The (X,Y,Z) sizes of the Voxels <<
521 G4BestUnit(sizeOfVoxelAlongX, "Length") < <<
522 G4BestUnit(sizeOfVoxelAlongY, "Length") < <<
523 G4BestUnit(sizeOfVoxelAlongZ, "Length") << <<
524 <<
525 G4cout << "The number of Voxels along (X,Y <<
526 numberOfVoxelsAlongX << ',' << <<
527 numberOfVoxelsAlongY <<',' << <<
528 numberOfVoxelsAlongZ << ')' << G4endl; <<
529 } <<
530 435
>> 436 G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
>> 437 G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
>> 438 G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
>> 439 G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
>> 440
>> 441 G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
>> 442 G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
>> 443 G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
>> 444 G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
>> 445
>> 446 G4cout << "Displacement between Phantom and World is: ";
>> 447 G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
>> 448 "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
>> 449 "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
>> 450
>> 451 G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
>> 452 G4BestUnit(sizeOfVoxelAlongX, "Length") << ',' <<
>> 453 G4BestUnit(sizeOfVoxelAlongY, "Length") << ',' <<
>> 454 G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
>> 455
>> 456 G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
>> 457 numberOfVoxelsAlongX << ',' <<
>> 458 numberOfVoxelsAlongY <<',' <<
>> 459 numberOfVoxelsAlongZ << ')' << G4endl;
531 460
>> 461 }
532 462