ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HadrontherapyDetectorConstruction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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"
55 
56 #include "HadrontherapyMatrix.hh"
57 
58 #include "HadrontherapyRBE.hh"
59 #include "G4SystemOfUnits.hh"
60 
61 #include <cmath>
62 
63 
64 
68 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
69 detectorSD(0), detectorROGeometry(0), matrix(0),
70 phantom(0), detector(0),
71 phantomLogicalVolume(0), detectorLogicalVolume(0),
72 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
73 aRegion(0)
74 {
75 
76 
77  /* NOTE! that the HadrontherapyDetectorConstruction class
78  * does NOT inherit from G4VUserDetectorConstruction G4 class
79  * So the Construct() mandatory virtual method is inside another geometric class
80  * like the passiveProtonBeamLIne, ...
81  */
82 
83  // Messenger to change parameters of the phantom/detector geometry
85 
86  // Default detector voxels size
87  // 200 slabs along the beam direction (X)
88  sizeOfVoxelAlongX = 200 *um;
89  sizeOfVoxelAlongY = 4 *cm;
90  sizeOfVoxelAlongZ = 4 *cm;
91 
92  // Define here the material of the water phantom and of the detector
93  SetPhantomMaterial("G4_WATER");
94  // Construct geometry (messenger commands)
95  // SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
96  SetDetectorSize(4. *cm, 4. *cm, 4. *cm);
97  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
98 
99  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
102  //GetDetectorToWorldPosition();
103 
104  // Write virtual parameters to the real ones and check for consistency
105  UpdateGeometry();
106 
107 
108 
109 }
110 
113 {
114  delete detectorROGeometry;
115  delete matrix;
116  delete detectorMessenger;
117 }
118 
121 {
122  return instance;
123 }
124 
126 // ConstructPhantom() is the method that construct a water box (called phantom
127 // (or water phantom) in the usual Medical physicists slang).
128 // A water phantom can be considered a good approximation of a an human body.
130 {
131  // Definition of the solid volume of the Phantom
132  phantom = new G4Box("Phantom",
133  phantomSizeX/2,
134  phantomSizeY/2,
135  phantomSizeZ/2);
136 
137  // Definition of the logical volume of the Phantom
140  "phantomLog", 0, 0, 0);
141 
142  // Definition of the physics volume of the Phantom
145  "phantomPhys",
147  motherPhys,
148  false,
149  0);
150 
151  // Visualisation attributes of the phantom
152  red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
153  red -> SetVisibility(true);
154  red -> SetForceSolid(true);
155  red -> SetForceWireframe(true);
156  phantomLogicalVolume -> SetVisAttributes(red);
157 }
158 
160 // ConstructDetector() is the method the reconstruct a detector region
161 // inside the water phantom. It is a volume, located inside the water phantom.
162 //
163 // **************************
164 // * water phantom *
165 // * *
166 // * *
167 // *--------------- *
168 // Beam * - *
169 // -----> * detector - *
170 // * - *
171 // *--------------- *
172 // * *
173 // * *
174 // * *
175 // **************************
176 //
177 // The detector can be dived in slices or voxelized
178 // and inside it different quantities (dose distribution, fluence distribution, LET, etc)
179 // can be stored.
181 
182 {
183  // Definition of the solid volume of the Detector
184  detector = new G4Box("Detector",
185 
186  phantomSizeX/2,
187 
188  phantomSizeY/2,
189 
190  phantomSizeZ/2);
191 
192  // Definition of the logic volume of the Phantom
195  "DetectorLog",
196  0,0,0);
197  // Definition of the physical volume of the Phantom
199  detectorPosition, // Setted by displacement
200  "DetectorPhys",
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  // **************
214  // Cut per Region
215  // **************
216  //
217  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
218  // required accuracy
219  if (!aRegion)
220  {
221  aRegion = new G4Region("DetectorLog");
222  detectorLogicalVolume -> SetRegion(aRegion);
224  }
225 }
226 
231  detectorToWorldPosition)
232 {
233  RO->Initialize(detectorToWorldPosition,
234  detectorSizeX/2,
235  detectorSizeY/2,
236  detectorSizeZ/2,
240 }
242 {
243 
244  //Virtual plane
246  NewSource= Varbool;
247  if(NewSource == true)
248  {
249  // std::cout<<"trr"<<std::endl;
251 
252  solidVirtualLayer = new G4Box("VirtualLayer",
253  1.*um,
254  20.*cm,
255  40.*cm);
256 
259  airNist,
260  "VirtualLayer");
261 
263  "VirtualLayer",
265  motherPhys,
266  false,
267  0);
268 
269  logicVirtualLayer -> SetVisAttributes(skyBlue);
270  }
271 
272 
273 
274 
275 }
276 
277 
280 {
281  // Check phantom/detector sizes & relative position
282  if (!IsInside(detectorSizeX,
285  phantomSizeX,
286  phantomSizeY,
287  phantomSizeZ,
289  ))
290  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
291 
292  // Check Detector sizes respect to the voxel ones
293 
295  G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
296  }
298  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
299  }
301  G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
302  }
303 }
304 
307 {
308 
309  if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
310  {
311  phantomMaterial = pMat;
312  detectorMaterial = pMat;
314  {
317 
318  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
320  G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
321  }
322  }
323  else
324  {
325  G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
326  " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
327  G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
328  return false;
329  }
330 
331  return true;
332 }
335 {
336  if (sizeX > 0.) phantomSizeX = sizeX;
337  if (sizeY > 0.) phantomSizeY = sizeY;
338  if (sizeZ > 0.) phantomSizeZ = sizeZ;
339 }
340 
343 {
344  if (sizeX > 0.) {detectorSizeX = sizeX;}
345  if (sizeY > 0.) {detectorSizeY = sizeY;}
346  if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
348 }
349 
352 {
353  if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
354  if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
355  if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
356 }
357 
360 {
362 }
363 
366 {
368 }
369 
371 {
372 
375 
376 }
379 {
380  /*
381  * Check parameters consistency
382  */
383  ParametersCheck();
384 
385  G4GeometryManager::GetInstance() -> OpenGeometry();
386  if (phantom)
387  {
388  phantom -> SetXHalfLength(phantomSizeX/2);
389  phantom -> SetYHalfLength(phantomSizeY/2);
390  phantom -> SetZHalfLength(phantomSizeZ/2);
391 
392  phantomPhysicalVolume -> SetTranslation(phantomPosition);
393  }
394  else ConstructPhantom();
395 
396 
397  // Get the center of the detector
399  if (detector)
400  {
401 
402  detector -> SetXHalfLength(detectorSizeX/2);
403  detector -> SetYHalfLength(detectorSizeY/2);
404  detector -> SetZHalfLength(detectorSizeZ/2);
405 
406  detectorPhysicalVolume -> SetTranslation(detectorPosition);
407  }
408  else ConstructDetector();
409 
410  //std::cout<<NewSource<<std::endl;
411  /*if(NewSource)
412  {
413  std::cout<<"via"<<std::endl;
414  }*/
415 
416 
417  // std::cout<<"i"<<std::endl;
418  // std::cout<<VirtualLayerPosition<<std::endl;
419  // physVirtualLayer->SetTranslation(VirtualLayerPosition);
420 
421 
422 
423 
424 
425  // Round to nearest integer number of voxel
426 
434 
436 
438 
439  //Set parameters, either for the Construct() or for the UpdateROGeometry()
441  detectorSizeX/2,
442  detectorSizeY/2,
443  detectorSizeZ/2,
447 
448  //This method below has an effect only if the RO geometry is already built.
449  RO->UpdateROGeometry();
450 
451 
452 
454  massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
455  // This will clear the existing matrix (together with all data inside it)!
459  massOfVoxel);
460 
461 
462  // Initialize RBE
464 
465  // Comment out the line below if let calculation is not needed!
466  // Initialize LET with energy of primaries and clear data inside
467  if ( (let = HadrontherapyLet::GetInstance(this)) )
468  {
470  }
471 
472 
473  // Initialize analysis
474  // Inform the kernel about the new geometry
476  G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
477 
478  PrintParameters();
479 
480  // CheckOverlaps();
481 }
482 
484 //Check of the geometry
487 {
489  G4cout << thePVStore->size() << " physical volumes are defined" << G4endl;
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]->CheckOverlaps(res,tol,fCheckOverlapsVerbosity) | overlapFlag;
496  overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,true) | overlapFlag; }
497  if (overlapFlag)
498  G4cout << "Check: there are overlapping volumes" << G4endl;
499 }
500 
503 {
504 
505  G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
506  G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
507  G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
508  G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
509 
510  G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
511  G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
512  G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
513  G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
514 
515  G4cout << "Displacement between Phantom and World is: ";
516  G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
517  "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
518  "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
519 
520  G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
521  G4BestUnit(sizeOfVoxelAlongX, "Length") << ',' <<
522  G4BestUnit(sizeOfVoxelAlongY, "Length") << ',' <<
523  G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
524 
525  G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
526  numberOfVoxelsAlongX << ',' <<
527  numberOfVoxelsAlongY <<',' <<
528  numberOfVoxelsAlongZ << ')' << G4endl;
529 }
530 
531