diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 063745e757816..db00f0465e645 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -108,6 +108,16 @@ o2_add_executable(g4-determine-unknown-pdg-properties SOURCES g4DetermineUnknownPdgProperties.cxx PUBLIC_LINK_LIBRARIES O2::O2TrivialMC MC::Geant4 MC::Geant4VMC O2::SimConfig) +# Reads a geometry file and a field and reports where the two do not fit together. +# It needs neither the simulation engines nor any detector code -- the medium flags +# it reasons about are GSTMED parameters stored in the geometry file itself -- so it +# deliberately does not link internal::allsim. +o2_add_executable(geometry-doctor + COMPONENT_NAME sim + SOURCES o2sim_geometry_doctor.cxx + PUBLIC_LINK_LIBRARIES O2::Field ROOT::Geom Boost::program_options + nlohmann_json::nlohmann_json) + o2_add_executable(mctracks-proxy COMPONENT_NAME sim SOURCES o2sim_mctracks_proxy.cxx diff --git a/run/geometry-doctor-anchors.json b/run/geometry-doctor-anchors.json new file mode 100644 index 0000000000000..3c8bef47740e4 --- /dev/null +++ b/run/geometry-doctor-anchors.json @@ -0,0 +1,15 @@ +{ + "comment": "Volumes of the ALICE Run 3 geometry whose relation to the field is known independently of this tool, from the geometry of the LHC machine elements around the RB24 beam pipe. Pass this file to o2-sim-geometry-doctor --verify-anchors to check a run against them. The five copper tubes sit 14 to 63 cm clear of the machine band in z and must come out field-free; RB243 and RB243A are assemblies, so they carry no material and nothing can be flagged on them; caveRB24 must NOT come out field-free, because what its daughters leave over is a sliver of cave air between an oval beam pipe and a circular field cylinder that carries more than 13 kG; and RB24B1BellowM is placed in the solenoid fringe while carrying a field-free medium, which is what the reverse audit exists to find.", + "anchors": [ + {"volume": "voRB242CuTube", "expect": "OUT", "placements": 2}, + {"volume": "voRB243CuTube", "expect": "OUT", "placements": 2}, + {"volume": "voRB24CuTubeF", "expect": "OUT", "placements": 7}, + {"volume": "voRB24cCuTube", "expect": "OUT", "placements": 1}, + {"volume": "voRB24cCuTubeM", "expect": "OUT", "placements": 1}, + {"volume": "RB243", "expect": "ASSEMBLY", "placements": 1}, + {"volume": "RB243A", "expect": "ASSEMBLY", "placements": 2}, + {"volume": "caveRB24", "expect": "NOT_OUT", "min_max_B_kG": 13.0}, + {"volume": "voRB24CuTubeM", "expect": "IN", "min_max_B_kG": 4.0}, + {"volume": "RB24B1BellowM", "expect": "REVERSE_AUDIT_FLAGGED", "min_max_B_kG": 4.0} + ] +} diff --git a/run/o2sim_geometry_doctor.cxx b/run/o2sim_geometry_doctor.cxx new file mode 100644 index 0000000000000..da9072f5b9238 --- /dev/null +++ b/run/o2sim_geometry_doctor.cxx @@ -0,0 +1,1775 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file o2sim_geometry_doctor.cxx +/// \brief Audits a placed geometry against the magnetic field it will be transported in. +/// +/// The tool reads a geometry file and a magnetic field, and answers two questions +/// that cannot be answered by looking at either one alone: +/// +/// * where does the geometry fail to express the field's structure -- a logical +/// volume placed both inside and outside the field, or a mother whose own +/// material straddles the boundary, so that no per-medium field flag can be +/// right for all of its placements; +/// * 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 detects and proposes. It never modifies a geometry. +/// +/// The field enters in two ways. A support model (an outer bound on where |B| +/// exceeds a threshold, maximised over phi) supplies the geometric argument: a +/// placement whose extent misses every band of the model is field-free, and that +/// is the only way this tool ever concludes "field-free". Sampling enters only in +/// the opposite direction, to disprove -- finding field inside a volume settles +/// the question, finding none does not. +/// +/// Reading the geometry needs no detector code: a medium's ifield flag and its +/// sensitivity are GSTMED parameters 1 and 0 of the TGeoMedium, and are recovered +/// from the file itself. + +#include "Field/MagneticField.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace bpo = boost::program_options; +using json = nlohmann::json; + +namespace +{ + +// --------------------------------------------------------------------------- +// small helpers +// --------------------------------------------------------------------------- + +/// printf into a std::string, for the tables this tool prints and stores. +template +std::string form(const char* fmt, Args... args) +{ + char buffer[4096]; + std::snprintf(buffer, sizeof(buffer), fmt, args...); + return std::string(buffer); +} + +/// Writes every line to stdout as it is produced and keeps a copy for the report file. +class Report +{ + public: + void operator()(const std::string& line) + { + std::cout << line << '\n'; + mLines.push_back(line); + } + void write(const std::string& path) const + { + std::ofstream out(path); + for (const auto& line : mLines) { + out << line << '\n'; + } + } + + private: + std::vector mLines; +}; + +/// Progress and diagnostics go to stderr so that the report on stdout stays a report. +void progress(const std::string& line) { std::cerr << line << std::endl; } + +// --------------------------------------------------------------------------- +// the field +// --------------------------------------------------------------------------- + +// A serialized MagneticField keeps its measured map in a transient member, so a +// reader has to call CreateField() again to get a usable field back. That call +// carries a trap: it feeds mMultipicativeFactorSolenoid/Dipole back through +// setters which negate under the LHC polarity convention, so CreateField() is not +// idempotent and the second call inverts the polarity of both the measured map +// and the machine compensators. |B| is untouched, which is exactly why the flip +// survives any magnitude-based check -- a 200k-point round trip on |B| passes +// while every field vector points the wrong way. +// +// A field file therefore has to carry reference field VECTORS taken from the live +// object at write time. The loader re-evaluates them, repairs a pure global flip, +// and refuses the file if what it gets back is anything other than what was +// written. A file without those probes cannot be verified and is refused too; +// --field-current builds the field from scratch instead. + +constexpr const char* kFieldObjectKey = "MagneticField"; +constexpr const char* kFieldProbeKey = "ReferenceProbes"; + +/// 0 = reproduces the probes, 1 = exactly negated, -1 = neither. +int compareToProbes(o2::field::MagneticField* field, const TVectorD& probes) +{ + const int n = probes.GetNrows() / 6; + bool same = true; + bool flipped = true; + for (int i = 0; i < n; ++i) { + double x[3] = {probes[6 * i], probes[6 * i + 1], probes[6 * i + 2]}; + double b[3] = {0., 0., 0.}; + field->Field(x, b); + for (int k = 0; k < 3; ++k) { + const double want = probes[6 * i + 3 + k]; + same = same && (b[k] == want); + flipped = flipped && (b[k] == -want); + } + } + return same ? 0 : (flipped ? 1 : -1); +} + +o2::field::MagneticField* loadFieldFromFile(const std::string& path) +{ + TFile* file = TFile::Open(path.c_str()); + if (file == nullptr || file->IsZombie()) { + progress("error: cannot open field file " + path); + return nullptr; + } + auto* field = dynamic_cast(file->Get(kFieldObjectKey)); + auto* probes = dynamic_cast(file->Get(kFieldProbeKey)); + if (field == nullptr) { + progress(form("error: no '%s' object in %s", kFieldObjectKey, path.c_str())); + return nullptr; + } + if (probes == nullptr) { + progress(form( + "error: no '%s' in %s -- the field cannot be verified against what was written, " + "and a silently inverted field is exactly what this check exists to catch. " + "Use --field-current instead.", + kFieldProbeKey, path.c_str())); + return nullptr; + } + const TVectorD reference(*probes); + file->Close(); + delete file; + + // Reload the parameterisation from this file, whatever path was stored at write + // time, so that the file is self-contained and relocatable. + field->setDataFileName(path.c_str()); + field->CreateField(); + + int comparison = compareToProbes(field, reference); + if (comparison == 1) { + field->setFactorSolenoid(-field->getFactorSolenoid()); + field->setFactorDipole(-field->getFactorDipole()); + comparison = compareToProbes(field, reference); + if (comparison == 0) { + progress("field: polarity flip from the non-idempotent CreateField() detected and repaired"); + } + } + if (comparison != 0) { + progress(form( + "error: the field reloaded from %s does not reproduce its own reference probes; " + "refusing to hand back a field that is not the one written", + path.c_str())); + return nullptr; + } + progress(form("field: %s verified against %d reference probe vectors", path.c_str(), reference.GetNrows() / 6)); + return field; +} + +double fieldMag(o2::field::MagneticField* field, double x, double y, double z) +{ + double point[3] = {x, y, z}; + double b[3] = {0., 0., 0.}; + field->Field(point, b); + return std::sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]); +} + +double fieldMagCyl(o2::field::MagneticField* field, double r, double phi, double z) +{ + return fieldMag(field, r * std::cos(phi), r * std::sin(phi), z); +} + +// --------------------------------------------------------------------------- +// the field-support model +// --------------------------------------------------------------------------- + +// The model states, per threshold, a list of z-bands each carrying the radial +// intervals in which |B| exceeds that threshold, maximised over phi. It is an +// OUTER bound: a point outside every band has |B| <= threshold. That direction is +// what lets a sampled quantity support a geometric argument about a volume. +// +// Two features of the real field dictate the sampling, and both were found by +// this model getting them wrong first: +// +// * the LHC machine elements are hard cylinders with a discontinuous edge (the +// A-side compensator aperture is exactly r < 4.0 cm), so every threshold +// crossing is bisected rather than left on the grid; +// +// * the measured map's coverage is a BOX in (x, y), not a cylinder. Between the +// box's inscribed and corner radius the field survives only in ~2 degree wedges +// at the four corners -- at r = 194.5, z = -797.9 it is 8.2 kG at phi = 132.8 +// degrees and exactly zero at 127.5 and 135. A model sampling 16 phi values +// steps straight over those wedges and declares 8 kG of dipole field +// unsupported. phi sampling is therefore bounded by ARC LENGTH, so the angular +// resolution follows the feature size at every radius. +// +// Interval edges are stored as best estimates with a separate edge uncertainty +// rather than pre-inflated, because the study this tool comes from turns on a +// volume whose inner radius is exactly the field boundary: it is separated from +// the field by exactly zero, and no amount of inflation may promote that to a +// margin. + +constexpr int kMinPhiSamples = 24; +constexpr double kPhiArcStep = 3.0; // cm, the azimuthal sampling bound +constexpr double kBisectionTol = 0.01; // cm, also the model's edge uncertainty +constexpr double kZStepCoarse = 1.0; // cm +constexpr double kZStepRefine = 0.1; // cm, used wherever the radial structure changes +constexpr double kScanRMax = 2100.; // cm, outer reach of the scan +constexpr double kScanRMaxFine = 900.; // cm, beyond this the radial grid is coarse +constexpr double kTightMargin = 0.05; // cm, for boundaries that are analytically hard +constexpr long kViolationScanPoints = 400000; + +struct Interval { + double lo, hi; +}; + +struct Band { + double zlo, zhi; + std::vector iv; +}; + +struct Model { + double thresholdKG = 0.; + std::vector bands; +}; + +struct Support { + std::vector models; + double edgeUncertainty = kBisectionTol; + double marginStrict = 5.0; + double marginTight = kTightMargin; + double zmin = 0., zmax = 0., rmax = 0.; + std::string parameterisation; + + /// Separation in the (z,r) half-plane between a placement's extent and the + /// support; 0 means they intersect. The support is a union over phi, so the + /// nearest support point is reachable at the volume's own phi and this planar + /// distance is the true 3-D distance. + double separation(int t, double vzmin, double vzmax, double vrmin, double vrmax) const + { + double best = 1e30; + for (const auto& band : models[t].bands) { + const double dz = std::max(0., std::max(band.zlo - vzmax, vzmin - band.zhi)); + for (const auto& iv : band.iv) { + const double lo = iv.lo - edgeUncertainty; + const double hi = iv.hi + edgeUncertainty; + const double dr = std::max(0., std::max(lo - vrmax, vrmin - hi)); + best = std::min(best, std::sqrt(dz * dz + dr * dr)); + if (best <= 0.) { + return 0.; + } + } + } + return (best > 1e29) ? 1e30 : best; + } + + bool supportAt(int t, double z, double r) const { return separation(t, z, z, r, r) <= 0.; } + + /// How deeply the extent penetrates the expanded support, 0 if it does not. + /// This separates "overlaps by less than the model can resolve" -- a volume + /// sharing a boundary surface with the field and nothing else -- from + /// "genuinely reaches into the field". + double penetration(int t, double vzmin, double vzmax, double vrmin, double vrmax) const + { + double worst = 0.; + for (const auto& band : models[t].bands) { + if (band.zlo > vzmax || band.zhi < vzmin) { + continue; + } + for (const auto& iv : band.iv) { + const double lo = iv.lo - edgeUncertainty; + const double hi = iv.hi + edgeUncertainty; + if (lo > vrmax || hi < vrmin) { + continue; + } + worst = std::max(worst, std::min(hi - vrmin, vrmax - lo)); + } + } + return worst; + } + + /// Is the whole extent inside the support? Then the placement is in the field + /// everywhere and there is nothing left to disprove. This is what makes the tool + /// tractable: most of the ~10^5 passive placements sit deep inside the solenoid, + /// and sampling each of them to rediscover that would dominate the runtime for + /// no information. It is a coverage walk rather than a single-band lookup, + /// because bands are cut wherever the radial structure changes and a long volume + /// crosses many of them. + bool coveredBySupport(int t, double vzmin, double vzmax, double vrmin, double vrmax) const + { + std::vector> covering; + for (const auto& band : models[t].bands) { + if (band.zhi < vzmin || band.zlo > vzmax) { + continue; + } + for (const auto& iv : band.iv) { + if (iv.lo - edgeUncertainty <= vrmin && iv.hi + edgeUncertainty >= vrmax) { + covering.push_back({band.zlo, band.zhi}); + break; + } + } + } + std::sort(covering.begin(), covering.end()); + double frontier = vzmin; + for (const auto& segment : covering) { + if (segment.first > frontier + 1e-9) { + return false; + } + frontier = std::max(frontier, segment.second); + if (frontier >= vzmax) { + return true; + } + } + return frontier >= vzmax; + } + + /// Outside its domain the model makes no claim, and the answer is not OUT. + bool inDomain(double vzmin, double vzmax, double vrmax) const + { + return vzmin >= zmin && vzmax <= zmax && vrmax <= rmax; + } +}; + +int phiSamplesAt(double r) +{ + if (r <= 0.) { + return 1; + } + return std::max(kMinPhiSamples, (int)std::ceil(2 * M_PI * r / kPhiArcStep)); +} + +double maxFieldOverPhi(o2::field::MagneticField* field, double r, double z) +{ + if (r == 0.) { + return fieldMag(field, 0., 0., z); + } + const int n = phiSamplesAt(r); + double worst = 0.; + for (int i = 0; i < n; ++i) { + worst = std::max(worst, fieldMagCyl(field, r, 2 * M_PI * i / n, z)); + } + return worst; +} + +std::vector radialGrid() +{ + std::vector grid; + for (double r = 0.; r < 20.; r += 0.1) { + grid.push_back(r); + } + for (double r = 20.; r < 100.; r += 1.0) { + grid.push_back(r); + } + for (double r = 100.; r < 800.; r += 2.0) { + grid.push_back(r); + } + for (double r = 800.; r <= kScanRMaxFine; r += 10.0) { + grid.push_back(r); + } + // A coarse extension so that an unexpected far feature is not invisible by construction. + for (double r = kScanRMaxFine + 25.; r <= kScanRMax; r += 25.0) { + grid.push_back(r); + } + return grid; +} + +/// Crossing radius between rOut (|B| <= threshold) and rIn (|B| > threshold). +/// The two may be given in either radial order: assuming rIn < rOut turns every +/// falling edge into a grid midpoint, which understates the support radius, and +/// understating support is the one direction this model may not err in. +double bisectCrossing(o2::field::MagneticField* field, double z, double rOut, double rIn, double threshold) +{ + for (int i = 0; i < 60 && std::fabs(rOut - rIn) > kBisectionTol; ++i) { + const double middle = 0.5 * (rOut + rIn); + if (maxFieldOverPhi(field, middle, z) > threshold) { + rIn = middle; + } else { + rOut = middle; + } + } + return 0.5 * (rOut + rIn); +} + +struct Slice { + double z = 0., zlo = 0., zhi = 0.; + std::vector> iv; // one per threshold +}; + +Slice sliceAt(o2::field::MagneticField* field, double z, const std::vector& thresholds) +{ + Slice slice; + slice.z = z; + slice.iv.resize(thresholds.size()); + const std::vector grid = radialGrid(); + std::vector b(grid.size()); + for (size_t i = 0; i < grid.size(); ++i) { + b[i] = maxFieldOverPhi(field, grid[i], z); + } + for (size_t t = 0; t < thresholds.size(); ++t) { + bool open = false; + Interval current{0., 0.}; + for (size_t i = 0; i < grid.size(); ++i) { + const bool above = b[i] > thresholds[t]; + if (above && !open) { + current.lo = (i > 0) ? std::max(0., bisectCrossing(field, z, grid[i - 1], grid[i], thresholds[t])) : grid[i]; + open = true; + } else if (!above && open) { + current.hi = bisectCrossing(field, z, grid[i], grid[i - 1], thresholds[t]); + slice.iv[t].push_back(current); + open = false; + } + } + if (open) { + current.hi = grid.back(); + slice.iv[t].push_back(current); + } + } + return slice; +} + +void mergeIntervals(std::vector& into, const std::vector& from) +{ + into.insert(into.end(), from.begin(), from.end()); + if (into.empty()) { + return; + } + std::sort(into.begin(), into.end(), [](const Interval& a, const Interval& b) { return a.lo < b.lo; }); + std::vector merged{into.front()}; + for (size_t i = 1; i < into.size(); ++i) { + if (into[i].lo <= merged.back().hi + 1e-9) { + merged.back().hi = std::max(merged.back().hi, into[i].hi); + } else { + merged.push_back(into[i]); + } + } + into.swap(merged); +} + +bool sameStructure(const std::vector& a, const std::vector& b) +{ + if (a.size() != b.size()) { + return false; + } + for (size_t i = 0; i < a.size(); ++i) { + if (std::fabs(a[i].lo - b[i].lo) > 0.02 || std::fabs(a[i].hi - b[i].hi) > 0.02) { + return false; + } + } + return true; +} + +Support buildSupport(o2::field::MagneticField* field, const std::vector& thresholds, double zmin, double zmax) +{ + progress(form( + "support: scanning z %.0f..%.0f, dz %.1f cm refined to %.1f, r to %.0f cm, " + "phi by arc length <= %.1f cm (%d samples at r=200)", + zmin, zmax, kZStepCoarse, kZStepRefine, kScanRMax, kPhiArcStep, phiSamplesAt(200.))); + + std::vector slices; + Slice previous; + bool havePrevious = false; + for (double z = zmin; z <= zmax + 1e-9; z += kZStepCoarse) { + Slice slice = sliceAt(field, z, thresholds); + bool changed = false; + for (size_t t = 0; havePrevious && t < thresholds.size(); ++t) { + changed = changed || !sameStructure(previous.iv[t], slice.iv[t]); + } + if (changed) { + for (double zz = previous.z + kZStepRefine; zz < z - 1e-9; zz += kZStepRefine) { + slices.push_back(sliceAt(field, zz, thresholds)); + } + } + slices.push_back(slice); + previous = slice; + havePrevious = true; + if (std::fmod(z - zmin, 500.) < kZStepCoarse / 2) { + progress(form("support: ... z = %.0f", z)); + } + } + + for (size_t i = 0; i < slices.size(); ++i) { + const double zPrev = (i == 0) ? slices[i].z - kZStepCoarse : slices[i - 1].z; + const double zNext = (i + 1 == slices.size()) ? slices[i].z + kZStepCoarse : slices[i + 1].z; + slices[i].zlo = 0.5 * (zPrev + slices[i].z); + slices[i].zhi = 0.5 * (slices[i].z + zNext); + } + + // A band carries the union of the intervals of ITS OWN slices, extended by half + // a coarse step at each end so that neighbouring bands overlap and a point + // between two samples is claimed by both. Unioning a band with its bracketing + // slices instead is sound only for a band one step wide: applied to a merged + // band it smears a neighbour's support across the whole length, and the field- + // free windows this tool exists to find are precisely empty bands between two + // populated ones. Conservatism has to stay local or it destroys them. + Support support; + support.models.resize(thresholds.size()); + for (size_t t = 0; t < thresholds.size(); ++t) { + support.models[t].thresholdKG = thresholds[t]; + size_t i = 0; + while (i < slices.size()) { + size_t j = i; + while (j + 1 < slices.size() && sameStructure(slices[j].iv[t], slices[j + 1].iv[t])) { + ++j; + } + Band band; + band.zlo = slices[i].zlo - 0.5 * kZStepCoarse; + band.zhi = slices[j].zhi + 0.5 * kZStepCoarse; + for (size_t k = i; k <= j; ++k) { + mergeIntervals(band.iv, slices[k].iv[t]); + } + if (!band.iv.empty()) { + support.models[t].bands.push_back(band); + } + i = j + 1; + } + } + support.zmin = zmin; + support.zmax = zmax; + support.rmax = kScanRMax; + support.parameterisation = field->getParameterName(); + return support; +} + +/// The model claims to be an outer bound. Test that against the field itself: +/// where the model reports no support, the real field must not exceed the +/// threshold. This is sampling used in the one direction it is allowed, to +/// disprove, and it is the only check that catches a model which is merely +/// plausible. A violated bound invalidates every OUT verdict downstream, so it is +/// fatal rather than a warning. +bool violationScan(o2::field::MagneticField* field, const Support& support, Report& report) +{ + TRandom3 random(10001); + std::vector violations(support.models.size(), 0); + std::vector worst(support.models.size(), 0.); + std::vector worstR(support.models.size(), 0.); + std::vector worstZ(support.models.size(), 0.); + for (long i = 0; i < kViolationScanPoints; ++i) { + // Half the points anywhere in the domain, half near the axis where the + // machine elements are narrow enough to hide between uniform samples. + const bool nearAxis = (i % 2 == 1); + const double z = nearAxis ? random.Uniform(std::max(support.zmin, -2200.), std::min(support.zmax, 2200.)) + : random.Uniform(support.zmin, support.zmax); + const double r = random.Uniform(0., nearAxis ? 20. : support.rmax); + const double b = fieldMagCyl(field, r, random.Uniform(0., 2 * M_PI), z); + for (size_t t = 0; t < support.models.size(); ++t) { + if (b > support.models[t].thresholdKG && !support.supportAt(t, z, r)) { + ++violations[t]; + if (b > worst[t]) { + worst[t] = b; + worstR[t] = r; + worstZ[t] = z; + } + } + } + } + bool ok = true; + for (size_t t = 0; t < support.models.size(); ++t) { + report(form(" outer bound at %6.1f G: %ld / %ld sampled points with field outside every band%s", + support.models[t].thresholdKG * 1000., violations[t], kViolationScanPoints, + violations[t] == 0 ? " (bound holds)" : " <-- THE MODEL IS NOT AN OUTER BOUND")); + if (violations[t] != 0) { + report(form(" worst: |B| = %.4f kG at r = %.3f, z = %.3f", worst[t], worstR[t], worstZ[t])); + ok = false; + } + } + return ok; +} + +json supportToJson(const Support& support, const std::string& fieldSource) +{ + const std::time_t now = std::time(nullptr); + char stamp[64]; + std::strftime(stamp, sizeof(stamp), "%Y-%m-%dT%H:%M:%S", std::gmtime(&now)); + + json out; + out["schema"] = "o2-sim-geometry-doctor/field_support/1"; + out["generated_utc"] = stamp; + out["field_source"] = fieldSource; + out["parameterisation"] = support.parameterisation; + out["units"] = "kGauss, cm"; + out["semantics"] = + "Outer bound on the support of |B|, maximised over phi. A point outside every band, " + "after expanding intervals by edge_uncertainty_cm, has |B| <= threshold."; + out["resolution"] = {{"dz_coarse", kZStepCoarse}, + {"dz_refine", kZStepRefine}, + {"dr_near_axis", 0.1}, + {"phi_arc_step_cm", kPhiArcStep}, + {"phi_min_samples", kMinPhiSamples}, + {"bisection_tol_cm", kBisectionTol}, + {"edge_uncertainty_cm", support.edgeUncertainty}}; + out["domain"] = {{"zmin", support.zmin}, {"zmax", support.zmax}, {"rmax", support.rmax}}; + out["recommended_margins_cm"] = {{"strict", support.marginStrict}, {"tight", support.marginTight}}; + out["models"] = json::array(); + for (const auto& model : support.models) { + json m; + m["threshold_kG"] = model.thresholdKG; + m["threshold_gauss"] = model.thresholdKG * 1000.; + m["n_bands"] = model.bands.size(); + m["bands"] = json::array(); + for (const auto& band : model.bands) { + json b; + b["zlo"] = band.zlo; + b["zhi"] = band.zhi; + b["iv"] = json::array(); + for (const auto& iv : band.iv) { + b["iv"].push_back(json::array({iv.lo, iv.hi})); + } + m["bands"].push_back(b); + } + out["models"].push_back(m); + } + return out; +} + +/// An empty or truncated model would declare the whole geometry field-free, which +/// is the worst direction to fail in, so every problem here is fatal. +bool supportFromJson(const json& in, Support& support) +{ + try { + support.edgeUncertainty = in.at("resolution").at("edge_uncertainty_cm").get(); + support.marginStrict = in.at("recommended_margins_cm").at("strict").get(); + support.marginTight = in.at("recommended_margins_cm").at("tight").get(); + support.zmin = in.at("domain").at("zmin").get(); + support.zmax = in.at("domain").at("zmax").get(); + support.rmax = in.at("domain").at("rmax").get(); + support.parameterisation = in.value("parameterisation", std::string()); + for (const auto& m : in.at("models")) { + Model model; + model.thresholdKG = m.at("threshold_kG").get(); + for (const auto& b : m.at("bands")) { + Band band; + band.zlo = b.at("zlo").get(); + band.zhi = b.at("zhi").get(); + for (const auto& iv : b.at("iv")) { + band.iv.push_back({iv.at(0).get(), iv.at(1).get()}); + } + model.bands.push_back(band); + } + support.models.push_back(model); + } + } catch (const std::exception& e) { + progress(std::string("error: cannot read the support model: ") + e.what()); + return false; + } + if (support.models.empty() || support.models.front().bands.empty()) { + progress( + "error: the support model is empty -- refusing to continue, since an empty model " + "would declare the whole geometry field-free"); + return false; + } + return true; +} + +// --------------------------------------------------------------------------- +// the placement table +// --------------------------------------------------------------------------- + +constexpr int kMaxDepth = 14; +constexpr size_t kMaxRows = 400000; +constexpr double kSampleDr = 0.05; // cm, radial step of the disproof scan +constexpr double kSampleArc = 0.5; // cm, azimuthal step of the disproof scan +constexpr double kSampleDzMax = 2.0; // cm +constexpr long kMaxSamplesPerRow = 4000000; + +struct Row { + std::string path, lv, medium, mother, shape; + std::string effectiveMother; ///< nearest non-assembly ancestor: whose material really surrounds this + std::string verdict = "UNCLASSIFIED"; + int ifield = -1; + int copyNo = 0, nDaughters = 0, depth = 0; + bool sensitive = false, assembly = false, approximateExtent = false, resolved = true; + double zmin = 0., zmax = 0., rmin = 0., rmax = 0.; + double separation = -1., penetration = 0.; + double maxB = -1., minB = -1.; + long nSampled = 0; + bool wholeVolumeSampled = false; + double wholeMaxB = -1., wholeMinB = -1.; + TGeoNode* node = nullptr; + TGeoHMatrix matrix; +}; + +bool isOutFamily(const std::string& verdict) +{ + return verdict == "OUT" || verdict == "OUT_TIGHT" || verdict == "OUT_BOUNDARY"; +} + +/// Everything that is not an OUT-family verdict counts as in-field, including +/// UNKNOWN and OUTSIDE_DOMAIN. Defaulting the unclear cases to "has field" is the +/// only safe direction: the cost of being wrong the other way is straight-line +/// transport through a real field. +bool isInFamily(const std::string& verdict) +{ + return verdict == "IN" || verdict == "IN_COVERED" || verdict == "UNKNOWN" || verdict == "OUTSIDE_DOMAIN"; +} + +/// Radial extent from the SHAPE, not from bounding-box corners: the corners of a +/// box around an on-axis tube all sit at the same radius and report rmin == rmax. +/// Two shapes need care beyond that. TGeoPgon's Rmin/Rmax are apothems, so the +/// circumscribed radius is Rmax / cos(pi/n). And TGeoEltu MUST be tested before +/// TGeoTube: it derives from TGeoTube and reuses fRmin/fRmax to hold its two +/// semi-axes, so the tube branch makes a solid elliptical pipe report an inner +/// radius of 3.175 cm, as though it avoided the beam axis it plainly contains. +/// That is the dangerous direction -- it promotes an in-field volume to field-free. +bool shapeRadii(TGeoShape* shape, double& rmin, double& rmax) +{ + if (auto* pgon = dynamic_cast(shape)) { + rmin = 1e30; + rmax = -1e30; + for (int i = 0; i < pgon->GetNz(); ++i) { + rmin = std::min(rmin, pgon->GetRmin(i)); + rmax = std::max(rmax, pgon->GetRmax(i)); + } + const double edges = pgon->GetNedges() > 2 ? pgon->GetNedges() : 3; + rmax /= std::cos(M_PI / edges); + return true; + } + if (auto* pcon = dynamic_cast(shape)) { + rmin = 1e30; + rmax = -1e30; + for (int i = 0; i < pcon->GetNz(); ++i) { + rmin = std::min(rmin, pcon->GetRmin(i)); + rmax = std::max(rmax, pcon->GetRmax(i)); + } + return true; + } + if (auto* cone = dynamic_cast(shape)) { + rmin = std::min(cone->GetRmin1(), cone->GetRmin2()); + rmax = std::max(cone->GetRmax1(), cone->GetRmax2()); + return true; + } + if (auto* eltu = dynamic_cast(shape)) { + rmin = 0.; // solid, so it contains the axis + rmax = std::max(eltu->GetA(), eltu->GetB()); + return true; + } + if (auto* tube = dynamic_cast(shape)) { + rmin = tube->GetRmin(); + rmax = tube->GetRmax(); + return true; + } + return false; +} + +bool zPreserving(const TGeoHMatrix& m) +{ + const Double_t* r = m.GetRotationMatrix(); + return std::fabs(std::fabs(r[8]) - 1.) < 1e-9 && std::fabs(r[2]) < 1e-9 && std::fabs(r[5]) < 1e-9 && + std::fabs(r[6]) < 1e-9 && std::fabs(r[7]) < 1e-9; +} + +void extentOf(TGeoNode* node, const TGeoHMatrix& matrix, Row& row) +{ + TGeoShape* shape = node->GetVolume()->GetShape(); + auto* box = dynamic_cast(shape); + if (box == nullptr) { + // No bounding box means no analytic extent; claim everything, which classifies + // as in-field and never as OUT. + row.zmin = row.zmax = 0.; + row.rmin = 0.; + row.rmax = 1e9; + row.approximateExtent = true; + return; + } + + row.zmin = 1e30; + row.zmax = -1e30; + double boxRmax = 0.; + const double dx = box->GetDX(), dy = box->GetDY(), dz = box->GetDZ(); + const double* origin = box->GetOrigin(); + for (int i = 0; i < 8; ++i) { + double local[3] = {origin[0] + ((i & 1) ? dx : -dx), origin[1] + ((i & 2) ? dy : -dy), + origin[2] + ((i & 4) ? dz : -dz)}; + double global[3]; + matrix.LocalToMaster(local, global); + row.zmin = std::min(row.zmin, global[2]); + row.zmax = std::max(row.zmax, global[2]); + boxRmax = std::max(boxRmax, std::hypot(global[0], global[1])); + } + + const double* translation = matrix.GetTranslation(); + const double offAxis = std::hypot(translation[0] + origin[0], translation[1] + origin[1]); + double localRmin = 0., localRmax = 0.; + if (zPreserving(matrix) && shapeRadii(shape, localRmin, localRmax)) { + row.rmin = std::max(0., localRmin - offAxis); + row.rmax = localRmax + offAxis; + row.approximateExtent = false; + } else { + row.rmax = boxRmax; + row.rmin = (offAxis <= std::hypot(dx, dy)) ? 0. : std::max(0., offAxis - std::hypot(dx, dy)); + row.approximateExtent = true; + } +} + +// --------------------------------------------------------------------------- +// the doctor +// --------------------------------------------------------------------------- + +struct ContainerProposal { + std::string mother, motherPath; + double zlo, zhi, rmax; + int nDaughters; + bool clearedByStrictMargin; + bool sensitive; +}; + +struct SharedVolume { + std::string lv; + std::vector out, in; +}; + +class Doctor +{ + public: + Doctor(o2::field::MagneticField* field, const Support& support) : mField(field), mSupport(support) + { + mThreshold = support.models.front().thresholdKG; + } + + void walk(TGeoNode* node) { walk(node, nullptr, TGeoHMatrix(), 0, "", ""); } + void classifyAll(); + void findFindings(); + + const std::vector& rows() const { return mRows; } + size_t nPruned() const { return mPruned; } + const std::vector& reverseAudit() const { return mReverse; } + const std::vector& sharedVolumes() const { return mShared; } + const std::vector& straddlingMothers() const { return mStraddling; } + const std::vector& containers() const { return mContainers; } + bool hasSensitive(TGeoVolume* volume); + + private: + void walk(TGeoNode* node, TGeoVolume* mother, const TGeoHMatrix& parent, int depth, const std::string& path, + const std::string& effectiveMother); + void classify(Row& row); + void disproofScan(Row& row); + void wholeVolumeScan(Row& row); + bool ownMaterialAt(const Row& row, const double* global) const; + + o2::field::MagneticField* mField; + const Support& mSupport; + double mThreshold; + std::vector mRows; + std::map mSensitiveCache; + size_t mPruned = 0; + + std::vector mReverse; + std::vector mShared; + std::vector mStraddling; + std::vector mContainers; +}; + +bool Doctor::hasSensitive(TGeoVolume* volume) +{ + auto cached = mSensitiveCache.find(volume); + if (cached != mSensitiveCache.end() && cached->second >= 0) { + return cached->second == 1; + } + mSensitiveCache[volume] = 0; // guard against re-entry + auto* medium = volume->GetMedium(); + bool found = medium != nullptr && medium->GetParam(0) != 0.; + for (int i = 0; i < volume->GetNdaughters() && !found; ++i) { + found = hasSensitive(volume->GetNode(i)->GetVolume()); + } + mSensitiveCache[volume] = found ? 1 : 0; + return found; +} + +/// A volume is structural if it is an assembly or made of its mother's medium +/// (barrel inside cave, caveRB24 inside cave). The walk descends through those and +/// through fully passive subtrees, and prunes detector modules: this tool is about +/// the passive geometry, and a proposal touching a sensitive subtree is refused by +/// default anyway, because sensitive volumes are resolved by name after +/// construction and re-parenting one silently produces zero hits. +bool isStructural(TGeoNode* node, TGeoVolume* mother) +{ + if (mother == nullptr) { + return true; + } + TGeoVolume* volume = node->GetVolume(); + if (volume->IsAssembly()) { + return true; + } + auto* mine = volume->GetMedium(); + auto* theirs = mother->GetMedium(); + return mine != nullptr && theirs != nullptr && std::strcmp(mine->GetName(), theirs->GetName()) == 0; +} + +void Doctor::walk(TGeoNode* node, TGeoVolume* mother, const TGeoHMatrix& parent, int depth, const std::string& path, + const std::string& effectiveMother) +{ + if (mRows.size() >= kMaxRows || depth > kMaxDepth) { + return; + } + TGeoHMatrix here = parent; + here.Multiply(node->GetMatrix()); + TGeoVolume* volume = node->GetVolume(); + const std::string myPath = path + "/" + volume->GetName() + "_" + std::to_string(node->GetNumber()); + + Row row; + row.path = myPath; + row.lv = volume->GetName(); + row.mother = mother != nullptr ? mother->GetName() : ""; + row.effectiveMother = effectiveMother; + row.shape = volume->GetShape()->ClassName(); + row.copyNo = node->GetNumber(); + row.nDaughters = volume->GetNdaughters(); + row.depth = depth; + row.assembly = volume->IsAssembly(); + auto* medium = volume->GetMedium(); + row.medium = medium != nullptr ? medium->GetName() : "(none)"; + row.ifield = medium != nullptr ? (int)medium->GetParam(1) : -1; + row.sensitive = medium != nullptr && medium->GetParam(0) != 0.; + row.node = node; + row.matrix = here; + extentOf(node, here, row); + mRows.push_back(row); + + if (hasSensitive(volume) && !isStructural(node, mother)) { + ++mPruned; + return; + } + for (int i = 0; i < node->GetNdaughters(); ++i) { + walk(node->GetDaughter(i), volume, here, depth + 1, myPath, row.assembly ? effectiveMother : myPath); + } +} + +// A mother's shape includes the space its daughters occupy, but its MATERIAL is +// only what they leave over. That distinction is the whole point of these scans: +// caveRB24's shape reaches the beam axis, and what its daughters leave over is a +// sliver of cave air at r ~ 3-4 cm, between an oval beam pipe and the circular +// field cylinder, carrying 13.2 kG. +bool insideAnyDaughter(TGeoVolume* volume, const double* local) +{ + for (int i = 0; i < volume->GetNdaughters(); ++i) { + TGeoNode* daughter = volume->GetNode(i); + double inDaughter[3]; + daughter->GetMatrix()->MasterToLocal(local, inDaughter); + if (!daughter->GetVolume()->GetShape()->Contains(inDaughter)) { + continue; + } + if (daughter->GetVolume()->IsAssembly()) { + if (insideAnyDaughter(daughter->GetVolume(), inDaughter)) { + return true; + } + } else { + return true; + } + } + return false; +} + +bool Doctor::ownMaterialAt(const Row& row, const double* global) const +{ + double local[3]; + row.matrix.MasterToLocal(global, local); + if (!row.node->GetVolume()->GetShape()->Contains(local)) { + return false; + } + return !insideAnyDaughter(row.node->GetVolume(), local); +} + +double thinnestDaughter(TGeoVolume* volume) +{ + double thinnest = 1e30; + for (int i = 0; i < volume->GetNdaughters(); ++i) { + auto* box = dynamic_cast(volume->GetNode(i)->GetVolume()->GetShape()); + if (box != nullptr) { + thinnest = std::min(thinnest, 2 * std::min(box->GetDX(), std::min(box->GetDY(), box->GetDZ()))); + } + } + return thinnest; +} + +/// Sampling restricted to where an answer can be: the intersection of the +/// placement's extent with the support bands. Everywhere else the field is below +/// threshold by construction, so sampling it would cost time and tell nothing. +void Doctor::disproofScan(Row& row) +{ + const double thinnest = thinnestDaughter(row.node->GetVolume()); + row.resolved = (row.nDaughters == 0) || (kSampleDr <= std::max(0.1, thinnest)); + long budget = kMaxSamplesPerRow; + for (const auto& band : mSupport.models.front().bands) { + const double z0 = std::max(row.zmin, band.zlo); + const double z1 = std::min(row.zmax, band.zhi); + if (z1 < z0) { + continue; + } + for (const auto& iv : band.iv) { + const double r0 = std::max(row.rmin, iv.lo - mSupport.edgeUncertainty); + const double r1 = std::min(row.rmax, iv.hi + mSupport.edgeUncertainty); + if (r1 < r0) { + continue; + } + const double dz = std::min(kSampleDzMax, std::max(0.25, (z1 - z0) / 200.)); + for (double z = z0; z <= z1 + 1e-9 && budget > 0; z += dz) { + for (double r = r0; r <= r1 + 1e-9 && budget > 0; r += kSampleDr) { + const int nphi = (r <= 0.) ? 1 : std::max(8, (int)std::ceil(2 * M_PI * r / kSampleArc)); + for (int i = 0; i < nphi && budget > 0; ++i) { + const double phi = 2 * M_PI * i / nphi; + double global[3] = {r * std::cos(phi), r * std::sin(phi), z}; + --budget; + if (!ownMaterialAt(row, global)) { + continue; + } + ++row.nSampled; + const double b = fieldMag(mField, global[0], global[1], global[2]); + row.maxB = std::max(row.maxB, b); + row.minB = (row.minB < 0.) ? b : std::min(row.minB, b); + } + } + } + } + } + if (budget <= 0) { + row.resolved = false; + } +} + +/// A coarse scan over the whole placement, to see whether its own material spans +/// both sides of the predicate. Cheap, and it is what makes "the mother straddles +/// the boundary" a measurement rather than an assertion. +void Doctor::wholeVolumeScan(Row& row) +{ + row.wholeVolumeSampled = true; + const double dz = std::max(1.0, (row.zmax - row.zmin) / 300.); + const double dr = std::max(0.5, (row.rmax - row.rmin) / 200.); + long budget = 2000000; + for (double z = row.zmin; z <= row.zmax && budget > 0; z += dz) { + for (double r = row.rmin; r <= row.rmax && budget > 0; r += dr) { + const int nphi = (r <= 0.) ? 1 : std::max(8, (int)std::ceil(2 * M_PI * r / std::max(2.0, dr))); + for (int i = 0; i < nphi && budget > 0; ++i) { + const double phi = 2 * M_PI * i / nphi; + double global[3] = {r * std::cos(phi), r * std::sin(phi), z}; + --budget; + if (!ownMaterialAt(row, global)) { + continue; + } + const double b = fieldMag(mField, global[0], global[1], global[2]); + row.wholeMaxB = std::max(row.wholeMaxB, b); + row.wholeMinB = (row.wholeMinB < 0.) ? b : std::min(row.wholeMinB, b); + } + } + } +} + +void Doctor::classify(Row& row) +{ + if (row.assembly) { + row.verdict = "ASSEMBLY"; // virtual: no material, so nothing to flag + return; + } + if (!mSupport.inDomain(row.zmin, row.zmax, row.rmax)) { + row.verdict = "OUTSIDE_DOMAIN"; + return; + } + row.separation = mSupport.separation(0, row.zmin, row.zmax, row.rmin, row.rmax); + if (row.separation >= mSupport.marginStrict) { + row.verdict = "OUT"; + return; + } + if (row.separation >= mSupport.marginTight) { + row.verdict = "OUT_TIGHT"; + return; + } + row.penetration = mSupport.penetration(0, row.zmin, row.zmax, row.rmin, row.rmax); + if (mSupport.coveredBySupport(0, row.zmin, row.zmax, row.rmin, row.rmax)) { + row.verdict = "IN_COVERED"; // in the field by containment, nothing to disprove + return; + } + disproofScan(row); + if (row.maxB > mThreshold) { + row.verdict = "IN"; + } else if (row.separation > 0. || row.penetration <= 2 * mSupport.edgeUncertainty) { + // Touches the support only to within what the model can resolve. Field-freeness + // then rests on an analytic claim about a hard boundary, so it is emitted for a + // human to confirm and never treated as established. + row.verdict = "OUT_BOUNDARY"; + } else { + row.verdict = "UNKNOWN"; + } +} + +void Doctor::classifyAll() +{ + size_t done = 0; + for (auto& row : mRows) { + classify(row); + if (++done % 50000 == 0) { + progress(form("classify: %zu / %zu placements", done, mRows.size())); + } + } + for (auto& row : mRows) { + // Only a straddler can be heterogeneous: a placement wholly inside the support + // is uniformly in the field, one that is OUT is uniformly out of it. + if (row.nDaughters > 0 && !row.assembly && + (row.verdict == "IN" || row.verdict == "UNKNOWN" || row.verdict == "OUT_BOUNDARY")) { + wholeVolumeScan(row); + } + } +} + +void Doctor::findFindings() +{ + // --- the reverse audit: ifield == 0 media that reach into the field --------- + // An ifield == 0 medium tells the transport engine to move in a straight line. + // Where that is wrong it is a physics bug, so the finding must carry the field + // actually present in the volume's own material, measured. Rows short-circuited + // as IN_COVERED were never sampled, so they are sampled now. + for (auto& row : mRows) { + if (row.assembly || row.ifield != 0 || row.medium == "dummy" || row.medium == "(none)") { + continue; + } + if (isInFamily(row.verdict)) { + if (!row.wholeVolumeSampled) { + wholeVolumeScan(row); + } + mReverse.push_back(&row); + } + } + + // --- shared logical volumes placed on both sides (signature C) -------------- + std::map> byVolume; + for (const auto& row : mRows) { + if (!row.assembly) { + byVolume[row.lv].push_back(&row); + } + } + for (const auto& entry : byVolume) { + SharedVolume shared; + shared.lv = entry.first; + for (const auto* row : entry.second) { + if (isOutFamily(row->verdict)) { + shared.out.push_back(row); + } else if (isInFamily(row->verdict)) { + shared.in.push_back(row); + } + } + if (!shared.out.empty() && !shared.in.empty()) { + mShared.push_back(shared); + } + } + std::sort(mShared.begin(), mShared.end(), + [](const SharedVolume& a, const SharedVolume& b) { return a.out.size() > b.out.size(); }); + + // --- mothers whose own material straddles the predicate (signature D) ------- + for (const auto& row : mRows) { + if (row.nDaughters > 0 && !row.assembly && row.wholeVolumeSampled && row.wholeMaxB > mThreshold && + row.wholeMinB <= mThreshold) { + mStraddling.push_back(&row); + } + } + + // --- missing containers ----------------------------------------------------- + // Where a container belongs is not where the daughters have a spatial gap: a + // beam pipe is a continuous chain of volumes with no gap wider than a flange, so + // gap-based clustering returns one cluster spanning everything and proposes + // nothing. What separates the daughters is the PREDICATE. A maximal run of + // consecutive field-free daughters is exactly a group that cannot be given its + // own medium today, because the only thing enclosing it is the mother's air, + // which is not field-free along its whole length. + for (const auto* mother : mStraddling) { + std::vector kids; + for (const auto& row : mRows) { + if (row.effectiveMother == mother->path && !row.assembly) { + kids.push_back(&row); + } + } + std::sort(kids.begin(), kids.end(), [](const Row* a, const Row* b) { return a->zmin < b->zmin; }); + + size_t i = 0; + while (i < kids.size()) { + if (!isOutFamily(kids[i]->verdict)) { + ++i; + continue; + } + size_t j = i; + double zlo = kids[i]->zmin, zhi = kids[i]->zmax, rmax = kids[i]->rmax; + while (j + 1 < kids.size() && isOutFamily(kids[j + 1]->verdict)) { + ++j; + zlo = std::min(zlo, kids[j]->zmin); + zhi = std::max(zhi, kids[j]->zmax); + rmax = std::max(rmax, kids[j]->rmax); + } + const int n = (int)(j - i + 1); + + // A container must not swallow anything outside its run. Rejecting any run + // with an intruder is too blunt: a 0.05 cm overhang from a neighbouring pipe + // section would kill a 51-daughter proposal. The right answer to a 0.05 cm + // overhang is to move the container's edge, so the edges are clamped past any + // overhanging neighbour and only a residual overlap makes the run unusable. + double clampedLo = zlo, clampedHi = zhi; + std::vector intruders; + for (size_t k = 0; k < kids.size(); ++k) { + if (k >= i && k <= j) { + continue; + } + const Row* other = kids[k]; + if (other->zmax <= zlo + 0.01 || other->zmin >= zhi - 0.01 || other->rmin >= rmax - 0.01) { + continue; + } + intruders.push_back(other); + } + for (const Row* other : intruders) { + if (other->zmin <= clampedLo + 1e-9 && other->zmax > clampedLo) { + clampedLo = other->zmax; + } + if (other->zmax >= clampedHi - 1e-9 && other->zmin < clampedHi) { + clampedHi = other->zmin; + } + } + bool residual = false; + for (const Row* other : intruders) { + residual = residual || (other->zmax > clampedLo + 0.01 && other->zmin < clampedHi - 0.01 && + other->rmin < rmax - 0.01); + } + + const double separation = + (clampedHi > clampedLo) ? mSupport.separation(0, clampedLo, clampedHi, 0., rmax) : -1.; + if (n >= 2 && separation >= mSupport.marginTight && !residual) { + ContainerProposal proposal; + proposal.mother = mother->lv; + proposal.motherPath = mother->path; + proposal.zlo = clampedLo; + proposal.zhi = clampedHi; + proposal.rmax = rmax; + proposal.nDaughters = n; + proposal.clearedByStrictMargin = separation >= mSupport.marginStrict; + proposal.sensitive = false; + for (size_t k = i; k <= j; ++k) { + proposal.sensitive = proposal.sensitive || hasSensitive(kids[k]->node->GetVolume()); + } + mContainers.push_back(proposal); + } + i = j + 1; + } + } +} + +// --------------------------------------------------------------------------- +// outputs +// --------------------------------------------------------------------------- + +void writePlacementCsv(const std::vector& rows, const std::string& path) +{ + std::FILE* out = std::fopen(path.c_str(), "w"); + if (out == nullptr) { + progress("error: cannot write " + path); + return; + } + std::fprintf(out, + "path,lv,medium,ifield,shape,mother,copy,ndaughters,sensitive,assembly,approx," + "zmin,zmax,rmin,rmax,verdict,separation_cm,penetration_cm,maxB_kG,nsampled\n"); + for (const auto& row : rows) { + std::fprintf(out, "%s,%s,%s,%d,%s,%s,%d,%d,%d,%d,%d,%.3f,%.3f,%.3f,%.3f,%s,%.4f,%.4f,%.6f,%ld\n", row.path.c_str(), + row.lv.c_str(), row.medium.c_str(), row.ifield, row.shape.c_str(), row.mother.c_str(), row.copyNo, + row.nDaughters, (int)row.sensitive, (int)row.assembly, (int)row.approximateExtent, row.zmin, row.zmax, + row.rmin, row.rmax, row.verdict.c_str(), row.separation, row.penetration, row.maxB, row.nSampled); + } + std::fclose(out); +} + +json placementJson(const Row& row) +{ + return json{{"path", row.path}, + {"copy", row.copyNo}, + {"z", json::array({row.zmin, row.zmax})}, + {"r", json::array({row.rmin, row.rmax})}, + {"verdict", row.verdict}, + {"separation_cm", row.separation}, + {"max_B_kG", std::max(row.maxB, row.wholeMaxB)}}; +} + +json proposalsToJson(Doctor& doctor, const Support& support, const std::string& geometryFile, + const std::string& fieldSource) +{ + const std::time_t now = std::time(nullptr); + char stamp[64]; + std::strftime(stamp, sizeof(stamp), "%Y-%m-%dT%H:%M:%S", std::gmtime(&now)); + + json out; + out["schema"] = "o2-sim-geometry-doctor/proposals/1"; + out["generated_utc"] = stamp; + out["geometry"] = geometryFile; + out["field_source"] = fieldSource; + out["threshold_kG"] = support.models.front().thresholdKG; + out["margins_cm"] = {{"strict", support.marginStrict}, + {"tight", support.marginTight}, + {"edge_uncertainty", support.edgeUncertainty}}; + out["proposals"] = json::array(); + + for (const auto& shared : doctor.sharedVolumes()) { + bool refused = false; + for (const auto* row : shared.out) { + refused = refused || row->sensitive || doctor.hasSensitive(row->node->GetVolume()); + } + json entry; + entry["signature"] = "shared-volume"; + entry["action"] = "split the logical volume, so that its field-free placements can carry a field-free medium"; + entry["logical_volume"] = shared.lv; + entry["n_out"] = shared.out.size(); + entry["n_in"] = shared.in.size(); + entry["status"] = refused ? "refused by default (sensitive path)" : "proposed"; + entry["out_placements"] = json::array(); + for (const auto* row : shared.out) { + entry["out_placements"].push_back(placementJson(*row)); + } + entry["in_placements"] = json::array(); + for (const auto* row : shared.in) { + entry["in_placements"].push_back(placementJson(*row)); + } + out["proposals"].push_back(entry); + } + + for (const auto* row : doctor.reverseAudit()) { + const double maxB = std::max(row->maxB, row->wholeMaxB); + json entry; + entry["signature"] = "reverse-audit"; + entry["action"] = maxB > support.models.front().thresholdKG + ? "the field-free medium assignment is wrong: straight-line transport inside real field" + : "the field-free medium reaches field support but no field was found in its own material, review"; + entry["path"] = row->path; + entry["logical_volume"] = row->lv; + entry["medium"] = row->medium; + entry["copy"] = row->copyNo; + entry["z"] = json::array({row->zmin, row->zmax}); + entry["r"] = json::array({row->rmin, row->rmax}); + entry["min_B_kG"] = row->wholeMinB; + entry["max_B_kG"] = maxB; + entry["verdict"] = row->verdict; + out["proposals"].push_back(entry); + } + + for (const auto& container : doctor.containers()) { + json entry; + entry["signature"] = "missing-container"; + entry["action"] = "insert a container with a field-free medium and re-parent the cluster into it"; + entry["mother"] = container.mother; + entry["mother_path"] = container.motherPath; + entry["z"] = json::array({container.zlo, container.zhi}); + entry["rmax"] = container.rmax; + entry["n_daughters"] = container.nDaughters; + entry["status"] = container.sensitive ? "refused by default (sensitive path)" : "proposed"; + entry["clearance"] = container.clearedByStrictMargin ? "strict margin" : "tight margin"; + out["proposals"].push_back(entry); + } + + for (const auto* row : doctor.straddlingMothers()) { + json entry; + entry["signature"] = "heterogeneous-mother"; + entry["action"] = + "the mother's own material spans both sides of the predicate, so no per-medium flag can " + "express it; it needs a container"; + entry["path"] = row->path; + entry["logical_volume"] = row->lv; + entry["medium"] = row->medium; + entry["min_B_kG"] = row->wholeMinB; + entry["max_B_kG"] = row->wholeMaxB; + entry["n_daughters"] = row->nDaughters; + out["proposals"].push_back(entry); + } + return out; +} + +// --------------------------------------------------------------------------- +// the anchor self-check +// --------------------------------------------------------------------------- + +// A unit test for this tool would need a placed ALICE geometry and a field map, +// neither of which belongs in the repository, so the regression gate is instead a +// file of expectations that any ALICE Run 3 geometry must satisfy, checked against +// a real run with --verify-anchors. See run/geometry-doctor-anchors.json. +// +// Supported expectations: OUT (every placement of the volume is in the OUT +// family), IN, NOT_OUT, ASSEMBLY, and REVERSE_AUDIT_FLAGGED. An anchor may also +// require a placement count and a lower bound on the field found in the volume's +// own material. + +bool verifyAnchors(const std::string& path, Doctor& doctor, Report& report) +{ + std::ifstream in(path); + if (!in) { + report(" cannot open the anchor file " + path); + return false; + } + json anchors; + try { + in >> anchors; + } catch (const std::exception& e) { + report(std::string(" cannot parse the anchor file: ") + e.what()); + return false; + } + + std::set flaggedByReverseAudit; + for (const auto* row : doctor.reverseAudit()) { + flaggedByReverseAudit.insert(row->lv); + } + + bool allPassed = true; + report(form(" %-22s %-22s %-10s %s", "volume", "expected", "verdict", "evidence")); + for (const auto& anchor : anchors.at("anchors")) { + const auto volume = anchor.at("volume").get(); + const auto expected = anchor.at("expect").get(); + + int placements = 0; + int failures = 0; + double worstSeparation = 1e30; + double bestField = -1.; + std::string reported; + for (const auto& row : doctor.rows()) { + if (row.lv != volume) { + continue; + } + ++placements; + bool ok = true; + if (expected == "OUT") { + ok = isOutFamily(row.verdict) || row.verdict == "ASSEMBLY"; + } else if (expected == "IN") { + ok = row.verdict == "IN"; + } else if (expected == "NOT_OUT") { + ok = !isOutFamily(row.verdict); + } else if (expected == "ASSEMBLY") { + ok = row.verdict == "ASSEMBLY"; + } else if (expected == "REVERSE_AUDIT_FLAGGED") { + ok = true; // decided below, on the volume rather than the placement + } else { + report(" unknown expectation '" + expected + "' for " + volume); + return false; + } + failures += ok ? 0 : 1; + bestField = std::max(bestField, std::max(row.maxB, row.wholeMaxB)); + if (row.separation < worstSeparation) { + worstSeparation = row.separation; + reported = row.verdict; + } + } + + if (expected == "REVERSE_AUDIT_FLAGGED") { + failures = flaggedByReverseAudit.count(volume) > 0 ? 0 : 1; + reported = failures == 0 ? "flagged" : "not flagged"; + } + if (placements == 0) { + failures = 1; + reported = "not placed"; + } + if (anchor.contains("placements") && placements != anchor.at("placements").get()) { + failures += 1; + } + if (anchor.contains("min_max_B_kG") && bestField < anchor.at("min_max_B_kG").get()) { + failures += 1; + } + + allPassed = allPassed && failures == 0; + report(form(" %-22s %-22s %-10s %d placement(s), separation %.3f cm, max|B| %.4f kG %s", volume.c_str(), + expected.c_str(), reported.c_str(), placements, worstSeparation > 1e29 ? -1. : worstSeparation, + bestField, failures == 0 ? "PASS" : "FAIL")); + } + return allPassed; +} + +// --------------------------------------------------------------------------- + +struct Options { + std::string geometryFile; + std::string fieldFile; + int fieldCurrent = 0; + std::string supportFile; + std::string anchorFile; + std::vector thresholdsGauss; + double margin = 5.0; + std::string outputPrefix = "geometry-doctor"; +}; + +} // namespace + +int main(int argc, char** argv) +{ + Options options; + bpo::options_description description( + "Audits a placed geometry against the magnetic field it will be transported " + "in, and reports where the two do not fit together.\n\nOptions"); + description.add_options() // + ("help,h", "print this help message") // + ("geometry-file", bpo::value(&options.geometryFile)->required(), // + "the geometry to audit, e.g. o2sim_geometry.root") // + ("field-file", bpo::value(&options.fieldFile), // + "a serialized MagneticField carrying reference probe vectors") // + ("field-current", bpo::value(&options.fieldCurrent), // + "build the nominal field for this L3 current instead, e.g. -5") // + ("support-file", bpo::value(&options.supportFile), // + "field-support model cache: read it if it exists, otherwise write it") // + ("threshold", bpo::value>(&options.thresholdsGauss)->composing(), // + "field threshold in Gauss, repeatable; the lowest one decides the verdicts (default 1 and 10)") // + ("margin", bpo::value(&options.margin)->default_value(5.0), // + "clearance in cm a placement must keep from the field support to be called field-free") // + ("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"); + + bpo::variables_map arguments; + try { + bpo::store(bpo::parse_command_line(argc, argv, description), arguments); + if (arguments.count("help") != 0u) { + std::cout << description << '\n'; + return 0; + } + bpo::notify(arguments); + } catch (const bpo::error& e) { + std::cerr << "error: " << e.what() << "\n\n" + << description << '\n'; + return 1; + } + + const bool haveFieldFile = arguments.count("field-file") != 0u; + const bool haveFieldCurrent = arguments.count("field-current") != 0u; + if (haveFieldFile == haveFieldCurrent) { + std::cerr << "error: give exactly one of --field-file and --field-current\n"; + return 1; + } + if (options.thresholdsGauss.empty()) { + options.thresholdsGauss = {1., 10.}; + } + std::sort(options.thresholdsGauss.begin(), options.thresholdsGauss.end()); + std::vector thresholds; // kGauss, as the field itself reports + for (double gauss : options.thresholdsGauss) { + thresholds.push_back(gauss * 1e-3); + } + + const std::string fieldSource = + haveFieldFile ? options.fieldFile : form("createNominalField(%d)", options.fieldCurrent); + o2::field::MagneticField* field = + haveFieldFile ? loadFieldFromFile(options.fieldFile) : o2::field::MagneticField::createNominalField(options.fieldCurrent); + if (field == nullptr) { + std::cerr << "error: no usable magnetic field\n"; + return 1; + } + + Report report; + report("ALICE simulation geometry doctor"); + report(""); + report(" geometry : " + options.geometryFile); + report(" field : " + fieldSource + ", parameterisation " + field->getParameterName()); + + // --- the field-support model ------------------------------------------------ + Support support; + bool supportFromCache = false; + if (!options.supportFile.empty()) { + std::ifstream cache(options.supportFile); + if (cache) { + json cached; + try { + cache >> cached; + } catch (const std::exception& e) { + std::cerr << "error: cannot parse " << options.supportFile << ": " << e.what() << '\n'; + return 1; + } + if (!supportFromJson(cached, support)) { + return 1; + } + supportFromCache = true; + } + } + + if (supportFromCache) { + if (support.models.size() != thresholds.size()) { + std::cerr << "error: " << options.supportFile << " carries " << support.models.size() + << " thresholds but " << thresholds.size() << " were requested\n"; + return 1; + } + for (size_t t = 0; t < thresholds.size(); ++t) { + if (std::fabs(support.models[t].thresholdKG - thresholds[t]) > 1e-9) { + std::cerr << "error: " << options.supportFile << " was built for a different threshold (" + << support.models[t].thresholdKG * 1000. << " G against " << thresholds[t] * 1000. << " G)\n"; + return 1; + } + } + if (!support.parameterisation.empty() && support.parameterisation != field->getParameterName()) { + std::cerr << "error: " << options.supportFile << " was built for parameterisation " + << support.parameterisation << ", not " << field->getParameterName() << '\n'; + return 1; + } + report(" support model : " + options.supportFile + " (cached)"); + } else { + support = buildSupport(field, thresholds, -3000., 3000.); + if (!options.supportFile.empty()) { + std::ofstream out(options.supportFile); + out << supportToJson(support, fieldSource).dump(1, '\t') << '\n'; + report(" support model : built and written to " + options.supportFile); + } else { + report(" support model : built for this run"); + } + } + support.marginStrict = options.margin; + + std::string bandCounts; + for (const auto& model : support.models) { + bandCounts += form("%s%.1f G -> %zu bands", bandCounts.empty() ? "" : ", ", model.thresholdKG * 1000., + model.bands.size()); + } + report(" " + bandCounts); + report(form(" margins: strict %.2f cm, tight %.2f cm, edge uncertainty %.2f cm", + support.marginStrict, support.marginTight, support.edgeUncertainty)); + report(""); + + // The model is only worth anything if it really is an outer bound, and only the + // field itself can say so. + report("outer-bound check"); + if (!violationScan(field, support, report)) { + report(""); + report("The support model is not an outer bound on this field, so no placement can be called"); + report("field-free from it. Refusing to classify."); + report.write(options.outputPrefix + "-report.txt"); + return 1; + } + report(""); + + // --- the geometry ----------------------------------------------------------- + TGeoManager::Import(options.geometryFile.c_str()); + if (gGeoManager == nullptr) { + std::cerr << "error: no TGeoManager in " << options.geometryFile << '\n'; + return 1; + } + report(form(" volumes : %d, media %d", gGeoManager->GetListOfVolumes()->GetEntries(), + gGeoManager->GetListOfMedia()->GetEntries())); + + Doctor doctor(field, support); + doctor.walk(gGeoManager->GetTopNode()); + report(form(" placements : %zu classified, %zu detector subtrees pruned", doctor.rows().size(), + doctor.nPruned())); + report(""); + + progress("classify: sampling the field inside every placement that reaches the support"); + doctor.classifyAll(); + doctor.findFindings(); + + std::map verdicts; + for (const auto& row : doctor.rows()) { + ++verdicts[row.verdict]; + } + report("verdicts"); + for (const auto& verdict : verdicts) { + report(form(" %-16s %7d", verdict.first.c_str(), verdict.second)); + } + report(""); + + // --- the reverse audit ------------------------------------------------------ + const double threshold = support.models.front().thresholdKG; + std::map> reverseByVolume; + for (const auto* row : doctor.reverseAudit()) { + auto& entry = reverseByVolume[row->lv + " [" + row->medium + "]"]; + ++entry.first; + entry.second = std::max(entry.second, std::max(row->maxB, row->wholeMaxB)); + } + int inRealField = 0; + for (const auto& entry : reverseByVolume) { + inRealField += entry.second.second > threshold ? 1 : 0; + } + report(form("reverse audit: %zu placements carry a field-free medium yet reach into the field support", + doctor.reverseAudit().size())); + report(form(" %zu logical volumes, %d of them with real field in their own material", + reverseByVolume.size(), inRealField)); + report(form(" %-46s %11s %16s", "volume [medium]", "placements", "max |B| [kG]")); + for (const auto& entry : reverseByVolume) { + report(form(" %-46s %11d %16.4f%s", entry.first.c_str(), entry.second.first, entry.second.second, + entry.second.second > threshold ? " <-- straight-line transport in real field" : "")); + } + report(""); + + // --- the forward findings --------------------------------------------------- + report(form("shared volumes: %zu logical volumes are placed both out of and into the field", + doctor.sharedVolumes().size())); + report(form(" %-28s %8s %8s %s", "logical volume", "out", "in", "status")); + for (size_t i = 0; i < doctor.sharedVolumes().size() && i < 20; ++i) { + const auto& shared = doctor.sharedVolumes()[i]; + bool refused = false; + bool approximate = false; + for (const auto* row : shared.out) { + refused = refused || row->sensitive || doctor.hasSensitive(row->node->GetVolume()); + approximate = approximate || row->approximateExtent; + } + report(form(" %-28s %8zu %8zu %s%s", shared.lv.c_str(), shared.out.size(), shared.in.size(), + refused ? "refused by default (sensitive)" : "proposed", + approximate ? " [extent approximate]" : "")); + } + if (doctor.sharedVolumes().size() > 20) { + report(form(" ... and %zu more, all of them in the proposals file", doctor.sharedVolumes().size() - 20)); + } + report(""); + + report(form("heterogeneous mothers: %zu whose own material straddles the predicate", + doctor.straddlingMothers().size())); + for (size_t i = 0; i < doctor.straddlingMothers().size() && i < 10; ++i) { + const auto* row = doctor.straddlingMothers()[i]; + report(form(" %-40s |B| in own material %.4g .. %.4f kG, %d daughters", row->lv.c_str(), row->wholeMinB, + row->wholeMaxB, row->nDaughters)); + } + report(""); + + report(form("missing containers: %zu daughter clusters lie wholly on the field-free side", + doctor.containers().size())); + for (const auto& container : doctor.containers()) { + report(form(" in %-14s z %9.2f .. %9.2f rmax %7.2f %3d daughters %s%s", container.mother.c_str(), + container.zlo, container.zhi, container.rmax, container.nDaughters, + container.clearedByStrictMargin ? "clear by the strict margin" : "clear by the tight margin", + container.sensitive ? " [refused: sensitive]" : "")); + } + report(""); + + // --- the anchor self-check -------------------------------------------------- + bool anchorsPassed = true; + if (!options.anchorFile.empty()) { + report("anchors"); + anchorsPassed = verifyAnchors(options.anchorFile, doctor, report); + report(anchorsPassed ? " all anchors reproduced" : " ANCHORS FAILED"); + report(""); + } + + // --- outputs ---------------------------------------------------------------- + const std::string proposalsPath = options.outputPrefix + "-proposals.json"; + const std::string tablePath = options.outputPrefix + "-placements.csv"; + const std::string reportPath = options.outputPrefix + "-report.txt"; + std::ofstream proposals(proposalsPath); + proposals << proposalsToJson(doctor, support, options.geometryFile, fieldSource).dump(1, '\t') << '\n'; + writePlacementCsv(doctor.rows(), tablePath); + report("wrote " + proposalsPath + ", " + tablePath + " and " + reportPath); + report.write(reportPath); + + return anchorsPassed ? 0 : 2; +}