Skip to content
Snippets Groups Projects
Commit 44a56feb authored by Johannes Junggeburth's avatar Johannes Junggeburth :dog2: Committed by Vakhtang Tsulaia
Browse files

GeoMaterialHelper - Minor clean

parent 75034a8a
Branches
Tags
1 merge request!404GeoMaterialHelper - Minor clean
Pipeline #9971236 passed
Stage: step-A
Show whitespace changes
Inline Side-by-side
GeoModelTools/GeoModelValidation/CMakeLists.txt
+ 3
− 1
View file @ 44a56feb
......@@ -8,7 +8,9 @@ file( GLOB HEADERS GeoModelValidation/*.h )
# Create the library.
add_library( GeoModelValidation SHARED ${HEADERS} ${SOURCES} )
target_link_libraries( GeoModelValidation PUBLIC Eigen3::Eigen GeoModelCore::GeoGenericFunctions GeoModelCore::GeoModelKernel GeoModelIO::GeoModelDBManager GeoModelIO::GeoModelWrite GeoModelIO::GeoModelRead )
target_link_libraries( GeoModelValidation PUBLIC Eigen3::Eigen GeoModelCore::GeoGenericFunctions GeoModelCore::GeoModelKernel
GeoModelCore::GeoModelHelpers GeoModelIO::GeoModelDBManager
GeoModelIO::GeoModelWrite GeoModelIO::GeoModelRead )
target_include_directories( GeoModelValidation PUBLIC
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}> )
......
GeoModelTools/GeoModelValidation/GeoModelValidation/GeoMaterialHelper.h
+ 11
− 12
View file @ 44a56feb
......@@ -10,9 +10,9 @@
#define GEOMATERIALHELPER_H
#include <GeoModelKernel/GeoDefinitions.h>
#include <GeoModelKernel/GeoMaterial.h>
#include <GeoModelKernel/GeoVPhysVol.h>
class GeoMaterial;
class GeoVPhysVol;
namespace GeoModelTools {
......@@ -24,30 +24,29 @@ namespace GeoModelTools {
@author sarka.todorova@cern.ch
*/
typedef std::pair< const GeoMaterial*, double > MaterialComponent;
using GeoMaterialPtr = GeoIntrusivePtr<const GeoMaterial>;
using MaterialComponent = std::pair<GeoMaterialPtr, double>;
class GeoMaterialHelper {
public:
/** Default constructor*/
GeoMaterialHelper()
{}
GeoMaterialHelper() = default;
/** Destructor*/
virtual ~GeoMaterialHelper(){}
virtual ~GeoMaterialHelper() = default;
/** Evaluate mass ( of a GeoModel tree/branch/volume ) */
float evaluateMass(const GeoVPhysVol* gv, bool inclusive = true) const;
double evaluateMass(PVConstLink gv, bool inclusive = true) const;
/** Collect and blend material : return blended material (the client takes the ownership) and effective volume */
MaterialComponent collectMaterial(const GeoVPhysVol* gv) const;
MaterialComponent collectMaterial(PVConstLink gv) const;
/** hardcoded dummy materials : TODO : find generic criterium ( density ? radiation length ? ) */
bool dummy_material(const GeoVPhysVol*) const;
bool dummy_material(PVConstLink gv) const;
private:
/** Internal recursive loop over material components */
void collectMaterialContent( const GeoVPhysVol* gv, std::vector< MaterialComponent >& materialContent ) const;
void collectMaterialContent( PVConstLink gv, std::vector< MaterialComponent >& materialContent) const;
};
......
GeoModelTools/GeoModelValidation/src/GeoMaterialHelper.cxx
+ 23
− 34
View file @ 44a56feb
/*
Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
///////////////////////////////////////////////////////////////////
......@@ -16,69 +16,66 @@
#include <string>
float GeoModelTools::GeoMaterialHelper::evaluateMass( const GeoVPhysVol* gv, bool inclusive) const {
double GeoModelTools::GeoMaterialHelper::evaluateMass( const PVConstLink gv, bool inclusive) const {
const GeoLogVol* lv = gv->getLogVol();
double weight{0.}, motherVolume{0.};
float weight = 0.;
double motherVolume = 0.;
// skip volume calculation for dummy material configuration
if ( !this->dummy_material(gv) ) motherVolume = lv->getShape()->volume();
if ( !dummy_material(gv)) motherVolume = lv->getShape()->volume();
// subtract children volumes
//GeoVolumeVec_t vols = geoGetVolumes(gv);
//for (const GeoVolumeVec_t::value_type& p : vols) {
unsigned int nc = gv->getNChildVols();
for (unsigned int ic = 0; ic < nc; ic++) {
const GeoVPhysVol* cv = &(*(gv->getChildVol(ic)));
const PVConstLink cv = gv->getChildVol(ic);
double childVolume = cv->getLogVol()->getShape()->volume() ;
// skip volume calculation for dummy material configuration
if (!this->dummy_material(gv) ) motherVolume -= childVolume; // child volumes fully confined by construction
if (!dummy_material(gv) ) motherVolume -= childVolume; // child volumes fully confined by construction
if (inclusive) weight += evaluateMass(cv, inclusive);
}
// skip mass calculation for dummy material configuration
if (!this->dummy_material(gv)) weight += motherVolume * lv->getMaterial()->getDensity();
if (!dummy_material(gv)) weight += motherVolume * lv->getMaterial()->getDensity();
return weight;
}
GeoModelTools::MaterialComponent GeoModelTools::GeoMaterialHelper::collectMaterial( const GeoVPhysVol* gv) const {
GeoModelTools::MaterialComponent GeoModelTools::GeoMaterialHelper::collectMaterial(PVConstLink gv) const {
// iterative loop over substructure
std::vector< GeoModelTools::MaterialComponent > materialContent ;
collectMaterialContent( gv, materialContent );
// material blending
double volume = 0.;
double mass = 0.;
double volume{0.}, mass{0.};
for ( auto matInput : materialContent ) {
volume += matInput.second;
mass += matInput.first->getDensity() * matInput.second;
}
GeoMaterial* gmMatTmp=new GeoMaterial("MaterialBlend", mass / volume);
auto gmMatTmp= make_intrusive<GeoMaterial>("MaterialBlend", mass / volume);
for ( auto matInput : materialContent ) {
gmMatTmp->add( matInput.first, matInput.first->getDensity() * matInput.second / mass );
}
gmMatTmp->lock();
return std::pair< const GeoMaterial*, double > ( gmMatTmp, volume );
return std::make_pair( gmMatTmp, volume );
}
void GeoModelTools::GeoMaterialHelper::collectMaterialContent( const GeoVPhysVol* gv, std::vector< GeoModelTools::MaterialComponent >& materialContent ) const {
void GeoModelTools::GeoMaterialHelper::collectMaterialContent(const PVConstLink gv,
std::vector< GeoModelTools::MaterialComponent >& materialContent ) const {
const GeoLogVol* lv = gv->getLogVol();
double motherVolume = 0.;
double motherVolume{0.};
// skip volume calculation for dummy material configuration
if (!this->dummy_material(gv) ) motherVolume = lv->getShape()->volume();
if (!dummy_material(gv) ) motherVolume = lv->getShape()->volume();
// subtract children volumes from mother volume & collect their material content
//GeoVolumeVec_t vols = geoGetVolumes(gv);
......@@ -88,28 +85,20 @@ void GeoModelTools::GeoMaterialHelper::collectMaterialContent( const GeoVPhysVol
const GeoVPhysVol* cv = &(*(gv->getChildVol(ic)));
double childVolume = cv->getLogVol()->getShape()->volume() ;
// skip for dummy material configuration
if (!this->dummy_material(gv) ) motherVolume -= childVolume; // child volumes fully confined by construction
if (!dummy_material(gv) ) motherVolume -= childVolume; // child volumes fully confined by construction
collectMaterialContent( cv, materialContent );
}
// skip for dummy material configuration
if (!this->dummy_material(gv) ) materialContent.push_back( std::pair<const GeoMaterial*, double> ( lv->getMaterial(), motherVolume ) );
if (!dummy_material(gv) ) materialContent.emplace_back(lv->getMaterial(), motherVolume);
return;
}
bool GeoModelTools::GeoMaterialHelper::dummy_material(const GeoVPhysVol* vol) const {
bool dummyMat = false;
std::string mat = vol->getLogVol()->getMaterial()->getName();
if ( mat=="special::Ether" ) dummyMat = true;
if ( mat=="WorldLog:Air" ) dummyMat = true;
if ( mat=="std::Air" ) dummyMat = true;
if ( mat=="Air" ) dummyMat = true;
return dummyMat;
bool GeoModelTools::GeoMaterialHelper::dummy_material(const PVConstLink vol) const {
static const std::array<std::string, 4> dummyMatNames{"special::Ether","WorldLog:Air", "std::Air" , "Air"};
const std::string& matName = vol->getLogVol()->getMaterial()->getName();
return std::find(dummyMatNames.begin(), dummyMatNames.end(),matName ) != dummyMatNames.end();
}
GeoModelTools/GeoModelValidation/src/GeoPhysVolHelper.cxx
+ 10
− 175
View file @ 44a56feb
......@@ -34,6 +34,9 @@
#include "GeoModelKernel/GeoPgon.h"
#include "GeoModelKernel/GeoPara.h"
#include "GeoModelHelpers/GeoShapeSorter.h"
#include "GeoModelHelpers/TransformSorter.h"
#include "GeoModelKernel/Units.h"
#define SYSTEM_OF_UNITS GeoModelKernelUnits
......@@ -539,172 +542,9 @@ void GeoModelTools::GeoPhysVolHelper::printChildren(const GeoVPhysVol* pv, int l
}
}
bool GeoModelTools::GeoPhysVolHelper::compareShapes( const GeoShape* sh1, const GeoShape* sh2, float tol ) const {
if (sh1->type() != sh2->type()) return false;
if ( sh1->type()=="Pgon") {
const GeoPgon* pgon1 = dynamic_cast<const GeoPgon*>(sh1);
const GeoPgon* pgon2 = dynamic_cast<const GeoPgon*>(sh2);
if (!pgon1 || !pgon2) return false;
if (pgon1->getNPlanes() != pgon2->getNPlanes()) return false;
if (pgon1->getNSides() != pgon2->getNSides()) return false;
if (fabs(pgon1->getSPhi() - pgon2->getSPhi()) > tol) return false;
if (fabs(pgon1->getDPhi() - pgon2->getDPhi()) > tol) return false;
return true;
} else if (sh1->type()=="Trd") {
const GeoTrd* trd1 = dynamic_cast<const GeoTrd*> (sh1);
const GeoTrd* trd2 = dynamic_cast<const GeoTrd*> (sh2);
if ((trd1->getXHalfLength1() - trd2->getXHalfLength1())>tol) return false;
if ((trd1->getXHalfLength2() - trd2->getXHalfLength2())>tol) return false;
if ((trd1->getYHalfLength1() - trd2->getYHalfLength1())>tol) return false;
if ((trd1->getYHalfLength2() - trd2->getYHalfLength2())>tol) return false;
if ((trd1->getZHalfLength() - trd2->getZHalfLength())>tol) return false;
return true;
} else if ( sh1->type()=="Box") {
const GeoBox* box1 = dynamic_cast<const GeoBox*> (sh1);
const GeoBox* box2 = dynamic_cast<const GeoBox*> (sh2);
if (fabs(box1->getXHalfLength() - box2->getXHalfLength()) > tol) return false;
if (fabs(box1->getYHalfLength() - box2->getYHalfLength()) > tol) return false;
if (fabs(box1->getZHalfLength() - box2->getZHalfLength()) > tol) return false;
return true;
} else if ( sh1->type() == "Tube" ) {
const GeoTube* tube1=dynamic_cast<const GeoTube*> (sh1);
const GeoTube* tube2=dynamic_cast<const GeoTube*> (sh2);
if ( fabs(tube1->getRMin() - tube2->getRMin()) > tol) return false;
if ( fabs(tube1->getRMax() - tube2->getRMax()) > tol) return false;
if ( fabs(tube1->getZHalfLength() - tube2->getZHalfLength()) > tol) return false;
return true;
} else if ( sh1->type() == "Tubs" ) {
const GeoTubs* tubs1=dynamic_cast<const GeoTubs*> (sh1);
const GeoTubs* tubs2=dynamic_cast<const GeoTubs*> (sh2);
if ( fabs(tubs1->getRMin() - tubs2->getRMin()) >tol) return false;
if ( fabs(tubs1->getRMax() - tubs2->getRMax()) > tol) return false;
if ( fabs(tubs1->getZHalfLength() - tubs2->getZHalfLength()) > tol) return false;
if ( fabs(tubs1->getSPhi() - tubs2->getSPhi()) > tol) return false;
if ( fabs(tubs1->getDPhi() - tubs2->getDPhi()) > tol) return false;
return true;
} else if ( sh1->type() == "Cons" ) {
const GeoCons* cons1=dynamic_cast<const GeoCons*> (sh1);
const GeoCons* cons2=dynamic_cast<const GeoCons*> (sh2);
if ( fabs(cons1->getRMin1() - cons2->getRMin1()) > tol) return false;
if ( fabs(cons1->getRMin2() - cons2->getRMin2()) >tol) return false;
if ( fabs(cons1->getRMax1() - cons2->getRMax1()) > tol) return false;
if ( fabs(cons1->getRMax2() - cons2->getRMax2()) > tol) return false;
if ( fabs(cons1->getDZ() - cons2->getDZ()) > tol) return false;
if ( fabs(cons1->getSPhi() - cons2->getSPhi()) > tol) return false;
if ( fabs(cons1->getDPhi() - cons2->getDPhi()) > tol) return false;
return true;
} else if ( sh1->type()=="SimplePolygonBrep") {
const GeoSimplePolygonBrep* spb1 = dynamic_cast<const GeoSimplePolygonBrep*> (sh1);
const GeoSimplePolygonBrep* spb2 = dynamic_cast<const GeoSimplePolygonBrep*> (sh2);
if (!spb1 || !spb2) return false;
unsigned int nv1 = spb1->getNVertices();
unsigned int nv2 = spb2->getNVertices();
if (nv1 != nv2) return false;
if (fabs(spb1->getDZ() - spb2->getDZ()) > tol) return false;
for (unsigned int iv = 0; iv < nv1; iv++) {
if (fabs(spb1->getXVertex(iv) - spb2->getXVertex(iv)) > tol) return false;
if (fabs(spb1->getYVertex(iv) - spb2->getYVertex(iv)) > tol) return false;
}
return true;
} else if ( sh1->type()=="Pcon") {
const GeoPcon* pc1 = dynamic_cast<const GeoPcon*> (sh1);
const GeoPcon* pc2 = dynamic_cast<const GeoPcon*> (sh2);
if (!pc1 || !pc2) return false;
if ( fabs(pc1->getSPhi() - pc2->getSPhi()) > tol) return false;
if ( fabs(pc1->getDPhi() - pc2->getDPhi()) > tol) return false;
unsigned int nv1 = pc1->getNPlanes();
unsigned int nv2 = pc2->getNPlanes();
if (nv1 != nv2) return false;
for (unsigned int iv = 0; iv < nv1; iv++) {
if (fabs(pc1->getZPlane(iv) - pc2->getZPlane(iv)) > tol) return false;
if (fabs(pc1->getRMinPlane(iv) - pc2->getRMinPlane(iv)) > tol) return false;
if (fabs(pc1->getRMaxPlane(iv) - pc2->getRMaxPlane(iv)) > tol) return false;
}
return true;
} else if ( sh1->type()=="Subtraction") {
const GeoShapeSubtraction* sub1 = dynamic_cast<const GeoShapeSubtraction*> (sh1);
const GeoShapeSubtraction* sub2 = dynamic_cast<const GeoShapeSubtraction*> (sh2);
if (!sub1 || !sub2) return false;
if (!compareShapes(sub1->getOpA(),sub2->getOpA(),tol)) return false;
if (!compareShapes(sub1->getOpB(),sub2->getOpB(),tol)) return false;
return true;
} else if ( sh1->type()=="Union") {
const GeoShapeUnion* sub1 = dynamic_cast<const GeoShapeUnion*> (sh1);
const GeoShapeUnion* sub2 = dynamic_cast<const GeoShapeUnion*> (sh2);
if (!sub1 || !sub2) return false;
if (!compareShapes(sub1->getOpA(),sub2->getOpA(),tol)) return false;
if (!compareShapes(sub1->getOpB(),sub2->getOpB(),tol)) return false;
return true;
} else if ( sh1->type()=="Shift") {
const GeoShapeShift* shift1 = dynamic_cast<const GeoShapeShift*> (sh1);
const GeoShapeShift* shift2 = dynamic_cast<const GeoShapeShift*> (sh2);
if (!shift1 || !shift2) return false;
if (!compareShapes(shift1->getOp(),shift2->getOp(),tol)) return false;
const GeoTrf::Transform3D& transf1 = shift1->getX();
const GeoTrf::Transform3D& transf2 = shift2->getX();
if ((transf1.translation()-transf2.translation()).norm()>tol) return false;
if (!identity_check(transf1.rotation()*transf2.rotation().inverse(), tol)) return false;
return true;
}
std::cout <<"unknown shape to compare:"<<sh1->type()<< std::endl;
return false;
}
bool GeoModelTools::GeoPhysVolHelper::identity_check(GeoTrf::RotationMatrix3D rotation, float tol) const {
if (fabs(rotation(0,1))>tol) return false;
if (fabs(rotation(0,2))>tol) return false;
if (fabs(rotation(1,2))>tol) return false;
return true;
bool GeoModelTools::GeoPhysVolHelper::compareShapes( const GeoShape* sh1, const GeoShape* sh2, float /*tol*/ ) const {
GeoShapeSorter sorter{};
return sorter.compare(sh1, sh2) != 0;
}
void GeoModelTools::GeoPhysVolHelper::decodeShape(const GeoShape* sh) const {
......@@ -763,15 +603,10 @@ void GeoModelTools::GeoPhysVolHelper::decodeShape(const GeoShape* sh) const {
}
int GeoModelTools::GeoPhysVolHelper::compareTransforms(GeoTrf::Transform3D tr_test, GeoTrf::Transform3D tr_ref, float tolerance) const {
int trcode = 0; // assuming identical
if ((tr_test.translation()-tr_ref.translation()).norm()>tolerance) trcode = 6;
if (!identity_check(tr_test.rotation()*tr_ref.rotation().inverse(), tolerance)) trcode = 7;
return trcode;
int GeoModelTools::GeoPhysVolHelper::compareTransforms(GeoTrf::Transform3D tr_test, GeoTrf::Transform3D tr_ref,
float /*tolerance*/) const {
GeoTrf::TransformSorter sorter{};
return sorter.compare(tr_test, tr_ref);
}
void GeoModelTools::GeoPhysVolHelper::printTranslationDiff(GeoTrf::Transform3D tr_test, GeoTrf::Transform3D tr_ref, float tolerance) const {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment