Geant4 Cross Reference |
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 // $Id$ 26 // 27 // 27 /// \file medical/DICOM/src/DicomDetectorCons 28 /// \file medical/DICOM/src/DicomDetectorConstruction.cc 28 /// \brief Implementation of the DicomDetector 29 /// \brief Implementation of the DicomDetectorConstruction class 29 // 30 // 30 31 31 #include "DicomDetectorConstruction.hh" << 32 #include "globals.hh" 32 << 33 #include "CLHEP/Units/SystemOfUnits.h" << 34 #include "DicomHandler.hh" << 35 #include "DicomPhantomZSliceHeader.hh" << 36 33 37 #include "G4Box.hh" 34 #include "G4Box.hh" 38 #include "G4Element.hh" << 39 #include "G4LogicalVolume.hh" 35 #include "G4LogicalVolume.hh" 40 #include "G4Material.hh" << 36 #include "G4VPhysicalVolume.hh" 41 #include "G4NistManager.hh" << 42 #include "G4PVPlacement.hh" 37 #include "G4PVPlacement.hh" 43 #include "G4PhysicalConstants.hh" << 38 #include "G4Material.hh" >> 39 #include "G4Element.hh" 44 #include "G4UIcommand.hh" 40 #include "G4UIcommand.hh" 45 #include "G4VPhysicalVolume.hh" << 41 #include "G4PhysicalConstants.hh" 46 #include "globals.hh" << 42 #include "G4NistManager.hh" >> 43 // #include "G4SystemOfUnits.hh" >> 44 #include "CLHEP/Units/SystemOfUnits.h" 47 45 >> 46 #include "DicomDetectorConstruction.hh" >> 47 #include "DicomPhantomZSliceHeader.hh" >> 48 >> 49 #include "DicomRunAction.hh" >> 50 #include "DicomRun.hh" 48 #ifdef G4_DCMTK 51 #ifdef G4_DCMTK 49 # include "DicomFileMgr.hh" << 52 #include "DicomFileMgr.hh" 50 #endif 53 #endif 51 #include "G4VisAttributes.hh" 54 #include "G4VisAttributes.hh" 52 55 >> 56 using CLHEP::m; 53 using CLHEP::cm3; 57 using CLHEP::cm3; >> 58 using CLHEP::mole; 54 using CLHEP::g; 59 using CLHEP::g; 55 using CLHEP::m; << 56 using CLHEP::mg; 60 using CLHEP::mg; 57 using CLHEP::mole; << 58 using CLHEP::perCent; 61 using CLHEP::perCent; 59 62 60 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.. 61 DicomDetectorConstruction::DicomDetectorConstr 64 DicomDetectorConstruction::DicomDetectorConstruction() 62 : G4VUserDetectorConstruction(), << 65 : G4VUserDetectorConstruction(), 63 fAir(0), << 66 fAir(0), 64 67 65 fWorld_solid(0), << 68 fWorld_solid(0), 66 fWorld_logic(0), << 69 fWorld_logic(0), 67 fWorld_phys(0), << 70 fWorld_phys(0), 68 << 71 69 fContainer_solid(0), << 72 fContainer_solid(0), 70 fContainer_logic(0), << 73 fContainer_logic(0), 71 fContainer_phys(0), << 74 fContainer_phys(0), 72 << 75 73 fNoFiles(0), << 76 fNoFiles(0), 74 fMateIDs(0), << 77 fMateIDs(0), 75 << 78 76 fZSliceHeaderMerged(0), << 79 fZSliceHeaderMerged(0), 77 << 80 78 fNoVoxelsX(0), << 81 fNVoxelX(0), 79 fNoVoxelsY(0), << 82 fNVoxelY(0), 80 fNoVoxelsZ(0), << 83 fNVoxelZ(0), 81 fVoxelHalfDimX(0), << 84 fVoxelHalfDimX(0), 82 fVoxelHalfDimY(0), << 85 fVoxelHalfDimY(0), 83 fVoxelHalfDimZ(0), << 86 fVoxelHalfDimZ(0), 84 87 85 fConstructed(false) << 88 fConstructed(false) 86 {} << 89 { >> 90 >> 91 } 87 92 88 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo... 89 DicomDetectorConstruction::~DicomDetectorConst << 94 DicomDetectorConstruction::~DicomDetectorConstruction() >> 95 { >> 96 } 90 97 91 //....oooOO0OOooo........oooOO0OOooo........oo 98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 92 G4VPhysicalVolume* DicomDetectorConstruction:: 99 G4VPhysicalVolume* DicomDetectorConstruction::Construct() 93 { 100 { 94 if (!fConstructed || fWorld_phys == 0) { << 101 if(!fConstructed || fWorld_phys == 0) { 95 fConstructed = true; 102 fConstructed = true; 96 InitialisationOfMaterials(); 103 InitialisationOfMaterials(); 97 104 98 //----- Build world 105 //----- Build world 99 G4double worldXDimension = 1. * m; << 106 G4double worldXDimension = 1.*m; 100 G4double worldYDimension = 1. * m; << 107 G4double worldYDimension = 1.*m; 101 G4double worldZDimension = 1. * m; << 108 G4double worldZDimension = 1.*m; 102 << 109 103 fWorld_solid = new G4Box("WorldSolid", wor << 110 fWorld_solid = new G4Box( "WorldSolid", 104 << 111 worldXDimension, 105 fWorld_logic = new G4LogicalVolume(fWorld_ << 112 worldYDimension, 106 << 113 worldZDimension ); 107 fWorld_phys = new G4PVPlacement(0, G4Three << 114 108 << 115 fWorld_logic = new G4LogicalVolume( fWorld_solid, 109 fWorld_logic->SetVisAttributes(G4VisAttrib << 116 fAir, >> 117 "WorldLogical", >> 118 0, 0, 0 ); >> 119 >> 120 fWorld_phys = new G4PVPlacement( 0, >> 121 G4ThreeVector(0,0,0), >> 122 "World", >> 123 fWorld_logic, >> 124 0, >> 125 false, >> 126 0 ); 110 127 111 #ifdef G4_DCMTK 128 #ifdef G4_DCMTK 112 ReadPhantomDataNew(); 129 ReadPhantomDataNew(); 113 ConstructPhantomContainerNew(); 130 ConstructPhantomContainerNew(); 114 #else 131 #else 115 ReadPhantomData(); 132 ReadPhantomData(); 116 ConstructPhantomContainer(); 133 ConstructPhantomContainer(); 117 #endif 134 #endif 118 << 135 119 ConstructPhantom(); 136 ConstructPhantom(); 120 } 137 } 121 return fWorld_phys; << 138 return fWorld_phys; 122 } 139 } 123 140 124 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........ 125 void DicomDetectorConstruction::Initialisation 142 void DicomDetectorConstruction::InitialisationOfMaterials() 126 { 143 { 127 // Creating elements : << 144 // Creating elements : 128 G4double z, a, density; << 145 G4double z, a, density; 129 G4String name, symbol; << 146 G4String name, symbol; 130 << 147 131 G4Element* elC = new G4Element(name = "Carbo << 148 G4Element* elC = new G4Element( name = "Carbon", 132 G4Element* elH = new G4Element(name = "Hydro << 149 symbol = "C", 133 G4Element* elN = new G4Element(name = "Nitro << 150 z = 6.0, a = 12.011 * g/mole ); 134 G4Element* elO = new G4Element(name = "Oxyge << 151 G4Element* elH = new G4Element( name = "Hydrogen", 135 G4Element* elNa = << 152 symbol = "H", 136 new G4Element(name = "Sodium", symbol = "N << 153 z = 1.0, a = 1.008 * g/mole ); 137 G4Element* elMg = << 154 G4Element* elN = new G4Element( name = "Nitrogen", 138 new G4Element(name = "Magnesium", symbol = << 155 symbol = "N", 139 G4Element* elP = << 156 z = 7.0, a = 14.007 * g/mole ); 140 new G4Element(name = "Phosphorus", symbol << 157 G4Element* elO = new G4Element( name = "Oxygen", 141 G4Element* elS = new G4Element(name = "Sulfu << 158 symbol = "O", 142 G4Element* elCl = << 159 z = 8.0, a = 16.00 * g/mole ); 143 new G4Element(name = "Chlorine", symbol = << 160 G4Element* elNa = new G4Element( name = "Sodium", 144 G4Element* elK = << 161 symbol = "Na", 145 new G4Element(name = "Potassium", symbol = << 162 z= 11.0, a = 22.98977* g/mole ); 146 << 163 G4Element* elMg = new G4Element( name = "Magnesium", 147 G4Element* elFe = new G4Element(name = "Iron << 164 symbol = "Mg", 148 << 165 z = 12.0, a = 24.3050* g/mole ); 149 G4Element* elCa = new G4Element(name = "Calc << 166 G4Element* elP = new G4Element( name = "Phosphorus", 150 << 167 symbol = "P", 151 G4Element* elZn = new G4Element(name = "Zinc << 168 z = 15.0, a = 30.973976* g/mole ); 152 << 169 G4Element* elS = new G4Element( name = "Sulfur", 153 // Creating Materials : << 170 symbol = "S", 154 G4int numberofElements; << 171 z = 16.0,a = 32.065* g/mole ); 155 << 172 G4Element* elCl = new G4Element( name = "Chlorine", 156 // Air << 173 symbol = "P", 157 fAir = new G4Material("Air", 1.290 * mg / cm << 174 z = 17.0, a = 35.453* g/mole ); 158 fAir->AddElement(elN, 0.7); << 175 G4Element* elK = new G4Element( name = "Potassium", 159 fAir->AddElement(elO, 0.3); << 176 symbol = "P", 160 << 177 z = 19.0, a = 30.0983* g/mole ); 161 // Soft tissue (ICRP - NIST) << 178 162 G4Material* softTissue = new G4Material("Sof << 179 G4Element* elFe = new G4Element( name = "Iron", 163 softTissue->AddElement(elH, 10.4472 * perCen << 180 symbol = "Fe", 164 softTissue->AddElement(elC, 23.219 * perCent << 181 z = 26, a = 56.845* g/mole ); 165 softTissue->AddElement(elN, 2.488 * perCent) << 182 166 softTissue->AddElement(elO, 63.0238 * perCen << 183 G4Element* elCa = new G4Element( name="Calcium", 167 softTissue->AddElement(elNa, 0.113 * perCent << 184 symbol = "Ca", 168 softTissue->AddElement(elMg, 0.0113 * perCen << 185 z = 20.0, a = 40.078* g/mole ); 169 softTissue->AddElement(elP, 0.113 * perCent) << 186 170 softTissue->AddElement(elS, 0.199 * perCent) << 187 G4Element* elZn = new G4Element( name = "Zinc", 171 softTissue->AddElement(elCl, 0.134 * perCent << 188 symbol = "Zn", 172 softTissue->AddElement(elK, 0.199 * perCent) << 189 z = 30.0,a = 65.382* g/mole ); 173 softTissue->AddElement(elCa, 0.023 * perCent << 190 174 softTissue->AddElement(elFe, 0.005 * perCent << 191 // Creating Materials : 175 softTissue->AddElement(elZn, 0.003 * perCent << 192 G4int numberofElements; 176 << 193 177 // Lung Inhale << 194 // Air 178 G4Material* lunginhale = << 195 fAir = new G4Material( "Air", 179 new G4Material("LungInhale", density = 0.2 << 196 1.290*mg/cm3, 180 lunginhale->AddElement(elH, 0.103); << 197 numberofElements = 2 ); 181 lunginhale->AddElement(elC, 0.105); << 198 fAir->AddElement(elN, 0.7); 182 lunginhale->AddElement(elN, 0.031); << 199 fAir->AddElement(elO, 0.3); 183 lunginhale->AddElement(elO, 0.749); << 200 184 lunginhale->AddElement(elNa, 0.002); << 201 185 lunginhale->AddElement(elP, 0.002); << 202 // Soft tissue (ICRP - NIST) 186 lunginhale->AddElement(elS, 0.003); << 203 G4Material* softTissue = new G4Material ("SoftTissue", 1.00*g/cm3, 187 lunginhale->AddElement(elCl, 0.002); << 204 numberofElements = 13); 188 lunginhale->AddElement(elK, 0.003); << 205 softTissue->AddElement(elH, 10.4472*perCent); 189 << 206 softTissue->AddElement(elC, 23.219*perCent); 190 // Lung exhale << 207 softTissue->AddElement(elN, 2.488*perCent); 191 G4Material* lungexhale = << 208 softTissue->AddElement(elO, 63.0238*perCent); 192 new G4Material("LungExhale", density = 0.5 << 209 softTissue->AddElement(elNa, 0.113*perCent); 193 lungexhale->AddElement(elH, 0.103); << 210 softTissue->AddElement(elMg, 0.0113*perCent); 194 lungexhale->AddElement(elC, 0.105); << 211 softTissue->AddElement(elP, 0.113*perCent); 195 lungexhale->AddElement(elN, 0.031); << 212 softTissue->AddElement(elS, 0.199*perCent); 196 lungexhale->AddElement(elO, 0.749); << 213 softTissue->AddElement(elCl, 0.134*perCent); 197 lungexhale->AddElement(elNa, 0.002); << 214 softTissue->AddElement(elK, 0.199*perCent); 198 lungexhale->AddElement(elP, 0.002); << 215 softTissue->AddElement(elCa, 0.023*perCent); 199 lungexhale->AddElement(elS, 0.003); << 216 softTissue->AddElement(elFe, 0.005*perCent); 200 lungexhale->AddElement(elCl, 0.002); << 217 softTissue->AddElement(elZn, 0.003*perCent); 201 lungexhale->AddElement(elK, 0.003); << 218 202 << 219 // Lung Inhale 203 // Adipose tissue << 220 G4Material* lunginhale = new G4Material( "LungInhale", 204 G4Material* adiposeTissue = << 221 density = 0.217*g/cm3, 205 new G4Material("AdiposeTissue", density = << 222 numberofElements = 9); 206 adiposeTissue->AddElement(elH, 0.114); << 223 lunginhale->AddElement(elH,0.103); 207 adiposeTissue->AddElement(elC, 0.598); << 224 lunginhale->AddElement(elC,0.105); 208 adiposeTissue->AddElement(elN, 0.007); << 225 lunginhale->AddElement(elN,0.031); 209 adiposeTissue->AddElement(elO, 0.278); << 226 lunginhale->AddElement(elO,0.749); 210 adiposeTissue->AddElement(elNa, 0.001); << 227 lunginhale->AddElement(elNa,0.002); 211 adiposeTissue->AddElement(elS, 0.001); << 228 lunginhale->AddElement(elP,0.002); 212 adiposeTissue->AddElement(elCl, 0.001); << 229 lunginhale->AddElement(elS,0.003); 213 << 230 lunginhale->AddElement(elCl,0.002); 214 // Brain (ICRP - NIST) << 231 lunginhale->AddElement(elK,0.003); 215 G4Material* brainTissue = new G4Material("Br << 232 216 brainTissue->AddElement(elH, 11.0667 * perCe << 233 // Lung exhale 217 brainTissue->AddElement(elC, 12.542 * perCen << 234 G4Material* lungexhale = new G4Material( "LungExhale", 218 brainTissue->AddElement(elN, 1.328 * perCent << 235 density = 0.508*g/cm3, 219 brainTissue->AddElement(elO, 73.7723 * perCe << 236 numberofElements = 9 ); 220 brainTissue->AddElement(elNa, 0.1840 * perCe << 237 lungexhale->AddElement(elH,0.103); 221 brainTissue->AddElement(elMg, 0.015 * perCen << 238 lungexhale->AddElement(elC,0.105); 222 brainTissue->AddElement(elP, 0.356 * perCent << 239 lungexhale->AddElement(elN,0.031); 223 brainTissue->AddElement(elS, 0.177 * perCent << 240 lungexhale->AddElement(elO,0.749); 224 brainTissue->AddElement(elCl, 0.236 * perCen << 241 lungexhale->AddElement(elNa,0.002); 225 brainTissue->AddElement(elK, 0.31 * perCent) << 242 lungexhale->AddElement(elP,0.002); 226 brainTissue->AddElement(elCa, 0.009 * perCen << 243 lungexhale->AddElement(elS,0.003); 227 brainTissue->AddElement(elFe, 0.005 * perCen << 244 lungexhale->AddElement(elCl,0.002); 228 brainTissue->AddElement(elZn, 0.001 * perCen << 245 lungexhale->AddElement(elK,0.003); 229 << 246 230 // Breast << 247 // Adipose tissue 231 G4Material* breast = new G4Material("Breast" << 248 G4Material* adiposeTissue = new G4Material( "AdiposeTissue", 232 breast->AddElement(elH, 0.109); << 249 density = 0.967*g/cm3, 233 breast->AddElement(elC, 0.506); << 250 numberofElements = 7); 234 breast->AddElement(elN, 0.023); << 251 adiposeTissue->AddElement(elH,0.114); 235 breast->AddElement(elO, 0.358); << 252 adiposeTissue->AddElement(elC,0.598); 236 breast->AddElement(elNa, 0.001); << 253 adiposeTissue->AddElement(elN,0.007); 237 breast->AddElement(elP, 0.001); << 254 adiposeTissue->AddElement(elO,0.278); 238 breast->AddElement(elS, 0.001); << 255 adiposeTissue->AddElement(elNa,0.001); 239 breast->AddElement(elCl, 0.001); << 256 adiposeTissue->AddElement(elS,0.001); 240 << 257 adiposeTissue->AddElement(elCl,0.001); 241 // Spinal Disc << 258 242 G4Material* spinalDisc = new G4Material("Spi << 259 // Brain (ICRP - NIST) 243 spinalDisc->AddElement(elH, 9.60 * perCent); << 260 G4Material* brainTissue = new G4Material ("BrainTissue", 1.03 * g/cm3, 244 spinalDisc->AddElement(elC, 9.90 * perCent); << 261 numberofElements = 13); 245 spinalDisc->AddElement(elN, 2.20 * perCent); << 262 brainTissue->AddElement(elH, 11.0667*perCent); 246 spinalDisc->AddElement(elO, 74.40 * perCent) << 263 brainTissue->AddElement(elC, 12.542*perCent); 247 spinalDisc->AddElement(elNa, 0.50 * perCent) << 264 brainTissue->AddElement(elN, 1.328*perCent); 248 spinalDisc->AddElement(elP, 2.20 * perCent); << 265 brainTissue->AddElement(elO, 73.7723*perCent); 249 spinalDisc->AddElement(elS, 0.90 * perCent); << 266 brainTissue->AddElement(elNa, 0.1840*perCent); 250 spinalDisc->AddElement(elCl, 0.30 * perCent) << 267 brainTissue->AddElement(elMg, 0.015*perCent); 251 << 268 brainTissue->AddElement(elP, 0.356*perCent); 252 // Water << 269 brainTissue->AddElement(elS, 0.177*perCent); 253 G4Material* water = new G4Material("Water", << 270 brainTissue->AddElement(elCl, 0.236*perCent); 254 water->AddElement(elH, 0.112); << 271 brainTissue->AddElement(elK, 0.31*perCent); 255 water->AddElement(elO, 0.888); << 272 brainTissue->AddElement(elCa, 0.009*perCent); 256 << 273 brainTissue->AddElement(elFe, 0.005*perCent); 257 // Muscle << 274 brainTissue->AddElement(elZn, 0.001*perCent); 258 G4Material* muscle = new G4Material("Muscle" << 275 259 muscle->AddElement(elH, 0.102); << 276 260 muscle->AddElement(elC, 0.143); << 277 // Breast 261 muscle->AddElement(elN, 0.034); << 278 G4Material* breast = new G4Material( "Breast", 262 muscle->AddElement(elO, 0.710); << 279 density = 0.990*g/cm3, 263 muscle->AddElement(elNa, 0.001); << 280 numberofElements = 8 ); 264 muscle->AddElement(elP, 0.002); << 281 breast->AddElement(elH,0.109); 265 muscle->AddElement(elS, 0.003); << 282 breast->AddElement(elC,0.506); 266 muscle->AddElement(elCl, 0.001); << 283 breast->AddElement(elN,0.023); 267 muscle->AddElement(elK, 0.004); << 284 breast->AddElement(elO,0.358); 268 << 285 breast->AddElement(elNa,0.001); 269 // Liver << 286 breast->AddElement(elP,0.001); 270 G4Material* liver = new G4Material("Liver", << 287 breast->AddElement(elS,0.001); 271 liver->AddElement(elH, 0.102); << 288 breast->AddElement(elCl,0.001); 272 liver->AddElement(elC, 0.139); << 289 273 liver->AddElement(elN, 0.030); << 290 // Spinal Disc 274 liver->AddElement(elO, 0.716); << 291 G4Material* spinalDisc = new G4Material ("SpinalDisc", 1.10 * g/cm3, 275 liver->AddElement(elNa, 0.002); << 292 numberofElements = 8); 276 liver->AddElement(elP, 0.003); << 293 spinalDisc->AddElement(elH, 9.60*perCent); 277 liver->AddElement(elS, 0.003); << 294 spinalDisc->AddElement(elC, 9.90*perCent); 278 liver->AddElement(elCl, 0.002); << 295 spinalDisc->AddElement(elN, 2.20*perCent); 279 liver->AddElement(elK, 0.003); << 296 spinalDisc->AddElement(elO, 74.40*perCent); 280 << 297 spinalDisc->AddElement(elNa, 0.50*perCent); 281 // Tooth Dentin << 298 spinalDisc->AddElement(elP, 2.20*perCent); 282 G4Material* toothDentin = new G4Material("To << 299 spinalDisc->AddElement(elS, 0.90*perCent); 283 toothDentin->AddElement(elH, 2.67 * perCent) << 300 spinalDisc->AddElement(elCl, 0.30*perCent); 284 toothDentin->AddElement(elC, 12.77 * perCent << 301 285 toothDentin->AddElement(elN, 4.27 * perCent) << 302 286 toothDentin->AddElement(elO, 40.40 * perCent << 303 // Water 287 toothDentin->AddElement(elNa, 0.65 * perCent << 304 G4Material* water = new G4Material( "Water", 288 toothDentin->AddElement(elMg, 0.59 * perCent << 305 density = 1.0*g/cm3, 289 toothDentin->AddElement(elP, 11.86 * perCent << 306 numberofElements = 2 ); 290 toothDentin->AddElement(elCl, 0.04 * perCent << 307 water->AddElement(elH,0.112); 291 toothDentin->AddElement(elCa, 26.74 * perCen << 308 water->AddElement(elO,0.888); 292 toothDentin->AddElement(elZn, 0.01 * perCent << 309 293 << 310 // Muscle 294 // Trabecular Bone << 311 G4Material* muscle = new G4Material( "Muscle", 295 G4Material* trabecularBone = << 312 density = 1.061*g/cm3, 296 new G4Material("TrabecularBone", density = << 313 numberofElements = 9 ); 297 trabecularBone->AddElement(elH, 0.085); << 314 muscle->AddElement(elH,0.102); 298 trabecularBone->AddElement(elC, 0.404); << 315 muscle->AddElement(elC,0.143); 299 trabecularBone->AddElement(elN, 0.058); << 316 muscle->AddElement(elN,0.034); 300 trabecularBone->AddElement(elO, 0.367); << 317 muscle->AddElement(elO,0.710); 301 trabecularBone->AddElement(elNa, 0.001); << 318 muscle->AddElement(elNa,0.001); 302 trabecularBone->AddElement(elMg, 0.001); << 319 muscle->AddElement(elP,0.002); 303 trabecularBone->AddElement(elP, 0.034); << 320 muscle->AddElement(elS,0.003); 304 trabecularBone->AddElement(elS, 0.002); << 321 muscle->AddElement(elCl,0.001); 305 trabecularBone->AddElement(elCl, 0.002); << 322 muscle->AddElement(elK,0.004); 306 trabecularBone->AddElement(elK, 0.001); << 323 307 trabecularBone->AddElement(elCa, 0.044); << 324 // Liver 308 trabecularBone->AddElement(elFe, 0.001); << 325 G4Material* liver = new G4Material( "Liver", 309 << 326 density = 1.071*g/cm3, 310 // Trabecular bone used in the DICOM Head << 327 numberofElements = 9); 311 << 328 liver->AddElement(elH,0.102); 312 G4Material* trabecularBone_head = << 329 liver->AddElement(elC,0.139); 313 new G4Material("TrabecularBone_HEAD", 1.18 << 330 liver->AddElement(elN,0.030); 314 trabecularBone_head->AddElement(elH, 8.50 * << 331 liver->AddElement(elO,0.716); 315 trabecularBone_head->AddElement(elC, 40.40 * << 332 liver->AddElement(elNa,0.002); 316 trabecularBone_head->AddElement(elN, 2.80 * << 333 liver->AddElement(elP,0.003); 317 trabecularBone_head->AddElement(elO, 36.70 * << 334 liver->AddElement(elS,0.003); 318 trabecularBone_head->AddElement(elNa, 0.10 * << 335 liver->AddElement(elCl,0.002); 319 trabecularBone_head->AddElement(elMg, 0.10 * << 336 liver->AddElement(elK,0.003); 320 trabecularBone_head->AddElement(elP, 3.40 * << 337 321 trabecularBone_head->AddElement(elS, 0.20 * << 338 // Tooth Dentin 322 trabecularBone_head->AddElement(elCl, 0.20 * << 339 G4Material* toothDentin = new G4Material ("ToothDentin", 2.14 * g/cm3, 323 trabecularBone_head->AddElement(elK, 0.10 * << 340 numberofElements = 10); 324 trabecularBone_head->AddElement(elCa, 7.40 * << 341 toothDentin->AddElement(elH, 2.67*perCent); 325 trabecularBone_head->AddElement(elFe, 0.10 * << 342 toothDentin->AddElement(elC, 12.77*perCent); 326 << 343 toothDentin->AddElement(elN, 4.27*perCent); 327 // Dense Bone << 344 toothDentin->AddElement(elO, 40.40*perCent); 328 G4Material* denseBone = << 345 toothDentin->AddElement(elNa, 0.65*perCent); 329 new G4Material("DenseBone", density = 1.57 << 346 toothDentin->AddElement(elMg, 0.59*perCent); 330 denseBone->AddElement(elH, 0.056); << 347 toothDentin->AddElement(elP, 11.86*perCent); 331 denseBone->AddElement(elC, 0.235); << 348 toothDentin->AddElement(elCl, 0.04*perCent); 332 denseBone->AddElement(elN, 0.050); << 349 toothDentin->AddElement(elCa, 26.74*perCent); 333 denseBone->AddElement(elO, 0.434); << 350 toothDentin->AddElement(elZn, 0.01*perCent); 334 denseBone->AddElement(elNa, 0.001); << 351 335 denseBone->AddElement(elMg, 0.001); << 352 336 denseBone->AddElement(elP, 0.072); << 353 // Trabecular Bone 337 denseBone->AddElement(elS, 0.003); << 354 G4Material* trabecularBone = new G4Material( "TrabecularBone", density = 1.159*g/cm3, 338 denseBone->AddElement(elCl, 0.001); << 355 numberofElements = 12 ); 339 denseBone->AddElement(elK, 0.001); << 356 trabecularBone->AddElement(elH,0.085); 340 denseBone->AddElement(elCa, 0.146); << 357 trabecularBone->AddElement(elC,0.404); 341 << 358 trabecularBone->AddElement(elN,0.058); 342 // Cortical Bone (ICRP - NIST) << 359 trabecularBone->AddElement(elO,0.367); 343 G4Material* corticalBone = new G4Material("C << 360 trabecularBone->AddElement(elNa,0.001); 344 corticalBone->AddElement(elH, 4.7234 * perCe << 361 trabecularBone->AddElement(elMg,0.001); 345 corticalBone->AddElement(elC, 14.4330 * perC << 362 trabecularBone->AddElement(elP,0.034); 346 corticalBone->AddElement(elN, 4.199 * perCen << 363 trabecularBone->AddElement(elS,0.002); 347 corticalBone->AddElement(elO, 44.6096 * perC << 364 trabecularBone->AddElement(elCl,0.002); 348 corticalBone->AddElement(elMg, 0.22 * perCen << 365 trabecularBone->AddElement(elK,0.001); 349 corticalBone->AddElement(elP, 10.497 * perCe << 366 trabecularBone->AddElement(elCa,0.044); 350 corticalBone->AddElement(elS, 0.315 * perCen << 367 trabecularBone->AddElement(elFe,0.001); 351 corticalBone->AddElement(elCa, 20.993 * perC << 368 352 corticalBone->AddElement(elZn, 0.01 * perCen << 369 // Trabecular bone used in the DICOM Head 353 << 370 354 // Tooth enamel << 371 G4Material* trabecularBone_head = new G4Material ("TrabecularBone_HEAD", 355 G4Material* toothEnamel = new G4Material("To << 372 1.18 * g/cm3, 356 toothEnamel->AddElement(elH, 0.95 * perCent) << 373 numberofElements = 12); 357 toothEnamel->AddElement(elC, 1.11 * perCent) << 374 trabecularBone_head->AddElement(elH, 8.50*perCent); 358 toothEnamel->AddElement(elN, 0.23 * perCent) << 375 trabecularBone_head->AddElement(elC, 40.40*perCent); 359 toothEnamel->AddElement(elO, 41.66 * perCent << 376 trabecularBone_head->AddElement(elN, 2.80*perCent); 360 toothEnamel->AddElement(elNa, 0.79 * perCent << 377 trabecularBone_head->AddElement(elO, 36.70*perCent); 361 toothEnamel->AddElement(elMg, 0.23 * perCent << 378 trabecularBone_head->AddElement(elNa, 0.10*perCent); 362 toothEnamel->AddElement(elP, 18.71 * perCent << 379 trabecularBone_head->AddElement(elMg, 0.10*perCent); 363 toothEnamel->AddElement(elCl, 0.34 * perCent << 380 trabecularBone_head->AddElement(elP, 3.40*perCent); 364 toothEnamel->AddElement(elCa, 35.97 * perCen << 381 trabecularBone_head->AddElement(elS, 0.20*perCent); 365 toothEnamel->AddElement(elZn, 0.02 * perCent << 382 trabecularBone_head->AddElement(elCl, 0.20*perCent); >> 383 trabecularBone_head->AddElement(elK, 0.10*perCent); >> 384 trabecularBone_head->AddElement(elCa, 7.40*perCent); >> 385 trabecularBone_head->AddElement(elFe, 0.10*perCent); >> 386 >> 387 // Dense Bone >> 388 G4Material* denseBone = new G4Material( "DenseBone", >> 389 density = 1.575*g/cm3, >> 390 numberofElements = 11 ); >> 391 denseBone->AddElement(elH,0.056); >> 392 denseBone->AddElement(elC,0.235); >> 393 denseBone->AddElement(elN,0.050); >> 394 denseBone->AddElement(elO,0.434); >> 395 denseBone->AddElement(elNa,0.001); >> 396 denseBone->AddElement(elMg,0.001); >> 397 denseBone->AddElement(elP,0.072); >> 398 denseBone->AddElement(elS,0.003); >> 399 denseBone->AddElement(elCl,0.001); >> 400 denseBone->AddElement(elK,0.001); >> 401 denseBone->AddElement(elCa,0.146); >> 402 >> 403 // Cortical Bone (ICRP - NIST) >> 404 G4Material* corticalBone = new G4Material ("CorticalBone", 1.85 * g/cm3, >> 405 numberofElements = 9); >> 406 corticalBone->AddElement(elH, 4.7234*perCent); >> 407 corticalBone->AddElement(elC, 14.4330*perCent); >> 408 corticalBone->AddElement(elN, 4.199*perCent); >> 409 corticalBone->AddElement(elO, 44.6096*perCent); >> 410 corticalBone->AddElement(elMg, 0.22*perCent); >> 411 corticalBone->AddElement(elP, 10.497*perCent); >> 412 corticalBone->AddElement(elS, 0.315*perCent); >> 413 corticalBone->AddElement(elCa, 20.993*perCent); >> 414 corticalBone->AddElement(elZn, 0.01*perCent); >> 415 >> 416 >> 417 // Tooth enamel >> 418 G4Material* toothEnamel = new G4Material ("ToothEnamel", 2.89 * g/cm3, >> 419 numberofElements = 10); >> 420 toothEnamel->AddElement(elH, 0.95*perCent); >> 421 toothEnamel->AddElement(elC, 1.11*perCent); >> 422 toothEnamel->AddElement(elN, 0.23*perCent); >> 423 toothEnamel->AddElement(elO,41.66*perCent); >> 424 toothEnamel->AddElement(elNa, 0.79*perCent); >> 425 toothEnamel->AddElement(elMg, 0.23*perCent); >> 426 toothEnamel->AddElement(elP, 18.71*perCent); >> 427 toothEnamel->AddElement(elCl, 0.34*perCent); >> 428 toothEnamel->AddElement(elCa, 35.97*perCent); >> 429 toothEnamel->AddElement(elZn, 0.02*perCent); 366 430 367 #ifdef DICOM_USE_HEAD 431 #ifdef DICOM_USE_HEAD 368 //----- Put the materials in a vector HEAD P 432 //----- Put the materials in a vector HEAD PHANTOM 369 fOriginalMaterials.push_back(fAir); // 0.00 << 433 fOriginalMaterials.push_back(fAir); //0.00129 g/cm3 370 fOriginalMaterials.push_back(softTissue); / << 434 fOriginalMaterials.push_back(softTissue); // 1.055 g/cm3 371 fOriginalMaterials.push_back(brainTissue); << 435 fOriginalMaterials.push_back(brainTissue); // 1.07 g/cm3 372 fOriginalMaterials.push_back(spinalDisc); / << 436 fOriginalMaterials.push_back(spinalDisc); // 1.10 g/cm3 373 fOriginalMaterials.push_back(trabecularBone_ << 437 fOriginalMaterials.push_back(trabecularBone_head); // 1.13 g/cm3 374 fOriginalMaterials.push_back(toothDentin); << 438 fOriginalMaterials.push_back(toothDentin); // 1.66 g/cm3 375 fOriginalMaterials.push_back(corticalBone); 439 fOriginalMaterials.push_back(corticalBone); // 1.75 g/cm3 376 fOriginalMaterials.push_back(toothEnamel); << 440 fOriginalMaterials.push_back(toothEnamel); // 2.04 g/cm3 377 G4cout << "The materials of the DICOM Head h << 441 G4cout << "The materials of the DICOM Head have been used" << G4endl; 378 #else 442 #else 379 fOriginalMaterials.push_back(fAir); // rho << 443 fOriginalMaterials.push_back(fAir); // rho = 0.00129 380 fOriginalMaterials.push_back(lunginhale); / << 444 fOriginalMaterials.push_back(lunginhale); // rho = 0.217 381 fOriginalMaterials.push_back(lungexhale); / << 445 fOriginalMaterials.push_back(lungexhale); // rho = 0.508 382 fOriginalMaterials.push_back(adiposeTissue); << 446 fOriginalMaterials.push_back(adiposeTissue); // rho = 0.967 383 fOriginalMaterials.push_back(breast); // rh << 447 fOriginalMaterials.push_back(breast ); // rho = 0.990 384 fOriginalMaterials.push_back(water); // rho << 448 fOriginalMaterials.push_back(water); // rho = 1.018 385 fOriginalMaterials.push_back(muscle); // rh << 449 fOriginalMaterials.push_back(muscle); // rho = 1.061 386 fOriginalMaterials.push_back(liver); // rho << 450 fOriginalMaterials.push_back(liver); // rho = 1.071 387 fOriginalMaterials.push_back(trabecularBone) << 451 fOriginalMaterials.push_back(trabecularBone); // rho = 1.159 - HEAD PHANTOM 388 fOriginalMaterials.push_back(denseBone); // << 452 fOriginalMaterials.push_back(denseBone); // rho = 1.575 389 G4cout << "Default materials of the DICOM Ex << 453 G4cout << "The default materials of the DICOM Extended examples have been used" >> 454 << G4endl; 390 #endif 455 #endif 391 } 456 } 392 457 393 //....oooOO0OOooo........oooOO0OOooo........oo 458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 394 void DicomDetectorConstruction::ReadPhantomDat 459 void DicomDetectorConstruction::ReadPhantomDataNew() 395 { 460 { 396 #ifdef G4_DCMTK 461 #ifdef G4_DCMTK 397 G4String fileName = DicomFileMgr::GetInstanc 462 G4String fileName = DicomFileMgr::GetInstance()->GetFileOutName(); 398 << 463 399 std::ifstream fin(fileName); 464 std::ifstream fin(fileName); 400 std::vector<G4String> wl; 465 std::vector<G4String> wl; 401 G4int nMaterials; 466 G4int nMaterials; 402 fin >> nMaterials; 467 fin >> nMaterials; 403 G4String mateName; 468 G4String mateName; 404 G4int nmate; 469 G4int nmate; 405 for (G4int ii = 0; ii < nMaterials; ++ii) { << 470 for( G4int ii = 0; ii < nMaterials; ii++ ){ 406 fin >> nmate; 471 fin >> nmate; 407 fin >> mateName; 472 fin >> mateName; 408 if (mateName[0] == '"' && mateName[mateNam << 473 if( mateName[0] == '"' && mateName[mateName.length()-1] == '"' ) { 409 mateName = mateName.substr(1, mateName.l << 474 mateName = mateName.substr(1,mateName.length()-2); 410 } 475 } 411 G4cout << "GmReadPhantomG4Geometry::ReadPh << 476 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData reading nmate " >> 477 << ii << " = " << nmate 412 << " mate " << mateName << G4endl; 478 << " mate " << mateName << G4endl; 413 if (ii != nmate) << 479 if( ii != nmate ) 414 G4Exception("GmReadPhantomG4Geometry::Re << 480 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", 415 FatalErrorInArgument, << 481 "Wrong argument", 416 "Material number should be i << 482 FatalErrorInArgument, >> 483 "Material number should be in increasing order:wrong material number"); 417 484 418 G4Material* mate = 0; 485 G4Material* mate = 0; 419 const G4MaterialTable* matTab = G4Material 486 const G4MaterialTable* matTab = G4Material::GetMaterialTable(); 420 for (auto matite = matTab->cbegin(); matit << 487 std::vector<G4Material*>::const_iterator matite; 421 if ((*matite)->GetName() == mateName) { << 488 for( matite = matTab->begin(); matite != matTab->end(); ++matite ) { >> 489 if( (*matite)->GetName() == mateName ) { 422 mate = *matite; 490 mate = *matite; 423 } 491 } 424 } 492 } 425 if (mate == 0) { << 493 if( mate == 0 ) { 426 mate = G4NistManager::Instance()->FindOr 494 mate = G4NistManager::Instance()->FindOrBuildMaterial(mateName); 427 } 495 } 428 if (!mate) << 496 if( !mate ) G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", 429 G4Exception("GmReadPhantomG4Geometry::Re << 497 "Wrong argument", 430 FatalErrorInArgument, ("Mate << 498 FatalErrorInArgument, 431 fPhantomMaterialsOriginal[nmate] = mate; << 499 ("Material not found" + mateName).c_str()); >> 500 thePhantomMaterialsOriginal[nmate] = mate; 432 } 501 } 433 502 434 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxels << 503 fin >> fNVoxelX >> fNVoxelY >> fNVoxelZ; 435 G4cout << "GmReadPhantomG4Geometry::ReadPhan << 504 G4cout << "GmReadPhantomG4Geometry::ReadPhantomData fNVoxel X/Y/Z " 436 << fNoVoxelsY << " " << fNoVoxelsZ << << 505 << fNVoxelX << " " >> 506 << fNVoxelY << " " << fNVoxelZ << G4endl; 437 fin >> fMinX >> fMaxX; 507 fin >> fMinX >> fMaxX; 438 fin >> fMinY >> fMaxY; 508 fin >> fMinY >> fMaxY; 439 fin >> fMinZ >> fMaxZ; 509 fin >> fMinZ >> fMaxZ; 440 fVoxelHalfDimX = (fMaxX - fMinX) / fNoVoxels << 510 fVoxelHalfDimX = (fMaxX-fMinX)/fNVoxelX/2.; 441 fVoxelHalfDimY = (fMaxY - fMinY) / fNoVoxels << 511 fVoxelHalfDimY = (fMaxY-fMinY)/fNVoxelY/2.; 442 fVoxelHalfDimZ = (fMaxZ - fMinZ) / fNoVoxels << 512 fVoxelHalfDimZ = (fMaxZ-fMinZ)/fNVoxelZ/2.; 443 # ifdef G4VERBOSE << 513 #ifdef G4VERBOSE 444 G4cout << " Extension in X " << fMinX << " " << 514 G4cout << " Extension in X " << fMinX << " " << fMaxX << G4endl 445 << " " << fMaxY << G4endl << " Extens << 515 << " Extension in Y " << fMinY << " " << fMaxY << G4endl 446 # endif << 516 << " Extension in Z " << fMinZ << " " << fMaxZ << G4endl; 447 << 517 #endif 448 fMateIDs = new size_t[fNoVoxelsX * fNoVoxels << 518 449 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) { << 519 fMateIDs = new size_t[fNVoxelX*fNVoxelY*fNVoxelZ]; 450 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) << 520 for( G4int iz = 0; iz < fNVoxelZ; iz++ ) { 451 for (G4int ix = 0; ix < fNoVoxelsX; ++ix << 521 for( G4int iy = 0; iy < fNVoxelY; iy++ ) { >> 522 for( G4int ix = 0; ix < fNVoxelX; ix++ ) { 452 G4int mateID; 523 G4int mateID; 453 fin >> mateID; << 524 fin >> mateID; 454 G4int nnew = ix + (iy)*fNoVoxelsX + (i << 525 G4int nnew = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY; 455 if (mateID < 0 || mateID >= nMaterials << 526 if( mateID < 0 || mateID >= nMaterials ) { 456 G4Exception("GmReadPhantomG4Geometry << 527 G4Exception("GmReadPhantomG4Geometry::ReadPhantomData", >> 528 "Wrong index in phantom file", 457 FatalException, 529 FatalException, 458 G4String("It should be b 530 G4String("It should be between 0 and " 459 + G4UIcommand:: << 531 + G4UIcommand::ConvertToString(nMaterials-1) 460 + G4UIcommand:: << 532 + ", while it is " 461 .c_str()); << 533 + G4UIcommand::ConvertToString(mateID)).c_str()); 462 } 534 } 463 fMateIDs[nnew] = mateID; 535 fMateIDs[nnew] = mateID; 464 } 536 } 465 } 537 } 466 } 538 } 467 539 468 ReadVoxelDensities(fin); << 540 ReadVoxelDensities( fin ); 469 541 470 fin.close(); 542 fin.close(); 471 #endif 543 #endif 472 } << 473 544 474 //....oooOO0OOooo........oooOO0OOooo........oo << 545 } 475 void DicomDetectorConstruction::ReadVoxelDensi << 546 void DicomDetectorConstruction::ReadVoxelDensities( std::ifstream& fin ) 476 { 547 { 477 G4String stemp; 548 G4String stemp; 478 std::map<G4int, std::pair<G4double, G4double << 549 std::map<G4int, std::pair<G4double,G4double> > densiMinMax; 479 std::map<G4int, std::pair<G4double, G4double << 550 std::map<G4int, std::pair<G4double,G4double> >::iterator mpite; 480 for (G4int ii = 0; ii < G4int(fPhantomMateri << 551 for( size_t ii = 0; ii < thePhantomMaterialsOriginal.size(); ii++ ){ 481 densiMinMax[ii] = std::pair<G4double, G4do << 552 densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX); 482 } 553 } 483 554 484 char* part = std::getenv("DICOM_CHANGE_MATER << 555 char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" ); 485 G4double densityDiff = -1.; 556 G4double densityDiff = -1.; 486 if (part) densityDiff = G4UIcommand::Convert << 557 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part); 487 558 488 std::map<G4int, G4double> densityDiffs; << 559 std::map<G4int,G4double> densityDiffs; 489 for (G4int ii = 0; ii < G4int(fPhantomMateri << 560 for( size_t ii = 0; ii < thePhantomMaterialsOriginal.size(); ii++ ){ 490 densityDiffs[ii] = densityDiff; // curren << 561 densityDiffs[ii] = densityDiff; //currently all materials with same step 491 } 562 } 492 // densityDiffs[0] = 0.0001; //air 563 // densityDiffs[0] = 0.0001; //air 493 564 494 //--- Calculate the average material density 565 //--- Calculate the average material density for each material/density bin 495 std::map<std::pair<G4Material*, G4int>, matI << 566 std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens; 496 567 497 //---- Read the material densities 568 //---- Read the material densities 498 G4double dens; 569 G4double dens; 499 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) { << 570 for( G4int iz = 0; iz < fNVoxelZ; iz++ ) { 500 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) << 571 for( G4int iy = 0; iy < fNVoxelY; iy++ ) { 501 for (G4int ix = 0; ix < fNoVoxelsX; ++ix << 572 for( G4int ix = 0; ix < fNVoxelX; ix++ ) { 502 fin >> dens; << 573 fin >> dens; 503 G4int copyNo = ix + (iy)*fNoVoxelsX + << 574 // G4cout << ix << " " << iy << " " << iz 504 << 575 // << " density " << dens << G4endl; 505 if (densityDiff != -1.) continue; << 576 G4int copyNo = ix + (iy)*fNVoxelX + (iz)*fNVoxelX*fNVoxelY; 506 << 577 507 //--- store the minimum and maximum de << 578 if( densityDiff != -1. ) continue; 508 mpite = densiMinMax.find(G4int(fMateID << 579 509 if (dens < (*mpite).second.first) (*mp << 580 //--- store the minimum and maximum density for each material (just for printin 510 if (dens > (*mpite).second.second) (*m << 581 mpite = densiMinMax.find( fMateIDs[copyNo] ); >> 582 if( dens < (*mpite).second.first ) (*mpite).second.first = dens; >> 583 if( dens > (*mpite).second.second ) (*mpite).second.second = dens; 511 //--- Get material from original list 584 //--- Get material from original list of material in file 512 G4int mateID = G4int(fMateIDs[copyNo]) << 585 int mateID = fMateIDs[copyNo]; 513 std::map<G4int, G4Material*>::const_it << 586 std::map<G4int,G4Material*>::const_iterator imite = 514 << 587 thePhantomMaterialsOriginal.find(mateID); >> 588 // G4cout << copyNo << " mateID " << mateID << G4endl; 515 //--- Check if density is equal to the 589 //--- Check if density is equal to the original material density 516 if (std::fabs(dens - (*imite).second-> << 590 if(std::fabs(dens - (*imite).second->GetDensity()/CLHEP::g*CLHEP::cm3 ) < 1.e-9 ) continue; 517 continue; << 591 >> 592 //--- Build material name with thePhantomMaterialsOriginal name + density >> 593 //float densityBin = densityDiffs[mateID] * (G4int(dens/densityDiffs[mateID])+0.5); >> 594 G4int densityBin = (G4int(dens/densityDiffs[mateID])); 518 595 519 //--- Build material name with fPhanto << 596 G4String mateName = (*imite).second->GetName()+G4UIcommand::ConvertToString(densityBin); 520 G4int densityBin = (G4int(dens / densi << 521 << 522 G4String mateName = (*imite).second->G << 523 //--- Look if it is the first voxel wi 597 //--- Look if it is the first voxel with this material/densityBin 524 std::pair<G4Material*, G4int> matdens( << 598 std::pair<G4Material*,G4int> matdens((*imite).second, densityBin ); 525 599 526 auto mppite = newMateDens.find(matdens << 600 std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite = 527 if (mppite != newMateDens.cend()) { << 601 newMateDens.find( matdens ); >> 602 if( mppite != newMateDens.end() ){ 528 matInfo* mi = (*mppite).second; 603 matInfo* mi = (*mppite).second; 529 mi->fSumdens += dens; 604 mi->fSumdens += dens; 530 mi->fNvoxels++; 605 mi->fNvoxels++; 531 fMateIDs[copyNo] = fPhantomMaterials << 606 fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId; 532 } << 607 } else { 533 else { << 534 matInfo* mi = new matInfo; 608 matInfo* mi = new matInfo; 535 mi->fSumdens = dens; 609 mi->fSumdens = dens; 536 mi->fNvoxels = 1; 610 mi->fNvoxels = 1; 537 mi->fId = G4int(newMateDens.size() + << 611 mi->fId = newMateDens.size()+1; 538 newMateDens[matdens] = mi; 612 newMateDens[matdens] = mi; 539 fMateIDs[copyNo] = fPhantomMaterials << 613 fMateIDs[copyNo] = thePhantomMaterialsOriginal.size()-1 + mi->fId; 540 } 614 } 541 } 615 } 542 } 616 } 543 } 617 } 544 618 545 if (densityDiff != -1.) { << 619 if( densityDiff != -1. ) { 546 for (mpite = densiMinMax.begin(); mpite != << 620 for( mpite = densiMinMax.begin(); mpite != densiMinMax.end(); mpite++ ){ 547 #ifdef G4VERBOSE 621 #ifdef G4VERBOSE 548 G4cout << "DicomDetectorConstruction::Re << 622 G4cout << "DicomDetectorConstruction::ReadVoxelDensities ORIG MATERIALS DENSITY " 549 << " ORIG MATERIALS DENSITY " << << 623 << (*mpite).first << " MIN " << (*mpite).second.first << " MAX " 550 << " MAX " << (*mpite).second.sec << 624 << (*mpite).second.second << G4endl; 551 #endif 625 #endif 552 } 626 } 553 } 627 } 554 628 555 //----- Build the list of phantom materials 629 //----- Build the list of phantom materials that go to Parameterisation 556 //--- Add original materials 630 //--- Add original materials 557 for (auto mimite = fPhantomMaterialsOriginal << 631 std::map<G4int,G4Material*>::const_iterator mimite; 558 ++mimite) << 632 for( mimite = thePhantomMaterialsOriginal.begin(); mimite != thePhantomMaterialsOriginal.end(); 559 { << 633 mimite++ ){ 560 fMaterials.push_back((*mimite).second); << 634 fMaterials.push_back( (*mimite).second ); 561 } 635 } 562 // << 636 // 563 //---- Build and add new materials 637 //---- Build and add new materials 564 for (auto mppite = newMateDens.cbegin(); mpp << 638 std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite; 565 G4double averdens = (*mppite).second->fSum << 639 for( mppite= newMateDens.begin(); mppite != newMateDens.end(); mppite++ ){ 566 G4double saverdens = G4int(1000.001 * aver << 640 G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels; >> 641 G4double saverdens = G4int(1000.001*averdens)/1000.; 567 #ifndef G4VERBOSE 642 #ifndef G4VERBOSE 568 G4cout << "DicomDetectorConstruction::Read << 643 G4cout << "DicomDetectorConstruction::ReadVoxelDensities AVER DENS " 569 << saverdens << " -> " << G4int(100 << 644 << averdens << " -> " 570 << " " << G4int(1000 * averdens) / << 645 << saverdens << " -> " << G4int(1000*averdens) << " " >> 646 << G4int(1000*averdens)/1000 >> 647 << " " << G4int(1000*averdens)/1000. << G4endl; 571 #endif 648 #endif 572 649 573 G4String mateName = << 650 G4String mateName = ((*mppite).first).first->GetName() + "_" 574 ((*mppite).first).first->GetName() + "_" << 651 + G4UIcommand::ConvertToString(saverdens); 575 fMaterials.push_back( << 652 fMaterials.push_back( BuildMaterialWithChangingDensity( 576 BuildMaterialWithChangingDensity((*mppit << 653 (*mppite).first.first, averdens, mateName ) ); 577 } 654 } >> 655 578 } 656 } 579 657 580 //....oooOO0OOooo........oooOO0OOooo........oo 658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....... 581 void DicomDetectorConstruction::ReadPhantomDat 659 void DicomDetectorConstruction::ReadPhantomData() 582 { 660 { 583 G4String dataFile = DicomHandler::GetDicomDa << 661 #ifdef DICOM_USE_HEAD >> 662 G4String path = getenv("DICOM_PATH"); >> 663 G4String dataFile = path+"/Data.dat"; >> 664 #else >> 665 G4String dataFile = "Data.dat"; >> 666 #endif 584 std::ifstream finDF(dataFile.c_str()); 667 std::ifstream finDF(dataFile.c_str()); 585 G4String fname; 668 G4String fname; 586 669 587 if (finDF.good() != 1) { << 670 if(finDF.good() != 1 ) 588 G4String descript = "Problem reading data << 671 { 589 G4Exception(" DicomDetectorConstruction::R << 672 G4String descript = "Problem reading data file: "+dataFile; >> 673 G4Exception(" DicomDetectorConstruction::ReadPhantomData"," ", >> 674 FatalException,descript); 590 } 675 } 591 676 592 G4int compression; 677 G4int compression; 593 finDF >> compression; // not used here << 678 finDF >> compression; // not used here 594 finDF >> fNoFiles; 679 finDF >> fNoFiles; 595 680 596 for (G4int i = 0; i < fNoFiles; ++i) { << 681 for(G4int i = 0; i < fNoFiles; i++ ) { 597 finDF >> fname; << 598 682 599 //--- Read one data file << 683 finDF >> fname; 600 fname += ".g4dcm"; << 601 684 602 ReadPhantomDataFile(fname); << 685 //--- Read one data file 603 } << 686 fname += ".g4dcm"; >> 687 >> 688 ReadPhantomDataFile(fname); 604 689 605 //----- Merge data headers << 606 MergeZSliceHeaders(); << 607 finDF.close(); << 608 } 690 } 609 691 >> 692 //----- Merge data headers >> 693 MergeZSliceHeaders(); >> 694 finDF.close(); >> 695 } 610 //....oooOO0OOooo........oooOO0OOooo........oo 696 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........ 611 void DicomDetectorConstruction::ReadPhantomDat 697 void DicomDetectorConstruction::ReadPhantomDataFile(const G4String& fname) 612 { 698 { 613 G4cout << " DicomDetectorConstruction::ReadP << 699 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " 614 << G4endl; // GDEB << 700 << fname << G4endl; //GDEB 615 701 616 #ifdef G4VERBOSE 702 #ifdef G4VERBOSE 617 G4cout << " DicomDetectorConstruction::ReadP << 703 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " >> 704 << fname << G4endl; 618 #endif 705 #endif 619 << 706 620 std::ifstream fin(fname.c_str(), std::ios_ba 707 std::ifstream fin(fname.c_str(), std::ios_base::in); 621 if (!fin.is_open()) { << 708 if( !fin.is_open() ) { 622 G4Exception("DicomDetectorConstruction::Re << 709 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", 623 G4String("File not found " + f << 710 "", >> 711 FatalErrorInArgument, >> 712 G4String("File not found " + fname ).c_str()); 624 } 713 } 625 //----- Define density differences (maximum 714 //----- Define density differences (maximum density difference to create 626 // a new material) 715 // a new material) 627 char* part = std::getenv("DICOM_CHANGE_MATER << 716 char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" ); 628 G4double densityDiff = -1.; 717 G4double densityDiff = -1.; 629 if (part) densityDiff = G4UIcommand::Convert << 718 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part); 630 if (densityDiff != -1.) { << 719 if( densityDiff != -1. ) { 631 for (unsigned int ii = 0; ii < fOriginalMa << 720 for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){ 632 fDensityDiffs[ii] = densityDiff; // cur << 721 fDensityDiffs[ii] = densityDiff; //currently all materials with 633 // same difference 722 // same difference 634 } 723 } 635 } << 724 }else { 636 else { << 725 if( fMaterials.size() == 0 ) { // do it only for first slice 637 if (fMaterials.size() == 0) { // do it on << 726 for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){ 638 for (unsigned int ii = 0; ii < fOriginal << 727 fMaterials.push_back( fOriginalMaterials[ii] ); 639 fMaterials.push_back(fOriginalMaterial << 640 } 728 } 641 } 729 } 642 } 730 } 643 << 731 644 //----- Read data header 732 //----- Read data header 645 DicomPhantomZSliceHeader* sliceHeader = new << 733 DicomPhantomZSliceHeader* sliceHeader = new DicomPhantomZSliceHeader( fin ); 646 fZSliceHeaders.push_back(sliceHeader); << 734 fZSliceHeaders.push_back( sliceHeader ); 647 << 735 648 //----- Read material indices 736 //----- Read material indices 649 G4int nVoxels = sliceHeader->GetNoVoxels(); 737 G4int nVoxels = sliceHeader->GetNoVoxels(); 650 << 738 651 //--- If first slice, initiliaze fMateIDs 739 //--- If first slice, initiliaze fMateIDs 652 if (fZSliceHeaders.size() == 1) { << 740 if( fZSliceHeaders.size() == 1 ) { 653 // fMateIDs = new unsigned int[fNoFiles*nV << 741 //fMateIDs = new unsigned int[fNoFiles*nVoxels]; 654 fMateIDs = new size_t[fNoFiles * nVoxels]; << 742 fMateIDs = new size_t[fNoFiles*nVoxels]; >> 743 655 } 744 } 656 << 745 657 unsigned int mateID; 746 unsigned int mateID; 658 // number of voxels from previously read sli 747 // number of voxels from previously read slices 659 G4int voxelCopyNo = G4int((fZSliceHeaders.si << 748 G4int voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels; 660 for (G4int ii = 0; ii < nVoxels; ++ii, voxel << 749 for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){ 661 fin >> mateID; 750 fin >> mateID; 662 fMateIDs[voxelCopyNo] = mateID; 751 fMateIDs[voxelCopyNo] = mateID; 663 } 752 } 664 << 753 665 //----- Read material densities and build ne 754 //----- Read material densities and build new materials if two voxels have 666 // same material but its density is in a di << 755 // same material but its density is in a different density interval 667 // (size of density intervals defined by den 756 // (size of density intervals defined by densityDiff) 668 G4double density; 757 G4double density; 669 // number of voxels from previously read sli 758 // number of voxels from previously read slices 670 voxelCopyNo = G4int((fZSliceHeaders.size() - << 759 voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels; 671 for (G4int ii = 0; ii < nVoxels; ++ii, voxel << 760 for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){ 672 fin >> density; 761 fin >> density; 673 << 762 674 //-- Get material from list of original ma 763 //-- Get material from list of original materials 675 mateID = unsigned(fMateIDs[voxelCopyNo]); << 764 mateID = fMateIDs[voxelCopyNo]; 676 G4Material* mateOrig = fOriginalMaterials[ << 765 G4Material* mateOrig = fOriginalMaterials[mateID]; 677 << 766 678 //-- Get density bin: middle point of the 767 //-- Get density bin: middle point of the bin in which the current 679 // density is included 768 // density is included 680 G4String newMateName = mateOrig->GetName() 769 G4String newMateName = mateOrig->GetName(); 681 G4float densityBin = 0.; << 770 float densityBin = 0.; 682 if (densityDiff != -1.) { << 771 if( densityDiff != -1.) { 683 densityBin = G4float(fDensityDiffs[mateI << 772 densityBin = fDensityDiffs[mateID] * >> 773 (G4int(density/fDensityDiffs[mateID])+0.5); 684 //-- Build the new material name 774 //-- Build the new material name 685 newMateName += G4UIcommand::ConvertToStr 775 newMateName += G4UIcommand::ConvertToString(densityBin); 686 } 776 } 687 << 777 688 //-- Look if a material with this name is 778 //-- Look if a material with this name is already created 689 // (because a previous voxel was already 779 // (because a previous voxel was already in this density bin) 690 unsigned int im; 780 unsigned int im; 691 for (im = 0; im < fMaterials.size(); ++im) << 781 for( im = 0; im < fMaterials.size(); im++ ){ 692 if (fMaterials[im]->GetName() == newMate << 782 if( fMaterials[im]->GetName() == newMateName ) { 693 break; 783 break; 694 } 784 } 695 } 785 } 696 //-- If material is already created use in 786 //-- If material is already created use index of this material 697 if (im != fMaterials.size()) { << 787 if( im != fMaterials.size() ) { 698 fMateIDs[voxelCopyNo] = im; 788 fMateIDs[voxelCopyNo] = im; 699 //-- else, create the material 789 //-- else, create the material 700 } << 790 } else { 701 else { << 791 if( densityDiff != -1.) { 702 if (densityDiff != -1.) { << 792 fMaterials.push_back( BuildMaterialWithChangingDensity( mateOrig, 703 fMaterials.push_back(BuildMaterialWith << 793 densityBin, newMateName ) ); 704 fMateIDs[voxelCopyNo] = fMaterials.siz << 794 fMateIDs[voxelCopyNo] = fMaterials.size()-1; 705 } << 795 } else { 706 else { << 796 G4cerr << " im " << im << " < " << fMaterials.size() << " name " 707 G4cerr << " im " << im << " < " << fMa << 797 << newMateName << G4endl; 708 G4Exception("DicomDetectorConstruction << 798 G4Exception("DicomDetectorConstruction::ReadPhantomDataFile", 709 "Wrong index in material") << 799 "", >> 800 FatalErrorInArgument, >> 801 "Wrong index in material"); //it should never reach here 710 } 802 } 711 } 803 } 712 } 804 } >> 805 713 } 806 } 714 807 715 //....oooOO0OOooo........oooOO0OOooo........oo 808 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 716 void DicomDetectorConstruction::MergeZSliceHea 809 void DicomDetectorConstruction::MergeZSliceHeaders() 717 { 810 { 718 //----- Images must have the same dimension 811 //----- Images must have the same dimension ... 719 fZSliceHeaderMerged = new DicomPhantomZSlice << 812 fZSliceHeaderMerged = new DicomPhantomZSliceHeader( *fZSliceHeaders[0] ); 720 for (unsigned int ii = 1; ii < fZSliceHeader << 813 for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ii++ ) { 721 *fZSliceHeaderMerged += *fZSliceHeaders[ii 814 *fZSliceHeaderMerged += *fZSliceHeaders[ii]; 722 } << 815 }; >> 816 723 } 817 } 724 818 725 //....oooOO0OOooo........oooOO0OOooo........oo 819 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 726 G4Material* DicomDetectorConstruction::BuildMa << 820 G4Material* DicomDetectorConstruction::BuildMaterialWithChangingDensity( 727 << 821 const G4Material* origMate, float density, G4String newMateName ) 728 << 729 { 822 { 730 //----- Copy original material, but with new 823 //----- Copy original material, but with new density 731 G4int nelem = G4int(origMate->GetNumberOfEle << 824 G4int nelem = origMate->GetNumberOfElements(); 732 G4Material* mate = << 825 G4Material* mate = new G4Material( newMateName, density*g/cm3, nelem, 733 new G4Material(newMateName, density * g / << 826 kStateUndefined, STP_Temperature ); 734 << 827 735 for (G4int ii = 0; ii < nelem; ++ii) { << 828 for( G4int ii = 0; ii < nelem; ii++ ){ 736 G4double frac = origMate->GetFractionVecto 829 G4double frac = origMate->GetFractionVector()[ii]; 737 G4Element* elem = const_cast<G4Element*>(o 830 G4Element* elem = const_cast<G4Element*>(origMate->GetElement(ii)); 738 mate->AddElement(elem, frac); << 831 mate->AddElement( elem, frac ); 739 } 832 } 740 << 833 741 return mate; 834 return mate; 742 } 835 } 743 836 744 //....oooOO0OOooo........oooOO0OOooo........oo 837 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 745 void DicomDetectorConstruction::ConstructPhant 838 void DicomDetectorConstruction::ConstructPhantomContainer() 746 { 839 { 747 //---- Extract number of voxels and voxel di 840 //---- Extract number of voxels and voxel dimensions 748 fNoVoxelsX = fZSliceHeaderMerged->GetNoVoxel << 841 fNVoxelX = fZSliceHeaderMerged->GetNoVoxelX(); 749 fNoVoxelsY = fZSliceHeaderMerged->GetNoVoxel << 842 fNVoxelY = fZSliceHeaderMerged->GetNoVoxelY(); 750 fNoVoxelsZ = fZSliceHeaderMerged->GetNoVoxel << 843 fNVoxelZ = fZSliceHeaderMerged->GetNoVoxelZ(); 751 << 844 752 fVoxelHalfDimX = fZSliceHeaderMerged->GetVox 845 fVoxelHalfDimX = fZSliceHeaderMerged->GetVoxelHalfX(); 753 fVoxelHalfDimY = fZSliceHeaderMerged->GetVox 846 fVoxelHalfDimY = fZSliceHeaderMerged->GetVoxelHalfY(); 754 fVoxelHalfDimZ = fZSliceHeaderMerged->GetVox 847 fVoxelHalfDimZ = fZSliceHeaderMerged->GetVoxelHalfZ(); 755 #ifdef G4VERBOSE 848 #ifdef G4VERBOSE 756 G4cout << " fNoVoxelsX " << fNoVoxelsX << " << 849 G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX 757 G4cout << " fNoVoxelsY " << fNoVoxelsY << " << 850 <<G4endl; 758 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " << 851 G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY 759 G4cout << " totalPixels " << fNoVoxelsX * fN << 852 <<G4endl; >> 853 G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ >> 854 <<G4endl; >> 855 G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl; 760 #endif 856 #endif 761 << 857 762 //----- Define the volume that contains all 858 //----- Define the volume that contains all the voxels 763 fContainer_solid = new G4Box("phantomContain << 859 fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX, 764 fNoVoxelsY * fV << 860 fNVoxelY*fVoxelHalfDimY, >> 861 fNVoxelZ*fVoxelHalfDimZ); 765 fContainer_logic = 862 fContainer_logic = 766 new G4LogicalVolume(fContainer_solid, << 863 new G4LogicalVolume( fContainer_solid, 767 // the material is not << 864 //the material is not important, it will be fully filled by the voxels 768 fMaterials[0], "phanto << 865 fMaterials[0], >> 866 "phantomContainer", >> 867 0, 0, 0 ); 769 //--- Place it on the world 868 //--- Place it on the world 770 G4double fOffsetX = (fZSliceHeaderMerged->Ge << 869 G4double fOffsetX = (fZSliceHeaderMerged->GetMaxX() + 771 G4double fOffsetY = (fZSliceHeaderMerged->Ge << 870 fZSliceHeaderMerged->GetMinX() ) /2.; 772 G4double fOffsetZ = (fZSliceHeaderMerged->Ge << 871 G4double fOffsetY = (fZSliceHeaderMerged->GetMaxY() + 773 G4ThreeVector posCentreVoxels(fOffsetX, fOff << 872 fZSliceHeaderMerged->GetMinY() ) /2.; >> 873 G4double fOffsetZ = (fZSliceHeaderMerged->GetMaxZ() + >> 874 fZSliceHeaderMerged->GetMinZ() ) /2.; >> 875 G4ThreeVector posCentreVoxels(fOffsetX,fOffsetY,fOffsetZ); 774 #ifdef G4VERBOSE 876 #ifdef G4VERBOSE 775 G4cout << " placing voxel container volume a 877 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl; 776 #endif 878 #endif 777 fContainer_phys = new G4PVPlacement(0, // r << 879 fContainer_phys = 778 posCentr << 880 new G4PVPlacement(0, // rotation 779 fContain << 881 posCentreVoxels, 780 "phantom << 882 fContainer_logic, // The logic volume 781 fWorld_l << 883 "phantomContainer", // Name 782 false, << 884 fWorld_logic, // Mother 783 1); // << 885 false, // No op. bool. >> 886 1); // Copy number >> 887 >> 888 //fContainer_logic->SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.))); 784 } 889 } 785 890 786 //....oooOO0OOooo........oooOO0OOooo........oo 891 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 787 void DicomDetectorConstruction::ConstructPhant 892 void DicomDetectorConstruction::ConstructPhantomContainerNew() 788 { 893 { 789 #ifdef G4_DCMTK 894 #ifdef G4_DCMTK 790 //---- Extract number of voxels and voxel di 895 //---- Extract number of voxels and voxel dimensions 791 # ifdef G4VERBOSE << 896 #ifdef G4VERBOSE 792 G4cout << " fNoVoxelsX " << fNoVoxelsX << " << 897 G4cout << " fNVoxelX " << fNVoxelX << " fVoxelHalfDimX " << fVoxelHalfDimX 793 G4cout << " fNoVoxelsY " << fNoVoxelsY << " << 898 <<G4endl; 794 G4cout << " fNoVoxelsZ " << fNoVoxelsZ << " << 899 G4cout << " fNVoxelY " << fNVoxelY << " fVoxelHalfDimY " << fVoxelHalfDimY 795 G4cout << " totalPixels " << fNoVoxelsX * fN << 900 <<G4endl; 796 # endif << 901 G4cout << " fNVoxelZ " << fNVoxelZ << " fVoxelHalfDimZ " << fVoxelHalfDimZ 797 << 902 <<G4endl; >> 903 G4cout << " totalPixels " << fNVoxelX*fNVoxelY*fNVoxelZ << G4endl; >> 904 #endif >> 905 798 //----- Define the volume that contains all 906 //----- Define the volume that contains all the voxels 799 fContainer_solid = new G4Box("phantomContain << 907 fContainer_solid = new G4Box("phantomContainer",fNVoxelX*fVoxelHalfDimX, 800 fNoVoxelsY * fV << 908 fNVoxelY*fVoxelHalfDimY, >> 909 fNVoxelZ*fVoxelHalfDimZ); 801 fContainer_logic = 910 fContainer_logic = 802 new G4LogicalVolume(fContainer_solid, << 911 new G4LogicalVolume( fContainer_solid, 803 // the material is not << 912 //the material is not important, it will be fully filled by the voxels 804 fMaterials[0], "phanto << 913 fMaterials[0], >> 914 "phantomContainer", >> 915 0, 0, 0 ); 805 916 806 G4ThreeVector posCentreVoxels((fMinX + fMaxX << 917 G4ThreeVector posCentreVoxels((fMinX+fMaxX)/2.,(fMinY+fMaxY)/2.,(fMinZ+fMaxZ)/2.); 807 # ifdef G4VERBOSE << 918 #ifdef G4VERBOSE 808 G4cout << " placing voxel container volume a 919 G4cout << " placing voxel container volume at " << posCentreVoxels << G4endl; 809 # endif << 920 #endif 810 fContainer_phys = new G4PVPlacement(0, // r << 921 fContainer_phys = 811 posCentr << 922 new G4PVPlacement(0, // rotation 812 fContain << 923 posCentreVoxels, 813 "phantom << 924 fContainer_logic, // The logic volume 814 fWorld_l << 925 "phantomContainer", // Name 815 false, << 926 fWorld_logic, // Mother 816 1); // << 927 false, // No op. bool. >> 928 1); // Copy number >> 929 >> 930 //fContainer_logic->SetVisAttributes(new G4VisAttributes(G4Colour(1.,0.,0.))); 817 #endif 931 #endif 818 } 932 } 819 933 >> 934 #include "G4SDManager.hh" 820 #include "G4MultiFunctionalDetector.hh" 935 #include "G4MultiFunctionalDetector.hh" 821 #include "G4PSDoseDeposit.hh" 936 #include "G4PSDoseDeposit.hh" 822 #include "G4PSDoseDeposit3D.hh" 937 #include "G4PSDoseDeposit3D.hh" 823 #include "G4SDManager.hh" << 824 938 825 //....oooOO0OOooo........oooOO0OOooo........oo 939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo. 826 void DicomDetectorConstruction::SetScorer(G4Lo 940 void DicomDetectorConstruction::SetScorer(G4LogicalVolume* voxel_logic) 827 { 941 { >> 942 828 #ifdef G4VERBOSE 943 #ifdef G4VERBOSE 829 G4cout << "\t SET SCORER : " << voxel_logic- 944 G4cout << "\t SET SCORER : " << voxel_logic->GetName() << G4endl; 830 #endif 945 #endif 831 << 946 832 fScorers.insert(voxel_logic); 947 fScorers.insert(voxel_logic); >> 948 833 } 949 } 834 950 835 //....oooOO0OOooo........oooOO0OOooo........oo 951 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 836 952 837 void DicomDetectorConstruction::ConstructSDand 953 void DicomDetectorConstruction::ConstructSDandField() 838 { 954 { >> 955 839 #ifdef G4VERBOSE 956 #ifdef G4VERBOSE 840 G4cout << "\t CONSTRUCT SD AND FIELD" << G4e 957 G4cout << "\t CONSTRUCT SD AND FIELD" << G4endl; 841 #endif 958 #endif 842 << 959 843 // G4SDManager* SDman = G4SDManager::GetSDMp << 960 //G4SDManager* SDman = G4SDManager::GetSDMpointer(); 844 << 961 845 // SDman->SetVerboseLevel(1); << 962 //SDman->SetVerboseLevel(1); 846 << 963 847 // 964 // 848 // Sensitive Detector Name 965 // Sensitive Detector Name 849 G4String concreteSDname = "phantomSD"; 966 G4String concreteSDname = "phantomSD"; 850 std::vector<G4String> scorer_names; 967 std::vector<G4String> scorer_names; 851 scorer_names.push_back(concreteSDname); 968 scorer_names.push_back(concreteSDname); 852 //------------------------ 969 //------------------------ 853 // MultiFunctionalDetector 970 // MultiFunctionalDetector 854 //------------------------ 971 //------------------------ 855 // 972 // 856 // Define MultiFunctionalDetector with name. 973 // Define MultiFunctionalDetector with name. 857 // declare MFDet as a MultiFunctionalDetecto 974 // declare MFDet as a MultiFunctionalDetector scorer 858 G4MultiFunctionalDetector* MFDet = new G4Mul << 975 G4MultiFunctionalDetector* MFDet = 859 G4SDManager::GetSDMpointer()->AddNewDetector << 976 new G4MultiFunctionalDetector(concreteSDname); 860 char* nest = std::getenv("DICOM_NESTED_PARAM << 977 G4SDManager::GetSDMpointer()->AddNewDetector( MFDet ); 861 if (nest && G4String(nest) == "1") { << 978 //G4VPrimitiveScorer* dosedep = new G4PSDoseDeposit("DoseDeposit"); 862 G4VPrimitiveScorer* dosedep = new G4PSDose << 979 G4VPrimitiveScorer* dosedep = 863 "DoseDeposit", fNoVoxelsZ, fNoVoxelsY, f << 980 new G4PSDoseDeposit3D("DoseDeposit", fNVoxelX, fNVoxelY, fNVoxelZ); 864 // - the last 3 arguments correspond to th << 981 MFDet->RegisterPrimitive(dosedep); 865 MFDet->RegisterPrimitive(dosedep); << 982 866 } << 983 for(std::set<G4LogicalVolume*>::iterator ite = fScorers.begin(); 867 else { << 984 ite != fScorers.end(); ++ite) { 868 G4VPrimitiveScorer* dosedep = new G4PSDose << 869 MFDet->RegisterPrimitive(dosedep); << 870 } << 871 << 872 for (auto ite = fScorers.cbegin(); ite != fS << 873 SetSensitiveDetector(*ite, MFDet); 985 SetSensitiveDetector(*ite, MFDet); 874 } 986 } >> 987 >> 988 /*if(DicomRunAction::Instance()->GetDicomRun()) { >> 989 DicomRunAction::Instance()->GetDicomRun()->ConstructMFD(scorer_names); >> 990 }*/ >> 991 875 } 992 } 876 993