diff --git a/Detectors/Passive/src/Cave.cxx b/Detectors/Passive/src/Cave.cxx index 208084a335ab5..bc28c852fca95 100644 --- a/Detectors/Passive/src/Cave.cxx +++ b/Detectors/Passive/src/Cave.cxx @@ -85,7 +85,7 @@ void Cave::ConstructGeometry() shCaveTR1->DefineSection(0, -706. - 8.6, 0., 790.5); shCaveTR1->DefineSection(1, 707. + 7.6, 0., 790.5); TGeoTube* shCaveTR2 = new TGeoTube("shCaveTR2", 0., 150., 110.); - TGeoTube* shCaveTR3 = new TGeoTube("shCaveTR3", 0., 80., 75.); + TGeoTube* shCaveTR3 = new TGeoTube("shCaveTR3", 0., 105., 75.); TGeoTranslation* transCaveTR2 = new TGeoTranslation("transTR2", 0, 30., -505. - 110.); TGeoTranslation* transCaveTR3 = new TGeoTranslation("transTR3", 0, 30., 714.6 + 75.); diff --git a/Detectors/ZDC/simulation/include/ZDCSimulation/Detector.h b/Detectors/ZDC/simulation/include/ZDCSimulation/Detector.h index 40e4babe8d760..2b9df85aa42a9 100644 --- a/Detectors/ZDC/simulation/include/ZDCSimulation/Detector.h +++ b/Detectors/ZDC/simulation/include/ZDCSimulation/Detector.h @@ -122,6 +122,10 @@ class Detector : public o2::base::DetImpl void createMagnets(); void createDetectors(); + /// Whether the ZEM calorimeters are built: FoCal occupies the same space, so not + /// when it is active. + static Bool_t withZEM(); + // determine detector; sector/tower and impact coordinates given volumename and position void getDetIDandSecID(TString const& volname, math_utils::Vector3D const& x, math_utils::Vector3D& xDet, int& detector, int& sector) const; diff --git a/Detectors/ZDC/simulation/src/Detector.cxx b/Detectors/ZDC/simulation/src/Detector.cxx index b8b81379a4dff..427f0a8ee387c 100644 --- a/Detectors/ZDC/simulation/src/Detector.cxx +++ b/Detectors/ZDC/simulation/src/Detector.cxx @@ -30,6 +30,8 @@ #include #include #include "ZDCSimulation/ZDCSimParam.h" +#include "SimConfig/SimConfig.h" +#include #ifdef ZDC_FASTSIM_ONNX #include "Utils.h" // for normal_distribution() #include "FastSimulations.h" // for fastsim module @@ -257,7 +259,7 @@ void Detector::defineSensitiveVolumes() AddSensitiveVolume(vol); mZEMVolID = vol->GetNumber(); AddSensitiveVolume(gGeoManager->GetVolume("ZEMF")); - } else { + } else if (withZEM()) { LOG(fatal) << "can't find volume ZEM"; } } @@ -2302,6 +2304,10 @@ void Detector::createDetectors() // ------------------------------------------------------------------------------- // -> EM calorimeter (ZEM) + if (!withZEM()) { + LOG(warning) << "ZDC: FoCal is active, so the ZEM calorimeters are not built -- they occupy the same space"; + return; + } int32_t irotzem1, irotzem2; double rangzem1[6] = {0., 0., 90., 90., -90., 0.}; double rangzem2[6] = {180., 0., 90., 45. + 90., 90., 45.}; @@ -2354,16 +2360,19 @@ void Detector::createDetectors() TVirtualMC::GetMC()->Gspos("ZEV1", 1, "ZETR", -zemVoidLayer[0] + zemTranLength, 0., 0., 0, "ONLY"); // --- Positioning the ZEM into the ZDC - rotation for 90 degrees - // NB -> ZEM is positioned in cave volume - const float z0 = 1313.3475 + 75.; // center of caveRB24 mother volume - TVirtualMC::GetMC()->Gspos("ZEM ", 1, "caveRB24", -Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + Geometry::ZEMDIMENSION[0] - z0, irotzem1, "ONLY"); + // NB -> ZEM is positioned in the barrel volume + // The ZEM calorimeters and their supports occupy 750 < z < 860 cm, a region that + // belongs to the barrel mother volume (caveRB24 starts only at z = 864.6 cm). + // The barrel is placed at y = -30 cm in the cave, hence the y0 offset. + const float y0 = 30.; + TVirtualMC::GetMC()->Gspos("ZEM ", 1, "barrel", -Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + Geometry::ZEMDIMENSION[0], irotzem1, "ONLY"); // Second EM ZDC (same side w.r.t. IP, just on the other side w.r.t. beam pipe) - TVirtualMC::GetMC()->Gspos("ZEM ", 2, "caveRB24", Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + Geometry::ZEMDIMENSION[0] - z0, irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEM ", 2, "barrel", Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + Geometry::ZEMDIMENSION[0], irotzem1, "ONLY"); // --- Adding last slice at the end of the EM calorimeter float zLastSlice = Geometry::ZEMPOSITION[2] + zemPbSlice[0] + 2 * Geometry::ZEMDIMENSION[0]; - TVirtualMC::GetMC()->Gspos("ZEL2", 1, "caveRB24", Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1], zLastSlice - z0, irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEL2", 1, "barrel", Geometry::ZEMPOSITION[0], Geometry::ZEMPOSITION[1] + y0, zLastSlice, irotzem1, "ONLY"); // ------------------------------------------------------------------------------- // -> ZEM supports @@ -2376,31 +2385,31 @@ void Detector::createDetectors() // Bridge TVirtualMC::GetMC()->Gsvolu("ZESH", "BOX ", getMediumID(kAl), const_cast(zemSupport1), 3); float ybridge = Geometry::ZEMPOSITION[1] - Geometry::ZEMDIMENSION[1] - 2. * 2. * zemSupportBox[3 + 1] - 5. - zemSupport1[1]; - TVirtualMC::GetMC()->Gspos("ZESH", 1, "caveRB24", Geometry::ZEMPOSITION[0], ybridge, zbox - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZESH", 2, "caveRB24", -Geometry::ZEMPOSITION[0], ybridge, zbox - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESH", 1, "barrel", Geometry::ZEMPOSITION[0], ybridge + y0, zbox, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESH", 2, "barrel", -Geometry::ZEMPOSITION[0], ybridge + y0, zbox, 0, "ONLY"); // TVirtualMC::GetMC()->Gsvolu("ZESV", "BOX ", getMediumID(kAl), const_cast(zemSupport2), 3); - TVirtualMC::GetMC()->Gspos("ZESV", 1, "caveRB24", Geometry::ZEMPOSITION[0] - zemSupportBox[0] + zemSupport2[0], ybox - zemSupportBox[1] - zemSupport2[1], zbox - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZESV", 2, "caveRB24", Geometry::ZEMPOSITION[0] + zemSupportBox[0] - zemSupport2[0], ybox - zemSupportBox[1] - zemSupport2[1], zbox - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZESV", 3, "caveRB24", -(Geometry::ZEMPOSITION[0] - zemSupportBox[0] + zemSupport2[0]), ybox - zemSupportBox[1] - zemSupport2[1], zbox - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZESV", 4, "caveRB24", -(Geometry::ZEMPOSITION[0] + zemSupportBox[0] - zemSupport2[0]), ybox - zemSupportBox[1] - zemSupport2[1], zbox - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESV", 1, "barrel", Geometry::ZEMPOSITION[0] - zemSupportBox[0] + zemSupport2[0], ybox - zemSupportBox[1] - zemSupport2[1] + y0, zbox, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESV", 2, "barrel", Geometry::ZEMPOSITION[0] + zemSupportBox[0] - zemSupport2[0], ybox - zemSupportBox[1] - zemSupport2[1] + y0, zbox, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESV", 3, "barrel", -(Geometry::ZEMPOSITION[0] - zemSupportBox[0] + zemSupport2[0]), ybox - zemSupportBox[1] - zemSupport2[1] + y0, zbox, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZESV", 4, "barrel", -(Geometry::ZEMPOSITION[0] + zemSupportBox[0] - zemSupport2[0]), ybox - zemSupportBox[1] - zemSupport2[1] + y0, zbox, 0, "ONLY"); // Table TVirtualMC::GetMC()->Gsvolu("ZETA", "BOX ", getMediumID(kAl), const_cast(zemSupportTable), 3); float ytable = ybridge - zemSupport1[1] - zemSupportTable[1]; - TVirtualMC::GetMC()->Gspos("ZETA", 1, "caveRB24", 0.0, ytable, zbox - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZETA", 2, "caveRB24", 0.0, ytable - 13. + 2. * zemSupportTable[1], zbox - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZETA", 1, "barrel", 0.0, ytable + y0, zbox, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZETA", 2, "barrel", 0.0, ytable - 13. + 2. * zemSupportTable[1] + y0, zbox, 0, "ONLY"); //Screens around ZEM TVirtualMC::GetMC()->Gsvolu("ZEFL", "BOX ", getMediumID(kAl), const_cast(zemSupport3), 3); - TVirtualMC::GetMC()->Gspos("ZEFL", 1, "caveRB24", Geometry::ZEMPOSITION[0], -Geometry::ZEMDIMENSION[1] - zemSupport3[1], zSupport + zemSupport3[2] - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEFL", 2, "caveRB24", -Geometry::ZEMPOSITION[0], -Geometry::ZEMDIMENSION[1] - zemSupport3[1], zSupport + zemSupport3[2] - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEFL", 1, "barrel", Geometry::ZEMPOSITION[0], -Geometry::ZEMDIMENSION[1] - zemSupport3[1] + y0, zSupport + zemSupport3[2], 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEFL", 2, "barrel", -Geometry::ZEMPOSITION[0], -Geometry::ZEMDIMENSION[1] - zemSupport3[1] + y0, zSupport + zemSupport3[2], 0, "ONLY"); TVirtualMC::GetMC()->Gsvolu("ZELA", "PARA", getMediumID(kAl), const_cast(zemSupport4), 6); - TVirtualMC::GetMC()->Gspos("ZELA", 1, "caveRB24", Geometry::ZEMPOSITION[0] - Geometry::ZEMDIMENSION[2] - zemSupport4[2], Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + zemSupport4[0] - z0, irotzem1, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZELA", 2, "caveRB24", Geometry::ZEMPOSITION[0] + Geometry::ZEMDIMENSION[2] + zemSupport4[2], Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + zemSupport4[0] - z0, irotzem1, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZELA", 3, "caveRB24", -(Geometry::ZEMPOSITION[0] - Geometry::ZEMDIMENSION[2] - zemSupport4[2]), Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + zemSupport4[0] - z0, irotzem1, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZELA", 4, "caveRB24", -(Geometry::ZEMPOSITION[0] + Geometry::ZEMDIMENSION[2] + zemSupport4[2]), Geometry::ZEMPOSITION[1], Geometry::ZEMPOSITION[2] + zemSupport4[0] - z0, irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZELA", 1, "barrel", Geometry::ZEMPOSITION[0] - Geometry::ZEMDIMENSION[2] - zemSupport4[2], Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + zemSupport4[0], irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZELA", 2, "barrel", Geometry::ZEMPOSITION[0] + Geometry::ZEMDIMENSION[2] + zemSupport4[2], Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + zemSupport4[0], irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZELA", 3, "barrel", -(Geometry::ZEMPOSITION[0] - Geometry::ZEMDIMENSION[2] - zemSupport4[2]), Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + zemSupport4[0], irotzem1, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZELA", 4, "barrel", -(Geometry::ZEMPOSITION[0] + Geometry::ZEMDIMENSION[2] + zemSupport4[2]), Geometry::ZEMPOSITION[1] + y0, Geometry::ZEMPOSITION[2] + zemSupport4[0], irotzem1, "ONLY"); // Containers for ZEM calorimeters TVirtualMC::GetMC()->Gsvolu("ZEW1", "BOX ", getMediumID(kAl), const_cast(zemWallH), 3); @@ -2410,22 +2419,37 @@ void Detector::createDetectors() // float yh1 = Geometry::ZEMPOSITION[1] - Geometry::ZEMDIMENSION[1] - 2 * zemSupport3[1] - zemWallH[1]; float zh1 = zSupport + zemWallH[2]; - TVirtualMC::GetMC()->Gspos("ZEW1", 1, "caveRB24", Geometry::ZEMPOSITION[0], yh1, zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW1", 2, "caveRB24", Geometry::ZEMPOSITION[0], yh1 + 2 * zemSupportBox[1], zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW1", 3, "caveRB24", -Geometry::ZEMPOSITION[0], yh1, zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW1", 4, "caveRB24", -Geometry::ZEMPOSITION[0], yh1 + 2 * zemSupportBox[1], zh1 - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW1", 1, "barrel", Geometry::ZEMPOSITION[0], yh1 + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW1", 2, "barrel", Geometry::ZEMPOSITION[0], yh1 + 2 * zemSupportBox[1] + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW1", 3, "barrel", -Geometry::ZEMPOSITION[0], yh1 + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW1", 4, "barrel", -Geometry::ZEMPOSITION[0], yh1 + 2 * zemSupportBox[1] + y0, zh1, 0, "ONLY"); // - TVirtualMC::GetMC()->Gspos("ZEW2", 1, "caveRB24", Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1], zSupport - zemWallVfwd[2] - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW3", 1, "caveRB24", Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1], zSupport + 2 * zemWallH[2] - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW2", 2, "caveRB24", -Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1], zSupport - zemWallVfwd[2] - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW3", 2, "caveRB24", -Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1], zSupport + 2 * zemWallH[2] - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW2", 1, "barrel", Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1] + y0, zSupport - zemWallVfwd[2], 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW3", 1, "barrel", Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1] + y0, zSupport + 2 * zemWallH[2], 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW2", 2, "barrel", -Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1] + y0, zSupport - zemWallVfwd[2], 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW3", 2, "barrel", -Geometry::ZEMPOSITION[0], yh1 + zemSupportBox[1] + y0, zSupport + 2 * zemWallH[2], 0, "ONLY"); // float xl1 = Geometry::ZEMPOSITION[0] - Geometry::ZEMDIMENSION[2] - 2. * zemSupport4[2] - zemWallVside[0]; float xl2 = Geometry::ZEMPOSITION[0] + Geometry::ZEMDIMENSION[2] + 2. * zemSupport4[2] + zemWallVside[0]; - TVirtualMC::GetMC()->Gspos("ZEW4", 1, "caveRB24", xl1, yh1 + zemSupportBox[1], zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW4", 2, "caveRB24", xl2, yh1 + zemSupportBox[1], zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW4", 3, "caveRB24", -xl1, yh1 + zemSupportBox[1], zh1 - z0, 0, "ONLY"); - TVirtualMC::GetMC()->Gspos("ZEW4", 4, "caveRB24", -xl2, yh1 + zemSupportBox[1], zh1 - z0, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW4", 1, "barrel", xl1, yh1 + zemSupportBox[1] + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW4", 2, "barrel", xl2, yh1 + zemSupportBox[1] + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW4", 3, "barrel", -xl1, yh1 + zemSupportBox[1] + y0, zh1, 0, "ONLY"); + TVirtualMC::GetMC()->Gspos("ZEW4", 4, "barrel", -xl2, yh1 + zemSupportBox[1] + y0, zh1, 0, "ONLY"); +} + +//_____________________________________________________________________________ +/// FoCal occupies the space the ZEM calorimeters sit in -- every point of ZEM lies +/// inside the FOCAL box -- so the two cannot both be built. Which of them the Run 4 +/// layout keeps, and where a Run 4 ZEM would go, is for the ZDC and FoCal groups to +/// settle; until they do, FoCal wins, which is what the geometry did by accident +/// while ZEM was placed outside its mother and unreachable. +/// +/// An empty active-module list (a geometry built outside o2-sim) means Run 3, and +/// ZEM is built. +Bool_t Detector::withZEM() +{ + const auto& modules = o2::conf::SimConfig::Instance().getActiveModules(); + return std::find(modules.begin(), modules.end(), "FOC") == modules.end(); } //_____________________________________________________________________________ diff --git a/run/o2sim_geometry_doctor.cxx b/run/o2sim_geometry_doctor.cxx index da9072f5b9238..1664580766276 100644 --- a/run/o2sim_geometry_doctor.cxx +++ b/run/o2sim_geometry_doctor.cxx @@ -22,6 +22,12 @@ /// * where does it already express it wrongly -- a medium with ifield == 0, /// which asks the transport engine for straight lines, sitting in real field. /// +/// It also answers one question about the geometry alone, and therefore runs that +/// part without a field under --reachability-only: does every placement actually +/// occupy the space it was built in? A daughter outside its mother, or one shadowed +/// by an overlapping sibling, is never reached by the navigator, carries no material +/// and produces no hits, and nothing in the construction code says so. +/// /// It detects and proposes. It never modifies a geometry. /// /// The field enters in two ways. A support model (an outer bound on where |B| @@ -697,6 +703,166 @@ bool supportFromJson(const json& in, Support& support) return true; } +// --------------------------------------------------------------------------- +// reachability +// --------------------------------------------------------------------------- + +/// A placement can be perfectly built and still occupy no space. TGeo never +/// descends into a daughter lying outside its mother, and an overlapping sibling +/// can shadow one that does not; either way the volume carries no material, takes +/// no steps and produces no hits, and nothing in the construction code complains. +/// Only the navigator can settle it, so ask it: draw points inside a placement's +/// own shape and check that FindNode() comes back through that placement. +/// +/// One representative path per node object, which is the granularity at which the +/// common defect lives -- a daughter outside its mother is a property of the node, +/// not of the path that reaches it. A node whose mother is itself placed many +/// times is therefore sampled once, in the first of those placements. +struct Reach { + std::string medium, mother, worstPath; + long sampled = 0; + double fraction = 1.; +}; + +constexpr int kReachRejectionTries = 400; + +class ReachAudit +{ + public: + explicit ReachAudit(int samples) : mSamples(samples) {} + + void walk(TGeoNode* node) { walk(node, TGeoHMatrix(), ""); } + + const std::vector& dead() const { return mDead; } + const std::vector& partial() const { return mPartial; } + long nodesVisited() const { return mVisited; } + long nodesSampled() const { return mSampled; } + long nodesUnsampleable() const { return mUnsampleable; } + + private: + void walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path); + bool samplePoint(TGeoShape* shape, double* local); + + int mSamples; + long mVisited = 0, mSampled = 0, mUnsampleable = 0; + TRandom3 mRandom{20260901}; + std::set mSeen; + std::vector mDead, mPartial; +}; + +/// Rejection sampling against the shape itself. A TGeoCompositeShape inherits +/// TGeoBBox, so its DX/DY/DZ describe a box that still contains the holes and +/// subtractions -- only Contains() knows the difference. GetOrigin() matters too: +/// the box need not be centred on the local origin. +bool ReachAudit::samplePoint(TGeoShape* shape, double* local) +{ + auto* box = dynamic_cast(shape); + if (box == nullptr) { + return false; + } + const double* origin = box->GetOrigin(); + for (int attempt = 0; attempt < kReachRejectionTries; ++attempt) { + local[0] = origin[0] + box->GetDX() * (2. * mRandom.Rndm() - 1.); + local[1] = origin[1] + box->GetDY() * (2. * mRandom.Rndm() - 1.); + local[2] = origin[2] + box->GetDZ() * (2. * mRandom.Rndm() - 1.); + if (shape->Contains(local)) { + return true; + } + } + return false; +} + +void ReachAudit::walk(TGeoNode* node, const TGeoHMatrix& parent, const std::string& path) +{ + if (!mSeen.insert(node).second) { + return; // this node object, and therefore its whole subtree, is already covered + } + TGeoHMatrix here = parent; + here.Multiply(node->GetMatrix()); + TGeoVolume* volume = node->GetVolume(); + const std::string myPath = path + "/" + node->GetName(); + ++mVisited; + + // an assembly is expanded away at closure, so FindNode never returns one + if (!volume->IsAssembly()) { + int drawn = 0, reached = 0; + for (int i = 0; i < mSamples; ++i) { + double local[3], global[3]; + if (!samplePoint(volume->GetShape(), local)) { + break; + } + ++drawn; + here.LocalToMaster(local, global); + if (gGeoManager->FindNode(global[0], global[1], global[2]) == nullptr) { + continue; + } + const std::string found = gGeoManager->GetPath(); + // reached if the navigator's own path passes through this placement + if (found.rfind(myPath, 0) == 0) { + ++reached; + } + } + if (drawn == 0) { + ++mUnsampleable; // a sliver too thin for the rejection budget; says nothing + } else { + ++mSampled; + const double fraction = double(reached) / drawn; + if (fraction < 0.999) { + auto* medium = volume->GetMedium(); + Reach entry; + entry.medium = medium != nullptr ? medium->GetName() : "(none)"; + entry.mother = node->GetMotherVolume() != nullptr ? node->GetMotherVolume()->GetName() : "-"; + entry.worstPath = myPath; + entry.sampled = drawn; + entry.fraction = fraction; + (fraction == 0. ? mDead : mPartial).push_back(entry); + } + } + } + + for (int i = 0; i < node->GetNdaughters(); ++i) { + walk(node->GetDaughter(i), here, myPath); + } +} + +/// Prints the audit and returns how many placements the navigator cannot reach at all. +long reportReachability(int samples, Report& report) +{ + if (samples <= 0) { + return 0; + } + progress("reachability: asking the navigator to find every placement from inside its own shape"); + ReachAudit audit(samples); + audit.walk(gGeoManager->GetTopNode()); + + report(form("reachability: %ld node objects visited, %ld sampled, %ld too thin to sample", + audit.nodesVisited(), audit.nodesSampled(), audit.nodesUnsampleable())); + report(form(" %ld placements the navigator never reaches, %zu it reaches only in part", + (long)audit.dead().size(), audit.partial().size())); + if (!audit.dead().empty()) { + report(" unreachable -- these carry no material and produce no hits:"); + report(form(" %-12s %-18s %10s %s", "mother", "medium", "sampled", "path")); + for (const auto& entry : audit.dead()) { + report(form(" %-12s %-18s %10ld %s", entry.mother.c_str(), entry.medium.c_str(), entry.sampled, + entry.worstPath.c_str())); + } + } + for (size_t i = 0; i < audit.partial().size() && i < 20; ++i) { + const auto& entry = audit.partial()[i]; + if (i == 0) { + report(" partially shadowed -- an overlapping sibling or an extruding placement:"); + report(form(" %-12s %-18s %8s %s", "mother", "medium", "reached", "path")); + } + report(form(" %-12s %-18s %7.1f%% %s", entry.mother.c_str(), entry.medium.c_str(), + 100. * entry.fraction, entry.worstPath.c_str())); + } + if (audit.partial().size() > 20) { + report(form(" ... and %zu more", audit.partial().size() - 20)); + } + report(""); + return (long)audit.dead().size(); +} + // --------------------------------------------------------------------------- // the placement table // --------------------------------------------------------------------------- @@ -1513,6 +1679,8 @@ struct Options { std::string anchorFile; std::vector thresholdsGauss; double margin = 5.0; + int reachSamples = 32; + bool reachabilityOnly = false; std::string outputPrefix = "geometry-doctor"; }; @@ -1541,7 +1709,11 @@ int main(int argc, char** argv) ("output-prefix", bpo::value(&options.outputPrefix)->default_value("geometry-doctor"), // "prefix for the report, the proposals and the placement table") // ("verify-anchors", bpo::value(&options.anchorFile), // - "check the classification against known-good volumes listed in this JSON file"); + "check the classification against known-good volumes listed in this JSON file") // + ("reachability-samples", bpo::value(&options.reachSamples)->default_value(32), // + "points drawn inside each placement for the reachability audit; 0 disables it") // + ("reachability-only", bpo::bool_switch(&options.reachabilityOnly), // + "run only the reachability audit, which needs no magnetic field"); bpo::variables_map arguments; try { @@ -1559,6 +1731,29 @@ int main(int argc, char** argv) const bool haveFieldFile = arguments.count("field-file") != 0u; const bool haveFieldCurrent = arguments.count("field-current") != 0u; + + // The reachability audit is a question about the geometry alone, so it is the one + // part of this tool that can run without a field. + if (options.reachabilityOnly) { + TGeoManager::Import(options.geometryFile.c_str()); + if (gGeoManager == nullptr) { + std::cerr << "error: no TGeoManager in " << options.geometryFile << '\n'; + return 1; + } + Report report; + report("ALICE simulation geometry doctor -- reachability audit"); + report(""); + report(" geometry : " + options.geometryFile); + report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), + gGeoManager->GetListOfMedia()->GetEntries())); + report(""); + const long dead = reportReachability(options.reachSamples, report); + const std::string reportPath = options.outputPrefix + "-report.txt"; + report("wrote " + reportPath); + report.write(reportPath); + return dead == 0 ? 0 : 3; + } + if (haveFieldFile == haveFieldCurrent) { std::cerr << "error: give exactly one of --field-file and --field-current\n"; return 1; @@ -1668,6 +1863,8 @@ int main(int argc, char** argv) } report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), gGeoManager->GetListOfMedia()->GetEntries())); + report(""); + reportReachability(options.reachSamples, report); Doctor doctor(field, support); doctor.walk(gGeoManager->GetTopNode());