| Geant4 Cross Reference |
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"
48 #include "HadrontherapyDetectorConstruction.hh"
49 #include "HadrontherapyDetectorROGeometry.hh"
50 #include "HadrontherapyDetectorMessenger.hh"
51 #include "HadrontherapyDetectorSD.hh"
52 #include "HadrontherapyMatrix.hh"
53 #include "HadrontherapyLet.hh"
54 #include "PassiveProtonBeamLine.hh"
55 #include "BESTPassiveProtonBeamLine.hh"
56 #include "HadrontherapyMatrix.hh"
57
58 #include "HadrontherapyRBE.hh"
59 #include "G4SystemOfUnits.hh"
60
61 #include <cmath>
62
63
64
65 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::instance = 0;
66 /////////////////////////////////////////////////////////////////////////////
67 HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
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
84 detectorMessenger = new HadrontherapyDetectorMessenger(this);
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));
100 SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 0. *cm, 0. *cm));
101 SetDetectorPosition();
102 //GetDetectorToWorldPosition();
103
104 // Write virtual parameters to the real ones and check for consistency
105 UpdateGeometry();
106
107
108
109 }
110
111 /////////////////////////////////////////////////////////////////////////////
112 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
113 {
114 delete detectorROGeometry;
115 delete matrix;
116 delete detectorMessenger;
117 }
118
119 /////////////////////////////////////////////////////////////////////////////
120 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::GetInstance()
121 {
122 return instance;
123 }
124
125 /////////////////////////////////////////////////////////////////////////////
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.
129 void HadrontherapyDetectorConstruction::ConstructPhantom()
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
138 phantomLogicalVolume = new G4LogicalVolume(phantom,
139 phantomMaterial,
140 "phantomLog", 0, 0, 0);
141
142 // Definition of the physics volume of the Phantom
143 phantomPhysicalVolume = new G4PVPlacement(0,
144 phantomPosition,
145 "phantomPhys",
146 phantomLogicalVolume,
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
159 /////////////////////////////////////////////////////////////////////////////
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.
180 void HadrontherapyDetectorConstruction::ConstructDetector()
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
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 // **************
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);
223 aRegion->AddRootLogicalVolume( detectorLogicalVolume );
224 }
225 }
226
227 ///////////////////////////////////////////////////////////////////////
228 void HadrontherapyDetectorConstruction::InitializeDetectorROGeometry(
229 HadrontherapyDetectorROGeometry* RO,
230 G4ThreeVector
231 detectorToWorldPosition)
232 {
233 RO->Initialize(detectorToWorldPosition,
234 detectorSizeX/2,
235 detectorSizeY/2,
236 detectorSizeZ/2,
237 numberOfVoxelsAlongX,
238 numberOfVoxelsAlongY,
239 numberOfVoxelsAlongZ);
240 }
241 void HadrontherapyDetectorConstruction::VirtualLayer(G4bool Varbool)
242 {
243
244 //Virtual plane
245 VirtualLayerPosition = G4ThreeVector(0*cm,0*cm,0*cm);
246 NewSource= Varbool;
247 if(NewSource == true)
248 {
249 // std::cout<<"trr"<<std::endl;
250 G4Material* airNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
251
252 solidVirtualLayer = new G4Box("VirtualLayer",
253 1.*um,
254 20.*cm,
255 40.*cm);
256
257 logicVirtualLayer = new G4LogicalVolume(
258 solidVirtualLayer,
259 airNist,
260 "VirtualLayer");
261
262 physVirtualLayer= new G4PVPlacement(0,VirtualLayerPosition,
263 "VirtualLayer",
264 logicVirtualLayer,
265 motherPhys,
266 false,
267 0);
268
269 logicVirtualLayer -> SetVisAttributes(skyBlue);
270 }
271
272
273
274
275 }
276
277
278 ///////////////////////////////////////////////////////////////////////
279 void HadrontherapyDetectorConstruction::ParametersCheck()
280 {
281 // Check phantom/detector sizes & relative position
282 if (!IsInside(detectorSizeX,
283 detectorSizeY,
284 detectorSizeZ,
285 phantomSizeX,
286 phantomSizeY,
287 phantomSizeZ,
288 detectorToPhantomPosition
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
294 if ( detectorSizeX < sizeOfVoxelAlongX) {
295 G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error: Detector X size must be bigger or equal than that of Voxel X!");
296 }
297 if ( detectorSizeY < sizeOfVoxelAlongY) {
298 G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error: Detector Y size must be bigger or equal than that of Voxel Y!");
299 }
300 if ( detectorSizeZ < sizeOfVoxelAlongZ) {
301 G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error: Detector Z size must be bigger or equal than that of Voxel Z!");
302 }
303 }
304
305 ///////////////////////////////////////////////////////////////////////
306 G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
307 {
308
309 if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
310 {
311 phantomMaterial = pMat;
312 detectorMaterial = pMat;
313 if (detectorLogicalVolume && phantomLogicalVolume)
314 {
315 detectorLogicalVolume -> SetMaterial(pMat);
316 phantomLogicalVolume -> SetMaterial(pMat);
317
318 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
319 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
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 }
333 /////////////////////////////////////////////////////////////////////////////
334 void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
335 {
336 if (sizeX > 0.) phantomSizeX = sizeX;
337 if (sizeY > 0.) phantomSizeY = sizeY;
338 if (sizeZ > 0.) phantomSizeZ = sizeZ;
339 }
340
341 /////////////////////////////////////////////////////////////////////////////
342 void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
343 {
344 if (sizeX > 0.) {detectorSizeX = sizeX;}
345 if (sizeY > 0.) {detectorSizeY = sizeY;}
346 if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
347 SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
348 }
349
350 /////////////////////////////////////////////////////////////////////////////
351 void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
352 {
353 if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
354 if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
355 if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
356 }
357
358 ///////////////////////////////////////////////////////////////////////
359 void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
360 {
361 phantomPosition = pos;
362 }
363
364 /////////////////////////////////////////////////////////////////////////////
365 void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
366 {
367 detectorToPhantomPosition = displ;
368 }
369
370 void HadrontherapyDetectorConstruction::SetVirtualLayerPosition(G4ThreeVector position)
371 {
372
373 VirtualLayerPosition = position;
374 physVirtualLayer->SetTranslation(VirtualLayerPosition);
375
376 }
377 /////////////////////////////////////////////////////////////////////////////
378 void HadrontherapyDetectorConstruction::UpdateGeometry()
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
398 SetDetectorPosition();
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
427 numberOfVoxelsAlongX = G4lrint(detectorSizeX / sizeOfVoxelAlongX);
428 sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
429 numberOfVoxelsAlongY = G4lrint(detectorSizeY / sizeOfVoxelAlongY);
430 sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
431 numberOfVoxelsAlongZ = G4lrint(detectorSizeZ / sizeOfVoxelAlongZ);
432 sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
433 PassiveProtonBeamLine *ppbl= (PassiveProtonBeamLine*)
434
435 G4RunManager::GetRunManager()->GetUserDetectorConstruction();
436
437 HadrontherapyDetectorROGeometry* RO = (HadrontherapyDetectorROGeometry*) ppbl->GetParallelWorld(0);
438
439 //Set parameters, either for the Construct() or for the UpdateROGeometry()
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 the RO geometry is already built.
449 RO->UpdateROGeometry();
450
451
452
453 volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
454 massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
455 // This will clear the existing matrix (together with all data inside it)!
456 matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
457 numberOfVoxelsAlongY,
458 numberOfVoxelsAlongZ,
459 massOfVoxel);
460
461
462 // Initialize RBE
463 HadrontherapyRBE::CreateInstance(numberOfVoxelsAlongX, numberOfVoxelsAlongY, numberOfVoxelsAlongZ, massOfVoxel);
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 {
469 HadrontherapyLet::GetInstance() -> Initialize();
470 }
471
472
473 // Initialize analysis
474 // Inform the kernel about the new geometry
475 G4RunManager::GetRunManager() -> GeometryHasBeenModified();
476 G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
477
478 PrintParameters();
479
480 // CheckOverlaps();
481 }
482
483 /////////////////////////////////////////////////////////////////////////////
484 //Check of the geometry
485 /////////////////////////////////////////////////////////////////////////////
486 void HadrontherapyDetectorConstruction::CheckOverlaps()
487 {
488 G4PhysicalVolumeStore* thePVStore = G4PhysicalVolumeStore::GetInstance();
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
501 /////////////////////////////////////////////////////////////////////////////
502 void HadrontherapyDetectorConstruction::PrintParameters()
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
532