Skip to content

Adds ParameterHit info to KalSeed so it can be regrown #1853

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
bonventre wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Adds ParameterHit info to KalSeed so it can be regrown #1853

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
ff8931a
Add ParameterHit info to KalSeed and allow regrowing
bonventre Jun 8, 2026
0dbe116
bug comparing signed to unsigned
bonventre Jun 8, 2026
b518083
paramhit utility function
bonventre Jun 8, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions CosmicReco/src/LineFinder_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -190,8 +190,8 @@ int LineFinder::findLine(const ComboHitCollection& shC, std::vector<StrawHitInde
Straw const& strawk = tracker->getStraw(shC[kloc].strawId());
TwoLinePCA pca( strawk.getMidPoint(), strawk.getDirection(),
ipos, newdir);
double dist = (pca.point1()-strawk.getMidPoint()).mag();
if (pca.dca() < _maxDOCA && dist < strawk.halfLength()){
double dist = (pca.point1()-strawk.getMidPoint()).dot(strawk.getDirection());
if (pca.dca() < _maxDOCA && fabs(dist) < strawk.halfLength()){
count += 1;
ll += pow(dist-shC[kloc].wireDist(),2)/shC[kloc].wireVar();
}
Expand Down
33 changes: 31 additions & 2 deletions Mu2eKinKal/inc/KKFit.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,9 @@ namespace mu2e {
using KKCALOHIT = KKCaloHit<KTRAJ>;
using KKCALOHITPTR = std::shared_ptr<KKCALOHIT>;
using KKCALOHITCOL = std::vector<KKCALOHITPTR>;
using PARAMHIT = KinKal::ParameterHit<KTRAJ>;
using PARAMHITPTR = std::shared_ptr<PARAMHIT>;
using PARAMHITCOL = std::vector<PARAMHITPTR>;
// using KKPANELHIT = KKPanelHit<KTRAJ>;
using MEAS = KinKal::Hit<KTRAJ>;
using MEASPTR = std::shared_ptr<MEAS>;
Expand All @@ -99,10 +102,12 @@ namespace mu2e {
bool makeStrawHits(Tracker const& tracker,StrawResponse const& strawresponse, BFieldMap const& kkbf,
PTRAJ const& ptraj, ComboHitCollection const& chcol, StrawHitIndexCollection const& strawHitIdxs,
KKSTRAWHITCOL& hits, KKSTRAWXINGCOL& exings) const;
void makeSeedParamHit(KTRAJ const& seedtraj, std::vector<double> const& paramconstraints, PARAMHITCOL& paramhits) const;
// regrow KKTrack components from a KalSeed
bool regrowComponents(KalSeed const& kseed, ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap,
Tracker const& tracker,Calorimeter const& calo, StrawResponse const& strawresponse,BFieldMap const& kkbf,
PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const;
PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, PARAMHITCOL& paramhits,
KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const;
std::shared_ptr<SensorLine> caloAxis(CaloCluster const& cluster, Calorimeter const& calo) const; // should come from CaloCluster TODO
bool makeCaloHit(CCPtr const& cluster, Calorimeter const& calo, PTRAJ const& pktraj, KKCALOHITCOL& hits) const;
// extend a track with a new configuration, optionally searching for and adding hits and straw material
Expand Down Expand Up @@ -251,11 +256,29 @@ namespace mu2e {
return ngood >= minNStrawHits_;
}

template <class KTRAJ> void KKFit<KTRAJ>::makeSeedParamHit(KTRAJ const& seedtraj, std::vector<double> const& paramconstraints, PARAMHITCOL& paramhits) const {
std::array<bool,KinKal::NParams()> mask = {false};
KinKal::Parameters cparams = seedtraj.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;
}
if (paramconstraints[ipar] > 0){
mask[ipar] = true;
cparams.covariance()[ipar][ipar] = paramconstraints[ipar]*paramconstraints[ipar];
}else{
cparams.covariance()[ipar][ipar] = 1.0; // otherwise inversion fails
}
}
paramhits.push_back(std::make_shared<PARAMHIT>(seedtraj.range().mid(),seedtraj,cparams,mask));
}

template <class KTRAJ> bool KKFit<KTRAJ>::regrowComponents(KalSeed const& kseed, // primary event input
ComboHitCollection const& chcol, mu2e::IndexMap const& strawindexmap, // ancillary event input
Tracker const& tracker,Calorimeter const& calo, // geometries
StrawResponse const& strawresponse,BFieldMap const& kkbf, // other conditions
PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const { // return values
PTRAJPTR& ptraj, KKSTRAWHITCOL& strawhits, KKCALOHITCOL& calohits, PARAMHITCOL& paramhits,
KKSTRAWXINGCOL& exings, DOMAINCOL& domains) const { // return values
GeomHandle<mu2e::KKMaterial> kkmat_h;
auto const& smat = kkmat_h->strawMaterial();
unsigned ngood(0), nactive(0), nsactive(0);
Expand All @@ -282,6 +305,9 @@ namespace mu2e {
if(kseed.caloCluster()){
makeCaloHit(kseed.caloCluster(),calo,*ptraj,calohits);
}
for (auto const& paramhit : kseed.paramHits()){
paramhits.push_back(std::make_shared<PARAMHIT>(paramhit.time(),*ptraj,paramhit.params(),paramhit.pmask()));
}
if(matcorr_){
// add Straw Xings for straws without hits
for(auto const& sx : kseed.straws()){
Expand Down Expand Up @@ -724,6 +750,9 @@ namespace mu2e {
calohit->unbiasedClosestApproach().tpData(),
ctres,ca.particleTraj().momentum3(ca.particleToca()));
}
for (auto const& paramhit : kktrk.paramHits()){
kseed._paramhits.emplace_back(paramhit->time(),paramhit->constraintParameters(),paramhit->constraintMask());
}
kseed._straws.reserve(kktrk.strawXings().size());
for(auto const& sxing : kktrk.strawXings()) {
// create and fill the flag
Expand Down
56 changes: 21 additions & 35 deletions Mu2eKinKal/inc/KKTrack.hh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,9 @@ namespace mu2e {
using KKCALOHIT = KKCaloHit<KTRAJ>;
using KKCALOHITPTR = std::shared_ptr<KKCALOHIT>;
using KKCALOHITCOL = std::vector<KKCALOHITPTR>;
using PARAMHIT = KinKal::ParameterHit<KTRAJ>;
using PARAMHITPTR = std::shared_ptr<PARAMHIT>;
using PARAMHITCOL = std::vector<PARAMHITPTR>;
using MEAS = KinKal::Hit<KTRAJ>;
using MEASPTR = std::shared_ptr<MEAS>;
using MEASCOL = std::vector<MEASPTR>;
Expand All @@ -60,10 +63,10 @@ namespace mu2e {
using DOMAINCOL = std::set<DOMAINPTR>;
// construct from configuration, fit environment, and hits and materials
KKTrack(Config const& config, BFieldMap const& bfield, KTRAJ const& seedtraj, PDGCode::type tpart, KKSTRAWHITCLUSTERER const& shclusterer,
KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits, std::array<double, KinKal::NParams()> constraints = {0});
KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits, PARAMHITCOL const& paramhits);
// construct from regrown constituants
KKTrack(Config const& config, BFieldMap const& bfield, PDGCode::type tpart,
PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits, DOMAINCOL& domains);
PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits, PARAMHITCOL const& paramhits, DOMAINCOL& domains);
// copy constructor
KKTrack(KKTrack<KTRAJ> const& rhs, CloneContext& context): KinKal::Track<KTRAJ>(rhs, context),
tpart_(rhs.fitParticle()),
Expand Down Expand Up @@ -118,7 +121,7 @@ namespace mu2e {

// extend the track according to new configuration, hits, and/or exings
void extendTrack(Config const& config,
KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings, KKCALOHITCOL const& calohits );
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 IPA Xing
Expand All @@ -135,6 +138,7 @@ namespace mu2e {
KKSTRAWHITCOL const& strawHits() const { return strawhits_; }
KKSTRAWHITCLUSTERCOL const& strawHitClusters() const { return strawhitclusters_; }
KKSTRAWXINGCOL const& strawXings() const { return strawxings_; }
PARAMHITCOL const& paramHits() const { return paramhits_; }
KKIPAXINGCOL const& IPAXings() const { return ipaxings_; }
KKSTXINGCOL const& STXings() const { return stxings_; }
KKCRVXINGCOL const& CRVXings() const { return crvxings_; }
Expand All @@ -153,13 +157,14 @@ namespace mu2e {
KKSTRAWHITCOL strawhits_; // straw hits used in this fit
KKSTRAWXINGCOL strawxings_; // straw material crossings used in this fit
KKCALOHITCOL calohits_; // calo hits used in this fit
PARAMHITCOL paramhits_;
KKIPAXINGCOL ipaxings_; // ipa material crossings used in extrapolation
KKSTXINGCOL stxings_; // stopping target material crossings used in extrapolation
KKCRVXINGCOL crvxings_; // crv crossings using in extrapolation
KKINTERCOL inters_; // other recorded intersections
KKSTRAWHITCLUSTERCOL strawhitclusters_; // straw hit clusters used in this fit
// utility function to convert to generic types
void convertTypes( KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings,KKCALOHITCOL const& calohits,
void convertTypes( KKSTRAWHITCOL const& strawhits, KKSTRAWXINGCOL const& strawxings,KKCALOHITCOL const& calohits, PARAMHITCOL const& paramhits,
MEASCOL& hits, EXINGCOL& exings);
// add hits to clusters
void addHitClusters(KKSTRAWHITCOL const& strawhits,KKSTRAWXINGCOL const& strawxings, MEASCOL& hits);
Expand All @@ -170,9 +175,9 @@ namespace mu2e {
KKSTRAWHITCOL const& strawhits,
KKSTRAWXINGCOL const& strawxings,
KKCALOHITCOL const& calohits,
std::array<double, KinKal::NParams()> constraints) :
PARAMHITCOL const& paramhits) :
KinKal::Track<KTRAJ>(config,bfield), tpart_(tpart), shclusterer_(shclusterer),
strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits) {
strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits), paramhits_(paramhits) {
MEASCOL hits; // polymorphic container of hits
EXINGCOL exings; // polymorphic container of detector element crossings
// add the hits to clusters, as required
Expand All @@ -183,44 +188,22 @@ namespace mu2e {
shcluster->print(std::cout,this->config().plevel_);
}
}
convertTypes(strawhits_, strawxings_, calohits_, hits, exings);
convertTypes(strawhits_, strawxings_, calohits_, paramhits_, hits, exings);

std::array<bool,KinKal::NParams()> mask = {false};
bool constraining = false;
for (size_t i=0;i<KinKal::NParams();i++){
if (constraints[i] > 0){
mask[i] = true;
constraining = true;
}
}
if (constraining){
KinKal::Parameters cparams = seedtraj.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;
}
}
for(size_t ipar=0; ipar < KinKal::NParams(); ipar++){
if (mask[ipar])
cparams.covariance()[ipar][ipar] = constraints[ipar]*constraints[ipar];
else
cparams.covariance()[ipar][ipar] = 1.0; // otherwise inversion fails
}
hits.push_back(std::make_shared<KinKal::ParameterHit<KTRAJ>>(seedtraj.range().mid(),seedtraj,cparams,mask));
}
// now fit these
this->fit(hits, exings, seedtraj);
}


template <class KTRAJ> KKTrack<KTRAJ>::KKTrack(Config const& config, BFieldMap const& bfield, PDGCode::type tpart,
PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits, DOMAINCOL& domains) :
PKTRAJPTR& fittraj, KKSTRAWHITCOL& strawhits, KKSTRAWXINGCOL& strawxings, KKCALOHITCOL& calohits,
PARAMHITCOL const& paramhits, DOMAINCOL& domains) :
KinKal::Track<KTRAJ>(config,bfield), tpart_(tpart),shclusterer_(StrawIdMask::none,0,0.0),
strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits) {
strawhits_(strawhits), strawxings_(strawxings), calohits_(calohits), paramhits_(paramhits) {
// convert types
MEASCOL hits; // polymorphic container of hits
EXINGCOL exings; // polymorphic container of detector element crossings
convertTypes(strawhits_, strawxings_, calohits_, hits, exings);
convertTypes(strawhits_, strawxings_, calohits_, paramhits_, hits, exings);
this->fit(hits,exings,domains,fittraj);
// add ParameterHit, intersect, extrapolate, ... TODO
}
Expand Down Expand Up @@ -265,11 +248,13 @@ namespace mu2e {
KKSTRAWHITCOL const& strawhits,
KKSTRAWXINGCOL const& strawxings,
KKCALOHITCOL const& calohits,
PARAMHITCOL const& paramhits,
MEASCOL& hits, EXINGCOL& exings) {
hits.reserve(strawhits_.size() + calohits_.size());
hits.reserve(strawhits_.size() + calohits_.size()+paramhits_.size());
exings.reserve(strawxings_.size());
for(auto const& strawhit : strawhits)hits.emplace_back(std::static_pointer_cast<MEAS>(strawhit));
for(auto const& calohit : calohits)hits.emplace_back(std::static_pointer_cast<MEAS>(calohit));
for(auto const& paramhit : paramhits)hits.emplace_back(std::static_pointer_cast<MEAS>(paramhit));
for(auto const& strawxing : strawxings)exings.emplace_back(std::static_pointer_cast<EXING>(strawxing));
}

Expand All @@ -294,7 +279,8 @@ namespace mu2e {
}
if(nhit != strawhits_.size()+strawhits.size()) std::cout << "cluster hit sum doesn't match " << nhit << " " << strawhits_.size() << std::endl;
}
convertTypes(strawhits,strawxings,calohits,hits,exings);
PARAMHITCOL paramhits;
convertTypes(strawhits,strawxings,calohits,paramhits,hits,exings);
this->extend(config,hits,exings);
// store the new hits
strawhits_.reserve(strawhits_.size()+strawhits.size());
Expand Down
29 changes: 23 additions & 6 deletions Mu2eKinKal/src/CentralHelixFit_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,9 @@ namespace mu2e {
using KKCALOHIT = KKCaloHit<KTRAJ>;
using KKCALOHITPTR = std::shared_ptr<KKCALOHIT>;
using KKCALOHITCOL = std::vector<KKCALOHITPTR>;
using PARAMHIT = KinKal::ParameterHit<KTRAJ>;
using PARAMHITPTR = std::shared_ptr<PARAMHIT>;
using PARAMHITCOL = std::vector<PARAMHITPTR>;
using KKFIT = KKFit<KTRAJ>;
using KinKal::VEC3;
using KinKal::DMAT;
Expand All @@ -113,7 +116,7 @@ namespace mu2e {
fhicl::OptionalAtom<double> fixedBField { Name("ConstantBField"), Comment("Constant BField value") };
fhicl::Atom<double> seedMom { Name("SeedMomentum"), Comment("Seed momentum") };
fhicl::Atom<int> seedCharge { Name("SeedCharge"), Comment("Seed charge MAGNITUDE, in electron charge units") };
fhicl::Sequence<float> parconst { Name("ParameterConstraints"), Comment("External constraint on parameters to seed values (rms, various units)") };
fhicl::Sequence<double> parconst { Name("ParameterConstraints"), Comment("External constraint on parameters to seed values (rms, various units)") };
fhicl::Sequence<std::string> sampleSurfaces { Name("SampleSurfaces"), Comment("When creating the KalSeed, sample the fit at these surfaces") };
fhicl::Atom<bool> sampleInRange { Name("SampleInRange"), Comment("Require sample times to be inside the fit trajectory time range") };
fhicl::Atom<bool> sampleInBounds { Name("SampleInBounds"), Comment("Require sample intersection point be inside surface bounds (within tolerance)") };
Expand Down Expand Up @@ -165,14 +168,15 @@ namespace mu2e {
bool fixedfield_; //
double seedMom_;
int seedCharge_;
std::vector<double> paramconstraints_;
double intertol_; // surface intersection tolerance (mm)
double sampletbuff_; // simple time buffer; replace this with extrapolation TODO
bool useFitCharge_; // Set the PDG particle to agree with the fit charge
double minCenterRho_; // min center distance to z axis
bool sampleinrange_, sampleinbounds_; // require samples to be in range or on surface
SurfaceIdCollection ssids_;
KinKalGeom::SurfacePairCollection surfacess_to_sample_; // surfaces to sample the fit
std::array<double,KinKal::NParams()> paramconstraints_;
bool constraining_;
};

CentralHelixFit::CentralHelixFit(const Parameters& settings) : art::EDProducer{settings},
Expand All @@ -189,6 +193,7 @@ namespace mu2e {
fixedfield_(false),
seedMom_(settings().modSettings().seedMom()),
seedCharge_(settings().modSettings().seedCharge()),
paramconstraints_(settings().modSettings().parconst()),
intertol_(settings().modSettings().interTol()),
sampletbuff_(settings().modSettings().sampleTBuff()),
useFitCharge_(settings().modSettings().useFitCharge()),
Expand All @@ -208,9 +213,15 @@ namespace mu2e {
for(size_t ipar=0;ipar < seederrors.size(); ++ipar){
seedcov_[ipar][ipar] = seederrors[ipar]*seederrors[ipar];
}
auto const& parerrors = settings().modSettings().parconst();
for (size_t ipar=0;ipar<KinKal::NParams();ipar++)
paramconstraints_[ipar] = parerrors[ipar];
constraining_ = false;
if (paramconstraints_.size() == KinKal::NParams()){
for (size_t ipar=0;ipar<KinKal::NParams();ipar++){
if (paramconstraints_[ipar] > 0)
constraining_ = true;
}
}else if (paramconstraints_.size() > 0){
throw cet::exception("RECO")<<"mu2e::CentralHelixFit: Parameter constraint configuration error"<< endl;
}

if(print_ > 0) std::cout << "Fit " << config_ << "Extension " << exconfig_;
double bz(0.0);
Expand Down Expand Up @@ -330,8 +341,14 @@ namespace mu2e {

try {
seedtraj.range() = kkfit_.range(strawhits,calohits,strawxings);
// create parameter constraint, use seedtrajectory parameters
PARAMHITCOL paramhits;
if (constraining_){
kkfit_.makeSeedParamHit(seedtraj,paramconstraints_,paramhits);
}

// create and fit the track
auto ktrk = make_unique<KKTRK>(config_,*kkbf_,seedtraj,fitpart,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramconstraints_);
auto ktrk = make_unique<KKTRK>(config_,*kkbf_,seedtraj,fitpart,kkfit_.strawHitClusterer(),strawhits,strawxings,calohits,paramhits);
// Check the fit
auto goodfit = goodFit(*ktrk);
// if we have an extension schedule, extend.
Expand Down
Loading