| 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 //GEANT4 - Depth-of-Interaction enabled Positron emission tomography (PET) advanced example
27
28 //Authors and contributors
29
30 // Author list to be updated, with names of co-authors and contributors from National Institute of Radiological Sciences (NIRS)
31
32 // Abdella M. Ahmed (1, 2), Andrew Chacon (1, 2), Harley Rutherford (1, 2),
33 // Hideaki Tashima (3), Go Akamatsu (3), Akram Mohammadi (3), Eiji Yoshida (3), Taiga Yamaya (3)
34 // Susanna Guatelli (2), and Mitra Safavi-Naeini (1, 2)
35
36 // (1) Australian Nuclear Science and Technology Organisation, Australia
37 // (2) University of Wollongong, Australia
38 // (3) National Institute of Radiological Sciences, Japan
39
40
41
42 //A whole-body PET scanner with depth-of-intercation (DOI) detector is
43 //defined in the doiPETDetectorConstruction class and varius types phantoms
44 //are also defined based on the NEMA NU-2 2012 standard protocol for
45 //performance evaluation PET scanners.
46 //The scanner's geometrical specifications (like number of crystals, number of rings ...) are defined as a global variable in the doiPETGlobalParameters.hh
47 //If the phantom is other there the phantoms defined, exception will be thrown. If you want to run the simulation without phantom (for any reason),
48 //you can comment out ConstructPhantom(world_logicalV) method.
49
50 #include "doiPETDetectorConstruction.hh"
51 #include "doiPETAnalysis.hh"
52 #include "doiPETDetectorConstructionMessenger.hh"
53 #include "doiPETGlobalParameters.hh"
54
55 #include "G4NistManager.hh"
56 #include "G4Box.hh"
57 #include "G4Tubs.hh"
58 #include "G4LogicalVolume.hh"
59 #include "G4PVPlacement.hh"
60 #include "G4RotationMatrix.hh"
61 #include "G4Transform3D.hh"
62 #include "G4SDManager.hh"
63 #include "G4MultiFunctionalDetector.hh"
64 #include "G4VPrimitiveScorer.hh"
65 #include "G4PSEnergyDeposit.hh"
66 #include "G4PSDoseDeposit.hh"
67 #include "G4VisAttributes.hh"
68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh"
70
71 //
72 #include "G4SubtractionSolid.hh"
73 #include "G4UnionSolid.hh"
74 #include "G4IntersectionSolid.hh"
75 //
76 #include "G4Element.hh"
77 #include "G4Material.hh"
78 #include "G4Sphere.hh"
79 #include "G4Colour.hh"
80 #include "G4RunManager.hh"
81 #include "G4StateManager.hh"
82 #include "G4GeometryManager.hh"
83 #include "G4PhysicalVolumeStore.hh"
84 #include "G4LogicalVolumeStore.hh"
85 #include "G4SolidStore.hh"
86
87
88 ////////////////////////////////////////////////////////////////////////////////////////
89
90 doiPETDetectorConstruction::doiPETDetectorConstruction()
91 : G4VUserDetectorConstruction(),
92 fCheckOverlaps(true)
93 {
94 DefineMaterials();
95 fDetectorMessenger = new doiPETDetectorConstructionMessenger(this);
96
97 //Initialize variables. The following variables can be changed via the run.mac file.
98 PhantomType = "NEMA_imageQualityPhantom_wholeBody";//"Phantom_sensitivity";// "NEMA_phantom_NECR";//
99 phantomRadius = 100 * mm;
100 phantomLength = 700 * mm;
101 phantomPosition = G4ThreeVector(0,0,0);
102 numOfSleeves = 5;
103 }
104
105 //////////////////////////////////////////////////////////////////////////////////////
106
107 doiPETDetectorConstruction::~doiPETDetectorConstruction()
108 {
109 delete fDetectorMessenger;
110 }
111
112 /////////////////////////////////// Define Materials /////////////////////////////////////////////////////
113
114 void doiPETDetectorConstruction::DefineMaterials()
115 {
116 G4NistManager* nist = G4NistManager::Instance();
117
118 //Define air
119 air = nist->FindOrBuildMaterial("G4_AIR");
120
121 //Define PMMA
122 pmma = nist->FindOrBuildMaterial("G4_PLEXIGLASS"); //Default 1.19 g/cm3
123
124 //Define water
125 water = nist->FindOrBuildMaterial("G4_WATER");
126
127 //Defining polyethylene from NIST and modifying the density
128 polyethylene = nist->BuildMaterialWithNewDensity("polyethylene","G4_POLYETHYLENE",0.959*g/cm3);
129 polyethylene->GetIonisation()->SetMeanExcitationEnergy(56*eV);
130
131 //polyethylene_NEMA, defualt density 0.94 g/cm3, and excitation energy = 57.4 eV
132 polyethylene_NEMA = nist->FindOrBuildMaterial("G4_POLYETHYLENE");
133
134
135 //Define expanded polystyrene by modifiying the density to mimic lung phantom used in phantom experiment
136 polystyrene = nist->BuildMaterialWithNewDensity( "polystyrene","G4_POLYSTYRENE",0.3*g/cm3);
137
138 isotopes = false;
139
140 //Defile Aluminum material for the detetor cover
141 Aluminum = nist->FindOrBuildMaterial("G4_Al", isotopes);
142
143 //Define elements for the GSO crystal (scintillator)
144 O = nist->FindOrBuildElement("O" , isotopes);
145 Si = nist->FindOrBuildElement("Si", isotopes);
146 Gd = nist->FindOrBuildElement("Gd", isotopes);
147
148
149 //define GSO crystal for PET detector
150 GSO = new G4Material("GSO", 6.7*g/cm3, 3);
151 GSO->AddElement(Gd, 2);
152 GSO->AddElement(Si, 1);
153 GSO->AddElement(O, 5);
154 crystalMaterial = nist->FindOrBuildMaterial("GSO");
155 }
156
157 ////////////////////////////////////////////////////////////////////////////////////////
158
159 G4VPhysicalVolume* doiPETDetectorConstruction::Construct()
160 {
161 //size of the world
162 worldSizeX = 2 * m;//0.5 *mm
163 worldSizeY = 2 * m;//0.5*mm
164 worldSizeZ = 4. * m;
165
166
167 // Define a solid shape to describe the world volume
168 G4Box* solid_world = new G4Box("world", worldSizeX/2., worldSizeY/2., worldSizeZ/2.);
169
170 // Define a logical volume for the world volume
171 world_logicalV = new G4LogicalVolume(solid_world, air, "world_logicalV", 0,0,0);
172
173
174 // Define the physical world volume
175 world_physicalV = new G4PVPlacement(0,G4ThreeVector(),world_logicalV, "world_physicalV", 0, false, 0, fCheckOverlaps);
176 world_logicalV->SetVisAttributes (G4VisAttributes::GetInvisible());
177
178 //NOTE!!!
179 //The scanner specification (like size and number of crystals in each detctor) are given in the "doiPETGlobalParameters.hh" header file.
180
181 //Each block detector is identified with its unique number, blockIndex.
182 blockIndex = 0;
183
184 //Each crystal is identified with its unique number, crystalIndex
185 crystalIndex = 0;
186
187 //This is to place the phantom.
188
189 ConstructPhantom(world_logicalV);
190
191 //Define air volume (box) to fill the detector block. Crystal elements (scintillators) is then placed.
192 sizeOfAirBox_DOI = (numberOfCrystal_DOI * sizeOfCrystal_DOI) + (numberOfCrystal_DOI - 1)*crystalGap_DOI;
193 sizeOfAirBox_axial = (numberOfCrystal_axial * sizeOfCrystal_axial) + (numberOfCrystal_axial - 1)*crystalGap_axial;
194 sizeOfAirBox_tangential = (numberOfCrystal_tangential * sizeOfCrystal_tangential) + (numberOfCrystal_tangential - 1)*crystalGap_tangential;
195
196 //This will pass the size of the detector block so that the position of the PMT will be calculated with respect to the axis of the detector block.
197 //see doiPETAnalysis.cc
198 pAnalysis = doiPETAnalysis::GetInstance();
199 pAnalysis->GetSizeOfDetector(sizeOfAirBox_DOI,sizeOfAirBox_tangential, sizeOfAirBox_axial);
200
201 G4cout<<"size of crytal element: "<<sizeOfCrystal_tangential<<" "<<sizeOfCrystal_axial<<" "<<sizeOfCrystal_DOI<<G4endl;
202 G4cout<<"Size of detector block (without Al cover): "<<sizeOfAirBox_tangential<<" "<<sizeOfAirBox_axial<<" "<<sizeOfAirBox_DOI<<G4endl;
203
204
205
206 //Define the size of the detector block.
207 sizeOfBlockDetector_DOI = sizeOfAirBox_DOI + AluminumCoverThickness;
208 sizeOfBlockDetector_axial = sizeOfAirBox_axial + AluminumCoverThickness;
209 sizeOfBlockDetector_tangential = sizeOfAirBox_tangential + AluminumCoverThickness;
210
211
212
213 //Define solid shape for the detector block
214 G4Box* blockDetector = new G4Box("blockDetector",sizeOfBlockDetector_DOI/2,sizeOfBlockDetector_tangential/2,sizeOfBlockDetector_axial/2);
215
216 //Define the logical volume for the detector block
217 blockDetector_logicalV = new G4LogicalVolume(blockDetector,Aluminum,"blockDetector_logicalV", 0,0,0);
218
219
220
221 //Define air (box) inside the detector block. Crystal elements will be placed in it.
222 G4Box* airBox = new G4Box("airBox", sizeOfAirBox_DOI/2, sizeOfAirBox_tangential/2,sizeOfAirBox_axial/2);
223
224 //Define the logical volume
225 airBox_logicalV = new G4LogicalVolume(airBox,air,"airBox_logicalV", 0,0,0);
226
227 //Define its physical volume and place it inside the detector block
228 airBox_physicalV = new G4PVPlacement (0,G4ThreeVector(0,0,0),airBox_logicalV,"airBox_physicalV", blockDetector_logicalV,false,0,fCheckOverlaps);
229
230
231 ///////////////////////////////////////// Arrange the PET ring and place the PET detectors in the ring(s) ////////////////////////////////////
232
233 for(G4int Ring = 0; Ring< numberOfRings; Ring++)
234 {
235 //place the detectors in a ring along the axial direction. Note that the ring gap between two adjcent rings is measured from scintillator to scintillator. It does not include the Aluminum thickness cover.
236
237 detectorPositionZ = (Ring-((G4double)numberOfRings)/2 + 0.5)*(sizeOfBlockDetector_axial + ringGap - AluminumCoverThickness);
238
239 for(G4int i = 0; i<numberOfDetector_perRing; i++)
240 {
241 //The azimuthal angle to arrange the detectors in a ring
242 thetaDetector = (double)(i*twopi/numberOfDetector_perRing);
243
244 //The radius of the scanner is measured from opposing crystal (scintillator) faces. It does not include the Aluminum thickness cover.
245 detectorPositionX = (scannerRadius + sizeOfBlockDetector_DOI/2 - AluminumCoverThickness/2)*std::cos(thetaDetector);
246 detectorPositionY = (scannerRadius + sizeOfBlockDetector_DOI/2 - AluminumCoverThickness/2)*std::sin(thetaDetector);
247
248 //Define the rotation matrix for correct placement of detetors
249 G4RotationMatrix rotm_PET = G4RotationMatrix();
250 rotm_PET.rotateZ(thetaDetector);
251 G4ThreeVector uz_PET = G4ThreeVector(detectorPositionX,detectorPositionY,detectorPositionZ);
252 G4Transform3D transform = G4Transform3D(rotm_PET,uz_PET);
253
254 //Define the physical volume of the detectors.
255 blockDetector_physicalV = new G4PVPlacement (transform,blockDetector_logicalV,"blockDetector_physicalV", world_logicalV,false,blockIndex,fCheckOverlaps);
256 blockIndex++;
257 //G4cout<<Ring<<" "<<detectorPositionX- ((sizeOfBlockDetector_DOI - AluminumCoverThickness)/2)*cos(thetaDetector)<<" "<<detectorPositionY- ((sizeOfBlockDetector_DOI- AluminumCoverThickness)/2)*sin(thetaDetector)<<" "<<detectorPositionZ<<G4endl;
258
259 }
260 }
261
262 //Define the solid crystal
263 G4VSolid* CrystalSolid = new G4Box("Crystal", sizeOfCrystal_DOI/2., sizeOfCrystal_tangential/2., sizeOfCrystal_axial/2.);
264
265 //Define the local volume of the crystal
266 crystal_logicalV = new G4LogicalVolume(CrystalSolid,crystalMaterial,"Crystal_logicalV", 0,0,0);
267
268 //Place the crystals inside the detectors and give them a unique number with crystalIndex.
269 for(G4int i_DOI = 0; i_DOI<numberOfCrystal_DOI; i_DOI++){
270 crystalPositionX=(i_DOI-((G4double)numberOfCrystal_DOI)/2 + 0.5)*(sizeOfCrystal_DOI + crystalGap_DOI);
271 for(G4int i_axial=0; i_axial< numberOfCrystal_axial;i_axial++){
272 crystalPositionZ = (i_axial-((G4double)numberOfCrystal_axial)/2 + 0.5)*(sizeOfCrystal_axial + crystalGap_axial);
273 for(G4int i_tan=0; i_tan<numberOfCrystal_tangential;i_tan++){
274 crystalPositionY=(i_tan-((G4double)numberOfCrystal_tangential)/2 + 0.5)*(sizeOfCrystal_tangential + crystalGap_tangential);
275
276 //G4cout<<crystalIndex<<" "<<crystalPositionX<<" "<<crystalPositionY<<" "<<crystalPositionZ<<G4endl;
277 //place the crystal inside the block detector.
278 crystal_physicalV = new G4PVPlacement (0, G4ThreeVector (crystalPositionX,crystalPositionY,crystalPositionZ), crystal_logicalV, "Crystal_physicalV", airBox_logicalV,false,crystalIndex/*,fCheckOverlaps*/);
279 crystalIndex++;
280 }
281 }
282 }
283
284 //****************** Visualization *****************************//
285
286 //visualization for the block detector
287 G4VisAttributes* blockDetectorVisAtt;
288 blockDetectorVisAtt = new G4VisAttributes(G4Colour(1,1.0,1.0));
289 blockDetectorVisAtt->SetVisibility (true);
290 //blockDetectorVisAtt->SetForceWireframe (true);
291 blockDetector_logicalV->SetVisAttributes (blockDetectorVisAtt);
292 //blockDetector_logicalV->SetVisAttributes (G4VisAttributes::GetInvisible());
293
294 //visualization for the the box filled with air
295 G4VisAttributes* airBoxVisAtt;
296 airBoxVisAtt = new G4VisAttributes(G4Colour(1,1.0,1.0));
297 airBoxVisAtt->SetVisibility (true);
298 airBoxVisAtt->SetForceWireframe (true);
299 airBox_logicalV->SetVisAttributes (airBoxVisAtt);
300 airBox_logicalV->SetVisAttributes (G4VisAttributes::GetInvisible());
301
302 //visualization for the crystal
303 G4VisAttributes* crystalVisAtt;
304 crystalVisAtt = new G4VisAttributes(G4Colour(0.88,0.55,1.0));
305 //crystalVisAtt->SetVisibility (true);
306 crystalVisAtt->SetForceWireframe (true);
307 crystal_logicalV->SetVisAttributes (crystalVisAtt);
308 crystal_logicalV->SetVisAttributes (G4VisAttributes::GetInvisible());
309
310
311 //always return the physical World
312 return world_physicalV;
313 }
314
315
316 /////////////////////////////////////////////// Construct Phantom //////////////////////////////////////////////
317
318 void doiPETDetectorConstruction::ConstructPhantom(G4LogicalVolume* worldLogical)
319 {
320 G4cout<<"------------------================- " <<PhantomType<<" Phantom has been used -==================-----------\n"<<G4endl;
321
322 //The following phantoms are defined based on NEMA NU-2 Standards Publication (Performance Measurements of Positron Emission Tomographs)
323 if(PhantomType == "NEMA_imageQualityPhantom_wholeBody"){
324
325 //The body phantom (also called image quality phantom) can presisely be created by the G4UninSolid method from (1) one big half cylinder (with a diameter of 150 mm),
326 //(2) two small full (or quarter) cylinders (with 80mm), and (3) a box with a height width of 140 mm and height of 80mm.
327 //All the material to surround (hold) the water is PMMA.
328
329 //offset position for the body phantom. The center of the body phantom is mover by 35 mm in the +y direction
330 yOffsetBodyPhantom = 35*mm;
331
332 //The body phantom is moved by 70 mm in the z direction so that the sphere (cold and hot spheres) are coplanar with the middle slice of the scanner
333 zOffsetBodyPhantom =70*mm;
334
335 //length (in the z-direction) of the body phantom.
336 lengthOfBodyPhantom = 180*mm; //Interior length ( = 180m mm) + wallthickness (= 3mm)
337
338 //outer radius of the body phantom. (inner radius + wallthinkness)
339 radiusOfBodyPhantom = 147 * mm;
340
341 //wall thickness of the surrounding pmma phantom
342 wallThicknessOfBodyPhantom = 3 * mm;
343
344 //
345 radiusOfSmallcyl = 77 * mm;
346 boxWidth = 140*mm;// along the x-axis
347 boxHeight = 77 * mm; // along the y-axis
348
349
350 /************************* Surrounding PMMA ***********************************/
351
352 //define a half cylinder
353 G4Tubs* halfcyl_pmma = new G4Tubs("halfcyl_pmma", 0*mm,(radiusOfBodyPhantom + wallThicknessOfBodyPhantom),(lengthOfBodyPhantom + wallThicknessOfBodyPhantom)/2,0*deg,180*deg);
354
355 //define two full small cylinder. (It can be also be a quarter or half).
356 G4Tubs* cyl1_pmma = new G4Tubs("cyl_pmma", 0*mm,radiusOfSmallcyl + wallThicknessOfBodyPhantom, (lengthOfBodyPhantom + wallThicknessOfBodyPhantom)/2,0*deg,360*deg);
357 G4Tubs* cyl2_pmma = new G4Tubs("cyl_pmma", 0*mm,radiusOfSmallcyl + wallThicknessOfBodyPhantom, (lengthOfBodyPhantom + wallThicknessOfBodyPhantom)/2,0*deg,360*deg);
358
359 //define a box
360 G4Box* box_pmma = new G4Box("box_pmma", boxWidth/2, (boxHeight + wallThicknessOfBodyPhantom)/2, (lengthOfBodyPhantom + wallThicknessOfBodyPhantom)/2);
361
362 //a translation vector for the small cylinder with respect to the half cylinder at position (70 mm,0,0)
363 G4ThreeVector translation_cyl1 = G4ThreeVector(boxWidth/2, 0, 0);
364
365 //a translation vector for the small cylinder with respect to the half cylinder at position (-70 mm,0,0)
366 G4ThreeVector translation_cyl2 = G4ThreeVector(-boxWidth/2, 0, 0);
367
368 //translation for the box with respect to half cylinder
369 G4ThreeVector translation_box = G4ThreeVector(0,-(boxHeight + wallThicknessOfBodyPhantom)/2, 0);
370
371 //define union1_pmma by uniting the solids of halfcyl and cyl1
372 G4UnionSolid* union1_pmma = new G4UnionSolid("union1",halfcyl_pmma,cyl1_pmma,0,translation_cyl1);
373
374 //define union2_pmma by uniting the solids of union1_pmma and cyl2
375 G4UnionSolid* union2_pmma = new G4UnionSolid("union2",union1_pmma,cyl2_pmma,0,translation_cyl2);
376
377 //********** Now define the solid volume of the surrounding PMMA body phantom *********//
378 //define phantom by uniting the solids of union2_pmma and box_pmma
379 G4UnionSolid* phantom = new G4UnionSolid("phantom",union2_pmma,box_pmma,0,translation_box);
380
381 //Define the logical volume of the body phantom
382 phantom_logicalV = new G4LogicalVolume(phantom, pmma, "phantom_logicalV",0,0,0);
383
384 //Define the physical volume of the body phantom
385 phantom_physicalV = new G4PVPlacement(0, G4ThreeVector(0,-yOffsetBodyPhantom,-(lengthOfBodyPhantom + 2*wallThicknessOfBodyPhantom)/2 + zOffsetBodyPhantom), phantom_logicalV, "pmma_phantom_physicalV", worldLogical, false, 0, fCheckOverlaps);
386
387
388
389 //****************************** Water inside the PMMA phantom ********************************/
390
391 //define a half cylinder
392 G4Tubs* halfcyl_water = new G4Tubs("halfcyl_water", 0*mm,radiusOfBodyPhantom,lengthOfBodyPhantom/2,0*deg,180*deg);
393
394 //define a full small cylinder phantom. (It can be also be a quarter or half).
395 G4Tubs* cyl1_water = new G4Tubs("cyl_water", 0*mm,radiusOfSmallcyl, lengthOfBodyPhantom/2,0*deg,360*deg);
396 G4Tubs* cyl2_water = new G4Tubs("cyl_water", 0*mm,radiusOfSmallcyl, lengthOfBodyPhantom/2,0*deg,360*deg);
397
398 //define a box
399 G4Box* box_water = new G4Box("box_water", boxWidth/2, boxHeight/2, lengthOfBodyPhantom/2);
400
401 //a translation vector for the small cylinder with respect to the half cylinder at position (70 mm,0,0)
402 G4ThreeVector translation_cyl1_water = G4ThreeVector(boxWidth/2, 0, 0);
403
404 //a translation vector for the small cylinder with respect to the half cylinder at position (-70 mm,0,0)
405 G4ThreeVector translation_cyl2_water = G4ThreeVector(-boxWidth/2, 0, 0);
406
407 //translation for the box with respect to half cylinder
408 G4ThreeVector translation_box_water = G4ThreeVector(0,-boxHeight/2, 0);
409
410 //define union1_water by uniting the solids of halfcyl and cyl1
411 G4UnionSolid* union1_water = new G4UnionSolid("union1",halfcyl_water,cyl1_water,0,translation_cyl1_water);
412
413 //define union2_water by uniting the solids of union1_water and cyl2
414 G4UnionSolid* union2_water = new G4UnionSolid("union2",union1_water,cyl2_water,0,translation_cyl2_water);
415
416 //********** Now define the solid volume of the body phantom to be filled with water *********//
417 //define phantom by uniting the solids of union2_water and box_water
418 G4UnionSolid* phantom_water = new G4UnionSolid("phantom_water",union2_water,box_water,0,translation_box_water);
419
420 //Define the logical volume of the body phantom
421 water_logicalV = new G4LogicalVolume(phantom_water, water, "phantom_logicalV",0,0,0);
422
423 //Define the physical volume of the body phantom
424 water_physicalV = new G4PVPlacement(0, G4ThreeVector(0,0,0), water_logicalV, "phantom_physicalV", phantom_logicalV, false, 0, fCheckOverlaps);
425
426 //
427 /////////////////////////////// lung phantom /////////////////////////////////////////////////////
428
429 //A lung phantom (with low density material) is inserted in the body phantom.
430 G4double wallThiknessOfLungPhantom = 4 *mm;
431 radiusOfLungPhantom = 21 *mm;
432
433 //define surrounding pmma phantom for the lung
434 //Define the solid shape for the lung phantom
435 G4Tubs* phantom_lungPMMA = new G4Tubs("Phantom_lung", 0*mm, radiusOfLungPhantom + wallThiknessOfLungPhantom,lengthOfBodyPhantom/2,0*deg,360*deg);
436
437 //Define the logical volume of for the lung phantom
438 lung_logicalV_PMMA = new G4LogicalVolume(phantom_lungPMMA, pmma, "coldRegion_logicalV",0,0,0);
439
440 //Define the physical volume for the lung phantom. The center of the lung phantom is moved (in the y-axis) by the same distance as that of the body phantom.
441 lung_physicalVPMMA = new G4PVPlacement(0, G4ThreeVector(0,yOffsetBodyPhantom,0), lung_logicalV_PMMA, "coldRegion_physicalV", water_logicalV, false, 0, fCheckOverlaps);
442
443
444 //Define the solid shape for the lung phantom
445 G4Tubs* phantom_lung = new G4Tubs("Phantom_lung", 0*mm,radiusOfLungPhantom,lengthOfBodyPhantom/2,0*deg,360*deg);
446
447 //Define the logical volume of for the lung phantom
448 lung_logicalV = new G4LogicalVolume(phantom_lung, polystyrene, "coldRegion_logicalV",0,0,0);
449
450 //Define the physical volume for the lung phantom and place it in the phantom_lungPMMA.
451 lung_physicalV = new G4PVPlacement(0, G4ThreeVector(0,0,0), lung_logicalV, "coldRegion_physicalV", lung_logicalV_PMMA, false, 0, fCheckOverlaps);
452 //lung_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
453
454
455 ////////////////////////////////// Test phantom ////////////////////////////////////////////////////////
456 //This phantom has the same characteristics with that of NECR phantom except its axial position.
457
458 hieghtOfTestPhantom = 705* mm;
459 diameterOfTestPhantom = 203 * mm;
460 G4double zOffset_testPhantom = 0*mm;//this is to make some gap between the phantoms if necessary
461
462 //Define the solid shape for the test phantom
463 G4Tubs* phantom_test = new G4Tubs("Phantom", 0*mm,diameterOfTestPhantom/2,hieghtOfTestPhantom/2,0*deg,360*deg);
464
465 //Define the logical volume of the test phantom
466 test_logicalV = new G4LogicalVolume(phantom_test, polyethylene_NEMA, "phantom_logicalV",0,0,0);
467
468 //Define the physical volume of the test phantom. The test phantom is placed next to the body phantom.
469 test_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,hieghtOfTestPhantom/2+zOffsetBodyPhantom + zOffset_testPhantom), test_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
470 //test_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
471
472 ////////////////////////////////// Six Spherical phantoms for hot and cold lesions (sources) placed inside the body phatom ///////////////////////////////////
473
474 //Define diameter from the center of the body phantom to the center of the spheres
475 distanceFromCenter = 114.4*mm;
476
477 //Define total number of spherical phantoms (both hot and cold spheres)
478 numberOfSpheres = 6;
479
480 //Define the wall thickness of the spherical phantoms. It should be less than or equal to 1 mm according to the NEMA NU-2 protocol.
481 sphereWallThickness = 1*mm;
482
483 //The centers of the sphere phantoms are placed 68 mm from the endplate of the body phantom so that the plane through the centers of the spheres is coplanar with to the middile slice of the scanner
484 zOffsetSpherePhantom = lengthOfBodyPhantom/2 + wallThicknessOfBodyPhantom- zOffsetBodyPhantom;
485
486 //Place the spherical phantoms in the body phantom
487 for(G4int i = 0; i < numberOfSpheres; i++){
488 if(i == 0) sphereDiameter = 37*mm;//cold sphere
489 if(i == 1) sphereDiameter = 10*mm;// hot sphere
490 if(i == 2) sphereDiameter = 13*mm;// hot sphere
491 if(i == 3) sphereDiameter = 17*mm;// hot sphere
492 if(i == 4) sphereDiameter = 22*mm;// hot sphere
493 if(i == 5) sphereDiameter = 28*mm;// cold sphere
494
495 spherePositionX = distanceFromCenter/2*std::cos(((G4double)i*twopi/numberOfSpheres));
496 spherePositionY = distanceFromCenter/2*std::sin(((G4double)i*twopi/numberOfSpheres));
497 //zOffsetBodyPhantom = 0;
498
499 //hot sphere phantoms
500 if(i>0 && i <5){
501
502 //Surrounding PMMA phantom for the hot sphere
503 //Define the solid shape for the surrounding PMMA phantom for the hot sphere pahntom
504 G4Sphere* hotSpherePMMA = new G4Sphere("hotSphere", 0., sphereDiameter/2 + sphereWallThickness, 0*deg, 360*deg, 0*deg, 180*deg);
505
506 //Deifne the logical volume of the surrounding PMMA for the hot sphere phantom
507 hotSpherePMMA_logicalV = new G4LogicalVolume(hotSpherePMMA, pmma , "hotSphere_logicalV",0,0,0);
508
509 //Deifne the physical volume of the surrounding PMMA for the hot sphere phantom
510 hotSpherePMMA_physicalV = new G4PVPlacement(0, G4ThreeVector(spherePositionX,spherePositionY+yOffsetBodyPhantom,zOffsetSpherePhantom), hotSpherePMMA_logicalV, "hotSphere_physicalV", water_logicalV, false, i,fCheckOverlaps);
511
512
513 //Defining and placing the water in the surrounding PMMA for hot sphere
514 //Define the solid shape of the water phantom for the hot sphere phantom
515 G4Sphere* hotSphereWater = new G4Sphere("hotSphere", 0., sphereDiameter/2, 0*deg, 360*deg, 0*deg, 180*deg);
516
517 //Define the logical volume of the water phatom for the hot sphere
518 hotSphereWater_logicalV = new G4LogicalVolume(hotSphereWater, water , "hotSphere_logicalV",0,0,0);
519
520 //Define the physical volume of the water phatom for the hot sphere
521 hotSphereWater_physicalV = new G4PVPlacement(0, G4ThreeVector(0,0,0), hotSphereWater_logicalV, "phantom_physicalV", hotSpherePMMA_logicalV, false, i,fCheckOverlaps);
522 //G4cout<<"Hot sphere "<<i<<" is placed. "<< " sphereDiameter = " <<sphereDiameter <<" "<< "Position " <<spherePositionX<<" "<< spherePositionY<<", "<<sphereDiameter<<G4endl;
523 //hotSpherePMMA_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
524 //hotSphereWater_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
525 }
526
527 //Cold sphere phantoms
528 if(i==0 || i==5){
529 //Surrounding PMMA phantom for the cold sphere
530 //Define the solid shape for the surrounding PMMA phantom for the cold sphere pahntom
531 G4Sphere* coldSpherePMMA = new G4Sphere("coldSphere", 0., sphereDiameter/2 + sphereWallThickness, 0*deg, 360*deg, 0*deg, 180*deg);
532
533 //Deifne the logical volume of the surrounding PMMA for the cold sphere phantom
534 coldSpherePMMA_logicalV = new G4LogicalVolume(coldSpherePMMA, pmma , "coldRegion_logicalV",0,0,0);
535
536 //Deifne the physical volume of the surrounding PMMA for the cold sphere phantom
537 coldSpherePMMA_physicalV = new G4PVPlacement(0, G4ThreeVector(spherePositionX,spherePositionY+yOffsetBodyPhantom,zOffsetSpherePhantom), coldSpherePMMA_logicalV, "coldRegion_physicalV", water_logicalV, false, i,fCheckOverlaps);
538
539
540 //Defining and placing the water in the surrounding PMMA for cold sphere
541 //Define the solid shape of the water phantom for the cold sphere phantom
542 G4Sphere* coldSphereWater = new G4Sphere("coldSphere", 0., sphereDiameter/2, 0*deg, 360*deg, 0*deg, 180*deg);
543
544 //Define the logical volume of the water phatom for the cold sphere
545 coldSphereWater_logicalV = new G4LogicalVolume(coldSphereWater, water , "coldRegion_logicalV",0,0,0);
546
547 //Define the physical volume of the water phatom for the cold sphere
548 coldSphereWater_physicalV = new G4PVPlacement(0, G4ThreeVector(0,0,0), coldSphereWater_logicalV, "coldRegion_physicalV", coldSpherePMMA_logicalV, false, i,fCheckOverlaps);
549 //G4cout<<"Cold sphere "<<i<<" is placed. "<< " sphereDiameter = " <<sphereDiameter <<" "<< "Position " <<spherePositionX<<" "<< spherePositionY<<" "<<sphereDiameter<<G4endl;
550 //coldSpherePMMA_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
551 //coldSphereWater_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
552 }
553 }
554 }
555
556 else if(PhantomType == "NEMA_imageQualityPhantom_smallAnimal"){
557 //The follwoing is for NEMA NU-2 image quality phantom for small animal
558 //To see the details of the phantom, please see: http://www.qrm.de/content/pdf/QRM-MicroPET-IQ.pdf
559
560 //Outside radius of PMMMA
561 phantomRadius = 16.75*mm;// dia=33.5*mm;
562
563 //Outside length of PMMA
564 phantomLength = 63*mm;
565
566 //Dimension of water phantom to be filled with activity (hot region).
567 waterPhantomRadius = 15*mm;
568 waterPhantomLength = 30*mm;
569
570
571 //Dimension of rod phantom (hot region)
572 rodPhantomLength = 20*mm;
573
574 //There are five rod phantoms with different diameters
575 numberOfRods = 5;
576
577 //Distance from the center of the cylinder to the center of the rods
578 distanceFromCenter = 7*mm;
579
580 //surrounding PMMA phantom.
581 //Define PMMA solid
582 G4Tubs* phantom = new G4Tubs("Phantom", 0*mm,phantomRadius,phantomLength/2,0*deg,360*deg);
583
584 //Define PMMA logical volume
585 phantom_logicalV = new G4LogicalVolume(phantom, pmma, "phantom_logicalV",0,0,0);
586
587 //Define PMMA physical volume
588 phantom_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,0), phantom_logicalV, "pmma_phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
589
590 //The rods are placed at one end of the surrounding PMMA cylinderical phantom
591 rodPositionZ = -waterPhantomLength/2;
592
593 for(G4int i = 0; i < numberOfRods; i++){
594
595 if(i == 0) rodDiameter = 1 * mm;
596 if(i == 1) rodDiameter = 2 * mm;
597 if(i == 2) rodDiameter = 3 * mm;
598 if(i == 3) rodDiameter = 4 * mm;
599 if(i == 4) rodDiameter = 5 * mm;
600
601 rodPositionX = distanceFromCenter*std::cos(((G4double)i*twopi/numberOfRods));
602 rodPositionY = distanceFromCenter*std::sin(((G4double)i*twopi/numberOfRods));
603
604 //Define rod phantom
605 G4Tubs* rod_phantom = new G4Tubs("Phantom", 0*mm,rodDiameter/2,rodPhantomLength/2,0*deg,360*deg);
606
607 //Define rod phantom logical volume
608 rod_phantom_logicalV = new G4LogicalVolume(rod_phantom, water, "rod_phantom_logicalV",0,0,0);
609
610 //Define rod phantom physical volume
611 rod_phantom_physicalV = new G4PVPlacement(0,G4ThreeVector(rodPositionX,rodPositionY,rodPositionZ), rod_phantom_logicalV, "phantom_physicalV", phantom_logicalV, false, 0, fCheckOverlaps);
612 }
613
614 //Out dimensions of the surrounding PMMA for cold region chambers
615 chamberPhantomLength = 15*mm;
616 chamberDiameter = 10*mm;
617
618 //Wall thickness of the surrounding PMMA phantom for the cold region chambers
619 wallThicknessOfChamber = 1*mm;
620
621 //The centers of the cold chambers is
622 chamberPositionX = distanceFromCenter;
623 chamberPositionY = 0*mm;
624 chamberPositionZ = waterPhantomLength/2 - chamberPhantomLength/2;
625
626 //hot region filled with water
627 G4Tubs* water_phantom = new G4Tubs("Phantom", 0*mm,waterPhantomRadius,waterPhantomLength/2,0*deg,360*deg);
628
629 waterPhantom_logicalV = new G4LogicalVolume(water_phantom, water, "waterPhantom_logicalV",0,0,0);
630
631 //place the phantom at one end of the PMMA phantom
632 WaterPhantom_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,rodPhantomLength/2), waterPhantom_logicalV, "phantom_physicalV", phantom_logicalV, false, 0,fCheckOverlaps);
633 //waterPhantom_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
634
635 //define the surrounding PMMA chamber
636 G4Tubs* chamberPMMA = new G4Tubs("chamber_phantom1", 0*mm,chamberDiameter/2,chamberPhantomLength/2,0*deg,360*deg);
637
638 //define the logical volume of the surrounding PMMA chamber
639 chamberPMMA_logicalV = new G4LogicalVolume(chamberPMMA, pmma, "chamberPMMA_logicalV",0,0,0);
640
641 //define the physical volume of the surrounding PMMA chamber
642 chamberPMMA_physicalV = new G4PVPlacement(0,G4ThreeVector(chamberPositionX, chamberPositionY, chamberPositionZ), chamberPMMA_logicalV, "pmma_phantom_physicalV", waterPhantom_logicalV, false, 0,fCheckOverlaps);
643
644
645 //Two cold region chambers: one of them filled with (non-radioactive) water and the other with air and is placed inside a hot region
646 //chamber one filled with water (cold region)
647 //Define the cold region chamber phantom to be filled with water
648 G4Tubs* chamberWater = new G4Tubs("chamber_phantom1", 0*mm,chamberDiameter/2 - wallThicknessOfChamber,chamberPhantomLength/2 - wallThicknessOfChamber,0*deg,360*deg);
649
650 //Define the logical volume of the cold region phantom
651 chamberWater_logicalV = new G4LogicalVolume(chamberWater, water, "chamber_phantom_logicalV",0,0,0);
652
653 //Define the physical volume of the cold region phantom
654 chamberWater_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,0), chamberWater_logicalV, "coldRegion_physicalV", chamberPMMA_logicalV, false, 0,fCheckOverlaps);
655
656 //chamber2 filled with air
657 //Place the surrounding chamber at another location (at -chamberPositionX)
658 chamberPMMA_physicalV = new G4PVPlacement(0,G4ThreeVector(-chamberPositionX, chamberPositionY, chamberPositionZ), chamberPMMA_logicalV, "pmma_phantom_physicalV", waterPhantom_logicalV, false, 0,fCheckOverlaps);
659
660 //Define the logical volume of the cold region phantom to be filled with air
661 G4Tubs* chamberAir = new G4Tubs("chamber_phantom1", 0*mm,chamberDiameter/2 - wallThicknessOfChamber,chamberPhantomLength/2 - wallThicknessOfChamber,0*deg,360*deg);
662
663 //Define its logical volume
664 chamberAir_logicalV = new G4LogicalVolume(chamberAir, air, "chamber_phantom_logicalV",0,0,0);
665
666 //Define the physical volume of air-filled chamber
667 chamberAir_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,0), chamberAir_logicalV, "coldRegion_physicalV", chamberPMMA_logicalV, false, 0,fCheckOverlaps);
668
669 }
670
671 else if(PhantomType == "NEMA_phantom_NECR"){
672 ////////////////////////////// Phantom for NECR //////////////////////////////////////
673 //The phantom is 203 mm in diameter and 700 mm in height and is placed at the center of the AFOV
674
675 //Define its solid shape
676 G4Tubs* phantom = new G4Tubs("Phantom", 0*mm,phantomRadius,phantomLength/2,0*deg,360*deg);
677
678 //Define its logical volume
679 phantom_logicalV = new G4LogicalVolume(phantom, polyethylene_NEMA, "phantom_logicalV",0,0,0);
680
681 //Define its physical volume
682 phantom_physicalV = new G4PVPlacement(0,G4ThreeVector(0,0,0), phantom_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
683 }
684 else if(PhantomType == "Phantom_sensitivity"){
685
686 ///////////////////////////////////////// Sensitiviy phantom /////////////////////////////////////////////
687 //There are six diffrent sizes of the sensitivity phantom which differe in their diameter.
688 //The size of the phantom can be changed via the macro file by changing the number of sleeves
689 //The length of sensitivity phantom is 700 mm
690 //(3.9, 6.4), (7.0, 9.5), (10.2, 12.7), (13.4, 15.9), (16.6, 19.1)
691
692 G4double Rmin, Rmax;
693 G4cout<<"Total number of sleeves : "<<numOfSleeves<<G4endl;
694
695 for(G4int i=0; i<numOfSleeves; i++){
696 if(i==0){
697 Rmin = 3.9/2 * mm;
698 Rmax = 3.2 * mm;
699 }
700 if(i==1){
701 Rmin = 7.0/2 * mm;
702 Rmax = 4.75 * mm;
703 }
704 if(i==2){
705 Rmin = 10.2/2 * mm;
706 Rmax = 6.35*mm;
707 }
708 if(i==3){
709 Rmin = 13.4/2 * mm;
710 Rmax = 7.95 *mm;
711 }
712 if(i==4){
713 Rmin = 16.6/2 * mm;
714 Rmax = 9.55 * mm;
715
716 }
717 G4cout<<"Sleeve "<<i+1 <<" is placed, Rmin = "<<Rmin<< " mm, Rmax = "<<Rmax<< " mm, "<< "length = " <<phantomLength <<" mm."<<G4endl;
718
719 //Concenric aluminum sleeves
720 G4Tubs* phantom = new G4Tubs("Phantom_sensitivity", Rmin,Rmax,phantomLength/2,0*deg,360*deg);
721
722 //Define its logical volume. The Matrial of the tubes are Aluminum
723 phantom_logicalV = new G4LogicalVolume(phantom, Aluminum, "phantom_logicalV",0,0,0);
724
725 //Define its physical volume
726 phantom_physicalV = new G4PVPlacement(0, phantomPosition, phantom_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
727 }
728
729 //Define the inner most polyethylene (PE) tube with a diameter of 3 mm. This size will not be changed. This is the phantom that holds the source.
730 G4Tubs* phantomPE = new G4Tubs("Phantom_sensitivity", 0 *mm,1.5*mm,phantomLength/2,0*deg,360*deg);
731
732 //Define its logical volume. The Matrial of the tubes are polyethylene
733 phantomPE_logicalV = new G4LogicalVolume(phantomPE, polyethylene_NEMA, "phantom_logicalV",0,0,0);
734
735 //Define its physical volume
736 phantomPE_physicalV = new G4PVPlacement(0, phantomPosition, phantomPE_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
737
738
739 G4cout<<"---------------Phantom dimension: "<<"radius (mm) = "<<phantomRadius<<" "<<"Length (mm) = "<<phantomLength<<G4endl;
740 G4cout<<"---------------Phantom position (mm) : "<<phantomPosition<<G4endl;
741 }
742 else if(PhantomType == "Phantom_spatialResolution"){
743 ////////////////// phantom for point source (For spatial resolution evaluation) /////////////////////////////////
744 //The position of the point source phantom can be changed in the macro file.
745 //It has cylindrical shape and its radius and length can be set in the macro file
746
747 //Define its solid shape
748 G4Tubs* phantom = new G4Tubs("Phantom_point", 0., phantomRadius,phantomLength/2,0*deg,360*deg);
749
750 //Define its logical volume
751 phantom_logicalV = new G4LogicalVolume(phantom, polyethylene_NEMA , "phantom_logicalV",0,0,0);
752
753 //place the phantom
754 phantomPosition = G4ThreeVector(phantomPosition.x(),phantomPosition.y(),phantomPosition.z());
755
756 //Define its physical volume
757 phantom_physicalV = new G4PVPlacement(0, phantomPosition, phantom_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
758 G4cout<<"---------------Phantom dimension: "<<"radius (mm) = "<<phantomRadius<<" "<<", Length (mm) = "<<phantomLength<<G4endl;
759 G4cout<<"---------------Phantom position: "<<phantomPosition<< "mm"<<G4endl;
760 }
761
762 else if (PhantomType == "Normalization")
763 {
764 //The normalization phantom is a hallow cylindrical phantom to represent a (rotating) line source. The source is confined in the phantom.
765 //The thickness of the phantom is 3 mm, and its diameter is 350 mm.
766
767 //Define the solid shape.
768 G4Tubs* phantom = new G4Tubs("Phantom", phantomRadius - 3*mm, phantomRadius,phantomLength/2,0*deg,360*deg);
769
770 //Define its logical volume.
771 phantom_logicalV = new G4LogicalVolume(phantom, polyethylene_NEMA, "phantom_logicalV",0,0,0);
772
773 //Define its physical volume
774 phantom_physicalV = new G4PVPlacement(0, phantomPosition, phantom_logicalV, "phantom_physicalV", worldLogical, false, 0,fCheckOverlaps);
775
776
777 G4cout<<"---------------Phantom dimension: "<<"radius (mm) = "<<phantomRadius<<" "<<"Length = "<<phantomLength<<G4endl;
778 G4cout<<"---------------Phantom position : "<<phantomPosition<< "mm"<<G4endl;
779
780 }
781 //-------------------------------------================- No Phantom Set -==================------------------
782 else
783 {
784 G4cerr << "****************************************\nERROR: Phantom Remains: " << PhantomType <<"\n****************************************" << G4endl;
785 exit(0);
786
787 }
788 ///////////Visualization of phantoms///////////
789
790 G4VisAttributes* phantomVisAtt;
791 phantomVisAtt = new G4VisAttributes(G4Colour(0.88,0.55,1.0));
792 phantomVisAtt->SetVisibility (true);
793 phantomVisAtt->SetForceWireframe (true);
794 phantom_logicalV->SetVisAttributes (phantomVisAtt);
795 //phantom_logicalV->SetVisAttributes (G4VisAttributes::Invisible);
796 }
797
798 //Change the type of the phantom via .mac file
799 void doiPETDetectorConstruction::ChangePhantom(G4String NewPhantomtype)
800 {
801 PhantomType = NewPhantomtype;
802 }
803
804 //Change position of the phantom via .mac file
805 void doiPETDetectorConstruction::SetPhantomPosition(G4ThreeVector NewphantomPosition)
806 {
807 phantomPosition = NewphantomPosition;
808 }
809
810 //Change the radius of the phantom via .mac file
811 void doiPETDetectorConstruction::SetPhantomRadius(G4double newPhantomRadius){
812 phantomRadius = newPhantomRadius;
813 }
814
815 //Change the length of the phantom via .mac file
816 void doiPETDetectorConstruction::SetPhantomLength(G4double newPhantomLength){
817 phantomLength = newPhantomLength;
818 }
819 //
820 void doiPETDetectorConstruction::SetNumberOfSleeves(G4int NewNumOfSleeves){
821 numOfSleeves = NewNumOfSleeves;
822 //G4cout<<"numOfSleeves "<<numOfSleeves<<G4endl;
823 }
824