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