Geant4  10.01.p02
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 // Hadrontherapy development is supported by the
27 // Italian Institute for Nuclear Physics (INFN)
28 // in the framework of the MC-INFN Project
29 
30 #include "G4UnitsTable.hh"
31 #include "G4SDManager.hh"
32 #include "G4RunManager.hh"
33 #include "G4GeometryManager.hh"
34 #include "G4SolidStore.hh"
35 #include "G4PhysicalVolumeStore.hh"
36 #include "G4LogicalVolumeStore.hh"
37 #include "G4Box.hh"
38 #include "G4LogicalVolume.hh"
39 #include "G4ThreeVector.hh"
40 #include "G4PVPlacement.hh"
41 #include "globals.hh"
42 #include "G4Transform3D.hh"
43 #include "G4RotationMatrix.hh"
44 #include "G4Colour.hh"
45 #include "G4UserLimits.hh"
46 #include "G4UnitsTable.hh"
47 #include "G4VisAttributes.hh"
48 #include "G4NistManager.hh"
53 #include "HadrontherapyMatrix.hh"
54 #include "HadrontherapyLet.hh"
55 #include "PassiveProtonBeamLine.hh"
57 #include "G4SystemOfUnits.hh"
58 
59 #include <cmath>
60 
64 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
65 detectorSD(0), detectorROGeometry(0), matrix(0),
66 phantom(0), detector(0),
67 phantomLogicalVolume(0), detectorLogicalVolume(0),
68 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
69 aRegion(0)
70 {
72 
73  /* NOTE! that the HadrontherapyDetectorConstruction class
74  * does NOT inherit from G4VUserDetectorConstruction G4 class
75  * So the Construct() mandatory virtual method is inside another geometric class
76  * like the passiveProtonBeamLIne, ...
77  */
78 
79  // Messenger to change parameters of the phantom/detector geometry
81 
82  // Default detector voxels size
83  // 200 slabs along the beam direction (X)
84  sizeOfVoxelAlongX = 200 *um;
85  sizeOfVoxelAlongY = 4 *cm;
86  sizeOfVoxelAlongZ = 4 *cm;
87 
88  // Define here the material of the water phantom and of the detector
89  SetPhantomMaterial("G4_WATER");
90  // Construct geometry (messenger commands)
91  SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
92  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
93  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
96  //GetDetectorToWorldPosition();
97 
98  // Write virtual parameters to the real ones and check for consistency
100 }
101 
104 {
105  delete detectorROGeometry;
106  delete matrix;
107  delete detectorMessenger;
108 }
109 
112 {
113  return instance;
114 }
115 
117 // ConstructPhantom() is the method that construct a water box (called phantom
118 // (or water phantom) in the usual Medical physicists slang).
119 // A water phantom can be considered a good approximation of a an human body.
121 {
122  // Definition of the solid volume of the Phantom
123  phantom = new G4Box("Phantom",
124  phantomSizeX/2,
125  phantomSizeY/2,
126  phantomSizeZ/2);
127 
128  // Definition of the logical volume of the Phantom
131  "phantomLog", 0, 0, 0);
132 
133  // Definition of the physics volume of the Phantom
136  "phantomPhys",
138  motherPhys,
139  false,
140  0);
141 
142  // Visualisation attributes of the phantom
143  red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
144  red -> SetVisibility(true);
145  red -> SetForceSolid(true);
146  red -> SetForceWireframe(true);
147  phantomLogicalVolume -> SetVisAttributes(red);
148 }
149 
151 // ConstructDetector() is the method the reconstruct a detector region
152 // inside the water phantom. It is a volume, located inside the water phantom.
153 //
154 // **************************
155 // * water phantom *
156 // * *
157 // * *
158 // *--------------- *
159 // Beam * - *
160 // -----> * detector - *
161 // * - *
162 // *--------------- *
163 // * *
164 // * *
165 // * *
166 // **************************
167 //
168 // The detector can be dived in slices or voxelized
169 // and inside it different quantities (dose distribution, fluence distribution, LET, etc)
170 // can be stored.
172 {
173  // Definition of the solid volume of the Detector
174  detector = new G4Box("Detector",
175  detectorSizeX/2,
176  detectorSizeY/2,
177  detectorSizeZ/2);
178 
179  // Definition of the logic volume of the Phantom
182  "DetectorLog",
183  0,0,0);
184  // Definition of the physical volume of the Phantom
186  detectorPosition, // Setted by displacement
187  "DetectorPhys",
190  false,0);
191 
192  // Visualisation attributes of the detector
193  skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. , 235/255. ));
194  skyBlue -> SetVisibility(true);
195  skyBlue -> SetForceSolid(true);
196  //skyBlue -> SetForceWireframe(true);
197  detectorLogicalVolume -> SetVisAttributes(skyBlue);
198 
199  // **************
200  // Cut per Region
201  // **************
202  //
203  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
204  // required accuracy
205  if (!aRegion)
206  {
207  aRegion = new G4Region("DetectorLog");
208  detectorLogicalVolume -> SetRegion(aRegion);
209  aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
210  }
211 }
212 
217  detectorToWorldPosition)
218 {
219  RO->Initialize(detectorToWorldPosition,
220  detectorSizeX/2,
221  detectorSizeY/2,
222  detectorSizeZ/2,
226 }
227 
230 {
231  // Check phantom/detector sizes & relative position
232  if (!IsInside(detectorSizeX,
235  phantomSizeX,
236  phantomSizeY,
237  phantomSizeZ,
239  ))
240  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
241 
242  // Check Detector sizes respect to the voxel ones
243 
245  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
246  }
248  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
249  }
251  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
252  }
253 }
254 
257 {
258 
259  if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
260  {
261  phantomMaterial = pMat;
262  detectorMaterial = pMat;
264  {
265  detectorLogicalVolume -> SetMaterial(pMat);
266  phantomLogicalVolume -> SetMaterial(pMat);
267 
268  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
270  G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
271  }
272  }
273  else
274  {
275  G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
276  " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
277  G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
278  return false;
279  }
280 
281  return true;
282 }
285 {
286  if (sizeX > 0.) phantomSizeX = sizeX;
287  if (sizeY > 0.) phantomSizeY = sizeY;
288  if (sizeZ > 0.) phantomSizeZ = sizeZ;
289 }
290 
293 {
294  if (sizeX > 0.) {detectorSizeX = sizeX;}
295  if (sizeY > 0.) {detectorSizeY = sizeY;}
296  if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
298 }
299 
302 {
303  if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
304  if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
305  if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
306 }
307 
310 {
312 }
313 
316 {
318 }
319 
322 {
323  /*
324  * Check parameters consistency
325  */
326  ParametersCheck();
327 
328  G4GeometryManager::GetInstance() -> OpenGeometry();
329  if (phantom)
330  {
331  phantom -> SetXHalfLength(phantomSizeX/2);
332  phantom -> SetYHalfLength(phantomSizeY/2);
333  phantom -> SetZHalfLength(phantomSizeZ/2);
334  phantomPhysicalVolume -> SetTranslation(phantomPosition);
335  }
336  else ConstructPhantom();
337 
338  // Get the center of the detector
340  if (detector)
341  {
342  detector -> SetXHalfLength(detectorSizeX/2);
343  detector -> SetYHalfLength(detectorSizeY/2);
344  detector -> SetZHalfLength(detectorSizeZ/2);
345  detectorPhysicalVolume -> SetTranslation(detectorPosition);
346  }
347  else ConstructDetector();
348 
349  // Round to nearest integer number of voxel
350 
358 
360 
362 
363  //Set parameters, either for the Construct() or for the UpdateROGeometry()
365  detectorSizeX/2,
366  detectorSizeY/2,
367  detectorSizeZ/2,
371 
372  //This method below has an effect only if the RO geometry is already built.
373  RO->UpdateROGeometry();
374 
375 
376 
378  massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
379  // This will clear the existing matrix (together with all data inside it)!
383  massOfVoxel);
384 
385 
386  // Comment out the line below if let calculation is not needed!
387  // Initialize LET with energy of primaries and clear data inside
388  if ( (let = HadrontherapyLet::GetInstance(this)) )
389  {
391  }
392 
393 
394  // Initialize analysis
395 #ifdef G4ANALYSIS_USE_ROOT
397  analysis -> flush(); // Finalize the root file
398  analysis -> book();
399 #endif
400  // Inform the kernel about the new geometry
402  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
403 
404  PrintParameters();
405  // CheckOverlaps();
406 }
407 
409 //Check of the geometry
412 {
414  G4cout << thePVStore->size() << " physical volumes are defined" << G4endl;
415  G4bool overlapFlag = false;
416  G4int res=1000;
417  G4double tol=0.; //tolerance
418  for (size_t i=0;i<thePVStore->size();i++)
419  {
420  //overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,fCheckOverlapsVerbosity) | overlapFlag;
421  overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,true) | overlapFlag; }
422  if (overlapFlag)
423  G4cout << "Check: there are overlapping volumes" << G4endl;
424 }
425 
428 {
429 
430  G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
431  G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
432  G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
433  G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
434 
435  G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
436  G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
437  G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
438  G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
439 
440  G4cout << "Displacement between Phantom and World is: ";
441  G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
442  "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
443  "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
444 
445  G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
446  G4BestUnit(sizeOfVoxelAlongX, "Length") << ',' <<
447  G4BestUnit(sizeOfVoxelAlongY, "Length") << ',' <<
448  G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
449 
450  G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
451  numberOfVoxelsAlongX << ',' <<
452  numberOfVoxelsAlongY <<',' <<
453  numberOfVoxelsAlongZ << ')' << G4endl;
454 }
455 
static const double cm
Definition: G4SIunits.hh:106
static HadrontherapyAnalysisManager * GetInstance()
Get the pointer to the analysis manager.
void SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
static HadrontherapyDetectorConstruction * instance
HadrontherapyDetectorMessenger * detectorMessenger
CLHEP::Hep3Vector G4ThreeVector
Definition: G4Box.hh:64
static HadrontherapyLet * GetInstance()
const G4VUserDetectorConstruction * GetUserDetectorConstruction() const
void Initialize()
Definition: errprop.cc:101
void Initialize(G4ThreeVector detectorPos, G4double detectorDimX, G4double detectorDimY, G4double detectorDimZ, G4int numberOfVoxelsX, G4int numberOfVoxelsY, G4int numberOfVoxelsZ)
static HadrontherapyDetectorConstruction * GetInstance()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
static G4PhysicalVolumeStore * GetInstance()
static HadrontherapyMatrix * GetInstance()
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)
HadrontherapyDetectorROGeometry * detectorROGeometry
G4VUserParallelWorld * GetParallelWorld(G4int i) const
static G4GeometryManager * GetInstance()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
int G4lrint(double ad)
Definition: templates.hh:163
A class for connecting the simulation to an analysis package.
#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)
static const G4double pos