Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:21:44 

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 //
0027 /// \file  medical/DICOM/src/DicomDetectorConstruction.cc
0028 /// \brief Implementation of the DicomDetectorConstruction class
0029 //
0030 
0031 #include "DicomDetectorConstruction.hh"
0032 
0033 #include "CLHEP/Units/SystemOfUnits.h"
0034 #include "DicomHandler.hh"
0035 #include "DicomPhantomZSliceHeader.hh"
0036 
0037 #include "G4Box.hh"
0038 #include "G4Element.hh"
0039 #include "G4LogicalVolume.hh"
0040 #include "G4Material.hh"
0041 #include "G4NistManager.hh"
0042 #include "G4PVPlacement.hh"
0043 #include "G4PhysicalConstants.hh"
0044 #include "G4UIcommand.hh"
0045 #include "G4VPhysicalVolume.hh"
0046 #include "globals.hh"
0047 
0048 #ifdef G4_DCMTK
0049 #  include "DicomFileMgr.hh"
0050 #endif
0051 #include "G4VisAttributes.hh"
0052 
0053 using CLHEP::cm3;
0054 using CLHEP::g;
0055 using CLHEP::m;
0056 using CLHEP::mg;
0057 using CLHEP::mole;
0058 using CLHEP::perCent;
0059 
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..
0061 DicomDetectorConstruction::DicomDetectorConstruction()
0062   : G4VUserDetectorConstruction(),
0063     fAir(0),
0064 
0065     fWorld_solid(0),
0066     fWorld_logic(0),
0067     fWorld_phys(0),
0068 
0069     fContainer_solid(0),
0070     fContainer_logic(0),
0071     fContainer_phys(0),
0072 
0073     fNoFiles(0),
0074     fMateIDs(0),
0075 
0076     fZSliceHeaderMerged(0),
0077 
0078     fNoVoxelsX(0),
0079     fNoVoxelsY(0),
0080     fNoVoxelsZ(0),
0081     fVoxelHalfDimX(0),
0082     fVoxelHalfDimY(0),
0083     fVoxelHalfDimZ(0),
0084 
0085     fConstructed(false)
0086 {}
0087 
0088 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
0089 DicomDetectorConstruction::~DicomDetectorConstruction() {}
0090 
0091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0092 G4VPhysicalVolume* DicomDetectorConstruction::Construct()
0093 {
0094   if (!fConstructed || fWorld_phys == 0) {
0095     fConstructed = true;
0096     InitialisationOfMaterials();
0097 
0098     //----- Build world
0099     G4double worldXDimension = 1. * m;
0100     G4double worldYDimension = 1. * m;
0101     G4double worldZDimension = 1. * m;
0102 
0103     fWorld_solid = new G4Box("WorldSolid", worldXDimension, worldYDimension, worldZDimension);
0104 
0105     fWorld_logic = new G4LogicalVolume(fWorld_solid, fAir, "WorldLogical", 0, 0, 0);
0106 
0107     fWorld_phys = new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "World", fWorld_logic, 0, false, 0);
0108 
0109     fWorld_logic->SetVisAttributes(G4VisAttributes::GetInvisible());
0110 
0111 #ifdef G4_DCMTK
0112     ReadPhantomDataNew();
0113     ConstructPhantomContainerNew();
0114 #else
0115     ReadPhantomData();
0116     ConstructPhantomContainer();
0117 #endif
0118 
0119     ConstructPhantom();
0120   }
0121   return fWorld_phys;
0122 }
0123 
0124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
0125 void DicomDetectorConstruction::InitialisationOfMaterials()
0126 {
0127   // Creating elements :
0128   G4double z, a, density;
0129   G4String name, symbol;
0130 
0131   G4Element* elC = new G4Element(name = "Carbon", symbol = "C", z = 6.0, a = 12.011 * g / mole);
0132   G4Element* elH = new G4Element(name = "Hydrogen", symbol = "H", z = 1.0, a = 1.008 * g / mole);
0133   G4Element* elN = new G4Element(name = "Nitrogen", symbol = "N", z = 7.0, a = 14.007 * g / mole);
0134   G4Element* elO = new G4Element(name = "Oxygen", symbol = "O", z = 8.0, a = 16.00 * g / mole);
0135   G4Element* elNa =
0136     new G4Element(name = "Sodium", symbol = "Na", z = 11.0, a = 22.98977 * g / mole);
0137   G4Element* elMg =
0138     new G4Element(name = "Magnesium", symbol = "Mg", z = 12.0, a = 24.3050 * g / mole);
0139   G4Element* elP =
0140     new G4Element(name = "Phosphorus", symbol = "P", z = 15.0, a = 30.973976 * g / mole);
0141   G4Element* elS = new G4Element(name = "Sulfur", symbol = "S", z = 16.0, a = 32.065 * g / mole);
0142   G4Element* elCl =
0143     new G4Element(name = "Chlorine", symbol = "Cl", z = 17.0, a = 35.453 * g / mole);
0144   G4Element* elK =
0145     new G4Element(name = "Potassium", symbol = "K", z = 19.0, a = 30.0983 * g / mole);
0146 
0147   G4Element* elFe = new G4Element(name = "Iron", symbol = "Fe", z = 26, a = 56.845 * g / mole);
0148 
0149   G4Element* elCa = new G4Element(name = "Calcium", symbol = "Ca", z = 20.0, a = 40.078 * g / mole);
0150 
0151   G4Element* elZn = new G4Element(name = "Zinc", symbol = "Zn", z = 30.0, a = 65.382 * g / mole);
0152 
0153   // Creating Materials :
0154   G4int numberofElements;
0155 
0156   // Air
0157   fAir = new G4Material("Air", 1.290 * mg / cm3, numberofElements = 2);
0158   fAir->AddElement(elN, 0.7);
0159   fAir->AddElement(elO, 0.3);
0160 
0161   // Soft tissue (ICRP - NIST)
0162   G4Material* softTissue = new G4Material("SoftTissue", 1.00 * g / cm3, numberofElements = 13);
0163   softTissue->AddElement(elH, 10.4472 * perCent);
0164   softTissue->AddElement(elC, 23.219 * perCent);
0165   softTissue->AddElement(elN, 2.488 * perCent);
0166   softTissue->AddElement(elO, 63.0238 * perCent);
0167   softTissue->AddElement(elNa, 0.113 * perCent);
0168   softTissue->AddElement(elMg, 0.0113 * perCent);
0169   softTissue->AddElement(elP, 0.113 * perCent);
0170   softTissue->AddElement(elS, 0.199 * perCent);
0171   softTissue->AddElement(elCl, 0.134 * perCent);
0172   softTissue->AddElement(elK, 0.199 * perCent);
0173   softTissue->AddElement(elCa, 0.023 * perCent);
0174   softTissue->AddElement(elFe, 0.005 * perCent);
0175   softTissue->AddElement(elZn, 0.003 * perCent);
0176 
0177   //  Lung Inhale
0178   G4Material* lunginhale =
0179     new G4Material("LungInhale", density = 0.217 * g / cm3, numberofElements = 9);
0180   lunginhale->AddElement(elH, 0.103);
0181   lunginhale->AddElement(elC, 0.105);
0182   lunginhale->AddElement(elN, 0.031);
0183   lunginhale->AddElement(elO, 0.749);
0184   lunginhale->AddElement(elNa, 0.002);
0185   lunginhale->AddElement(elP, 0.002);
0186   lunginhale->AddElement(elS, 0.003);
0187   lunginhale->AddElement(elCl, 0.002);
0188   lunginhale->AddElement(elK, 0.003);
0189 
0190   // Lung exhale
0191   G4Material* lungexhale =
0192     new G4Material("LungExhale", density = 0.508 * g / cm3, numberofElements = 9);
0193   lungexhale->AddElement(elH, 0.103);
0194   lungexhale->AddElement(elC, 0.105);
0195   lungexhale->AddElement(elN, 0.031);
0196   lungexhale->AddElement(elO, 0.749);
0197   lungexhale->AddElement(elNa, 0.002);
0198   lungexhale->AddElement(elP, 0.002);
0199   lungexhale->AddElement(elS, 0.003);
0200   lungexhale->AddElement(elCl, 0.002);
0201   lungexhale->AddElement(elK, 0.003);
0202 
0203   // Adipose tissue
0204   G4Material* adiposeTissue =
0205     new G4Material("AdiposeTissue", density = 0.967 * g / cm3, numberofElements = 7);
0206   adiposeTissue->AddElement(elH, 0.114);
0207   adiposeTissue->AddElement(elC, 0.598);
0208   adiposeTissue->AddElement(elN, 0.007);
0209   adiposeTissue->AddElement(elO, 0.278);
0210   adiposeTissue->AddElement(elNa, 0.001);
0211   adiposeTissue->AddElement(elS, 0.001);
0212   adiposeTissue->AddElement(elCl, 0.001);
0213 
0214   // Brain (ICRP - NIST)
0215   G4Material* brainTissue = new G4Material("BrainTissue", 1.03 * g / cm3, numberofElements = 13);
0216   brainTissue->AddElement(elH, 11.0667 * perCent);
0217   brainTissue->AddElement(elC, 12.542 * perCent);
0218   brainTissue->AddElement(elN, 1.328 * perCent);
0219   brainTissue->AddElement(elO, 73.7723 * perCent);
0220   brainTissue->AddElement(elNa, 0.1840 * perCent);
0221   brainTissue->AddElement(elMg, 0.015 * perCent);
0222   brainTissue->AddElement(elP, 0.356 * perCent);
0223   brainTissue->AddElement(elS, 0.177 * perCent);
0224   brainTissue->AddElement(elCl, 0.236 * perCent);
0225   brainTissue->AddElement(elK, 0.31 * perCent);
0226   brainTissue->AddElement(elCa, 0.009 * perCent);
0227   brainTissue->AddElement(elFe, 0.005 * perCent);
0228   brainTissue->AddElement(elZn, 0.001 * perCent);
0229 
0230   // Breast
0231   G4Material* breast = new G4Material("Breast", density = 0.990 * g / cm3, numberofElements = 8);
0232   breast->AddElement(elH, 0.109);
0233   breast->AddElement(elC, 0.506);
0234   breast->AddElement(elN, 0.023);
0235   breast->AddElement(elO, 0.358);
0236   breast->AddElement(elNa, 0.001);
0237   breast->AddElement(elP, 0.001);
0238   breast->AddElement(elS, 0.001);
0239   breast->AddElement(elCl, 0.001);
0240 
0241   // Spinal Disc
0242   G4Material* spinalDisc = new G4Material("SpinalDisc", 1.10 * g / cm3, numberofElements = 8);
0243   spinalDisc->AddElement(elH, 9.60 * perCent);
0244   spinalDisc->AddElement(elC, 9.90 * perCent);
0245   spinalDisc->AddElement(elN, 2.20 * perCent);
0246   spinalDisc->AddElement(elO, 74.40 * perCent);
0247   spinalDisc->AddElement(elNa, 0.50 * perCent);
0248   spinalDisc->AddElement(elP, 2.20 * perCent);
0249   spinalDisc->AddElement(elS, 0.90 * perCent);
0250   spinalDisc->AddElement(elCl, 0.30 * perCent);
0251 
0252   // Water
0253   G4Material* water = new G4Material("Water", density = 1.0 * g / cm3, numberofElements = 2);
0254   water->AddElement(elH, 0.112);
0255   water->AddElement(elO, 0.888);
0256 
0257   // Muscle
0258   G4Material* muscle = new G4Material("Muscle", density = 1.061 * g / cm3, numberofElements = 9);
0259   muscle->AddElement(elH, 0.102);
0260   muscle->AddElement(elC, 0.143);
0261   muscle->AddElement(elN, 0.034);
0262   muscle->AddElement(elO, 0.710);
0263   muscle->AddElement(elNa, 0.001);
0264   muscle->AddElement(elP, 0.002);
0265   muscle->AddElement(elS, 0.003);
0266   muscle->AddElement(elCl, 0.001);
0267   muscle->AddElement(elK, 0.004);
0268 
0269   // Liver
0270   G4Material* liver = new G4Material("Liver", density = 1.071 * g / cm3, numberofElements = 9);
0271   liver->AddElement(elH, 0.102);
0272   liver->AddElement(elC, 0.139);
0273   liver->AddElement(elN, 0.030);
0274   liver->AddElement(elO, 0.716);
0275   liver->AddElement(elNa, 0.002);
0276   liver->AddElement(elP, 0.003);
0277   liver->AddElement(elS, 0.003);
0278   liver->AddElement(elCl, 0.002);
0279   liver->AddElement(elK, 0.003);
0280 
0281   // Tooth Dentin
0282   G4Material* toothDentin = new G4Material("ToothDentin", 2.14 * g / cm3, numberofElements = 10);
0283   toothDentin->AddElement(elH, 2.67 * perCent);
0284   toothDentin->AddElement(elC, 12.77 * perCent);
0285   toothDentin->AddElement(elN, 4.27 * perCent);
0286   toothDentin->AddElement(elO, 40.40 * perCent);
0287   toothDentin->AddElement(elNa, 0.65 * perCent);
0288   toothDentin->AddElement(elMg, 0.59 * perCent);
0289   toothDentin->AddElement(elP, 11.86 * perCent);
0290   toothDentin->AddElement(elCl, 0.04 * perCent);
0291   toothDentin->AddElement(elCa, 26.74 * perCent);
0292   toothDentin->AddElement(elZn, 0.01 * perCent);
0293 
0294   // Trabecular Bone
0295   G4Material* trabecularBone =
0296     new G4Material("TrabecularBone", density = 1.159 * g / cm3, numberofElements = 12);
0297   trabecularBone->AddElement(elH, 0.085);
0298   trabecularBone->AddElement(elC, 0.404);
0299   trabecularBone->AddElement(elN, 0.058);
0300   trabecularBone->AddElement(elO, 0.367);
0301   trabecularBone->AddElement(elNa, 0.001);
0302   trabecularBone->AddElement(elMg, 0.001);
0303   trabecularBone->AddElement(elP, 0.034);
0304   trabecularBone->AddElement(elS, 0.002);
0305   trabecularBone->AddElement(elCl, 0.002);
0306   trabecularBone->AddElement(elK, 0.001);
0307   trabecularBone->AddElement(elCa, 0.044);
0308   trabecularBone->AddElement(elFe, 0.001);
0309 
0310   // Trabecular bone used in the DICOM Head
0311 
0312   G4Material* trabecularBone_head =
0313     new G4Material("TrabecularBone_HEAD", 1.18 * g / cm3, numberofElements = 12);
0314   trabecularBone_head->AddElement(elH, 8.50 * perCent);
0315   trabecularBone_head->AddElement(elC, 40.40 * perCent);
0316   trabecularBone_head->AddElement(elN, 2.80 * perCent);
0317   trabecularBone_head->AddElement(elO, 36.70 * perCent);
0318   trabecularBone_head->AddElement(elNa, 0.10 * perCent);
0319   trabecularBone_head->AddElement(elMg, 0.10 * perCent);
0320   trabecularBone_head->AddElement(elP, 3.40 * perCent);
0321   trabecularBone_head->AddElement(elS, 0.20 * perCent);
0322   trabecularBone_head->AddElement(elCl, 0.20 * perCent);
0323   trabecularBone_head->AddElement(elK, 0.10 * perCent);
0324   trabecularBone_head->AddElement(elCa, 7.40 * perCent);
0325   trabecularBone_head->AddElement(elFe, 0.10 * perCent);
0326 
0327   // Dense Bone
0328   G4Material* denseBone =
0329     new G4Material("DenseBone", density = 1.575 * g / cm3, numberofElements = 11);
0330   denseBone->AddElement(elH, 0.056);
0331   denseBone->AddElement(elC, 0.235);
0332   denseBone->AddElement(elN, 0.050);
0333   denseBone->AddElement(elO, 0.434);
0334   denseBone->AddElement(elNa, 0.001);
0335   denseBone->AddElement(elMg, 0.001);
0336   denseBone->AddElement(elP, 0.072);
0337   denseBone->AddElement(elS, 0.003);
0338   denseBone->AddElement(elCl, 0.001);
0339   denseBone->AddElement(elK, 0.001);
0340   denseBone->AddElement(elCa, 0.146);
0341 
0342   // Cortical Bone (ICRP - NIST)
0343   G4Material* corticalBone = new G4Material("CorticalBone", 1.85 * g / cm3, numberofElements = 9);
0344   corticalBone->AddElement(elH, 4.7234 * perCent);
0345   corticalBone->AddElement(elC, 14.4330 * perCent);
0346   corticalBone->AddElement(elN, 4.199 * perCent);
0347   corticalBone->AddElement(elO, 44.6096 * perCent);
0348   corticalBone->AddElement(elMg, 0.22 * perCent);
0349   corticalBone->AddElement(elP, 10.497 * perCent);
0350   corticalBone->AddElement(elS, 0.315 * perCent);
0351   corticalBone->AddElement(elCa, 20.993 * perCent);
0352   corticalBone->AddElement(elZn, 0.01 * perCent);
0353 
0354   // Tooth enamel
0355   G4Material* toothEnamel = new G4Material("ToothEnamel", 2.89 * g / cm3, numberofElements = 10);
0356   toothEnamel->AddElement(elH, 0.95 * perCent);
0357   toothEnamel->AddElement(elC, 1.11 * perCent);
0358   toothEnamel->AddElement(elN, 0.23 * perCent);
0359   toothEnamel->AddElement(elO, 41.66 * perCent);
0360   toothEnamel->AddElement(elNa, 0.79 * perCent);
0361   toothEnamel->AddElement(elMg, 0.23 * perCent);
0362   toothEnamel->AddElement(elP, 18.71 * perCent);
0363   toothEnamel->AddElement(elCl, 0.34 * perCent);
0364   toothEnamel->AddElement(elCa, 35.97 * perCent);
0365   toothEnamel->AddElement(elZn, 0.02 * perCent);
0366 
0367 #ifdef DICOM_USE_HEAD
0368   //----- Put the materials in a vector HEAD PHANTOM
0369   fOriginalMaterials.push_back(fAir);  // 0.00129 g/cm3
0370   fOriginalMaterials.push_back(softTissue);  // 1.055 g/cm3
0371   fOriginalMaterials.push_back(brainTissue);  // 1.07 g/cm3
0372   fOriginalMaterials.push_back(spinalDisc);  // 1.10 g/cm3
0373   fOriginalMaterials.push_back(trabecularBone_head);  // 1.13 g/cm3
0374   fOriginalMaterials.push_back(toothDentin);  // 1.66 g/cm3
0375   fOriginalMaterials.push_back(corticalBone);  // 1.75 g/cm3
0376   fOriginalMaterials.push_back(toothEnamel);  // 2.04 g/cm3
0377   G4cout << "The materials of the DICOM Head have been used" << G4endl;
0378 #else
0379   fOriginalMaterials.push_back(fAir);  // rho = 0.00129
0380   fOriginalMaterials.push_back(lunginhale);  // rho = 0.217
0381   fOriginalMaterials.push_back(lungexhale);  // rho = 0.508
0382   fOriginalMaterials.push_back(adiposeTissue);  // rho = 0.967
0383   fOriginalMaterials.push_back(breast);  // rho = 0.990
0384   fOriginalMaterials.push_back(water);  // rho = 1.018
0385   fOriginalMaterials.push_back(muscle);  // rho = 1.061
0386   fOriginalMaterials.push_back(liver);  // rho = 1.071
0387   fOriginalMaterials.push_back(trabecularBone);  // rho = 1.159 - HEAD PHANTOM
0388   fOriginalMaterials.push_back(denseBone);  // rho = 1.575
0389   G4cout << "Default materials of the DICOM Extended examples have been used" << G4endl;
0390 #endif
0391 }
0392 
0393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0394 void DicomDetectorConstruction::ReadPhantomDataNew()
0395 {
0396 #ifdef G4_DCMTK
0397   G4String fileName = DicomFileMgr::GetInstance()->GetFileOutName();
0398 
0399   std::ifstream fin(fileName);
0400   std::vector<G4String> wl;
0401   G4int nMaterials;
0402   fin >> nMaterials;
0403   G4String mateName;
0404   G4int nmate;
0405   for (G4int ii = 0; ii < nMaterials; ++ii) {
0406     fin >> nmate;
0407     fin >> mateName;
0408     if (mateName[0] == '"' && mateName[mateName.length() - 1] == '"') {
0409       mateName = mateName.substr(1, mateName.length() - 2);
0410     }
0411     G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate " << ii << " = " << nmate
0412            << " mate " << mateName << G4endl;
0413     if (ii != nmate)
0414       G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong argument",
0415                   FatalErrorInArgument,
0416                   "Material number should be in increasing order:wrong material number");
0417 
0418     G4Material* mate = 0;
0419     const G4MaterialTable* matTab = G4Material::GetMaterialTable();
0420     for (auto matite = matTab->cbegin(); matite != matTab->cend(); ++matite) {
0421       if ((*matite)->GetName() == mateName) {
0422         mate = *matite;
0423       }
0424     }
0425     if (mate == 0) {
0426       mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName);
0427     }
0428     if (!mate)
0429       G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong argument",
0430                   FatalErrorInArgument, ("Material not found" + mateName).c_str());
0431     fPhantomMaterialsOriginal[nmate] = mate;
0432   }
0433 
0434   fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
0435   G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNoVoxels X/Y/Z " << fNoVoxelsX << " "
0436          << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
0437   fin >> fMinX >> fMaxX;
0438   fin >> fMinY >> fMaxY;
0439   fin >> fMinZ >> fMaxZ;
0440   fVoxelHalfDimX = (fMaxX - fMinX) / fNoVoxelsX / 2.;
0441   fVoxelHalfDimY = (fMaxY - fMinY) / fNoVoxelsY / 2.;
0442   fVoxelHalfDimZ = (fMaxZ - fMinZ) / fNoVoxelsZ / 2.;
0443 #  ifdef G4VERBOSE
0444   G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl << " Extension in Y " << fMinY
0445          << " " << fMaxY << G4endl << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl;
0446 #  endif
0447 
0448   fMateIDs = new size_t[fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ];
0449   for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {
0450     for (G4int iy = 0; iy < fNoVoxelsY; ++iy) {
0451       for (G4int ix = 0; ix < fNoVoxelsX; ++ix) {
0452         G4int mateID;
0453         fin >> mateID;
0454         G4int nnew = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX * fNoVoxelsY;
0455         if (mateID < 0 || mateID >= nMaterials) {
0456           G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", "Wrong index in phantom file",
0457                       FatalException,
0458                       G4String("It should be between 0 and "
0459                                + G4UIcommand::ConvertToString(nMaterials - 1) + ", while it is "
0460                                + G4UIcommand::ConvertToString(mateID))
0461                         .c_str());
0462         }
0463         fMateIDs[nnew] = mateID;
0464       }
0465     }
0466   }
0467 
0468   ReadVoxelDensities(fin);
0469 
0470   fin.close();
0471 #endif
0472 }
0473 
0474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0475 void DicomDetectorConstruction::ReadVoxelDensities(std::ifstream& fin)
0476 {
0477   G4String stemp;
0478   std::map<G4int, std::pair<G4double, G4double>> densiMinMax;
0479   std::map<G4int, std::pair<G4double, G4double>>::iterator mpite;
0480   for (G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii) {
0481     densiMinMax[ii] = std::pair<G4double, G4double>(DBL_MAX, -DBL_MAX);
0482   }
0483 
0484   char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
0485   G4double densityDiff = -1.;
0486   if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
0487 
0488   std::map<G4int, G4double> densityDiffs;
0489   for (G4int ii = 0; ii < G4int(fPhantomMaterialsOriginal.size()); ++ii) {
0490     densityDiffs[ii] = densityDiff;  // currently all materials with same step
0491   }
0492   //  densityDiffs[0] = 0.0001; //air
0493 
0494   //--- Calculate the average material density for each material/density bin
0495   std::map<std::pair<G4Material*, G4int>, matInfo*> newMateDens;
0496 
0497   //---- Read the material densities
0498   G4double dens;
0499   for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {
0500     for (G4int iy = 0; iy < fNoVoxelsY; ++iy) {
0501       for (G4int ix = 0; ix < fNoVoxelsX; ++ix) {
0502         fin >> dens;
0503         G4int copyNo = ix + (iy)*fNoVoxelsX + (iz)*fNoVoxelsX * fNoVoxelsY;
0504 
0505         if (densityDiff != -1.) continue;
0506 
0507         //--- store the minimum and maximum density for each material
0508         mpite = densiMinMax.find(G4int(fMateIDs[copyNo]));
0509         if (dens < (*mpite).second.first) (*mpite).second.first = dens;
0510         if (dens > (*mpite).second.second) (*mpite).second.second = dens;
0511         //--- Get material from original list of material in file
0512         G4int mateID = G4int(fMateIDs[copyNo]);
0513         std::map<G4int, G4Material*>::const_iterator imite = fPhantomMaterialsOriginal.find(mateID);
0514 
0515         //--- Check if density is equal to the original material density
0516         if (std::fabs(dens - (*imite).second->GetDensity() / CLHEP::g * CLHEP::cm3) < 1.e-9)
0517           continue;
0518 
0519         //--- Build material name with fPhantomMaterialsOriginal name+density
0520         G4int densityBin = (G4int(dens / densityDiffs[mateID]));
0521 
0522         G4String mateName = (*imite).second->GetName() + G4UIcommand::ConvertToString(densityBin);
0523         //--- Look if it is the first voxel with this material/densityBin
0524         std::pair<G4Material*, G4int> matdens((*imite).second, densityBin);
0525 
0526         auto mppite = newMateDens.find(matdens);
0527         if (mppite != newMateDens.cend()) {
0528           matInfo* mi = (*mppite).second;
0529           mi->fSumdens += dens;
0530           mi->fNvoxels++;
0531           fMateIDs[copyNo] = fPhantomMaterialsOriginal.size() - 1 + mi->fId;
0532         }
0533         else {
0534           matInfo* mi = new matInfo;
0535           mi->fSumdens = dens;
0536           mi->fNvoxels = 1;
0537           mi->fId = G4int(newMateDens.size() + 1);
0538           newMateDens[matdens] = mi;
0539           fMateIDs[copyNo] = fPhantomMaterialsOriginal.size() - 1 + mi->fId;
0540         }
0541       }
0542     }
0543   }
0544 
0545   if (densityDiff != -1.) {
0546     for (mpite = densiMinMax.begin(); mpite != densiMinMax.end(); ++mpite) {
0547 #ifdef G4VERBOSE
0548       G4cout << "DicomDetectorConstruction::ReadVoxelDensities"
0549              << " ORIG MATERIALS DENSITY " << (*mpite).first << " MIN " << (*mpite).second.first
0550              << " MAX " << (*mpite).second.second << G4endl;
0551 #endif
0552     }
0553   }
0554 
0555   //----- Build the list of phantom materials that go to Parameterisation
0556   //--- Add original materials
0557   for (auto mimite = fPhantomMaterialsOriginal.cbegin(); mimite != fPhantomMaterialsOriginal.cend();
0558        ++mimite)
0559   {
0560     fMaterials.push_back((*mimite).second);
0561   }
0562   //
0563   //---- Build and add new materials
0564   for (auto mppite = newMateDens.cbegin(); mppite != newMateDens.cend(); ++mppite) {
0565     G4double averdens = (*mppite).second->fSumdens / (*mppite).second->fNvoxels;
0566     G4double saverdens = G4int(1000.001 * averdens) / 1000.;
0567 #ifndef G4VERBOSE
0568     G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS " << averdens << " -> "
0569            << saverdens << " -> " << G4int(1000 * averdens) << " " << G4int(1000 * averdens) / 1000
0570            << " " << G4int(1000 * averdens) / 1000. << G4endl;
0571 #endif
0572 
0573     G4String mateName =
0574       ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
0575     fMaterials.push_back(
0576       BuildMaterialWithChangingDensity((*mppite).first.first, G4float(averdens), mateName));
0577   }
0578 }
0579 
0580 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.......
0581 void DicomDetectorConstruction::ReadPhantomData()
0582 {
0583   G4String dataFile = DicomHandler::GetDicomDataFile();
0584   std::ifstream finDF(dataFile.c_str());
0585   G4String fname;
0586 
0587   if (finDF.good() != 1) {
0588     G4String descript = "Problem reading data file: " + dataFile;
0589     G4Exception(" DicomDetectorConstruction::ReadPhantomData", " ", FatalException, descript);
0590   }
0591 
0592   G4int compression;
0593   finDF >> compression;  // not used here
0594   finDF >> fNoFiles;
0595 
0596   for (G4int i = 0; i < fNoFiles; ++i) {
0597     finDF >> fname;
0598 
0599     //--- Read one data file
0600     fname += ".g4dcm";
0601 
0602     ReadPhantomDataFile(fname);
0603   }
0604 
0605   //----- Merge data headers
0606   MergeZSliceHeaders();
0607   finDF.close();
0608 }
0609 
0610 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........
0611 void DicomDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
0612 {
0613   G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname
0614          << G4endl;  // GDEB
0615 
0616 #ifdef G4VERBOSE
0617   G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
0618 #endif
0619 
0620   std::ifstream fin(fname.c_str(), std::ios_base::in);
0621   if (!fin.is_open()) {
0622     G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", "", FatalErrorInArgument,
0623                 G4String("File not found " + fname).c_str());
0624   }
0625   //----- Define density differences (maximum density difference to create
0626   // a new material)
0627   char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
0628   G4double densityDiff = -1.;
0629   if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
0630   if (densityDiff != -1.) {
0631     for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
0632       fDensityDiffs[ii] = densityDiff;  // currently all materials with
0633       // same difference
0634     }
0635   }
0636   else {
0637     if (fMaterials.size() == 0) {  // do it only for first slice
0638       for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
0639         fMaterials.push_back(fOriginalMaterials[ii]);
0640       }
0641     }
0642   }
0643 
0644   //----- Read data header
0645   DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader(fin);
0646   fZSliceHeaders.push_back(sliceHeader);
0647 
0648   //----- Read material indices
0649   G4int nVoxels = sliceHeader->GetNoVoxels();
0650 
0651   //--- If first slice, initiliaze fMateIDs
0652   if (fZSliceHeaders.size() == 1) {
0653     // fMateIDs = new unsigned int[fNoFiles*nVoxels];
0654     fMateIDs = new size_t[fNoFiles * nVoxels];
0655   }
0656 
0657   unsigned int mateID;
0658   // number of voxels from previously read slices
0659   G4int voxelCopyNo = G4int((fZSliceHeaders.size() - 1) * nVoxels);
0660   for (G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++) {
0661     fin >> mateID;
0662     fMateIDs[voxelCopyNo] = mateID;
0663   }
0664 
0665   //----- Read material densities and build new materials if two voxels have
0666   //  same material but its density is in a different density interval
0667   // (size of density intervals defined by densityDiff)
0668   G4double density;
0669   // number of voxels from previously read slices
0670   voxelCopyNo = G4int((fZSliceHeaders.size() - 1) * nVoxels);
0671   for (G4int ii = 0; ii < nVoxels; ++ii, voxelCopyNo++) {
0672     fin >> density;
0673 
0674     //-- Get material from list of original materials
0675     mateID = unsigned(fMateIDs[voxelCopyNo]);
0676     G4Material* mateOrig = fOriginalMaterials[mateID];
0677 
0678     //-- Get density bin: middle point of the bin in which the current
0679     // density is included
0680     G4String newMateName = mateOrig->GetName();
0681     G4float densityBin = 0.;
0682     if (densityDiff != -1.) {
0683       densityBin = G4float(fDensityDiffs[mateID]) * (G4int(density / fDensityDiffs[mateID]) + 0.5);
0684       //-- Build the new material name
0685       newMateName += G4UIcommand::ConvertToString(densityBin);
0686     }
0687 
0688     //-- Look if a material with this name is already created
0689     //  (because a previous voxel was already in this density bin)
0690     unsigned int im;
0691     for (im = 0; im < fMaterials.size(); ++im) {
0692       if (fMaterials[im]->GetName() == newMateName) {
0693         break;
0694       }
0695     }
0696     //-- If material is already created use index of this material
0697     if (im != fMaterials.size()) {
0698       fMateIDs[voxelCopyNo] = im;
0699       //-- else, create the material
0700     }
0701     else {
0702       if (densityDiff != -1.) {
0703         fMaterials.push_back(BuildMaterialWithChangingDensity(mateOrig, densityBin, newMateName));
0704         fMateIDs[voxelCopyNo] = fMaterials.size() - 1;
0705       }
0706       else {
0707         G4cerr << " im " << im << " < " << fMaterials.size() << " name " << newMateName << G4endl;
0708         G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", "", FatalErrorInArgument,
0709                     "Wrong index in material");  // it should never reach here
0710       }
0711     }
0712   }
0713 }
0714 
0715 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0716 void DicomDetectorConstruction::MergeZSliceHeaders()
0717 {
0718   //----- Images must have the same dimension ...
0719   fZSliceHeaderMerged = new DicomPhantomZSliceHeader(*fZSliceHeaders[0]);
0720   for (unsigned int ii = 1; ii < fZSliceHeaders.size(); ++ii) {
0721     *fZSliceHeaderMerged += *fZSliceHeaders[ii];
0722   }
0723 }
0724 
0725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0726 G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity(const G4Material* origMate,
0727                                                                         G4float density,
0728                                                                         G4String newMateName)
0729 {
0730   //----- Copy original material, but with new density
0731   G4int nelem = G4int(origMate->GetNumberOfElements());
0732   G4Material* mate =
0733     new G4Material(newMateName, density * g / cm3, nelem, kStateUndefined, STP_Temperature);
0734 
0735   for (G4int ii = 0; ii < nelem; ++ii) {
0736     G4double frac = origMate->GetFractionVector()[ii];
0737     G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii));
0738     mate->AddElement(elem, frac);
0739   }
0740 
0741   return mate;
0742 }
0743 
0744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0745 void DicomDetectorConstruction::ConstructPhantomContainer()
0746 {
0747   //---- Extract number of voxels and voxel dimensions
0748   fNoVoxelsX = fZSliceHeaderMerged->GetNoVoxelsX();
0749   fNoVoxelsY = fZSliceHeaderMerged->GetNoVoxelsY();
0750   fNoVoxelsZ = fZSliceHeaderMerged->GetNoVoxelsZ();
0751 
0752   fVoxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX();
0753   fVoxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY();
0754   fVoxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ();
0755 #ifdef G4VERBOSE
0756   G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX << G4endl;
0757   G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY << G4endl;
0758   G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ << G4endl;
0759   G4cout << " totalPixels " << fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ << G4endl;
0760 #endif
0761 
0762   //----- Define the volume that contains all the voxels
0763   fContainer_solid = new G4Box("phantomContainer", fNoVoxelsX * fVoxelHalfDimX,
0764                                fNoVoxelsY * fVoxelHalfDimY, fNoVoxelsZ * fVoxelHalfDimZ);
0765   fContainer_logic =
0766     new G4LogicalVolume(fContainer_solid,
0767                         // the material is not important, it will be fully filled by the voxels
0768                         fMaterials[0], "phantomContainer", 0, 0, 0);
0769   //--- Place it on the world
0770   G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() + fZSliceHeaderMerged->GetMinX()) / 2.;
0771   G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() + fZSliceHeaderMerged->GetMinY()) / 2.;
0772   G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() + fZSliceHeaderMerged->GetMinZ()) / 2.;
0773   G4ThreeVector posCentreVoxels(fOffsetX, fOffsetY, fOffsetZ);
0774 #ifdef G4VERBOSE
0775   G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
0776 #endif
0777   fContainer_phys = new G4PVPlacement(0,  // rotation
0778                                       posCentreVoxels,
0779                                       fContainer_logic,  // The logic volume
0780                                       "phantomContainer",  // Name
0781                                       fWorld_logic,  // Mother
0782                                       false,  // No op. bool.
0783                                       1);  // Copy number
0784 }
0785 
0786 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0787 void DicomDetectorConstruction::ConstructPhantomContainerNew()
0788 {
0789 #ifdef G4_DCMTK
0790   //---- Extract number of voxels and voxel dimensions
0791 #  ifdef G4VERBOSE
0792   G4cout << " fNoVoxelsX " << fNoVoxelsX << " fVoxelHalfDimX " << fVoxelHalfDimX << G4endl;
0793   G4cout << " fNoVoxelsY " << fNoVoxelsY << " fVoxelHalfDimY " << fVoxelHalfDimY << G4endl;
0794   G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ << G4endl;
0795   G4cout << " totalPixels " << fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ << G4endl;
0796 #  endif
0797 
0798   //----- Define the volume that contains all the voxels
0799   fContainer_solid = new G4Box("phantomContainer", fNoVoxelsX * fVoxelHalfDimX,
0800                                fNoVoxelsY * fVoxelHalfDimY, fNoVoxelsZ * fVoxelHalfDimZ);
0801   fContainer_logic =
0802     new G4LogicalVolume(fContainer_solid,
0803                         // the material is not important, it will be fully filled by the voxels
0804                         fMaterials[0], "phantomContainer", 0, 0, 0);
0805 
0806   G4ThreeVector posCentreVoxels((fMinX + fMaxX) / 2., (fMinY + fMaxY) / 2., (fMinZ + fMaxZ) / 2.);
0807 #  ifdef G4VERBOSE
0808   G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl;
0809 #  endif
0810   fContainer_phys = new G4PVPlacement(0,  // rotation
0811                                       posCentreVoxels,
0812                                       fContainer_logic,  // The logic volume
0813                                       "phantomContainer",  // Name
0814                                       fWorld_logic,  // Mother
0815                                       false,  // No op. bool.
0816                                       1);  // Copy number
0817 #endif
0818 }
0819 
0820 #include "G4MultiFunctionalDetector.hh"
0821 #include "G4PSDoseDeposit.hh"
0822 #include "G4PSDoseDeposit3D.hh"
0823 #include "G4SDManager.hh"
0824 
0825 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.
0826 void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic)
0827 {
0828 #ifdef G4VERBOSE
0829   G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl;
0830 #endif
0831 
0832   fScorers.insert(voxel_logic);
0833 }
0834 
0835 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
0836 
0837 void DicomDetectorConstruction::ConstructSDandField()
0838 {
0839 #ifdef G4VERBOSE
0840   G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl;
0841 #endif
0842 
0843   // G4SDManager* SDman = G4SDManager::GetSDMpointer();
0844 
0845   // SDman->SetVerboseLevel(1);
0846 
0847   //
0848   // Sensitive Detector Name
0849   G4String concreteSDname = "phantomSD";
0850   std::vector<G4String> scorer_names;
0851   scorer_names.push_back(concreteSDname);
0852   //------------------------
0853   // MultiFunctionalDetector
0854   //------------------------
0855   //
0856   // Define MultiFunctionalDetector with name.
0857   // declare MFDet as a MultiFunctionalDetector scorer
0858   G4MultiFunctionalDetector* MFDet = new G4MultiFunctionalDetector(concreteSDname);
0859   G4SDManager::GetSDMpointer()->AddNewDetector(MFDet);
0860   char* nest = std::getenv("DICOM_NESTED_PARAM");
0861   if (nest && G4String(nest) == "1") {
0862     G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit3D(
0863       "DoseDeposit", fNoVoxelsZ, fNoVoxelsY, fNoVoxelsX, 0, 2, 1);  // nested param replica indexing
0864     // - the last 3 arguments correspond to the replica depth for Z, Y and X respectively
0865     MFDet->RegisterPrimitive(dosedep);
0866   }
0867   else {
0868     G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit");  // regular geometry
0869     MFDet->RegisterPrimitive(dosedep);
0870   }
0871 
0872   for (auto ite = fScorers.cbegin(); ite != fScorers.cend(); ++ite) {
0873     SetSensitiveDetector(*ite, MFDet);
0874   }
0875 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!