Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-02-23 09:21:45 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file medical/DICOM/src/DicomPartialDetectorConstruction.cc
0028 /// \brief Implementation of the DicomPartialDetectorConstruction class
0029 //
0030 
0031 #include "DicomPartialDetectorConstruction.hh"
0032 
0033 #include "G4Box.hh"
0034 #include "G4Colour.hh"
0035 #include "G4Element.hh"
0036 #include "G4LogicalVolume.hh"
0037 #include "G4Material.hh"
0038 #include "G4NistManager.hh"
0039 #include "G4PVParameterised.hh"
0040 #include "G4PVPlacement.hh"
0041 #include "G4PartialPhantomParameterisation.hh"
0042 #include "G4SystemOfUnits.hh"
0043 #include "G4Tubs.hh"
0044 #include "G4UIcommand.hh"
0045 #include "G4VPhysicalVolume.hh"
0046 #include "G4VisAttributes.hh"
0047 #include "G4ios.hh"
0048 #include "G4tgbMaterialMgr.hh"
0049 #include "G4tgbVolumeMgr.hh"
0050 #include "G4tgrMessenger.hh"
0051 #include "globals.hh"
0052 
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054 DicomPartialDetectorConstruction::DicomPartialDetectorConstruction()
0055   : DicomDetectorConstruction(),
0056     fPartialPhantomParam(0),
0057     fNoVoxels(0),
0058     fDimX(0),
0059     fDimY(0),
0060     fDimZ(0),
0061     fOffsetX(0),
0062     fOffsetY(0),
0063     fOffsetZ(0)
0064 {}
0065 
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067 DicomPartialDetectorConstruction::~DicomPartialDetectorConstruction() {}
0068 
0069 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0070 G4VPhysicalVolume* DicomPartialDetectorConstruction::Construct()
0071 {
0072   InitialisationOfMaterials();
0073 
0074   //----- Build world
0075   G4double worldXDimension = 1. * m;
0076   G4double worldYDimension = 1. * m;
0077   G4double worldZDimension = 1. * m;
0078 
0079   G4Box* world_box = new G4Box("WorldSolid", worldXDimension, worldYDimension, worldZDimension);
0080 
0081   fWorld_logic = new G4LogicalVolume(world_box, fAir, "WorldLogical", 0, 0, 0);
0082 
0083   G4VPhysicalVolume* world_phys =
0084     new G4PVPlacement(0, G4ThreeVector(0, 0, 0), "World", fWorld_logic, 0, false, 0);
0085 
0086   ReadPhantomData();
0087 
0088   ConstructPhantomContainer();
0089   ConstructPhantom();
0090 
0091   return world_phys;
0092 }
0093 
0094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0095 void DicomPartialDetectorConstruction::ReadPhantomData()
0096 {
0097   G4String dataFile = "Data.dat";
0098   std::ifstream finDF(dataFile.c_str());
0099   G4String fname;
0100   if (finDF.good() != 1) {
0101     G4Exception(" DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
0102                 " Problem reading data file: Data.dat");
0103   }
0104 
0105   G4int compression;
0106   finDF >> compression;  // not used here
0107 
0108   G4int numFiles;
0109   finDF >> numFiles;  // only 1 file supported for the moment
0110   if (numFiles != 1) {
0111     G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
0112                 "More than 1 DICOM file is not supported");
0113   }
0114   for (G4int i = 0; i < numFiles; i++) {
0115     finDF >> fname;
0116     //--- Read one data file
0117     fname += ".g4pdcm";
0118     ReadPhantomDataFile(fname);
0119   }
0120 
0121   finDF.close();
0122 }
0123 
0124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0125 void DicomPartialDetectorConstruction::ReadPhantomDataFile(const G4String& fname)
0126 {
0127   //  G4String filename = "phantom.g4pdcm";
0128 #ifdef G4VERBOSE
0129   G4cout << " DicomDetectorConstruction::ReadPhantomDataFile opening file " << fname << G4endl;
0130 #endif
0131   std::ifstream fin(fname.c_str(), std::ios_base::in);
0132   if (!fin.is_open()) {
0133     G4Exception("DicomPartialDetectorConstruction::ReadPhantomDataFile", "", FatalException,
0134                 G4String("File not found " + fname).c_str());
0135   }
0136   G4int nMaterials;
0137   fin >> nMaterials;
0138   G4String stemp;
0139   G4int nmate;
0140   for (G4int ii = 0; ii < nMaterials; ii++) {
0141     fin >> nmate >> stemp;
0142     G4cout << "DicomPartialDetectorConstruction::ReadPhantomData reading nmate " << ii << " = "
0143            << nmate << " mate " << stemp << G4endl;
0144     if (ii != nmate)
0145       G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalErrorInArgument,
0146                   "Material number should be in increasing order: \
0147                           wrong material number ");
0148   }
0149 
0150   fin >> fNoVoxelsX >> fNoVoxelsY >> fNoVoxelsZ;
0151   G4cout << "DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z " << fNoVoxelsX << " "
0152          << fNoVoxelsY << " " << fNoVoxelsZ << G4endl;
0153   fin >> fOffsetX >> fDimX;
0154   fDimX = (fDimX - fOffsetX) / fNoVoxelsX;
0155   fin >> fOffsetY >> fDimY;
0156   fDimY = (fDimY - fOffsetY) / fNoVoxelsY;
0157   fin >> fOffsetZ >> fDimZ;
0158   fDimZ = (fDimZ - fOffsetZ) / fNoVoxelsZ;
0159   G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimX " << fDimX << " fOffsetX "
0160          << fOffsetX << G4endl;
0161   G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimY " << fDimY << " fOffsetY "
0162          << fOffsetY << G4endl;
0163   G4cout << "DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ " << fDimZ << " fOffsetZ "
0164          << fOffsetZ << G4endl;
0165 
0166   //--- Read voxels that are filled
0167   fNoVoxels = 0;
0168   //  G4bool* isFilled = new G4bool[fNoVoxelsX*fNoVoxelsY*fNoVoxelsZ];
0169   //  fFilledIDs = new size_t[fNoVoxelsZ*fNoVoxelsY+1];
0170   //?  fFilledIDs.insert(0);
0171   G4int ifxmin1, ifxmax1;
0172   for (G4int iz = 0; iz < fNoVoxelsZ; iz++) {
0173     std::map<G4int, G4int> ifmin, ifmax;
0174     for (G4int iy = 0; iy < fNoVoxelsY; iy++) {
0175       fin >> ifxmin1 >> ifxmax1;
0176       // check coherence ...
0177 
0178       ifmin[iy] = ifxmin1;
0179       ifmax[iy] = ifxmax1;
0180       G4int ifxdiff = ifxmax1 - ifxmin1 + 1;
0181       if (ifxmax1 == -1 && ifxmin1 == -1) ifxdiff = 0;
0182       fFilledIDs.insert(std::pair<G4int, G4int>(ifxdiff + fNoVoxels - 1, ifxmin1));
0183       G4cout << "DicomPartialDetectorConstruction::ReadPhantomData insert "
0184              << " FilledIDs " << ifxdiff + fNoVoxels - 1 << " min " << ifxmin1
0185              << " N= " << fFilledIDs.size() << G4endl;
0186       // filledIDs[iz*fNoVoxelsY+iy+1] = ifxmax1-ifxmin1+1;
0187       for (G4int ix = 0; ix < fNoVoxelsX; ix++) {
0188         if (ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1)) {
0189           fNoVoxels++;
0190         }
0191         else {
0192         }
0193       }
0194     }
0195     fFilledMins[iz] = ifmin;
0196     fFilledMaxs[iz] = ifmax;
0197   }
0198 
0199   //--- Read material IDs
0200   G4int mateID1;
0201   fMateIDs = new size_t[fNoVoxelsX * fNoVoxelsY * fNoVoxelsZ];
0202   G4int copyNo = 0;
0203   for (G4int iz = 0; iz < fNoVoxelsZ; iz++) {
0204     std::map<G4int, G4int> ifmin = fFilledMins[iz];
0205     std::map<G4int, G4int> ifmax = fFilledMaxs[iz];
0206     for (G4int iy = 0; iy < fNoVoxelsY; iy++) {
0207       ifxmin1 = ifmin[iy];
0208       ifxmax1 = ifmax[iy];
0209       for (G4int ix = 0; ix < fNoVoxelsX; ix++) {
0210         if (ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1)) {
0211           fin >> mateID1;
0212           if (mateID1 < 0 || mateID1 >= nMaterials) {
0213             G4Exception("DicomPartialDetectorConstruction::ReadPhantomData", "", FatalException,
0214                         G4String("Wrong index in phantom file: It should be between 0 and "
0215                                  + G4UIcommand::ConvertToString(nMaterials - 1) + ", while it is "
0216                                  + G4UIcommand::ConvertToString(mateID1))
0217                           .c_str());
0218           }
0219           fMateIDs[copyNo] = mateID1;
0220           copyNo++;
0221         }
0222       }
0223     }
0224   }
0225 
0226   ReadVoxelDensitiesPartial(fin);
0227 
0228   fin.close();
0229 }
0230 
0231 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0232 void DicomPartialDetectorConstruction::ReadVoxelDensitiesPartial(std::ifstream& fin)
0233 {
0234   std::map<G4int, std::pair<G4double, G4double>> densiMinMax;
0235   std::map<G4int, std::pair<G4double, G4double>>::iterator mpite;
0236   for (G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii) {
0237     densiMinMax[ii] = std::pair<G4double, G4double>(DBL_MAX, -DBL_MAX);
0238   }
0239 
0240   //----- Define density differences (maximum density difference to create
0241   // a new material)
0242   char* part = std::getenv("DICOM_CHANGE_MATERIAL_DENSITY");
0243   G4double densityDiff = -1.;
0244   if (part) densityDiff = G4UIcommand::ConvertToDouble(part);
0245   std::map<G4int, G4double> densityDiffs;
0246   if (densityDiff != -1.) {
0247     for (G4int ii = 0; ii < G4int(fOriginalMaterials.size()); ++ii) {
0248       densityDiffs[ii] = densityDiff;  // currently all materials
0249       // with same step
0250     }
0251   }
0252   else {
0253     if (fMaterials.size() == 0) {  // do it only for first slice
0254       for (unsigned int ii = 0; ii < fOriginalMaterials.size(); ++ii) {
0255         fMaterials.push_back(fOriginalMaterials[ii]);
0256       }
0257     }
0258   }
0259   //  densityDiffs[0] = 0.0001; //fAir
0260 
0261   //--- Calculate the average material density for each material/density bin
0262   std::map<std::pair<G4Material*, G4int>, matInfo*> newMateDens;
0263 
0264   G4double dens1;
0265   G4int ifxmin1, ifxmax1;
0266 
0267   //---- Read the material densities
0268   G4int copyNo = 0;
0269   for (G4int iz = 0; iz < fNoVoxelsZ; ++iz) {
0270     std::map<G4int, G4int> ifmin = fFilledMins[iz];
0271     std::map<G4int, G4int> ifmax = fFilledMaxs[iz];
0272     for (G4int iy = 0; iy < fNoVoxelsY; ++iy) {
0273       ifxmin1 = ifmin[iy];
0274       ifxmax1 = ifmax[iy];
0275       for (G4int ix = 0; ix < fNoVoxelsX; ++ix) {
0276         if (ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1)) {
0277           fin >> dens1;
0278           //--- store the minimum and maximum density for each material
0279           //(just for printing)
0280           mpite = densiMinMax.find(G4int(fMateIDs[copyNo]));
0281           if (dens1 < (*mpite).second.first) (*mpite).second.first = dens1;
0282           if (dens1 > (*mpite).second.second) (*mpite).second.second = dens1;
0283 
0284           //--- Get material from original list of material in file
0285           G4int mateID = G4int(fMateIDs[copyNo]);
0286           G4Material* mate = fOriginalMaterials[mateID];
0287           //        G4cout << copyNo << " mateID " << mateID << G4endl;
0288           //--- Check if density is equal to the original material density
0289           if (std::fabs(dens1 - mate->GetDensity() / g * cm3) < 1.e-9) continue;
0290 
0291           //--- Build material name with fOriginalMaterials name + density
0292           //        float densityBin = densityDiffs[mateID] *
0293           //        (G4int(dens1/densityDiffs[mateID])+0.5);
0294           G4String newMateName = mate->GetName();
0295 
0296           G4int densityBin = 0;
0297           //        G4cout << " densityBin " << densityBin << " "
0298           //        << dens1 << " " <<densityDiffs[mateID] << G4endl;
0299 
0300           if (densityDiff != -1.) {
0301             densityBin = (G4int(dens1 / densityDiffs[mateID]));
0302             newMateName = mate->GetName() + G4UIcommand::ConvertToString(densityBin);
0303           }
0304 
0305           //--- Look if it is the first voxel with this material/densityBin
0306           std::pair<G4Material*, G4int> matdens(mate, densityBin);
0307 
0308           auto mppite = newMateDens.find(matdens);
0309           if (mppite != newMateDens.cend()) {
0310             matInfo* mi = (*mppite).second;
0311             mi->fSumdens += dens1;
0312             mi->fNvoxels++;
0313             fMateIDs[copyNo] = fOriginalMaterials.size() - 1 + mi->fId;
0314             //  G4cout << copyNo << " mat new again "
0315             //<< fOriginalMaterials.size()-1 + mi->id << " " << mi->id << G4endl;
0316           }
0317           else {
0318             matInfo* mi = new matInfo;
0319             mi->fSumdens = dens1;
0320             mi->fNvoxels = 1;
0321             mi->fId = G4int(newMateDens.size() + 1);
0322             newMateDens[matdens] = mi;
0323             fMateIDs[copyNo] = fOriginalMaterials.size() - 1 + mi->fId;
0324             //          G4cout << copyNo << " mat new first "
0325             //          << fOriginalMaterials.size()-1 + mi->fId << G4endl;
0326           }
0327           copyNo++;
0328           //        G4cout << ix << " " << iy << " " << iz
0329           //  << " filling fMateIDs " << copyNo << " = " << atoi(cid)-1 << G4endl;
0330           //        fMateIDs[copyNo] = atoi(cid)-1;
0331         }
0332       }
0333     }
0334   }
0335 
0336   //----- Build the list of phantom materials that go to Parameterisation
0337   //--- Add original materials
0338   for (auto mimite = fOriginalMaterials.cbegin(); mimite != fOriginalMaterials.cend(); ++mimite) {
0339     fPhantomMaterials.push_back((*mimite));
0340   }
0341   //
0342   //---- Build and add new materials
0343   for (auto mppite = newMateDens.cbegin(); mppite != newMateDens.cend(); ++mppite) {
0344     G4double averdens = (*mppite).second->fSumdens / (*mppite).second->fNvoxels;
0345     G4double saverdens = G4int(1000.001 * averdens) / 1000.;
0346     G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
0347            << saverdens << " -> " << G4int(1000 * averdens) << " " << G4int(1000 * averdens) / 1000
0348            << " " << G4int(1000 * averdens) / 1000. << G4endl;
0349 
0350     G4cout << "GmReadPhantomGeometry::ReadVoxelDensities AVER DENS " << averdens << " -> "
0351            << saverdens << " -> " << G4UIcommand::ConvertToString(saverdens) << " Nvoxels "
0352            << (*mppite).second->fNvoxels << G4endl;
0353     G4String mateName =
0354       ((*mppite).first).first->GetName() + "_" + G4UIcommand::ConvertToString(saverdens);
0355     fPhantomMaterials.push_back(
0356       BuildMaterialWithChangingDensity((*mppite).first.first, G4float(averdens), mateName));
0357   }
0358 }
0359 
0360 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0361 void DicomPartialDetectorConstruction::ConstructPhantomContainer()
0362 {
0363   // build a clinder that encompass the partial phantom built in the example
0364   //  /dicom/intersectWithUserVolume 0. 0. 0. 0.*deg 0. 0. TUBE 0. 200. 100.
0365   //---- Extract number of voxels and voxel dimensions
0366 
0367   G4String contname = "PhantomContainer";
0368 
0369   //----- Define the volume that contains all the voxels
0370   G4Tubs* container_tube = new G4Tubs(contname, 0., 200., 100., 0., 360 * deg);
0371 
0372   fContainer_logic = new G4LogicalVolume(container_tube, fAir, contname, 0, 0, 0);
0373 
0374   G4ThreeVector posCentreVoxels(0., 0., 0.);
0375   // G4cout << " PhantomContainer posCentreVoxels " << posCentreVoxels << G4endl;
0376   G4RotationMatrix* rotm = new G4RotationMatrix;
0377 
0378   fContainer_phys = new G4PVPlacement(rotm,  // rotation
0379                                       posCentreVoxels,
0380                                       fContainer_logic,  // The logic volume
0381                                       contname,  // Name
0382                                       fWorld_logic,  // Mother
0383                                       false,  // No op. bool.
0384                                       1);  // Copy number
0385 }
0386 
0387 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0388 void DicomPartialDetectorConstruction::ConstructPhantom()
0389 {
0390   G4String OptimAxis = "kXAxis";
0391   G4bool bSkipEqualMaterials = 0;
0392   G4int regStructureID = 2;
0393 
0394   G4ThreeVector posCentreVoxels(0., 0., 0.);
0395 
0396   //----- Build phantom
0397   G4String voxelName = "phantom";
0398   G4Box* phantom_solid = new G4Box(voxelName, fDimX / 2., fDimY / 2., fDimZ / 2.);
0399   G4LogicalVolume* phantom_logic = new G4LogicalVolume(phantom_solid, fAir, voxelName, 0, 0, 0);
0400   G4bool pVis = 0;
0401   if (!pVis) {
0402     G4VisAttributes* visAtt = new G4VisAttributes(FALSE);
0403     phantom_logic->SetVisAttributes(visAtt);
0404   }
0405 
0406   G4double theSmartless = 0.5;
0407   //  fContainer_logic->SetSmartless( theSmartless );
0408   phantom_logic->SetSmartless(theSmartless);
0409 
0410   fPartialPhantomParam = new G4PartialPhantomParameterisation();
0411 
0412   fPartialPhantomParam->SetMaterials(fPhantomMaterials);
0413   fPartialPhantomParam->SetVoxelDimensions(fDimX / 2., fDimY / 2., fDimZ / 2.);
0414   fPartialPhantomParam->SetNoVoxels(fNoVoxelsX, fNoVoxelsY, fNoVoxelsZ);
0415   fPartialPhantomParam->SetMaterialIndices(fMateIDs);
0416 
0417   fPartialPhantomParam->SetFilledIDs(fFilledIDs);
0418 
0419   fPartialPhantomParam->SetFilledMins(fFilledMins);
0420 
0421   fPartialPhantomParam->BuildContainerWalls();
0422 
0423   //  G4cout << " Number of Materials " << fPhantomMaterials.size() << G4endl;
0424   //  G4cout << " SetMaterialIndices(0) " << fMateIDs[0] << G4endl;
0425 
0426   G4PVParameterised* phantom_phys = 0;
0427   if (OptimAxis == "kUndefined") {
0428     phantom_phys = new G4PVParameterised(voxelName, phantom_logic, fContainer_logic, kUndefined,
0429                                          fNoVoxels, fPartialPhantomParam);
0430   }
0431   else if (OptimAxis == "kXAxis") {
0432     //    G4cout << " optim kX " << G4endl;
0433     phantom_phys = new G4PVParameterised(voxelName, phantom_logic, fContainer_logic, kXAxis,
0434                                          fNoVoxels, fPartialPhantomParam);
0435     fPartialPhantomParam->SetSkipEqualMaterials(bSkipEqualMaterials);
0436   }
0437   else {
0438     G4Exception("GmReadPhantomGeometry::ConstructPhantom", "", FatalErrorInArgument,
0439                 G4String("Wrong argument to parameter \
0440          GmReadPhantomGeometry:Phantom:OptimAxis: \
0441          Only allowed 'kUndefined' or 'kXAxis', it is: "
0442                          + OptimAxis)
0443                   .c_str());
0444   }
0445 
0446   phantom_phys->SetRegularStructureId(regStructureID);
0447 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!