From 3c8934a128bbc3215c8a7d693cdf0ac5fcc66499 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Sat, 11 Jul 2026 16:52:46 -0500 Subject: [PATCH 1/6] Update KinKal track extrapolation to include scattering in passive materials (DS and concrete). --- Mu2eKinKal/fcl/prolog.fcl | 3 + Mu2eKinKal/inc/ExtrapolateCRV.hh | 108 -------------- Mu2eKinKal/inc/ExtrapolateCRVRegion.hh | 159 ++++++++++++++++++++ Mu2eKinKal/inc/ExtrapolateCylinders.hh | 123 +++++++++++++++ Mu2eKinKal/inc/ExtrapolatePlanes.hh | 109 ++++++++++++++ Mu2eKinKal/inc/ExtrapolateToR.hh | 75 ++++++++++ Mu2eKinKal/inc/ExtrapolateToZ.hh | 15 +- Mu2eKinKal/inc/KKExtrap.hh | 198 ++++++++++++++++++++++--- Mu2eKinKal/inc/KKFit.hh | 21 ++- Mu2eKinKal/inc/KKFitSettings.hh | 6 + Mu2eKinKal/inc/KKTrack.hh | 44 +++++- Mu2eKinKal/src/KKBField.cc | 13 +- 12 files changed, 736 insertions(+), 138 deletions(-) delete mode 100644 Mu2eKinKal/inc/ExtrapolateCRV.hh create mode 100644 Mu2eKinKal/inc/ExtrapolateCRVRegion.hh create mode 100644 Mu2eKinKal/inc/ExtrapolateCylinders.hh create mode 100644 Mu2eKinKal/inc/ExtrapolatePlanes.hh create mode 100644 Mu2eKinKal/inc/ExtrapolateToR.hh diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 9d4acebcc1..53ab2f1536 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -39,6 +39,7 @@ Mu2eKinKal : { SampleSurfaces : [] # specific to the fit type SaveHitCalibInfo : false SaveDomains : false + SaveMaterialXings : true # save generic passive material Xings (filled during extrapolation) in the KalSeed } CHSEEDFIT: { @@ -401,6 +402,7 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns + MinVLong : 265.0 # mm/ns; near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction ToTrackerEnds : false Upstream : false @@ -414,6 +416,7 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns + MinVLong : 265.0 # mm/ns; near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction ToTrackerEnds : false Upstream : false diff --git a/Mu2eKinKal/inc/ExtrapolateCRV.hh b/Mu2eKinKal/inc/ExtrapolateCRV.hh deleted file mode 100644 index 9ba607f6ae..0000000000 --- a/Mu2eKinKal/inc/ExtrapolateCRV.hh +++ /dev/null @@ -1,108 +0,0 @@ -// predicate to extrapolate to CRV -#ifndef Mu2eKinKal_ExtrapolateCRV_hh -#define Mu2eKinKal_ExtrapolateCRV_hh -#include "Offline/Mu2eKinKal/inc/KKTrack.hh" -#include "KinKal/General/TimeDir.hh" -#include "KinKal/General/TimeRange.hh" -#include "KinKal/Geometry/Intersection.hh" -#include "KinKal/Geometry/Rectangle.hh" -#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" -#include "Offline/KinKalGeom/inc/CRV.hh" -#include "cetlib_except/exception.h" -#include -#include -#include -namespace mu2e { - using KinKal::TimeDir; - using KinKal::TimeRange; - using KKGeom::CRV; - using KinKal::Rectangle; - using KinKal::Intersection; - class ExtrapolateCRV { - public: - using CRVSV = std::vector; - struct CRVI { - Intersection inter_; - double whw_ = 0; - int isect_ = -1; - CRVI(){} - CRVI(Intersection const& inter,double whw, int isect) : inter_(inter),whw_(whw),isect_(isect) {} - }; - using CRVIV = std::vector; - - struct sortIntersections { // sort in time direction - TimeDir tdir_; - bool operator () (CRVI const& inter1, CRVI const& inter2) { - return tdir_ == TimeDir::forwards ? inter1.inter_.time_ < inter2.inter_.time_ : inter1.inter_.time_ > inter2.inter_.time_; - } - sortIntersections(TimeDir tdir) : tdir_(tdir) {} - }; - - ExtrapolateCRV(double maxdt, double maxdtstep, double dptol, double intertol, double minv, CRV const& crv, int debug=0) : - maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), debug_(debug), sectors_(crv.sectors()) {} - // interface for extrapolation - double maxDt() const { return maxDt_; } - double maxDtStep() const { return maxDtStep_; } - double dpTolerance() const { return dptol_; } - double interTolerance() const { return intertol_; } - // CRV specific - auto const& sectors() const { return sectors_; } - auto const& sector(size_t isect) const { return sectors_[isect]; } - auto const& intersections() const { return inters_; } - int debug() const { return debug_; } - void reset() const { inters_.clear(); } - // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time - template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; - private: - double maxDt_ = -1; // maximum extrapolation time - double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration - double dptol_ = 1e10; // fractional momentum tolerance - double intertol_ = 1e10; // intersection tolerance (mm) - double minvnorm_ = 1e-5; // minimum vel normal (outwards) to plane - int debug_ = 0; // debug level - CRVSV sectors_; // sectors cache - mutable CRVIV inters_; // cache of most recent intersections - }; - - template bool ExtrapolateCRV::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { - // we are answering the question: did the segment last added to this extrapolated track cross a CRV sector or not? - // if so, stop extrapolating (for now). If not, and if we're still heading towards at least 1 CRV sector , keep going. - // cache the intersection if it's found, so it can be used without recomputing it. - bool retval(false); - reset(); - auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); - // Make sure the range is positive definite - static const double epsilon(1.0e-6); - auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; - auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); - TimeRange trange(stime,etime,false); - // test from the start of the end piece - auto pos = ktraj.position3(stime); - auto vel = ktraj.velocity(stime); - if(debug_ > 4)std::cout << "CRV extrap tdir " << tdir << " stime " << stime << " etime " << etime << " vel " << vel << " pos " << pos << std::endl; - for(size_t isect = 0; isect < sectors_.size(); ++isect ){ - auto const& sector = sectors_[isect]; - double normvel = vel.Dot(sector.sector_->normal())*timeDirSign(tdir); // sign by extrapolation direction - double sdist = (sector.sector_->center()-pos).Dot(vel)*timeDirSign(tdir); - if(debug_ > 4)std::cout << "CRV extrap normvel " << normvel << " sdist " << sdist << std::endl; - // stop if horizontal or plane is past the current point. - if(fabs(normvel) < minvnorm_ || sdist < 0 )continue; - // try to intersect - auto newinter = KinKal::intersect(fittraj,*sector.sector_,trange,intertol_,tdir); - if(debug_ > 3)std::cout << "CRV " << newinter << std::endl; - if(newinter.good()){ - inters_.emplace_back(newinter,sector.whw_,(int)isect); - if(debug_ > 1)std::cout << "Good CRV " << newinter << " sector " << sector.sname_ << std::endl; - } else if ( newinter.onsurface_ && newinter.inbounds_) { // inbounds might be too strict for CentralHelix tracks, will need to check TODO - retval |= trange.beyond(newinter.time_,tdir); - // there's a potential intersection, but the trajectory hasn't gotten there yet. Tell the track to keep extending - if(trange.beyond(newinter.time_,tdir) && debug_ > 2)std::cout << "Potential CRV " << newinter << std::endl; - } - } - // sort intersections in the time direction - sortIntersections isort(tdir); - std::sort(inters_.begin(),inters_.end(),isort); - return retval; - } -} -#endif diff --git a/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh b/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh new file mode 100644 index 0000000000..5aee23aed7 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh @@ -0,0 +1,159 @@ +// Predicate to extrapolate through the CRV region, crossing BOTH the CRV sector planes AND the +// per-sector aluminium strongback (passive) planes in ONE interleaved, time-ordered pass. +// +// WHY THIS EXISTS +// Each CRV module has a thin (~12.7 mm) Al strongback on the tracker-facing side of its scintillator +// stack. The sector (veto) plane is modeled at the stack center, so the strongback plane sits ~60 mm +// (about a stack half-thickness) toward the tracker from it; going outward the surfaces alternate +// strongback -> sector for every module. KinKal xings can only be ADDED at the trajectory front/back (an interior xing cannot +// be inserted), so these surfaces must be crossed in strict spatial/time order. The previous code +// crossed them in two SEPARATE passes -- all strongbacks (ExtrapolatePlanes), then all sectors +// (ExtrapolateCRV). In a STACKED-sector geometry (the extracted KPP modules EX/T1/T2 sit ~150 mm +// apart in y) the strongback pass advanced the trajectory front PAST the lower sectors, leaving +// them stranded in the interior where they could no longer be added -- CRV_T1 was lost on every +// track, and EX/T2 were under-counted. Searching both surface types together, this predicate stops +// the extrapolation at the NEXT crossing of EITHER type, so the KKExtrap loop adds them strictly +// outward (T1_SB -> T1 -> EX_SB -> EX -> ...) and nothing is stranded. +// +// RUN-2 SAFETY: run-2 geometry presents at most one CRV sector (and its one strongback) per +// extrapolation direction, so the outward order is unchanged (strongback then sector) and no +// stacking occurs. Concrete passive planes are NOT handled here (they stay in +// extrapolatePassiveMaterial); only the CRV strongbacks are merged with the sectors. Per-surface +// grazing floors are preserved exactly: sectors keep the z-normal end-cap floor (minvlong), the +// strongback planes keep minvnorm, identical to the originals. +#ifndef Mu2eKinKal_ExtrapolateCRVRegion_hh +#define Mu2eKinKal_ExtrapolateCRVRegion_hh +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/Rectangle.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/CRV.hh" +#include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include +#include +#include +#include +namespace mu2e { + using KinKal::TimeDir; + using KinKal::TimeRange; + using KKGeom::CRV; + using KinKal::Intersection; + using KinKal::timeDirSign; + class ExtrapolateCRVRegion { + public: + using CRVSV = std::vector; + using MaterialPlane = KinKalGeom::PassiveMaterialPlane; + using MaterialPlaneCollection = std::vector; + // a crossing of either a CRV sector or a strongback passive plane + struct RegionXing { + Intersection inter_; + bool isSector_ = false; + int sectorIdx_ = -1; // index into sectors_ (sector crossings) + double whw_ = 0; // sector module half-width (sector crossings) + MaterialPlane const* plane_ = nullptr; // the strongback plane (plane crossings) + RegionXing() {} + RegionXing(Intersection const& inter, int isect, double whw) : + inter_(inter), isSector_(true), sectorIdx_(isect), whw_(whw) {} + RegionXing(Intersection const& inter, MaterialPlane const& plane) : + inter_(inter), isSector_(false), plane_(&plane) {} + }; + using RegionXingCollection = std::vector; + struct sortXings { // sort in the extrapolation time direction + TimeDir tdir_; + bool operator () (RegionXing const& a, RegionXing const& b) { + return tdir_ == TimeDir::forwards ? a.inter_.time_ < b.inter_.time_ : a.inter_.time_ > b.inter_.time_; + } + sortXings(TimeDir tdir) : tdir_(tdir) {} + }; + + ExtrapolateCRVRegion(double maxdt, double maxdtstep, double dptol, double intertol, + double minv, double minvlong, CRV const& crv, MaterialPlaneCollection const& planes, int debug=0) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), + minvnorm_(minv), minvlong_(minvlong), debug_(debug), sectors_(crv.sectors()), planes_(planes) { + // per-sector grazing floor: stricter (minvlong) for z-normal end-cap planes, looser (minvnorm) + // for the oblique-cosmic top/side planes -- identical to ExtrapolateCRV. + sectorminv_.reserve(sectors_.size()); + for(auto const& sector : sectors_) + sectorminv_.push_back( fabs(sector.sector_->normal().Z()) > 0.9 ? minvlong_ : minvnorm_ ); + } + double maxDt() const { return maxDt_; } + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + auto const& sectors() const { return sectors_; } + auto const& sector(size_t isect) const { return sectors_[isect]; } + auto const& intersections() const { return inters_; } + int debug() const { return debug_; } + void reset() const { inters_.clear(); } + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + private: + double maxDt_ = -1; + double maxDtStep_ = -1; + double dptol_ = 1e10; + double intertol_ = 1e10; + double minvnorm_ = 1e-5; + double minvlong_ = 1e-5; + int debug_ = 0; + CRVSV sectors_; + MaterialPlaneCollection const& planes_; + std::vector sectorminv_; + mutable RegionXingCollection inters_; + }; + + template bool ExtrapolateCRVRegion::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { + // Stop the extrapolation at the NEXT crossing of any CRV sector OR strongback plane. The + // KKExtrap loop adds whichever crossing(s) this piece produced (time-sorted, outermost first) + // then re-extrapolates, so the surfaces are crossed strictly outward and none is left interior. + bool retval(false); + reset(); + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + static const double epsilon(1.0e-6); + static const double coincidenttol(1.0e-3); // ns; crossings closer than this are the same physical point (a seam) + auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; + auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + TimeRange trange(stime,etime,false); + auto pos = ktraj.position3(stime); + auto vel = ktraj.velocity(stime); + double piece_span = fabs(stime - etime); + // helper: has a crossing at this time already been cached (a coplanar seam)? + auto coincident = [&](double t){ for(auto const& p : inters_) if(fabs(p.inter_.time_ - t) < coincidenttol) return true; return false; }; + // --- CRV sectors --- + for(size_t isect = 0; isect < sectors_.size(); ++isect ){ + auto const& sector = sectors_[isect]; + double normvel = vel.Dot(sector.sector_->normal())*timeDirSign(tdir); + double signed_perp = (sector.sector_->center()-pos).Dot(sector.sector_->normal()); + double time_to_sector = signed_perp / normvel; + // skip grazing, or planes already behind this piece's evaluation point (already crossed/added) + if(fabs(normvel) < sectorminv_[isect] || time_to_sector < 0 )continue; + auto newinter = KinKal::intersect(fittraj,*sector.sector_,trange,intertol_,tdir); + if(newinter.good()){ + if(!coincident(newinter.time_)) inters_.emplace_back(newinter,(int)isect,sector.whw_); + } else if(newinter.onsurface_ && newinter.inbounds_) { + retval |= trange.beyond(newinter.time_,tdir); // crossing just beyond this piece: keep extrapolating + } else if(!newinter.onsurface_ && tdir == TimeDir::backwards) { + if(time_to_sector > piece_span) retval = true; // far sector beyond this (short) piece: keep going + } + } + // --- strongback passive planes (minvnorm floor, like ExtrapolatePlanes) --- + for(auto const& plane : planes_){ + double normvel = vel.Dot(plane.surface_->normal())*timeDirSign(tdir); + double signed_perp = (plane.surface_->center()-pos).Dot(plane.surface_->normal()); + double time_to_plane = signed_perp / normvel; + if(fabs(normvel) < minvnorm_ || time_to_plane < 0 )continue; + auto newinter = KinKal::intersect(fittraj,*plane.surface_,trange,intertol_,tdir); + if(newinter.good()){ + if(!coincident(newinter.time_)) inters_.emplace_back(newinter,plane); + } else if(newinter.onsurface_ && newinter.inbounds_) { + retval |= trange.beyond(newinter.time_,tdir); + } else if(!newinter.onsurface_ && tdir == TimeDir::backwards) { + if(time_to_plane > piece_span) retval = true; + } + } + // sort the crossings in the time direction so the KKExtrap loop adds them outermost-first + std::sort(inters_.begin(),inters_.end(),sortXings(tdir)); + if(!inters_.empty()) return false; // stop: this piece crossed the CRV region; record now + return retval; + } +} +#endif diff --git a/Mu2eKinKal/inc/ExtrapolateCylinders.hh b/Mu2eKinKal/inc/ExtrapolateCylinders.hh new file mode 100644 index 0000000000..f24483d3bc --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateCylinders.hh @@ -0,0 +1,123 @@ +// Predicate to extrapolate to the next intersection with passive cylindrical material shells. +#ifndef Mu2eKinKal_ExtrapolateCylinders_hh +#define Mu2eKinKal_ExtrapolateCylinders_hh + +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/DetectorSolenoid.hh" + +#include +#include +#include +#include +#include + +namespace mu2e { + using KinKal::Intersection; + using KinKal::TimeDir; + using KinKal::TimeRange; + using KinKal::timeDirSign; + + class ExtrapolateCylinders { + public: + using MaterialCylinder = KKGeom::DetectorSolenoid::MaterialCylinder; + using MaterialCylinderCollection = KKGeom::DetectorSolenoid::MaterialCylinderCollection; + struct CylinderIntersection { + Intersection inter_; + MaterialCylinder const* cylinder_ = nullptr; + CylinderIntersection() {} + CylinderIntersection(Intersection const& inter, MaterialCylinder const& cylinder) : + inter_(inter), cylinder_(&cylinder) {} + }; + using CylinderIntersectionCollection = std::vector; + + struct SortIntersections { + TimeDir tdir_; + bool operator () (CylinderIntersection const& lhs, CylinderIntersection const& rhs) { + return tdir_ == TimeDir::forwards ? lhs.inter_.time_ < rhs.inter_.time_ : lhs.inter_.time_ > rhs.inter_.time_; + } + SortIntersections(TimeDir tdir) : tdir_(tdir) {} + }; + + ExtrapolateCylinders(double maxdt, double maxdtstep, double dptol, double intertol, + double minv, MaterialCylinderCollection const& cylinders, int debug=0) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), cylinders_(cylinders), debug_(debug) { + for(auto const& cylinder : cylinders_) { + maxRadius_ = std::max(maxRadius_, cylinder.surface_->radius() + 0.5*cylinder.thickness_); + } + } + + double maxDt() const { return maxDt_; } + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + auto const& intersections() const { return inters_; } + int debug() const { return debug_; } + void reset() const { inters_.clear(); } + + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + + private: + double maxDt_ = -1; + double maxDtStep_ = -1; + double dptol_ = 1e10; + double intertol_ = 1e10; + double minvnorm_ = 1e-5; + double minShellCosine_ = 0.1; + MaterialCylinderCollection const& cylinders_; + int debug_ = 0; + double maxRadius_ = 0.0; + mutable CylinderIntersectionCollection inters_; + }; + + template bool ExtrapolateCylinders::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { + reset(); + if(cylinders_.empty()) return false; + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + static const double epsilon(1.0e-6); + if(ktraj.range().range() <= epsilon) return true; + double stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; + double etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + // March the shell search across the end piece in BOUNDED sub-windows (capped by maxDtStep) from the + // near end outward, rather than one intersect over the whole piece. The thin concentric shells are + // closely spaced (~960-1315mm), so a single intersect over a piece that spans several of them can + // step over the inner ones (they were under-found, innermost worst). Marching brackets the nearest + // crossings and records them HERE, with no overstep; the KKExtrap loop banks/de-dups and re-extends. + double const sgn = (etime >= stime) ? 1.0 : -1.0; + double const dtwin = (maxDtStep_ > 0.0) ? maxDtStep_ : std::fabs(etime - stime); + double t0 = stime; + while(sgn*(etime - t0) > epsilon) { + // stop once we are outside the shells AND moving further out: any shell is behind us + auto p0 = ktraj.position3(t0); + double const rho0 = std::sqrt(p0.X()*p0.X() + p0.Y()*p0.Y()); + auto v0 = ktraj.velocity(t0); + double const vrho0 = rho0 > 0.0 ? (p0.X()*v0.X() + p0.Y()*v0.Y())/rho0 : 0.0; + if(rho0 > maxRadius_ && vrho0*timeDirSign(tdir) > 0.0) return false; + double const t1 = t0 + sgn*std::min(dtwin, std::fabs(etime - t0)); + TimeRange win = (sgn > 0.0) ? TimeRange(t0,t1,false) : TimeRange(t1,t0,false); + for(auto const& cylinder : cylinders_) { + auto newinter = KinKal::intersect(fittraj,*cylinder.surface_,win,intertol_,tdir); + if(debug_ > 3)std::cout << cylinder.sid_ << " " << newinter << std::endl; + if(newinter.good()) { + double const normvel = std::fabs(newinter.norm_.Dot(newinter.pdir_)); + // thin-shell approximation: skip near-tangent crossings (unbounded path through a finite shell) + if(normvel > std::max(minvnorm_,minShellCosine_)) inters_.emplace_back(newinter,cylinder); + else if(debug_ > 1) std::cout << "Skipping grazing DS material intersection " << cylinder.sid_ + << " normvel " << normvel << std::endl; + } + } + if(!inters_.empty()) { + SortIntersections isort(tdir); + std::sort(inters_.begin(),inters_.end(),isort); + return false; // nearest shell crossing(s) found and recorded here -- no overstep + } + t0 = t1; + } + // reached the piece end with no crossing while still approaching the shells: extend and retry + return true; + } +} + +#endif diff --git a/Mu2eKinKal/inc/ExtrapolatePlanes.hh b/Mu2eKinKal/inc/ExtrapolatePlanes.hh new file mode 100644 index 0000000000..05510595b8 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolatePlanes.hh @@ -0,0 +1,109 @@ +// Predicate to extrapolate to the next intersection with passive planar material shells. +#ifndef Mu2eKinKal_ExtrapolatePlanes_hh +#define Mu2eKinKal_ExtrapolatePlanes_hh + +#include "KinKal/General/TimeDir.hh" +#include "KinKal/General/TimeRange.hh" +#include "KinKal/Geometry/Intersection.hh" +#include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "Offline/KinKalGeom/inc/KinKalGeom.hh" + +#include +#include +#include +#include + +namespace mu2e { + using KinKal::Intersection; + using KinKal::TimeDir; + using KinKal::TimeRange; + using KinKal::timeDirSign; + + class ExtrapolatePlanes { + public: + using MaterialPlane = KinKalGeom::PassiveMaterialPlane; + using MaterialPlaneCollection = KinKalGeom::PassiveMaterialPlaneCollection; + struct PlaneIntersection { + Intersection inter_; + MaterialPlane const* plane_ = nullptr; + PlaneIntersection() {} + PlaneIntersection(Intersection const& inter, MaterialPlane const& plane) : + inter_(inter), plane_(&plane) {} + }; + using PlaneIntersectionCollection = std::vector; + + struct SortIntersections { + TimeDir tdir_; + bool operator () (PlaneIntersection const& lhs, PlaneIntersection const& rhs) { + return tdir_ == TimeDir::forwards ? lhs.inter_.time_ < rhs.inter_.time_ : lhs.inter_.time_ > rhs.inter_.time_; + } + SortIntersections(TimeDir tdir) : tdir_(tdir) {} + }; + + ExtrapolatePlanes(double maxdt, double maxdtstep, double dptol, double intertol, + double minv, MaterialPlaneCollection const& planes, int debug=0) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), planes_(planes), debug_(debug) {} + + double maxDt() const { return maxDt_; } + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + auto const& intersections() const { return inters_; } + int debug() const { return debug_; } + void reset() const { inters_.clear(); } + + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + + private: + double maxDt_ = -1; + double maxDtStep_ = -1; + double dptol_ = 1e10; + double intertol_ = 1e10; + double minvnorm_ = 1e-5; + MaterialPlaneCollection const& planes_; + int debug_ = 0; + mutable PlaneIntersectionCollection inters_; + }; + + template bool ExtrapolatePlanes::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { + reset(); + if(planes_.empty()) return false; + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + static const double epsilon(1.0e-6); + if(ktraj.range().range() <= epsilon) return true; + auto stime = tdir == TimeDir::forwards ? ktraj.range().begin()+epsilon : ktraj.range().end()-epsilon; + auto etime = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + TimeRange trange(stime,etime,false); + auto pos = ktraj.position3(stime); + auto vel = ktraj.velocity(stime); + double piece_span = std::fabs(stime - etime); + bool keepGoing = false; + for(auto const& plane : planes_) { + double normvel = vel.Dot(plane.surface_->normal())*timeDirSign(tdir); + double pdist = (plane.surface_->center()-pos).Dot(plane.surface_->normal()); + if(debug_ > 4)std::cout << "Passive plane extrap " << plane.sid_ << " normvel " << normvel << " pdist " << pdist << std::endl; + if(std::fabs(normvel) < minvnorm_) continue; + if(pdist*normvel < 0.0) continue; + auto newinter = KinKal::intersect(fittraj,*plane.surface_,trange,intertol_,tdir); + if(debug_ > 3)std::cout << plane.sid_ << " " << newinter << std::endl; + if(newinter.good()) { + inters_.emplace_back(newinter,plane); + } else if(newinter.onsurface_ && newinter.inbounds_) { + keepGoing |= trange.beyond(newinter.time_,tdir); + } else if(!newinter.onsurface_ && tdir == TimeDir::backwards) { + // a far plane beyond this (short) piece: the range-bounded stepIntersect can't reach it, so + // use the analytic ray estimate (pdist/normvel) to keep marching outward (mirrors + // ExtrapolateCRVRegion). This lets a tight KinematicLine-tail seed reach the concrete plane + // instead of overshooting it (and the CRV) to maxDt. + if(pdist/normvel > piece_span) keepGoing = true; + } + } + SortIntersections isort(tdir); + std::sort(inters_.begin(),inters_.end(),isort); + // stop as soon as a crossing is cached so the KKExtrap loop records it HERE (no overshoot) + if(!inters_.empty()) return false; + return keepGoing; + } +} + +#endif diff --git a/Mu2eKinKal/inc/ExtrapolateToR.hh b/Mu2eKinKal/inc/ExtrapolateToR.hh new file mode 100644 index 0000000000..74f51ea7dd --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateToR.hh @@ -0,0 +1,75 @@ +// predicate for extrapolation to a given transverse radius rho = sqrt(x^2+y^2), bounded in z. +// This is the radial analog of ExtrapolateToZ: it targets a constant-rho (cylindrical) surface rather +// than a constant-z plane. That is the natural tracker exit for a CentralHelix (cosmic) track, which +// traverses the tracker radially (mostly in y) and leaves through the cylindrical side, NOT through the +// z-end planes. The [minz,maxz] bounds keep a near-horizontal (large v_z) cosmic from chasing the target +// radius far down the beamline, past the tracker ends -- the symmetric guard to ExtrapolateToZ's behaviour +// for near-vertical cosmics. +#ifndef Mu2eKinKal_ExtrapolateToR_hh +#define Mu2eKinKal_ExtrapolateToR_hh +#include "KinKal/Trajectory/ParticleTrajectory.hh" +#include "KinKal/General/TimeDir.hh" +#include +#include +#include +#include "cetlib_except/exception.h" +namespace mu2e { + using KinKal::TimeDir; + using KinKal::timeDirSign; + class ExtrapolateToR { + public: + ExtrapolateToR(double maxdt, double maxdtstep, double tol, double rval, double minz, double maxz, int debug) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), rval_(rval), minz_(minz), maxz_(maxz), debug_(debug) {} + // interface for extrapolation + double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } // tolerance on fractional momentum change + double rVal() const { return rval_; } + double minZ() const { return minz_; } + double maxZ() const { return maxz_; } + int debug() const { return debug_; } + // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + private: + double maxDt_ = -1; // maximum extrapolation time + double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration + double dptol_ = 1e10; // fractional momentum tolerance in BField domain + double rval_ = 0; // target transverse radius (mm) + double minz_ = -std::numeric_limits::max(); // minimum z allowed during extrapolation (mm) + double maxz_ = std::numeric_limits::max(); // maximum z allowed during extrapolation (mm) + int debug_ = 0; // debug level + }; + + template bool ExtrapolateToR::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + auto time = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + auto vel = ktraj.velocity(time); + auto pos = ktraj.position3(time); + double rho = std::sqrt(pos.X()*pos.X() + pos.Y()*pos.Y()); + double z = pos.Z(); + // radial velocity, signed by the extrapolation direction (positive = heading outward) + double rvel = (rho > 0.0 ? (pos.X()*vel.X() + pos.Y()*vel.Y())/rho : 0.0)*timeDirSign(tdir); + double vz = vel.Z()*timeDirSign(tdir); // z velocity, signed by the extrapolation direction + if(debug_ > 2)std::cout << "R extrap start time " << time << " rho " << rho << " rvel " << rvel << " z " << z + << " vz " << vz << " rval " << rval_ << " z[" << minz_ << "," << maxz_ << "]" << std::endl; + // already outside the z bounds, or already at/beyond the target radius: stop + if(z > maxz_ || z < minz_) return false; + if(rho >= rval_) return false; + // heading inward: this direction will never reach the outer target radius + if(rvel <= 0.0) return false; + // Decide GEOMETRICALLY which boundary the track reaches first (the BField domain step out here is set + // by the field tolerance, not MaxDtStep, and can be large -- so a predicate that only stops AFTER the + // step overshoots the z bound by a whole domain). A near-axial cosmic (large |v_z|, small v_rho) + // reaches the z end before the target radius; extending it toward the radius is exactly what overshoots + // the z end, carrying rho far past the DS shells and stranding them as interior pieces for the material + // pass. In that case stop at the tracker boundary now and let the DS-material pass (no z bound, marches + // shell-by-shell) extend it. Otherwise the track exits radially: keep going to the target radius. + double time_to_rval = (rval_ - rho)/rvel; + double time_to_zbound = std::numeric_limits::max(); + if(vz > 0.0) time_to_zbound = (maxz_ - z)/vz; + else if(vz < 0.0) time_to_zbound = (z - minz_)/(-vz); + if(time_to_zbound < time_to_rval) return false; // exits via the z end: stop at the tracker boundary + return true; // radial exit: keep extending toward the target radius + } +} +#endif diff --git a/Mu2eKinKal/inc/ExtrapolateToZ.hh b/Mu2eKinKal/inc/ExtrapolateToZ.hh index f4a4fc513f..ffcf9d73a0 100644 --- a/Mu2eKinKal/inc/ExtrapolateToZ.hh +++ b/Mu2eKinKal/inc/ExtrapolateToZ.hh @@ -5,13 +5,16 @@ #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/General/TimeDir.hh" #include +#include #include "cetlib_except/exception.h" namespace mu2e { using KinKal::TimeDir; using KinKal::timeDirSign; class ExtrapolateToZ { public: - ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double zval,int debug) : maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), zval_(zval), debug_(debug) {} + ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double zval,int debug, + double maxradius = std::numeric_limits::max()) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), zval_(zval), debug_(debug), maxradius_(maxradius) {} // interface for extrapolation double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement double maxDtStep() const { return maxDtStep_; } @@ -26,6 +29,9 @@ namespace mu2e { double dptol_ = 1e10; // fractional momentum tolerance in BField domain double zval_ = 0; // z value targeted int debug_ = 0; // debug level + double maxradius_ = std::numeric_limits::max(); // max transverse radius rho=sqrt(x^2+y^2) (mm) + // the track may reach; stops near-vertical (small-v_z) cosmics from chasing the + // constant-z target plane out to large rho. Default = no limit. }; template bool ExtrapolateToZ::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { @@ -37,7 +43,12 @@ namespace mu2e { double zval = pos.Z(); auto const& bnom = ktraj.bnom(time); double zref = vel.R()*fabs(sin(bnom.Theta())); - if(debug_ > 2)std::cout << "Z extrap start time " << time << " z " << zval << " zvel " << zvel << " zref " << zref << std::endl; + double rho = std::sqrt(pos.X()*pos.X() + pos.Y()*pos.Y()); + if(debug_ > 2)std::cout << "Z extrap start time " << time << " z " << zval << " zvel " << zvel << " zref " << zref << " rho " << rho << std::endl; + // stop if we have reached the maximum allowed transverse radius: a near-vertical (small-v_z) track + // would otherwise keep chasing the constant-z target plane out to large rho (e.g. a cosmic running + // metres in y), pre-extending the trajectory past the DS before the DS-material pass. Default = no limit. + if(rho > maxradius_) return false; // if z velocity is unreliable, continue if(fabs(zvel) < zref) return true; // stop if we're heading away from the target z diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index c4f32c80d4..829d2e20d2 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -10,13 +10,18 @@ #include "Offline/KinKalGeom/inc/Tracker.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateToR.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateCRV.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateCylinders.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolatePlanes.hh" #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/GeometryService/inc/GeomHandle.hh" +#include +#include #include "Offline/GeometryService/inc/GeometryService.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" namespace mu2e { @@ -38,15 +43,17 @@ namespace mu2e { template bool extrapolateST(KKTrack& ktrk,TimeDir trkdir) const; template bool extrapolateTracker(KKTrack& ktrk,TimeDir tdir) const; template bool extrapolateTSDA(KKTrack& ktrk,TimeDir tdir) const; + template void extrapolateDSMaterial(KKTrack& ktrk,TimeDir tdir) const; + template void extrapolatePassiveMaterial(KKTrack& ktrk,TimeDir tdir) const; template void extrapolateCRV(KKTrack& ktrk,TimeDir tdir) const; template void extrapolateOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; template void toTrackerEnds(KKTrack& ktrk) const; - // TODO add DS and shielding material private: int debug_; - double btol_, intertol_, maxdt_, maxdtstep_, minv_; - bool backToTracker_, extrapolateOPA_, toTrackerEnds_, upstream_, toCRV_; + double btol_, intertol_, maxdt_, maxdtstep_, minv_, minvlong_, maxtrackerendsradius_; + double mintrackerendsz_, maxtrackerendsz_; + bool backToTracker_, extrapolateOPA_, toTrackerEnds_, toTrackerEndsRadius_, upstream_, toCRV_; double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO }; @@ -58,9 +65,14 @@ namespace mu2e { maxdt_(extrapconfig.MaxDt()), maxdtstep_(extrapconfig.MaxDtStep()), minv_(extrapconfig.MinV()), + minvlong_(extrapconfig.MinVLong()), + maxtrackerendsradius_(extrapconfig.MaxTrackerEndsRadius()), + mintrackerendsz_(extrapconfig.MinTrackerEndsZ()), + maxtrackerendsz_(extrapconfig.MaxTrackerEndsZ()), backToTracker_(extrapconfig.BackToTracker()), extrapolateOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), + toTrackerEndsRadius_(extrapconfig.ToTrackerEndsRadius()), upstream_(extrapconfig.Upstream()), toCRV_(extrapconfig.ToCRV()) {} @@ -103,7 +115,11 @@ namespace mu2e { } // try CRV in both directions if(toCRV_){ + extrapolateDSMaterial(ktrk,TimeDir::backwards); + extrapolatePassiveMaterial(ktrk,TimeDir::backwards); extrapolateCRV(ktrk,TimeDir::backwards); + extrapolateDSMaterial(ktrk,TimeDir::forwards); + extrapolatePassiveMaterial(ktrk,TimeDir::forwards); extrapolateCRV(ktrk,TimeDir::forwards); } } @@ -111,8 +127,38 @@ namespace mu2e { template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { GeomHandle kkg_h; - ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_); - ExtrapolateToZ trackerBack(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->back().center().Z(),debug_); + // CentralHelix (cosmic) mode: extrapolate to the tracker OUTER CYLINDER (constant rho), not the z-end + // planes. A cosmic traverses the tracker radially (mostly in y), so the outer cylinder is its true exit; + // ExtrapolateToZ would chase the constant-z plane out to large rho. The z bounds keep near-horizontal + // cosmics from chasing the radius past the tracker ends. + if(toTrackerEndsRadius_){ + auto const& outer = kkg_h->tracker()->outer(); + ExtrapolateToR toOuter(maxdt_,maxdtstep_,btol_,outer.radius(),mintrackerendsz_,maxtrackerendsz_,debug_); + auto const& ftraj = ktrk.fitTraj(); + // the cosmic crosses the outer cylinder on both ends; extend in both time directions + auto tofwd = ktrk.extrapolate(TimeDir::forwards,toOuter); + auto tobwd = ktrk.extrapolate(TimeDir::backwards,toOuter); + static const SurfaceId tt_outer("TT_Outer"); + static const SurfaceId tt_mid("TT_Mid"); + // middle reference, recorded only if the track actually crosses it + auto midinter = KinKal::intersect(ftraj,kkg_h->tracker()->middle(),ftraj.range(),intertol_); + if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); + // outer-cylinder crossings (tracker entry and exit) from the respective end pieces + if(tofwd){ + auto const& fhel = ftraj.back(); + auto inter = KinKal::intersect(fhel,outer,fhel.range(),intertol_,TimeDir::forwards); + if(inter.good()) ktrk.addIntersection(tt_outer,inter); + } + if(tobwd){ + auto const& bhel = ftraj.front(); + auto inter = KinKal::intersect(bhel,outer,bhel.range(),intertol_,TimeDir::backwards); + if(inter.good()) ktrk.addIntersection(tt_outer,inter); + } + return; + } + + ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_,maxtrackerendsradius_); + ExtrapolateToZ trackerBack(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->back().center().Z(),debug_,maxtrackerendsradius_); // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't // reflect inside the tracker volume auto const& ftraj = ktrk.fitTraj(); @@ -284,32 +330,142 @@ namespace mu2e { } //extrapolate to the CRV sectors. This only makes sense for KKLine or CentralHelix + template void KKExtrap::extrapolateDSMaterial(KKTrack& ktrk,TimeDir tdir) const { + using KKMATCYLXING = KKShellXing; + GeomHandle kkg_h; + GeomHandle kkmat_h; + auto const& cylinders = kkg_h->DS()->materialCylinders(); + if(cylinders.empty()) return; + auto extrapCylinders = ExtrapolateCylinders(maxdt_,maxdtstep_,btol_,intertol_,minv_,cylinders,debug_); + auto const& ftraj = ktrk.fitTraj(); + using MaterialCylinder = ExtrapolateCylinders::MaterialCylinder; + std::unordered_set banked; // bank each cylinder once (avoid re-found duplicates) + do { + ktrk.extrapolate(tdir,extrapCylinders); + if(debug_ > 0) std::cout << "Found " << extrapCylinders.intersections().size() << " DS material intersections " << std::endl; + // Bank every crossing found in this pass, in strict tracker-inward time order (descending time + // for backward, ascending for forward). The DS material shells are concentric and closely spaced, + // so a single extrapolation step overshoots the whole cluster and the trajectory piece spans all + // of them; adding the crossings in this order makes each prepend/append shrink the piece down to + // exactly the next crossing, so every one stays within the addable front/back range and none is + // stranded as interior. (intersections() is not reliably time-sorted, so sort a local copy.) The + // 'banked' set skips a crossing re-found on a later step; once a pass adds nothing new, we are done. + auto inters = extrapCylinders.intersections(); + std::sort(inters.begin(),inters.end(),[tdir](auto const& a,auto const& b){ + return tdir == TimeDir::forwards ? a.inter_.time_ < b.inter_.time_ : a.inter_.time_ > b.inter_.time_; }); + bool newbank = false; + for(auto const& inter : inters) { + if(!banked.insert(inter.cylinder_).second) continue; + newbank = true; + auto const& cylinder = *inter.cylinder_; + auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + auto matxingptr = std::make_shared(cylinder.surface_,cylinder.sid_,*kkmat_h->material(cylinder.material_), + inter.inter_,reftrajptr,cylinder.thickness_,extrapCylinders.interTolerance(),kkmat_h->applyBetheCorrection()); + ktrk.addMaterialCylXing(matxingptr,tdir); + } + if(!newbank) break; + } while(extrapCylinders.intersections().size()>0); + } + + template void KKExtrap::extrapolatePassiveMaterial(KKTrack& ktrk,TimeDir tdir) const { + using KKMATRECXING = KKShellXing; // passive plane (Rectangle or Annulus) + using MaterialPlane = ExtrapolatePlanes::MaterialPlane; + GeomHandle kkg_h; + GeomHandle kkmat_h; + auto const& planes = kkg_h->passiveMaterialPlanes(); + if(planes.empty()) return; + auto extrapPlanes = ExtrapolatePlanes(maxdt_,maxdtstep_,btol_,intertol_,minv_,planes,debug_); + auto const& ftraj = ktrk.fitTraj(); + + // Normal path: process each plane at most once. + // Stop before calling ktrk.extrapolate when all planes have been processed so + // we don't advance the trajectory past the last plane (which would overshoot the + // CRV and prevent extrapolateCRV from finding it). + std::unordered_set processed; + // The CRV strongback planes are crossed together with the CRV sectors in extrapolateCRV + // (interleaved, time-ordered): processing all strongbacks HERE advances the front past the + // lower sectors of a stacked geometry and strands them. Exclude the strongbacks from this + // (concrete-only for the extracted geometry) pass by pre-marking them processed. + SurfaceId const sbid(SurfaceIdEnum::CRV_StrongBack); + for(auto const& plane : planes) if(plane.sid_.id() == sbid.id()) processed.insert(&plane); + do { + if(processed.size() >= planes.size()) break; + ktrk.extrapolate(tdir,extrapPlanes); + if(debug_ > 0) std::cout << "Found " << extrapPlanes.intersections().size() << " passive material plane intersections " << std::endl; + bool found = false; + for(auto const& inter : extrapPlanes.intersections()) { + if(processed.count(inter.plane_) > 0) continue; + auto const& plane = *inter.plane_; + auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); + auto matxingptr = std::make_shared(plane.surface_,plane.sid_,*kkmat_h->material(plane.material_), + inter.inter_,reftrajptr,plane.thickness_,extrapPlanes.interTolerance(),kkmat_h->applyBetheCorrection()); + ktrk.addMaterialPlaneXing(matxingptr,tdir); + processed.insert(inter.plane_); + found = true; + break; + } + if(!found) break; + } while(extrapPlanes.intersections().size()>0); + } + template void KKExtrap::extrapolateCRV(KKTrack& ktrk,TimeDir tdir) const { using KKCRVXING = KKShellXing; - using KKCRVXINGPTR = std::shared_ptr; + using KKMATRECXING = KKShellXing; // strongback (passive) plane + using MaterialPlane = ExtrapolateCRVRegion::MaterialPlane; GeomHandle kkg_h; GeomHandle kkmat_h; - // extrapolate to the extracted CRV. Loop to cover multiple intersections - auto extrapCRV = ExtrapolateCRV(maxdt_,maxdtstep_,btol_,intertol_,minv_,*kkg_h->CRV(),debug_); - if(debug_ > 5){std::cout << "Extrapolating to CRV with " << extrapCRV.sectors().size() << " sectors" << std::endl; - for(auto const& sector : kkg_h->CRV()->sectors()) { + // The CRV sector planes and their aluminium strongbacks are spatially interleaved (the sector plane + // is modeled at the scintillator-stack center; its tracker-facing strongback plane sits ~60 mm -- + // a stack half-thickness -- toward the tracker), so they MUST be crossed in ONE outward, + // time-ordered pass -- a xing can only be added at the trajectory front/back, never the interior. + // Crossing all strongbacks first (the old extrapolatePassiveMaterial pass) advanced the front past + // the lower sectors of a stacked geometry (extracted EX/T1/T2 ~150 mm apart in y) and stranded them + // (CRV_T1 lost on every track). Hand both surface types to ExtrapolateCRVRegion, which stops at the + // next crossing of either; add them one at a time, innermost first, re-extrapolating between each + // so every crossing is added at the current front (mirrors extrapolatePassiveMaterial). Run-2 sees + // <=1 sector per direction, so the order (strongback then sector) is unchanged. + ExtrapolateCRVRegion::MaterialPlaneCollection strongbacks; + SurfaceId const sbid(SurfaceIdEnum::CRV_StrongBack); + for(auto const& plane : kkg_h->passiveMaterialPlanes()) + if(plane.sid_.id() == sbid.id()) strongbacks.push_back(plane); + auto extrapCRV = ExtrapolateCRVRegion(maxdt_,maxdtstep_,btol_,intertol_,minv_,minvlong_,*kkg_h->CRV(),strongbacks,debug_); + if(debug_ > 5){std::cout << "Extrapolating to CRV with " << extrapCRV.sectors().size() << " sectors + " + << strongbacks.size() << " strongbacks" << std::endl; + for(auto const& sector : kkg_h->CRV()->sectors()) std::cout << sector.sname_ << " position " << sector.sector_->center() << " halfwidth " << sector.whw_ << std::endl; - } } auto const& ftraj = ktrk.fitTraj(); + std::unordered_set processedSectors; + std::unordered_set processedPlanes; do { - // iterate until we no longer hit CRV modules ktrk.extrapolate(tdir,extrapCRV); - if(debug_ > 0) std::cout << "Found " << extrapCRV.intersections().size() << " CRV intersections " << std::endl; - for(auto const& inter : extrapCRV.intersections()){ - // we have a good intersection. Use this to create a Shell material Xing + if(debug_ > 0) std::cout << "Found " << extrapCRV.intersections().size() << " CRV-region intersections " << std::endl; + bool added = false; + // intersections are time-sorted (innermost first for this tdir): add the nearest un-processed + // one, then re-extrapolate so the next is again at the front. + for(auto const& xing : extrapCRV.intersections()){ auto const& reftrajptr = tdir == TimeDir::backwards ? ftraj.frontPtr() : ftraj.backPtr(); - auto crvxingptr = std::make_shared(extrapCRV.sector((size_t)inter.isect_).sector_, - SurfaceId(kkg_h->CRV()->sectorName(inter.isect_)),*kkmat_h->CRVMaterial(),inter.inter_,reftrajptr, - 2*inter.whw_, extrapCRV.interTolerance(),kkmat_h->applyBetheCorrection()); - ktrk.addCRVXing(crvxingptr,tdir); - if(debug_ > 1) std::cout << "Good CRV " << inter.inter_ << ftraj.range() << std::endl; + if(xing.isSector_){ + if(processedSectors.count(xing.sectorIdx_) > 0) continue; + auto crvxingptr = std::make_shared(extrapCRV.sector((size_t)xing.sectorIdx_).sector_, + SurfaceId(kkg_h->CRV()->sectorName(xing.sectorIdx_)),*kkmat_h->CRVMaterial(),xing.inter_,reftrajptr, + 2*xing.whw_, extrapCRV.interTolerance(),kkmat_h->applyBetheCorrection()); + ktrk.addCRVXing(crvxingptr,tdir); + processedSectors.insert(xing.sectorIdx_); + if(debug_ > 1) std::cout << "Good CRV sector " << xing.sectorIdx_ << " " << xing.inter_ << std::endl; + } else { + if(processedPlanes.count(xing.plane_) > 0) continue; + auto const& plane = *xing.plane_; + auto matxingptr = std::make_shared(plane.surface_,plane.sid_,*kkmat_h->material(plane.material_), + xing.inter_,reftrajptr,plane.thickness_,extrapCRV.interTolerance(),kkmat_h->applyBetheCorrection()); + ktrk.addMaterialPlaneXing(matxingptr,tdir); + processedPlanes.insert(xing.plane_); + if(debug_ > 1) std::cout << "Good CRV strongback " << plane.sid_ << " " << xing.inter_ << std::endl; + } + added = true; + break; } + if(!added) break; } while(extrapCRV.intersections().size()>0); } diff --git a/Mu2eKinKal/inc/KKFit.hh b/Mu2eKinKal/inc/KKFit.hh index 3b4dee0c0e..1dfa09a039 100644 --- a/Mu2eKinKal/inc/KKFit.hh +++ b/Mu2eKinKal/inc/KKFit.hh @@ -164,6 +164,7 @@ namespace mu2e { double addStrawMinDz_; int strawNBuffer_; bool saveHitCalib_; + bool savematxings_; // save generic (post-extrapolation) passive material Xings in the KalSeed SurfaceIdCollection ssids_; }; @@ -198,7 +199,8 @@ namespace mu2e { skipStrawCheck_(fitconfig.skipStrawCheck()), addStrawMinDz_(fitconfig.addStrawMinDz()), strawNBuffer_(fitconfig.strawNBuffer()), - saveHitCalib_(fitconfig.saveHitCalib()) + saveHitCalib_(fitconfig.saveHitCalib()), + savematxings_(fitconfig.saveMatXings()) { if (fitconfig.saveTraj() == "T0") { savetraj_ = t0seg; @@ -890,6 +892,23 @@ namespace mu2e { crvxing->materialEffects(dmom,paramomvar,perpmomvar); inters.emplace_back(ktraj.stateEstimate(crvxing->time()),XYZVectorF(ktraj.bnom()),crvxing->surfaceId(),crvxing->intersection(),dmom); } + // generic passive material Xings are only filled during extrapolation; saving them in the KalSeed is optional + if(savematxings_){ + for(auto const& matxing : kktrk.materialCylXings()){ + double stime = matxing->time() - epsilon; + auto const& ktraj = ptraj.nearestPiece(stime); + double dmom,paramomvar,perpmomvar; + matxing->materialEffects(dmom,paramomvar,perpmomvar); + inters.emplace_back(ktraj.stateEstimate(matxing->time()),XYZVectorF(ktraj.bnom()),matxing->surfaceId(),matxing->intersection(),dmom); + } + for(auto const& matxing : kktrk.materialPlaneXings()){ + double stime = matxing->time() - epsilon; + auto const& ktraj = ptraj.nearestPiece(stime); + double dmom,paramomvar,perpmomvar; + matxing->materialEffects(dmom,paramomvar,perpmomvar); + inters.emplace_back(ktraj.stateEstimate(matxing->time()),XYZVectorF(ktraj.bnom()),matxing->surfaceId(),matxing->intersection(),dmom); + } + } // record other intersections saved in the track for(auto const& interpair : kktrk.intersections()) { auto const& sid = std::get<0>(interpair); diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 32959710c8..4f00d6952b 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -98,6 +98,7 @@ namespace mu2e { fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; fhicl::Atom sampleInRange { Name("SampleInRange"), Comment("Require sample times to be inside the fit trajectory time range") }; fhicl::Atom sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") }; + fhicl::Atom saveMatXings { Name("SaveMaterialXings"), Comment("Save generic passive material Xings (filled during extrapolation) in the KalSeed intersections"), true }; }; // struct for configuring a KinKal fit module struct KKModuleConfig { @@ -117,6 +118,11 @@ namespace mu2e { fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit (ns)") }; fhicl::Atom MaxDtStep { Name("MaxDtStep"), Comment("Maximum time step when extrapolating a fit (ns)") }; fhicl::Atom MinV { Name("MinV"), Comment("Minimum velocity perp to the surface extrapolate a straight track fit") }; + fhicl::Atom MinVLong { Name("MinVLong"), Comment("Minimum velocity perp to a z-normal CRV end-cap plane (mm/ns); suppresses grazing-cosmic fakes"), 265.0f }; + fhicl::Atom MaxTrackerEndsRadius { Name("MaxTrackerEndsRadius"), Comment("Max transverse radius rho=sqrt(x^2+y^2) (mm) a track may reach in the toTrackerEnds step"), 1.0e10f }; + fhicl::Atom ToTrackerEndsRadius { Name("ToTrackerEndsRadius"), Comment("In toTrackerEnds, extrapolate to the tracker outer cylinder (TT_Outer) rather than the z-end planes"), false }; + fhicl::Atom MaxTrackerEndsZ { Name("MaxTrackerEndsZ"), Comment("Max z (mm) reachable in toTrackerEnds (used only when ToTrackerEndsRadius=true)"), 1.0e10f }; + fhicl::Atom MinTrackerEndsZ { Name("MinTrackerEndsZ"), Comment("Min z (mm) reachable in toTrackerEnds (used only when ToTrackerEndsRadius=true)"), -1.0e10f }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; diff --git a/Mu2eKinKal/inc/KKTrack.hh b/Mu2eKinKal/inc/KKTrack.hh index 20f772ba59..3f6412726d 100644 --- a/Mu2eKinKal/inc/KKTrack.hh +++ b/Mu2eKinKal/inc/KKTrack.hh @@ -54,6 +54,12 @@ namespace mu2e { using KKCRVXING = KKShellXing; using KKCRVXINGPTR = std::shared_ptr; using KKCRVXINGCOL = std::vector; + using KKMATCYLXING = KKShellXing; + using KKMATCYLXINGPTR = std::shared_ptr; + using KKMATCYLXINGCOL = std::vector; + using KKMATRECXING = KKShellXing; // passive plane: Rectangle or Annulus + using KKMATRECXINGPTR = std::shared_ptr; + using KKMATRECXINGCOL = std::vector; using KKINTER = std::tuple; using KKINTERCOL = std::vector; using TRACK = KinKal::Track; @@ -73,8 +79,8 @@ namespace mu2e { shclusterer_(rhs.strawHitClusterer()), strawhits_(rhs.strawHits()), strawxings_(rhs.strawXings()), - inters_(rhs.intersections()), - calohits_(rhs.caloHits()){ + calohits_(rhs.caloHits()), + inters_(rhs.intersections()){ // hits and crossings were reallocated into the base copy; here, // we propagate references to those reallocations to the subclass this->remap_pointer_collection(strawhits_, this->hits(), context); @@ -85,6 +91,8 @@ namespace mu2e { this->clone_pointer_collection(ipaxings_, rhs.IPAXings(), context); this->clone_pointer_collection(stxings_, rhs.STXings(), context); this->clone_pointer_collection(crvxings_, rhs.CRVXings(), context); + this->clone_pointer_collection(matcylxings_, rhs.materialCylXings(), context); + this->clone_pointer_collection(matrecxings_, rhs.materialPlaneXings(), context); this->clone_pointer_collection(strawhitclusters_, rhs.strawHitClusters(), context); } @@ -124,6 +132,10 @@ namespace mu2e { KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits); // extend the track to cover a new set of material Xings. This will reuse the existing config object void extendTrack(EXINGCOL const& xings); + // add passive cylindrical material Xing + void addMaterialCylXing(KKMATCYLXINGPTR const& matxing,TimeDir const& tdir); + // add passive rectangular material Xing + void addMaterialPlaneXing(KKMATRECXINGPTR const& matxing,TimeDir const& tdir); // add IPA Xing void addIPAXing(KKIPAXINGPTR const& ipaxing,TimeDir const& tdir); // add ST Xing @@ -142,6 +154,8 @@ namespace mu2e { KKIPAXINGCOL const& IPAXings() const { return ipaxings_; } KKSTXINGCOL const& STXings() const { return stxings_; } KKCRVXINGCOL const& CRVXings() const { return crvxings_; } + KKMATCYLXINGCOL const& materialCylXings() const { return matcylxings_; } + KKMATRECXINGCOL const& materialPlaneXings() const { return matrecxings_; } KKINTERCOL const& intersections() const { return inters_; } KKCALOHITCOL const& caloHits() const { return calohits_; } void printFit(std::ostream& ost=std::cout,int detail=0) const; @@ -161,6 +175,8 @@ namespace mu2e { KKIPAXINGCOL ipaxings_; // ipa material crossings used in extrapolation KKSTXINGCOL stxings_; // stopping target material crossings used in extrapolation KKCRVXINGCOL crvxings_; // crv crossings using in extrapolation + KKMATCYLXINGCOL matcylxings_; // passive cylindrical material crossings used in extrapolation + KKMATRECXINGCOL matrecxings_; // passive planar material crossings used in extrapolation KKINTERCOL inters_; // other recorded intersections KKSTRAWHITCLUSTERCOL strawhitclusters_; // straw hit clusters used in this fit // utility function to convert to generic types @@ -296,6 +312,30 @@ namespace mu2e { this->extend(this->config(),nohits,exings); } + template void KKTrack::addMaterialCylXing(KKMATCYLXINGPTR const& matxingptr,TimeDir const& tdir) { + // convert to a generic Xing + std::shared_ptr> exptr = std::static_pointer_cast>(matxingptr); + // extrapolate the fit through this xing + if(!this->extrapolate(exptr,tdir))throw cet::exception("RECO")<<"mu2e::KKTrack: cylindrical material shell extrapolation failure" + << " sid " << matxingptr->surfaceId() << " time " << matxingptr->time() + << " normDot " << matxingptr->intersection().norm_.Dot(matxingptr->intersection().pdir_) + << " active " << matxingptr->active() << std::endl; + // store the xing + matcylxings_.push_back(matxingptr); + } + + template void KKTrack::addMaterialPlaneXing(KKMATRECXINGPTR const& matxingptr,TimeDir const& tdir) { + // convert to a generic Xing + std::shared_ptr> exptr = std::static_pointer_cast>(matxingptr); + // extrapolate the fit through this xing + if(!this->extrapolate(exptr,tdir))throw cet::exception("RECO")<<"mu2e::KKTrack: passive material plane extrapolation failure" + << " sid " << matxingptr->surfaceId() << " time " << matxingptr->time() + << " normDot " << matxingptr->intersection().norm_.Dot(matxingptr->intersection().pdir_) + << " active " << matxingptr->active() << std::endl; + // store the xing + matrecxings_.push_back(matxingptr); + } + template void KKTrack::printFit(std::ostream& ost,int printlevel) const { TRACK::print(ost,0); ost << "Fit with " << strawhits_.size() << " StrawHits and " << calohits_.size() << " CaloHits and " << strawxings_.size() << " Straw Xings" << std::endl; diff --git a/Mu2eKinKal/src/KKBField.cc b/Mu2eKinKal/src/KKBField.cc index 5623968684..08dc93c57a 100644 --- a/Mu2eKinKal/src/KKBField.cc +++ b/Mu2eKinKal/src/KKBField.cc @@ -12,9 +12,12 @@ namespace mu2e { // = bfmgr_.getBField(vpoint_mu2e); if(bfmgr_.getBFieldWithStatus(vpoint_mu2e,field)) return VEC3(field); -// see if there's no maps; that says this is the no-field case -// FIXME need to deal with case when outside all maps - static const VEC3 nullfield(0.0,0.0,0.0); + // Outside all field maps the real field is ~0 (detector hall, far from the solenoids), but a + // CentralHelix evaluated at IDENTICALLY zero field is singular (omega->0 -> 1/omega -> NaN). A + // bfcorr=true fit's domain walk can transiently sample just outside the maps, so return a tiny + // NON-zero axial placeholder instead: physically negligible, but it keeps omega finite so the + // fit can't NaN. + static const VEC3 nullfield(0.0,0.0,1.0e-3); return nullfield; } @@ -34,7 +37,9 @@ namespace mu2e { // numerical derivatives for now: TODO! VEC3 KKBField::fieldDeriv(VEC3 const& position, VEC3 const& velocity) const { static double dt(0.1); // 100 psec, ~3cm. this is arbitrary and should be set according to the field sampling TODO - if(!inRange(position))throw std::runtime_error("Position out of range"); + // fieldVect is defined everywhere (tiny non-zero placeholder outside the maps), so the derivative + // can always be evaluated; out of the maps it is ~0. No out-of-range throw: a bfcorr=true fit whose + // domain walk steps just outside the maps must not be aborted. VEC3 start = fieldVect(position); VEC3 end = fieldVect(position + velocity*dt); return (end-start)/dt; From adeba8e0f664b69676728099c22e61b6d428dc9b Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Sat, 11 Jul 2026 17:08:51 -0500 Subject: [PATCH 2/6] Make passive material xing off by default. --- Mu2eKinKal/fcl/prolog.fcl | 2 +- Mu2eKinKal/inc/KKFitSettings.hh | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 53ab2f1536..df42b8f6f7 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -39,7 +39,7 @@ Mu2eKinKal : { SampleSurfaces : [] # specific to the fit type SaveHitCalibInfo : false SaveDomains : false - SaveMaterialXings : true # save generic passive material Xings (filled during extrapolation) in the KalSeed + SaveMaterialXings : false # save generic passive material Xings (filled during extrapolation) in the KalSeed } CHSEEDFIT: { diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 4f00d6952b..f4e23ffafc 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -98,7 +98,7 @@ namespace mu2e { fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; fhicl::Atom sampleInRange { Name("SampleInRange"), Comment("Require sample times to be inside the fit trajectory time range") }; fhicl::Atom sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") }; - fhicl::Atom saveMatXings { Name("SaveMaterialXings"), Comment("Save generic passive material Xings (filled during extrapolation) in the KalSeed intersections"), true }; + fhicl::Atom saveMatXings { Name("SaveMaterialXings"), Comment("Save generic passive material Xings (filled during extrapolation) in the KalSeed intersections"), false }; }; // struct for configuring a KinKal fit module struct KKModuleConfig { From e0873d0d67b981b5d28e91339e35469a690bed7f Mon Sep 17 00:00:00 2001 From: robmina Date: Mon, 13 Jul 2026 11:20:43 -0400 Subject: [PATCH 3/6] Add clarifying comments to address AI review issues. --- Mu2eKinKal/fcl/prolog.fcl | 4 ++-- Mu2eKinKal/inc/ExtrapolateCylinders.hh | 3 ++- Mu2eKinKal/inc/KKExtrap.hh | 20 +++++++++++++++++--- Mu2eKinKal/src/KKBField.cc | 11 +++++++++-- 4 files changed, 30 insertions(+), 8 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index df42b8f6f7..3cd6c01466 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -402,7 +402,7 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns - MinVLong : 265.0 # mm/ns; near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics + MinVLong : 265.0 # mm/ns (= 0.884c -> cos(incidence)>~0.88, i.e. within ~28deg of normal); near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction ToTrackerEnds : false Upstream : false @@ -416,7 +416,7 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns - MinVLong : 265.0 # mm/ns; near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics + MinVLong : 265.0 # mm/ns (= 0.884c -> cos(incidence)>~0.88, i.e. within ~28deg of normal); near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction ToTrackerEnds : false Upstream : false diff --git a/Mu2eKinKal/inc/ExtrapolateCylinders.hh b/Mu2eKinKal/inc/ExtrapolateCylinders.hh index f24483d3bc..d9086f3867 100644 --- a/Mu2eKinKal/inc/ExtrapolateCylinders.hh +++ b/Mu2eKinKal/inc/ExtrapolateCylinders.hh @@ -65,7 +65,8 @@ namespace mu2e { double dptol_ = 1e10; double intertol_ = 1e10; double minvnorm_ = 1e-5; - double minShellCosine_ = 0.1; + double minShellCosine_ = 0.1; // min |cos(incidence)| to accept a shell crossing (~within 84 deg of normal); + // below this the thin-shell path length blows up (near-tangent), so skip it MaterialCylinderCollection const& cylinders_; int debug_ = 0; double maxRadius_ = 0.0; diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 829d2e20d2..3fa36e27ca 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -54,8 +54,13 @@ namespace mu2e { double btol_, intertol_, maxdt_, maxdtstep_, minv_, minvlong_, maxtrackerendsradius_; double mintrackerendsz_, maxtrackerendsz_; bool backToTracker_, extrapolateOPA_, toTrackerEnds_, toTrackerEndsRadius_, upstream_, toCRV_; - double ipathick_ = 0.511; // ipa thickness: should come from geometry service TODO - double stthick_ = 0.1056; // st foil thickness: should come from geometry service TODO + // IPA/ST thicknesses are nominal literals, consistent with the rest of KinKalGeom, whose IPA + // (single Cylinder) and ST (uniform Annulus foils) surfaces are themselves hard-coded nominal + // placeholders (see KinKalGeomMaker::makeDS/makeTarget). Sourcing these from the real geometry + // service is part of the broader "build KinKalGeom from the geometry service" work, not a standalone + // fix -- doing it for the thicknesses alone while the surfaces stay nominal would be inconsistent. TODO + double ipathick_ = 0.511; // ipa thickness (mm) + double stthick_ = 0.1056; // st foil thickness (mm) }; KKExtrap::KKExtrap(KKExtrapConfig const& extrapconfig) : @@ -339,6 +344,12 @@ namespace mu2e { auto extrapCylinders = ExtrapolateCylinders(maxdt_,maxdtstep_,btol_,intertol_,minv_,cylinders,debug_); auto const& ftraj = ktrk.fitTraj(); using MaterialCylinder = ExtrapolateCylinders::MaterialCylinder; + // NB: this bank-all-sorted idiom differs deliberately from the one-crossing-at-a-time idiom in + // extrapolatePassiveMaterial / extrapolateCRV below. It is safe HERE because the DS shells are + // CONCENTRIC and closely spaced: one step overshoots the whole cluster and the resulting trajectory + // piece spans all of them, so all crossings can be added in a single sorted pass. The planar surfaces + // are STACKED and spatially separated, so adding more than the nearest crossing before re-extrapolating + // advances the front past the farther planes and strands them -- hence they must go one at a time. std::unordered_set banked; // bank each cylinder once (avoid re-found duplicates) do { ktrk.extrapolate(tdir,extrapCylinders); @@ -377,7 +388,10 @@ namespace mu2e { auto extrapPlanes = ExtrapolatePlanes(maxdt_,maxdtstep_,btol_,intertol_,minv_,planes,debug_); auto const& ftraj = ktrk.fitTraj(); - // Normal path: process each plane at most once. + // Normal path: process each plane at most once, ONE crossing per extrapolate pass (add nearest, + // re-extrapolate, repeat). Unlike extrapolateDSMaterial's bank-all-sorted pass, the passive planes + // are stacked and spatially separated, so adding more than the nearest crossing before re-extrapolating + // would advance the front past the farther planes and strand them (see the cross-reference note there). // Stop before calling ktrk.extrapolate when all planes have been processed so // we don't advance the trajectory past the last plane (which would overshoot the // CRV and prevent extrapolateCRV from finding it). diff --git a/Mu2eKinKal/src/KKBField.cc b/Mu2eKinKal/src/KKBField.cc index 08dc93c57a..b69eb3e7c9 100644 --- a/Mu2eKinKal/src/KKBField.cc +++ b/Mu2eKinKal/src/KKBField.cc @@ -16,7 +16,11 @@ namespace mu2e { // CentralHelix evaluated at IDENTICALLY zero field is singular (omega->0 -> 1/omega -> NaN). A // bfcorr=true fit's domain walk can transiently sample just outside the maps, so return a tiny // NON-zero axial placeholder instead: physically negligible, but it keeps omega finite so the - // fit can't NaN. + // fit can't NaN. 1e-3 T is negligible for every fit type: the implied curvature radius + // (R = p/(0.3 B) ~ 3 km for a 1 GeV/c track) dwarfs the ~10 m out-of-map path we extrapolate, and + // KKLine (straight) ignores the field entirely -- only CentralHelix's 1/omega ever sees it. + // TODO(KinKal): remove this placeholder once the upstream bfcorr singularity at B==0 is fixed; + // then out-of-map queries can return an honest (0,0,0). static const VEC3 nullfield(0.0,0.0,1.0e-3); return nullfield; } @@ -39,7 +43,10 @@ namespace mu2e { static double dt(0.1); // 100 psec, ~3cm. this is arbitrary and should be set according to the field sampling TODO // fieldVect is defined everywhere (tiny non-zero placeholder outside the maps), so the derivative // can always be evaluated; out of the maps it is ~0. No out-of-range throw: a bfcorr=true fit whose - // domain walk steps just outside the maps must not be aborted. + // domain walk steps just outside the maps must not be aborted. Dropping the throw does not hide a + // genuinely broken position -- one that is nonsensical still fails downstream as a bad surface + // intersection / extrapolation failure, which is caught there; the throw here only aborted the + // legitimate just-outside-the-maps case. VEC3 start = fieldVect(position); VEC3 end = fieldVect(position + velocity*dt); return (end-start)/dt; From 10ae3ba4220c8d66faf310bc329c5d3380fe5539 Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Mon, 13 Jul 2026 16:30:12 -0500 Subject: [PATCH 4/6] Replace ExtrapolateToR with ExtrapolateToTrackerPerimeter, remove now-unused max-R configuration for ExtrapolateToZ. --- Mu2eKinKal/inc/ExtrapolateToR.hh | 75 ---------- .../inc/ExtrapolateToTrackerPerimeter.hh | 76 ++++++++++ Mu2eKinKal/inc/ExtrapolateToZ.hh | 15 +- Mu2eKinKal/inc/KKExtrap.hh | 134 ++++++++---------- Mu2eKinKal/inc/KKFitSettings.hh | 4 - 5 files changed, 137 insertions(+), 167 deletions(-) delete mode 100644 Mu2eKinKal/inc/ExtrapolateToR.hh create mode 100644 Mu2eKinKal/inc/ExtrapolateToTrackerPerimeter.hh diff --git a/Mu2eKinKal/inc/ExtrapolateToR.hh b/Mu2eKinKal/inc/ExtrapolateToR.hh deleted file mode 100644 index 74f51ea7dd..0000000000 --- a/Mu2eKinKal/inc/ExtrapolateToR.hh +++ /dev/null @@ -1,75 +0,0 @@ -// predicate for extrapolation to a given transverse radius rho = sqrt(x^2+y^2), bounded in z. -// This is the radial analog of ExtrapolateToZ: it targets a constant-rho (cylindrical) surface rather -// than a constant-z plane. That is the natural tracker exit for a CentralHelix (cosmic) track, which -// traverses the tracker radially (mostly in y) and leaves through the cylindrical side, NOT through the -// z-end planes. The [minz,maxz] bounds keep a near-horizontal (large v_z) cosmic from chasing the target -// radius far down the beamline, past the tracker ends -- the symmetric guard to ExtrapolateToZ's behaviour -// for near-vertical cosmics. -#ifndef Mu2eKinKal_ExtrapolateToR_hh -#define Mu2eKinKal_ExtrapolateToR_hh -#include "KinKal/Trajectory/ParticleTrajectory.hh" -#include "KinKal/General/TimeDir.hh" -#include -#include -#include -#include "cetlib_except/exception.h" -namespace mu2e { - using KinKal::TimeDir; - using KinKal::timeDirSign; - class ExtrapolateToR { - public: - ExtrapolateToR(double maxdt, double maxdtstep, double tol, double rval, double minz, double maxz, int debug) : - maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), rval_(rval), minz_(minz), maxz_(maxz), debug_(debug) {} - // interface for extrapolation - double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement - double maxDtStep() const { return maxDtStep_; } - double dpTolerance() const { return dptol_; } // tolerance on fractional momentum change - double rVal() const { return rval_; } - double minZ() const { return minz_; } - double maxZ() const { return maxz_; } - int debug() const { return debug_; } - // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time - template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; - private: - double maxDt_ = -1; // maximum extrapolation time - double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration - double dptol_ = 1e10; // fractional momentum tolerance in BField domain - double rval_ = 0; // target transverse radius (mm) - double minz_ = -std::numeric_limits::max(); // minimum z allowed during extrapolation (mm) - double maxz_ = std::numeric_limits::max(); // maximum z allowed during extrapolation (mm) - int debug_ = 0; // debug level - }; - - template bool ExtrapolateToR::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { - auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); - auto time = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); - auto vel = ktraj.velocity(time); - auto pos = ktraj.position3(time); - double rho = std::sqrt(pos.X()*pos.X() + pos.Y()*pos.Y()); - double z = pos.Z(); - // radial velocity, signed by the extrapolation direction (positive = heading outward) - double rvel = (rho > 0.0 ? (pos.X()*vel.X() + pos.Y()*vel.Y())/rho : 0.0)*timeDirSign(tdir); - double vz = vel.Z()*timeDirSign(tdir); // z velocity, signed by the extrapolation direction - if(debug_ > 2)std::cout << "R extrap start time " << time << " rho " << rho << " rvel " << rvel << " z " << z - << " vz " << vz << " rval " << rval_ << " z[" << minz_ << "," << maxz_ << "]" << std::endl; - // already outside the z bounds, or already at/beyond the target radius: stop - if(z > maxz_ || z < minz_) return false; - if(rho >= rval_) return false; - // heading inward: this direction will never reach the outer target radius - if(rvel <= 0.0) return false; - // Decide GEOMETRICALLY which boundary the track reaches first (the BField domain step out here is set - // by the field tolerance, not MaxDtStep, and can be large -- so a predicate that only stops AFTER the - // step overshoots the z bound by a whole domain). A near-axial cosmic (large |v_z|, small v_rho) - // reaches the z end before the target radius; extending it toward the radius is exactly what overshoots - // the z end, carrying rho far past the DS shells and stranding them as interior pieces for the material - // pass. In that case stop at the tracker boundary now and let the DS-material pass (no z bound, marches - // shell-by-shell) extend it. Otherwise the track exits radially: keep going to the target radius. - double time_to_rval = (rval_ - rho)/rvel; - double time_to_zbound = std::numeric_limits::max(); - if(vz > 0.0) time_to_zbound = (maxz_ - z)/vz; - else if(vz < 0.0) time_to_zbound = (z - minz_)/(-vz); - if(time_to_zbound < time_to_rval) return false; // exits via the z end: stop at the tracker boundary - return true; // radial exit: keep extending toward the target radius - } -} -#endif diff --git a/Mu2eKinKal/inc/ExtrapolateToTrackerPerimeter.hh b/Mu2eKinKal/inc/ExtrapolateToTrackerPerimeter.hh new file mode 100644 index 0000000000..797ae5f32a --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateToTrackerPerimeter.hh @@ -0,0 +1,76 @@ +// Predicate for extrapolation to the tracker PERIMETER: the closed cylinder bounding the tracker +// active volume, {rho = sqrt(x^2+y^2) <= R} intersected with {zmin <= z <= zmax}. The track is +// extended until it reaches the FIRST bounding face it meets -- the curved outer wall (constant rho) +// OR either z-end plane (constant z) -- whichever comes first for that track. +// +// This is the direction-agnostic generalization of the two former predicates: ExtrapolateToZ (a +// downstream track exits through a constant-z end plane) and ExtrapolateToR (a cosmic traverses the +// tracker radially -- mostly in y -- and exits through the curved outer wall). One predicate handles +// both, choosing the face per track from its own kinematics, so no configuration switch is needed. +// +// All three bounds (R, zmin, zmax) come from the tracker geometry (KinKalGeom::tracker(), itself +// sourced from GeometryService), so the tracker perimeter is defined in exactly one place -- no +// hand-entered radius/z limits that must be kept in sync with the geometry. +#ifndef Mu2eKinKal_ExtrapolateToTrackerPerimeter_hh +#define Mu2eKinKal_ExtrapolateToTrackerPerimeter_hh +#include "KinKal/Trajectory/ParticleTrajectory.hh" +#include "KinKal/General/TimeDir.hh" +#include +#include +#include +#include "cetlib_except/exception.h" +namespace mu2e { + using KinKal::TimeDir; + using KinKal::timeDirSign; + class ExtrapolateToTrackerPerimeter { + public: + ExtrapolateToTrackerPerimeter(double maxdt, double maxdtstep, double tol, + double rval, double zmin, double zmax, int debug) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), rval_(rval), zmin_(zmin), zmax_(zmax), debug_(debug) {} + // interface for extrapolation + double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } // tolerance on fractional momentum change + double rVal() const { return rval_; } + double zMin() const { return zmin_; } + double zMax() const { return zmax_; } + int debug() const { return debug_; } + // extrapolation predicate: the track will be extrapolated until this predicate returns false, subject to the maximum time + template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; + private: + double maxDt_ = -1; // maximum extrapolation time + double maxDtStep_ = -1; // maximum extrapolation time step in a single iteration + double dptol_ = 1e10; // fractional momentum tolerance in BField domain + double rval_ = 0; // outer wall transverse radius (mm) + double zmin_ = -std::numeric_limits::max(); // upstream z-end plane (mm) + double zmax_ = std::numeric_limits::max(); // downstream z-end plane (mm) + int debug_ = 0; // debug level + }; + + template bool ExtrapolateToTrackerPerimeter::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { + auto const& ktraj = tdir == TimeDir::forwards ? fittraj.back() : fittraj.front(); + auto time = tdir == TimeDir::forwards ? ktraj.range().end() : ktraj.range().begin(); + auto vel = ktraj.velocity(time); + auto pos = ktraj.position3(time); + double rho = std::sqrt(pos.X()*pos.X() + pos.Y()*pos.Y()); + double z = pos.Z(); + // velocities signed by the extrapolation direction (positive = heading toward that boundary) + double rvel = (rho > 0.0 ? (pos.X()*vel.X() + pos.Y()*vel.Y())/rho : 0.0)*timeDirSign(tdir); // outward radial + double vz = vel.Z()*timeDirSign(tdir); + if(debug_ > 2)std::cout << "Perimeter extrap start time " << time << " rho " << rho << " rvel " << rvel + << " z " << z << " vz " << vz << " R " << rval_ << " z[" << zmin_ << "," << zmax_ << "]" << std::endl; + // already at/outside any bounding face: stop (the recording step finds the precise crossing) + if(rho >= rval_) return false; + if(z <= zmin_ || z >= zmax_) return false; + // time to reach each face we are actually approaching (INF if heading away from it) + static const double INF = std::numeric_limits::max(); + double t_r = rvel > 0.0 ? (rval_ - rho)/rvel : INF; // outer wall + double t_zlo = vz < 0.0 ? (z - zmin_)/(-vz) : INF; // upstream z-end + double t_zhi = vz > 0.0 ? (zmax_ - z)/vz : INF; // downstream z-end + // inside the volume but heading toward no face (e.g. a stalled step): nothing more to reach + if(t_r == INF && t_zlo == INF && t_zhi == INF) return false; + // still strictly inside and closing on a face: keep extending + return true; + } +} +#endif diff --git a/Mu2eKinKal/inc/ExtrapolateToZ.hh b/Mu2eKinKal/inc/ExtrapolateToZ.hh index ffcf9d73a0..f4a4fc513f 100644 --- a/Mu2eKinKal/inc/ExtrapolateToZ.hh +++ b/Mu2eKinKal/inc/ExtrapolateToZ.hh @@ -5,16 +5,13 @@ #include "KinKal/Trajectory/ParticleTrajectory.hh" #include "KinKal/General/TimeDir.hh" #include -#include #include "cetlib_except/exception.h" namespace mu2e { using KinKal::TimeDir; using KinKal::timeDirSign; class ExtrapolateToZ { public: - ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double zval,int debug, - double maxradius = std::numeric_limits::max()) : - maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), zval_(zval), debug_(debug), maxradius_(maxradius) {} + ExtrapolateToZ(double maxdt, double maxdtstep, double tol, double zval,int debug) : maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(tol), zval_(zval), debug_(debug) {} // interface for extrapolation double maxDt() const { return maxDt_; } // maximum time to extend the track, WRT the time of the first(last) measurement double maxDtStep() const { return maxDtStep_; } @@ -29,9 +26,6 @@ namespace mu2e { double dptol_ = 1e10; // fractional momentum tolerance in BField domain double zval_ = 0; // z value targeted int debug_ = 0; // debug level - double maxradius_ = std::numeric_limits::max(); // max transverse radius rho=sqrt(x^2+y^2) (mm) - // the track may reach; stops near-vertical (small-v_z) cosmics from chasing the - // constant-z target plane out to large rho. Default = no limit. }; template bool ExtrapolateToZ::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, KinKal::TimeDir tdir) const { @@ -43,12 +37,7 @@ namespace mu2e { double zval = pos.Z(); auto const& bnom = ktraj.bnom(time); double zref = vel.R()*fabs(sin(bnom.Theta())); - double rho = std::sqrt(pos.X()*pos.X() + pos.Y()*pos.Y()); - if(debug_ > 2)std::cout << "Z extrap start time " << time << " z " << zval << " zvel " << zvel << " zref " << zref << " rho " << rho << std::endl; - // stop if we have reached the maximum allowed transverse radius: a near-vertical (small-v_z) track - // would otherwise keep chasing the constant-z target plane out to large rho (e.g. a cosmic running - // metres in y), pre-extending the trajectory past the DS before the DS-material pass. Default = no limit. - if(rho > maxradius_) return false; + if(debug_ > 2)std::cout << "Z extrap start time " << time << " z " << zval << " zvel " << zvel << " zref " << zref << std::endl; // if z velocity is unreliable, continue if(fabs(zvel) < zref) return true; // stop if we're heading away from the target z diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 3fa36e27ca..43e1150bfa 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -10,7 +10,7 @@ #include "Offline/KinKalGeom/inc/Tracker.hh" #include "Offline/Mu2eKinKal/inc/KKTrack.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateToZ.hh" -#include "Offline/Mu2eKinKal/inc/ExtrapolateToR.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateToTrackerPerimeter.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateIPA.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateST.hh" #include "Offline/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh" @@ -22,6 +22,10 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include #include +#include +#include +#include +#include #include "Offline/GeometryService/inc/GeometryService.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" namespace mu2e { @@ -47,13 +51,12 @@ namespace mu2e { template void extrapolatePassiveMaterial(KKTrack& ktrk,TimeDir tdir) const; template void extrapolateCRV(KKTrack& ktrk,TimeDir tdir) const; template void extrapolateOPA(KKTrack& ktrk, double tstart, TimeDir tdir) const; - template void toTrackerEnds(KKTrack& ktrk) const; + template void toTrackerPerimeter(KKTrack& ktrk) const; private: int debug_; - double btol_, intertol_, maxdt_, maxdtstep_, minv_, minvlong_, maxtrackerendsradius_; - double mintrackerendsz_, maxtrackerendsz_; - bool backToTracker_, extrapolateOPA_, toTrackerEnds_, toTrackerEndsRadius_, upstream_, toCRV_; + double btol_, intertol_, maxdt_, maxdtstep_, minv_, minvlong_; + bool backToTracker_, extrapolateOPA_, toTrackerEnds_, upstream_, toCRV_; // IPA/ST thicknesses are nominal literals, consistent with the rest of KinKalGeom, whose IPA // (single Cylinder) and ST (uniform Annulus foils) surfaces are themselves hard-coded nominal // placeholders (see KinKalGeomMaker::makeDS/makeTarget). Sourcing these from the real geometry @@ -71,13 +74,9 @@ namespace mu2e { maxdtstep_(extrapconfig.MaxDtStep()), minv_(extrapconfig.MinV()), minvlong_(extrapconfig.MinVLong()), - maxtrackerendsradius_(extrapconfig.MaxTrackerEndsRadius()), - mintrackerendsz_(extrapconfig.MinTrackerEndsZ()), - maxtrackerendsz_(extrapconfig.MaxTrackerEndsZ()), backToTracker_(extrapconfig.BackToTracker()), extrapolateOPA_(extrapconfig.ToOPA()), toTrackerEnds_(extrapconfig.ToTrackerEnds()), - toTrackerEndsRadius_(extrapconfig.ToTrackerEndsRadius()), upstream_(extrapconfig.Upstream()), toCRV_(extrapconfig.ToCRV()) {} @@ -88,7 +87,7 @@ namespace mu2e { auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; - if(toTrackerEnds_)toTrackerEnds(ktrk); + if(toTrackerEnds_)toTrackerPerimeter(ktrk); if(upstream_){ double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); // extrapolate through the IPA in this direction. @@ -129,77 +128,62 @@ namespace mu2e { } } - template void KKExtrap::toTrackerEnds(KKTrack& ktrk) const { + // Extrapolate a fit to the tracker PERIMETER (the closed cylinder bounding the tracker active volume) + // and record the bounding-surface intersection the track actually reaches. One code path serves every + // track type: a downstream track exits through a z-end plane (TT_Front/TT_Back), a cosmic exits through + // the outer wall (TT_Outer). The perimeter (R, z-ends) is taken entirely from KinKalGeom::tracker(), + // which is built from GeometryService -- there are no hand-entered radius/z limits to keep in sync. + template void KKExtrap::toTrackerPerimeter(KKTrack& ktrk) const { GeomHandle kkg_h; - - // CentralHelix (cosmic) mode: extrapolate to the tracker OUTER CYLINDER (constant rho), not the z-end - // planes. A cosmic traverses the tracker radially (mostly in y), so the outer cylinder is its true exit; - // ExtrapolateToZ would chase the constant-z plane out to large rho. The z bounds keep near-horizontal - // cosmics from chasing the radius past the tracker ends. - if(toTrackerEndsRadius_){ - auto const& outer = kkg_h->tracker()->outer(); - ExtrapolateToR toOuter(maxdt_,maxdtstep_,btol_,outer.radius(),mintrackerendsz_,maxtrackerendsz_,debug_); - auto const& ftraj = ktrk.fitTraj(); - // the cosmic crosses the outer cylinder on both ends; extend in both time directions - auto tofwd = ktrk.extrapolate(TimeDir::forwards,toOuter); - auto tobwd = ktrk.extrapolate(TimeDir::backwards,toOuter); - static const SurfaceId tt_outer("TT_Outer"); - static const SurfaceId tt_mid("TT_Mid"); - // middle reference, recorded only if the track actually crosses it - auto midinter = KinKal::intersect(ftraj,kkg_h->tracker()->middle(),ftraj.range(),intertol_); - if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); - // outer-cylinder crossings (tracker entry and exit) from the respective end pieces - if(tofwd){ - auto const& fhel = ftraj.back(); - auto inter = KinKal::intersect(fhel,outer,fhel.range(),intertol_,TimeDir::forwards); - if(inter.good()) ktrk.addIntersection(tt_outer,inter); - } - if(tobwd){ - auto const& bhel = ftraj.front(); - auto inter = KinKal::intersect(bhel,outer,bhel.range(),intertol_,TimeDir::backwards); - if(inter.good()) ktrk.addIntersection(tt_outer,inter); - } - return; - } - - ExtrapolateToZ trackerFront(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->front().center().Z(),debug_,maxtrackerendsradius_); - ExtrapolateToZ trackerBack(maxdt_,maxdtstep_,btol_,kkg_h->tracker()->back().center().Z(),debug_,maxtrackerendsradius_); - // time direction to reach the bounding surfaces from the active region depends on the z momentum. This calculation assumes the particle doesn't - // reflect inside the tracker volume + auto const& trk = *kkg_h->tracker(); + double const R = trk.outer().radius(); + double const zf = trk.front().center().Z(); + double const zb = trk.back().center().Z(); + double const zmin = std::min(zf,zb); // upstream z-end plane (do not assume front 0) ? TimeDir::backwards : TimeDir::forwards; - TimeDir backtdir = (dir0.Z() > 0) ? TimeDir::forwards : TimeDir::backwards; - auto tofront = ktrk.extrapolate(fronttdir,trackerFront); - auto toback = ktrk.extrapolate(backtdir,trackerBack); - // record the standard tracker intersections + static const SurfaceId tt_outer("TT_Outer"); static const SurfaceId tt_front("TT_Front"); - static const SurfaceId tt_mid("TT_Mid"); static const SurfaceId tt_back("TT_Back"); + static const SurfaceId tt_mid("TT_Mid"); - // start with the middle - auto midinter = KinKal::intersect(ftraj,kkg_h->tracker()->middle(),ftraj.range(),intertol_); + // middle reference, recorded once if the track crosses it + auto midinter = KinKal::intersect(ftraj,trk.middle(),ftraj.range(),intertol_); if(midinter.good()) ktrk.addIntersection(tt_mid,midinter); - if(tofront){ - // check the front piece first; that is usually correct - // track extrapolation to the front succeeded, but the intersection failed. Use the last trajectory to force an intersection - auto fhel = fronttdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); - auto frontinter = KinKal::intersect(fhel,kkg_h->tracker()->front(),fhel.range(),intertol_,fronttdir); - if(!frontinter.good()){ - // start from the middle - TimeRange frange = ftraj.range(); - if(midinter.good())frange = fronttdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - frontinter = KinKal::intersect(ftraj,kkg_h->tracker()->front(),frange,intertol_,fronttdir); - } - if(frontinter.good()) ktrk.addIntersection(tt_front,frontinter); - } - if(toback){ - // start from the middle - TimeRange brange = ftraj.range(); - if(midinter.good())brange = backtdir == TimeDir::forwards ? TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); - auto backinter = KinKal::intersect(ftraj,kkg_h->tracker()->back(),brange,intertol_,backtdir); - if(backinter.good())ktrk.addIntersection(tt_back,backinter); - } + + // The track meets the perimeter at BOTH ends (tracker entry and exit); extend each time direction and + // record the face it reaches. For a given direction the exit face is the bounding surface crossed + // first, so among the good candidate intersections we keep the earliest one in that direction. Disk + // (z-end) intersections can fail at the edge, so on failure retry over the middle->end sub-range using + // the full trajectory (the outer wall does not need this fallback). + auto recordFace = [&](TimeDir tdir){ + if(!ktrk.extrapolate(tdir,toPerim)) return; + auto const& end = tdir == TimeDir::forwards ? ftraj.back() : ftraj.front(); + double const sgn = KinKal::timeDirSign(tdir); + double const t0 = tdir == TimeDir::forwards ? end.range().begin() : end.range().end(); + std::optional> best; + double bestdt = std::numeric_limits::max(); + auto consider = [&](SurfaceId const& sid, auto const& surf, bool midFallback){ + auto inter = KinKal::intersect(end,surf,end.range(),intertol_,tdir); + if(!inter.good() && midFallback){ + TimeRange rng = ftraj.range(); + if(midinter.good()) rng = tdir == TimeDir::forwards ? + TimeRange(midinter.time_,ftraj.range().end()) : TimeRange(ftraj.range().begin(),midinter.time_); + inter = KinKal::intersect(ftraj,surf,rng,intertol_,tdir); + } + if(inter.good()){ + double const dt = sgn*(inter.time_ - t0); // time from the tracker interior out to this face + if(dt >= 0.0 && dt < bestdt){ bestdt = dt; best.emplace(sid,inter); } + } + }; + consider(tt_outer, trk.outer(), false); + consider(tt_front, trk.front(), true); + consider(tt_back, trk.back(), true); + if(best) ktrk.addIntersection(best->first,best->second); + }; + recordFace(TimeDir::backwards); + recordFace(TimeDir::forwards); } template bool KKExtrap::extrapolateIPA(KKTrack& ktrk,TimeDir tdir) const { diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index f4e23ffafc..a61bd059ab 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -119,10 +119,6 @@ namespace mu2e { fhicl::Atom MaxDtStep { Name("MaxDtStep"), Comment("Maximum time step when extrapolating a fit (ns)") }; fhicl::Atom MinV { Name("MinV"), Comment("Minimum velocity perp to the surface extrapolate a straight track fit") }; fhicl::Atom MinVLong { Name("MinVLong"), Comment("Minimum velocity perp to a z-normal CRV end-cap plane (mm/ns); suppresses grazing-cosmic fakes"), 265.0f }; - fhicl::Atom MaxTrackerEndsRadius { Name("MaxTrackerEndsRadius"), Comment("Max transverse radius rho=sqrt(x^2+y^2) (mm) a track may reach in the toTrackerEnds step"), 1.0e10f }; - fhicl::Atom ToTrackerEndsRadius { Name("ToTrackerEndsRadius"), Comment("In toTrackerEnds, extrapolate to the tracker outer cylinder (TT_Outer) rather than the z-end planes"), false }; - fhicl::Atom MaxTrackerEndsZ { Name("MaxTrackerEndsZ"), Comment("Max z (mm) reachable in toTrackerEnds (used only when ToTrackerEndsRadius=true)"), 1.0e10f }; - fhicl::Atom MinTrackerEndsZ { Name("MinTrackerEndsZ"), Comment("Min z (mm) reachable in toTrackerEnds (used only when ToTrackerEndsRadius=true)"), -1.0e10f }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; From 664b09ddd4deec42a7fdf03ee35deb5de8cb5f2f Mon Sep 17 00:00:00 2001 From: robmina Date: Mon, 13 Jul 2026 19:34:58 -0400 Subject: [PATCH 5/6] Revert KKBField.cc to upstream main. --- Mu2eKinKal/src/KKBField.cc | 20 ++++---------------- 1 file changed, 4 insertions(+), 16 deletions(-) diff --git a/Mu2eKinKal/src/KKBField.cc b/Mu2eKinKal/src/KKBField.cc index b69eb3e7c9..5623968684 100644 --- a/Mu2eKinKal/src/KKBField.cc +++ b/Mu2eKinKal/src/KKBField.cc @@ -12,16 +12,9 @@ namespace mu2e { // = bfmgr_.getBField(vpoint_mu2e); if(bfmgr_.getBFieldWithStatus(vpoint_mu2e,field)) return VEC3(field); - // Outside all field maps the real field is ~0 (detector hall, far from the solenoids), but a - // CentralHelix evaluated at IDENTICALLY zero field is singular (omega->0 -> 1/omega -> NaN). A - // bfcorr=true fit's domain walk can transiently sample just outside the maps, so return a tiny - // NON-zero axial placeholder instead: physically negligible, but it keeps omega finite so the - // fit can't NaN. 1e-3 T is negligible for every fit type: the implied curvature radius - // (R = p/(0.3 B) ~ 3 km for a 1 GeV/c track) dwarfs the ~10 m out-of-map path we extrapolate, and - // KKLine (straight) ignores the field entirely -- only CentralHelix's 1/omega ever sees it. - // TODO(KinKal): remove this placeholder once the upstream bfcorr singularity at B==0 is fixed; - // then out-of-map queries can return an honest (0,0,0). - static const VEC3 nullfield(0.0,0.0,1.0e-3); +// see if there's no maps; that says this is the no-field case +// FIXME need to deal with case when outside all maps + static const VEC3 nullfield(0.0,0.0,0.0); return nullfield; } @@ -41,12 +34,7 @@ namespace mu2e { // numerical derivatives for now: TODO! VEC3 KKBField::fieldDeriv(VEC3 const& position, VEC3 const& velocity) const { static double dt(0.1); // 100 psec, ~3cm. this is arbitrary and should be set according to the field sampling TODO - // fieldVect is defined everywhere (tiny non-zero placeholder outside the maps), so the derivative - // can always be evaluated; out of the maps it is ~0. No out-of-range throw: a bfcorr=true fit whose - // domain walk steps just outside the maps must not be aborted. Dropping the throw does not hide a - // genuinely broken position -- one that is nonsensical still fails downstream as a bad surface - // intersection / extrapolation failure, which is caught there; the throw here only aborted the - // legitimate just-outside-the-maps case. + if(!inRange(position))throw std::runtime_error("Position out of range"); VEC3 start = fieldVect(position); VEC3 end = fieldVect(position + velocity*dt); return (end-start)/dt; From a7602904c7d459861b8707bc2fdfcbd1ae5fdeb8 Mon Sep 17 00:00:00 2001 From: robmina Date: Thu, 16 Jul 2026 11:37:59 -0400 Subject: [PATCH 6/6] Remove MinVLong, consolidate toTrackerEnds -> toTrackerPerimeter, common ExtrapolateSurface template. --- Mu2eKinKal/fcl/prolog.fcl | 10 ++--- Mu2eKinKal/inc/ExtrapolateCRVRegion.hh | 24 +++++------- Mu2eKinKal/inc/ExtrapolateCylinders.hh | 46 +++++++---------------- Mu2eKinKal/inc/ExtrapolatePlanes.hh | 46 +++++++---------------- Mu2eKinKal/inc/ExtrapolateSurface.hh | 51 ++++++++++++++++++++++++++ Mu2eKinKal/inc/KKExtrap.hh | 11 +++--- Mu2eKinKal/inc/KKFitSettings.hh | 3 +- 7 files changed, 96 insertions(+), 95 deletions(-) create mode 100644 Mu2eKinKal/inc/ExtrapolateSurface.hh diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 3cd6c01466..77417868b2 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -377,7 +377,7 @@ Mu2eKinKal : { MaxDtStep : 1.0 # (ns) MinV : 0.0 # mm/ns BCorrTolerance : 1e-4 # momemntum fraction - ToTrackerEnds : true + ToTrackerPerimeter : true Upstream : true BackToTracker : false ToOPA : true @@ -390,7 +390,7 @@ Mu2eKinKal : { MaxDtStep : 10.0 # (ns) MinV : 0.0 # mm/ns BCorrTolerance : 1e-2 # momemntum fraction - ToTrackerEnds : true + ToTrackerPerimeter : true Upstream : false BackToTracker : false ToOPA : false @@ -402,9 +402,8 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns - MinVLong : 265.0 # mm/ns (= 0.884c -> cos(incidence)>~0.88, i.e. within ~28deg of normal); near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction - ToTrackerEnds : false + ToTrackerPerimeter : false Upstream : false BackToTracker : false ToOPA : false @@ -416,9 +415,8 @@ Mu2eKinKal : { MaxDt : 200.0 # (ns) MaxDtStep : 1.0 # (ns) MinV : 1e-5 # mm/ns - MinVLong : 265.0 # mm/ns (= 0.884c -> cos(incidence)>~0.88, i.e. within ~28deg of normal); near-normal-incidence floor for CRV end-cap (z-normal U/D/E) planes -> suppresses grazing fakes from soft transverse cosmics BCorrTolerance : 1e-2 # momemntum fraction - ToTrackerEnds : false + ToTrackerPerimeter : false Upstream : false BackToTracker : false ToOPA : false diff --git a/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh b/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh index 5aee23aed7..5cef2f85a1 100644 --- a/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh +++ b/Mu2eKinKal/inc/ExtrapolateCRVRegion.hh @@ -18,9 +18,10 @@ // RUN-2 SAFETY: run-2 geometry presents at most one CRV sector (and its one strongback) per // extrapolation direction, so the outward order is unchanged (strongback then sector) and no // stacking occurs. Concrete passive planes are NOT handled here (they stay in -// extrapolatePassiveMaterial); only the CRV strongbacks are merged with the sectors. Per-surface -// grazing floors are preserved exactly: sectors keep the z-normal end-cap floor (minvlong), the -// strongback planes keep minvnorm, identical to the originals. +// extrapolatePassiveMaterial); only the CRV strongbacks are merged with the sectors. Sectors and +// strongback planes share the same minvnorm floor (a near-zero guard against the normvel divide in +// time_to_sector, not a physics cut): grazing incidence is not pre-filtered by a velocity threshold, +// it is left to fail at the actual surface intersection (newinter.good()), same as any other surface. #ifndef Mu2eKinKal_ExtrapolateCRVRegion_hh #define Mu2eKinKal_ExtrapolateCRVRegion_hh #include "KinKal/General/TimeDir.hh" @@ -68,15 +69,9 @@ namespace mu2e { }; ExtrapolateCRVRegion(double maxdt, double maxdtstep, double dptol, double intertol, - double minv, double minvlong, CRV const& crv, MaterialPlaneCollection const& planes, int debug=0) : + double minv, CRV const& crv, MaterialPlaneCollection const& planes, int debug=0) : maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), - minvnorm_(minv), minvlong_(minvlong), debug_(debug), sectors_(crv.sectors()), planes_(planes) { - // per-sector grazing floor: stricter (minvlong) for z-normal end-cap planes, looser (minvnorm) - // for the oblique-cosmic top/side planes -- identical to ExtrapolateCRV. - sectorminv_.reserve(sectors_.size()); - for(auto const& sector : sectors_) - sectorminv_.push_back( fabs(sector.sector_->normal().Z()) > 0.9 ? minvlong_ : minvnorm_ ); - } + minvnorm_(minv), debug_(debug), sectors_(crv.sectors()), planes_(planes) {} double maxDt() const { return maxDt_; } double maxDtStep() const { return maxDtStep_; } double dpTolerance() const { return dptol_; } @@ -93,11 +88,9 @@ namespace mu2e { double dptol_ = 1e10; double intertol_ = 1e10; double minvnorm_ = 1e-5; - double minvlong_ = 1e-5; int debug_ = 0; CRVSV sectors_; MaterialPlaneCollection const& planes_; - std::vector sectorminv_; mutable RegionXingCollection inters_; }; @@ -124,8 +117,9 @@ namespace mu2e { double normvel = vel.Dot(sector.sector_->normal())*timeDirSign(tdir); double signed_perp = (sector.sector_->center()-pos).Dot(sector.sector_->normal()); double time_to_sector = signed_perp / normvel; - // skip grazing, or planes already behind this piece's evaluation point (already crossed/added) - if(fabs(normvel) < sectorminv_[isect] || time_to_sector < 0 )continue; + // skip planes already behind this piece's evaluation point (already crossed/added); grazing + // incidence is not pre-filtered here, it is left to fail at the intersection attempt below + if(fabs(normvel) < minvnorm_ || time_to_sector < 0 )continue; auto newinter = KinKal::intersect(fittraj,*sector.sector_,trange,intertol_,tdir); if(newinter.good()){ if(!coincident(newinter.time_)) inters_.emplace_back(newinter,(int)isect,sector.whw_); diff --git a/Mu2eKinKal/inc/ExtrapolateCylinders.hh b/Mu2eKinKal/inc/ExtrapolateCylinders.hh index d9086f3867..f64b5e6f03 100644 --- a/Mu2eKinKal/inc/ExtrapolateCylinders.hh +++ b/Mu2eKinKal/inc/ExtrapolateCylinders.hh @@ -7,6 +7,7 @@ #include "KinKal/Geometry/Intersection.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/KinKalGeom/inc/DetectorSolenoid.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateSurface.hh" #include #include @@ -20,57 +21,36 @@ namespace mu2e { using KinKal::TimeRange; using KinKal::timeDirSign; - class ExtrapolateCylinders { + struct CylinderIntersection { + using MaterialCylinder = KKGeom::DetectorSolenoid::MaterialCylinder; + Intersection inter_; + MaterialCylinder const* cylinder_ = nullptr; + CylinderIntersection() {} + CylinderIntersection(Intersection const& inter, MaterialCylinder const& cylinder) : + inter_(inter), cylinder_(&cylinder) {} + }; + + class ExtrapolateCylinders : public ExtrapolateSurface { public: + using Base = ExtrapolateSurface; using MaterialCylinder = KKGeom::DetectorSolenoid::MaterialCylinder; using MaterialCylinderCollection = KKGeom::DetectorSolenoid::MaterialCylinderCollection; - struct CylinderIntersection { - Intersection inter_; - MaterialCylinder const* cylinder_ = nullptr; - CylinderIntersection() {} - CylinderIntersection(Intersection const& inter, MaterialCylinder const& cylinder) : - inter_(inter), cylinder_(&cylinder) {} - }; - using CylinderIntersectionCollection = std::vector; - - struct SortIntersections { - TimeDir tdir_; - bool operator () (CylinderIntersection const& lhs, CylinderIntersection const& rhs) { - return tdir_ == TimeDir::forwards ? lhs.inter_.time_ < rhs.inter_.time_ : lhs.inter_.time_ > rhs.inter_.time_; - } - SortIntersections(TimeDir tdir) : tdir_(tdir) {} - }; ExtrapolateCylinders(double maxdt, double maxdtstep, double dptol, double intertol, double minv, MaterialCylinderCollection const& cylinders, int debug=0) : - maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), cylinders_(cylinders), debug_(debug) { + Base(maxdt,maxdtstep,dptol,intertol,minv,debug), cylinders_(cylinders) { for(auto const& cylinder : cylinders_) { maxRadius_ = std::max(maxRadius_, cylinder.surface_->radius() + 0.5*cylinder.thickness_); } } - double maxDt() const { return maxDt_; } - double maxDtStep() const { return maxDtStep_; } - double dpTolerance() const { return dptol_; } - double interTolerance() const { return intertol_; } - auto const& intersections() const { return inters_; } - int debug() const { return debug_; } - void reset() const { inters_.clear(); } - template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; private: - double maxDt_ = -1; - double maxDtStep_ = -1; - double dptol_ = 1e10; - double intertol_ = 1e10; - double minvnorm_ = 1e-5; double minShellCosine_ = 0.1; // min |cos(incidence)| to accept a shell crossing (~within 84 deg of normal); // below this the thin-shell path length blows up (near-tangent), so skip it MaterialCylinderCollection const& cylinders_; - int debug_ = 0; double maxRadius_ = 0.0; - mutable CylinderIntersectionCollection inters_; }; template bool ExtrapolateCylinders::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { diff --git a/Mu2eKinKal/inc/ExtrapolatePlanes.hh b/Mu2eKinKal/inc/ExtrapolatePlanes.hh index 05510595b8..d064853a99 100644 --- a/Mu2eKinKal/inc/ExtrapolatePlanes.hh +++ b/Mu2eKinKal/inc/ExtrapolatePlanes.hh @@ -7,6 +7,7 @@ #include "KinKal/Geometry/Intersection.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +#include "Offline/Mu2eKinKal/inc/ExtrapolateSurface.hh" #include #include @@ -19,50 +20,29 @@ namespace mu2e { using KinKal::TimeRange; using KinKal::timeDirSign; - class ExtrapolatePlanes { + struct PlaneIntersection { + using MaterialPlane = KinKalGeom::PassiveMaterialPlane; + Intersection inter_; + MaterialPlane const* plane_ = nullptr; + PlaneIntersection() {} + PlaneIntersection(Intersection const& inter, MaterialPlane const& plane) : + inter_(inter), plane_(&plane) {} + }; + + class ExtrapolatePlanes : public ExtrapolateSurface { public: + using Base = ExtrapolateSurface; using MaterialPlane = KinKalGeom::PassiveMaterialPlane; using MaterialPlaneCollection = KinKalGeom::PassiveMaterialPlaneCollection; - struct PlaneIntersection { - Intersection inter_; - MaterialPlane const* plane_ = nullptr; - PlaneIntersection() {} - PlaneIntersection(Intersection const& inter, MaterialPlane const& plane) : - inter_(inter), plane_(&plane) {} - }; - using PlaneIntersectionCollection = std::vector; - - struct SortIntersections { - TimeDir tdir_; - bool operator () (PlaneIntersection const& lhs, PlaneIntersection const& rhs) { - return tdir_ == TimeDir::forwards ? lhs.inter_.time_ < rhs.inter_.time_ : lhs.inter_.time_ > rhs.inter_.time_; - } - SortIntersections(TimeDir tdir) : tdir_(tdir) {} - }; ExtrapolatePlanes(double maxdt, double maxdtstep, double dptol, double intertol, double minv, MaterialPlaneCollection const& planes, int debug=0) : - maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), planes_(planes), debug_(debug) {} - - double maxDt() const { return maxDt_; } - double maxDtStep() const { return maxDtStep_; } - double dpTolerance() const { return dptol_; } - double interTolerance() const { return intertol_; } - auto const& intersections() const { return inters_; } - int debug() const { return debug_; } - void reset() const { inters_.clear(); } + Base(maxdt,maxdtstep,dptol,intertol,minv,debug), planes_(planes) {} template bool needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const; private: - double maxDt_ = -1; - double maxDtStep_ = -1; - double dptol_ = 1e10; - double intertol_ = 1e10; - double minvnorm_ = 1e-5; MaterialPlaneCollection const& planes_; - int debug_ = 0; - mutable PlaneIntersectionCollection inters_; }; template bool ExtrapolatePlanes::needsExtrapolation(KinKal::ParticleTrajectory const& fittraj, TimeDir tdir) const { diff --git a/Mu2eKinKal/inc/ExtrapolateSurface.hh b/Mu2eKinKal/inc/ExtrapolateSurface.hh new file mode 100644 index 0000000000..a01d5e0e43 --- /dev/null +++ b/Mu2eKinKal/inc/ExtrapolateSurface.hh @@ -0,0 +1,51 @@ +// Common scaffolding shared by the surface-crossing extrapolation predicates (ExtrapolateCylinders, +// ExtrapolatePlanes): dt/tolerance bookkeeping, intersection storage, and time-ordering. Each +// predicate's needsExtrapolation search is geometry-specific and NOT shared here: ExtrapolateCylinders +// marches in bounded sub-windows to avoid stepping over closely-spaced concentric shells, while +// ExtrapolatePlanes gates directly on normal velocity/distance over the whole piece. Forcing those two +// search strategies into one shared body would obscure both, so only the bookkeeping below is common. +#ifndef Mu2eKinKal_ExtrapolateSurface_hh +#define Mu2eKinKal_ExtrapolateSurface_hh + +#include "KinKal/General/TimeDir.hh" + +#include + +namespace mu2e { + using KinKal::TimeDir; + + template class ExtrapolateSurface { + public: + using SurfaceIntersectionCollection = std::vector; + + struct SortIntersections { + TimeDir tdir_; + bool operator () (SurfaceIntersection const& lhs, SurfaceIntersection const& rhs) { + return tdir_ == TimeDir::forwards ? lhs.inter_.time_ < rhs.inter_.time_ : lhs.inter_.time_ > rhs.inter_.time_; + } + SortIntersections(TimeDir tdir) : tdir_(tdir) {} + }; + + double maxDt() const { return maxDt_; } + double maxDtStep() const { return maxDtStep_; } + double dpTolerance() const { return dptol_; } + double interTolerance() const { return intertol_; } + auto const& intersections() const { return inters_; } + int debug() const { return debug_; } + void reset() const { inters_.clear(); } + + protected: + ExtrapolateSurface(double maxdt, double maxdtstep, double dptol, double intertol, double minv, int debug) : + maxDt_(maxdt), maxDtStep_(maxdtstep), dptol_(dptol), intertol_(intertol), minvnorm_(minv), debug_(debug) {} + + double maxDt_ = -1; + double maxDtStep_ = -1; + double dptol_ = 1e10; + double intertol_ = 1e10; + double minvnorm_ = 1e-5; + int debug_ = 0; + mutable SurfaceIntersectionCollection inters_; + }; +} + +#endif diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index 43e1150bfa..6491674ada 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -55,8 +55,8 @@ namespace mu2e { private: int debug_; - double btol_, intertol_, maxdt_, maxdtstep_, minv_, minvlong_; - bool backToTracker_, extrapolateOPA_, toTrackerEnds_, upstream_, toCRV_; + double btol_, intertol_, maxdt_, maxdtstep_, minv_; + bool backToTracker_, extrapolateOPA_, toTrackerPerimeter_, upstream_, toCRV_; // IPA/ST thicknesses are nominal literals, consistent with the rest of KinKalGeom, whose IPA // (single Cylinder) and ST (uniform Annulus foils) surfaces are themselves hard-coded nominal // placeholders (see KinKalGeomMaker::makeDS/makeTarget). Sourcing these from the real geometry @@ -73,10 +73,9 @@ namespace mu2e { maxdt_(extrapconfig.MaxDt()), maxdtstep_(extrapconfig.MaxDtStep()), minv_(extrapconfig.MinV()), - minvlong_(extrapconfig.MinVLong()), backToTracker_(extrapconfig.BackToTracker()), extrapolateOPA_(extrapconfig.ToOPA()), - toTrackerEnds_(extrapconfig.ToTrackerEnds()), + toTrackerPerimeter_(extrapconfig.ToTrackerPerimeter()), upstream_(extrapconfig.Upstream()), toCRV_(extrapconfig.ToCRV()) {} @@ -87,7 +86,7 @@ namespace mu2e { auto const& ftraj = ktrk.fitTraj(); auto dir0 = ftraj.direction(ftraj.t0()); TimeDir tdir = (dir0.Z() > 0) ? TimeDir::backwards : TimeDir::forwards; - if(toTrackerEnds_)toTrackerPerimeter(ktrk); + if(toTrackerPerimeter_)toTrackerPerimeter(ktrk); if(upstream_){ double starttime = tdir == TimeDir::forwards ? ftraj.range().end() : ftraj.range().begin(); // extrapolate through the IPA in this direction. @@ -426,7 +425,7 @@ namespace mu2e { SurfaceId const sbid(SurfaceIdEnum::CRV_StrongBack); for(auto const& plane : kkg_h->passiveMaterialPlanes()) if(plane.sid_.id() == sbid.id()) strongbacks.push_back(plane); - auto extrapCRV = ExtrapolateCRVRegion(maxdt_,maxdtstep_,btol_,intertol_,minv_,minvlong_,*kkg_h->CRV(),strongbacks,debug_); + auto extrapCRV = ExtrapolateCRVRegion(maxdt_,maxdtstep_,btol_,intertol_,minv_,*kkg_h->CRV(),strongbacks,debug_); if(debug_ > 5){std::cout << "Extrapolating to CRV with " << extrapCRV.sectors().size() << " sectors + " << strongbacks.size() << " strongbacks" << std::endl; for(auto const& sector : kkg_h->CRV()->sectors()) diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index a61bd059ab..60ec284281 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -118,9 +118,8 @@ namespace mu2e { fhicl::Atom MaxDt { Name("MaxDt"), Comment("Maximum time to extrapolate a fit (ns)") }; fhicl::Atom MaxDtStep { Name("MaxDtStep"), Comment("Maximum time step when extrapolating a fit (ns)") }; fhicl::Atom MinV { Name("MinV"), Comment("Minimum velocity perp to the surface extrapolate a straight track fit") }; - fhicl::Atom MinVLong { Name("MinVLong"), Comment("Minimum velocity perp to a z-normal CRV end-cap plane (mm/ns); suppresses grazing-cosmic fakes"), 265.0f }; fhicl::Atom BackToTracker { Name("BackToTracker"), Comment("Extrapolate reflecting tracks back to the tracker") }; - fhicl::Atom ToTrackerEnds { Name("ToTrackerEnds"), Comment("Extrapolate tracks to the tracker ends") }; + fhicl::Atom ToTrackerPerimeter { Name("ToTrackerPerimeter"), Comment("Extrapolate tracks to the tracker perimeter (z-ends and outer wall)") }; fhicl::Atom Upstream { Name("Upstream"), Comment("Extrapolate tracks upstream") }; fhicl::Atom ToOPA { Name("ToOPA"), Comment("Test tracks for intersection with the OPA") }; fhicl::Atom ToCRV { Name("ToCRV"), Comment("Extrapolate tracks to the CRV modules") };