/* Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration */ #include "ActsGeometry/ActsCaloTrackingVolumeBuilder.h" #include "ActsInterop/Logger.h" #include "StoreGate/ReadHandle.h" #include "CaloDetDescr/CaloDetDescrManager.h" #include "CaloDetDescr/CaloDetDescrElement.h" #include "CaloDetDescrUtils/CaloDetDescrBuilder.h" #include "Acts/Geometry/TrackingVolume.hpp" #include "Acts/Geometry/VolumeBounds.hpp" #include "Acts/Geometry/GlueVolumesDescriptor.hpp" // ACTS #include "Acts/Geometry/GenericCuboidVolumeBounds.hpp" #include "Acts/Geometry/CutoutCylinderVolumeBounds.hpp" #include "Acts/Geometry/CylinderVolumeBounds.hpp" //#include "Acts/Utilities/IVisualization.hpp" #include "Acts/Utilities/Helpers.hpp" #include "Acts/Geometry/TrackingVolumeArrayCreator.hpp" #include "Acts/Utilities/BinnedArrayXD.hpp" #include "Acts/Geometry/CylinderVolumeHelper.hpp" #include <fstream> using Box = Acts::Volume::BoundingBox; // shortcut using CVBBV = Acts::CylinderVolumeBounds::BoundValues; using CCVBBV = Acts::CutoutCylinderVolumeBounds::BoundValues; ActsCaloTrackingVolumeBuilder::ActsCaloTrackingVolumeBuilder(const std::string& type, const std::string& name, const IInterface* parent) : base_class(type, name, parent) { } StatusCode ActsCaloTrackingVolumeBuilder::initialize() { m_caloMgr = detStore()->tryConstRetrieve<CaloDetDescrManager>(caloMgrStaticKey); if(!m_caloMgr) { std::unique_ptr<CaloDetDescrManager> caloMgrPtr = buildCaloDetDescrNoAlign(serviceLocator() , Athena::getMessageSvc()); ATH_CHECK(detStore()->record(std::move(caloMgrPtr), caloMgrStaticKey)); ATH_CHECK(detStore()->retrieve(m_caloMgr, caloMgrStaticKey)); } return StatusCode::SUCCESS; } std::shared_ptr<Acts::TrackingVolume> ActsCaloTrackingVolumeBuilder::trackingVolume( const Acts::GeometryContext& gctx, std::shared_ptr<const Acts::TrackingVolume> insideVolume, std::shared_ptr<const Acts::VolumeBounds> /*outsideBounds*/) const { // generate the calo cell volume description std::vector<std::unique_ptr<Acts::Volume>> cells; cells = cellFactory(); ATH_MSG_VERBOSE("Collected " << cells.size() << " calo cells"); // we need to turn the cells into boundary boxes std::vector<std::unique_ptr<Box>> boxStore; std::vector<Box*> prims; for (const auto& cell : cells) { boxStore.push_back( std::make_unique<Box>(cell->boundingBox({0.1, 0.1, 0.1}))); prims.push_back(boxStore.back().get()); } ATH_MSG_VERBOSE("Generated Bounding Boxes"); ATH_MSG_VERBOSE("Figure out dimensions of wrapping volume"); std::shared_ptr<Acts::CutoutCylinderVolumeBounds> caloVolBounds = makeCaloVolumeBounds(boxStore, insideVolume); // build a BVH octree for the bounding boxes // but only AFTER we've built the calo volume bounds // Box* top; // top = Acts::make_octree(boxStore, prims, 1, 0.1); // Create Tracking Volume that coutains the Calo // This needs to own the Abstract Volumes (cells), but we // need to up-cast them to Volume, since that's what TrackingVolume can own std::vector<std::unique_ptr<const Acts::Volume>> cellVols; cellVols.reserve(cells.size()); for(auto& cell : cells) { std::unique_ptr<const Acts::Volume> up; // release, up-cast, then immediately pack again up = std::unique_ptr<const Acts::Volume>(dynamic_cast<const Acts::Volume*>(cell.release())); cellVols.push_back(std::move(up)); } // This was removed in https://github.com/acts-project/acts/pull/3029 // To be reimplemented using new geometry model instead of explicit TrackingVolume content (void)gctx; // suppress compiler warning throw std::runtime_error{"Calo building for ACTS currently disabled"}; /***** TODO START ***** std::shared_ptr<Acts::TrackingVolume> calo; // = Acts::TrackingVolume::create(Acts::Transform3::Identity(), // caloVolBounds, // std::move(boxStore), // std::move(cellVols), // top, // nullptr, // no material for now // "Calo"); // We need to interglue all the volumes together std::shared_ptr<Acts::TrackingVolume> mutInsideVolume = std::const_pointer_cast<Acts::TrackingVolume>(insideVolume); auto idBounds = dynamic_cast<const Acts::CylinderVolumeBounds*>(&insideVolume->volumeBounds()); if (idBounds == nullptr) { // protection against nullptr ATH_MSG_ERROR("Unable to dynamic cast volume bounds to Acts::CylinderVolumeBounds"); throw std::runtime_error("Error casting to CylinderVolumeBounds"); } // we want gap volumes at pos and neg xy face, and at outer cyl cover // which will include the solenoid area auto trackingVolumeArrayCreator = std::make_shared<const Acts::TrackingVolumeArrayCreator>( Acts::TrackingVolumeArrayCreator::Config{}, makeActsAthenaLogger(this, std::string("TrkVolArrCrtr"), std::string("ActsTGSvc"))); Acts::CylinderVolumeHelper::Config cvhCfg; cvhCfg.trackingVolumeArrayCreator = trackingVolumeArrayCreator; Acts::CylinderVolumeHelper cvh(cvhCfg, makeActsAthenaLogger(this, std::string("ACaloTrkVB"), std::string("CylVolHlp"))); std::vector<double> lPos = {}; std::vector<std::shared_ptr<Acts::TrackingVolume>> noVolumes; ATH_MSG_VERBOSE("Creating gap volume to extend ID"); // positive xy gap auto idGapPosXY = cvh.createGapTrackingVolume(gctx, noVolumes, nullptr, idBounds->get(CVBBV::eMinR), idBounds->get(CVBBV::eMaxR), idBounds->get(CVBBV::eHalfLengthZ), caloVolBounds->get(CCVBBV::eHalfLengthZcutout), lPos, false, "ID::PositiveGap" ); // negative xy gap auto idGapNegXY = cvh.createGapTrackingVolume(gctx, noVolumes, nullptr, idBounds->get(CVBBV::eMinR), idBounds->get(CVBBV::eMaxR), -caloVolBounds->get(CCVBBV::eHalfLengthZcutout), -idBounds->get(CVBBV::eHalfLengthZ), lPos, false, "ID::NegativeGap" ); // outer cover gap auto idGapCylOuter = cvh.createGapTrackingVolume(gctx, noVolumes, nullptr, idBounds->get(CVBBV::eMaxR), caloVolBounds->get(CCVBBV::eMedR), -caloVolBounds->get(CCVBBV::eHalfLengthZcutout), +caloVolBounds->get(CCVBBV::eHalfLengthZcutout), lPos, false, "ID::CylOutGap" ); ATH_MSG_VERBOSE("Create container volume to contain ID and gap volumes"); auto idContainerZ = cvh.createContainerTrackingVolume(gctx, {idGapNegXY, mutInsideVolume, idGapPosXY}); auto idContainer = cvh.createContainerTrackingVolume(gctx, {idContainerZ, idGapCylOuter}); ATH_MSG_VERBOSE("Begin volume glueing"); const Acts::GlueVolumesDescriptor& gvd = idContainer->glueVolumesDescriptor(); // let's see what the GVD says is on the inner cover of the ID const auto& tVolArrPos = gvd.glueVolumes(Acts::positiveFaceXY); const auto& tVolArrNeg = gvd.glueVolumes(Acts::negativeFaceXY); std::cout << "POSITIVE: " << std::endl; for(const auto& subvol : tVolArrPos->arrayObjects()) { std::cout << subvol->volumeName() << std::endl; std::cout << *subvol << std::endl; } std::cout << "NEGATIVE: " << std::endl; for(const auto& subvol : tVolArrNeg->arrayObjects()) { std::cout << subvol->volumeName() << std::endl; std::cout << *subvol << std::endl; } using BoundarySurface = Acts::BoundarySurfaceT<Acts::TrackingVolume>; // Glue outer radial cover of ID to inner cover of Calo cutout auto idOutVolArray = gvd.glueVolumes(Acts::tubeOuterCover); // Attach that volume array to the calo inner cover ATH_MSG_VERBOSE("Glueing " << calo->volumeName() << " inner cover to " << idOutVolArray->arrayObjects().size() << " volumes"); std::const_pointer_cast<BoundarySurface>(calo->boundarySurfaces().at(Acts::tubeInnerCover)) ->attachVolumeArray(idOutVolArray, Acts::Direction::Backward); // Loop through the array and attach their boundary surfaces to the calo for(const auto& idVol : idOutVolArray->arrayObjects()){ ATH_MSG_VERBOSE("Glueing outer cover of " << idVol->volumeName() << " to inner cover of " << calo->volumeName()); std::const_pointer_cast<BoundarySurface>(idVol->boundarySurfaces().at(Acts::tubeOuterCover)) ->attachVolume(calo.get(), Acts::Direction::Forward); } // Glue positive XY face of ID to inner positive XY face of Calo. // ID has multiple, Calo has only one auto idPosXYVolArray = gvd.glueVolumes(Acts::positiveFaceXY); ATH_MSG_VERBOSE("Glueing " << calo->volumeName() << " positive inner cutout disc to " << idPosXYVolArray->arrayObjects().size() << " volumes"); std::const_pointer_cast<BoundarySurface>(calo->boundarySurfaces().at(Acts::index5)) ->attachVolumeArray(idPosXYVolArray, Acts::Direction::Backward); // Other way round, attach ID volumes to calo for(const auto& idVol : idPosXYVolArray->arrayObjects()){ ATH_MSG_VERBOSE("Glueing positive XY face of " << idVol->volumeName() << " to positive inner coutout disc of " << calo->volumeName()); std::const_pointer_cast<BoundarySurface>(idVol->boundarySurfaces().at(Acts::positiveFaceXY)) ->attachVolume(calo.get(), Acts::Direction::Forward); } // Glue negative XY face of ID to inner negative XY face of Calo. // ID has multiple, Calo has only one auto idNegXYVolArray = gvd.glueVolumes(Acts::negativeFaceXY); ATH_MSG_VERBOSE("Glueing " << calo->volumeName() << " negative inner cutout disc to " << idNegXYVolArray->arrayObjects().size() << " volumes"); std::const_pointer_cast<BoundarySurface>(calo->boundarySurfaces().at(Acts::index4)) ->attachVolumeArray(idNegXYVolArray, Acts::Direction::Forward); // Other way round, attach ID volumes to calo for(const auto& idVol : idNegXYVolArray->arrayObjects()){ ATH_MSG_VERBOSE("Glueing negative XY face of " << idVol->volumeName() << " to negative inner coutout disc of " << calo->volumeName()); std::const_pointer_cast<BoundarySurface>(idVol->boundarySurfaces().at(Acts::negativeFaceXY)) ->attachVolume(calo.get(), Acts::Direction::Backward); } // For navigational purposes we need to create three pseudo container cylinders. // One for the middle, which includes the central part of the Calo and the ID, and // two that fit cleanly at the +- XY face that then only include the Calo // Construct track vol array for use in positive and negative pseudocontainer. // This will only contain the calo double caloRMin = caloVolBounds->get(CCVBBV::eMinR); double caloRMed = caloVolBounds->get(CCVBBV::eMedR); double caloRMax = caloVolBounds->get(CCVBBV::eMaxR); double caloDZ1 = caloVolBounds->get(CCVBBV::eHalfLengthZ); double caloDZ2 = caloVolBounds->get(CCVBBV::eHalfLengthZcutout); Acts::Vector3 caloChokeRPos = {caloRMin + (caloRMax - caloRMin)/2., 0, 0}; std::vector<Acts::TrackingVolumeOrderPosition> tVolOrdPosNeg; tVolOrdPosNeg.push_back(std::make_pair(calo, caloChokeRPos)); std::vector<float> posNegBoundaries = {float(caloRMin), float(caloRMax)}; auto binUtilityPosNeg = std::make_unique<const Acts::BinUtility>(posNegBoundaries, Acts::open, Acts::BinningValue::binR); auto tVolArrPosNeg = std::make_shared<const Acts::BinnedArrayXD<Acts::TrackingVolumePtr>>( tVolOrdPosNeg, std::move(binUtilityPosNeg)); double chokeZOffset = caloDZ2 + (caloDZ1 - caloDZ2)/2.; Acts::Transform3 posTrf(Acts::Translation3(Acts::Vector3::UnitZ() * chokeZOffset)); Acts::Transform3 negTrf(Acts::Translation3(Acts::Vector3::UnitZ()* -1 *chokeZOffset)); auto posNegCylBounds = std::make_shared<Acts::CylinderVolumeBounds>( caloRMin, caloRMax, (caloDZ1 - caloDZ2) / 2.); // they share the same bounds and tvol array auto posContainer = std::make_shared<Acts::TrackingVolume>( posTrf, posNegCylBounds, nullptr, nullptr, tVolArrPosNeg, Acts::MutableTrackingVolumeVector{}); ATH_MSG_VERBOSE("Built positive container " << *posContainer); ATH_MSG_VERBOSE(" - containing: " << calo->volumeName()); auto negContainer = std::make_shared<Acts::TrackingVolume>( negTrf, posNegCylBounds, nullptr, nullptr, tVolArrPosNeg, Acts::MutableTrackingVolumeVector{}); ATH_MSG_VERBOSE("Built negative container " << *negContainer); ATH_MSG_VERBOSE(" - containing: " << calo->volumeName()); // now build the central pseudocontainer std::vector<Acts::TrackingVolumeOrderPosition> tVolOrderedCtr; tVolOrderedCtr.push_back(std::make_pair(idContainer, Acts::Vector3(caloRMed / 2., 0, 0))); tVolOrderedCtr.push_back(std::make_pair(calo, Acts::Vector3(caloRMed + (caloRMax- caloRMed) / 2., 0, 0))); std::vector<float> ctrBoundaries = {0, float(caloRMed), float(caloRMax)}; auto binUtilityCtr = std::make_unique<const Acts::BinUtility>( ctrBoundaries, Acts::open, Acts::BinningValue::binR); auto tVolArrCtr = std::make_shared<const Acts::BinnedArrayXD<Acts::TrackingVolumePtr>>( tVolOrderedCtr, std::move(binUtilityCtr)); auto ctrContainer = std::make_shared<Acts::TrackingVolume>(Acts::Transform3::Identity(), std::make_shared<Acts::CylinderVolumeBounds>( caloRMin, caloRMax, caloDZ2), nullptr, nullptr, tVolArrCtr, Acts::MutableTrackingVolumeVector{} ); ATH_MSG_VERBOSE("Built central container " << *ctrContainer); ATH_MSG_VERBOSE("- containing: " << idContainer->volumeName() << ", " << calo->volumeName()); // and now combine those together into another one Acts::TrackingVolumeArrayCreator tvac{Acts::TrackingVolumeArrayCreator::Config{}}; auto mainContainer = std::make_shared<Acts::TrackingVolume>(Acts::Transform3::Identity(), std::make_shared<Acts::CylinderVolumeBounds>( caloRMin, caloRMax, caloDZ1), nullptr, nullptr, tvac.trackingVolumeArray(gctx, {negContainer, ctrContainer, posContainer}, Acts::BinningValue::binZ), Acts::MutableTrackingVolumeVector{} ); ATH_MSG_VERBOSE("Built main container: " << *mainContainer); return mainContainer; ***** TODO END *****/ } std::shared_ptr<Acts::CutoutCylinderVolumeBounds> ActsCaloTrackingVolumeBuilder::makeCaloVolumeBounds(const std::vector<std::unique_ptr<Box>>& boxStore, std::shared_ptr<const Acts::TrackingVolume> insideVolume) const { using namespace Acts::VectorHelpers; // determine the dimensions of the double rmin_at_center = std::numeric_limits<double>::max(); double rmin_at_choke = std::numeric_limits<double>::max(); double rmax = std::numeric_limits<double>::lowest(); double zmin = std::numeric_limits<double>::max(); double zmax = std::numeric_limits<double>::lowest(); double cutout_zmin_abs = std::numeric_limits<double>::max(); // We need to figure out what the size of the inner cutout cylinder is // so we can make sure everything worked fine! // We check what the min radius at small z, and then we turn it around and // check z bounds at lower radii. for (const auto& box : boxStore) { Acts::Vector3 vmin = box->min().cast<double>(); Acts::Vector3 vmax = box->max().cast<double>(); double vrmin = perp(vmin); double vrmax = perp(vmax); rmin_at_choke = std::min(rmin_at_choke, std::min(vrmin, vrmax)); rmax = std::max(rmax, std::max(vrmin, vrmax)); zmin = std::min(zmin, std::min(vmin.z(), vmax.z())); zmax = std::max(zmax, std::max(vmin.z(), vmax.z())); if (std::abs(vmin.z()) < 100) { rmin_at_center = std::min(vrmin, rmin_at_center); } if (std::abs(vmax.z()) < 100) { rmin_at_center = std::min(vrmax, rmin_at_center); } } for (const auto& box : boxStore) { Acts::Vector3 vmin = box->min().cast<double>(); Acts::Vector3 vmax = box->max().cast<double>(); double vrmin = perp(vmin); double vrmax = perp(vmax); if (vrmin < rmin_at_center * 0.9) { cutout_zmin_abs = std::min(cutout_zmin_abs, std::abs(vmin.z())); } if (vrmax < rmin_at_center * 0.9) { cutout_zmin_abs = std::min(cutout_zmin_abs, std::abs(vmax.z())); } } double dz1 = (zmax - zmin) / 2.; double dz2 = cutout_zmin_abs; // envelopes double envZ = 5; double envR = 5; dz1 += envZ; dz2 -= envZ; rmax += envR; if(rmin_at_choke > envR) rmin_at_choke -= envR; rmin_at_center -= envR; ATH_MSG_VERBOSE("rmin_at_center: " << rmin_at_center << " rmin at choke: " << rmin_at_choke << " rmax: " << rmax << " zmin: " << zmin << " zmax: " << zmax << " coutout_zmin_abs: " << cutout_zmin_abs); // Ok now let's analyse what we're wrapping the calo around: the ID // The ID will have to be built already. // We need to figure out the dimensions of the ID. // Assuming the wrapping volume is a cylinder. auto idCylBds = dynamic_cast<const Acts::CylinderVolumeBounds*>(&insideVolume->volumeBounds()); if (idCylBds == nullptr) { // protection against nullptr ATH_MSG_ERROR("Unable to dynamic cast volume bounds to Acts::CylinderVolumeBounds"); throw std::runtime_error("Error casting to CylinderVolumeBounds"); } double idRMax = idCylBds->get(CVBBV::eMaxR); double idRMin = idCylBds->get(CVBBV::eMinR); double idHlZ = idCylBds->get(CVBBV::eHalfLengthZ); ATH_MSG_VERBOSE("ID volume bounds:\n" << *idCylBds); ATH_MSG_VERBOSE("Inside volume transform: \n" << insideVolume->transform().matrix()); if (!insideVolume->transform().isApprox(Acts::Transform3::Identity())) { ATH_MSG_VERBOSE("Inside volume transform is not unity."); // transformation matrix is NOT unity. Let's check: // - Rotation is approximate unity // - Translation is only along z axis const auto& trf = insideVolume->transform(); Acts::RotationMatrix3 rot = trf.rotation(); bool unityRot = rot.isApprox(Acts::RotationMatrix3::Identity()); ATH_MSG_VERBOSE("\n" << rot); // dot product with Z axis is about 1 => ok const Acts::Vector3 trl = trf.translation(); bool transZOnly = std::abs(1 - std::abs(Acts::Vector3::UnitZ().dot(trl.normalized()))) < 1e-6; ATH_MSG_VERBOSE("TRL "<< trl.transpose()); ATH_MSG_VERBOSE("TRL "<< trl.normalized().dot(Acts::Vector3::UnitZ())); if(!unityRot || !transZOnly) { ATH_MSG_ERROR("The ID appears to be shifted from the origin. I cannot handle this."); ATH_MSG_ERROR("(I'm not building the Calo!)"); throw std::runtime_error("Error building calo"); } else { ATH_MSG_VERBOSE("Checked: non unitarity is ONLY due to shift along z axis: that's ok"); double prevIdHlZ = idHlZ; idHlZ += std::abs(trl.z()); ATH_MSG_VERBOSE("Modifying effective half length of ID cylinder: " << prevIdHlZ << " => " << idHlZ); } } // make sure we can fit the ID inside the calo cutout if (idRMax > rmin_at_center || idHlZ > dz2 || (idRMin > rmin_at_choke && idRMin != 0.)) { ATH_MSG_ERROR("Cannot fit ID inside the Calo"); ATH_MSG_ERROR("This can be because the ID overlaps into the calo volume"); ATH_MSG_ERROR("Or because the Calo choke radius is SMALLER than the ID inner radius"); ATH_MSG_ERROR("Currently, I can only make the choke radius smaller, I can not make it larger"); ATH_MSG_ERROR("nor can I manipulate the ID volume bounds at this point."); ATH_MSG_ERROR("ID rMax: " << idRMax << " Calo rMin@center: " << rmin_at_center); ATH_MSG_ERROR("ID hlZ: " << idHlZ << " Calo inner Z hl: " << dz2); ATH_MSG_ERROR("ID rMin: " << idRMin << " Calo rMin@choke: " << rmin_at_choke); ATH_MSG_ERROR("(I'm not building the Calo!)"); throw std::runtime_error("Error building calo"); } // Let's harmonize the sizes, so we have a exact wrap of the ID // Choke is now exactly as wide as space inside the ID. // We can fit the beam pipe in there later. rmin_at_choke = idRMin; std::shared_ptr<Acts::CutoutCylinderVolumeBounds> volBds = nullptr; volBds = std::make_shared<Acts::CutoutCylinderVolumeBounds>( rmin_at_choke, rmin_at_center, rmax, dz1, dz2); ATH_MSG_VERBOSE(*volBds); return volBds; } namespace { double eta_to_theta(double eta) { return 2 * std::atan(std::exp(-eta)); } } Acts::Volume ActsCaloTrackingVolumeBuilder::build_endcap(double z, double dz, double eta, double deta, double phi, double dphi) const { double eta_max = eta + deta * 0.5; double eta_min = eta - deta * 0.5; double theta_max = eta_to_theta(eta_max); double theta = eta_to_theta(eta); double theta_min = eta_to_theta(eta_min); double phi_max = phi + dphi * 0.5; double phi_min = phi - dphi * 0.5; double z_min = z - dz; double z_max = z + dz; double r_min, r_max; // inner face r_min = std::tan(theta_min) * z_min; r_max = std::tan(theta_max) * z_min; Acts::Vector3 p1(r_min * std::cos(phi_min), r_min * std::sin(phi_min), z_min); Acts::Vector3 p2(r_min * std::cos(phi_max), r_min * std::sin(phi_max), z_min); Acts::Vector3 p3(r_max * std::cos(phi_max), r_max * std::sin(phi_max), z_min); Acts::Vector3 p4(r_max * std::cos(phi_min), r_max * std::sin(phi_min), z_min); // outer face r_min = std::tan(theta_min) * z_max; r_max = std::tan(theta_max) * z_max; Acts::Vector3 p5(r_min * std::cos(phi_min), r_min * std::sin(phi_min), z_max); Acts::Vector3 p6(r_min * std::cos(phi_max), r_min * std::sin(phi_max), z_max); Acts::Vector3 p7(r_max * std::cos(phi_max), r_max * std::sin(phi_max), z_max); Acts::Vector3 p8(r_max * std::cos(phi_min), r_max * std::sin(phi_min), z_max); double r_mid = std::tan(theta) * z_min; Acts::Vector3 center; center.x() = r_mid * std::cos(phi); center.y() = r_mid * std::sin(phi); center.z() = z; Acts::Transform3 glob2vol = Acts::Transform3::Identity(); glob2vol *= Acts::AngleAxis3(-phi, Acts::Vector3::UnitZ()); glob2vol *= Acts::AngleAxis3( -theta, Acts::Vector3::UnitZ().cross(center).normalized()); glob2vol *= Acts::Translation3(-(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.); p1 = glob2vol * p1; p2 = glob2vol * p2; p3 = glob2vol * p3; p4 = glob2vol * p4; p5 = glob2vol * p5; p6 = glob2vol * p6; p7 = glob2vol * p7; p8 = glob2vol * p8; auto globalToLocal = glob2vol.inverse(); auto cubo = std::make_shared<Acts::GenericCuboidVolumeBounds>( std::array<Acts::Vector3, 8>({{p1, p2, p3, p4, p5, p6, p7, p8}})); Acts::Volume vol(globalToLocal, std::move(cubo)); return vol; } Acts::Volume ActsCaloTrackingVolumeBuilder::build_barrel(double r, double dr, double eta, double deta, double phi, double dphi) const { // std::cout << "build barrel" << std::endl; double eta_max = eta + deta * 0.5; double eta_min = eta - deta * 0.5; double theta = eta_to_theta(eta); double theta_max = eta_to_theta(eta_max); double theta_min = eta_to_theta(eta_min); double phi_max = phi + dphi * 0.5; double phi_min = phi - dphi * 0.5; double r_min = r - dr; double r_max = r + dr; double z_min, z_max; // inner face z_min = r_min / std::tan(theta_min); z_max = r_min / std::tan(theta_max); Acts::Vector3 p1(r_min * std::cos(phi_min), r_min * std::sin(phi_min), z_min); Acts::Vector3 p2(r_min * std::cos(phi_min), r_min * std::sin(phi_min), z_max); Acts::Vector3 p3(r_min * std::cos(phi_max), r_min * std::sin(phi_max), z_max); Acts::Vector3 p4(r_min * std::cos(phi_max), r_min * std::sin(phi_max), z_min); // outer face z_min = r_max / std::tan(theta_min); z_max = r_max / std::tan(theta_max); Acts::Vector3 p5(r_max * std::cos(phi_min), r_max * std::sin(phi_min), z_min); Acts::Vector3 p6(r_max * std::cos(phi_min), r_max * std::sin(phi_min), z_max); Acts::Vector3 p7(r_max * std::cos(phi_max), r_max * std::sin(phi_max), z_max); Acts::Vector3 p8(r_max * std::cos(phi_max), r_max * std::sin(phi_max), z_min); Acts::Vector3 center; center.x() = r * std::cos(phi); center.y() = r * std::sin(phi); center.z() = r / std::tan(theta); Acts::Transform3 glob2vol = Acts::Transform3::Identity(); glob2vol *= Acts::AngleAxis3(-phi, Acts::Vector3::UnitZ()); glob2vol *= Acts::AngleAxis3( -theta, Acts::Vector3::UnitZ().cross(center).normalized()); glob2vol *= Acts::Translation3(-(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.); p1 = glob2vol * p1; p2 = glob2vol * p2; p3 = glob2vol * p3; p4 = glob2vol * p4; p5 = glob2vol * p5; p6 = glob2vol * p6; p7 = glob2vol * p7; p8 = glob2vol * p8; auto globalToLocal = glob2vol.inverse(); auto cubo = std::make_shared<Acts::GenericCuboidVolumeBounds>( std::array<Acts::Vector3, 8>({{p1, p2, p3, p4, p5, p6, p7, p8}})); Acts::Volume vol(globalToLocal, std::move(cubo)); return vol; } Acts::Volume ActsCaloTrackingVolumeBuilder::build_box(double x, double dx, double y, double dy, double z, double dz) const { // std::cout << "build box" << std::endl; double x_min, x_max, y_min, y_max, z_min, z_max; x_min = x - dx; x_max = x + dx; y_min = y - dy; y_max = y + dy; z_min = z - dz; z_max = z + dz; // inner face Acts::Vector3 p1(x_min, y_min, z_min); Acts::Vector3 p2(x_min, y_max, z_min); Acts::Vector3 p3(x_max, y_max, z_min); Acts::Vector3 p4(x_max, y_min, z_min); // outer face Acts::Vector3 p5(x_min, y_min, z_max); Acts::Vector3 p6(x_min, y_max, z_max); Acts::Vector3 p7(x_max, y_max, z_max); Acts::Vector3 p8(x_max, y_min, z_max); Acts::Transform3 glob2vol = Acts::Transform3::Identity(); glob2vol *= Acts::Translation3(-(p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8) / 8.); p1 = glob2vol * p1; p2 = glob2vol * p2; p3 = glob2vol * p3; p4 = glob2vol * p4; p5 = glob2vol * p5; p6 = glob2vol * p6; p7 = glob2vol * p7; p8 = glob2vol * p8; auto globalToLocal = glob2vol.inverse(); auto cubo = std::make_shared<Acts::GenericCuboidVolumeBounds>( std::array<Acts::Vector3, 8>({{p1, p2, p3, p4, p5, p6, p7, p8}})); Acts::Volume vol(globalToLocal, std::move(cubo)); return vol; } std::vector<std::unique_ptr<Acts::Volume>> ActsCaloTrackingVolumeBuilder::cellFactory() const { //Acts::ply_helper<double> ply_lar; //Acts::ply_helper<double> ply_tile; //Acts::ply_helper<double> ply_fcal; //float x, y, z, r, phi_raw, eta_raw, dphi, deta, dr, dx, dy, dz; float z, dz, eta_raw, deta, phi_raw, dphi, r, dr, x, y, dx, dy; size_t calosample; float scale; // storage of cells we will produce std::vector<std::unique_ptr<Acts::Volume>> cells; cells.reserve(m_caloMgr->element_size()); // about 180k for(auto it = m_caloMgr->element_begin();it < m_caloMgr->element_end();++it) { const CaloDetDescrElement* cde = *it; z = cde->z(); dz = cde->dz(); eta_raw = cde->eta_raw(); deta = cde->deta(); phi_raw = cde->phi_raw(); dphi = cde->dphi(); r = cde->r(); dr = cde->dr(); x = cde->x(); y = cde->y(); dx = cde->dx(); dy = cde->dy(); calosample = cde->getSampling(); scale = 1.; if (calosample >= 12 && calosample <= 20) { scale = 0.5; } //Acts::ply_helper<double>* ply; //if (calosample <= 11) { //ply = &ply_lar; //} else if (calosample <= 20) { //ply = &ply_tile; //} else { //ply = &ply_fcal; //} switch (calosample) { case 4: case 5: case 6: case 7: case 8: case 9: case 10: case 11: case 17: dz *= scale; cells.push_back(std::make_unique<Acts::Volume>( build_endcap(z, dz, eta_raw, deta, phi_raw, dphi))); break; case 0: case 1: case 2: case 3: case 12: case 13: case 14: case 15: case 16: case 18: case 19: case 20: dr *= scale; cells.push_back(std::make_unique<Acts::Volume>( build_barrel(r, dr, eta_raw, deta, phi_raw, dphi))); break; case 21: case 22: case 23: scale = 1.; dx *= scale; dy *= scale; // dz *= scale; cells.push_back(std::make_unique<Acts::Volume>( build_box(x, dx, y, dy, z, dz))); break; default: std::stringstream ss; ss << "Unkown calo sample " << calosample; std::runtime_error(ss.str()); } //cells.back()->boundingBox({0.1, 0.1, 0.1}).draw(*ply); //auto cvb = dynamic_cast<const //Acts::GenericCuboidVolumeBounds*>(&cells.back()->volumeBounds()); //cvb->draw(*ply, cells.back()->transform()); } //std::ofstream os("lar.ply"); //os << ply_lar << std::flush; //os.close(); //os = std::ofstream("tile.ply"); //os << ply_tile << std::flush; //os.close(); //os = std::ofstream("fcal.ply"); //os << ply_fcal << std::flush; //os.close(); return cells; }