diff --git a/CMakeLists.txt b/CMakeLists.txt index 44b0d0a..45151cd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -5,7 +5,7 @@ project(k4ActsTracking) find_package(ROOT COMPONENTS RIO Tree) find_package(EDM4HEP) find_package(k4FWCore) -find_package(Acts) +find_package(Acts REQUIRED COMPONENTS Core PluginDD4hep) find_package(DD4hep) #--------------------------------------------------------------- diff --git a/k4ActsTracking/CMakeLists.txt b/k4ActsTracking/CMakeLists.txt index 20b7edd..3a86731 100644 --- a/k4ActsTracking/CMakeLists.txt +++ b/k4ActsTracking/CMakeLists.txt @@ -8,9 +8,10 @@ gaudi_add_module(k4ActsTrackingPlugins file(GLOB _dd4hep_plugin_sources src/dd4hepplugin/*.cpp) add_dd4hep_plugin(k4ActsDD4hepPlugin SHARED ${_dd4hep_plugin_sources}) target_link_libraries(k4ActsDD4hepPlugin DD4hep::DDCore DD4hep::DDRec DD4hep::DDG4 DD4hep::DDParsers) +target_link_libraries(k4ActsDD4hepPlugin ActsCore ActsPluginDD4hep) -SET(compile_flags -Wformat-security -Wno-long-long -Wdeprecated -fdiagnostics-color=auto -Wall -Wextra -pedantic) +SET(compile_flags -Wformat-security -Wno-long-long -Wdeprecated -Wall -Wextra -pedantic) target_compile_options(k4ActsTrackingPlugins PUBLIC ${compile_flags}) target_compile_options(k4ActsDD4hepPlugin PUBLIC ${compile_flags}) diff --git a/k4ActsTracking/options/run_acts_plugin.py b/k4ActsTracking/options/run_acts_plugin.py index 4cd4970..5d26916 100644 --- a/k4ActsTracking/options/run_acts_plugin.py +++ b/k4ActsTracking/options/run_acts_plugin.py @@ -3,19 +3,19 @@ from Gaudi.Configuration import * from Configurables import ActsAlg, GeoSvc + algList = [] geosvc = GeoSvc("GeoSvc") # using this instead of CLIC_o3_v14 because it is much faster to instantiate -geosvc.detectors = [os.environ["LCGEO"]+"/CLIC/compact/CLIC_o2_v04/CLIC_o2_v04.xml"] +geosvc.detectors = [os.environ["LCGEO"] + "/CLIC/compact/CLIC_o2_v04/CLIC_o2_v04.xml"] +# geosvc.detectors = [os.environ["LCGEO"] + "/CLIC/compact/CLIC_o3_v14/CLIC_o3_v14.xml"] -a = ActsAlg("MyActsAlg") +a = ActsAlg("MyActsAlg", OutputLevel=VERBOSE) algList.append(a) from Configurables import ApplicationMgr -ApplicationMgr( TopAlg = algList, - EvtSel = 'NONE', - EvtMax = 2, - ExtSvc = [geosvc], - OutputLevel=INFO - ) + +ApplicationMgr( + TopAlg=algList, EvtSel="NONE", EvtMax=2, ExtSvc=[geosvc], OutputLevel=INFO +) diff --git a/k4ActsTracking/src/dd4hepplugin/actsdd4hepplugin.cpp b/k4ActsTracking/src/dd4hepplugin/actsdd4hepplugin.cpp index 27ef2cc..7f26e8e 100644 --- a/k4ActsTracking/src/dd4hepplugin/actsdd4hepplugin.cpp +++ b/k4ActsTracking/src/dd4hepplugin/actsdd4hepplugin.cpp @@ -9,6 +9,14 @@ #include #include +#include +#include +#include +#include +#include +#include + +#include #include using dd4hep::DetElement; @@ -20,6 +28,177 @@ namespace { * */ + std::pair, double> convertShapePlane(const TGeoShape& shape, + const std::string& axes, + Acts::LoggerWrapper logger) { + ACTS_VERBOSE("Converting TGeoShape to plane with local axes: " << axes); + + std::shared_ptr bounds; + double thickness; + if (const auto* trapezoid1 = dynamic_cast(&shape); trapezoid1 != nullptr) { + ACTS_VERBOSE("Have TGeoTrd1"); + throw std::runtime_error{"TGeoTrd1 not implemented yet"}; + } else if (const auto* trapezoid2 = dynamic_cast(&shape); trapezoid2 != nullptr) { + ACTS_VERBOSE("Have TGeoTrd2"); + + double dx1 = trapezoid2->GetDx1(); + double dx2 = trapezoid2->GetDx2(); + double dy1 = trapezoid2->GetDy1(); + double dy2 = trapezoid2->GetDy2(); + double dz = trapezoid2->GetDz(); + ACTS_VERBOSE(dx1 << " " << dx2 << " " << dy1 << " " << dy2 << " " << dz); + + if (axes == "zxy") { + if (std::abs(dy1 - dy2) > Acts::s_epsilon) { + throw std::runtime_error{"Inconsistent axes"}; + } + bounds = std::make_shared(dx1, dx2, dz); + thickness = dy1; + } else { + throw std::runtime_error{"Axes " + axes + " not yet implemented"}; + } + + return {bounds, thickness}; + } else if (const auto* box = dynamic_cast(&shape); box != nullptr) { + // this should always be true + ACTS_VERBOSE("Have Box"); + double dx = box->GetDX(); + double dy = box->GetDY(); + double dz = box->GetDZ(); + ACTS_VERBOSE(dx << " " << dy << " " << dz); + if (axes == "xyz") { + bounds = std::make_shared(dx, dy); + thickness = dz; + } else if (axes == "yzx") { + bounds = std::make_shared(dy, dz); + thickness = dx; + } else { + throw std::runtime_error{"Axes " + axes + " not yet implemented"}; + } + return {bounds, thickness}; + + } else { + throw std::runtime_error{"Unable to convert TGeoShape"}; + } + } + + // std::shared_ptr + std::shared_ptr createDetectorElement(const dd4hep::rec::Surface& srf) { + auto _logger = Acts::getDefaultLogger("DDRec2Acts", Acts::Logging::VERBOSE); + Acts::LoggerWrapper logger{*_logger}; + + ACTS_VERBOSE("Converting dd4hep::rec::Surface: checking SurfaceType"); + + auto getAxes = [&]() -> std::pair { + ACTS_VERBOSE("Interrogating surface for local coordinate system"); + + auto _origin = srf.globalToLocal(srf.origin()); + Acts::Vector2 origin{_origin.u(), _origin.v()}; + ACTS_VERBOSE(" - origin: " << origin.transpose() << " global " << srf.origin()); + + auto _normal = srf.normal(srf.origin()); + Acts::Vector3 normal{_normal.x(), _normal.y(), _normal.z()}; + + // ACTS_VERBOSE(" -> u.v = " << localU.dot(localV)); + + const TGeoNode& placement = *srf.detElement().placement().ptr(); + const TGeoMatrix& geoTransform = srf.detElement().nominal().worldTransformation(); + + Acts::Transform3 transform; + transform.translation() << geoTransform.GetTranslation()[0], geoTransform.GetTranslation()[1], + geoTransform.GetTranslation()[2]; + + const auto* rot = geoTransform.GetRotationMatrix(); + + std::array columns{Acts::Vector3{rot[0], rot[3], rot[6]}, Acts::Vector3{rot[1], rot[4], rot[7]}, + Acts::Vector3{rot[2], rot[5], rot[8]}}; + + std::array axes{'x', 'y', 'z'}; + + auto rotate = [](auto& a) { std::rotate(std::rbegin(a), std::rbegin(a) + 1, std::rend(a)); }; + + ACTS_VERBOSE("trafo: \n" << transform.matrix()); + ACTS_VERBOSE("origin loc: " << (transform * Acts::Vector3::Zero()).transpose()); + + bool found = false; + + for (size_t i = 0; i < 3; i++) { + transform.matrix().block<3, 1>(0, 0) = columns[0]; + transform.matrix().block<3, 1>(0, 1) = columns[1]; + transform.matrix().block<3, 1>(0, 2) = columns[2]; + + Acts::Vector3 normalLoc = transform.inverse().linear() * normal; + ACTS_VERBOSE("Axes: " << axes[0] << axes[1] << axes[2]); + ACTS_VERBOSE("normal global: " << normal.transpose()); + ACTS_VERBOSE(" -> local: " << normalLoc.transpose()); + + // ACTS_VERBOSE(normalLoc.dot(Acts::Vector3::UnitZ()) - 1); + + if (std::abs(normalLoc.dot(Acts::Vector3::UnitZ()) - 1) < Acts::s_epsilon) { + ACTS_VERBOSE("Normal accepted"); + found = true; + break; + } + + rotate(axes); + rotate(columns); + } + + if (!found) { + ACTS_ERROR("Unable to find permutation of axes for correct surface normal"); + throw std::runtime_error("Unable to find permutation of axis for correct surface normal"); + } + + auto _globalU = srf.u(srf.origin()); + Acts::Vector3 globalU{_globalU.x(), _globalU.y(), _globalU.z()}; + Acts::Vector3 localU = transform.inverse().linear() * globalU; + ACTS_VERBOSE(" - u: global " << globalU.transpose() << " local " << localU.transpose()); + auto _globalV = srf.v(srf.origin()); + Acts::Vector3 globalV{_globalV.x(), _globalV.y(), _globalV.z()}; + Acts::Vector3 localV = transform.inverse().linear() * globalV; + ACTS_VERBOSE(" - v: global " << globalV.transpose() << " local " << localV.transpose()); + + double localRotAngle = -std::acos(localU.normalized().dot(Acts::Vector3::UnitX())); + ACTS_VERBOSE("Local rot angle: " << (localRotAngle)); + + transform = transform * Acts::AngleAxis3{localRotAngle, Acts::Vector3::UnitZ()}; + + std::string saxes{axes.data(), 3}; + // ACTS_VERBOSE(saxes); + + return {transform, saxes}; + + // sanity check + + // ACTS_VERBOSE((transform.linear() * Acts::Vector3::UnitZ()).transpose()); + }; + + if (srf.type().isCylinder()) { + ACTS_VERBOSE("Surface is Cylinder"); + } else if (srf.type().isPlane()) { + ACTS_VERBOSE("Is 'Plane' (might be disc)"); + + auto [transform, axes] = getAxes(); + auto [bounds, thickness] = + convertShapePlane(*srf.detElement().placement().ptr()->GetVolume()->GetShape(), axes, logger); + + std::stringstream ss; + bounds->toStream(ss); + ACTS_VERBOSE("Created bounds: " << ss.str()); + ACTS_VERBOSE("Thickness: " << thickness); + + auto surface = Acts::Surface::makeShared(transform, bounds); + return surface; + + } else { + ACTS_ERROR("Surface type " << srf.type() << " currently unsupported"); + throw std::runtime_error("Unsupported surface type"); + } + + return nullptr; + // return std::make_shared(srf.detElement(), "Xz"); + } + static long addActsExtensions(dd4hep::Detector& description, int, char**) { const std::string LOG_SOURCE("AddActsExtensions"); printout(PrintLevel::INFO, LOG_SOURCE, "Running plugin"); @@ -27,6 +206,8 @@ namespace { //Getting the surface manager dd4hep::rec::SurfaceManager& surfMan = *description.extension(); + std::map, std::vector>> layers{}; + dd4hep::rec::SurfaceHelper ds(description.world()); dd4hep::rec::SurfaceList const& detSL = ds.surfaceList(); int counter = 0; @@ -38,28 +219,54 @@ namespace { continue; dd4hep::BitFieldCoder* cellIDcoder = nullptr; + + // std::cout << ddsurf->type() << std::endl; + // printout(PrintLevel::INFO, LOG_SOURCE, "%u", ddsurf->type()); if (ddsurf->type().isSensitive() and not ddsurf->type().isHelper()) { // surface is sensitive, assuming it has a readout auto detectorName = ddsurf->detElement().type(); try { auto theDetector = dd4hep::SensitiveDetector(description.sensitiveDetector(detectorName)); cellIDcoder = theDetector.readout().idSpec().decoder(); - printout(PrintLevel::DEBUG, LOG_SOURCE, "The encoding string for detector %s is %s", detectorName.c_str(), - cellIDcoder->fieldDescription().c_str()); + // printout(PrintLevel::DEBUG, LOG_SOURCE, "The encoding string for detector %s is %s", detectorName.c_str(), + // cellIDcoder->fieldDescription().c_str()); } catch (...) { printout(PrintLevel::INFO, LOG_SOURCE, "Not readout for detector %s", detectorName.c_str()); continue; } } else { + printout(PrintLevel::INFO, LOG_SOURCE, "Not sensitive"); continue; } counter++; + int system = cellIDcoder->get(ddsurf->id(), "system"); + // if (system != 0) { + // continue; + // } + std::string path = ddsurf->detElement().path(); auto layerNumber = cellIDcoder->get(ddsurf->id(), "layer"); printout(PrintLevel::INFO, LOG_SOURCE, "Found Surface at %s in layer %d at radius %3.2f mm", path.c_str(), layerNumber, ddsurf->origin().rho()); + for (const auto& field : cellIDcoder->fields()) { + // std::cout << "field: " << field.name() << std::endl; + std::cout << field.name() << ": " << cellIDcoder->get(ddsurf->id(), field.name()) << std::endl; + } + std::cout << "---" << std::endl; + + std::shared_ptr surface = createDetectorElement(*ddsurf); + + std::tuple key{ + system, cellIDcoder->get(ddsurf->id(), "layer"), + // cellIDcoder->get(ddsurf->id(), "barrel"), + }; + + layers[key].push_back(surface); + + // break; + //fixme: replace this with the ACTS type to be used as an extension dd4hep::rec::DoubleParameters* para = nullptr; try { // use existing map, or create a new one @@ -70,11 +277,25 @@ namespace { para->doubleParameters["SortingPolicy"] = 321; ddsurf->detElement().addExtension(para); } - if (counter > 1000) { - break; - } + // if (counter > 1000) { + // break; + // } } // for all surfaces + Acts::GeometryContext gctx; + for (const auto& [key, elements] : layers) { + Acts::ObjVisualization3D vis; + const auto [system, layer] = key; + for (const auto& element : elements) { + Acts::GeometryView3D::drawSurface(vis, *element, gctx); + } + std::stringstream ss; + ss << "obj/" << system << "_" << layer << "" + << ".obj"; + std::ofstream ofs{ss.str()}; + vis.write(ofs); + } + return 42; }