Extend poly edges

GeoModelCore/GeoModelHelpers/src/GeoShapeUtils.cxx

@@ -28,6 +28,10 @@
 
 #include <vector>
 
+namespace{
+    constexpr unsigned tubeApproxEdges = 8;
+}
+
 std::pair<const GeoShape* , const GeoShape*> getOps(const GeoShape* composed) {
     if (composed->typeID() == GeoShapeUnion::getClassTypeID()) {
         const GeoShapeUnion* unionShape = dynamic_cast<const GeoShapeUnion*>(composed);
@@ -48,7 +52,7 @@ std::pair<const GeoShape* , const GeoShape*> getOps(const GeoShape* composed) {
 unsigned int countComposedShapes(const GeoShape* shape) {
     std::pair<const GeoShape*, const GeoShape*> ops = getOps(shape);
     return 1 + (ops.first ? countComposedShapes(ops.first) : 0) +
-                (ops.second ? countComposedShapes(ops.second) : 0);
+               (ops.second ? countComposedShapes(ops.second) : 0);
 }
 
 GeoIntrusivePtr<const GeoShape> compressShift(const GeoShape* shift) {    
@@ -229,6 +233,46 @@ std::vector<GeoTrf::Vector3D> getPolyShapeEdges(const GeoShape* shape,
             
             }
         }
+    } else if (shape->typeID() == GeoTube::getClassTypeID()) {
+        constexpr double stepLength = 2.*M_PI / tubeApproxEdges;
+        const GeoTube* tube = static_cast<const GeoTube*>(shape);
+        for (const double z :{-1., 1.}){
+            for (unsigned int e = 0; e < tubeApproxEdges; ++e) {
+                edgePoints.emplace_back(tube->getRMax()*std::cos(stepLength*e),
+                                        tube->getRMax()*std::sin(stepLength*e),
+                                        z*tube->getZHalfLength());
+            }
+        }
+    } else if (shape->typeID() == GeoTubs::getClassTypeID()) {
+        const GeoTubs* tubs = static_cast<const GeoTubs*>(shape);
+        const double stepSize = tubs->getDPhi() / tubeApproxEdges;
+        for (const double z : {-1., 1.}) {
+            for (unsigned e = 0; e <= tubeApproxEdges; ++e) {
+                edgePoints.emplace_back(tubs->getRMax()*std::cos(tubs->getSPhi() + stepSize*e), 
+                                        tubs->getRMax()*std::sin(tubs->getSPhi() + stepSize*e), 
+                                        z*tubs->getZHalfLength());
+            }
+        }
+    } else if (shape->typeID() == GeoEllipticalTube::getClassTypeID()) {
+        const GeoEllipticalTube* tube = static_cast<const GeoEllipticalTube*>(shape);
+        constexpr double stepSize = 2.*M_PI / tubeApproxEdges;
+        for (const double z : {-1.,1.}) {
+            for (unsigned e = 0; e<tubeApproxEdges; ++e){
+                edgePoints.emplace_back(tube->getXHalfLength()*std::cos(stepSize*e),
+                                        tube->getYHalfLength()*std::sin(stepSize*e),
+                                        z* tube->getZHalfLength());
+            }
+        }
+    } else if (shape->typeID() == GeoTorus::getClassTypeID()) {
+        constexpr double stepSize = 2.*M_PI / tubeApproxEdges;
+        const GeoTorus* torus = static_cast<const GeoTorus*>(shape);
+        for (unsigned int donut =0; donut < tubeApproxEdges; ++donut) {
+            for (unsigned int slice =0; slice < tubeApproxEdges; ++slice) {
+                edgePoints.emplace_back((torus->getRTor() + torus->getRMax()* std::cos(slice *stepSize))*std::cos(donut*stepSize),
+                                        (torus->getRTor() + torus->getRMax()* std::cos(slice *stepSize))*std::sin(donut*stepSize),
+                                         torus->getRMax()* std::sin(slice*stepSize));
+            }
+        }
     } else {
         THROW_EXCEPTION("The shape "<<shape->type()<<" is not supported. Please add it to the list");
     }

GeoModelCore/GeoModelHelpers/src/MaterialManager.cxx

@@ -33,15 +33,14 @@ namespace {
 const GeoMaterial* MaterialManager::MaterialFactory::get() const { return m_material; }
 MaterialManager::MaterialFactory::~MaterialFactory() {
     if (!m_material) return;
-    // std::cout<<"MaterialFactory() -- finalize " <<m_material->getName()<<std::endl;
+
     if (m_components.empty()) {
       m_material->lock();
       return;
     }
     const double inv_totalFraction = 1. / m_totFraction;
     for(const auto& [element, fraction] : m_components) {
-      m_material->add(element, 
-                      fraction * element->getA() * inv_totalFraction);
+      m_material->add(element, fraction * element->getA() * inv_totalFraction);
     }
     m_material->lock();
 }
@@ -56,8 +55,7 @@ void MaterialManager::MaterialFactory::addComponent(const ConstMaterialPtr& mat,
         addComponent(elePtr, mat->getFraction(ele) * fraction);
     }
 }
-void MaterialManager::MaterialFactory::addComponent(const ConstElementPtr& ele, 
-                                                    double fraction) {
+void MaterialManager::MaterialFactory::addComponent(const ConstElementPtr& ele, double fraction) {
     m_components.emplace_back(std::make_pair(ele, fraction));
     m_totFraction += fraction;
 }
@@ -135,7 +133,7 @@ void MaterialManager::printAll() const {
 }
 
 void MaterialManager::addElement(const std::string &name, const std::string &symbol, double z, double a) {
-    GeoIntrusivePtr<GeoElement> newElement{new GeoElement(name,symbol,z,a*atomicDensity)};
+    GeoIntrusivePtr<GeoElement> newElement{make_intrusive<GeoElement>(name,symbol,z,a*atomicDensity)};
     addElement(newElement);
 }
 
@@ -158,7 +156,7 @@ void MaterialManager::setMaterialNamespace(const std::string &name) {
 void MaterialManager::lockMaterial() { m_factory.reset(); }
 
 void MaterialManager::addMaterial(const std::string &name, double density) {
-    MaterialPtr newMat{new GeoMaterial(name, density * volDensity)};
+    MaterialPtr newMat{make_intrusive<GeoMaterial>(name, density * volDensity)};
     addMaterial(newMat);
 }
 
@@ -200,14 +198,14 @@ void MaterialManager::addMatComponent(const std::string &name, double fraction)
 }
 void MaterialManager::buildSpecialMaterials() {
   // Ether  
-  GeoIntrusivePtr<GeoElement> ethElement{new GeoElement("Ether","ET",500.0,0.0)};
+  GeoIntrusivePtr<GeoElement> ethElement{make_intrusive<GeoElement>("Ether","ET",500.0,0.0)};
   m_elements.insert(std::make_pair("Ether",ethElement));
-  GeoIntrusivePtr<GeoMaterial> ether{new GeoMaterial("special::Ether",0.0)};
+  GeoIntrusivePtr<GeoMaterial> ether{make_intrusive<GeoMaterial>("special::Ether",0.0)};
   ether->add(ethElement,1.);
   ether->lock();
   m_materials.insert(std::make_pair("special::Ether", ether));
   // HyperUranium
-  GeoIntrusivePtr<GeoMaterial> hu{new GeoMaterial("special::HyperUranium",0.0)};
+  GeoIntrusivePtr<GeoMaterial> hu{make_intrusive<GeoMaterial>("special::HyperUranium",0.0)};
   hu->add(ethElement,1.);
   hu->lock();
   m_materials.insert(std::make_pair("special::HyperUranium", hu));