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