Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc

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

Diff markup

Differences between /examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc (Version 11.3.0) and /examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc (Version 10.4.p3)


  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