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