Context Navigation

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

Last change on this file since 1272 was 1230, checked in by garnier, 15 years ago
update to geant4.9.3
File size: 15.7 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	#include "G4SDManager.hh"
31	#include "G4RunManager.hh"
32	#include "G4Box.hh"
33	#include "G4LogicalVolume.hh"
34	#include "G4ThreeVector.hh"
35	#include "G4PVPlacement.hh"
36	#include "globals.hh"
37	#include "G4Transform3D.hh"
38	#include "G4RotationMatrix.hh"
39	#include "G4Colour.hh"
40	#include "G4UserLimits.hh"
41	#include "G4UnitsTable.hh"
42	#include "G4VisAttributes.hh"
43	#include "G4NistManager.hh"
44	#include "HadrontherapyDetectorROGeometry.hh"
45	#include "HadrontherapyDetectorMessenger.hh"
46	#include "HadrontherapyDetectorSD.hh"
47	#include "HadrontherapyDetectorConstruction.hh"
48	#include "HadrontherapyMatrix.hh"
49
50	/////////////////////////////////////////////////////////////////////////////
51	HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
52	: motherPhys(physicalTreatmentRoom),
53	detectorSD(0), detectorROGeometry(0), matrix(0),
54	phantomPhysicalVolume(0),
55	detectorLogicalVolume(0), detectorPhysicalVolume(0),
56	phantomSizeX(20.cm), phantomSizeY(20.cm), phantomSizeZ(20.*cm), // Default half dimensions
57	detectorSizeX(2.cm), detectorSizeY(2.cm), detectorSizeZ(2.*cm),
58	phantomPosition(20.cm, 0.cm, 0.*cm),
59	detectorToPhantomPosition(0.cm,18.cm,18.*cm)// Default displacement of the detector respect to the phantom
60	{
61	// NOTE! that the HadrontherapyDetectorConstruction class
62	// does NOT inherit from G4VUserDetectorConstruction G4 class
63	// So the Construct() mandatory virtual method is inside another geometric class
64	// (like the passiveProtonBeamLIne, ...)
65
66	// Messenger to change parameters of the phantom/detector geometry
67	detectorMessenger = new HadrontherapyDetectorMessenger(this);
68
69	// Default detector voxels size
70	// 200 slabs along the beam direction (X)
71	sizeOfVoxelAlongX = 200 *um;
72	sizeOfVoxelAlongY = 2 * detectorSizeY;
73	sizeOfVoxelAlongZ = 2 * detectorSizeZ;
74
75	// Calculate (and eventually set) detector position by displacement, phantom size and detector size
76	SetDetectorPosition();
77
78	// Build phantom and associated detector
79	ConstructPhantom();
80	ConstructDetector();
81	// Set number of the detector voxels along X Y and Z directions.
82	// This will construct also the sensitive detector, the ROGeometry
83	// and the matrix where the energy deposited is collected!
84	SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
85	}
86
87	/////////////////////////////////////////////////////////////////////////////
88	HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
89	{
90	delete detectorROGeometry;// This should be safe in C++ even if the argument is a NULL pointer
91	delete matrix;
92	delete detectorMessenger;
93	}
94
95	void HadrontherapyDetectorConstruction::ConstructPhantom()
96	{
97	//----------------------------------------
98	// Phantom:
99	// A box used to approximate tissues
100	//----------------------------------------
101
102	G4bool isotopes = false;
103	G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
104	phantom = new G4Box("Phantom",phantomSizeX, phantomSizeY, phantomSizeZ);
105	phantomLogicalVolume = new G4LogicalVolume(phantom,
106	waterNist,
107	"phantomLog", 0, 0, 0);
108
109	phantomPhysicalVolume = new G4PVPlacement(0,
110	phantomPosition,
111	"phantomPhys",
112	phantomLogicalVolume,
113	motherPhys,
114	false,
115	0);
116
117	// Visualisation attributes of the phantom
118	red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
119	red -> SetVisibility(true);
120	red -> SetForceSolid(true);
121	//red -> SetForceWireframe(true);
122	phantomLogicalVolume -> SetVisAttributes(red);
123	}
124
125	/////////////////////////////////////////////////////////////////////////////
126	void HadrontherapyDetectorConstruction::ConstructDetector()
127	{
128	//-----------
129	// Detector
130	//-----------
131	G4bool isotopes = false;
132	G4Material* waterNist = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER", isotopes);
133	detector = new G4Box("Detector",detectorSizeX,detectorSizeY,detectorSizeZ);
134	detectorLogicalVolume = new G4LogicalVolume(detector,
135	waterNist,
136	"DetectorLog",
137	0,0,0);
138	// Detector is attached by default to the phantom face directly exposed to the beam
139	detectorPhysicalVolume = new G4PVPlacement(0,
140	detectorPosition, // Setted by displacement
141	"DetectorPhys",
142	detectorLogicalVolume,
143	phantomPhysicalVolume,
144	false,0);
145
146	// Visualisation attributes of the detector
147	skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. , 235/255. ));
148	skyBlue -> SetVisibility(true);
149	skyBlue -> SetForceSolid(true);
150	//skyBlue -> SetForceWireframe(true);
151	detectorLogicalVolume -> SetVisAttributes(skyBlue);
152
153	}
154	/////////////////////////////////////////////////////////////////////////////
155	void HadrontherapyDetectorConstruction::ConstructSensitiveDetector(G4ThreeVector detectorToWorldPosition)
156	{
157	// Install new Sensitive Detector and ROGeometry
158	delete detectorROGeometry; // this should be safe in C++ also if we have a NULL pointer
159
160	// Sensitive Detector and ReadOut geometry definition
161	G4SDManager* sensitiveDetectorManager = G4SDManager::GetSDMpointer();
162
163	G4String sensitiveDetectorName = "Detector";
164	if (!detectorSD)
165	{
166	// The sensitive detector is instantiated
167	detectorSD = new HadrontherapyDetectorSD(sensitiveDetectorName);
168	}
169	// The Read Out Geometry is instantiated
170	G4String ROGeometryName = "DetectorROGeometry";
171	detectorROGeometry = new HadrontherapyDetectorROGeometry(ROGeometryName,
172	detectorToWorldPosition,
173	detectorSizeX,
174	detectorSizeY,
175	detectorSizeZ,
176	numberOfVoxelsAlongX,
177	numberOfVoxelsAlongY,
178	numberOfVoxelsAlongZ);
179
180	G4cout << "Instantiating new Read Out Geometry \"" << ROGeometryName << "\""<< G4endl;
181	// This will invoke Build() HadrontherapyDetectorROGeometry virtual method
182	detectorROGeometry -> BuildROGeometry();
183	// Attach ROGeometry to SDetector
184	detectorSD -> SetROgeometry(detectorROGeometry);
185	//sensitiveDetectorManager -> Activate(sensitiveDetectorName, true);
186	if (!sensitiveDetectorManager -> FindSensitiveDetector(sensitiveDetectorName, false))
187	{
188	G4cout << "Registering new DetectorSD \"" << sensitiveDetectorName << "\""<< G4endl;
189	// Register user SD
190	sensitiveDetectorManager -> AddNewDetector(detectorSD);
191	// Attach SD to detector logical volume
192	detectorLogicalVolume -> SetSensitiveDetector(detectorSD);
193	}
194	}
195
196	/////////////////
197	// MESSENGERS //
198	////////////////
199	G4bool HadrontherapyDetectorConstruction::SetNumberOfVoxelBySize(G4double sizeX, G4double sizeY, G4double sizeZ)
200	{
201	// Only change positive dimensions
202	// XXX numberOfVoxels must be an integer, warn the user
203
204	if (sizeX > 0)
205	{
206	if (sizeX > 2*detectorSizeX)
207	{
208	G4cout << "WARNING: Voxel X size must be smaller or equal than that of detector X" << G4endl;
209	return false;
210	}
211	// Round to the nearest integer
212	numberOfVoxelsAlongX = lrint(2 * detectorSizeX / sizeX);
213	sizeOfVoxelAlongX = (2 * detectorSizeX / numberOfVoxelsAlongX );
214	if(sizeOfVoxelAlongX!=sizeX) G4cout << "Rounding " <<
215	G4BestUnit(sizeX, "Length") << " to " <<
216	G4BestUnit(sizeOfVoxelAlongX, "Length") << G4endl;
217	}
218
219	if (sizeY > 0)
220	{
221	if (sizeY > 2*detectorSizeY)
222	{
223	G4cout << "WARNING: Voxel Y size must be smaller or equal than that of detector Y" << G4endl;
224	return false;
225	}
226	numberOfVoxelsAlongY = lrint(2 * detectorSizeY / sizeY);
227	sizeOfVoxelAlongY = (2 * detectorSizeY / numberOfVoxelsAlongY );
228	if(sizeOfVoxelAlongY!=sizeY) G4cout << "Rounding " <<
229	G4BestUnit(sizeY, "Length") << " to " <<
230	G4BestUnit(sizeOfVoxelAlongY, "Length") << G4endl;
231	}
232	if (sizeZ > 0)
233	{
234	if (sizeZ > 2*detectorSizeZ)
235	{
236	G4cout << "WARNING: Voxel Z size must be smaller or equal than that of detector Z" << G4endl;
237	return false;
238	}
239	numberOfVoxelsAlongZ = lrint(2 * detectorSizeZ / sizeZ);
240	sizeOfVoxelAlongZ = (2 * detectorSizeZ / numberOfVoxelsAlongZ );
241	if(sizeOfVoxelAlongZ!=sizeZ) G4cout << "Rounding " <<
242	G4BestUnit(sizeZ, "Length") << " to " <<
243	G4BestUnit(sizeOfVoxelAlongZ, "Length") << G4endl;
244	}
245
246	G4cout << "The (X, Y, Z) sizes of the Voxels are: (" <<
247	G4BestUnit(sizeOfVoxelAlongX, "Length") << ", " <<
248	G4BestUnit(sizeOfVoxelAlongY, "Length") << ", " <<
249	G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
250
251	G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
252	numberOfVoxelsAlongX << ", " <<
253	numberOfVoxelsAlongY << ", " <<
254	numberOfVoxelsAlongZ << ')' << G4endl;
255
256	// This will clear the existing matrix (together with data inside it)!
257	matrix = HadrontherapyMatrix::getInstance(numberOfVoxelsAlongX,
258	numberOfVoxelsAlongY,
259	numberOfVoxelsAlongZ);
260
261	// Here construct the Sensitive Detector and Read Out Geometry
262	ConstructSensitiveDetector(GetDetectorToWorldPosition());
263	G4RunManager::GetRunManager() -> GeometryHasBeenModified();
264	return true;
265	}
266	/////////////////////////////////////////////////////////////////////////////
267	G4bool HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
268	{
269	// Check that the detector stay inside the phantom
270	if (sizeX > 0 && sizeX < sizeOfVoxelAlongX) {G4cout << "WARNING: Detector X size must be bigger than that of Voxel X" << G4endl; return false;}
271	if (sizeY > 0 && sizeY < sizeOfVoxelAlongY) {G4cout << "WARNING: Detector Y size must be bigger than that of Voxel Y" << G4endl; return false;}
272	if (sizeZ > 0 && sizeZ < sizeOfVoxelAlongZ) {G4cout << "WARNING: Detector Z size must be bigger than that of Voxel Z" << G4endl; return false;}
273
274	if (!IsInside(sizeX/2,
275	sizeY/2,
276	sizeZ/2,
277	phantomSizeX,
278	phantomSizeY,
279	phantomSizeZ,
280	detectorToPhantomPosition))
281	{return false;}
282	// Negative or null values mean don't change it!
283	if (sizeX > 0) {
284	detectorSizeX = sizeX/2;
285	detector -> SetXHalfLength(detectorSizeX);
286	}
287
288	if (sizeY > 0) {
289	detectorSizeY = sizeY/2;
290	detector -> SetYHalfLength(detectorSizeY);
291	}
292
293	if (sizeZ > 0) {
294	detectorSizeZ = sizeZ/2;
295	detector -> SetZHalfLength(detectorSizeZ);
296	}
297
298
299	G4cout << "The (X, Y, Z) dimensions of the detector are : (" <<
300	G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ", " <<
301	G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ", " <<
302	G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
303	// Adjust detector position
304	SetDetectorPosition();
305	// Adjust voxels number accordingly to new detector geometry
306	// Matrix will be re-instantiated!
307	// Voxels and ROGeometry must follow the detector!
308	SetNumberOfVoxelBySize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
309	G4RunManager::GetRunManager() -> GeometryHasBeenModified();
310	return true;
311	}
312
313	/////////////////////////////////////////////////////////////////////////////
314	G4bool HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
315	{
316
317	if (!IsInside(detectorSizeX,
318	detectorSizeY,
319	detectorSizeZ,
320	sizeX/2,//method parameters
321	sizeY/2,
322	sizeZ/2,
323	detectorToPhantomPosition
324	))
325	return false;
326
327	// Only change positive dimensions
328	if (sizeX > 0) {
329	phantomSizeX = sizeX/2;
330	phantom -> SetXHalfLength(phantomSizeX);
331	}
332	if (sizeY > 0) {
333	phantomSizeY = sizeY/2;
334	phantom -> SetYHalfLength(phantomSizeY);
335	}
336
337	if (sizeZ > 0) {
338	phantomSizeZ = sizeZ/2;
339	phantom -> SetZHalfLength(phantomSizeZ);
340	}
341
342
343	G4cout << "The (X, Y, Z) dimensions of the phantom are : (" <<
344	G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ", " <<
345	G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ", " <<
346	G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
347	//G4cout << '\n' << "Coordinate volume: " << phantomPhysicalVolume -> GetTranslation() << G4endl;
348	// Adjust detector position inside phantom
349	SetDetectorPosition();
350
351	ConstructSensitiveDetector(GetDetectorToWorldPosition());
352	G4RunManager::GetRunManager() -> GeometryHasBeenModified();
353	return true;
354	}
355	/////////////////////////////////////////////////////////////////////////////
356	G4bool HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector displacement)
357	{
358	// Set Phantom position respect to the World
359	// TODO check for overlap!
360	phantomPosition = displacement;
361	if (phantomPhysicalVolume)
362	{
363	phantomPhysicalVolume -> SetTranslation(phantomPosition);
364	G4cout << "Displacement between Phantom and World is: ";
365	G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") << ", " <<
366	"DY= "<< G4BestUnit(phantomPosition.getY(),"Length") << ", " <<
367	"DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
368
369	// Redraw ROGeometry!
370	ConstructSensitiveDetector(GetDetectorToWorldPosition());
371	G4RunManager::GetRunManager() -> GeometryHasBeenModified();
372	}
373	return true;
374	}
375
376	/////////////////////////////////////////////////////////////////////////////
377	G4bool HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displacement)
378	{
379	// Ignore negative values
380	if (displacement[0] < 0.) displacement[0] = detectorToPhantomPosition[0];
381	if (displacement[1] < 0.) displacement[1] = detectorToPhantomPosition[1];
382	if (displacement[2] < 0.) displacement[2] = detectorToPhantomPosition[2];
383
384	if (!IsInside(detectorSizeX,
385	detectorSizeY,
386	detectorSizeZ,
387	phantomSizeX,
388	phantomSizeY,
389	phantomSizeZ,
390	displacement // method parameter!
391	))
392	{return false;}
393	detectorToPhantomPosition = displacement;
394
395	// Adjust detector position inside phantom
396	SetDetectorPosition();
397
398	ConstructSensitiveDetector(GetDetectorToWorldPosition());
399	G4RunManager::GetRunManager() -> GeometryHasBeenModified();
400	return true;
401	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: