Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
HadrontherapyDetectorConstruction.cc
이 파일의 문서화 페이지로 가기
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 advanced example for Geant4
27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
28 
29 #include "G4UnitsTable.hh"
30 #include "G4SDManager.hh"
31 #include "G4RunManager.hh"
32 #include "G4GeometryManager.hh"
33 #include "G4SolidStore.hh"
34 #include "G4PhysicalVolumeStore.hh"
35 #include "G4LogicalVolumeStore.hh"
36 #include "G4Box.hh"
37 #include "G4LogicalVolume.hh"
38 #include "G4ThreeVector.hh"
39 #include "G4PVPlacement.hh"
40 #include "globals.hh"
41 #include "G4Transform3D.hh"
42 #include "G4RotationMatrix.hh"
43 #include "G4Colour.hh"
44 #include "G4UserLimits.hh"
45 #include "G4UnitsTable.hh"
46 #include "G4VisAttributes.hh"
47 #include "G4NistManager.hh"
52 #include "HadrontherapyMatrix.hh"
53 #include "HadrontherapyLet.hh"
54 #include "PassiveProtonBeamLine.hh"
56 #include "HadrontherapyMatrix.hh"
57 
58 #include "G4SystemOfUnits.hh"
59 
60 #include <cmath>
61 
62 
63 
67 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
68 detectorSD(0), detectorROGeometry(0), matrix(0),
69 phantom(0), detector(0),
70 phantomLogicalVolume(0), detectorLogicalVolume(0),
71 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
72 aRegion(0)
73 {
75 
76  /* NOTE! that the HadrontherapyDetectorConstruction class
77  * does NOT inherit from G4VUserDetectorConstruction G4 class
78  * So the Construct() mandatory virtual method is inside another geometric class
79  * like the passiveProtonBeamLIne, ...
80  */
81 
82  // Messenger to change parameters of the phantom/detector geometry
84 
85  // Default detector voxels size
86  // 200 slabs along the beam direction (X)
87  sizeOfVoxelAlongX = 200 *um;
88  sizeOfVoxelAlongY = 4 *cm;
89  sizeOfVoxelAlongZ = 4 *cm;
90 
91  // Define here the material of the water phantom and of the detector
92  SetPhantomMaterial("G4_WATER");
93  // Construct geometry (messenger commands)
94  // SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
95  SetDetectorSize(4. *cm, 4. *cm, 4. *cm);
96  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
97 
98  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
101  //GetDetectorToWorldPosition();
102 
103  // Write virtual parameters to the real ones and check for consistency
104  UpdateGeometry();
105 
106 }
107 
110 {
111  delete detectorROGeometry;
112  delete matrix;
113  delete detectorMessenger;
114 }
115 
118 {
119  return instance;
120 }
121 
123 // ConstructPhantom() is the method that construct a water box (called phantom
124 // (or water phantom) in the usual Medical physicists slang).
125 // A water phantom can be considered a good approximation of a an human body.
127 {
128  // Definition of the solid volume of the Phantom
129  phantom = new G4Box("Phantom",
130  phantomSizeX/2,
131  phantomSizeY/2,
132  phantomSizeZ/2);
133 
134  // Definition of the logical volume of the Phantom
137  "phantomLog", 0, 0, 0);
138 
139  // Definition of the physics volume of the Phantom
142  "phantomPhys",
144  motherPhys,
145  false,
146  0);
147 
148  // Visualisation attributes of the phantom
149  red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
150  red -> SetVisibility(true);
151  red -> SetForceSolid(true);
152  red -> SetForceWireframe(true);
153  phantomLogicalVolume -> SetVisAttributes(red);
154 }
155 
157 // ConstructDetector() is the method the reconstruct a detector region
158 // inside the water phantom. It is a volume, located inside the water phantom.
159 //
160 // **************************
161 // * water phantom *
162 // * *
163 // * *
164 // *--------------- *
165 // Beam * - *
166 // -----> * detector - *
167 // * - *
168 // *--------------- *
169 // * *
170 // * *
171 // * *
172 // **************************
173 //
174 // The detector can be dived in slices or voxelized
175 // and inside it different quantities (dose distribution, fluence distribution, LET, etc)
176 // can be stored.
178 
179 {
180  // Definition of the solid volume of the Detector
181  detector = new G4Box("Detector",
182 
183  phantomSizeX/2,
184 
185  phantomSizeY/2,
186 
187  phantomSizeZ/2);
188 
189  // Definition of the logic volume of the Phantom
192  "DetectorLog",
193  0,0,0);
194  // Definition of the physical volume of the Phantom
196  detectorPosition, // Setted by displacement
197  "DetectorPhys",
200  false,0);
201 
202  // Visualisation attributes of the detector
203  skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. , 235/255. ));
204  skyBlue -> SetVisibility(true);
205  skyBlue -> SetForceSolid(true);
206  //skyBlue -> SetForceWireframe(true);
207  detectorLogicalVolume -> SetVisAttributes(skyBlue);
208 
209  // **************
210  // Cut per Region
211  // **************
212  //
213  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
214  // required accuracy
215  if (!aRegion)
216  {
217  aRegion = new G4Region("DetectorLog");
218  detectorLogicalVolume -> SetRegion(aRegion);
220  }
221 }
222 
227  detectorToWorldPosition)
228 {
229  RO->Initialize(detectorToWorldPosition,
230  detectorSizeX/2,
231  detectorSizeY/2,
232  detectorSizeZ/2,
236 }
237 
240 {
241  // Check phantom/detector sizes & relative position
242  if (!IsInside(detectorSizeX,
245  phantomSizeX,
246  phantomSizeY,
247  phantomSizeZ,
249  ))
250  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
251 
252  // Check Detector sizes respect to the voxel ones
253 
255  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
256  }
258  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
259  }
261  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
262  }
263 }
264 
267 {
268 
269  if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
270  {
271  phantomMaterial = pMat;
272  detectorMaterial = pMat;
274  {
275  detectorLogicalVolume -> SetMaterial(pMat);
276  phantomLogicalVolume -> SetMaterial(pMat);
277 
278  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
280  G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
281  }
282  }
283  else
284  {
285  G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
286  " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
287  G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
288  return false;
289  }
290 
291  return true;
292 }
295 {
296  if (sizeX > 0.) phantomSizeX = sizeX;
297  if (sizeY > 0.) phantomSizeY = sizeY;
298  if (sizeZ > 0.) phantomSizeZ = sizeZ;
299 }
300 
303 {
304  if (sizeX > 0.) {detectorSizeX = sizeX;}
305  if (sizeY > 0.) {detectorSizeY = sizeY;}
306  if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
308 }
309 
312 {
313  if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
314  if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
315  if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
316 }
317 
320 {
322 }
323 
326 {
328 }
329 
332 {
333  /*
334  * Check parameters consistency
335  */
336  ParametersCheck();
337 
338  G4GeometryManager::GetInstance() -> OpenGeometry();
339  if (phantom)
340  {
341  phantom -> SetXHalfLength(phantomSizeX/2);
342  phantom -> SetYHalfLength(phantomSizeY/2);
343  phantom -> SetZHalfLength(phantomSizeZ/2);
344  phantomPhysicalVolume -> SetTranslation(phantomPosition);
345  }
346  else ConstructPhantom();
347 
348  // Get the center of the detector
350  if (detector)
351  {
352  detector -> SetXHalfLength(detectorSizeX/2);
353  detector -> SetYHalfLength(detectorSizeY/2);
354  detector -> SetZHalfLength(detectorSizeZ/2);
355  detectorPhysicalVolume -> SetTranslation(detectorPosition);
356  }
357  else ConstructDetector();
358 
359  // Round to nearest integer number of voxel
360 
368 
370 
372 
373  //Set parameters, either for the Construct() or for the UpdateROGeometry()
375  detectorSizeX/2,
376  detectorSizeY/2,
377  detectorSizeZ/2,
381 
382  //This method below has an effect only if the RO geometry is already built.
383  RO->UpdateROGeometry();
384 
385 
386 
388  massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
389  // This will clear the existing matrix (together with all data inside it)!
393  massOfVoxel);
394 
395 
396  // Comment out the line below if let calculation is not needed!
397  // Initialize LET with energy of primaries and clear data inside
398  if ( (let = HadrontherapyLet::GetInstance(this)) )
399  {
401  }
402 
403 
404  // Initialize analysis
405  // Inform the kernel about the new geometry
407  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
408 
409  PrintParameters();
410  // CheckOverlaps();
411 }
412 
414 //Check of the geometry
417 {
419  G4cout << thePVStore->size() << " physical volumes are defined" << G4endl;
420  G4bool overlapFlag = false;
421  G4int res=1000;
422  G4double tol=0.; //tolerance
423  for (size_t i=0;i<thePVStore->size();i++)
424  {
425  //overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,fCheckOverlapsVerbosity) | overlapFlag;
426  overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,true) | overlapFlag; }
427  if (overlapFlag)
428  G4cout << "Check: there are overlapping volumes" << G4endl;
429 }
430 
433 {
434 
435  G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
436  G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
437  G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
438  G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
439 
440  G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
441  G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
442  G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
443  G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
444 
445  G4cout << "Displacement between Phantom and World is: ";
446  G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
447  "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
448  "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
449 
450  G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
451  G4BestUnit(sizeOfVoxelAlongX, "Length") << ',' <<
452  G4BestUnit(sizeOfVoxelAlongY, "Length") << ',' <<
453  G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
454 
455  G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
456  numberOfVoxelsAlongX << ',' <<
457  numberOfVoxelsAlongY <<',' <<
458  numberOfVoxelsAlongZ << ')' << G4endl;
459 }
460 
void AddRootLogicalVolume(G4LogicalVolume *lv)
Definition: G4Region.cc:290
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:80
G4VUserParallelWorld * GetParallelWorld(G4int i) const
CLHEP::Hep3Vector G4ThreeVector
static const G4double pos
static HadrontherapyMatrix * GetInstance()
double getZ() const
void SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
#define G4endl
Definition: G4ios.hh:61
void Initialize()
Definition: errprop.cc:101
HadrontherapyDetectorMessenger * detectorMessenger
bool IsInside(G4double detectorX, G4double detectorY, G4double detectorZ, G4double phantomX, G4double phantomY, G4double phantomZ, G4ThreeVector pos)
void InitializeDetectorROGeometry(HadrontherapyDetectorROGeometry *, G4ThreeVector detectorToWorldPosition)
static HadrontherapyLet * GetInstance()
static constexpr double um
Definition: G4SIunits.hh:113
HadrontherapyDetectorROGeometry * detectorROGeometry
static HadrontherapyDetectorConstruction * GetInstance()
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static HadrontherapyDetectorConstruction * instance
static G4GeometryManager * GetInstance()
Definition: G4Box.hh:64
double getX() const
double getY() const
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.hh:65
static G4PhysicalVolumeStore * GetInstance()
int G4lrint(double ad)
Definition: templates.hh:151
int G4int
Definition: G4Types.hh:78
void SetDetectorToPhantomPosition(G4ThreeVector DetectorToPhantomPosition)
G4GLOB_DLL std::ostream G4cout
static constexpr double cm
Definition: G4SIunits.hh:119
void SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
void SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
static HadrontherapyAnalysisManager * GetInstance()
const G4VUserDetectorConstruction * GetUserDetectorConstruction() const
static G4NistManager * Instance()
void Initialize(G4ThreeVector detectorPos, G4double detectorDimX, G4double detectorDimY, G4double detectorDimZ, G4int numberOfVoxelsX, G4int numberOfVoxelsY, G4int numberOfVoxelsZ)