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