diff --git a/BFieldGeom/inc/BFieldConfig.hh b/BFieldGeom/inc/BFieldConfig.hh index b7100f5345..f01f9bd88b 100644 --- a/BFieldGeom/inc/BFieldConfig.hh +++ b/BFieldGeom/inc/BFieldConfig.hh @@ -64,9 +64,12 @@ namespace mu2e { bool flipBFieldMaps() const { return flipBFieldMaps_; } + // smallest |B| (T) KinKal may use for field-corrected transport; 0 disables its low-field protection + double minUsableField() const { return minUsableField_; } + private: BFieldConfig() - : scaleFactor_(1.), writeBinaries_(false), verbosityLevel_(1), flipBFieldMaps_(false) {} + : scaleFactor_(1.), writeBinaries_(false), verbosityLevel_(1), flipBFieldMaps_(false), minUsableField_(0.) {} // G4BL, PARAM or possible future types. BFMapType mapType_; @@ -91,6 +94,7 @@ namespace mu2e { bool writeBinaries_; int verbosityLevel_; bool flipBFieldMaps_; + double minUsableField_; }; } // namespace mu2e diff --git a/GeometryService/src/BFieldConfigMaker.cc b/GeometryService/src/BFieldConfigMaker.cc index 6149ead151..3b78fd0668 100644 --- a/GeometryService/src/BFieldConfigMaker.cc +++ b/GeometryService/src/BFieldConfigMaker.cc @@ -47,6 +47,7 @@ namespace mu2e { bfconf_->flipBFieldMaps_ = config.getBool("bfield.flipMaps", false); bfconf_->scaleFactor_ = config.getDouble("bfield.scaleFactor", 1.0); + bfconf_->minUsableField_ = config.getDouble("bfield.minUsableField", 0.0); bfconf_->dsFieldForm_ = BFieldConfig::DSFieldModel( config.getInt("detSolFieldForm", BFieldConfig::dsModelUniform)); diff --git a/Mu2eG4/geom/bfgeom_v01_minfield.txt b/Mu2eG4/geom/bfgeom_v01_minfield.txt new file mode 100644 index 0000000000..2e255171e4 --- /dev/null +++ b/Mu2eG4/geom/bfgeom_v01_minfield.txt @@ -0,0 +1,17 @@ +// +// bfgeom_v01 with KinKal low-field protection enabled at the null-field threshold. +// +// bfield.minUsableField is the smallest |B| (T) KinKal uses for field-corrected transport. Below it, and +// outside the maps, bfcorr extrapolation hands off to geometric free-particle continuation instead of +// building domains on a field it can't parameterize. 1e-6 T is KinKal's null-field value +// (BFieldMap::zeroField), so this only removes the degenerate samples. Reconstruction only: Geant4 ignores it. + +#include "Offline/Mu2eG4/geom/bfgeom_v01.txt" + +double bfield.minUsableField = 1.0e-6; + +// +// This tells emacs to view this file in c++ mode. +// Local Variables: +// mode:c++ +// End: diff --git a/Mu2eG4/geom/bfgeom_v01_minfield0p3.txt b/Mu2eG4/geom/bfgeom_v01_minfield0p3.txt new file mode 100644 index 0000000000..f59af60e7e --- /dev/null +++ b/Mu2eG4/geom/bfgeom_v01_minfield0p3.txt @@ -0,0 +1,16 @@ +// +// bfgeom_v01 with KinKal low-field protection at 0.3 T. +// +// bfield.minUsableField is the smallest |B| (T) KinKal uses for field-corrected transport; see +// bfgeom_v01_minfield.txt. 0.3 T hands off to geometric continuation while real bending remains, so it +// trades extrapolation accuracy for stability. Reconstruction only: Geant4 ignores it. + +#include "Offline/Mu2eG4/geom/bfgeom_v01.txt" + +double bfield.minUsableField = 0.3; + +// +// This tells emacs to view this file in c++ mode. +// Local Variables: +// mode:c++ +// End: diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index d1440a4d71..29e43eac39 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -53,6 +53,7 @@ Mu2eKinKal : { DivergenceGap : 20 # mm BFieldCorrection : false BCorrTolerance : 1e-2 # momemntum fraction + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded ProcessEnds : false MetaIterationSettings : [ # annealing temp, strawhit updater algorithm @@ -89,6 +90,7 @@ Mu2eKinKal : { BFieldCorrection : true ProcessEnds : false BCorrTolerance : 1e-2 # relative momemntum precision + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [ # annealing temp, strawhit updater algorithm [ 2.0, "CADSHU" ], @@ -128,6 +130,7 @@ Mu2eKinKal : { BFieldCorrection : false ProcessEnds : false BCorrTolerance : 1e-4 # momemntum fraction + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [ # annealing temp, strawhit updater algorithms [ 2.0, "CADSHU" ], @@ -194,6 +197,7 @@ Mu2eKinKal : { DivergenceGap : 10 # mm BFieldCorrection : true BCorrTolerance : 1e-2 # momemntum fraction + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded ProcessEnds : false MetaIterationSettings : [ # annealing temp, strawhit updater algorithm @@ -227,6 +231,7 @@ Mu2eKinKal : { BFieldCorrection : false ProcessEnds : false BCorrTolerance : 1e-4 # momemntum fraction + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [] CADSHUSettings : [] StrawXingUpdaterSettings : [] @@ -244,6 +249,7 @@ Mu2eKinKal : { BFieldCorrection : true ProcessEnds : false BCorrTolerance : 1e-2 # relative momemntum precision + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [ # annealing temp, strawhit updater algorithm [ 2.0, "CADSHU" ], @@ -283,6 +289,7 @@ Mu2eKinKal : { BFieldCorrection : true ProcessEnds : true BCorrTolerance : 1e-4 # momemntum fraction + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [ # annealing temp, strawhit updater algorithms [ 2.0, "CADSHU" ], @@ -358,6 +365,7 @@ Mu2eKinKal : { DivergenceGap : 100 # mm ProcessEnds : true BCorrTolerance : 1e-4 + MinDomainTimeStep : 0.0 # (ns); 0 leaves the BField domain step unbounded MetaIterationSettings : [ [ 0.0, "" ] ] CADSHUSettings : [ ] BkgANNSHUSettings : [ ] diff --git a/Mu2eKinKal/inc/KKBField.hh b/Mu2eKinKal/inc/KKBField.hh index f916223501..78b2383d66 100644 --- a/Mu2eKinKal/inc/KKBField.hh +++ b/Mu2eKinKal/inc/KKBField.hh @@ -17,7 +17,8 @@ namespace mu2e using Grad = ROOT::Math::SMatrix; // field gradient: ie dBi/d(x,y,z) // construct from BField object and system translator. // This should be a single BField map valid in the detector system, to avoid making continuous translations FIXME! - KKBField(BFieldManager const& bfmgr, DetectorSystem const& det) : bfmgr_(bfmgr), det_(det) {} + // The minimum usable field comes from the BField geometry (bfield.minUsableField), so every module in a job shares it. + KKBField(BFieldManager const& bfmgr, DetectorSystem const& det); virtual ~KKBField() {} // KinKal BField interface // return value of the field at a poin diff --git a/Mu2eKinKal/inc/KKExtrap.hh b/Mu2eKinKal/inc/KKExtrap.hh index a72489b0d6..13a1b179aa 100644 --- a/Mu2eKinKal/inc/KKExtrap.hh +++ b/Mu2eKinKal/inc/KKExtrap.hh @@ -19,6 +19,7 @@ #include "Offline/Mu2eKinKal/inc/KKShellXing.hh" #include "Offline/KinKalGeom/inc/KKMaterial.hh" #include "KinKal/Geometry/ParticleTrajectoryIntersect.hh" +#include "cetlib_except/exception.h" #include "Offline/GeometryService/inc/GeomHandle.hh" #include #include @@ -41,6 +42,8 @@ namespace mu2e { class KKExtrap { public: explicit KKExtrap(KKExtrapConfig const& exconfig); + // throw at configuration if the fit config's minimum domain step can't fit this extrapolation's maximum step + void checkConfig(KinKal::Config const& config) const; // extrapolation functions; these are templated on the type of trajectory template void extrapolate(KKTrack& ktrk) const; template bool extrapolateIPA(KKTrack& ktrk,TimeDir trkdir) const; @@ -80,6 +83,12 @@ namespace mu2e { toCRV_(extrapconfig.ToCRV()) {} + inline void KKExtrap::checkConfig(KinKal::Config const& config) const { + if(config.bfcorr_ && config.mindtstep_ > maxdtstep_) + throw cet::exception("RECO") << "mu2e::KKExtrap: MinDomainTimeStep " << config.mindtstep_ + << " ns exceeds ExtrapolationSettings.MaxDtStep " << maxdtstep_ << " ns" << std::endl; + } + template void KKExtrap::extrapolate(KKTrack& ktrk) const { GeomHandle kkg_h; // define the time direction according to the fit direction inside the tracker diff --git a/Mu2eKinKal/inc/KKFitSettings.hh b/Mu2eKinKal/inc/KKFitSettings.hh index 60ec284281..4395e5999d 100644 --- a/Mu2eKinKal/inc/KKFitSettings.hh +++ b/Mu2eKinKal/inc/KKFitSettings.hh @@ -4,6 +4,7 @@ // Structs for configuring the Mu2e KinKal fit // #include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/OptionalAtom.h" #include "fhiclcpp/types/Sequence.h" #include "fhiclcpp/types/OptionalSequence.h" #include "fhiclcpp/types/Tuple.h" @@ -35,6 +36,8 @@ namespace mu2e { fhicl::Atom bfieldCorr { Name("BFieldCorrection"), Comment("Apply correction for BField inhomogeneity") }; fhicl::Atom ends { Name("ProcessEnds"), Comment("Process purely passive sites at the time range ends") }; fhicl::Atom btol { Name("BCorrTolerance"), Comment("Tolerance on BField correction momentum fractional accuracy (dimensionless)") }; + fhicl::Atom minDomainTimeStep { Name("MinDomainTimeStep"), Comment("Minimum BField domain time step (ns); 0 leaves it unbounded"), 0.0 }; + fhicl::OptionalAtom domainMargin { Name("DomainMargin"), Comment("Maximum time a BField domain may extend beyond the active range (ns, >= 0); absent leaves domains unconfined") }; // Updater settings using MetaIterationSettings = fhicl::Sequence>; MetaIterationSettings miConfig { Name("MetaIterationSettings"), Comment("Temperature (dimensionless), StrawHitUpdater algorithm") }; diff --git a/Mu2eKinKal/src/CentralHelixFit_module.cc b/Mu2eKinKal/src/CentralHelixFit_module.cc index dd0aaca7de..11e4fb51a2 100644 --- a/Mu2eKinKal/src/CentralHelixFit_module.cc +++ b/Mu2eKinKal/src/CentralHelixFit_module.cc @@ -233,6 +233,8 @@ namespace mu2e { } // setup extrapolation if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); + // extrapolation runs under the extension config when there is one + if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_); // surfaces to sample; this interface is deprecatecd and should be replaced with extrapolation TODO for(auto const& sidname : settings().modSettings().sampleSurfaces()) { @@ -361,9 +363,13 @@ namespace mu2e { goodfit = goodFit(*ktrk); } // extrapolate as required - if(goodfit && extrap_)extrap_->extrapolate(*ktrk); + if(goodfit && extrap_){ + extrap_->extrapolate(*ktrk); + goodfit = goodFit(*ktrk); // extrapolation can mark the fit failing + } if(print_>1)ktrk->printFit(std::cout,print_); - if(goodfit || saveall_){ + // a fit that failed before its trajectory was built has nothing to save + if((goodfit || saveall_) && ktrk->hasTraj()){ TrkFitFlag fitflag; fitflag.merge(fitflag_); if(goodfit) diff --git a/Mu2eKinKal/src/KKBField.cc b/Mu2eKinKal/src/KKBField.cc index fee672de22..214f1acaa5 100644 --- a/Mu2eKinKal/src/KKBField.cc +++ b/Mu2eKinKal/src/KKBField.cc @@ -1,9 +1,14 @@ #include "Offline/Mu2eKinKal/inc/KKBField.hh" +#include "Offline/BFieldGeom/inc/BFieldConfig.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" #include "cetlib_except/exception.h" namespace mu2e { using Grad = ROOT::Math::SMatrix; using SVEC3 = KinKal::SVEC3; + KKBField::KKBField(BFieldManager const& bfmgr, DetectorSystem const& det) : + KinKal::BFieldMap(GeomHandle()->minUsableField()), bfmgr_(bfmgr), det_(det) {} + VEC3 KKBField::fieldVect(VEC3 const& position) const { // change coordinates to mu2e; the map should be native in detector coordinates FIXME! CLHEP::Hep3Vector vpoint(position.x(),position.y(),position.z()); diff --git a/Mu2eKinKal/src/KKFitSettings.cc b/Mu2eKinKal/src/KKFitSettings.cc index 9f25a90e75..a8d78bce14 100644 --- a/Mu2eKinKal/src/KKFitSettings.cc +++ b/Mu2eKinKal/src/KKFitSettings.cc @@ -2,6 +2,7 @@ #include "Offline/Mu2eKinKal/inc/WireHitState.hh" #include "Offline/Mu2eKinKal/inc/StrawHitUpdaters.hh" #include "Offline/GeneralUtilities/inc/splitLine.hh" +#include "cetlib_except/exception.h" #include namespace mu2e { @@ -22,6 +23,12 @@ namespace mu2e { config.bfcorr_ = fitconfig.bfieldCorr(); config.ends_ = fitconfig.ends(); config.tol_ = fitconfig.btol(); + config.mindtstep_ = fitconfig.minDomainTimeStep(); + double domainmargin(0.0); + if(fitconfig.domainMargin(domainmargin)){ // absent keeps KinKal's unconfined default + if(domainmargin < 0.0) throw cet::exception("RECO") << "mu2e::Mu2eKinKal::makeConfig: DomainMargin must be >= 0, got " << domainmargin << std::endl; + config.domainmargin_ = domainmargin; + } // create the updaters requested std::vector cadshusettings; std::vector driftannshusettings; diff --git a/Mu2eKinKal/src/KKTruthSeed_module.cc b/Mu2eKinKal/src/KKTruthSeed_module.cc index af9d62ea26..0e80cfb7ba 100644 --- a/Mu2eKinKal/src/KKTruthSeed_module.cc +++ b/Mu2eKinKal/src/KKTruthSeed_module.cc @@ -145,6 +145,8 @@ namespace mu2e { produces(); produces(); if(settings().extrapSettings()) extrap_ = std::make_unique(*settings().extrapSettings()); + // the truth track extrapolates with the extension config's bfcorr_ and mindtstep_ + if(extrap_) extrap_->checkConfig(exconfig_); if(truthSeedParamConstraints_.size() != KinKal::NParams()) throw cet::exception("RECO") << "mu2e::KKTruthSeed: TruthSeedParameterConstraints must have " << KinKal::NParams() << " entries" << std::endl; @@ -219,6 +221,14 @@ namespace mu2e { auto const htime = fpiece.range().mid(); VEC3 tpos, tmom; double ttime = 0.0; if(!trueStateAt(mctrajs, det, fpiece.position3(htime), mass, tpos, tmom, ttime)) continue; + // a helix needs a usable local field: outside the map, or where the helix can't represent the true state there + // (null field, circle centered on the axis), the constructor would refuse it + if constexpr(!isLine){ + VEC4 tpos4(tpos.X(),tpos.Y(),tpos.Z(),ttime); + MOM4 tmom4(tmom.X(),tmom.Y(),tmom.Z(),mass); + if(!kkbf_->usable(tpos) || + !KTRAJ::constructible(KinKal::ParticleState(tpos4,tmom4,tcharge),VEC3(0.0,0.0,kkbf_->fieldVect(tpos).Z()))) continue; + } // the KTRAJ nominal BField: the local field for a helix; a fixed non-zero z for a field-free line // (KinematicLine requires a non-zero bnom for interface consistency). KTRAJ tpiece = [&]{ @@ -289,6 +299,12 @@ namespace mu2e { // without re-running the full extension schedule (which diverges a hit-less truth-pinned fit). truthConfig.bfcorr_ = exconfig_.bfcorr_; truthConfig.tol_ = exconfig_.tol_; + truthConfig.mindtstep_ = exconfig_.mindtstep_; // floors the extrapolation step too + // a field-corrected truth fit needs the reco domains, as in RegrowLoopHelix + if(truthConfig.bfcorr_ && truthDomains.empty()){ + if(debug_ > 0) std::cout << "KKTruthSeed: no BField domains from " << recoseed.domainBounds().size() << " saved bounds; skipping" << std::endl; + return std::nullopt; + } auto ktrk_truth = std::make_unique(truthConfig,*kkbf_,recoseed.particle(),truthTraj, nostrawhits,nostrawxings,nocalohits,truthParamHits,truthDomains); if(!ktrk_truth->fitStatus().usable()){ diff --git a/Mu2eKinKal/src/KinematicLineFit_module.cc b/Mu2eKinKal/src/KinematicLineFit_module.cc index cbe9d71810..1ecc83e40c 100644 --- a/Mu2eKinKal/src/KinematicLineFit_module.cc +++ b/Mu2eKinKal/src/KinematicLineFit_module.cc @@ -219,6 +219,8 @@ namespace mu2e { } // setup extrapolation if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); + // extrapolation runs under the extension config when there is one + if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_); if(print_ > 0) std::cout << config_; } @@ -316,7 +318,8 @@ namespace mu2e { // extrapolate as required if(goodfit && extrap_) extrap_->extrapolate(*kktrk); bool save = goodFit(*kktrk); - if(save || saveall_){ + // a fit that failed before its trajectory was built has nothing to save + if((save || saveall_) && kktrk->hasTraj()){ TrkFitFlag fitflag(hptr->status()); fitflag.merge(TrkFitFlag::KKLine); sampleFit(*kktrk); diff --git a/Mu2eKinKal/src/LoopHelixFit_module.cc b/Mu2eKinKal/src/LoopHelixFit_module.cc index 6dbe153cbf..94249e6cb3 100644 --- a/Mu2eKinKal/src/LoopHelixFit_module.cc +++ b/Mu2eKinKal/src/LoopHelixFit_module.cc @@ -225,6 +225,8 @@ namespace mu2e { } // setup extrapolation if(settings().extrapSettings())extrap_ = make_unique(*settings().extrapSettings()); + // extrapolation runs under the extension config when there is one + if(extrap_) extrap_->checkConfig(exconfig_.schedule().size() > 0 ? exconfig_ : config_); // setup optional fit finalization; this just updates the internals, not the fit result itself if(settings().finalSettings()){ @@ -366,7 +368,8 @@ namespace mu2e { } if(print_>0) printf("[LoopHelixFit::%s] After extending the fit : goodFit = %o, fitcon = %.4f, nHits = %2lu, %lu calo-hits\n", __func__, goodfit, ktrk->fitStatus().chisq_.probability(), ktrk->strawHits().size(), ktrk->caloHits().size()); - if((!goodfit) && (! saveall_)) ktrk.reset(); + // a fit that failed before its trajectory was built has nothing to save + if(!ktrk->hasTraj() || ((!goodfit) && (! saveall_))) ktrk.reset(); return ktrk; } diff --git a/RecoDataProducts/inc/TrkStrawHitSeed.hh b/RecoDataProducts/inc/TrkStrawHitSeed.hh index dc56783c9c..1a1721506b 100644 --- a/RecoDataProducts/inc/TrkStrawHitSeed.hh +++ b/RecoDataProducts/inc/TrkStrawHitSeed.hh @@ -172,7 +172,7 @@ namespace mu2e { float _stoca =0; // reference sensor time of closest approach (TOCA) float _rdoca, _rdocavar =0; // reference (biased) DOCA from the track to the wire, signed by the angular momentum WRT the wire and the measurement end (and variance) float _rdt, _rtocavar =0; // reference (biased) time difference (and variance) at POCA - float _udoca, _udocavar =0; // unbiaed DOCA (and variance) + float _udoca, _udocavar =0; // unbiased DOCA (and variance) WRT the final fit trajectory; the hit-state updaters (ANNs) used the unbiased DOCA from the previous fit cycle. float _udt, _utocavar =0; //unbiased dt and variance float _rupos =0; // reference POCA position along the straw WRT the straw middle float _uupos =0; // unbiased POCA position along the straw WRT the straw middle diff --git a/scripts/build/SConscript b/scripts/build/SConscript index 78adb49141..85d02f1592 100644 --- a/scripts/build/SConscript +++ b/scripts/build/SConscript @@ -41,4 +41,3 @@ exe_env.AlwaysBuild(env.Alias('INCPATH', [], incpath)) # mode:python # End: # vi:syntax=python -