Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc (Version 11.3.0) and /examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc (Version 10.7.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // Hadrontherapy advanced example for Geant4       26 // Hadrontherapy advanced example for Geant4
 27 // See more at: https://twiki.cern.ch/twiki/bi     27 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
 28                                                    28 
 29 #include "G4UnitsTable.hh"                         29 #include "G4UnitsTable.hh"
 30 #include "G4SDManager.hh"                          30 #include "G4SDManager.hh"
 31 #include "G4RunManager.hh"                         31 #include "G4RunManager.hh"
 32 #include "G4GeometryManager.hh"                    32 #include "G4GeometryManager.hh"
 33 #include "G4SolidStore.hh"                         33 #include "G4SolidStore.hh"
 34 #include "G4PhysicalVolumeStore.hh"                34 #include "G4PhysicalVolumeStore.hh"
 35 #include "G4LogicalVolumeStore.hh"                 35 #include "G4LogicalVolumeStore.hh"
 36 #include "G4Box.hh"                                36 #include "G4Box.hh"
 37 #include "G4LogicalVolume.hh"                      37 #include "G4LogicalVolume.hh"
 38 #include "G4ThreeVector.hh"                        38 #include "G4ThreeVector.hh"
 39 #include "G4PVPlacement.hh"                        39 #include "G4PVPlacement.hh"
 40 #include "globals.hh"                              40 #include "globals.hh"
 41 #include "G4Transform3D.hh"                        41 #include "G4Transform3D.hh"
 42 #include "G4RotationMatrix.hh"                     42 #include "G4RotationMatrix.hh"
 43 #include "G4Colour.hh"                             43 #include "G4Colour.hh"
 44 #include "G4UserLimits.hh"                         44 #include "G4UserLimits.hh"
 45 #include "G4UnitsTable.hh"                         45 #include "G4UnitsTable.hh"
 46 #include "G4VisAttributes.hh"                      46 #include "G4VisAttributes.hh"
 47 #include "G4NistManager.hh"                        47 #include "G4NistManager.hh"
 48 #include "HadrontherapyDetectorConstruction.hh     48 #include "HadrontherapyDetectorConstruction.hh"
 49 #include "HadrontherapyDetectorROGeometry.hh"      49 #include "HadrontherapyDetectorROGeometry.hh"
 50 #include "HadrontherapyDetectorMessenger.hh"       50 #include "HadrontherapyDetectorMessenger.hh"
 51 #include "HadrontherapyDetectorSD.hh"              51 #include "HadrontherapyDetectorSD.hh"
 52 #include "HadrontherapyMatrix.hh"                  52 #include "HadrontherapyMatrix.hh"
 53 #include "HadrontherapyLet.hh"                     53 #include "HadrontherapyLet.hh"
 54 #include "PassiveProtonBeamLine.hh"                54 #include "PassiveProtonBeamLine.hh"
 55 #include "BESTPassiveProtonBeamLine.hh"            55 #include "BESTPassiveProtonBeamLine.hh"
 56 #include "HadrontherapyMatrix.hh"                  56 #include "HadrontherapyMatrix.hh"
 57                                                    57 
 58 #include "HadrontherapyRBE.hh"                     58 #include "HadrontherapyRBE.hh"
 59 #include "G4SystemOfUnits.hh"                      59 #include "G4SystemOfUnits.hh"
 60                                                    60 
 61 #include <cmath>                                   61 #include <cmath>
 62                                                    62 
 63                                                    63 
 64                                                    64 
 65 HadrontherapyDetectorConstruction* Hadronthera     65 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::instance = 0;
 66 //////////////////////////////////////////////     66 /////////////////////////////////////////////////////////////////////////////
 67 HadrontherapyDetectorConstruction::Hadronthera     67 HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
 68 : motherPhys(physicalTreatmentRoom), // pointe     68 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
 69 detectorSD(0), detectorROGeometry(0), matrix(0     69 detectorSD(0), detectorROGeometry(0), matrix(0),
 70 phantom(0), detector(0),                           70 phantom(0), detector(0),
 71 phantomLogicalVolume(0), detectorLogicalVolume     71 phantomLogicalVolume(0), detectorLogicalVolume(0),
 72 phantomPhysicalVolume(0), detectorPhysicalVolu     72 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
 73 aRegion(0)                                         73 aRegion(0)
 74 {                                                  74 {
 75                                                    75     
 76                                                    76     
 77     /* NOTE! that the HadrontherapyDetectorCon     77     /* NOTE! that the HadrontherapyDetectorConstruction class
 78      * does NOT inherit from G4VUserDetectorCo     78      * does NOT inherit from G4VUserDetectorConstruction G4 class
 79      * So the Construct() mandatory virtual me     79      * So the Construct() mandatory virtual method is inside another geometric class
 80      * like the passiveProtonBeamLIne, ...         80      * like the passiveProtonBeamLIne, ...
 81      */                                            81      */
 82                                                    82     
 83     // Messenger to change parameters of the p     83     // Messenger to change parameters of the phantom/detector geometry
 84     detectorMessenger = new HadrontherapyDetec     84     detectorMessenger = new HadrontherapyDetectorMessenger(this);
 85                                                    85     
 86     // Default detector voxels size                86     // Default detector voxels size
 87     // 200 slabs along the beam direction (X)      87     // 200 slabs along the beam direction (X)
 88     sizeOfVoxelAlongX = 200 *um;                   88     sizeOfVoxelAlongX = 200 *um;
 89     sizeOfVoxelAlongY = 4 *cm;                     89     sizeOfVoxelAlongY = 4 *cm;
 90     sizeOfVoxelAlongZ = 4 *cm;                     90     sizeOfVoxelAlongZ = 4 *cm;
 91                                                    91     
 92     // Define here the material of the water p     92     // Define here the material of the water phantom and of the detector
 93     SetPhantomMaterial("G4_WATER");                93     SetPhantomMaterial("G4_WATER");
 94     // Construct geometry (messenger commands)     94     // Construct geometry (messenger commands)
 95     // SetDetectorSize(4.*cm, 4.*cm, 4.*cm);       95     // SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
 96     SetDetectorSize(4. *cm, 4. *cm, 4. *cm);       96     SetDetectorSize(4. *cm, 4. *cm, 4. *cm);
 97     SetPhantomSize(40. *cm, 40. *cm, 40. *cm);     97     SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
 98                                                    98     
 99     SetPhantomPosition(G4ThreeVector(20. *cm,      99     SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
100     SetDetectorToPhantomPosition(G4ThreeVector    100     SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 0. *cm, 0. *cm));
101     SetDetectorPosition();                        101     SetDetectorPosition();
102     //GetDetectorToWorldPosition();               102     //GetDetectorToWorldPosition();
103                                                   103     
104     // Write virtual parameters to the real on    104     // Write virtual parameters to the real ones and check for consistency
105     UpdateGeometry();                             105     UpdateGeometry();
106                                                   106     
107                                                   107     
108                                                   108     
109 }                                                 109 }
110                                                   110 
111 //////////////////////////////////////////////    111 /////////////////////////////////////////////////////////////////////////////
112 HadrontherapyDetectorConstruction::~Hadronther    112 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
113 {                                                 113 {
114     delete detectorROGeometry;                    114     delete detectorROGeometry;
115     delete matrix;                                115     delete matrix;
116     delete detectorMessenger;                     116     delete detectorMessenger;
117 }                                                 117 }
118                                                   118 
119 //////////////////////////////////////////////    119 /////////////////////////////////////////////////////////////////////////////
120 HadrontherapyDetectorConstruction* Hadronthera    120 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::GetInstance()
121 {                                                 121 {
122     return instance;                              122     return instance;
123 }                                                 123 }
124                                                   124 
125 //////////////////////////////////////////////    125 /////////////////////////////////////////////////////////////////////////////
126 // ConstructPhantom() is the method that const    126 // ConstructPhantom() is the method that construct a water box (called phantom
127 // (or water phantom) in the usual Medical phy    127 // (or water phantom) in the usual Medical physicists slang).
128 // A water phantom can be considered a good ap    128 // A water phantom can be considered a good approximation of a an human body.
129 void HadrontherapyDetectorConstruction::Constr    129 void HadrontherapyDetectorConstruction::ConstructPhantom()
130 {                                                 130 {
131     // Definition of the solid volume of the P    131     // Definition of the solid volume of the Phantom
132     phantom = new G4Box("Phantom",                132     phantom = new G4Box("Phantom",
133                         phantomSizeX/2,           133                         phantomSizeX/2,
134                         phantomSizeY/2,           134                         phantomSizeY/2,
135                         phantomSizeZ/2);          135                         phantomSizeZ/2);
136                                                   136     
137     // Definition of the logical volume of the    137     // Definition of the logical volume of the Phantom
138     phantomLogicalVolume = new G4LogicalVolume    138     phantomLogicalVolume = new G4LogicalVolume(phantom,
139                                                   139                                                phantomMaterial,
140                                                   140                                                "phantomLog", 0, 0, 0);
141                                                   141     
142     // Definition of the physics volume of the    142     // Definition of the physics volume of the Phantom
143     phantomPhysicalVolume = new G4PVPlacement(    143     phantomPhysicalVolume = new G4PVPlacement(0,
144                                                   144                                               phantomPosition,
145                                                   145                                               "phantomPhys",
146                                                   146                                               phantomLogicalVolume,
147                                                   147                                               motherPhys,
148                                                   148                                               false,
149                                                   149                                               0);
150                                                   150     
151     // Visualisation attributes of the phantom    151     // Visualisation attributes of the phantom
152     red = new G4VisAttributes(G4Colour(255/255    152     red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
153     red -> SetVisibility(true);                   153     red -> SetVisibility(true);
154     red -> SetForceSolid(true);                   154     red -> SetForceSolid(true);
155     red -> SetForceWireframe(true);               155     red -> SetForceWireframe(true);
156     phantomLogicalVolume -> SetVisAttributes(r    156     phantomLogicalVolume -> SetVisAttributes(red);
157 }                                                 157 }
158                                                   158 
159 //////////////////////////////////////////////    159 /////////////////////////////////////////////////////////////////////////////
160 // ConstructDetector() is the method the recon    160 // ConstructDetector() is the method the reconstruct a detector region
161 // inside the water phantom. It is a volume, l    161 // inside the water phantom. It is a volume, located inside the water phantom.
162 //                                                162 //
163 //           **************************           163 //           **************************
164 //           *    water phantom       *           164 //           *    water phantom       *
165 //           *                        *           165 //           *                        *
166 //           *                        *           166 //           *                        *
167 //           *---------------         *           167 //           *---------------         *
168 //  Beam     *              -         *           168 //  Beam     *              -         *
169 //  ----->   * detector     -         *           169 //  ----->   * detector     -         *
170 //           *              -         *           170 //           *              -         *
171 //           *---------------         *           171 //           *---------------         *
172 //           *                        *           172 //           *                        *
173 //           *                        *           173 //           *                        *
174 //           *                        *           174 //           *                        *
175 //           **************************           175 //           **************************
176 //                                                176 //
177 // The detector can be dived in slices or voxe    177 // The detector can be dived in slices or voxelized
178 // and inside it different quantities (dose di    178 // and inside it different quantities (dose distribution, fluence distribution, LET, etc)
179 // can be stored.                                 179 // can be stored.
180 void HadrontherapyDetectorConstruction::Constr    180 void HadrontherapyDetectorConstruction::ConstructDetector()
181                                                   181 
182 {                                                 182 {
183     // Definition of the solid volume of the D    183     // Definition of the solid volume of the Detector
184     detector = new G4Box("Detector",              184     detector = new G4Box("Detector",
185                                                   185                          
186                          phantomSizeX/2,          186                          phantomSizeX/2,
187                                                   187                          
188                          phantomSizeY/2,          188                          phantomSizeY/2,
189                                                   189                          
190                          phantomSizeZ/2);         190                          phantomSizeZ/2);
191                                                   191     
192     // Definition of the logic volume of the P    192     // Definition of the logic volume of the Phantom
193     detectorLogicalVolume = new G4LogicalVolum    193     detectorLogicalVolume = new G4LogicalVolume(detector,
194                                                   194                                                 detectorMaterial,
195                                                   195                                                 "DetectorLog",
196                                                   196                                                 0,0,0);
197     // Definition of the physical volume of th    197     // Definition of the physical volume of the Phantom
198     detectorPhysicalVolume = new G4PVPlacement    198     detectorPhysicalVolume = new G4PVPlacement(0,
199                                                   199                                                detectorPosition, // Setted by displacement
200                                                   200                                                "DetectorPhys",
201                                                   201                                                detectorLogicalVolume,
202                                                   202                                                phantomPhysicalVolume,
203                                                   203                                                false,0);
204                                                   204     
205     // Visualisation attributes of the detecto    205     // Visualisation attributes of the detector
206     skyBlue = new G4VisAttributes( G4Colour(13    206     skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
207     skyBlue -> SetVisibility(true);               207     skyBlue -> SetVisibility(true);
208     skyBlue -> SetForceSolid(true);               208     skyBlue -> SetForceSolid(true);
209     //skyBlue -> SetForceWireframe(true);         209     //skyBlue -> SetForceWireframe(true);
210     detectorLogicalVolume -> SetVisAttributes(    210     detectorLogicalVolume -> SetVisAttributes(skyBlue);
211                                                   211     
212     // **************                             212     // **************
213     // **************                             213     // **************
214     // Cut per Region                             214     // Cut per Region
215     // **************                             215     // **************
216     //                                            216     //
217     // A smaller cut is fixed in the phantom t    217     // A smaller cut is fixed in the phantom to calculate the energy deposit with the
218     // required accuracy                          218     // required accuracy
219     if (!aRegion)                                 219     if (!aRegion)
220     {                                             220     {
221         aRegion = new G4Region("DetectorLog");    221         aRegion = new G4Region("DetectorLog");
222         detectorLogicalVolume -> SetRegion(aRe    222         detectorLogicalVolume -> SetRegion(aRegion);
223         aRegion->AddRootLogicalVolume( detecto    223         aRegion->AddRootLogicalVolume( detectorLogicalVolume );
224     }                                             224     }
225 }                                                 225 }
226                                                   226 
227 //////////////////////////////////////////////    227 ///////////////////////////////////////////////////////////////////////
228 void HadrontherapyDetectorConstruction::Initia    228 void HadrontherapyDetectorConstruction::InitializeDetectorROGeometry(
229                                                   229                                                                      HadrontherapyDetectorROGeometry* RO,
230                                                   230                                                                      G4ThreeVector
231                                                   231                                                                      detectorToWorldPosition)
232 {                                                 232 {
233     RO->Initialize(detectorToWorldPosition,       233     RO->Initialize(detectorToWorldPosition,
234                    detectorSizeX/2,               234                    detectorSizeX/2,
235                    detectorSizeY/2,               235                    detectorSizeY/2,
236                    detectorSizeZ/2,               236                    detectorSizeZ/2,
237                    numberOfVoxelsAlongX,          237                    numberOfVoxelsAlongX,
238                    numberOfVoxelsAlongY,          238                    numberOfVoxelsAlongY,
239                    numberOfVoxelsAlongZ);         239                    numberOfVoxelsAlongZ);
240 }                                                 240 }
241 void HadrontherapyDetectorConstruction::Virtua    241 void HadrontherapyDetectorConstruction::VirtualLayer(G4bool Varbool)
242 {                                                 242 {
243                                                   243    
244     //Virtual  plane                              244     //Virtual  plane
245     VirtualLayerPosition = G4ThreeVector(0*cm,    245     VirtualLayerPosition = G4ThreeVector(0*cm,0*cm,0*cm);
246     NewSource= Varbool;                           246     NewSource= Varbool;
247     if(NewSource == true)                         247     if(NewSource == true)
248     {                                             248     {
249        // std::cout<<"trr"<<std::endl;            249        // std::cout<<"trr"<<std::endl;
250         G4Material* airNist =  G4NistManager::    250         G4Material* airNist =  G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
251                                                   251         
252         solidVirtualLayer = new G4Box("Virtual    252         solidVirtualLayer = new G4Box("VirtualLayer",
253                                       1.*um,      253                                       1.*um,
254                                       20.*cm,     254                                       20.*cm,
255                                       40.*cm);    255                                       40.*cm);
256                                                   256         
257         logicVirtualLayer = new G4LogicalVolum    257         logicVirtualLayer = new G4LogicalVolume(
258                                                   258                                                 solidVirtualLayer,
259                                                   259                                                 airNist,
260                                                   260                                                 "VirtualLayer");
261                                                   261         
262         physVirtualLayer= new G4PVPlacement(0,    262         physVirtualLayer= new G4PVPlacement(0,VirtualLayerPosition,
263                                             "V    263                                             "VirtualLayer",
264                                             lo    264                                             logicVirtualLayer,
265                                             mo    265                                             motherPhys,
266                                             fa    266                                             false,
267                                             0)    267                                             0);
268                                                   268         
269         logicVirtualLayer -> SetVisAttributes(    269         logicVirtualLayer -> SetVisAttributes(skyBlue);
270     }                                             270     }
271                                                   271     
272                                                   272     
273                                                   273     
274                                                   274     
275 }                                                 275 }
276                                                   276 
277                                                   277 
278 //////////////////////////////////////////////    278 ///////////////////////////////////////////////////////////////////////
279 void  HadrontherapyDetectorConstruction::Param    279 void  HadrontherapyDetectorConstruction::ParametersCheck()
280 {                                                 280 {
281     // Check phantom/detector sizes & relative    281     // Check phantom/detector sizes & relative position
282     if (!IsInside(detectorSizeX,                  282     if (!IsInside(detectorSizeX,
283                   detectorSizeY,                  283                   detectorSizeY,
284                   detectorSizeZ,                  284                   detectorSizeZ,
285                   phantomSizeX,                   285                   phantomSizeX,
286                   phantomSizeY,                   286                   phantomSizeY,
287                   phantomSizeZ,                   287                   phantomSizeZ,
288                   detectorToPhantomPosition       288                   detectorToPhantomPosition
289                   ))                              289                   ))
290         G4Exception("HadrontherapyDetectorCons    290         G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
291                                                   291     
292     // Check Detector sizes respect to the vox    292     // Check Detector sizes respect to the voxel ones
293                                                   293     
294     if ( detectorSizeX < sizeOfVoxelAlongX) {     294     if ( detectorSizeX < sizeOfVoxelAlongX) {
295         G4Exception("HadrontherapyDetectorCons    295         G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error:  Detector X size must be bigger or equal than that of Voxel X!");
296     }                                             296     }
297     if ( detectorSizeY < sizeOfVoxelAlongY) {     297     if ( detectorSizeY < sizeOfVoxelAlongY) {
298         G4Exception(" HadrontherapyDetectorCon    298         G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error:  Detector Y size must be bigger or equal than that of Voxel Y!");
299     }                                             299     }
300     if ( detectorSizeZ < sizeOfVoxelAlongZ) {     300     if ( detectorSizeZ < sizeOfVoxelAlongZ) {
301         G4Exception(" HadrontherapyDetectorCon    301         G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error:  Detector Z size must be bigger or equal than that of Voxel Z!");
302     }                                             302     }
303 }                                                 303 }
304                                                   304 
305 //////////////////////////////////////////////    305 ///////////////////////////////////////////////////////////////////////
306 G4bool HadrontherapyDetectorConstruction::SetP    306 G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
307 {                                                 307 {
308                                                   308     
309     if (G4Material* pMat = G4NistManager::Inst    309     if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
310     {                                             310     {
311         phantomMaterial  = pMat;                  311         phantomMaterial  = pMat;
312         detectorMaterial = pMat;                  312         detectorMaterial = pMat;
313         if (detectorLogicalVolume && phantomLo    313         if (detectorLogicalVolume && phantomLogicalVolume)
314         {                                         314         {
315             detectorLogicalVolume -> SetMateri    315             detectorLogicalVolume -> SetMaterial(pMat);
316             phantomLogicalVolume ->  SetMateri    316             phantomLogicalVolume ->  SetMaterial(pMat);
317                                                   317             
318             G4RunManager::GetRunManager() -> P    318             G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
319             G4RunManager::GetRunManager() -> G    319             G4RunManager::GetRunManager() -> GeometryHasBeenModified();
320             G4cout << "The material of Phantom    320             G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
321         }                                         321         }
322     }                                             322     }
323     else                                          323     else
324     {                                             324     {
325         G4cout << "WARNING: material \"" << ma    325         G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
326         " table [located in $G4INSTALL/source/    326         " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
327         G4cout << "Use command \"/parameter/ni    327         G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
328         return false;                             328         return false;
329     }                                             329     }
330                                                   330     
331     return true;                                  331     return true;
332 }                                                 332 }
333 //////////////////////////////////////////////    333 /////////////////////////////////////////////////////////////////////////////
334 void HadrontherapyDetectorConstruction::SetPha    334 void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
335 {                                                 335 {
336     if (sizeX > 0.) phantomSizeX = sizeX;         336     if (sizeX > 0.) phantomSizeX = sizeX;
337     if (sizeY > 0.) phantomSizeY = sizeY;         337     if (sizeY > 0.) phantomSizeY = sizeY;
338     if (sizeZ > 0.) phantomSizeZ = sizeZ;         338     if (sizeZ > 0.) phantomSizeZ = sizeZ;
339 }                                                 339 }
340                                                   340 
341 //////////////////////////////////////////////    341 /////////////////////////////////////////////////////////////////////////////
342 void HadrontherapyDetectorConstruction::SetDet    342 void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
343 {                                                 343 {
344     if (sizeX > 0.) {detectorSizeX = sizeX;}      344     if (sizeX > 0.) {detectorSizeX = sizeX;}
345     if (sizeY > 0.) {detectorSizeY = sizeY;}      345     if (sizeY > 0.) {detectorSizeY = sizeY;}
346     if (sizeZ > 0.) {detectorSizeZ = sizeZ;}      346     if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
347     SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxe    347     SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
348 }                                                 348 }
349                                                   349 
350 //////////////////////////////////////////////    350 /////////////////////////////////////////////////////////////////////////////
351 void HadrontherapyDetectorConstruction::SetVox    351 void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
352 {                                                 352 {
353     if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX    353     if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
354     if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY    354     if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
355     if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ    355     if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
356 }                                                 356 }
357                                                   357 
358 //////////////////////////////////////////////    358 ///////////////////////////////////////////////////////////////////////
359 void HadrontherapyDetectorConstruction::SetPha    359 void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
360 {                                                 360 {
361     phantomPosition = pos;                        361     phantomPosition = pos;
362 }                                                 362 }
363                                                   363 
364 //////////////////////////////////////////////    364 /////////////////////////////////////////////////////////////////////////////
365 void HadrontherapyDetectorConstruction::SetDet    365 void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
366 {                                                 366 {
367     detectorToPhantomPosition = displ;            367     detectorToPhantomPosition = displ;
368 }                                                 368 }
369                                                   369 
370 void HadrontherapyDetectorConstruction::SetVir    370 void HadrontherapyDetectorConstruction::SetVirtualLayerPosition(G4ThreeVector position)
371 {                                                 371 {
372                                                   372     
373     VirtualLayerPosition = position;              373     VirtualLayerPosition = position;
374     physVirtualLayer->SetTranslation(VirtualLa    374     physVirtualLayer->SetTranslation(VirtualLayerPosition);
375                                                   375     
376 }                                                 376 }
377 //////////////////////////////////////////////    377 /////////////////////////////////////////////////////////////////////////////
378 void HadrontherapyDetectorConstruction::Update    378 void HadrontherapyDetectorConstruction::UpdateGeometry()
379 {                                                 379 {
380     /*                                            380     /*
381      * Check parameters consistency               381      * Check parameters consistency
382      */                                           382      */
383     ParametersCheck();                            383     ParametersCheck();
384                                                   384     
385     G4GeometryManager::GetInstance() -> OpenGe    385     G4GeometryManager::GetInstance() -> OpenGeometry();
386     if (phantom)                                  386     if (phantom)
387     {                                             387     {
388         phantom -> SetXHalfLength(phantomSizeX    388         phantom -> SetXHalfLength(phantomSizeX/2);
389         phantom -> SetYHalfLength(phantomSizeY    389         phantom -> SetYHalfLength(phantomSizeY/2);
390         phantom -> SetZHalfLength(phantomSizeZ    390         phantom -> SetZHalfLength(phantomSizeZ/2);
391                                                   391         
392         phantomPhysicalVolume -> SetTranslatio    392         phantomPhysicalVolume -> SetTranslation(phantomPosition);
393     }                                             393     }
394     else   ConstructPhantom();                    394     else   ConstructPhantom();
395                                                   395     
396                                                   396     
397     // Get the center of the detector             397     // Get the center of the detector
398     SetDetectorPosition();                        398     SetDetectorPosition();
399     if (detector)                                 399     if (detector)
400     {                                             400     {
401                                                   401         
402         detector -> SetXHalfLength(detectorSiz    402         detector -> SetXHalfLength(detectorSizeX/2);
403         detector -> SetYHalfLength(detectorSiz    403         detector -> SetYHalfLength(detectorSizeY/2);
404         detector -> SetZHalfLength(detectorSiz    404         detector -> SetZHalfLength(detectorSizeZ/2);
405                                                   405         
406         detectorPhysicalVolume -> SetTranslati    406         detectorPhysicalVolume -> SetTranslation(detectorPosition);
407     }                                             407     }
408     else    ConstructDetector();                  408     else    ConstructDetector();
409                                                   409     
410     //std::cout<<NewSource<<std::endl;            410     //std::cout<<NewSource<<std::endl;
411     /*if(NewSource)                               411     /*if(NewSource)
412      {                                            412      {
413      std::cout<<"via"<<std::endl;                 413      std::cout<<"via"<<std::endl;
414      }*/                                          414      }*/
415                                                   415     
416                                                   416     
417     // std::cout<<"i"<<std::endl;                 417     // std::cout<<"i"<<std::endl;
418     // std::cout<<VirtualLayerPosition<<std::e    418     // std::cout<<VirtualLayerPosition<<std::endl;
419     // physVirtualLayer->SetTranslation(Virtua    419     // physVirtualLayer->SetTranslation(VirtualLayerPosition);
420                                                   420     
421                                                   421     
422                                                   422     
423                                                   423     
424                                                   424     
425     // Round to nearest integer number of voxe    425     // Round to nearest integer number of voxel
426                                                   426     
427     numberOfVoxelsAlongX = G4lrint(detectorSiz    427     numberOfVoxelsAlongX = G4lrint(detectorSizeX / sizeOfVoxelAlongX);
428     sizeOfVoxelAlongX = ( detectorSizeX / numb    428     sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
429     numberOfVoxelsAlongY = G4lrint(detectorSiz    429     numberOfVoxelsAlongY = G4lrint(detectorSizeY / sizeOfVoxelAlongY);
430     sizeOfVoxelAlongY = ( detectorSizeY / numb    430     sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
431     numberOfVoxelsAlongZ = G4lrint(detectorSiz    431     numberOfVoxelsAlongZ = G4lrint(detectorSizeZ / sizeOfVoxelAlongZ);
432     sizeOfVoxelAlongZ = ( detectorSizeZ / numb    432     sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
433     PassiveProtonBeamLine *ppbl= (PassiveProto    433     PassiveProtonBeamLine *ppbl= (PassiveProtonBeamLine*)
434                                                   434     
435     G4RunManager::GetRunManager()->GetUserDete    435     G4RunManager::GetRunManager()->GetUserDetectorConstruction();
436                                                   436     
437     HadrontherapyDetectorROGeometry* RO = (Had    437     HadrontherapyDetectorROGeometry* RO = (HadrontherapyDetectorROGeometry*) ppbl->GetParallelWorld(0);
438                                                   438     
439     //Set parameters, either for the Construct    439     //Set parameters, either for the Construct() or for the UpdateROGeometry()
440     RO->Initialize(GetDetectorToWorldPosition(    440     RO->Initialize(GetDetectorToWorldPosition(),
441                    detectorSizeX/2,               441                    detectorSizeX/2,
442                    detectorSizeY/2,               442                    detectorSizeY/2,
443                    detectorSizeZ/2,               443                    detectorSizeZ/2,
444                    numberOfVoxelsAlongX,          444                    numberOfVoxelsAlongX,
445                    numberOfVoxelsAlongY,          445                    numberOfVoxelsAlongY,
446                    numberOfVoxelsAlongZ);         446                    numberOfVoxelsAlongZ);
447                                                   447     
448     //This method below has an effect only if     448     //This method below has an effect only if the RO geometry is already built.
449     RO->UpdateROGeometry();                       449     RO->UpdateROGeometry();
450                                                   450     
451                                                   451     
452                                                   452     
453     volumeOfVoxel = sizeOfVoxelAlongX * sizeOf    453     volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
454     massOfVoxel = detectorMaterial -> GetDensi    454     massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
455     //  This will clear the existing matrix (t    455     //  This will clear the existing matrix (together with all data inside it)!
456     matrix = HadrontherapyMatrix::GetInstance(    456     matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
457                                                   457                                               numberOfVoxelsAlongY,
458                                                   458                                               numberOfVoxelsAlongZ,
459                                                   459                                               massOfVoxel);
460                                                   460     
461                                                   461     
462     // Initialize RBE                             462     // Initialize RBE
463     HadrontherapyRBE::CreateInstance(numberOfV    463     HadrontherapyRBE::CreateInstance(numberOfVoxelsAlongX, numberOfVoxelsAlongY, numberOfVoxelsAlongZ, massOfVoxel);
464                                                   464     
465     // Comment out the line below if let calcu    465     // Comment out the line below if let calculation is not needed!
466     // Initialize LET with energy of primaries    466     // Initialize LET with energy of primaries and clear data inside
467     if ( (let = HadrontherapyLet::GetInstance(    467     if ( (let = HadrontherapyLet::GetInstance(this)) )
468     {                                             468     {
469         HadrontherapyLet::GetInstance() -> Ini    469         HadrontherapyLet::GetInstance() -> Initialize();
470     }                                             470     }
471                                                   471     
472                                                   472     
473     // Initialize analysis                        473     // Initialize analysis
474     // Inform the kernel about the new geometr    474     // Inform the kernel about the new geometry
475     G4RunManager::GetRunManager() -> GeometryH    475     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
476     G4RunManager::GetRunManager() -> PhysicsHa    476     G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
477                                                   477     
478     PrintParameters();                            478     PrintParameters();
479                                                   479     
480     // CheckOverlaps();                           480     // CheckOverlaps();
481 }                                                 481 }
482                                                   482 
483 //////////////////////////////////////////////    483 /////////////////////////////////////////////////////////////////////////////
484 //Check of the geometry                           484 //Check of the geometry
485 //////////////////////////////////////////////    485 /////////////////////////////////////////////////////////////////////////////
486 void HadrontherapyDetectorConstruction::CheckO    486 void HadrontherapyDetectorConstruction::CheckOverlaps()
487 {                                                 487 {
488     G4PhysicalVolumeStore* thePVStore = G4Phys    488     G4PhysicalVolumeStore* thePVStore = G4PhysicalVolumeStore::GetInstance();
489     G4cout << thePVStore->size() << " physical    489     G4cout << thePVStore->size() << " physical volumes are defined" << G4endl;
490     G4bool overlapFlag = false;                   490     G4bool overlapFlag = false;
491     G4int res=1000;                               491     G4int res=1000;
492     G4double tol=0.; //tolerance                  492     G4double tol=0.; //tolerance
493     for (size_t i=0;i<thePVStore->size();i++)     493     for (size_t i=0;i<thePVStore->size();i++)
494     {                                             494     {
495         //overlapFlag = (*thePVStore)[i]->Chec    495         //overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,fCheckOverlapsVerbosity) | overlapFlag;
496         overlapFlag = (*thePVStore)[i]->CheckO    496         overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,true) | overlapFlag;    }
497     if (overlapFlag)                              497     if (overlapFlag)
498         G4cout << "Check: there are overlappin    498         G4cout << "Check: there are overlapping volumes" << G4endl;
499 }                                                 499 }
500                                                   500 
501 //////////////////////////////////////////////    501 /////////////////////////////////////////////////////////////////////////////
502 void HadrontherapyDetectorConstruction::PrintP    502 void HadrontherapyDetectorConstruction::PrintParameters()
503 {                                                 503 {
504                                                   504     
505     G4cout << "The (X,Y,Z) dimensions of the p    505     G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
506     G4BestUnit( phantom -> GetXHalfLength()*2.    506     G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
507     G4BestUnit( phantom -> GetYHalfLength()*2.    507     G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
508     G4BestUnit( phantom -> GetZHalfLength()*2.    508     G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
509                                                   509     
510     G4cout << "The (X,Y,Z) dimensions of the d    510     G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
511     G4BestUnit( detector -> GetXHalfLength()*2    511     G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
512     G4BestUnit( detector -> GetYHalfLength()*2    512     G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
513     G4BestUnit( detector -> GetZHalfLength()*2    513     G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
514                                                   514     
515     G4cout << "Displacement between Phantom an    515     G4cout << "Displacement between Phantom and World is: ";
516     G4cout << "DX= "<< G4BestUnit(phantomPosit    516     G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
517     "DY= "<< G4BestUnit(phantomPosition.getY()    517     "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
518     "DZ= "<< G4BestUnit(phantomPosition.getZ()    518     "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
519                                                   519     
520     G4cout << "The (X,Y,Z) sizes of the Voxels    520     G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
521     G4BestUnit(sizeOfVoxelAlongX, "Length")  <    521     G4BestUnit(sizeOfVoxelAlongX, "Length")  << ',' <<
522     G4BestUnit(sizeOfVoxelAlongY, "Length")  <    522     G4BestUnit(sizeOfVoxelAlongY, "Length")  << ',' <<
523     G4BestUnit(sizeOfVoxelAlongZ, "Length") <<    523     G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
524                                                   524     
525     G4cout << "The number of Voxels along (X,Y    525     G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
526     numberOfVoxelsAlongX  << ',' <<               526     numberOfVoxelsAlongX  << ',' <<
527     numberOfVoxelsAlongY  <<','  <<               527     numberOfVoxelsAlongY  <<','  <<
528     numberOfVoxelsAlongZ  << ')' << G4endl;       528     numberOfVoxelsAlongZ  << ')' << G4endl;
529 }                                                 529 }
530                                                   530 
531                                                   531 
532                                                   532