| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // << 26 #include "globals.hh"
27 /// \file medical/DICOM/src/DicomPartialDetect <<
28 /// \brief Implementation of the DicomPartialD <<
29 // <<
30 27
31 #include "DicomPartialDetectorConstruction.hh" 28 #include "DicomPartialDetectorConstruction.hh"
32 29
33 #include "G4Box.hh" 30 #include "G4Box.hh"
34 #include "G4Colour.hh" <<
35 #include "G4Element.hh" <<
36 #include "G4LogicalVolume.hh" 31 #include "G4LogicalVolume.hh"
37 #include "G4Material.hh" <<
38 #include "G4NistManager.hh" <<
39 #include "G4PVParameterised.hh" <<
40 #include "G4PVPlacement.hh" <<
41 #include "G4PartialPhantomParameterisation.hh" <<
42 #include "G4SystemOfUnits.hh" <<
43 #include "G4Tubs.hh" <<
44 #include "G4UIcommand.hh" <<
45 #include "G4VPhysicalVolume.hh" 32 #include "G4VPhysicalVolume.hh"
>> 33 #include "G4PVPlacement.hh"
>> 34 #include "G4PVParameterised.hh"
>> 35 #include "G4Material.hh"
>> 36 #include "G4Element.hh"
46 #include "G4VisAttributes.hh" 37 #include "G4VisAttributes.hh"
>> 38 #include "G4Colour.hh"
47 #include "G4ios.hh" 39 #include "G4ios.hh"
48 #include "G4tgbMaterialMgr.hh" << 40 #include "G4PartialPhantomParameterisation.hh"
>> 41 #include "G4NistManager.hh"
>> 42 #include "G4UIcommand.hh"
>> 43 #include "G4Tubs.hh"
>> 44
49 #include "G4tgbVolumeMgr.hh" 45 #include "G4tgbVolumeMgr.hh"
>> 46 #include "G4tgbMaterialMgr.hh"
>> 47 #include "G4tgrMessenger.hh"
50 #include "G4tgrMessenger.hh" 48 #include "G4tgrMessenger.hh"
51 #include "globals.hh" <<
52 49
53 //....oooOO0OOooo........oooOO0OOooo........oo << 50 //---------------------------------------------------------------------------
54 DicomPartialDetectorConstruction::DicomPartial 51 DicomPartialDetectorConstruction::DicomPartialDetectorConstruction()
55 : DicomDetectorConstruction(), << 52 {
56 fPartialPhantomParam(0), << 53 }
57 fNoVoxels(0), << 54
58 fDimX(0), << 55 //---------------------------------------------------------------------------
59 fDimY(0), << 56 DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction()
60 fDimZ(0), << 57 {
61 fOffsetX(0), << 58 }
62 fOffsetY(0), <<
63 fOffsetZ(0) <<
64 {} <<
65 59
66 //....oooOO0OOooo........oooOO0OOooo........oo <<
67 DicomPartialDetectorConstruction::~DicomPartia <<
68 60
69 //....oooOO0OOooo........oooOO0OOooo........oo << 61 //-------------------------------------------------------------
70 G4VPhysicalVolume* DicomPartialDetectorConstru 62 G4VPhysicalVolume* DicomPartialDetectorConstruction::Construct()
71 { 63 {
72 InitialisationOfMaterials(); 64 InitialisationOfMaterials();
73 65
74 //----- Build world 66 //----- Build world
75 G4double worldXDimension = 1. * m; << 67 G4double worldXDimension = 1.*m;
76 G4double worldYDimension = 1. * m; << 68 G4double worldYDimension = 1.*m;
77 G4double worldZDimension = 1. * m; << 69 G4double worldZDimension = 1.*m;
78 << 70
79 G4Box* world_box = new G4Box("WorldSolid", w << 71 G4Box* world_solid = new G4Box( "WorldSolid",
80 << 72 worldXDimension,
81 fWorld_logic = new G4LogicalVolume(world_box << 73 worldYDimension,
82 << 74 worldZDimension );
83 G4VPhysicalVolume* world_phys = << 75
84 new G4PVPlacement(0, G4ThreeVector(0, 0, 0 << 76 world_logic = new G4LogicalVolume( world_solid,
>> 77 air,
>> 78 "WorldLogical",
>> 79 0, 0, 0 );
>> 80
>> 81 G4VPhysicalVolume* world_phys = new G4PVPlacement( 0,
>> 82 G4ThreeVector(0,0,0),
>> 83 "World",
>> 84 world_logic,
>> 85 0,
>> 86 false,
>> 87 0 );
85 88
86 ReadPhantomData(); 89 ReadPhantomData();
87 90
88 ConstructPhantomContainer(); 91 ConstructPhantomContainer();
89 ConstructPhantom(); 92 ConstructPhantom();
90 93
91 return world_phys; 94 return world_phys;
92 } 95 }
93 96
94 //....oooOO0OOooo........oooOO0OOooo........oo << 97
>> 98 //-------------------------------------------------------------
95 void DicomPartialDetectorConstruction::ReadPha 99 void DicomPartialDetectorConstruction::ReadPhantomData()
96 { 100 {
97 G4String dataFile = "Data.dat"; << 101 std::ifstream finDF("Data.dat");
98 std::ifstream finDF(dataFile.c_str()); <<
99 G4String fname; 102 G4String fname;
100 if (finDF.good() != 1) { << 103 if(finDF.good() != 1 ) {
101 G4Exception(" DicomPartialDetectorConstruc << 104 G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData. Problem reading data file: Data.dat");
102 " Problem reading data file: D <<
103 } 105 }
104 106
105 G4int compression; 107 G4int compression;
106 finDF >> compression; // not used here << 108 finDF >> compression; // not used here
107 109
108 G4int numFiles; << 110 G4int fNoFiles;
109 finDF >> numFiles; // only 1 file supported << 111 finDF >> fNoFiles; // only 1 file supported for the moment
110 if (numFiles != 1) { << 112 if( fNoFiles != 1 ) {
111 G4Exception("DicomPartialDetectorConstruct << 113 G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
112 "More than 1 DICOM file is not << 114 "More than 1 DICOM file is not supported",
>> 115 FatalErrorInArgument,
>> 116 "");
113 } 117 }
114 for (G4int i = 0; i < numFiles; i++) { << 118 for(G4int i = 0; i < fNoFiles; i++ ) {
115 finDF >> fname; 119 finDF >> fname;
116 //--- Read one data file 120 //--- Read one data file
117 fname += ".g4pdcm"; 121 fname += ".g4pdcm";
118 ReadPhantomDataFile(fname); 122 ReadPhantomDataFile(fname);
119 } 123 }
120 124
121 finDF.close(); 125 finDF.close();
>> 126
122 } 127 }
123 128
124 //....oooOO0OOooo........oooOO0OOooo........oo << 129
>> 130 //---------------------------------------------------------------------------
125 void DicomPartialDetectorConstruction::ReadPha 131 void DicomPartialDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
126 { 132 {
127 // G4String filename = "phantom.g4pdcm"; 133 // G4String filename = "phantom.g4pdcm";
128 #ifdef G4VERBOSE 134 #ifdef G4VERBOSE
129 G4cout << " DicomDetectorConstruction::ReadP 135 G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
130 #endif << 136 #endif
131 std::ifstream fin(fname.c_str(), std::ios_ba 137 std::ifstream fin(fname.c_str(), std::ios_base::in);
132 if (!fin.is_open()) { << 138 if( !fin.is_open() ) {
133 G4Exception("DicomPartialDetectorConstruct << 139 G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile. File not found " + fname );
134 G4String("File not found " + f <<
135 } 140 }
136 G4int nMaterials; << 141 G4int nMaterials;
137 fin >> nMaterials; 142 fin >> nMaterials;
138 G4String stemp; 143 G4String stemp;
139 G4int nmate; 144 G4int nmate;
140 for (G4int ii = 0; ii < nMaterials; ii++) { << 145 for( G4int ii = 0; ii < nMaterials; ii++ ){
141 fin >> nmate >> stemp; 146 fin >> nmate >> stemp;
142 G4cout << "DicomPartialDetectorConstructio << 147 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate " << ii << " = " << nmate << " mate " << stemp << G4endl;
143 << nmate << " mate " << stemp << G4 << 148 if( ii != nmate ) G4Exception("DicomPartialDetectorConstruction::ReadPhantomData material number should be in increasing order: wrong material number ");
144 if (ii != nmate) << 149 }
145 G4Exception("DicomPartialDetectorConstru << 150
146 "Material number should be i << 151 fin >> nVoxelX >> nVoxelY >> nVoxelZ;
147 wrong material numbe << 152 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z " << nVoxelX << " " << nVoxelY << " " << nVoxelZ << G4endl;
148 } << 153 fin >> offsetX >> dimX;
149 << 154 dimX = (dimX - offsetX)/nVoxelX;
150 fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxels << 155 fin >> offsetY >> dimY;
151 G4cout << "DicomPartialDetectorConstruction: << 156 dimY = (dimY - offsetY)/nVoxelY;
152 << fNoVoxelsY << " " << fNoVoxelsZ << << 157 fin >> offsetZ >> dimZ;
153 fin >> fOffsetX >> fDimX; << 158 dimZ = (dimZ - offsetZ)/nVoxelZ;
154 fDimX = (fDimX - fOffsetX) / fNoVoxelsX; << 159 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimX " << dimX << " offsetX " << offsetX << G4endl;
155 fin >> fOffsetY >> fDimY; << 160 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY " << dimY << " offsetY " << offsetY << G4endl;
156 fDimY = (fDimY - fOffsetY) / fNoVoxelsY; << 161 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ " << dimZ << " offsetZ " << offsetZ << G4endl;
157 fin >> fOffsetZ >> fDimZ; <<
158 fDimZ = (fDimZ - fOffsetZ) / fNoVoxelsZ; <<
159 G4cout << "DicomPartialDetectorConstruction: <<
160 << fOffsetX << G4endl; <<
161 G4cout << "DicomPartialDetectorConstruction: <<
162 << fOffsetY << G4endl; <<
163 G4cout << "DicomPartialDetectorConstruction: <<
164 << fOffsetZ << G4endl; <<
165 162
166 //--- Read voxels that are filled 163 //--- Read voxels that are filled
167 fNoVoxels = 0; << 164 theNVoxels = 0;
168 // G4bool* isFilled = new G4bool[fNoVoxelsX << 165 // G4bool* isFilled = new G4bool[nVoxelX*nVoxelY*nVoxelZ];
169 // fFilledIDs = new size_t[fNoVoxelsZ*fNoVo << 166 // theFilledIDs = new size_t[nVoxelZ*nVoxelY+1];
170 //? fFilledIDs.insert(0); << 167 //? theFilledIDs.insert(0);
171 G4int ifxmin1, ifxmax1; << 168 int ifxmin1, ifxmax1;
172 for (G4int iz = 0; iz < fNoVoxelsZ; iz++) { << 169 for( G4int iz = 0; iz < nVoxelZ; iz++ ) {
173 std::map<G4int, G4int> ifmin, ifmax; << 170 std::map< G4int, G4int > ifmin, ifmax;
174 for (G4int iy = 0; iy < fNoVoxelsY; iy++) << 171 for( G4int iy = 0; iy < nVoxelY; iy++ ) {
175 fin >> ifxmin1 >> ifxmax1; << 172 fin >> ifxmin1 >> ifxmax1;
176 // check coherence ... 173 // check coherence ...
177 174
178 ifmin[iy] = ifxmin1; 175 ifmin[iy] = ifxmin1;
179 ifmax[iy] = ifxmax1; 176 ifmax[iy] = ifxmax1;
180 G4int ifxdiff = ifxmax1 - ifxmin1 + 1; << 177 G4int ifxdiff = ifxmax1-ifxmin1+1;
181 if (ifxmax1 == -1 && ifxmin1 == -1) ifxd << 178 if( ifxmax1 == -1 && ifxmin1 == -1 ) ifxdiff = 0;
182 fFilledIDs.insert(std::pair<G4int, G4int << 179 // theFilledIDs.insert(std::multimap<G4int,size_t>::value_type(ifxdiff+theNVoxels, ifxmin1));
183 G4cout << "DicomPartialDetectorConstruct << 180 theFilledIDs.insert(std::pair<size_t,size_t>(ifxdiff+theNVoxels-1, ifxmin1));
184 << " FilledIDs " << ifxdiff + fNo << 181 G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert FilledIDs " << ifxdiff+theNVoxels-1 << " min " << ifxmin1 << " N= " << theFilledIDs.size() << G4endl;
185 << " N= " << fFilledIDs.size() << << 182 //filledIDs[iz*nVoxelY+iy+1] = ifxmax1-ifxmin1+1;
186 // filledIDs[iz*fNoVoxelsY+iy+1] = ifxma << 183 for( G4int ix = 0; ix < nVoxelX; ix++ ) {
187 for (G4int ix = 0; ix < fNoVoxelsX; ix++ << 184 // G4cout << "REGNAV DicomPartialDetectorConstruction::ReadPhantomData checking to add voxel iz " << iz << " iy " << iy << " ix " << ix << " min= " << ifxmin1 << " max= " << ifxmax1 << G4endl;
188 if (ix >= G4int(ifxmin1) && ix <= G4in << 185 // G4int copyNo = ix + (iy)*nVoxelX + (iz)*nVoxelX*nVoxelY;
189 fNoVoxels++; << 186 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
190 } << 187 // isFilled[copyNo] = true;
191 else { << 188 theNVoxels++;
192 } << 189 } else {
>> 190 // isFilled[copyNo] = false;
>> 191 }
193 } 192 }
194 } 193 }
195 fFilledMins[iz] = ifmin; << 194 theFilledMins[iz] = ifmin;
196 fFilledMaxs[iz] = ifmax; << 195 theFilledMaxs[iz] = ifmax;
197 } 196 }
198 197
199 //--- Read material IDs 198 //--- Read material IDs
200 G4int mateID1; 199 G4int mateID1;
201 fMateIDs = new size_t[fNoVoxelsX * fNoVoxels << 200 mateIDs = new size_t[nVoxelX*nVoxelY*nVoxelZ];
202 G4int copyNo = 0; 201 G4int copyNo = 0;
203 for (G4int iz = 0; iz < fNoVoxelsZ; iz++) { << 202 for( G4int iz = 0; iz < nVoxelZ; iz++ ) {
204 std::map<G4int, G4int> ifmin = fFilledMins << 203 std::map< G4int, G4int > ifmin = theFilledMins[iz];
205 std::map<G4int, G4int> ifmax = fFilledMaxs << 204 std::map< G4int, G4int > ifmax = theFilledMaxs[iz];
206 for (G4int iy = 0; iy < fNoVoxelsY; iy++) << 205 for( G4int iy = 0; iy < nVoxelY; iy++ ) {
207 ifxmin1 = ifmin[iy]; 206 ifxmin1 = ifmin[iy];
208 ifxmax1 = ifmax[iy]; 207 ifxmax1 = ifmax[iy];
209 for (G4int ix = 0; ix < fNoVoxelsX; ix++ << 208 for( G4int ix = 0; ix < nVoxelX; ix++ ) {
210 if (ix >= G4int(ifxmin1) && ix <= G4in << 209 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
211 fin >> mateID1; << 210 fin >> mateID1;
212 if (mateID1 < 0 || mateID1 >= nMater << 211 if( mateID1 < 0 || mateID1 >= nMaterials ) {
213 G4Exception("DicomPartialDetectorC << 212 G4Exception("DicomPartialDetectorConstruction::ReadPhantomData",
214 G4String("Wrong index << 213 "Wrong index in phantom file",
215 + G4UIcommand << 214 FatalException,
216 + G4UIcommand << 215 G4String("It should be between 0 and "
217 .c_str()); << 216 + G4UIcommand::ConvertToString(nMaterials-1)
218 } << 217 + ", while it is "
219 fMateIDs[copyNo] = mateID1; << 218 + G4UIcommand::ConvertToString(mateID1)).c_str());
220 copyNo++; << 219 }
221 } << 220 mateIDs[copyNo] = mateID1;
>> 221 copyNo++;
>> 222 }
222 } 223 }
223 } 224 }
224 } 225 }
225 226
226 ReadVoxelDensitiesPartial(fin); << 227 ReadVoxelDensitiesPartial( fin );
227 228
228 fin.close(); 229 fin.close();
229 } 230 }
230 231
231 //....oooOO0OOooo........oooOO0OOooo........oo << 232
232 void DicomPartialDetectorConstruction::ReadVox << 233 //------------------------------------------------------------------------
>> 234 void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial( std::ifstream& fin )
233 { 235 {
234 std::map<G4int, std::pair<G4double, G4double << 236 std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
235 std::map<G4int, std::pair<G4double, G4double << 237 std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
236 for (G4int ii = 0; ii < G4int(fOriginalMater << 238 for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
237 densiMinMax[ii] = std::pair<G4double, G4do << 239 densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
238 } 240 }
239 241
240 //----- Define density differences (maximum << 242 //----- Define density differences (maximum density difference to create a new material)
241 // a new material) << 243 char* part = getenv( "DICOM_CHANGE_MATERIAL_DENSITY" );
242 char* part = std::getenv("DICOM_CHANGE_MATER <<
243 G4double densityDiff = -1.; 244 G4double densityDiff = -1.;
244 if (part) densityDiff = G4UIcommand::Convert << 245 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
245 std::map<G4int, G4double> densityDiffs; << 246
246 if (densityDiff != -1.) { << 247 std::map<G4int,G4double> densityDiffs;
247 for (G4int ii = 0; ii < G4int(fOriginalMat << 248 if( densityDiff != -1. ) {
248 densityDiffs[ii] = densityDiff; // curr << 249 for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
249 // with same step << 250 densityDiffs[ii] = densityDiff; //currently all materials with same step
250 } 251 }
251 } << 252 }else {
252 else { << 253 if( fMaterials.size() == 0 ) { // do it only for first slice
253 if (fMaterials.size() == 0) { // do it on << 254 for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
254 for (unsigned int ii = 0; ii < fOriginal << 255 fMaterials.push_back( fOriginalMaterials[ii] );
255 fMaterials.push_back(fOriginalMaterial <<
256 } 256 }
257 } 257 }
258 } 258 }
259 // densityDiffs[0] = 0.0001; //fAir <<
260 259
261 //--- Calculate the average material density << 260 // densityDiffs[0] = 0.0001; //air
262 std::map<std::pair<G4Material*, G4int>, matI <<
263 261
>> 262 //--- Calculate the average material density for each material/density bin
>> 263 std::map< std::pair<G4Material*,G4int>, matInfo* > newMateDens;
>> 264
264 G4double dens1; 265 G4double dens1;
265 G4int ifxmin1, ifxmax1; << 266 size_t ifxmin1, ifxmax1;
266 267
267 //---- Read the material densities 268 //---- Read the material densities
268 G4int copyNo = 0; 269 G4int copyNo = 0;
269 for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) { << 270 for( G4int iz = 0; iz < nVoxelZ; iz++ ) {
270 std::map<G4int, G4int> ifmin = fFilledMins << 271 std::map< G4int, G4int > ifmin = theFilledMins[iz];
271 std::map<G4int, G4int> ifmax = fFilledMaxs << 272 std::map< G4int, G4int > ifmax = theFilledMaxs[iz];
272 for (G4int iy = 0; iy < fNoVoxelsY; ++iy) << 273 for( G4int iy = 0; iy < nVoxelY; iy++ ) {
273 ifxmin1 = ifmin[iy]; 274 ifxmin1 = ifmin[iy];
274 ifxmax1 = ifmax[iy]; 275 ifxmax1 = ifmax[iy];
275 for (G4int ix = 0; ix < fNoVoxelsX; ++ix << 276 for( G4int ix = 0; ix < nVoxelX; ix++ ) {
276 if (ix >= G4int(ifxmin1) && ix <= G4in << 277 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
277 fin >> dens1; << 278 fin >> dens1;
278 //--- store the minimum and maximum << 279 // G4cout << ix << " " << iy << " " << iz << " filling mateIDs " << copyNo << " = " << atoi(stemp.c_str())-1 << " " << stemp << G4endl;
279 //(just for printing) << 280
280 mpite = densiMinMax.find(G4int(fMate << 281 //--- store the minimum and maximum density for each material (just for printing)
281 if (dens1 < (*mpite).second.first) ( << 282 mpite = densiMinMax.find( mateIDs[copyNo] );
282 if (dens1 > (*mpite).second.second) << 283 if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
283 << 284 if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
284 //--- Get material from original lis << 285
285 G4int mateID = G4int(fMateIDs[copyNo << 286 //--- Get material from original list of material in file
286 G4Material* mate = fOriginalMaterial << 287 int mateID = mateIDs[copyNo];
287 // G4cout << copyNo << " mate << 288 G4Material* mate = fOriginalMaterials[mateID];
288 //--- Check if density is equal to t << 289 // G4cout << copyNo << " mateID " << mateID << G4endl;
289 if (std::fabs(dens1 - mate->GetDensi << 290 //--- Check if density is equal to the original material density
290 << 291 if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 ) continue;
291 //--- Build material name with fOrig << 292
292 // float densityBin = density << 293 //--- Build material name with fOriginalMaterials name + density
293 // (G4int(dens1/densityDiffs[ << 294 // float densityBin = densityDiffs[mateID] * (G4int(dens1/densityDiffs[mateID])+0.5);
294 G4String newMateName = mate->GetName << 295 G4String newMateName = mate->GetName();
295 << 296
296 G4int densityBin = 0; << 297 G4int densityBin = 0;
297 // G4cout << " densityBin " < << 298 // G4cout << " densityBin " << densityBin << " " << dens1 << " " <<densityDiffs[mateID] << G4endl;
298 // << dens1 << " " <<densityD << 299
299 << 300 if( densityDiff != -1.) {
300 if (densityDiff != -1.) { << 301 densityBin = (G4int(dens1/densityDiffs[mateID]));
301 densityBin = (G4int(dens1 / densit << 302 newMateName = mate->GetName()+G4UIcommand::ConvertToString(densityBin);
302 newMateName = mate->GetName() + G4 << 303 }
303 } << 304
304 << 305 //--- Look if it is the first voxel with this material/densityBin
305 //--- Look if it is the first voxel << 306 std::pair<G4Material*,G4int> matdens(mate, densityBin );
306 std::pair<G4Material*, G4int> matden << 307
307 << 308 std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite = newMateDens.find( matdens );
308 auto mppite = newMateDens.find(matde << 309 if( mppite != newMateDens.end() ){
309 if (mppite != newMateDens.cend()) { << 310 matInfo* mi = (*mppite).second;
310 matInfo* mi = (*mppite).second; << 311 mi->sumdens += dens1;
311 mi->fSumdens += dens1; << 312 mi->nvoxels++;
312 mi->fNvoxels++; << 313 mateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
313 fMateIDs[copyNo] = fOriginalMateri << 314 // G4cout << copyNo << " mat new again " << fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
314 // G4cout << copyNo << " mat new << 315 } else {
315 //<< fOriginalMaterials.size()-1 + << 316 matInfo* mi = new matInfo;
316 } << 317 mi->sumdens = dens1;
317 else { << 318 mi->nvoxels = 1;
318 matInfo* mi = new matInfo; << 319 mi->id = newMateDens.size()+1;
319 mi->fSumdens = dens1; << 320 newMateDens[matdens] = mi;
320 mi->fNvoxels = 1; << 321 mateIDs[copyNo] = fOriginalMaterials.size()-1 + mi->id;
321 mi->fId = G4int(newMateDens.size() << 322 // G4cout << copyNo << " mat new first " << fOriginalMaterials.size()-1 + mi->id << G4endl;
322 newMateDens[matdens] = mi; << 323 }
323 fMateIDs[copyNo] = fOriginalMateri << 324 copyNo++;
324 // G4cout << copyNo << " << 325 // G4cout << ix << " " << iy << " " << iz << " filling mateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
325 // << fOriginalMaterials. << 326 // mateIDs[copyNo] = atoi(cid)-1;
326 } << 327 }
327 copyNo++; <<
328 // G4cout << ix << " " << iy <<
329 // << " filling fMateIDs " << copyN <<
330 // fMateIDs[copyNo] = atoi(ci <<
331 } <<
332 } 328 }
333 } 329 }
334 } 330 }
335 331
336 //----- Build the list of phantom materials 332 //----- Build the list of phantom materials that go to Parameterisation
337 //--- Add original materials 333 //--- Add original materials
338 for (auto mimite = fOriginalMaterials.cbegin << 334 std::vector<G4Material*>::const_iterator mimite;
339 fPhantomMaterials.push_back((*mimite)); << 335 for( mimite = fOriginalMaterials.begin(); mimite != fOriginalMaterials.end(); mimite++ ){
>> 336 thePhantomMaterials.push_back( (*mimite) );
340 } 337 }
341 // << 338 //
342 //---- Build and add new materials 339 //---- Build and add new materials
343 for (auto mppite = newMateDens.cbegin(); mpp << 340 std::map< std::pair<G4Material*,G4int>, matInfo* >::iterator mppite;
344 G4double averdens = (*mppite).second->fSum << 341 for( mppite= newMateDens.begin(); mppite != newMateDens.end(); mppite++ ){
345 G4double saverdens = G4int(1000.001 * aver << 342 G4double averdens = (*mppite).second->sumdens/(*mppite).second->nvoxels;
346 G4cout << "GmReadPhantomGeometry::ReadVoxe << 343 G4double saverdens = G4int(1000.001*averdens)/1000.;
347 << saverdens << " -> " << G4int(100 << 344 G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> " << saverdens << " -> " << G4int(1000*averdens) << " " << G4int(1000*averdens)/1000 << " " << G4int(1000*averdens)/1000. << G4endl;
348 << " " << G4int(1000 * averdens) / << 345
349 << 346 G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> " << saverdens << " -> " << G4UIcommand::ConvertToString(saverdens) << " Nvoxels " <<(*mppite).second->nvoxels << G4endl;
350 G4cout << "GmReadPhantomGeometry::ReadVoxe << 347 G4String mateName = ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
351 << saverdens << " -> " << G4UIcomma << 348 thePhantomMaterials.push_back( BuildMaterialWithChangingDensity( (*mppite).first.first, averdens, mateName ) );
352 << (*mppite).second->fNvoxels << G4 <<
353 G4String mateName = <<
354 ((*mppite).first).first->GetName() + "_" <<
355 fPhantomMaterials.push_back( <<
356 BuildMaterialWithChangingDensity((*mppit <<
357 } 349 }
>> 350
358 } 351 }
359 352
360 //....oooOO0OOooo........oooOO0OOooo........oo << 353
>> 354 //-----------------------------------------------------------------------
361 void DicomPartialDetectorConstruction::Constru 355 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
362 { 356 {
363 // build a clinder that encompass the partia << 357 //build a clinder that encompass the partial phantom built in the example
364 // /dicom/intersectWithUserVolume 0. 0. 0. << 358 // /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
365 //---- Extract number of voxels and voxel di 359 //---- Extract number of voxels and voxel dimensions
366 360
367 G4String contname = "PhantomContainer"; << 361 G4String contname= "PhantomContainer";
368 362
369 //----- Define the volume that contains all 363 //----- Define the volume that contains all the voxels
370 G4Tubs* container_tube = new G4Tubs(contname << 364 G4Tubs* container_solid = new G4Tubs(contname,0., 200., 100., 0., 360*deg);
371 365
372 fContainer_logic = new G4LogicalVolume(conta << 366 container_logic =
>> 367 new G4LogicalVolume( container_solid,
>> 368 air,
>> 369 contname,
>> 370 0, 0, 0 );
373 371
374 G4ThreeVector posCentreVoxels(0., 0., 0.); << 372 G4ThreeVector posCentreVoxels(0.,0.,0.);
375 // G4cout << " PhantomContainer posCentreVox << 373 //G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
376 G4RotationMatrix* rotm = new G4RotationMatri 374 G4RotationMatrix* rotm = new G4RotationMatrix;
377 375
378 fContainer_phys = new G4PVPlacement(rotm, / << 376 container_phys =
379 posCentr << 377 new G4PVPlacement(rotm, // rotation
380 fContain << 378 posCentreVoxels,
381 contname << 379 container_logic, // The logic volume
382 fWorld_l << 380 contname, // Name
383 false, << 381 world_logic, // Mother
384 1); // << 382 false, // No op. bool.
>> 383 1); // Copy number
>> 384
385 } 385 }
386 386
387 //....oooOO0OOooo........oooOO0OOooo........oo << 387
>> 388 //---------------------------------------------------------------------------
388 void DicomPartialDetectorConstruction::Constru 389 void DicomPartialDetectorConstruction::ConstructPhantom()
389 { 390 {
390 G4String OptimAxis = "kXAxis"; 391 G4String OptimAxis = "kXAxis";
391 G4bool bSkipEqualMaterials = 0; 392 G4bool bSkipEqualMaterials = 0;
392 G4int regStructureID = 2; 393 G4int regStructureID = 2;
393 394
394 G4ThreeVector posCentreVoxels(0., 0., 0.); << 395 G4ThreeVector posCentreVoxels(0.,0.,0.);
395 396
396 //----- Build phantom 397 //----- Build phantom
397 G4String voxelName = "phantom"; 398 G4String voxelName = "phantom";
398 G4Box* phantom_solid = new G4Box(voxelName, << 399 G4Box* phantom_solid = new G4Box(voxelName, dimX/2., dimY/2., dimZ/2. );
399 G4LogicalVolume* phantom_logic = new G4Logic << 400 G4LogicalVolume* phantom_logic =
>> 401 new G4LogicalVolume( phantom_solid,
>> 402 air,
>> 403 voxelName,
>> 404 0, 0, 0 );
400 G4bool pVis = 0; 405 G4bool pVis = 0;
401 if (!pVis) { << 406 if( !pVis ) {
402 G4VisAttributes* visAtt = new G4VisAttribu << 407 G4VisAttributes* visAtt = new G4VisAttributes( FALSE );
403 phantom_logic->SetVisAttributes(visAtt); << 408 phantom_logic->SetVisAttributes( visAtt );
404 } 409 }
405 410
406 G4double theSmartless = 0.5; 411 G4double theSmartless = 0.5;
407 // fContainer_logic->SetSmartless( theSmart << 412 // container_logic->SetSmartless( theSmartless );
408 phantom_logic->SetSmartless(theSmartless); << 413 phantom_logic->SetSmartless( theSmartless );
409 414
410 fPartialPhantomParam = new G4PartialPhantomP << 415 thePartialPhantomParam = new G4PartialPhantomParameterisation();
411 416
412 fPartialPhantomParam->SetMaterials(fPhantomM << 417 thePartialPhantomParam->SetMaterials( thePhantomMaterials );
413 fPartialPhantomParam->SetVoxelDimensions(fDi << 418 thePartialPhantomParam->SetVoxelDimensions( dimX/2., dimY/2., dimZ/2. );
414 fPartialPhantomParam->SetNoVoxels(fNoVoxelsX << 419 thePartialPhantomParam->SetNoVoxel( nVoxelX, nVoxelY, nVoxelZ );
415 fPartialPhantomParam->SetMaterialIndices(fMa << 420 thePartialPhantomParam->SetMaterialIndices( mateIDs );
416 421
417 fPartialPhantomParam->SetFilledIDs(fFilledID << 422 thePartialPhantomParam->SetFilledIDs(theFilledIDs);
418 423
419 fPartialPhantomParam->SetFilledMins(fFilledM << 424 thePartialPhantomParam->SetFilledMins(theFilledMins);
420 425
421 fPartialPhantomParam->BuildContainerWalls(); << 426 thePartialPhantomParam->BuildContainerWalls();
422 427
423 // G4cout << " Number of Materials " << fPh << 428 // G4cout << " Number of Materials " << thePhantomMaterials.size() << G4endl;
424 // G4cout << " SetMaterialIndices(0) " << f << 429 // G4cout << " SetMaterialIndices(0) " << mateIDs[0] << G4endl;
425 430
426 G4PVParameterised* phantom_phys = 0; << 431 G4PVParameterised * phantom_phys = 0;
427 if (OptimAxis == "kUndefined") { << 432 if( OptimAxis == "kUndefined" ) {
428 phantom_phys = new G4PVParameterised(voxel << 433 phantom_phys = new G4PVParameterised(voxelName,phantom_logic,container_logic,
429 fNoVo << 434 kUndefined, theNVoxels, thePartialPhantomParam);
430 } << 435 } else if( OptimAxis == "kXAxis" ) {
431 else if (OptimAxis == "kXAxis") { <<
432 // G4cout << " optim kX " << G4endl; 436 // G4cout << " optim kX " << G4endl;
433 phantom_phys = new G4PVParameterised(voxel << 437 phantom_phys = new G4PVParameterised(voxelName,phantom_logic,container_logic,
434 fNoVo << 438 // kXAxis, nVoxelX*nVoxelY*nVoxelZ, thePartialPhantomParam);
435 fPartialPhantomParam->SetSkipEqualMaterial << 439 kXAxis, theNVoxels, thePartialPhantomParam);
436 } << 440 thePartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
437 else { << 441 } else {
438 G4Exception("GmReadPhantomGeometry::Constr << 442 G4Exception("GmReadPhantomGeometry::ConstructPhantom","Wrong argument to parameter GmReadPhantomGeometry:Phantom:OptimAxis",FatalErrorInArgument,(G4String("Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str()));
439 G4String("Wrong argument to pa <<
440 GmReadPhantomGeometry:Phantom:OptimAx <<
441 Only allowed 'kUndefined' or 'kXAxis' <<
442 + OptimAxis) <<
443 .c_str()); <<
444 } 443 }
445 444
446 phantom_phys->SetRegularStructureId(regStruc << 445 phantom_phys->SetRegularStructureId(regStructureID);
>> 446
>> 447 // G4cout << " GmReadPhantomGeometry phantom_phys " << phantom_phys << " trans " << phantom_phys->GetTranslation() << G4endl;
>> 448 //- thePartialPhantomParam->BuildContainerSolid(cont_phys);
>> 449
>> 450 // G4cout << " GmReadPhantomGeometry::constructPhantom ended " << G4endl;
>> 451
447 } 452 }
448 453