152 G4cout <<
" DicomDetectorConstruction::ReadPhantomDataFile opening file "
155 std::ifstream fin(fname.c_str(), std::ios_base::in);
156 if( !fin.is_open() ) {
157 G4Exception(
"DicomPartialDetectorConstruction::ReadPhantomDataFile",
160 G4String(
"File not found " + fname).c_str());
166 for( G4int ii = 0; ii < nMaterials; ii++ ){
167 fin >> nmate >> stemp;
168 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData reading nmate "
169 << ii <<
" = " << nmate <<
" mate " << stemp << G4endl;
171 G4Exception(
"DicomPartialDetectorConstruction::ReadPhantomData",
173 FatalErrorInArgument,
174 "Material number should be in increasing order: \
175 wrong material number ");
180 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData nVoxel X/Y/Z "
188 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData voxelDimX "
190 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData voxelDimY "
192 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData voxelDimZ "
200 G4int ifxmin1, ifxmax1;
202 std::map< G4int, G4int > ifmin, ifmax;
204 fin >> ifxmin1 >> ifxmax1;
209 G4int ifxdiff = ifxmax1-ifxmin1+1;
210 if( ifxmax1 == -1 && ifxmin1 == -1 ) ifxdiff = 0;
212 G4cout <<
"DicomPartialDetectorConstruction::ReadPhantomData insert "
213 <<
" FilledIDs " << ifxdiff+
fNoVoxels-1 <<
" min " << ifxmin1
217 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
238 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
240 if( mateID1 < 0 || mateID1 >= nMaterials ) {
241 G4Exception(
"DicomPartialDetectorConstruction::ReadPhantomData",
244 G4String(
"Wrong index in phantom file: It should be between 0 and "
245 + G4UIcommand::ConvertToString(nMaterials-1)
247 + G4UIcommand::ConvertToString(mateID1)).c_str());
265 std::map<G4int, std::pair<G4double,G4double> > densiMinMax;
266 std::map<G4int, std::pair<G4double,G4double> >::iterator mpite;
269 densiMinMax[ii] = std::pair<G4double,G4double>(DBL_MAX,-DBL_MAX);
274 char* part = std::getenv(
"DICOM_CHANGE_MATERIAL_DENSITY" );
275 G4double densityDiff = -1.;
276 if( part ) densityDiff = G4UIcommand::ConvertToDouble(part);
277 std::map<G4int,G4double> densityDiffs;
278 if( densityDiff != -1. )
282 densityDiffs[ii] = densityDiff;
298 std::map< std::pair<G4Material*,G4int>,
matInfo* > newMateDens;
301 G4int ifxmin1, ifxmax1;
315 if( ix >= G4int(ifxmin1) && ix <= G4int(ifxmax1) ) {
319 mpite = densiMinMax.find( G4int(
fMateIDs[copyNo]) );
320 if( dens1 < (*mpite).second.first ) (*mpite).second.first = dens1;
321 if( dens1 > (*mpite).second.second ) (*mpite).second.second = dens1;
324 G4int mateID = G4int(
fMateIDs[copyNo]);
328 if( std::fabs(dens1 - mate->GetDensity()/g*cm3 ) < 1.e-9 )
continue;
333 G4String newMateName = mate->GetName();
335 G4int densityBin = 0;
339 if( densityDiff != -1.) {
340 densityBin = (G4int(dens1/densityDiffs[mateID]));
342 mate->GetName()+G4UIcommand::ConvertToString(densityBin);
346 std::pair<G4Material*,G4int> matdens(mate, densityBin );
348 auto mppite = newMateDens.find( matdens );
349 if( mppite != newMateDens.cend() ){
350 matInfo* mi = (*mppite).second;
360 mi->
fId = G4int(newMateDens.size()+1);
361 newMateDens[matdens] = mi;
384 for(
auto mppite= newMateDens.cbegin();
385 mppite != newMateDens.cend(); ++mppite )
387 G4double averdens = (*mppite).second->fSumdens/(*mppite).second->fNvoxels;
388 G4double saverdens = G4int(1000.001*averdens)/1000.;
389 G4cout <<
"GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
390 << averdens <<
" -> " << saverdens <<
" -> "
391 << G4int(1000*averdens) <<
" " << G4int(1000*averdens)/1000 <<
" "
392 << G4int(1000*averdens)/1000. << G4endl;
394 G4cout <<
"GmReadPhantomGeometry::ReadVoxelDensities AVER DENS "
395 << averdens <<
" -> " << saverdens <<
" -> "
396 << G4UIcommand::ConvertToString(saverdens) <<
" Nvoxels "
397 <<(*mppite).second->fNvoxels << G4endl;
398 G4String mateName = ((*mppite).first).first->GetName() +
"_" +
399 G4UIcommand::ConvertToString(saverdens);
401 (*mppite).first.first,
402 G4float(averdens), mateName ) );
443 G4bool bSkipEqualMaterials = 0;
444 G4int regStructureID = 2;
446 G4ThreeVector posCentreVoxels(0.,0.,0.);
459 phantom_logic->SetVisAttributes( visAtt );
462 G4double theSmartless = 0.5;
464 phantom_logic->SetSmartless( theSmartless );
483 if( OptimAxis ==
"kUndefined" ) {
487 }
else if( OptimAxis ==
"kXAxis" ) {
494 G4Exception(
"GmReadPhantomGeometry::ConstructPhantom",
496 FatalErrorInArgument,
497 G4String(
"Wrong argument to parameter \
498 GmReadPhantomGeometry:Phantom:OptimAxis: \
499 Only allowed 'kUndefined' or 'kXAxis', it is: "+OptimAxis).c_str());
502 phantom_phys->SetRegularStructureId(regStructureID);