source: trunk/examples/advanced/medical_linac/src/MedLinacDetectorConstruction.cc @ 1282

Last change on this file since 1282 was 807, checked in by garnier, 16 years ago

update
File size: 18.1 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//
27// $Id: MedLinacDetectorConstruction.cc,v 1.12 2007/07/06 21:39:05 mpiergen Exp $
28//
29// Code developed by: M. Piergentili
30//
31//
32//------------------------------ beam line along z axis---------------------
33
34#include "MedLinacDetectorMessenger.hh"
35#include "MedLinacDetectorConstruction.hh"
36#include "MedLinacPhantomROGeometry.hh"
37#include "MedLinacPhantomSD.hh"
38#include "MedLinacVGeometryComponent.hh"
39#include "MedLinacHead.hh"
40#include "MedLinacDecorator.hh"
41#include "MedLinacTargetAndFilterDecorator.hh"
42#include "MedLinacMLCDecorator.hh"
43#include "G4CSGSolid.hh"
44#include "G4RotationMatrix.hh"
45#include "G4Transform3D.hh"
46#include "G4Material.hh"
47#include "G4MaterialPropertiesTable.hh"
48#include "G4MaterialTable.hh"
49#include "G4MaterialPropertyVector.hh"
50#include "G4Element.hh"
51#include "G4ElementTable.hh"
52#include "G4Box.hh"
53#include "G4Cons.hh"
54#include "G4Tubs.hh"
55#include "G4LogicalVolume.hh"
56#include "G4ThreeVector.hh"
57#include "G4PVPlacement.hh"
58#include "G4PVParameterised.hh"
59#include "globals.hh"
60#include "G4SDManager.hh"
61#include "G4RunManager.hh"
62#include "Randomize.hh"
63#include "G4TransportationManager.hh"
64#include "G4UserLimits.hh"
65#include "G4GeometryManager.hh"
66#include "G4BooleanSolid.hh"
67#include "G4SubtractionSolid.hh"
68#include "G4VSolid.hh"
69#include "G4Region.hh"
70
71#include "G4PhysicalVolumeStore.hh"
72#include "G4LogicalVolumeStore.hh"
73#include "G4SolidStore.hh"
74#include "G4RegionStore.hh"
75
76#include "G4VisAttributes.hh"
77#include "G4Colour.hh"
78
79#include "G4ios.hh"
80
81MedLinacDetectorConstruction::MedLinacDetectorConstruction(G4String SDName)
82  : phantomSD(0),phantomROGeometry(0),
83    experimentalHall_log(0),vacuumBlock_log(0),
84    JawY1_log(0),JawY2_log(0),
85    JawX1_log(0),JawX2_log(0),
86    Phantom_log(0),
87    experimentalHall_phys(0),vacuumBlock_phys(0),
88    JawY1_phys(0), JawY2_phys(0),
89    JawX1_phys(0), JawX2_phys(0),
90    Phantom_phys(0)
91{
92  phantomDim = 15.*cm;
93   //PhantomDimensionX = PhantomDim_x;
94   //PhantomDimensionY = PhantomDim_y;
95   //PhantomDimensionZ = PhantomDim_z;
96
97  //numberOfVoxelsAlongX=150;
98  //numberOfVoxelsAlongY=150;
99  //numberOfVoxelsAlongZ=150;
100
101  numberOfVoxelsAlongX=0;
102  numberOfVoxelsAlongY=0;
103  numberOfVoxelsAlongZ=0;
104  maxStep = 0;
105
106  sensitiveDetectorName = SDName;
107
108 // default parameter values
109
110  fieldX1 = -2.5*cm;
111  fieldX2 =  2.5*cm;
112  fieldY1 = -2.5*cm;
113  fieldY2 =  2.5*cm;
114
115
116  detectorMessenger = new MedLinacDetectorMessenger(this);
117  pHead = new MedLinacHead();
118  decorator = new MedLinacTargetAndFilterDecorator(pHead);
119  decorator1 = new MedLinacMLCDecorator(pHead);
120  ComputeDimVoxel();
121
122
123  }
124
125MedLinacDetectorConstruction* MedLinacDetectorConstruction::instance = 0;
126
127MedLinacDetectorConstruction* MedLinacDetectorConstruction::GetInstance(G4String SDName)
128{
129  if (instance == 0)
130    {
131      instance = new MedLinacDetectorConstruction(SDName);
132     
133    }
134  return instance;
135}
136
137MedLinacDetectorConstruction::~MedLinacDetectorConstruction()
138{
139  delete detectorMessenger;
140  delete pHead;
141  delete decorator;
142  delete decorator1;
143  }
144
145G4VPhysicalVolume* MedLinacDetectorConstruction::Construct()
146{
147  return ConstructGeom();
148}
149
150G4VPhysicalVolume* MedLinacDetectorConstruction::ConstructGeom ()
151{
152  // Cleanup old geometry
153
154  G4GeometryManager::GetInstance()->OpenGeometry();
155  G4PhysicalVolumeStore::GetInstance()->Clean();
156  G4LogicalVolumeStore::GetInstance()->Clean();
157  G4SolidStore::GetInstance()->Clean();
158
159  //    materials
160
161  G4double a;  // atomic mass
162  G4double z;  // atomic number
163  G4double density;
164  G4int ncomponents;
165  G4int natoms;
166  G4String name;
167  G4String symbol;
168  G4double fractionmass;                                           
169  G4double massOfMole;   
170  G4double pressure; 
171  G4double temperature; 
172
173  a = 14.00674*g/mole;
174  G4Element* elN = new G4Element(name="Nitrogen" ,symbol="N", z=7., a);
175
176  a = 16.00*g/mole;
177  G4Element* elO = new G4Element(name="Oxygen",symbol="O", z=8., a);
178
179  a = 1.00794*g/mole;
180  G4Element* elH = new G4Element(name="Hydrogen",symbol="H", z=1., a);
181 
182  density = 1.290*mg/cm3;
183  G4Material* Air = new G4Material(name="Air",density, ncomponents=2);
184  Air->AddElement(elN, fractionmass=70*perCent);
185  Air->AddElement(elO, fractionmass=30*perCent);
186
187  density = 1.000*g/cm3;
188  G4Material* H2O = new G4Material(name="Water",density, ncomponents=2);
189  H2O->AddElement(elH, natoms=2);
190  H2O->AddElement(elO, natoms=1);
191  H2O->GetIonisation()->SetMeanExcitationEnergy(75.0*eV);
192
193  //density = 19.3*g/cm3;
194  density = 18.*g/cm3;
195  a = 183.85*g/mole;
196  G4Material* W = new G4Material(name="Tungsten"  , z=74., a, density);
197
198  massOfMole = 1.008*g/mole;
199  density = 1.e-25*g/cm3;
200  temperature = 2.73*kelvin;
201  pressure = 3.e-18*pascal; 
202  G4Material* Vacuum = new G4Material("interGalactic", z=1.,massOfMole, 
203                                      density, kStateGas,temperature, pressure); 
204
205  //------------------------------ beam line along z axis------------------
206
207  //---------jaws position--------
208
209  //----1Y------------------------
210
211  G4double thetaY1 = std::atan(-fieldY1/(100.*cm));// in rad
212  G4double JawY1Pos_y = -((3.9*cm+28.*cm)*std::sin(thetaY1)+5.*cm*std::cos(thetaY1));
213  G4double JawY1Pos_z =115.*cm-((3.9*cm+28.*cm)*std::cos(thetaY1))+5.*cm*std::sin(thetaY1);
214  //----2Y-------------------------
215  G4double thetaY2 = std::atan(fieldY2/(100.*cm)); // in rad
216  G4double JawY2Pos_y = (3.9*cm+28.*cm)*std::sin(thetaY2)+5.*cm*std::cos(thetaY2);
217  G4double JawY2Pos_z = 115.*cm-(3.9*cm+28.*cm)*std::cos(thetaY2)+5.*cm*std::sin(thetaY2);
218  //----1X-------------------------
219  G4double thetaX1 = std::atan(-fieldX1/(100.*cm));
220  G4double JawX1Pos_x = -((36.7*cm+3.9*cm)*std::sin(thetaX1)+5.*cm*std::cos(thetaX1));
221  G4double JawX1Pos_z =115.*cm-(36.7*cm+3.9*cm)*std::cos(thetaX1)+5.*cm*std::sin(thetaX1);
222   //----2X-------------------------
223  G4double thetaX2 = std::atan(fieldX2/(100.*cm));
224  G4double JawX2Pos_x = (36.7*cm+3.9*cm)*std::sin(thetaX2)+5.*cm*std::cos(thetaX2);
225  G4double JawX2Pos_z = 115.*cm-(36.7*cm+3.9*cm)*std::cos(thetaX2)+5.*cm*std::sin(thetaX2);
226  //---------rotation jaw1Y--------
227
228  G4RotationMatrix* rotatejaw1Y=new G4RotationMatrix();
229  rotatejaw1Y->rotateX(thetaY1);
230
231//---------rotation jaw2Y--------
232
233  G4RotationMatrix* rotatejaw2Y=new G4RotationMatrix();
234  rotatejaw2Y->rotateX(-thetaY2);
235
236  //---------rotation jaw1X--------
237
238  G4RotationMatrix*  rotatejaw1X=new G4RotationMatrix();
239  rotatejaw1X->rotateY(-thetaX1);
240
241  //---------rotation jaw2X--------
242
243  G4RotationMatrix*  rotatejaw2X=new G4RotationMatrix();
244  rotatejaw2X->rotateY(thetaX2);
245
246  //---------------colors----------
247
248  G4Colour  white   (1.0, 1.0, 1.0);
249   G4Colour  magenta (1.0, 0.0, 1.0); 
250   G4Colour  lblue   (0.20, .50, .85);
251
252  //------------------------------------------------------ volumes
253
254  //------------------------------ experimental hall (world volume)
255  //------------------------------ beam line along z axis-------SSD=100cm-------------- 
256
257//------------------------------ experimental hall (world volume)
258  G4double expHall_x = 3.0*m;
259  G4double expHall_y = 3.0*m;
260  G4double expHall_z = 3.0*m;
261  G4Box* experimentalHall_box
262    = new G4Box("expHall_box",expHall_x,expHall_y,expHall_z);
263  experimentalHall_log = new G4LogicalVolume(experimentalHall_box,
264                                             Air,"expHall_log",0,0,0);
265  experimentalHall_phys = new G4PVPlacement(0,G4ThreeVector(),
266                                      "expHall",experimentalHall_log,0,false,0);
267
268
269 
270  //------------------------vacuum box------------------------
271
272  G4double vacudim_x = 9.*cm;
273  G4double vacudim_y = 9.*cm;
274  G4double vacudim_z = 8.755*cm;
275  G4Box* vacuumBlock_box = new G4Box("vacuumBlock_box",vacudim_x,
276                                     vacudim_y,vacudim_z);
277  vacuumBlock_log = new G4LogicalVolume(vacuumBlock_box,
278                                        Vacuum,"vacuumBlock_log",0,0,0);
279  G4double vacublockPos_x = 0.0*m;
280  G4double vacublockPos_y = 0.0*m;
281  G4double vacublockPos_z = 114.755*cm;
282  vacuumBlock_phys = new G4PVPlacement(0,
283             G4ThreeVector(vacublockPos_x,vacublockPos_y,vacublockPos_z),
284             "vacuBlock",vacuumBlock_log,experimentalHall_phys,false,0);
285
286
287  //----------------Jaw Y1 collimator---------------
288
289  G4double JawY1Dim_x = 10.0*cm;
290  G4double JawY1Dim_y = 5.0*cm;
291  G4double JawY1Dim_z = 3.9*cm;
292  G4Box* JawY1 = new G4Box("JawY1",JawY1Dim_x,
293                                  JawY1Dim_y,JawY1Dim_z);
294  JawY1_log = new G4LogicalVolume(JawY1,
295                                             W,"JawY1_log",0,0,0);
296
297
298
299  G4double JawY1Pos_x = 0.0*cm;
300  //G4double JawY1Pos_y = fieldY1;
301  //G4double JawY1Pos_z = 68.1*cm;
302  JawY1_phys = new G4PVPlacement(rotatejaw1Y,
303             G4ThreeVector(JawY1Pos_x,JawY1Pos_y,JawY1Pos_z),
304             "JawY1",JawY1_log,experimentalHall_phys,false,0);
305
306  //----------------Jaw Y2 collimator---------------
307
308  G4double JawY2Dim_x = 10.0*cm;
309  G4double JawY2Dim_y = 5.0*cm;
310  G4double JawY2Dim_z = 3.9*cm;
311  G4Box* JawY2 = new G4Box("JawY2",JawY2Dim_x,
312                                  JawY2Dim_y,JawY2Dim_z);
313  JawY2_log = new G4LogicalVolume(JawY2,
314                                             W,"JawY1_log",0,0,0);
315  G4double JawY2Pos_x = 0.0*cm;
316  //G4double JawY2Pos_y = fieldY2;
317  //G4double JawY2Pos_z = 68.1*cm;
318  JawY2_phys = new G4PVPlacement(rotatejaw2Y,
319             G4ThreeVector(JawY2Pos_x,JawY2Pos_y,JawY2Pos_z),
320             "JawY2",JawY2_log,experimentalHall_phys,false,0);
321
322
323 //----------------Jaw X1 collimator---------------
324
325  G4double JawX1Dim_x = 5.0*cm;
326  G4double JawX1Dim_y = 10.0*cm;
327  G4double JawX1Dim_z = 3.9*cm;
328  G4Box* JawX1 = new G4Box("JawX1",JawX1Dim_x,
329                                 JawX1Dim_y,JawX1Dim_z);
330  JawX1_log = new G4LogicalVolume(JawX1,
331                                             W,"JawX1_log",0,0,0);
332  //G4cout <<"_________________DetectorConstruction fieldX1 "<< fieldX1/cm << G4endl;
333
334  //G4double JawX1Pos_x = fieldX1;
335  G4double JawX1Pos_y = 0.0*cm;
336  //G4double JawX1Pos_z = 54.2*cm;
337  JawX1_phys = new G4PVPlacement(rotatejaw1X,
338             G4ThreeVector(JawX1Pos_x,JawX1Pos_y,JawX1Pos_z),
339             "JawX1",JawX1_log,experimentalHall_phys,false,0);
340
341  //----------------Jaw X2 collimator---------------
342
343  G4double JawX2Dim_x = 5.0*cm;
344  G4double JawX2Dim_y = 10.0*cm;
345  G4double JawX2Dim_z = 3.9*cm;
346  G4Box* JawX2 = new G4Box("JawX2",JawX2Dim_x,
347                                  JawX2Dim_y,JawX2Dim_z);
348  JawX2_log = new G4LogicalVolume(JawX2,
349                                             W,"JawX1_log",0,0,0);
350  //G4double JawX2Pos_x = fieldX2;
351  G4double JawX2Pos_y = 0.0*cm;
352  //G4double JawX2Pos_z = 54.2*cm;
353  JawX2_phys = new G4PVPlacement(rotatejaw2X,
354             G4ThreeVector(JawX2Pos_x,JawX2Pos_y,JawX2Pos_z),
355             "JawX2",JawX2_log,experimentalHall_phys,false,0);
356
357  //----------------Phantom---------
358  //phantomDim_x = 15.*cm;
359  //phantomDim_y = 15.*cm;
360  //phantomDim_z = 15.*cm;
361  phantomDim_x = phantomDim;                                                                                                             
362  phantomDim_y = phantomDim;                                                                                                             
363  phantomDim_z = phantomDim; 
364
365  G4Box* Phantom = new G4Box("Phantom",phantomDim_x,
366                                  phantomDim_y,phantomDim_z);
367  Phantom_log = new G4LogicalVolume(Phantom,
368                                             H2O,"Phantom_log",0,0,0);
369  G4double PhantomPos_x = 0.0*m;
370  G4double PhantomPos_y = 0.0*m;
371  G4double PhantomPos_z = 0.0*cm;
372  Phantom_phys = new G4PVPlacement(0,
373             G4ThreeVector(PhantomPos_x,PhantomPos_y,PhantomPos_z),
374             "Phantom_phys",Phantom_log,experimentalHall_phys,false,0);
375
376  numberOfVoxelsAlongX=numberOfVoxels;
377  numberOfVoxelsAlongY=numberOfVoxels;
378  numberOfVoxelsAlongZ=numberOfVoxels;
379
380  PrintParameters();   
381
382//--------- Visualization attributes -------------------------------
383
384
385 
386  G4VisAttributes* simpleH2OVisAtt= new G4VisAttributes(lblue);
387  simpleH2OVisAtt->SetVisibility(true);
388  simpleH2OVisAtt->SetForceSolid(true);
389   Phantom_log->SetVisAttributes(simpleH2OVisAtt);
390
391   G4VisAttributes* simpleVacuumWVisAtt= new G4VisAttributes(white);
392   simpleVacuumWVisAtt->SetVisibility(true);
393   simpleVacuumWVisAtt->SetForceWireframe(true);
394   vacuumBlock_log->SetVisAttributes(simpleVacuumWVisAtt);
395
396   G4VisAttributes* simpleTungstenSVisAtt= new G4VisAttributes(magenta);
397   simpleTungstenSVisAtt->SetVisibility(true);
398   simpleTungstenSVisAtt->SetForceSolid(true);
399   JawY1_log->SetVisAttributes(simpleTungstenSVisAtt);
400   JawY2_log->SetVisAttributes(simpleTungstenSVisAtt);
401   JawX1_log->SetVisAttributes(simpleTungstenSVisAtt);
402   JawX2_log->SetVisAttributes(simpleTungstenSVisAtt);
403
404 
405   G4VisAttributes* simpleWorldVisAtt= new G4VisAttributes(white);
406   simpleWorldVisAtt->SetVisibility(false);
407   experimentalHall_log ->SetVisAttributes(simpleWorldVisAtt);
408
409   //if you want not to see phantom, jaws e vacuum ...:
410
411   //Phantom_log->SetVisAttributes(simpleWorldVisAtt);
412   //JawY1_log->SetVisAttributes(simpleWorldVisAtt);
413   //JawY2_log->SetVisAttributes(simpleWorldVisAtt);
414   //JawX1_log->SetVisAttributes(simpleWorldVisAtt);
415   //JawX2_log->SetVisAttributes(simpleWorldVisAtt);
416   //vacuumBlock_log->SetVisAttributes(simpleWorldVisAtt);
417
418  //   Sets a max Step length in the detector
419   //G4double maxStep = 0.2*mm;
420   Phantom_log->SetUserLimits(new G4UserLimits(maxStep));
421
422   ConstructVolume();
423
424   ConstructSensitiveDetector();
425
426   //G4cout <<"????????????????????numberOfVoxels "<< numberOfVoxels <<G4endl ;
427   //G4cout <<"????????????????????numberOfVoxelsAlongX "<< numberOfVoxelsAlongX <<G4endl ;
428   //G4cout <<"????????????????????numberOfVoxelsAlongY "<< numberOfVoxelsAlongY <<G4endl ;
429   //G4cout <<"????????????????????numberOfVoxelsAlongZ "<< numberOfVoxelsAlongZ <<G4endl ;
430
431
432   //G4cout <<"????????????????????maxStep "<< maxStep/mm << " mm " <<G4endl ;
433
434
435
436
437   return experimentalHall_phys;
438}
439
440void MedLinacDetectorConstruction::PrintParameters()
441{ 
442  //G4cout <<"jaws1 x position "<< fieldX1/cm << " cm "<<G4endl ;
443  //G4cout <<"jaws2 x position "<< fieldX2/cm << " cm "<<G4endl ;
444  //G4cout <<"jaws1 y position "<< fieldY1/cm << " cm "<<G4endl ;
445  //G4cout <<"jaws2 y position "<< fieldY2/cm << " cm "<<G4endl ;
446  G4cout <<"************************phantom dimension "<< phantomDim_x/cm << " cm "<<G4endl ;
447  G4cout <<"************************phantom dimension "<< phantomDim/cm << " cm "<<G4endl;
448  //G4cout <<"************************numberOfVoxelsAlongX "<< numberOfVoxelsAlongX <<G4endl ;
449  //G4cout <<"************************numberOfVoxelsAlongY "<< numberOfVoxelsAlongY <<G4endl ;
450  //G4cout <<"************************numberOfVoxelsAlongZ "<< numberOfVoxelsAlongZ <<G4endl ;
451  //G4cout <<"************************maxStep "<< maxStep/mm << " mm " <<G4endl;
452}
453
454
455
456void MedLinacDetectorConstruction::SetJawX1Pos_x (G4double val)
457{ 
458  fieldX1 = val;
459}
460
461
462void MedLinacDetectorConstruction::SetJawX2Pos_x (G4double val)
463{
464  fieldX2 = val;
465}
466
467void MedLinacDetectorConstruction::SetJawY1Pos_y (G4double val)
468{
469  fieldY1 = val;
470}
471
472void MedLinacDetectorConstruction::SetJawY2Pos_y (G4double val)
473{
474  fieldY2 = val;
475  //G4cout <<"==============================DetectorConstruction JawY2  "<< fieldY2/cm<<"cm"<<G4endl;
476}
477
478void MedLinacDetectorConstruction::SetPhantomDim (G4double val)
479{
480  phantomDim = val;
481  G4cout <<"==============================phantomDim  "<< phantomDim/mm<<"mm"<<G4endl;
482}
483
484void MedLinacDetectorConstruction::SetNumberOfVoxels (G4int val)
485{
486  numberOfVoxels = val;
487  //G4cout <<"==============================numberOfVoxels  "<< numberOfVoxels<<G4endl;
488}
489
490void MedLinacDetectorConstruction::SetMaxStep (G4double val)
491{
492  maxStep = val;
493  //G4cout <<"==============================maxStep  "<< maxStep/mm<<"mm"<<G4endl;
494}
495void MedLinacDetectorConstruction::UpdateGeometry()
496{ 
497  //G4cout <<"@@@@@@@@@@@@@@@@@@@@@UpdateGeometry "<< G4endl;
498
499  G4RunManager::GetRunManager()->DefineWorldVolume(ConstructGeom());
500}
501
502
503
504void MedLinacDetectorConstruction::ConstructVolume()
505{
506  pHead ->ConstructComponent(experimentalHall_phys,vacuumBlock_phys);
507  decorator ->ConstructComponent( experimentalHall_phys,vacuumBlock_phys);
508  decorator1 ->ConstructComponent( experimentalHall_phys,vacuumBlock_phys);
509
510  }
511
512void  MedLinacDetectorConstruction::ConstructSensitiveDetector()
513// Sensitive Detector and ReadOut geometry definition
514{ 
515
516  G4SDManager* pSDManager = G4SDManager::GetSDMpointer();
517
518  if(!phantomSD)
519  {
520    phantomSD = new MedLinacPhantomSD (sensitiveDetectorName);
521
522    pSDManager->AddNewDetector(phantomSD);
523  }
524
525 G4String ROGeometryName = "PhantomROGeometry";
526    phantomROGeometry = new MedLinacPhantomROGeometry (ROGeometryName, phantomDim_x,phantomDim_y,phantomDim_z,
527                                                  numberOfVoxelsAlongX,numberOfVoxelsAlongY,numberOfVoxelsAlongZ);
528
529    phantomROGeometry->BuildROGeometry();
530    phantomSD->SetROgeometry(phantomROGeometry);
531
532    Phantom_log->SetSensitiveDetector(phantomSD);
533}
Note: See TracBrowser for help on using the repository browser.