Geant4_10
HadrontherapyDetectorConstruction.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // This is the *BASIC* version of Hadrontherapy, a Geant4-based application
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 
36 #include "G4UnitsTable.hh"
37 #include "G4SDManager.hh"
38 #include "G4RunManager.hh"
39 #include "G4GeometryManager.hh"
40 #include "G4SolidStore.hh"
41 #include "G4PhysicalVolumeStore.hh"
42 #include "G4LogicalVolumeStore.hh"
43 #include "G4Box.hh"
44 #include "G4LogicalVolume.hh"
45 #include "G4ThreeVector.hh"
46 #include "G4PVPlacement.hh"
47 #include "globals.hh"
48 #include "G4Transform3D.hh"
49 #include "G4RotationMatrix.hh"
50 #include "G4Colour.hh"
51 #include "G4UserLimits.hh"
52 #include "G4UnitsTable.hh"
53 #include "G4VisAttributes.hh"
54 #include "G4NistManager.hh"
55 
60 #include "HadrontherapyMatrix.hh"
61 #include "HadrontherapyLet.hh"
62 #include "PassiveProtonBeamLine.hh"
64 #include "G4SystemOfUnits.hh"
65 
66 #include <cmath>
67 
68 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::instance = 0;
71  : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
72  detectorSD(0), detectorROGeometry(0), matrix(0),
73  phantom(0), detector(0),
74  phantomLogicalVolume(0), detectorLogicalVolume(0),
75  phantomPhysicalVolume(0), detectorPhysicalVolume(0),
76  aRegion(0)
77 {
79 
80  /* NOTE! that the HadrontherapyDetectorConstruction class
81  * does NOT inherit from G4VUserDetectorConstruction G4 class
82  * So the Construct() mandatory virtual method is inside another geometric class
83  * like the passiveProtonBeamLIne, ...
84  */
85 
86  // Messenger to change parameters of the phantom/detector geometry
87  detectorMessenger = new HadrontherapyDetectorMessenger(this);
88 
89  // Default detector voxels size
90  // 200 slabs along the beam direction (X)
91  sizeOfVoxelAlongX = 200 *um;
92  sizeOfVoxelAlongY = 4 *cm;
93  sizeOfVoxelAlongZ = 4 *cm;
94 
95  // Define here the material of the water phantom and of the detector
96  SetPhantomMaterial("G4_WATER");
97  // Construct geometry (messenger commands)
98  SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
99  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
100  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
103 //GetDetectorToWorldPosition();
104 
105 
106  // Write virtual parameters to the real ones and check for consistency
107  UpdateGeometry();
108 }
109 
112 {
113  delete detectorROGeometry;
114  delete matrix;
115  delete detectorMessenger;
116 }
117 
118 
120 {
121 
122  return instance;
123 }
124 
125 
127 // ConstructPhantom() is the method that reconstuct a water box (called phantom
128 // (or water phantom) in the usual Medical physicists slang).
129 // A water phantom can be considered a good
130 // approximation of a an human body.
131 void HadrontherapyDetectorConstruction::ConstructPhantom()
132 {
133  // Definition of the solid volume of the Phantom
134  phantom = new G4Box("Phantom",
135  phantomSizeX/2,
136  phantomSizeY/2,
137  phantomSizeZ/2);
138 
139 // Definition of the logical volume of the Phantom
140  phantomLogicalVolume = new G4LogicalVolume(phantom,
141  phantomMaterial,
142  "phantomLog", 0, 0, 0);
143 
144  // Definition of the physics volume of the Phantom
145  phantomPhysicalVolume = new G4PVPlacement(0,
146  phantomPosition,
147  "phantomPhys",
148  phantomLogicalVolume,
149  motherPhys,
150  false,
151  0);
152 
153 // Visualisation attributes of the phantom
154  red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
155  red -> SetVisibility(true);
156  red -> SetForceSolid(true);
157  red -> SetForceWireframe(true);
158  phantomLogicalVolume -> SetVisAttributes(red);
159 }
160 
162 // ConstructDetector() it the method the reconstruct a detector region
163 // inside the water phantom. It is a volume, located inside the water phantom
164 // and with two coincident faces:
165 //
166 // **************************
167 // * water phantom *
168 // * *
169 // * *
170 // *--------------- *
171 // Beam * - *
172 // -----> * detector - *
173 // * - *
174 // *--------------- *
175 // * *
176 // * *
177 // * *
178 // **************************
179 //
180 // The detector is the volume that can be dived in slices or voxelized
181 // and in it we can collect a number of usefull information:
182 // dose distribution, fluence distribution, LET and so on
183 void HadrontherapyDetectorConstruction::ConstructDetector()
184 {
185 
186  // Definition of the solid volume of the Detector
187  detector = new G4Box("Detector",
188  detectorSizeX/2,
189  detectorSizeY/2,
190  detectorSizeZ/2);
191 
192  // Definition of the logic volume of the Phantom
193  detectorLogicalVolume = new G4LogicalVolume(detector,
194  detectorMaterial,
195  "DetectorLog",
196  0,0,0);
197 // Definition of the physical volume of the Phantom
198  detectorPhysicalVolume = new G4PVPlacement(0,
199  detectorPosition, // Setted by displacement
200  "DetectorPhys",
201  detectorLogicalVolume,
202  phantomPhysicalVolume,
203  false,0);
204 
205 // Visualisation attributes of the detector
206  skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. , 235/255. ));
207  skyBlue -> SetVisibility(true);
208  skyBlue -> SetForceSolid(true);
209  //skyBlue -> SetForceWireframe(true);
210  detectorLogicalVolume -> SetVisAttributes(skyBlue);
211 
212  // **************
213  // Cut per Region
214  // **************
215 
216  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
217  // required accuracy
218  if (!aRegion)
219  {
220  aRegion = new G4Region("DetectorLog");
221  detectorLogicalVolume -> SetRegion(aRegion);
222  aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
223  }
224 
225 }
226 
227 
229 
233  detectorToWorldPosition)
234 {
235  RO->Initialize(detectorToWorldPosition,
236  detectorSizeX/2,
237  detectorSizeY/2,
238  detectorSizeZ/2,
239  numberOfVoxelsAlongX,
240  numberOfVoxelsAlongY,
241  numberOfVoxelsAlongZ);
242 }
243 
244 
245 
246 
247 
249 
250 void HadrontherapyDetectorConstruction::ParametersCheck()
251 {
252  // Check phantom/detector sizes & relative position
253  if (!IsInside(detectorSizeX,
254  detectorSizeY,
255  detectorSizeZ,
256  phantomSizeX,
257  phantomSizeY,
258  phantomSizeZ,
259  detectorToPhantomPosition
260  ))
261  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
262 
263  // Check Detector sizes respect to the voxel ones
264 
265  if ( detectorSizeX < sizeOfVoxelAlongX) {
266  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
267  }
268  if ( detectorSizeY < sizeOfVoxelAlongY) {
269  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
270  }
271  if ( detectorSizeZ < sizeOfVoxelAlongZ) {
272  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
273  }
274 
275 }
277 // MESSENGERS //
279 
281 {
282 
283  if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
284  {
285  phantomMaterial = pMat;
286  detectorMaterial = pMat;
287  if (detectorLogicalVolume && phantomLogicalVolume)
288  {
289  detectorLogicalVolume -> SetMaterial(pMat);
290  phantomLogicalVolume -> SetMaterial(pMat);
291 
292  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
293  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
294  G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
295  }
296  }
297  else
298  {
299  G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
300  " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
301  G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
302  return false;
303  }
304 
305  return true;
306 }
309 {
310  if (sizeX > 0.) phantomSizeX = sizeX;
311  if (sizeY > 0.) phantomSizeY = sizeY;
312  if (sizeZ > 0.) phantomSizeZ = sizeZ;
313 }
317 {
318  if (sizeX > 0.) {detectorSizeX = sizeX;}
319  if (sizeY > 0.) {detectorSizeY = sizeY;}
320  if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
321  SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
322 }
324 
326 {
327  if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
328  if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
329  if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
330 }
332 {
333  phantomPosition = pos;
334 }
335 
338 {
339  detectorToPhantomPosition = displ;
340 }
343 {
344  /*
345  * Check parameters consistency
346  */
347  ParametersCheck();
348 
349  G4GeometryManager::GetInstance() -> OpenGeometry();
350  if (phantom)
351  {
352  phantom -> SetXHalfLength(phantomSizeX/2);
353  phantom -> SetYHalfLength(phantomSizeY/2);
354  phantom -> SetZHalfLength(phantomSizeZ/2);
355  phantomPhysicalVolume -> SetTranslation(phantomPosition);
356  }
357  else ConstructPhantom();
358 
359  // Get the center of the detector
361  if (detector)
362  {
363  detector -> SetXHalfLength(detectorSizeX/2);
364  detector -> SetYHalfLength(detectorSizeY/2);
365  detector -> SetZHalfLength(detectorSizeZ/2);
366  detectorPhysicalVolume -> SetTranslation(detectorPosition);
367  }
368  else ConstructDetector();
369 
370  // Round to nearest integer number of voxel
371 
372  numberOfVoxelsAlongX = G4lrint(detectorSizeX / sizeOfVoxelAlongX);
373 
374  sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
375 
376  numberOfVoxelsAlongY = G4lrint(detectorSizeY / sizeOfVoxelAlongY);
377 
378  sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
379 
380  numberOfVoxelsAlongZ = G4lrint(detectorSizeZ / sizeOfVoxelAlongZ);
381 
382  sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
383 
385 
387 
388 
390 
391  //Set parameters, either for the Construct() or for the UpdateROGeometry()
393  detectorSizeX/2,
394  detectorSizeY/2,
395  detectorSizeZ/2,
396  numberOfVoxelsAlongX,
397  numberOfVoxelsAlongY,
398  numberOfVoxelsAlongZ);
399 
400  //This method below has an effect only if the RO geometry is already built.
401  RO->UpdateROGeometry();
402 
403  volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
404  massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
405  // This will clear the existing matrix (together with all data inside it)!
406  matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
407  numberOfVoxelsAlongY,
408  numberOfVoxelsAlongZ,
409  massOfVoxel);
410 
411 
412 // Comment out the line below if let calculation is not needed!
413  // Initialize LET with energy of primaries and clear data inside
414  if ( (let = HadrontherapyLet::GetInstance(this)) )
415  {
417  }
418 
419 
420  // Initialize analysis
421 #ifdef G4ANALYSIS_USE_ROOT
423  analysis -> flush(); // Finalize the root file
424  analysis -> book();
425 #endif
426  // Inform the kernel about the new geometry
427  G4RunManager::GetRunManager() -> GeometryHasBeenModified();
428  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
429 
430  PrintParameters();
431 }
432 
434 {
435 
436  G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
437  G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
438  G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
439  G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
440 
441  G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
442  G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
443  G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
444  G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
445 
446  G4cout << "Displacement between Phantom and World is: ";
447  G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
448  "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
449  "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
450 
451  G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
452  G4BestUnit(sizeOfVoxelAlongX, "Length") << ',' <<
453  G4BestUnit(sizeOfVoxelAlongY, "Length") << ',' <<
454  G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
455 
456  G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
457  numberOfVoxelsAlongX << ',' <<
458  numberOfVoxelsAlongY <<',' <<
459  numberOfVoxelsAlongZ << ')' << G4endl;
460 
461 }
static HadrontherapyAnalysisManager * GetInstance()
void SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
CLHEP::Hep3Vector G4ThreeVector
Definition: G4Box.hh:63
static HadrontherapyLet * GetInstance()
const G4VUserDetectorConstruction * GetUserDetectorConstruction() const
void Initialize()
Definition: errprop.cc:96
void Initialize(G4ThreeVector detectorPos, G4double detectorDimX, G4double detectorDimY, G4double detectorDimZ, G4int numberOfVoxelsX, G4int numberOfVoxelsY, G4int numberOfVoxelsZ)
double getY() const
static HadrontherapyDetectorConstruction * GetInstance()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
static G4NistManager * Instance()
string material
Definition: eplot.py:19
static HadrontherapyMatrix * GetInstance()
double getX() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
void SetDetectorToPhantomPosition(G4ThreeVector DetectorToPhantomPosition)
bool IsInside(G4double detectorX, G4double detectorY, G4double detectorZ, G4double phantomX, G4double phantomY, G4double phantomZ, G4ThreeVector pos)
G4VUserParallelWorld * GetParallelWorld(G4int i) const
static G4GeometryManager * GetInstance()
def SetMaterial
Definition: EmPlot.py:25
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:74
int G4lrint(double ad)
Definition: templates.hh:163
double getZ() const
#define G4endl
Definition: G4ios.hh:61
void InitializeDetectorROGeometry(HadrontherapyDetectorROGeometry *, G4ThreeVector detectorToWorldPosition)
double G4double
Definition: G4Types.hh:76
void SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)