source: trunk/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc @ 1321

Last change on this file since 1321 was 1313, checked in by garnier, 14 years ago

geant4.9.4 beta rc0
File size: 16.3 KB
Line 
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// HadrontherapyDetectorConstruction.cc
27//
28// See more at: http://g4advancedexamples.lngs.infn.it/Examples/hadrontherapy
29
30
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"
49
50#include "HadrontherapyDetectorConstruction.hh"
51#include "HadrontherapyDetectorROGeometry.hh"
52#include "HadrontherapyDetectorMessenger.hh"
53#include "HadrontherapyDetectorSD.hh"
54#include "HadrontherapyMatrix.hh"
55#include "HadrontherapyAnalysisManager.hh"
56
57#include <cmath>
58
59/////////////////////////////////////////////////////////////////////////////
60HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
61  : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
62    detectorSD(0), detectorROGeometry(0), matrix(0), 
63    phantom(0), detector(0),
64    phantomLogicalVolume(0), detectorLogicalVolume(0), 
65    phantomPhysicalVolume(0), detectorPhysicalVolume(0),
66    aRegion(0)
67{
68  HadrontherapyAnalysisManager::GetInstance();
69
70  // NOTE! that the HadrontherapyDetectorConstruction class
71  // does NOT inherit from G4VUserDetectorConstruction G4 class
72  // So the Construct() mandatory virtual method is inside another geometric class
73  // (like the passiveProtonBeamLIne, ...)
74
75  // Messenger to change parameters of the phantom/detector geometry
76  detectorMessenger = new HadrontherapyDetectorMessenger(this);
77
78  // Default detector voxels size
79  // 200 slabs along the beam direction (X)
80  sizeOfVoxelAlongX = 200 *um; 
81  sizeOfVoxelAlongY = 4 *cm; 
82  sizeOfVoxelAlongZ = 4 *cm; 
83
84  // Define here the material of the water phantom and of the detector
85  SetPhantomMaterial("G4_WATER"); 
86  // Construct geometry (messenger commands)
87  SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
88  SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
89  SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
90  SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 18. *cm, 18. *cm));
91
92
93  // Write virtual parameters to the real ones and check for consistency     
94  UpdateGeometry();
95}
96
97/////////////////////////////////////////////////////////////////////////////
98HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
99{ 
100    delete detectorROGeometry; 
101    delete matrix; 
102    delete detectorMessenger;
103}
104
105/////////////////////////////////////////////////////////////////////////////
106// ConstructPhantom() is the method that reconstuct a water box (called phantom
107// (or water phantom) in the usual Medical physicists slang).
108// A water phantom can be considered a good
109// approximation of a an human body.
110void HadrontherapyDetectorConstruction::ConstructPhantom()
111{
112    // Definition of the solid volume of the Phantom
113    phantom = new G4Box("Phantom", 
114                        phantomSizeX/2, 
115                        phantomSizeY/2, 
116                        phantomSizeZ/2);
117   
118// Definition of the logical volume of the Phantom
119    phantomLogicalVolume = new G4LogicalVolume(phantom, 
120                                             phantomMaterial, 
121                                             "phantomLog", 0, 0, 0);
122 
123    // Definition of the physics volume of the Phantom
124    phantomPhysicalVolume = new G4PVPlacement(0,
125                                            phantomPosition,
126                                            "phantomPhys",
127                                            phantomLogicalVolume,
128                                            motherPhys,
129                                            false,
130                                            0);
131
132// Visualisation attributes of the phantom
133    red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
134    red -> SetVisibility(true);
135    red -> SetForceSolid(true);
136    //red -> SetForceWireframe(true);
137    phantomLogicalVolume -> SetVisAttributes(red);
138}
139
140/////////////////////////////////////////////////////////////////////////////
141// ConstructDetector() it the method the reconstruct a detector region
142// inside the water phantom. It is a volume, located inside the water phantom
143// and with two coincident faces:
144//
145//           **************************
146//           *   water phantom        *
147//           *                        *
148//           *                        *
149//           *---------------         *
150//  Beam     *              -         *
151//  ----->   * detector     -         *
152//           *              -         *
153//           *---------------         *
154//           *                        *
155//           *                        *
156//           *                        *
157//           **************************
158//
159// The detector is the volume that can be dived in slices or voxelized
160// and in it we can collect a number of usefull information:
161// dose distribution, fluence distribution, LET and so on
162void HadrontherapyDetectorConstruction::ConstructDetector()
163{
164
165    // Definition of the solid volume of the Detector
166    detector = new G4Box("Detector", 
167                         detectorSizeX/2, 
168                         detectorSizeY/2, 
169                         detectorSizeZ/2);
170   
171    // Definition of the logic volume of the Phantom
172    detectorLogicalVolume = new G4LogicalVolume(detector,
173                                                detectorMaterial,
174                                                "DetectorLog",
175                                                0,0,0);
176// Definition of the physical volume of the Phantom
177    detectorPhysicalVolume = new G4PVPlacement(0, 
178                                               detectorPosition, // Setted by displacement
179                                               "DetectorPhys", 
180                                               detectorLogicalVolume, 
181                                               phantomPhysicalVolume, 
182                                               false,0);
183
184// Visualisation attributes of the detector
185    skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
186    skyBlue -> SetVisibility(true);
187    skyBlue -> SetForceSolid(true);
188    //skyBlue -> SetForceWireframe(true);
189    detectorLogicalVolume -> SetVisAttributes(skyBlue);
190
191  // **************
192  // Cut per Region
193  // **************
194 
195  // A smaller cut is fixed in the phantom to calculate the energy deposit with the
196  // required accuracy
197    if (!aRegion)
198    {
199        aRegion = new G4Region("DetectorLog");
200        detectorLogicalVolume -> SetRegion(aRegion);
201        aRegion -> AddRootLogicalVolume(detectorLogicalVolume);
202    }
203
204}
205
206/////////////////////////////////////////////////////////////////////////////
207
208void  HadrontherapyDetectorConstruction::ConstructSensitiveDetector(G4ThreeVector detectorToWorldPosition)
209{ 
210    // Install new Sensitive Detector and ROGeometry
211    delete detectorROGeometry; // this should be safe in C++ also if we have a NULL pointer
212    //if (detectorSD) detectorSD->PrintAll();
213    //delete detectorSD;
214    // Sensitive Detector and ReadOut geometry definition
215    G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
216
217    static G4String sensitiveDetectorName = "Detector"; 
218    if (!detectorSD)
219        {
220            // The sensitive detector is instantiated
221            detectorSD = new HadrontherapyDetectorSD(sensitiveDetectorName);
222        }
223    // The Read Out Geometry is instantiated
224    static G4String ROGeometryName = "DetectorROGeometry";
225    detectorROGeometry = new HadrontherapyDetectorROGeometry(ROGeometryName,
226                                                            detectorToWorldPosition,
227                                                            detectorSizeX/2,
228                                                            detectorSizeY/2,
229                                                            detectorSizeZ/2,
230                                                            numberOfVoxelsAlongX,
231                                                            numberOfVoxelsAlongY,
232                                                            numberOfVoxelsAlongZ);
233
234    G4cout << "Instantiating new Read Out Geometry \"" << ROGeometryName << "\""<< G4endl;
235    // This will invoke Build() HadrontherapyDetectorROGeometry virtual method
236    detectorROGeometry -> BuildROGeometry();
237    // Attach ROGeometry to SDetector
238    detectorSD -> SetROgeometry(detectorROGeometry);
239    //sensitiveDetectorManager -> Activate(sensitiveDetectorName, true);
240    if (!sensitiveDetectorManager -> FindSensitiveDetector(sensitiveDetectorName, false))
241        {
242            G4cout << "Registering new DetectorSD \"" << sensitiveDetectorName << "\""<< G4endl;
243            // Register user SD
244            sensitiveDetectorManager -> AddNewDetector(detectorSD);
245            // Attach SD to detector logical volume
246            detectorLogicalVolume -> SetSensitiveDetector(detectorSD);
247        }
248}
249void  HadrontherapyDetectorConstruction::ParametersCheck()
250{
251    // Check phantom/detector sizes & relative position
252    if (!IsInside(detectorSizeX, 
253                detectorSizeY, 
254                detectorSizeZ,
255                phantomSizeX,
256                phantomSizeY,
257                phantomSizeZ,
258                detectorToPhantomPosition
259                ))
260        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector is not fully inside Phantom!");
261
262    // Check Detector sizes respect to the voxel ones
263
264    if ( detectorSizeX < sizeOfVoxelAlongX) {
265        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector X size must be bigger or equal than that of Voxel X");
266    }
267    if ( detectorSizeY < sizeOfVoxelAlongY) {
268        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Y size must be bigger or equal than that of Voxel Y");
269    }
270    if ( detectorSizeZ < sizeOfVoxelAlongZ) {
271        G4Exception("Error at HadrontherapyDetectorConstruction::ParametersCheck(). Detector Z size must be bigger or equal than that of Voxel Z");
272    }
273
274}
275/////////////////
276// MESSENGERS //
277////////////////
278
279G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
280{
281
282    if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
283    {
284        phantomMaterial  = pMat;
285        detectorMaterial = pMat;
286        if (detectorLogicalVolume && phantomLogicalVolume) 
287        {
288            detectorLogicalVolume -> SetMaterial(pMat); 
289            phantomLogicalVolume ->  SetMaterial(pMat);
290
291            G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
292            G4RunManager::GetRunManager() -> GeometryHasBeenModified();
293            G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
294        }
295    }
296    else
297    {
298        G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
299            " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl; 
300        G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl; 
301        return false;
302    }
303
304    return true;
305}
306/////////////////////////////////////////////////////////////////////////////
307void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
308{
309    if (sizeX > 0.) phantomSizeX = sizeX;
310    if (sizeY > 0.) phantomSizeY = sizeY;
311    if (sizeZ > 0.) phantomSizeZ = sizeZ;
312}
313/////////////////////////////////////////////////////////////////////////////
314/////////////////////////////////////////////////////////////////////////////
315void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
316{
317    if (sizeX > 0.) {detectorSizeX = sizeX;}
318    if (sizeY > 0.) {detectorSizeY = sizeY;}
319    if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
320    SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
321}
322/////////////////////////////////////////////////////////////////////////////
323
324void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
325{
326    if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
327    if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
328    if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
329}
330void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
331{
332    phantomPosition = pos;
333}
334
335/////////////////////////////////////////////////////////////////////////////
336void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
337{
338    detectorToPhantomPosition = displ;
339}
340/////////////////////////////////////////////////////////////////////////////
341void HadrontherapyDetectorConstruction::UpdateGeometry()
342{
343    ParametersCheck();
344
345    //G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
346    G4GeometryManager::GetInstance() -> OpenGeometry();
347    if (phantom)
348    {
349            phantom -> SetXHalfLength(phantomSizeX/2);
350            phantom -> SetYHalfLength(phantomSizeY/2);
351            phantom -> SetZHalfLength(phantomSizeZ/2);
352            phantomPhysicalVolume -> SetTranslation(phantomPosition);
353    }
354    else   ConstructPhantom();
355
356    // Get the center of the detector
357    SetDetectorPosition();
358    if (detector)
359    {
360                        detector -> SetXHalfLength(detectorSizeX/2);
361                        detector -> SetYHalfLength(detectorSizeY/2);
362                        detector -> SetZHalfLength(detectorSizeZ/2);
363                        detectorPhysicalVolume -> SetTranslation(detectorPosition);
364    }
365    else    ConstructDetector();
366   
367// Round to nearest integer number of voxel
368        numberOfVoxelsAlongX = lrint(detectorSizeX / sizeOfVoxelAlongX);
369        sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
370
371        numberOfVoxelsAlongY = lrint(detectorSizeY / sizeOfVoxelAlongY);
372        sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
373
374        numberOfVoxelsAlongZ = lrint(detectorSizeZ / sizeOfVoxelAlongZ);
375        sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
376
377    //G4cout << "*************** DetectorToWorldPosition " << GetDetectorToWorldPosition()/cm << "\n";
378    ConstructSensitiveDetector(GetDetectorToWorldPosition());
379
380    volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
381    massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
382    //  This will clear the existing matrix (together with all data inside it)!
383    matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX, 
384                                              numberOfVoxelsAlongY,
385                                              numberOfVoxelsAlongZ,
386                                              massOfVoxel);
387
388    // Inform the kernel about the new geometry
389    G4RunManager::GetRunManager() -> GeometryHasBeenModified();
390    G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
391
392    PrintParameters();
393}
394
395void HadrontherapyDetectorConstruction::PrintParameters()
396{
397
398    G4cout << "The (X,Y,Z) dimensions of the phantom are : (" << 
399        G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' << 
400        G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' << 
401        G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl; 
402   
403    G4cout << "The (X,Y,Z) dimensions of the detector are : (" << 
404        G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' << 
405        G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' << 
406        G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl; 
407
408    G4cout << "Displacement between Phantom and World is: "; 
409    G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") << 
410        "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") << 
411        "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
412
413    G4cout << "The (X,Y,Z) sizes of the Voxels are: (" << 
414        G4BestUnit(sizeOfVoxelAlongX, "Length")  << ',' << 
415        G4BestUnit(sizeOfVoxelAlongY, "Length")  << ',' << 
416        G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
417
418    G4cout << "The number of Voxels along (X,Y,Z) is: (" << 
419        numberOfVoxelsAlongX  << ',' <<
420        numberOfVoxelsAlongY  <<','  <<
421        numberOfVoxelsAlongZ  << ')' << G4endl;
422
423}
Note: See TracBrowser for help on using the repository browser.