Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /geant4/examples/advanced/hadrontherapy/src/HadrontherapyDetectorConstruction.cc was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Hadrontherapy advanced example for Geant4
0027 // See more at: https://twiki.cern.ch/twiki/bin/view/Geant4/AdvancedExamplesHadrontherapy
0028 
0029 #include "G4UnitsTable.hh"
0030 #include "G4SDManager.hh"
0031 #include "G4RunManager.hh"
0032 #include "G4GeometryManager.hh"
0033 #include "G4SolidStore.hh"
0034 #include "G4PhysicalVolumeStore.hh"
0035 #include "G4LogicalVolumeStore.hh"
0036 #include "G4Box.hh"
0037 #include "G4LogicalVolume.hh"
0038 #include "G4ThreeVector.hh"
0039 #include "G4PVPlacement.hh"
0040 #include "globals.hh"
0041 #include "G4Transform3D.hh"
0042 #include "G4RotationMatrix.hh"
0043 #include "G4Colour.hh"
0044 #include "G4UserLimits.hh"
0045 #include "G4UnitsTable.hh"
0046 #include "G4VisAttributes.hh"
0047 #include "G4NistManager.hh"
0048 #include "HadrontherapyDetectorConstruction.hh"
0049 #include "HadrontherapyDetectorROGeometry.hh"
0050 #include "HadrontherapyDetectorMessenger.hh"
0051 #include "HadrontherapyDetectorSD.hh"
0052 #include "HadrontherapyMatrix.hh"
0053 #include "HadrontherapyLet.hh"
0054 #include "PassiveProtonBeamLine.hh"
0055 #include "BESTPassiveProtonBeamLine.hh"
0056 #include "HadrontherapyMatrix.hh"
0057 
0058 #include "HadrontherapyRBE.hh"
0059 #include "G4SystemOfUnits.hh"
0060 
0061 #include <cmath>
0062 
0063 
0064 
0065 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::instance = 0;
0066 /////////////////////////////////////////////////////////////////////////////
0067 HadrontherapyDetectorConstruction::HadrontherapyDetectorConstruction(G4VPhysicalVolume* physicalTreatmentRoom)
0068 : motherPhys(physicalTreatmentRoom), // pointer to WORLD volume
0069 detectorSD(0), detectorROGeometry(0), matrix(0),
0070 phantom(0), detector(0),
0071 phantomLogicalVolume(0), detectorLogicalVolume(0),
0072 phantomPhysicalVolume(0), detectorPhysicalVolume(0),
0073 aRegion(0)
0074 {
0075     
0076     
0077     /* NOTE! that the HadrontherapyDetectorConstruction class
0078      * does NOT inherit from G4VUserDetectorConstruction G4 class
0079      * So the Construct() mandatory virtual method is inside another geometric class
0080      * like the passiveProtonBeamLIne, ...
0081      */
0082     
0083     // Messenger to change parameters of the phantom/detector geometry
0084     detectorMessenger = new HadrontherapyDetectorMessenger(this);
0085     
0086     // Default detector voxels size
0087     // 200 slabs along the beam direction (X)
0088     sizeOfVoxelAlongX = 200 *um;
0089     sizeOfVoxelAlongY = 4 *cm;
0090     sizeOfVoxelAlongZ = 4 *cm;
0091     
0092     // Define here the material of the water phantom and of the detector
0093     SetPhantomMaterial("G4_WATER");
0094     // Construct geometry (messenger commands)
0095     // SetDetectorSize(4.*cm, 4.*cm, 4.*cm);
0096     SetDetectorSize(4. *cm, 4. *cm, 4. *cm);
0097     SetPhantomSize(40. *cm, 40. *cm, 40. *cm);
0098     
0099     SetPhantomPosition(G4ThreeVector(20. *cm, 0. *cm, 0. *cm));
0100     SetDetectorToPhantomPosition(G4ThreeVector(0. *cm, 0. *cm, 0. *cm));
0101     SetDetectorPosition();
0102     //GetDetectorToWorldPosition();
0103     
0104     // Write virtual parameters to the real ones and check for consistency
0105     UpdateGeometry();
0106     
0107     
0108     
0109 }
0110 
0111 /////////////////////////////////////////////////////////////////////////////
0112 HadrontherapyDetectorConstruction::~HadrontherapyDetectorConstruction()
0113 {
0114     delete detectorROGeometry;
0115     delete matrix;
0116     delete detectorMessenger;
0117 }
0118 
0119 /////////////////////////////////////////////////////////////////////////////
0120 HadrontherapyDetectorConstruction* HadrontherapyDetectorConstruction::GetInstance()
0121 {
0122     return instance;
0123 }
0124 
0125 /////////////////////////////////////////////////////////////////////////////
0126 // ConstructPhantom() is the method that construct a water box (called phantom
0127 // (or water phantom) in the usual Medical physicists slang).
0128 // A water phantom can be considered a good approximation of a an human body.
0129 void HadrontherapyDetectorConstruction::ConstructPhantom()
0130 {
0131     // Definition of the solid volume of the Phantom
0132     phantom = new G4Box("Phantom",
0133                         phantomSizeX/2,
0134                         phantomSizeY/2,
0135                         phantomSizeZ/2);
0136     
0137     // Definition of the logical volume of the Phantom
0138     phantomLogicalVolume = new G4LogicalVolume(phantom,
0139                                                phantomMaterial,
0140                                                "phantomLog", 0, 0, 0);
0141     
0142     // Definition of the physics volume of the Phantom
0143     phantomPhysicalVolume = new G4PVPlacement(0,
0144                                               phantomPosition,
0145                                               "phantomPhys",
0146                                               phantomLogicalVolume,
0147                                               motherPhys,
0148                                               false,
0149                                               0);
0150     
0151     // Visualisation attributes of the phantom
0152     red = new G4VisAttributes(G4Colour(255/255., 0/255. ,0/255.));
0153     red -> SetVisibility(true);
0154     red -> SetForceSolid(true);
0155     red -> SetForceWireframe(true);
0156     phantomLogicalVolume -> SetVisAttributes(red);
0157 }
0158 
0159 /////////////////////////////////////////////////////////////////////////////
0160 // ConstructDetector() is the method the reconstruct a detector region
0161 // inside the water phantom. It is a volume, located inside the water phantom.
0162 //
0163 //           **************************
0164 //           *    water phantom       *
0165 //           *                        *
0166 //           *                        *
0167 //           *---------------         *
0168 //  Beam     *              -         *
0169 //  ----->   * detector     -         *
0170 //           *              -         *
0171 //           *---------------         *
0172 //           *                        *
0173 //           *                        *
0174 //           *                        *
0175 //           **************************
0176 //
0177 // The detector can be dived in slices or voxelized
0178 // and inside it different quantities (dose distribution, fluence distribution, LET, etc)
0179 // can be stored.
0180 void HadrontherapyDetectorConstruction::ConstructDetector()
0181 
0182 {
0183     // Definition of the solid volume of the Detector
0184     detector = new G4Box("Detector",
0185                          
0186                          phantomSizeX/2,
0187                          
0188                          phantomSizeY/2,
0189                          
0190                          phantomSizeZ/2);
0191     
0192     // Definition of the logic volume of the Phantom
0193     detectorLogicalVolume = new G4LogicalVolume(detector,
0194                                                 detectorMaterial,
0195                                                 "DetectorLog",
0196                                                 0,0,0);
0197     // Definition of the physical volume of the Phantom
0198     detectorPhysicalVolume = new G4PVPlacement(0,
0199                                                detectorPosition, // Setted by displacement
0200                                                "DetectorPhys",
0201                                                detectorLogicalVolume,
0202                                                phantomPhysicalVolume,
0203                                                false,0);
0204     
0205     // Visualisation attributes of the detector
0206     skyBlue = new G4VisAttributes( G4Colour(135/255. , 206/255. ,  235/255. ));
0207     skyBlue -> SetVisibility(true);
0208     skyBlue -> SetForceSolid(true);
0209     //skyBlue -> SetForceWireframe(true);
0210     detectorLogicalVolume -> SetVisAttributes(skyBlue);
0211     
0212     // **************
0213     // **************
0214     // Cut per Region
0215     // **************
0216     //
0217     // A smaller cut is fixed in the phantom to calculate the energy deposit with the
0218     // required accuracy
0219     if (!aRegion)
0220     {
0221         aRegion = new G4Region("DetectorLog");
0222         detectorLogicalVolume -> SetRegion(aRegion);
0223         aRegion->AddRootLogicalVolume( detectorLogicalVolume );
0224     }
0225 }
0226 
0227 ///////////////////////////////////////////////////////////////////////
0228 void HadrontherapyDetectorConstruction::InitializeDetectorROGeometry(
0229                                                                      HadrontherapyDetectorROGeometry* RO,
0230                                                                      G4ThreeVector
0231                                                                      detectorToWorldPosition)
0232 {
0233     RO->Initialize(detectorToWorldPosition,
0234                    detectorSizeX/2,
0235                    detectorSizeY/2,
0236                    detectorSizeZ/2,
0237                    numberOfVoxelsAlongX,
0238                    numberOfVoxelsAlongY,
0239                    numberOfVoxelsAlongZ);
0240 }
0241 void HadrontherapyDetectorConstruction::VirtualLayer(G4bool Varbool)
0242 {
0243    
0244     //Virtual  plane
0245     VirtualLayerPosition = G4ThreeVector(0*cm,0*cm,0*cm);
0246     NewSource= Varbool;
0247     if(NewSource == true)
0248     {
0249        // std::cout<<"trr"<<std::endl;
0250         G4Material* airNist =  G4NistManager::Instance()->FindOrBuildMaterial("G4_AIR");
0251         
0252         solidVirtualLayer = new G4Box("VirtualLayer",
0253                                       1.*um,
0254                                       20.*cm,
0255                                       40.*cm);
0256         
0257         logicVirtualLayer = new G4LogicalVolume(
0258                                                 solidVirtualLayer,
0259                                                 airNist,
0260                                                 "VirtualLayer");
0261         
0262         physVirtualLayer= new G4PVPlacement(0,VirtualLayerPosition,
0263                                             "VirtualLayer",
0264                                             logicVirtualLayer,
0265                                             motherPhys,
0266                                             false,
0267                                             0);
0268         
0269         logicVirtualLayer -> SetVisAttributes(skyBlue);
0270     }
0271     
0272     
0273     
0274     
0275 }
0276 
0277 
0278 ///////////////////////////////////////////////////////////////////////
0279 void  HadrontherapyDetectorConstruction::ParametersCheck()
0280 {
0281     // Check phantom/detector sizes & relative position
0282     if (!IsInside(detectorSizeX,
0283                   detectorSizeY,
0284                   detectorSizeZ,
0285                   phantomSizeX,
0286                   phantomSizeY,
0287                   phantomSizeZ,
0288                   detectorToPhantomPosition
0289                   ))
0290         G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0001", FatalException, "Error: Detector is not fully inside Phantom!");
0291     
0292     // Check Detector sizes respect to the voxel ones
0293     
0294     if ( detectorSizeX < sizeOfVoxelAlongX) {
0295         G4Exception("HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0002", FatalException, "Error:  Detector X size must be bigger or equal than that of Voxel X!");
0296     }
0297     if ( detectorSizeY < sizeOfVoxelAlongY) {
0298         G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0003", FatalException, "Error:  Detector Y size must be bigger or equal than that of Voxel Y!");
0299     }
0300     if ( detectorSizeZ < sizeOfVoxelAlongZ) {
0301         G4Exception(" HadrontherapyDetectorConstruction::ParametersCheck()", "Hadrontherapy0004", FatalException, "Error:  Detector Z size must be bigger or equal than that of Voxel Z!");
0302     }
0303 }
0304 
0305 ///////////////////////////////////////////////////////////////////////
0306 G4bool HadrontherapyDetectorConstruction::SetPhantomMaterial(G4String material)
0307 {
0308     
0309     if (G4Material* pMat = G4NistManager::Instance()->FindOrBuildMaterial(material, false) )
0310     {
0311         phantomMaterial  = pMat;
0312         detectorMaterial = pMat;
0313         if (detectorLogicalVolume && phantomLogicalVolume)
0314         {
0315             detectorLogicalVolume -> SetMaterial(pMat);
0316             phantomLogicalVolume ->  SetMaterial(pMat);
0317             
0318             G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0319             G4RunManager::GetRunManager() -> GeometryHasBeenModified();
0320             G4cout << "The material of Phantom/Detector has been changed to " << material << G4endl;
0321         }
0322     }
0323     else
0324     {
0325         G4cout << "WARNING: material \"" << material << "\" doesn't exist in NIST elements/materials"
0326         " table [located in $G4INSTALL/source/materials/src/G4NistMaterialBuilder.cc]" << G4endl;
0327         G4cout << "Use command \"/parameter/nist\" to see full materials list!" << G4endl;
0328         return false;
0329     }
0330     
0331     return true;
0332 }
0333 /////////////////////////////////////////////////////////////////////////////
0334 void HadrontherapyDetectorConstruction::SetPhantomSize(G4double sizeX, G4double sizeY, G4double sizeZ)
0335 {
0336     if (sizeX > 0.) phantomSizeX = sizeX;
0337     if (sizeY > 0.) phantomSizeY = sizeY;
0338     if (sizeZ > 0.) phantomSizeZ = sizeZ;
0339 }
0340 
0341 /////////////////////////////////////////////////////////////////////////////
0342 void HadrontherapyDetectorConstruction::SetDetectorSize(G4double sizeX, G4double sizeY, G4double sizeZ)
0343 {
0344     if (sizeX > 0.) {detectorSizeX = sizeX;}
0345     if (sizeY > 0.) {detectorSizeY = sizeY;}
0346     if (sizeZ > 0.) {detectorSizeZ = sizeZ;}
0347     SetVoxelSize(sizeOfVoxelAlongX, sizeOfVoxelAlongY, sizeOfVoxelAlongZ);
0348 }
0349 
0350 /////////////////////////////////////////////////////////////////////////////
0351 void HadrontherapyDetectorConstruction::SetVoxelSize(G4double sizeX, G4double sizeY, G4double sizeZ)
0352 {
0353     if (sizeX > 0.) {sizeOfVoxelAlongX = sizeX;}
0354     if (sizeY > 0.) {sizeOfVoxelAlongY = sizeY;}
0355     if (sizeZ > 0.) {sizeOfVoxelAlongZ = sizeZ;}
0356 }
0357 
0358 ///////////////////////////////////////////////////////////////////////
0359 void HadrontherapyDetectorConstruction::SetPhantomPosition(G4ThreeVector pos)
0360 {
0361     phantomPosition = pos;
0362 }
0363 
0364 /////////////////////////////////////////////////////////////////////////////
0365 void HadrontherapyDetectorConstruction::SetDetectorToPhantomPosition(G4ThreeVector displ)
0366 {
0367     detectorToPhantomPosition = displ;
0368 }
0369 
0370 void HadrontherapyDetectorConstruction::SetVirtualLayerPosition(G4ThreeVector position)
0371 {
0372     
0373     VirtualLayerPosition = position;
0374     physVirtualLayer->SetTranslation(VirtualLayerPosition);
0375     
0376 }
0377 /////////////////////////////////////////////////////////////////////////////
0378 void HadrontherapyDetectorConstruction::UpdateGeometry()
0379 {
0380     /*
0381      * Check parameters consistency
0382      */
0383     ParametersCheck();
0384     
0385     G4GeometryManager::GetInstance() -> OpenGeometry();
0386     if (phantom)
0387     {
0388         phantom -> SetXHalfLength(phantomSizeX/2);
0389         phantom -> SetYHalfLength(phantomSizeY/2);
0390         phantom -> SetZHalfLength(phantomSizeZ/2);
0391         
0392         phantomPhysicalVolume -> SetTranslation(phantomPosition);
0393     }
0394     else   ConstructPhantom();
0395     
0396     
0397     // Get the center of the detector
0398     SetDetectorPosition();
0399     if (detector)
0400     {
0401         
0402         detector -> SetXHalfLength(detectorSizeX/2);
0403         detector -> SetYHalfLength(detectorSizeY/2);
0404         detector -> SetZHalfLength(detectorSizeZ/2);
0405         
0406         detectorPhysicalVolume -> SetTranslation(detectorPosition);
0407     }
0408     else    ConstructDetector();
0409     
0410     //std::cout<<NewSource<<std::endl;
0411     /*if(NewSource)
0412      {
0413      std::cout<<"via"<<std::endl;
0414      }*/
0415     
0416     
0417     // std::cout<<"i"<<std::endl;
0418     // std::cout<<VirtualLayerPosition<<std::endl;
0419     // physVirtualLayer->SetTranslation(VirtualLayerPosition);
0420     
0421     
0422     
0423     
0424     
0425     // Round to nearest integer number of voxel
0426     
0427     numberOfVoxelsAlongX = G4lrint(detectorSizeX / sizeOfVoxelAlongX);
0428     sizeOfVoxelAlongX = ( detectorSizeX / numberOfVoxelsAlongX );
0429     numberOfVoxelsAlongY = G4lrint(detectorSizeY / sizeOfVoxelAlongY);
0430     sizeOfVoxelAlongY = ( detectorSizeY / numberOfVoxelsAlongY );
0431     numberOfVoxelsAlongZ = G4lrint(detectorSizeZ / sizeOfVoxelAlongZ);
0432     sizeOfVoxelAlongZ = ( detectorSizeZ / numberOfVoxelsAlongZ );
0433     PassiveProtonBeamLine *ppbl= (PassiveProtonBeamLine*)
0434     
0435     G4RunManager::GetRunManager()->GetUserDetectorConstruction();
0436     
0437     HadrontherapyDetectorROGeometry* RO = (HadrontherapyDetectorROGeometry*) ppbl->GetParallelWorld(0);
0438     
0439     //Set parameters, either for the Construct() or for the UpdateROGeometry()
0440     RO->Initialize(GetDetectorToWorldPosition(),
0441                    detectorSizeX/2,
0442                    detectorSizeY/2,
0443                    detectorSizeZ/2,
0444                    numberOfVoxelsAlongX,
0445                    numberOfVoxelsAlongY,
0446                    numberOfVoxelsAlongZ);
0447     
0448     //This method below has an effect only if the RO geometry is already built.
0449     RO->UpdateROGeometry();
0450     
0451     
0452     
0453     volumeOfVoxel = sizeOfVoxelAlongX * sizeOfVoxelAlongY * sizeOfVoxelAlongZ;
0454     massOfVoxel = detectorMaterial -> GetDensity() * volumeOfVoxel;
0455     //  This will clear the existing matrix (together with all data inside it)!
0456     matrix = HadrontherapyMatrix::GetInstance(numberOfVoxelsAlongX,
0457                                               numberOfVoxelsAlongY,
0458                                               numberOfVoxelsAlongZ,
0459                                               massOfVoxel);
0460     
0461     
0462     // Initialize RBE
0463     HadrontherapyRBE::CreateInstance(numberOfVoxelsAlongX, numberOfVoxelsAlongY, numberOfVoxelsAlongZ, massOfVoxel);
0464     
0465     // Comment out the line below if let calculation is not needed!
0466     // Initialize LET with energy of primaries and clear data inside
0467     if ( (let = HadrontherapyLet::GetInstance(this)) )
0468     {
0469         HadrontherapyLet::GetInstance() -> Initialize();
0470     }
0471     
0472     
0473     // Initialize analysis
0474     // Inform the kernel about the new geometry
0475     G4RunManager::GetRunManager() -> GeometryHasBeenModified();
0476     G4RunManager::GetRunManager() -> PhysicsHasBeenModified();
0477     
0478     PrintParameters();
0479     
0480     // CheckOverlaps();
0481 }
0482 
0483 /////////////////////////////////////////////////////////////////////////////
0484 //Check of the geometry
0485 /////////////////////////////////////////////////////////////////////////////
0486 void HadrontherapyDetectorConstruction::CheckOverlaps()
0487 {
0488     G4PhysicalVolumeStore* thePVStore = G4PhysicalVolumeStore::GetInstance();
0489     G4cout << thePVStore->size() << " physical volumes are defined" << G4endl;
0490     G4bool overlapFlag = false;
0491     G4int res=1000;
0492     G4double tol=0.; //tolerance
0493     for (size_t i=0;i<thePVStore->size();i++)
0494     {
0495         //overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,fCheckOverlapsVerbosity) | overlapFlag;
0496         overlapFlag = (*thePVStore)[i]->CheckOverlaps(res,tol,true) | overlapFlag;    }
0497     if (overlapFlag)
0498         G4cout << "Check: there are overlapping volumes" << G4endl;
0499 }
0500 
0501 /////////////////////////////////////////////////////////////////////////////
0502 void HadrontherapyDetectorConstruction::PrintParameters()
0503 {
0504     
0505     G4cout << "The (X,Y,Z) dimensions of the phantom are : (" <<
0506     G4BestUnit( phantom -> GetXHalfLength()*2., "Length") << ',' <<
0507     G4BestUnit( phantom -> GetYHalfLength()*2., "Length") << ',' <<
0508     G4BestUnit( phantom -> GetZHalfLength()*2., "Length") << ')' << G4endl;
0509     
0510     G4cout << "The (X,Y,Z) dimensions of the detector are : (" <<
0511     G4BestUnit( detector -> GetXHalfLength()*2., "Length") << ',' <<
0512     G4BestUnit( detector -> GetYHalfLength()*2., "Length") << ',' <<
0513     G4BestUnit( detector -> GetZHalfLength()*2., "Length") << ')' << G4endl;
0514     
0515     G4cout << "Displacement between Phantom and World is: ";
0516     G4cout << "DX= "<< G4BestUnit(phantomPosition.getX(),"Length") <<
0517     "DY= "<< G4BestUnit(phantomPosition.getY(),"Length") <<
0518     "DZ= "<< G4BestUnit(phantomPosition.getZ(),"Length") << G4endl;
0519     
0520     G4cout << "The (X,Y,Z) sizes of the Voxels are: (" <<
0521     G4BestUnit(sizeOfVoxelAlongX, "Length")  << ',' <<
0522     G4BestUnit(sizeOfVoxelAlongY, "Length")  << ',' <<
0523     G4BestUnit(sizeOfVoxelAlongZ, "Length") << ')' << G4endl;
0524     
0525     G4cout << "The number of Voxels along (X,Y,Z) is: (" <<
0526     numberOfVoxelsAlongX  << ',' <<
0527     numberOfVoxelsAlongY  <<','  <<
0528     numberOfVoxelsAlongZ  << ')' << G4endl;
0529 }
0530 
0531