| 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 // <<
27 /// \file medical/DICOM/src/DicomIntersectVolu <<
28 /// \brief Implementation of the DicomIntersec <<
29 // <<
30 <<
31 #include "DicomIntersectVolume.hh" 26 #include "DicomIntersectVolume.hh"
32 27
>> 28 #include "G4UIcmdWithAString.hh"
33 #include "G4LogicalVolume.hh" 29 #include "G4LogicalVolume.hh"
34 #include "G4LogicalVolumeStore.hh" 30 #include "G4LogicalVolumeStore.hh"
35 #include "G4Material.hh" <<
36 #include "G4PVParameterised.hh" <<
37 #include "G4PhantomParameterisation.hh" <<
38 #include "G4PhysicalVolumeStore.hh" <<
39 #include "G4SystemOfUnits.hh" <<
40 #include "G4UIcmdWithAString.hh" <<
41 #include "G4UIcommand.hh" <<
42 #include "G4VPhysicalVolume.hh" 31 #include "G4VPhysicalVolume.hh"
>> 32 #include "G4PhysicalVolumeStore.hh"
43 #include "G4VSolid.hh" 33 #include "G4VSolid.hh"
44 #include "G4tgbVolume.hh" <<
45 #include "G4tgrSolid.hh" 34 #include "G4tgrSolid.hh"
>> 35 #include "G4tgbVolume.hh"
>> 36 #include "G4Material.hh"
>> 37 #include "G4PhantomParameterisation.hh"
>> 38 #include "G4PVParameterised.hh"
>> 39 #include "G4UIcommand.hh"
46 40
47 //....oooOO0OOooo........oooOO0OOooo........oo << 41 //-----------------------------------------------------------------------
48 DicomIntersectVolume::DicomIntersectVolume() 42 DicomIntersectVolume::DicomIntersectVolume()
49 : G4UImessenger(), fG4VolumeCmd(0), fSolid(0 <<
50 { 43 {
51 fUserVolumeCmd = new G4UIcmdWithAString("/di << 44 theUserVolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
52 fUserVolumeCmd->SetGuidance( << 45 theUserVolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume and outputs the voxels that are totally inside the intersection as a new phantom file. It must have the parameters: POS_X POS_Y POS_Z ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
53 "Intersects a phantom with a user-defined << 46 theUserVolumeCmd->SetParameterName("choice",true);
54 "and outputs the voxels that are totally i << 47 theUserVolumeCmd->AvailableForStates(G4State_Idle);
55 " a new phantom file. It must have the par << 48
56 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_ << 49 theG4VolumeCmd = new G4UIcmdWithAString("/dicom/intersectWithUserVolume",this);
57 fUserVolumeCmd->SetParameterName("choice", t << 50 theG4VolumeCmd->SetGuidance("Intersects a phantom with a user-defined volume and outputs the voxels that are totally inside the intersection as a new phantom file. It must have the parameters: POS_X POS_Y POS_Z ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_1 (SOLID_PARAM_2 ...)");
58 fUserVolumeCmd->AvailableForStates(G4State_I << 51 theG4VolumeCmd->SetParameterName("choice",true);
59 << 52 theG4VolumeCmd->AvailableForStates(G4State_Idle);
60 fG4VolumeCmd = new G4UIcmdWithAString("/dico << 53
61 fG4VolumeCmd->SetGuidance( <<
62 "Intersects a phantom with a user-defined <<
63 " and outputs the voxels that are totally <<
64 " a new phantom file. It must have the par <<
65 "ANG_X ANG_Y ANG_Z SOLID_TYPE SOLID_PARAM_ <<
66 fG4VolumeCmd->SetParameterName("choice", tru <<
67 fG4VolumeCmd->AvailableForStates(G4State_Idl <<
68 } 54 }
69 55
70 //....oooOO0OOooo........oooOO0OOooo........oo << 56 //-----------------------------------------------------------------------
71 DicomIntersectVolume::~DicomIntersectVolume() 57 DicomIntersectVolume::~DicomIntersectVolume()
72 { 58 {
73 if (fUserVolumeCmd) delete fUserVolumeCmd; << 59 if (theUserVolumeCmd) delete theUserVolumeCmd;
74 if (fG4VolumeCmd) delete fG4VolumeCmd; << 60 if (theG4VolumeCmd) delete theG4VolumeCmd;
75 } 61 }
76 62
77 //....oooOO0OOooo........oooOO0OOooo........oo << 63 //---------------------------------------------------------------------
78 void DicomIntersectVolume::SetNewValue(G4UIcom << 64 void DicomIntersectVolume::SetNewValue(G4UIcommand * command,
79 { << 65 G4String newValues)
>> 66 {
80 G4AffineTransform theVolumeTransform; 67 G4AffineTransform theVolumeTransform;
81 68
82 if (command == fUserVolumeCmd) { << 69 if (command == theUserVolumeCmd) {
83 std::vector<G4String> params = GetWordsInS << 70
84 if (params.size() < 8) { << 71 std::vector<G4String> params = GetWordsInString( newValues );
85 G4Exception("DicomIntersectVolume::SetNe << 72 if( params.size() < 8 ) {
86 " There must be at least 8 p << 73 G4Exception("GmPositionUserVolumes::SetParams"," There must be at least 8 parameter: SOLID_TYPE POS_X POS_Y POS_Z ANG_X ANG_Y ANG_Z SOLID_PARAM_1 (SOLID_PARAM_2 ...)",FatalErrorInArgument,G4String("Number of parameters given = " + G4UIcommand::ConvertToString( G4int(params.size()) )).c_str());
87 "ANG_X ANG_Y ANG_Z SOLID_PAR << 74
88 FatalErrorInArgument, <<
89 G4String("Number of paramete <<
90 + G4UIcommand::Conv <<
91 .c_str()); <<
92 } 75 }
93 76
94 //----- Build G4VSolid << 77 //----- Build G4VSolid
95 BuildUserSolid(params); 78 BuildUserSolid(params);
96 79
97 //----- Calculate volume inverse 3D transf 80 //----- Calculate volume inverse 3D transform
98 G4ThreeVector pos = G4ThreeVector(G4UIcomm << 81 G4ThreeVector pos = G4ThreeVector( G4UIcommand::ConvertToDouble(params[0]), G4UIcommand::ConvertToDouble(params[1]), G4UIcommand::ConvertToDouble(params[2]) );
99 G4UIcomm <<
100 G4UIcomm <<
101 G4RotationMatrix* rotmat = new G4RotationM 82 G4RotationMatrix* rotmat = new G4RotationMatrix;
102 std::vector<G4double> angles; << 83 std::vector<double> angles;
103 rotmat->rotateX(G4UIcommand::ConvertToDoub << 84 rotmat->rotateX( G4UIcommand::ConvertToDouble(params[3]) );
104 rotmat->rotateY(G4UIcommand::ConvertToDoub << 85 rotmat->rotateY( G4UIcommand::ConvertToDouble(params[4]) );
105 rotmat->rotateY(G4UIcommand::ConvertToDoub << 86 rotmat->rotateY( G4UIcommand::ConvertToDouble(params[5]) );
106 theVolumeTransform = G4AffineTransform(rot << 87 theVolumeTransform = G4AffineTransform( rotmat, pos ).Invert();
107 } << 88
108 else if (command == fG4VolumeCmd) { << 89 } else if (command == theG4VolumeCmd) {
109 std::vector<G4String> params = GetWordsInS << 90 std::vector<G4String> params = GetWordsInString( newValues );
110 if (params.size() != 1) << 91 if( params.size() !=1 ) G4Exception("Command: "+ command->GetCommandPath() + "/" + command->GetCommandName() + " " + newValues + " needs 1 argument: VOLUME_NAME");
111 G4Exception("DicomIntersectVolume::SetNe <<
112 G4String("Command: " + comma <<
113 + " " + newValues + <<
114 .c_str()); <<
115 92
116 //----- Build G4VSolid << 93 //----- Build G4VSolid
117 BuildG4Solid(params); 94 BuildG4Solid(params);
118 95
119 //----- Calculate volume inverse 3D transf 96 //----- Calculate volume inverse 3D transform
120 G4VPhysicalVolume* pv = GetPhysicalVolumes << 97 G4VPhysicalVolume* pv = GetPhysicalVolumes( params[0], 1, 1)[0];
121 << 98
122 theVolumeTransform = G4AffineTransform(pv- << 99 theVolumeTransform = G4AffineTransform( pv->GetFrameRotation(), pv->GetFrameTranslation() );
123 } 100 }
124 101
125 //----- Calculate relative phantom - volume << 102 //----- Calculate relative phantom - volume 3D transform
126 G4PhantomParameterisation* thePhantomParam = 103 G4PhantomParameterisation* thePhantomParam = GetPhantomParam(true);
127 << 104
128 G4RotationMatrix* rotph = new G4RotationMatr << 105 G4RotationMatrix* rotph = new G4RotationMatrix(); // assumes the phantom mother is not rotated !!!
129 // assumes the phantom mother is not rotated << 106 G4AffineTransform thePhantomTransform( rotph, G4ThreeVector() ); // assumes the phantom mother is not translated !!!
130 G4AffineTransform thePhantomTransform(rotph, << 107
131 // assumes the phantom mother is not transla << 108 G4AffineTransform theTransform = theVolumeTransform*thePhantomTransform;
132 << 109
133 G4AffineTransform theTransform = theVolumeTr << 110 G4ThreeVector axisX( 1., 0., 0. );
134 << 111 G4ThreeVector axisY( 0., 1., 0. );
135 G4ThreeVector axisX(1., 0., 0.); << 112 G4ThreeVector axisZ( 0., 0., 1. );
136 G4ThreeVector axisY(0., 1., 0.); <<
137 G4ThreeVector axisZ(0., 0., 1.); <<
138 theTransform.ApplyAxisTransform(axisX); 113 theTransform.ApplyAxisTransform(axisX);
139 theTransform.ApplyAxisTransform(axisY); 114 theTransform.ApplyAxisTransform(axisY);
140 theTransform.ApplyAxisTransform(axisZ); 115 theTransform.ApplyAxisTransform(axisZ);
141 << 116
142 //----- Write phantom header 117 //----- Write phantom header
143 G4String thePhantomFileName = "phantom.g4pdc 118 G4String thePhantomFileName = "phantom.g4pdcm";
144 fout.open(thePhantomFileName); 119 fout.open(thePhantomFileName);
145 std::vector<G4Material*> materials = thePhan 120 std::vector<G4Material*> materials = thePhantomParam->GetMaterials();
146 fout << materials.size() << G4endl; 121 fout << materials.size() << G4endl;
147 for (unsigned int ii = 0; ii < materials.siz << 122 for( unsigned int ii = 0; ii < materials.size(); ii++ ) {
148 fout << ii << " " << materials[ii]->GetNam 123 fout << ii << " " << materials[ii]->GetName() << G4endl;
149 } 124 }
150 125
151 //----- Loop to pantom voxels 126 //----- Loop to pantom voxels
152 G4int nx = G4int(thePhantomParam->GetNoVoxel << 127 G4int nx = thePhantomParam->GetNoVoxelX();
153 G4int ny = G4int(thePhantomParam->GetNoVoxel << 128 G4int ny = thePhantomParam->GetNoVoxelY();
154 G4int nz = G4int(thePhantomParam->GetNoVoxel << 129 G4int nz = thePhantomParam->GetNoVoxelZ();
155 G4int nxy = nx * ny; << 130 G4int nxy = nx*ny;
156 fVoxelIsInside = new G4bool[nx * ny * nz]; << 131 theVoxelIsInside = new G4bool[nx*ny*nz];
157 G4double voxelHalfWidthX = thePhantomParam-> 132 G4double voxelHalfWidthX = thePhantomParam->GetVoxelHalfX();
158 G4double voxelHalfWidthY = thePhantomParam-> 133 G4double voxelHalfWidthY = thePhantomParam->GetVoxelHalfY();
159 G4double voxelHalfWidthZ = thePhantomParam-> 134 G4double voxelHalfWidthZ = thePhantomParam->GetVoxelHalfZ();
160 135
161 //----- Write phantom dimensions and limits 136 //----- Write phantom dimensions and limits
162 fout << nx << " " << ny << " " << nz << G4en 137 fout << nx << " " << ny << " " << nz << G4endl;
163 fout << -voxelHalfWidthX * nx + thePhantomTr << 138 fout << -voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << " "
164 << voxelHalfWidthX * nx + thePhantomTra << 139 << voxelHalfWidthX*nx+thePhantomTransform.NetTranslation().x() << G4endl;
165 fout << -voxelHalfWidthY * ny + thePhantomTr << 140 fout << -voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << " "
166 << voxelHalfWidthY * ny + thePhantomTra << 141 << voxelHalfWidthY*ny+thePhantomTransform.NetTranslation().y() << G4endl;
167 fout << -voxelHalfWidthZ * nz + thePhantomTr << 142 fout << -voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << " "
168 << voxelHalfWidthZ * nz + thePhantomTra << 143 << voxelHalfWidthZ*nz+thePhantomTransform.NetTranslation().z() << G4endl;
>> 144
>> 145 for( G4int iz = 0; iz < nz; iz++ ){
>> 146 for( G4int iy = 0; iy < ny; iy++) {
169 147
170 for (G4int iz = 0; iz < nz; ++iz) { <<
171 for (G4int iy = 0; iy < ny; ++iy) { <<
172 G4bool bPrevVoxelInside = true; 148 G4bool bPrevVoxelInside = true;
173 G4bool b1VoxelFoundInside = false; 149 G4bool b1VoxelFoundInside = false;
174 G4int firstVoxel = -1; 150 G4int firstVoxel = -1;
175 G4int lastVoxel = -1; 151 G4int lastVoxel = -1;
176 for (G4int ix = 0; ix < nx; ++ix) { << 152 for(G4int ix = 0; ix < nx; ix++ ){
177 G4ThreeVector voxelCentre((-nx + ix * << 153 G4ThreeVector voxelCentre( (-nx+ix*2+1)*voxelHalfWidthX, (-ny+iy*2+1)*voxelHalfWidthY, (-nz+iz*2+1)*voxelHalfWidthZ);
178 (-ny + iy * << 154 theTransform.ApplyPointTransform(voxelCentre);
179 (-nz + iz * << 155 G4bool bVoxelIsInside = true;
180 theTransform.ApplyPointTransform(voxel << 156 for( G4int ivx = -1; ivx <= 1; ivx+=2 ) {
181 G4bool bVoxelIsInside = true; << 157 for( G4int ivy = -1; ivy <= 1; ivy+=2 ){
182 for (G4int ivx = -1; ivx <= 1; ivx += << 158 for( G4int ivz = -1; ivz <= 1; ivz+=2 ) {
183 for (G4int ivy = -1; ivy <= 1; ivy + << 159 G4ThreeVector voxelPoint = voxelCentre + ivx*voxelHalfWidthX*axisX + ivy*voxelHalfWidthY*axisY + ivz*voxelHalfWidthZ*axisZ;
184 for (G4int ivz = -1; ivz <= 1; ivz << 160 if( theSolid->Inside( voxelPoint ) == kOutside ) {
185 G4ThreeVector voxelPoint = voxel << 161 bVoxelIsInside = false;
186 + ivy << 162 break;
187 + ivz << 163 } else {
188 if (fSolid->Inside(voxelPoint) = << 164 }
189 bVoxelIsInside = false; << 165 }
190 break; << 166 if( !bVoxelIsInside ) break;
191 } << 167 }
192 else { << 168 if( !bVoxelIsInside ) break;
193 } << 169 }
194 } << 170
195 if (!bVoxelIsInside) break; << 171 G4int copyNo = ix + nx*iy + nxy*iz;
196 } << 172 if( bVoxelIsInside ) {
197 if (!bVoxelIsInside) break; << 173 if( !bPrevVoxelInside ) {
198 } << 174 G4Exception("G4PhantomIntersectVolume",
199 << 175 "Volume cannot intersect phantom in discontiguous voxels, please use other voxel",
200 G4int copyNo = ix + nx * iy + nxy * iz << 176 FatalException,
201 if (bVoxelIsInside) { << 177 "");
202 if (!bPrevVoxelInside) { << 178 }
203 G4Exception("DicomIntersectVolume: << 179 if( !b1VoxelFoundInside ) {
204 "Volume cannot interse << 180 firstVoxel = ix;
205 "please use other voxe << 181 b1VoxelFoundInside = true;
206 } << 182 }
207 if (!b1VoxelFoundInside) { << 183 lastVoxel = ix;
208 firstVoxel = ix; << 184 theVoxelIsInside[copyNo] = true;
209 b1VoxelFoundInside = true; << 185 } else {
210 } << 186 theVoxelIsInside[copyNo] = false;
211 lastVoxel = ix; << 187 }
212 fVoxelIsInside[copyNo] = true; << 188
213 } <<
214 else { <<
215 fVoxelIsInside[copyNo] = false; <<
216 } <<
217 } 189 }
218 fout << firstVoxel << " " << lastVoxel < 190 fout << firstVoxel << " " << lastVoxel << G4endl;
219 } 191 }
220 } 192 }
221 193
222 //----- Now write the materials 194 //----- Now write the materials
223 for (G4int iz = 0; iz < nz; ++iz) { << 195 for( G4int iz = 0; iz < nz; iz++ ){
224 for (G4int iy = 0; iy < ny; ++iy) { << 196 for( G4int iy = 0; iy < ny; iy++) {
225 G4bool b1xFound = false; << 197 G4bool b1xFound = false;
226 for (G4int ix = 0; ix < nx; ++ix) { << 198 for(G4int ix = 0; ix < nx; ix++ ){
227 size_t copyNo = ix + ny * iy + nxy * i << 199 size_t copyNo = ix + ny*iy + nxy*iz;
228 if (fVoxelIsInside[copyNo]) { << 200 // fout << " iz " << iz << " i " << iy << " ix " << ix << G4endl;
229 fout << thePhantomParam->GetMaterial << 201 if( theVoxelIsInside[copyNo] ) {
230 b1xFound = true; << 202 fout << thePhantomParam->GetMaterialIndex(copyNo)<< " ";
231 } << 203 b1xFound = true;
>> 204 }
232 } 205 }
233 if (b1xFound) fout << G4endl; << 206 if(b1xFound ) fout << G4endl;
234 } 207 }
235 } 208 }
236 209
237 // Now write densities 210 // Now write densities
238 for (G4int iz = 0; iz < nz; ++iz) { << 211 for( G4int iz = 0; iz < nz; iz++ ){
239 for (G4int iy = 0; iy < ny; ++iy) { << 212 for( G4int iy = 0; iy < ny; iy++) {
240 G4bool b1xFound = false; << 213 G4bool b1xFound = false;
241 for (G4int ix = 0; ix < nx; ++ix) { << 214 for(G4int ix = 0; ix < nx; ix++ ){
242 size_t copyNo = ix + ny * iy + nxy * i << 215 size_t copyNo = ix + ny*iy + nxy*iz;
243 if (fVoxelIsInside[copyNo]) { << 216 if( theVoxelIsInside[copyNo] ) {
244 fout << thePhantomParam->GetMaterial << 217 fout << thePhantomParam->GetMaterial(copyNo)->GetDensity()/g*cm3 << " ";
245 b1xFound = true; << 218 b1xFound = true;
246 } << 219 }
247 } 220 }
248 if (b1xFound) fout << G4endl; << 221 if(b1xFound ) fout << G4endl;
249 } 222 }
250 } 223 }
251 } <<
252 224
253 //....oooOO0OOooo........oooOO0OOooo........oo << 225 }
254 void DicomIntersectVolume::BuildUserSolid(std: << 226
>> 227 //---------------------------------------------------------------------
>> 228 void DicomIntersectVolume::BuildUserSolid( std::vector<G4String> params )
255 { 229 {
256 for (G4int ii = 0; ii < 6; ++ii) << 230 for( G4int ii = 0; ii < 6; ii++ ) params.erase( params.begin() ); // take otu position and rotation angles
257 params.erase(params.begin()); << 231 params.insert( params.begin(), ":SOLID");
258 // take otu position and rotation angles << 232 params.insert( params.begin(), params[1] );
259 params.insert(params.begin(), ":SOLID"); <<
260 params.insert(params.begin(), params[1]); <<
261 G4tgrSolid* tgrSolid = new G4tgrSolid(params 233 G4tgrSolid* tgrSolid = new G4tgrSolid(params);
262 G4tgbVolume* tgbVolume = new G4tgbVolume(); << 234 G4tgbVolume* tgbVolume = new G4tgbVolume();
263 fSolid = tgbVolume->FindOrConstructG4Solid(t << 235 theSolid = tgbVolume->FindOrConstructG4Solid( tgrSolid );
>> 236
264 } 237 }
265 238
266 //....oooOO0OOooo........oooOO0OOooo........oo << 239 //---------------------------------------------------------------------
267 void DicomIntersectVolume::BuildG4Solid(std::v << 240 void DicomIntersectVolume::BuildG4Solid( std::vector<G4String> params )
268 { 241 {
269 fSolid = GetLogicalVolumes(params[0], 1, 1)[ << 242 theSolid = GetLogicalVolumes( params[0], 1, 1)[0]->GetSolid();
>> 243
270 } 244 }
271 245
272 //....oooOO0OOooo........oooOO0OOooo........oo << 246
>> 247
>> 248 //-----------------------------------------------------------------------
273 G4PhantomParameterisation* DicomIntersectVolum 249 G4PhantomParameterisation* DicomIntersectVolume::GetPhantomParam(G4bool bMustExist)
274 { 250 {
275 G4PhantomParameterisation* paramreg = 0; 251 G4PhantomParameterisation* paramreg = 0;
276 252
277 G4PhysicalVolumeStore* pvs = G4PhysicalVolum 253 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
278 for (auto cite = pvs->cbegin(); cite != pvs- << 254 std::vector<G4VPhysicalVolume*>::iterator cite;
279 if (IsPhantomVolume(*cite)) { << 255 for( cite = pvs->begin(); cite != pvs->end(); cite++ ) {
>> 256 // G4cout << " PV " << (*cite)->GetName() << " " << (*cite)->GetTranslation() << G4endl;
>> 257 if( IsPhantomVolume( *cite ) ) {
280 const G4PVParameterised* pvparam = stati 258 const G4PVParameterised* pvparam = static_cast<const G4PVParameterised*>(*cite);
281 G4VPVParameterisation* param = pvparam-> 259 G4VPVParameterisation* param = pvparam->GetParameterisation();
282 // if( static_cast<const G4PhantomPar 260 // if( static_cast<const G4PhantomParameterisation*>(param) ){
283 // if( static_cast<const G4PhantomPar 261 // if( static_cast<const G4PhantomParameterisation*>(param) ){
284 // G4cout << "G4PhantomParameterisa << 262 // G4cout << "G4PhantomParameterisation volume found " << (*cite)->GetName() << G4endl;
285 // << (*cite)->GetName() << G4endl; <<
286 paramreg = static_cast<G4PhantomParamete 263 paramreg = static_cast<G4PhantomParameterisation*>(param);
287 } 264 }
288 } 265 }
289 << 266
290 if (!paramreg && bMustExist) << 267 if( !paramreg && bMustExist ) G4Exception("DicomIntersectVolume::GetPhantomParam: No G4PhantomParameterisation found ");
291 G4Exception("DicomIntersectVolume::GetPhan << 268
292 " No G4PhantomParameterisation <<
293 <<
294 return paramreg; 269 return paramreg;
>> 270
295 } 271 }
296 272
297 //....oooOO0OOooo........oooOO0OOooo........oo << 273
298 std::vector<G4VPhysicalVolume*> DicomIntersect << 274
299 << 275 //-----------------------------------------------------------------------
>> 276 std::vector<G4VPhysicalVolume*> DicomIntersectVolume::GetPhysicalVolumes( const G4String& name, bool exists, G4int nVols )
300 { 277 {
301 std::vector<G4VPhysicalVolume*> vvolu; 278 std::vector<G4VPhysicalVolume*> vvolu;
302 std::string::size_type ial = name.rfind(":") 279 std::string::size_type ial = name.rfind(":");
303 G4String volname = ""; 280 G4String volname = "";
304 G4int volcopy = 0; 281 G4int volcopy = 0;
305 if (ial != std::string::npos) { << 282 if( ial != std::string::npos ) {
306 std::string::size_type ial2 = name.rfind(" << 283 std::string::size_type ial2 = name.rfind(":",ial-1);
307 if (ial2 != std::string::npos) { << 284 if( ial2 != std::string::npos ) {
308 G4Exception("DicomIntersectVolume::GetPh << 285 G4Exception("DicomIntersectVolume::GetPhysicalVolumes: name corresponds to a touchable " + name);
309 G4String("Name corresponds t << 286 }else {
310 } << 287 volname = name.substr( 0, ial );
311 else { << 288 volcopy = G4UIcommand::ConvertToInt( name.substr( ial+1, name.length() ).c_str() );
312 volname = name.substr(0, ial); <<
313 volcopy = G4UIcommand::ConvertToInt(name <<
314 } 289 }
315 } << 290 } else {
316 else { <<
317 volname = name; 291 volname = name;
318 volcopy = -1; 292 volcopy = -1;
319 } 293 }
320 294
321 G4PhysicalVolumeStore* pvs = G4PhysicalVolum 295 G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
322 for (auto citepv = pvs->cbegin(); citepv != << 296 std::vector<G4VPhysicalVolume*>::iterator citepv;
323 if (volname == (*citepv)->GetName() && ((* << 297 for( citepv = pvs->begin(); citepv != pvs->end(); citepv++ ) {
324 vvolu.push_back(*citepv); << 298 if( volname == (*citepv)->GetName()
>> 299 && ( (*citepv)->GetCopyNo() == volcopy || -1 == volcopy ) ){
>> 300 vvolu.push_back( *citepv );
325 } 301 }
326 } 302 }
327 303
328 //----- Check that at least one volume was f 304 //----- Check that at least one volume was found
329 if (vvolu.size() == 0) { << 305 if( vvolu.size() == 0 ) {
330 if (exists) { << 306 if(exists) {
331 G4Exception(" DicomIntersectVolume::GetP << 307 G4Exception(" DicomIntersectVolume::GetPhysicalVolumes: no physical volume found with name " + name);
332 G4String("No physical volume << 308 } else {
333 } << 309 G4cerr << "!!WARNING: DicomIntersectVolume::GetPhysicalVolumes: no physical volume found with name " << name << G4endl;
334 else { <<
335 G4cerr << "!!WARNING: DicomIntersectVolu <<
336 << " no physical volume found wit <<
337 } 310 }
338 } 311 }
339 312
340 if (nVols != -1 && G4int(vvolu.size()) != nV << 313 if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
341 G4Exception("DicomIntersectVolume::GetLogi 314 G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
342 "Wrong number of physical volu << 315 "Wrong number of physical volumes found",
343 ("Number of physical volumes " << 316 FatalErrorInArgument,
344 + ", requesting " + G4UIcomma << 317 ("Number of physical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size())) + ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
345 .c_str()); << 318 }
346 } <<
347 319
348 return vvolu; 320 return vvolu;
349 } 321 }
350 322
351 //....oooOO0OOooo........oooOO0OOooo........oo << 323
352 G4bool DicomIntersectVolume::IsPhantomVolume(G << 324 //-----------------------------------------------------------------------
>> 325 G4bool DicomIntersectVolume::IsPhantomVolume( G4VPhysicalVolume* pv )
353 { 326 {
354 EAxis axis; 327 EAxis axis;
355 G4int nReplicas; 328 G4int nReplicas;
356 G4double width, offset; << 329 G4double width,offset;
357 G4bool consuming; 330 G4bool consuming;
358 pv->GetReplicationData(axis, nReplicas, widt << 331 pv->GetReplicationData(axis,nReplicas,width,offset,consuming);
359 EVolume type = (consuming) ? kReplica : kPar 332 EVolume type = (consuming) ? kReplica : kParameterised;
360 if (type == kParameterised && pv->GetRegular << 333 if( type == kParameterised && pv->GetRegularStructureId() == 1 ) {
361 return TRUE; 334 return TRUE;
362 } << 335 } else {
363 else { <<
364 return FALSE; 336 return FALSE;
365 } 337 }
366 } <<
367 338
368 //....oooOO0OOooo........oooOO0OOooo........oo << 339 }
369 std::vector<G4LogicalVolume*> DicomIntersectVo << 340
370 << 341
>> 342 //-----------------------------------------------------------------------
>> 343 std::vector<G4LogicalVolume*> DicomIntersectVolume::GetLogicalVolumes( const G4String& name, bool exists, G4int nVols )
371 { 344 {
372 // G4cout << "GetLogicalVolumes " << name < 345 // G4cout << "GetLogicalVolumes " << name << " " << exists << G4endl;
373 std::vector<G4LogicalVolume*> vvolu; 346 std::vector<G4LogicalVolume*> vvolu;
374 G4int ial = G4int(name.rfind(":")); << 347 G4int ial = name.rfind(":");
375 if (ial != -1) { << 348 if( ial != -1 ) {
376 G4Exception("DicomIntersectVolume::GetLogi << 349 G4Exception("DicomIntersectVolume::GetLogicalVolumes: name corresponds to a touchable or physcal volume" + name);
377 G4String("Name corresponds to <<
378 } 350 }
379 351
380 G4LogicalVolumeStore* lvs = G4LogicalVolumeS 352 G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
381 for (auto citelv = lvs->cbegin(); citelv != << 353 std::vector<G4LogicalVolume*>::iterator citelv;
382 if (name == (*citelv)->GetName()) { << 354 for( citelv = lvs->begin(); citelv != lvs->end(); citelv++ ) {
383 vvolu.push_back(*citelv); << 355 if( name == (*citelv)->GetName() ) {
>> 356 vvolu.push_back( *citelv );
384 } 357 }
385 } 358 }
386 << 359
387 //----- Check that at least one volume was f 360 //----- Check that at least one volume was found
388 if (vvolu.size() == 0) { << 361 if( vvolu.size() == 0 ) {
389 if (exists) { << 362 if(exists) {
390 G4Exception("DicomIntersectVolume::GetLo << 363 G4Exception("DicomIntersectVolume::GetLogicalVolumes:","ERROR",FatalErrorInArgument,("no logical volume found with name " + name).c_str());
391 ("no logical volume found wi << 364 } else {
392 } << 365 G4Exception("DicomIntersectVolume::GetLogicalVolumes:","ERROR",JustWarning,("no logical volume found with name " + name).c_str());
393 else { <<
394 G4Exception("DicomIntersectVolume::GetLo <<
395 ("no logical volume found w <<
396 } 366 }
397 } 367 }
398 368
399 if (nVols != -1 && G4int(vvolu.size()) != nV << 369 if( nVols != -1 && G4int(vvolu.size()) != nVols ) {
400 G4Exception("DicomIntersectVolume::GetLogi << 370 G4Exception("DicomIntersectVolume::GetLogicalVolumes:",
401 FatalErrorInArgument, << 371 "Wrong number of logical volumes found",
402 ("Number of logical volumes " << 372 FatalErrorInArgument,
403 + ", requesting " + G4UIcomma << 373 ("Number of logical volumes " + G4UIcommand::ConvertToString(G4int(vvolu.size())) + ", requesting " + G4UIcommand::ConvertToString(nVols)).c_str());
404 .c_str()); << 374 }
405 } <<
406 375
407 return vvolu; 376 return vvolu;
>> 377
408 } 378 }
409 379
410 //....oooOO0OOooo........oooOO0OOooo........oo << 380
411 std::vector<G4String> DicomIntersectVolume::Ge << 381 //------------------------------------------------------------------------
>> 382 std::vector<G4String> DicomIntersectVolume::GetWordsInString( const G4String& stemp)
412 { 383 {
413 std::vector<G4String> wordlist; 384 std::vector<G4String> wordlist;
414 385
415 //---------- Read a line of file: 386 //---------- Read a line of file:
416 //----- Clear wordlist 387 //----- Clear wordlist
417 G4int ii; 388 G4int ii;
418 const char* cstr = stemp.c_str(); 389 const char* cstr = stemp.c_str();
419 G4int siz = G4int(stemp.length()); << 390 int siz = stemp.length();
420 G4int istart = 0; << 391 int istart = 0;
421 G4int nQuotes = 0; << 392 int nQuotes = 0;
422 G4bool lastIsBlank = false; << 393 bool lastIsBlank = false;
423 G4bool lastIsQuote = false; << 394 bool lastIsQuote = false;
424 for (ii = 0; ii < siz; ++ii) { << 395 for( ii = 0; ii < siz; ii++ ){
425 if (cstr[ii] == '\"') { << 396 if(cstr[ii] == '\"' ){
426 if (lastIsQuote) { << 397 if( lastIsQuote ){
427 G4Exception("GmGenUtils:GetWordsFromSt << 398 G4Exception("GmGenUtils:GetWordsFromString","ERROR",FatalException, ("There cannot be two quotes together " + stemp).c_str() );
428 ("There cannot be two quot <<
429 } <<
430 if (nQuotes % 2 == 1) { <<
431 // close word <<
432 wordlist.push_back(stemp.substr(istart <<
433 // G4cout << "GetWordsInString <<
434 //<< wordlist[wordlist.size()-1] << " <<
435 //<< " ii " << ii << G4endl; <<
436 } 399 }
437 else { << 400 if( nQuotes%2 == 1 ){
438 istart = ii + 1; << 401 //close word
>> 402 wordlist.push_back( stemp.substr(istart,ii-istart) );
>> 403 // G4cout << "GetWordsInString new word0 " << wordlist[wordlist.size()-1] << " istart " << istart << " ii " << ii << G4endl;
>> 404 } else {
>> 405 istart = ii+1;
439 } 406 }
440 ++nQuotes; << 407 nQuotes++;
441 lastIsQuote = true; 408 lastIsQuote = true;
442 lastIsBlank = false; 409 lastIsBlank = false;
443 } << 410 } else if(cstr[ii] == ' ' ){
444 else if (cstr[ii] == ' ') { << 411 // G4cout << "GetWordsInString blank nQuotes " << nQuotes << " lastIsBlank " << lastIsBlank << G4endl;
445 // G4cout << "GetWordsInStrin << 412 if( nQuotes%2 == 0 ){
446 //<< nQuotes << " lastIsBlank " << lastI << 413 if( !lastIsBlank && !lastIsQuote ) {
447 if (nQuotes % 2 == 0) { << 414 wordlist.push_back( stemp.substr(istart,ii-istart) );
448 if (!lastIsBlank && !lastIsQuote) { << 415 // G4cout << "GetWordsInString new word1 " << wordlist[wordlist.size()-1] << " lastIsBlank " << lastIsBlank << G4endl;
449 wordlist.push_back(stemp.substr(ista << 416 }
450 // G4cout << "GetWordsInStr << 417
451 //<< wordlist[wordlist.size()-1] << << 418 istart = ii+1;
452 //<< lastIsBlank << G4endl; << 419 lastIsQuote = false;
453 } << 420 lastIsBlank = true;
454 <<
455 istart = ii + 1; <<
456 lastIsQuote = false; <<
457 lastIsBlank = true; <<
458 } 421 }
459 } << 422 } else {
460 else { << 423 if( ii == siz-1 ) {
461 if (ii == siz - 1) { << 424 wordlist.push_back( stemp.substr(istart,ii-istart+1) );
462 wordlist.push_back(stemp.substr(istart << 425 // G4cout << "GetWordsInString new word2 " << wordlist[wordlist.size()-1] << " istart " << istart << G4endl;
463 // G4cout << "GetWordsI <<
464 //<< wordlist[wordlist.size()-1] << " <<
465 } 426 }
466 lastIsQuote = false; 427 lastIsQuote = false;
467 lastIsBlank = false; 428 lastIsBlank = false;
468 } 429 }
469 } 430 }
470 if (nQuotes % 2 == 1) { << 431 if( nQuotes%2 == 1 ) {
471 G4Exception("GmGenUtils:GetWordsFromString << 432 G4Exception("GmGenUtils:GetWordsFromString","ERROR",FatalException, ("unbalanced quotes in line " + stemp).c_str() );
472 ("unbalanced quotes in line " <<
473 } 433 }
474 434
475 return wordlist; 435 return wordlist;
476 } 436 }
>> 437
477 438