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