From 4c38bd55239620e831518ede0b7c284382525ceb Mon Sep 17 00:00:00 2001 From: Rob Mina Date: Wed, 1 Jul 2026 17:48:42 -0500 Subject: [PATCH] Truth-seeded track extrapolation to disentangle fit from extrapolation bias/uncertainty. --- Mu2eKinKal/fcl/prolog.fcl | 12 ++ Mu2eKinKal/src/CentralHelixFit_module.cc | 169 +++++++++++++++++++++- Mu2eKinKal/src/KinematicLineFit_module.cc | 135 ++++++++++++++++- 3 files changed, 314 insertions(+), 2 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 11fe39cf88..45330ec436 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -595,9 +595,21 @@ Mu2eKinKal.producers.KKCHSeedFitmuM.ModuleSettings.FitParticle : @local::Particl Mu2eKinKal.producers.KKCHSeedFitmuP.ModuleSettings.FitParticle : @local::Particle.muplus Mu2eKinKal.producers.KKCHmu.ModuleSettings.FitParticle : @local::Particle.muminus # charge floats in the fit, this just determines the mass +# CentralHelixFit can also extrapolate a copy of each fit re-seeded with the true muon state +# at the tracker, written as a parallel KKCHmu:TruthSeed KalSeedCollection, to isolate the +# extrapolation error from the fit error. To enable it, set both of: +# Mu2eKinKal.producers.KKCHmu.ModuleSettings.TruthSeedDiag : true +# Mu2eKinKal.producers.KKCHmu.ModuleSettings.MCTrajectories : "compressDigiMCs" + Mu2eKinKal.producers.KLSeedFit.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKLine.ModuleSettings.FitParticle : @local::Particle.muminus +# KinematicLineFit can also extrapolate a copy of each fit re-seeded with the true muon state +# (from the tracker VD crossing), written as a parallel KKLine:TruthSeed KalSeedCollection, +# to isolate the extrapolation error from the fit error. To enable it, set both of: +# Mu2eKinKal.producers.KKLine.ModuleSettings.TruthSeedDiag : true +# Mu2eKinKal.producers.KKLine.ModuleSettings.VDStepPointMCs : "compressDigiMCs:virtualdetector" + Mu2eKinKal.producers.KKDeSeedFit.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKDeSeedFit.FitDirection : @local::FitDir.downstream Mu2eKinKal.producers.KKDmuSeedFit.ModuleSettings.FitParticle : @local::Particle.muminus diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index ba96bd6135..c2df5db306 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -41,6 +41,10 @@ #include "Offline/RecoDataProducts/inc/CosmicTrackSeed.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +// MC (TruthSeedDiag only) +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/MCDataProducts/inc/MCTrajectory.hh" +#include "Offline/GeometryService/inc/DetectorSystem.hh" // KinKal #include "KinKal/Fit/Track.hh" #include "KinKal/Fit/Config.hh" @@ -72,6 +76,9 @@ #include #include #include +#include +#include +#include using KTRAJ= KinKal::CentralHelix; // this must come before HelixFit using namespace ROOT::Math; @@ -124,6 +131,13 @@ namespace mu2e { fhicl::Atom sampleTBuff { Name("SampleTimeBuffer"), Comment("Time buffer for sample intersections (nsec)") }; fhicl::Atom useFitCharge { Name("UseFitCharge"), Comment("Set the PDG particle according to the fit charge; otherwise reject fits that don't agree with the PDG particle charge") }; fhicl::Atom minCenterRho { Name("MinCenterRho"), Comment("Minimum transverse distance from the helix axis to the Z axis to consider the fit non-degenerate (mm)") }; + fhicl::Atom truthSeedDiag { Name("TruthSeedDiag"), Comment("Diagnostic: also extrapolate a copy of each fit re-seeded with the true muon state, isolating extrapolation error from fit error. Writes KalSeeds under instance 'TruthSeed'"), false }; + fhicl::OptionalAtom mcTrajectoryCollection { Name("MCTrajectories"), Comment("MCTrajectory map providing the true muon state at the fit reference point (required iff TruthSeedDiag)") }; + fhicl::Sequence truthSeedParamConstraints { + Name("TruthSeedParameterConstraints"), + Comment("TruthSeedDiag ParameterHit RMS constraints for CentralHelix parameters d0, phi0, omega, z0, tanDip, t0"), + std::vector{1.0e-3, 1.0e-6, 1.0e-9, 1.0e-3, 1.0e-6, 1.0e-3} + }; }; struct GlobalConfig { @@ -177,6 +191,9 @@ namespace mu2e { SurfaceIdCollection ssids_; KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit bool constraining_; + bool truthSeedDiag_; // diagnostic: extrapolate a truth-seeded copy of each fit + art::ProductToken mctrajcol_T_; // MC trajectories (truth seed) + std::vector truthSeedParamConstraints_; // ParameterHit RMS constraints for truth seeding }; CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings}, @@ -199,13 +216,31 @@ namespace mu2e { useFitCharge_(settings().modSettings().useFitCharge()), minCenterRho_(settings().modSettings().minCenterRho()), sampleinrange_(settings().modSettings().sampleInRange()), - sampleinbounds_(settings().modSettings().sampleInBounds()) + sampleinbounds_(settings().modSettings().sampleInBounds()), + truthSeedDiag_(settings().modSettings().truthSeedDiag()), + truthSeedParamConstraints_(settings().modSettings().truthSeedParamConstraints()) { // collection handling for(const auto& cseedtag : settings().modSettings().seedCollections()) { cseedCols_.emplace_back(consumes(cseedtag)); } produces(); produces(); // produces(); + // TruthSeedDiag: a parallel truth-seeded extrapolation, written under instance "TruthSeed" + if(truthSeedDiag_){ + art::InputTag mctrajtag; + if(!settings().modSettings().mcTrajectoryCollection(mctrajtag)) + throw cet::exception("RECO") << "mu2e::CentralHelixFit: TruthSeedDiag requires MCTrajectories" << endl; + if(truthSeedParamConstraints_.size() != KinKal::NParams()) + throw cet::exception("RECO") << "mu2e::CentralHelixFit: TruthSeedParameterConstraints must have " + << KinKal::NParams() << " entries" << endl; + for(auto const& sigma : truthSeedParamConstraints_){ + if(sigma <= 0.0) + throw cet::exception("RECO") << "mu2e::CentralHelixFit: TruthSeedParameterConstraints entries must be positive" << endl; + } + mctrajcol_T_ = consumes(mctrajtag); + produces("TruthSeed"); + produces("TruthSeed"); + } // build the initial seed covariance auto const& seederrors = settings().modSettings().seederrors(); if(seederrors.size() != KinKal::NParams()) @@ -269,6 +304,12 @@ namespace mu2e { // create output unique_ptr ktrkcol(new KKTRKCOL ); unique_ptr kkseedcol(new KalSeedCollection ); + // TruthSeedDiag parallel outputs + truth inputs + unique_ptr ktrkcol_truth(new KKTRKCOL ); + unique_ptr kkseedcol_truth(new KalSeedCollection ); + MCTrajectoryCollection const* mctrajs = nullptr; + GeomHandle det; + if(truthSeedDiag_) mctrajs = event.getValidHandle(mctrajcol_T_).product(); unsigned nseed(0); for (auto const& cseedtag : cseedCols_) { @@ -357,6 +398,128 @@ namespace mu2e { kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *ktrk ); goodfit = goodFit(*ktrk); } + // [TruthSeedDiag] Build a parallel track from the piecewise tracker truth trajectory, + // constrained to it by ParameterHits, and extrapolate that (via the PUBLIC KinKal + // reconstruction interface -- no fitTrajMutable substitution). The truth-seeded CRV + // residual isolates the extrapolation (B-field/material) error; differencing it from + // the fit-seeded residual gives the fit's contribution (the charge over-rotation). + if(truthSeedDiag_ && goodfit && extrap_ && mctrajs != nullptr){ + // The fit trajectory is PIECEWISE (one CentralHelix per domain) with dE/dx applied + // between pieces, so its |p| steps down through the tracker. A single truth helix has + // CONSTANT |p| and so cannot match the truth at every tracker surface for a low-|p| + // muon that sheds a large momentum fraction crossing the tracker. So we rebuild the + // truth seed with the SAME domain structure as the fit: one truth helix per fit piece, + // each anchored at the true muon state (MC trajectory) at that piece's location. This + // reproduces the truth |p| -- tracker energy loss included -- at every surface, giving a + // ~zero tracker residual at all momenta. The extrapolation then carries it to the CRV. + // NB: do NOT dereference the trajectory's SimParticle Ptr -- cosmic resampling/mixing + // leaves it dangling (ProductNotFound); the trajectory POINTS (pos/time/KE) are valid. + int tcharge = (ktrk->fitTraj().nearestPiece(ktrk->fitTraj().t0()).charge() > 0.0) ? 1 : -1; // from the fit + // truth muon state (pos, mom, time) from the MC trajectory point nearest a given location + auto trueStateAt = [&](VEC3 const& at, VEC3& tpos, VEC3& tmom, double& ttime)->bool{ + double bestd2 = std::numeric_limits::max(); bool found = false; + for(auto const& tjpair : *mctrajs){ + auto const& pts = tjpair.second.points(); + if(pts.size() < 2) continue; + for(size_t i = 0; i < pts.size(); ++i){ + auto pdh = det->toDetector(pts[i].pos()); + VEC3 pd(pdh.x(),pdh.y(),pdh.z()); + double d2 = (pd - at).Mag2(); + if(d2 >= bestd2) continue; + size_t j = (i + 1 < pts.size()) ? i + 1 : i - 1; // local tangent for the direction + auto qdh = det->toDetector(pts[j].pos()); + VEC3 qd(qdh.x(),qdh.y(),qdh.z()); + VEC3 dir = (j > i) ? (qd - pd) : (pd - qd); + double dn = sqrt(dir.Mag2()); + if(dn <= 0.0) continue; + dir /= dn; + double ke = pts[i].kineticEnergy(); + double pmag = std::sqrt(ke*ke + 2.0*ke*mass_); // |p| from kinetic energy + bestd2 = d2; found = true; tpos = pd; tmom = pmag*dir; ttime = pts[i].t(); + } + } + return found; + }; + // build the piecewise truth trajectory (one truth helix per fit domain) and remember each + // piece's reference time, to constrain it with a ParameterHit below. + PTRAJ truthpt; + std::vector truthHitTimes; + for(auto const& fpieceptr : ktrk->fitTraj().pieces()){ + auto const& fpiece = *fpieceptr; + auto const htime = fpiece.range().mid(); + VEC3 tpos, tmom; double ttime = 0.0; + if(!trueStateAt(fpiece.position3(htime), tpos, tmom, ttime)) continue; + auto tbnom = kkbf_->fieldVect(tpos); + KTRAJ tpiece(KinKal::VEC4(tpos.X(),tpos.Y(),tpos.Z(),ttime), + KinKal::MOM4(tmom.X(),tmom.Y(),tmom.Z(),mass_), + tcharge, tbnom.Z(), fpiece.range()); + tpiece.params() = KinKal::Parameters(tpiece.params().parameters(), seedcov_); + truthpt.append(tpiece); + truthHitTimes.emplace_back(htime); + } + if(!truthHitTimes.empty()){ + // Constrain each truth piece with a ParameterHit at its reference time (RMS from + // TruthSeedParameterConstraints). This pins the fit to the truth via the PUBLIC KinKal + // reconstruction interface -- no fitTrajMutable substitution -- at the cost of being a + // (tightly) constrained fit rather than an exact trajectory swap. + PARAMHITCOL truthParamHits; + PARAMHIT::PMASK truthMask; + truthMask.fill(true); + auto addTruthHit = [&](double hitTime, double covarianceScale) { + auto cparams = truthpt.nearestPiece(hitTime).params(); + for(size_t ipar = 0; ipar < KinKal::NParams(); ++ipar){ + for(size_t jpar = 0; jpar < KinKal::NParams(); ++jpar){ + cparams.covariance()[ipar][jpar] = 0.0; + } + auto const sigma = truthSeedParamConstraints_.at(ipar); + cparams.covariance()[ipar][ipar] = covarianceScale*sigma*sigma; + } + truthParamHits.push_back(std::make_shared( + hitTime, truthpt, cparams, truthMask)); + }; + // a single piece gives a zero-length fit range (lowNDOF); split it into two hits a + // sliver apart, with looser (x2) constraint, as the line tail does. + if(truthHitTimes.size() == 1){ + constexpr double truthHitDt = 1.0e-3; + addTruthHit(truthHitTimes.front() - 0.5*truthHitDt, 2.0); + addTruthHit(truthHitTimes.front() + 0.5*truthHitDt, 2.0); + } else { + for(auto const hitTime : truthHitTimes) addTruthHit(hitTime, 1.0); + } + // copy the fit's field domains so the truth fit/extrapolation see the same BField model + KKTRK::DOMAINCOL truthDomains; + for(auto const& domain : ktrk->domains()){ + truthDomains.emplace(std::make_shared(*domain)); + } + KKTRK::PKTRAJPTR truthTraj = std::make_unique(truthpt); + KKSTRAWHITCOL nostrawhits; KKSTRAWXINGCOL nostrawxings; KKCALOHITCOL nocalohits; + // one iteration, no convergence churn: the ParameterHits already pin the params + auto truthConfig = config_; + truthConfig.maxniter_ = 1; + truthConfig.divdchisq_ = std::numeric_limits::max(); + truthConfig.pdchisq_ = std::numeric_limits::max(); + truthConfig.divgap_ = std::numeric_limits::max(); + auto ktrk_truth = make_unique(truthConfig,*kkbf_,fitpart,truthTraj, + nostrawhits,nostrawxings,nocalohits,truthParamHits,truthDomains); + if(ktrk_truth->fitStatus().usable()){ + extrap_->extrapolate(*ktrk_truth); + TrkFitFlag tflag; tflag.merge(fitflag_); tflag.merge(TrkFitFlag::FitOK); + auto truthKSeed = kkfit_.createSeed(*ktrk_truth,tflag,*calo_h,*nominalTracker_h); + // EventNtuple and other consumers use the detector-hit list for bookkeeping even + // when hit details are disabled; carry the reco fit's hit metadata onto the truth + // KalSeed without adding the detector measurements to the truth-constrained fit. + auto recoMetadata = kkfit_.createSeed(*ktrk,tflag,*calo_h,*nominalTracker_h); + truthKSeed._hits = std::move(recoMetadata._hits); + truthKSeed._hitcalibs = std::move(recoMetadata._hitcalibs); + truthKSeed._chit = std::move(recoMetadata._chit); + kkseedcol_truth->push_back(std::move(truthKSeed)); + ktrkcol_truth->push_back(ktrk_truth.release()); + } else if(print_ > 0) { + std::cout << "CentralHelixFit TruthSeed ParameterHit fit unusable: " + << ktrk_truth->fitStatus() << std::endl; + } + } + } // extrapolate as required if(goodfit && extrap_)extrap_->extrapolate(*ktrk); if(print_>1)ktrk->printFit(std::cout,print_); @@ -385,6 +548,10 @@ namespace mu2e { if(print_ > 0) std::cout << "Fitted " << ktrkcol->size() << " tracks from " << nseed << " Seeds" << std::endl; event.put(move(ktrkcol)); event.put(move(kkseedcol)); + if(truthSeedDiag_){ + event.put(move(ktrkcol_truth),"TruthSeed"); + event.put(move(kkseedcol_truth),"TruthSeed"); + } } bool CentralHelixFit::goodFit(KKTRK const& ktrk) const { diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index 293240b649..42b5975dc6 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -37,6 +37,10 @@ #include "Offline/RecoDataProducts/inc/KKLine.hh" #include "Offline/DataProducts/inc/SurfaceId.hh" #include "Offline/KinKalGeom/inc/KinKalGeom.hh" +// MC (TruthSeedDiag only) +#include "Offline/MCDataProducts/inc/StepPointMC.hh" +#include "Offline/MCDataProducts/inc/SimParticle.hh" +#include "Offline/GeometryService/inc/DetectorSystem.hh" // KinKal #include "KinKal/Fit/Track.hh" #include "KinKal/Fit/Config.hh" @@ -63,6 +67,7 @@ #include #include #include +#include using namespace std; //using namespace KinKal; namespace mu2e { @@ -120,6 +125,13 @@ namespace mu2e { fhicl::Atom sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") }; fhicl::Atom interTol { Name("IntersectionTolerance"), Comment("Tolerance for surface intersections (mm)") }; fhicl::Atom sampleTBuff { Name("SampleTimeBuffer"), Comment("Time buffer for sample intersections (nsec)") }; + fhicl::Atom truthSeedDiag { Name("TruthSeedDiag"), Comment("Diagnostic: also extrapolate a copy of each fit re-seeded with the true muon state, isolating extrapolation error from fit error. Writes KalSeeds under instance 'TruthSeed'"), false }; + fhicl::OptionalAtom vdStepCollection { Name("VDStepPointMCs"), Comment("Virtual-detector StepPointMCs providing the true muon state at the tracker (required iff TruthSeedDiag)") }; + fhicl::Sequence truthSeedParamConstraints { + Name("TruthSeedParameterConstraints"), + Comment("TruthSeedDiag ParameterHit RMS constraints for KinematicLine parameters d0, phi0, z0, cost, t0, mom"), + std::vector{1.0e-3, 1.0e-6, 1.0e-3, 1.0e-6, 1.0e-3, 1.0e-3} + }; }; struct GlobalConfig { @@ -170,6 +182,9 @@ namespace mu2e { KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit Config config_; // initial fit configuration object Config exconfig_; // extension configuration object + bool truthSeedDiag_; // diagnostic: extrapolate a truth-seeded copy of each fit + art::ProductToken vdcol_T_; // VD StepPointMCs (truth seed) + std::vector truthSeedParamConstraints_; // ParameterHit RMS constraints for truth seeding }; KinematicLineFit::KinematicLineFit(const Parameters& settings) : art::EDProducer{settings}, @@ -187,7 +202,9 @@ namespace mu2e { sampleinrange_(settings().modSettings().sampleInRange()), sampleinbounds_(settings().modSettings().sampleInBounds()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), - exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())) + exconfig_(Mu2eKinKal::makeConfig(settings().extSettings())), + truthSeedDiag_(settings().modSettings().truthSeedDiag()), + truthSeedParamConstraints_(settings().modSettings().truthSeedParamConstraints()) { std::string fdir; if(settings().fitDirection(fdir))fdir_ = fdir; @@ -196,6 +213,22 @@ namespace mu2e { produces(); produces(); produces(); + // TruthSeedDiag: a parallel truth-seeded extrapolation, written under instance "TruthSeed" + if(truthSeedDiag_){ + art::InputTag vdtag; + if(!settings().modSettings().vdStepCollection(vdtag)) + throw cet::exception("RECO") << "mu2e::KinematicLineFit: TruthSeedDiag requires VDStepPointMCs" << endl; + if(truthSeedParamConstraints_.size() != KinKal::NParams()) + throw cet::exception("RECO") << "mu2e::KinematicLineFit: TruthSeedParameterConstraints must have " + << KinKal::NParams() << " entries" << endl; + for(auto const& sigma : truthSeedParamConstraints_){ + if(sigma <= 0.0) + throw cet::exception("RECO") << "mu2e::KinematicLineFit: TruthSeedParameterConstraints entries must be positive" << endl; + } + vdcol_T_ = consumes(vdtag); + produces("TruthSeed"); + produces("TruthSeed"); + } // build the initial seed covariance auto const& seederrors = settings().modSettings().seederrors(); if(seederrors.size() != KinKal::NParams()) @@ -251,6 +284,12 @@ namespace mu2e { unique_ptr kktrkcol(new KKTRKCOL ); unique_ptr kkseedcol(new KalSeedCollection ); //Needs to return a KalSeed unique_ptr kkseedassns(new KalLineAssns()); + // TruthSeedDiag parallel outputs + truth inputs + unique_ptr kktrkcol_truth(new KKTRKCOL ); + unique_ptr kkseedcol_truth(new KalSeedCollection ); + StepPointMCCollection const* vdsteps = nullptr; + GeomHandle det; + if(truthSeedDiag_) vdsteps = event.getValidHandle(vdcol_T_).product(); auto KalSeedCollectionPID = event.getProductID(); auto KalSeedCollectionGetter = event.productGetter(KalSeedCollectionPID); // find the track seed collections @@ -310,6 +349,96 @@ namespace mu2e { kkfit_.extendTrack(exconfig_,*kkbf_, *tracker,*strawresponse, chcol, *calo_h, cc_H, *kktrk ); } goodfit = goodFit(*kktrk); + // [TruthSeedDiag] Build a parallel track seeded with the TRUE muon state at the tracker + // (nearest VD crossing), constrained to it by ParameterHits, and extrapolate that -- via the + // PUBLIC KinKal reconstruction interface (no fitTrajMutable). The truth-seeded CRV residual + // isolates the extrapolation (material/scattering) error from the fit's contribution. + if(truthSeedDiag_ && goodfit && extrap_ && vdsteps != nullptr){ + double tt0 = kktrk->fitTraj().t0(); + auto refpos = kktrk->fitTraj().position3(tt0); + double bestd2 = std::numeric_limits::max(); + VEC3 tpos, tmom; double ttime = 0.0; int tcharge = 0; bool found = false; + for(auto const& vds : *vdsteps){ // nearest true muon VD crossing to the fit t0 point + auto const& sp = vds.simParticle(); + if(sp.isNull() || std::abs(static_cast(sp->pdgId())) != PDGCode::mu_minus) continue; + auto const hd = det->toDetector(vds.position()); + VEC3 dpos(hd.x(),hd.y(),hd.z()); + double d2 = (dpos-refpos).Mag2(); + if(d2 < bestd2){ + bestd2 = d2; found = true; ttime = vds.time(); + tpos = dpos; tmom = VEC3(vds.momentum().x(),vds.momentum().y(),vds.momentum().z()); + tcharge = (static_cast(sp->pdgId()) > 0) ? -1 : 1; // mu-(+13)->-1, mu+(-13)->+1 + } + } + if(found){ + try { + // field-off straight line: same non-zero nominal bfield (VEC3) as makeSeedTraj + VEC3 tbnom(0.0,0.0,0.001); + KTRAJ truthtraj(KinKal::VEC4(tpos.X(),tpos.Y(),tpos.Z(),ttime), + KinKal::MOM4(tmom.X(),tmom.Y(),tmom.Z(),mass_), + tcharge, tbnom, kktrk->fitTraj().range()); + truthtraj.params() = KinKal::Parameters(truthtraj.params().parameters(), seedcov_); + PTRAJ truthpt(truthtraj); + // constrain the line to the truth params via ParameterHits (public reco interface). A + // single piece is a zero-length fit range (lowNDOF), so pin TWO hits a sliver apart. + PARAMHITCOL truthParamHits; + PARAMHIT::PMASK truthMask; truthMask.fill(true); + auto addTruthHit = [&](double hitTime, double covarianceScale){ + auto cparams = truthpt.nearestPiece(hitTime).params(); + for(size_t ipar = 0; ipar < KinKal::NParams(); ++ipar){ + for(size_t jpar = 0; jpar < KinKal::NParams(); ++jpar) cparams.covariance()[ipar][jpar] = 0.0; + auto const sigma = truthSeedParamConstraints_.at(ipar); + cparams.covariance()[ipar][ipar] = covarianceScale*sigma*sigma; + } + truthParamHits.push_back(std::make_shared(hitTime, truthpt, cparams, truthMask)); + }; + constexpr double truthHitDt = 1.0e-3; + // Pin the constraints at an IN-RANGE time. The KKLine fit is field-off, so + // kktrk->domains() is empty and Track::fit's detrange collapses to just these hit + // times (Track.hh:202-214 only extend detrange to surviving domains). The VD + // crossing time `ttime` can fall outside the fit trajectory's time range; pinning + // the hits there leaves the single truth piece (range = fitTraj range) not + // overlapping detrange, so PiecewiseTrajectory::setRange(trim) erases every piece + // and then dereferences the empty deque -> SIGSEGV. The single-piece KinematicLine + // params are time-independent, so clamping the constraint time into the fit range + // pins exactly the same truth params. (CentralHelix avoids this via field domains + + // in-range piece-midpoint sampling.) + auto const& frange = kktrk->fitTraj().range(); + double tctr = std::max(frange.begin() + truthHitDt, std::min(ttime, frange.end() - truthHitDt)); + addTruthHit(tctr - 0.5*truthHitDt, 2.0); + addTruthHit(tctr + 0.5*truthHitDt, 2.0); + KKTRK::DOMAINCOL truthDomains; + for(auto const& domain : kktrk->domains()) truthDomains.emplace(std::make_shared(*domain)); + KKTRK::PKTRAJPTR truthTraj = std::make_unique(truthpt); + KKSTRAWHITCOL nostrawhits; KKSTRAWXINGCOL nostrawxings; KKCALOHITCOL nocalohits; + // one iteration, no convergence churn: the ParameterHits already pin the params + auto truthConfig = config_; + truthConfig.maxniter_ = 1; + truthConfig.divdchisq_ = std::numeric_limits::max(); + truthConfig.pdchisq_ = std::numeric_limits::max(); + truthConfig.divgap_ = std::numeric_limits::max(); + auto kktrk_truth = make_unique(truthConfig,*kkbf_,fpart_,truthTraj, + nostrawhits,nostrawxings,nocalohits,truthParamHits,truthDomains); + if(kktrk_truth->fitStatus().usable()){ + extrap_->extrapolate(*kktrk_truth); + TrkFitFlag tflag(hptr->status()); tflag.merge(TrkFitFlag::KKLine); + auto truthKSeed = kkfit_.createSeed(*kktrk_truth,tflag,*calo_h,*nominalTracker_h); + // carry the reco fit's hit metadata onto the truth KalSeed (the truth-constrained + // fit has no detector hits, but EventNtuple consumers expect the hit list) + auto recoMetadata = kkfit_.createSeed(*kktrk,tflag,*calo_h,*nominalTracker_h); + truthKSeed._hits = std::move(recoMetadata._hits); + truthKSeed._hitcalibs = std::move(recoMetadata._hitcalibs); + truthKSeed._chit = std::move(recoMetadata._chit); + kkseedcol_truth->push_back(std::move(truthKSeed)); + kktrkcol_truth->push_back(kktrk_truth.release()); + } else if(print_ > 0){ + std::cout << "KinematicLineFit TruthSeed ParameterHit fit unusable: " << kktrk_truth->fitStatus() << std::endl; + } + } catch (std::exception const& error) { + if(print_ > 0) std::cout << "KinematicLineFit TruthSeed Error " << error.what() << std::endl; + } + } + } // extrapolate as required if(goodfit && extrap_) extrap_->extrapolate(*kktrk); bool save = goodFit(*kktrk); @@ -334,6 +463,10 @@ namespace mu2e { event.put(move(kktrkcol)); event.put(move(kkseedcol)); event.put(move(kkseedassns)); + if(truthSeedDiag_){ + event.put(move(kktrkcol_truth),"TruthSeed"); + event.put(move(kkseedcol_truth),"TruthSeed"); + } } KTRAJ KinematicLineFit::makeSeedTraj(CosmicTrackSeed const& hseed) const {