From 8721501e8147b1a66922adc6033d2cd8a1e74c41 Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Tue, 23 Jun 2026 17:07:13 -0700 Subject: [PATCH 1/3] All regrowing from a Ptr collection --- Mu2eKinKal/fcl/prolog.fcl | 9 ++--- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 42 +++++++++++++++++------- 2 files changed, 33 insertions(+), 18 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 11fe39cf88..4de1b85b24 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -527,8 +527,8 @@ Mu2eKinKal : { # Additional width is needed here as straight tracks have worse time resolution MaxCaloClusterDt: 8 # save the full trajectory - SaveTrajectory: Full - SaveDomains : false # meaningless for Linefit + SaveTrajectory: Full + SaveDomains : false # meaningless for Linefit } FitSettings : @local::Mu2eKinKal.CHSEEDFIT ExtensionSettings : @local::Mu2eKinKal.CHDRIFTEXT @@ -608,7 +608,6 @@ Mu2eKinKal.producers.KKUeSeedFit.FitDirection : @local::FitDir.downstream Mu2eKinKal.producers.KKUmuSeedFit.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKUmuSeedFit.FitDirection : @local::FitDir.downstream -# save trajectories in the Detector region for downstream fits, just the T0 segment for the rest Mu2eKinKal.producers.KKDe.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKDe.FitDirection : @local::FitDir.downstream Mu2eKinKal.producers.KKDmu.ModuleSettings.FitParticle : @local::Particle.muminus @@ -617,9 +616,5 @@ Mu2eKinKal.producers.KKUe.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKUe.FitDirection : @local::FitDir.upstream Mu2eKinKal.producers.KKUmu.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKUmu.FitDirection : @local::FitDir.upstream -# propagate upstream fits back to the tracker -Mu2eKinKal.producers.KKUe.ExtrapolationSettings.BackToTracker : true -Mu2eKinKal.producers.KKUmu.ExtrapolationSettings.BackToTracker : true - # END_PROLOG diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index 085ed2166d..216b042026 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -84,9 +84,10 @@ namespace mu2e { using Comment = fhicl::Comment; struct RegrowLoopHelixConfig { fhicl::Atom debug{Name("debug"), Comment("Debug level"), 0}; - fhicl::Atom kalSeedCollection {Name("KalSeedCollection"), Comment("KalSeed collection to process") }; + fhicl::OptionalAtom kalSeedCollection {Name("KalSeedCollection"), Comment("KalSeed collection to process") }; + fhicl::OptionalAtom kalSeedPtrCollection {Name("KalSeedPtrCollection"), Comment("Collection of KalSeedPtrs to process") }; fhicl::Atom comboHitCollection {Name("ComboHitCollection"), Comment("Reduced ComboHit collection") }; - fhicl::Atom caloClusterCollection {Name("CaloClusterCollection"), Comment("CaloCluster collection ") }; + fhicl::OptionalAtom caloClusterCollection {Name("CaloClusterCollection"), Comment("CaloCluster collection ") }; fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; @@ -142,12 +143,14 @@ namespace mu2e { Config config_; // refit configuration object, containing the fit schedule KKFIT kkfit_; art::ProductToken kseedcol_T_; + art::ProductToken kseedptrcol_T_; art::ProductToken chcol_T_; art::ProductToken cccol_T_; art::ProductToken indexmap_T_; art::InputTag ksmca_T_; bool fillMCAssns_; bool extend_; + bool has_cccol_, has_kseedcol_, has_kseedptrcol_; std::unique_ptr extrap_; }; @@ -155,12 +158,13 @@ namespace mu2e { debug_(settings().debug()), config_(Mu2eKinKal::makeConfig(settings().fitSettings())), kkfit_(settings().kkfitSettings()), - kseedcol_T_(consumes(settings().kalSeedCollection())), chcol_T_(consumes(settings().comboHitCollection())), - cccol_T_(mayConsume(settings().caloClusterCollection())), indexmap_T_(consumes(settings().indexMap())), fillMCAssns_(settings().kalSeedMCAssns(ksmca_T_)), - extend_(settings().extend()) + extend_(settings().extend()), + has_cccol_(settings().caloClusterCollection()), + has_kseedcol_(settings().kalSeedCollection()), + has_kseedptrcol_(settings().kalSeedPtrCollection()) { produces(); produces(); @@ -169,6 +173,9 @@ namespace mu2e { consumes(ksmca_T_); produces (); } + if(has_kseedcol_)kseedcol_T_ = art::ProductToken(mayConsume(settings().kalSeedCollection().value())); + if(has_kseedptrcol_)kseedptrcol_T_ = art::ProductToken(mayConsume(settings().kalSeedPtrCollection().value())); + if(has_cccol_) cccol_T_ = art::ProductToken(mayConsume(settings().caloClusterCollection().value())); } void RegrowLoopHelix::beginRun(art::Run& run) @@ -186,8 +193,6 @@ namespace mu2e { GeomHandle nominalTracker_h; GeomHandle calo_h; // find input event data - auto kseed_H = event.getValidHandle(kseedcol_T_); - const auto& kseedcol = *kseed_H; auto ch_H = event.getValidHandle(chcol_T_); const auto& chcol = *ch_H; auto cc_H = event.getHandle(cccol_T_); @@ -207,7 +212,23 @@ namespace mu2e { ksmca = std::unique_ptr(new KalSeedMCAssns); } size_t iseed(0); - for (auto const& kseed : kseedcol) { + if(!(has_kseedcol_ || has_kseedptrcol_) || (has_kseedcol_ && has_kseedptrcol_)) + throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: exactly 1 of KalSeedCollection or KalSeedPtrCollection must be specified" << endl; + + KalSeedPtrCollection kseedptrs; + if(has_kseedptrcol_){ + auto kseedptr_H = event.getValidHandle(kseedptrcol_T_); + kseedptrs = *kseedptr_H; + } else if(has_kseedcol_){ + auto kseed_H = event.getValidHandle(kseedcol_T_); + const auto& kseedcol = *kseed_H; + for(size_t iks = 0; iks < kseedcol.size(); ++iks){ + kseedptrs.emplace_back(kseed_H,iks); + } + } + + for (auto const& kseedptr : kseedptrs) { + auto const& kseed = *kseedptr; if(!kseed.loopHelixFit())throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: passed KalSeed from non-LoopHelix fit " << endl; // regrow the components from the seed PKTRAJPTR trajptr = kseed.loopHelixFitTrajectory(); @@ -252,14 +273,13 @@ namespace mu2e { if(fillMCAssns_){ // find the MC assns auto ksmcai = (*ksmca_H)[iseed]; - auto origksp = art::Ptr(kseed_H,iseed); // test this is the right ptr - if(ksmcai.first != origksp)throw cet::exception("Reco")<<"mu2e::RegrowLoopHelix: wrong KalSeed ptr"<< std::endl; + if(ksmcai.first != kseedptr)throw cet::exception("Reco")<<"mu2e::RegrowLoopHelix: wrong KalSeed ptr"<< std::endl; auto mcseedp = ksmcai.second; auto rgksp = art::Ptr(KalSeedCollectionPID,rgkseedcol->size()-1,KalSeedCollectionGetter); ksmca->addSingle(rgksp,mcseedp); // add the original too - ksmca->addSingle(origksp,mcseedp); + ksmca->addSingle(kseedptr,mcseedp); } if(debug_ > 5)static_cast&>(ktrk->fitTraj()).print(std::cout,2); ktrkcol->push_back(ktrk.release()); From f8fceca0d25168d1febfb964901643188772106f Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Wed, 1 Jul 2026 22:12:33 -0700 Subject: [PATCH 2/3] Optionally regrow KalSeedMCs. Fix MCAssns in regrowing --- Mu2eKinKal/fcl/prolog.fcl | 24 ++++++++++- Mu2eKinKal/src/RegrowLoopHelix_module.cc | 54 ++++++++++++++++++------ TrkReco/src/SelectReflections_module.cc | 2 +- 3 files changed, 64 insertions(+), 16 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 4de1b85b24..12ae1391c6 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -485,8 +485,8 @@ Mu2eKinKal : { @table::Mu2eKinKal.KKFIT SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] # these are additional surfaces; surfaces used in extrapolation are also sampled # save trajectories in the Detector region - SaveDomains : true - SaveTrajectory : Detector + SaveDomains : true + SaveTrajectory : Detector } FitSettings : { @table::Mu2eKinKal.LHSEEDFIT @@ -573,6 +573,26 @@ Mu2eKinKal : { Mu2eKinKal : { @table::Mu2eKinKal + + RegrowLH : { + module_type : RegrowLoopHelix + CaloClusterCollection : "CaloClusterMaker" + StrawDigiIndexMap : "SelectReco:StrawDigiMap" + RefitSettings : { + @table::Mu2eKinKal.REGROW + BFieldCorrection : true + } + KKFitSettings: { + @table::Mu2eKinKal.KKFIT + SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] + AddHits : false + SaveDomains : true + SaveTrajectory : Detector + } + ExtrapolationSettings : @local::Mu2eKinKal.LHDRIFTXTRAP + Extend: false + } + producers : { # seed fits: these don't use drift information or BField corrections KLSeedFit: @local::Mu2eKinKal.KLSeedFit diff --git a/Mu2eKinKal/src/RegrowLoopHelix_module.cc b/Mu2eKinKal/src/RegrowLoopHelix_module.cc index 216b042026..4a84d6d19d 100644 --- a/Mu2eKinKal/src/RegrowLoopHelix_module.cc +++ b/Mu2eKinKal/src/RegrowLoopHelix_module.cc @@ -90,6 +90,7 @@ namespace mu2e { fhicl::OptionalAtom caloClusterCollection {Name("CaloClusterCollection"), Comment("CaloCluster collection ") }; fhicl::Atom indexMap {Name("StrawDigiIndexMap"), Comment("Map between original and reduced ComboHits") }; fhicl::OptionalAtom kalSeedMCAssns {Name("KalSeedMCAssns"), Comment("Association to KalSeedMC. If set, regrown KalSeeds will be associated with the same KalSeedMC as the original") }; + fhicl::OptionalAtom copyKalSeedMCs { Name("CopyKalSeedMCs"), Comment("If set, and KalSeedMCs are referenced, create a deep copy of the input") }; fhicl::Table kkfitSettings { Name("KKFitSettings") }; fhicl::Table fitSettings { Name("RefitSettings") }; fhicl::Atom extend {Name("Extend"), Comment("Extend the fit") }; @@ -148,7 +149,7 @@ namespace mu2e { art::ProductToken cccol_T_; art::ProductToken indexmap_T_; art::InputTag ksmca_T_; - bool fillMCAssns_; + bool fillMCAssns_, copyKalSeedMCs_; bool extend_; bool has_cccol_, has_kseedcol_, has_kseedptrcol_; std::unique_ptr extrap_; @@ -172,6 +173,12 @@ namespace mu2e { if( fillMCAssns_){ consumes(ksmca_T_); produces (); + // require KalSeedMC copy choice + if(!(settings().copyKalSeedMCs(copyKalSeedMCs_))) + throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: Specify if KalSeedMC Collection is deep copied or not" << endl; + if(copyKalSeedMCs_){ + produces(); + } } if(has_kseedcol_)kseedcol_T_ = art::ProductToken(mayConsume(settings().kalSeedCollection().value())); if(has_kseedptrcol_)kseedptrcol_T_ = art::ProductToken(mayConsume(settings().kalSeedPtrCollection().value())); @@ -200,18 +207,23 @@ namespace mu2e { const auto& indexmap = *indexmap_H; auto KalSeedCollectionPID = event.getProductID(); auto KalSeedCollectionGetter = event.productGetter(KalSeedCollectionPID); + art::ProductID KalSeedMCCollectionPID; art::Handle ksmca_H; // create outputs unique_ptr ktrkcol(new KKTRKCOL ); unique_ptr rgkseedcol(new KalSeedCollection ); - std::unique_ptr ksmca; + std::unique_ptr rgksmca; + std::unique_ptr rgksmcc; // deal with MC if(fillMCAssns_){ ksmca_H = event.getHandle(ksmca_T_); if(!ksmca_H)throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: No KalSeedMCAssns found" << endl; - ksmca = std::unique_ptr(new KalSeedMCAssns); + rgksmca = std::unique_ptr(new KalSeedMCAssns); + if(copyKalSeedMCs_){ + rgksmcc = std::unique_ptr(new KalSeedMCCollection); + KalSeedMCCollectionPID = event.getProductID(); + } } - size_t iseed(0); if(!(has_kseedcol_ || has_kseedptrcol_) || (has_kseedcol_ && has_kseedptrcol_)) throw cet::exception("RECO")<<"mu2e::RegrowLoopHelix: exactly 1 of KalSeedCollection or KalSeedPtrCollection must be specified" << endl; @@ -270,16 +282,32 @@ namespace mu2e { fitflag.merge(TrkFitFlag::Regrown); auto rgks = kkfit_.createSeed(*ktrk,fitflag,*calo_h,*nominalTracker_h); rgkseedcol->push_back(rgks); + auto rgksp = art::Ptr(KalSeedCollectionPID,rgkseedcol->size()-1,KalSeedCollectionGetter); if(fillMCAssns_){ // find the MC assns - auto ksmcai = (*ksmca_H)[iseed]; + auto ksmca = *ksmca_H; + auto ksmcai= ksmca.end(); + for(auto ksmci= ksmca.begin(); ksmci != ksmca.end(); ++ksmci){ + if(ksmci->first == kseedptr){ + ksmcai = ksmci; + break; + } + } // test this is the right ptr - if(ksmcai.first != kseedptr)throw cet::exception("Reco")<<"mu2e::RegrowLoopHelix: wrong KalSeed ptr"<< std::endl; - auto mcseedp = ksmcai.second; - auto rgksp = art::Ptr(KalSeedCollectionPID,rgkseedcol->size()-1,KalSeedCollectionGetter); - ksmca->addSingle(rgksp,mcseedp); - // add the original too - ksmca->addSingle(kseedptr,mcseedp); + if(ksmcai == ksmca.end())throw cet::exception("Reco")<<"mu2e::RegrowLoopHelix: can't find MC associated with KalSeed" << std::endl; + if(copyKalSeedMCs_){ + // deep-copy the KalSeedMC + rgksmcc->push_back(*ksmcai->second); + auto KalSeedMCCollectionGetter = event.productGetter(KalSeedMCCollectionPID); + auto mcseedp = art::Ptr(KalSeedMCCollectionPID,rgksmcc->size()-1,KalSeedMCCollectionGetter); + rgksmca->addSingle(rgksp,mcseedp); + } else { + // just reference the original KalSeedMC + auto mcseedp = ksmcai->second; + rgksmca->addSingle(rgksp,mcseedp); + // add the original too + rgksmca->addSingle(kseedptr,mcseedp); + } } if(debug_ > 5)static_cast&>(ktrk->fitTraj()).print(std::cout,2); ktrkcol->push_back(ktrk.release()); @@ -289,12 +317,12 @@ namespace mu2e { if(debug_ > 3)std::cout<< "original seed had " << kseed.hits().size() << " hits, NDOF " << kseed.nDOF()<< " fitcon " << kseed.fitConsistency() << " and status " << kseed.status() << std::endl; } } - ++iseed; } // store output event.put(move(ktrkcol)); event.put(move(rgkseedcol)); - if(fillMCAssns_)event.put(move(ksmca)); + if(fillMCAssns_)event.put(move(rgksmca)); + if(copyKalSeedMCs_)event.put(move(rgksmcc)); } void RegrowLoopHelix::endJob() diff --git a/TrkReco/src/SelectReflections_module.cc b/TrkReco/src/SelectReflections_module.cc index acdcd13058..f88190bed9 100644 --- a/TrkReco/src/SelectReflections_module.cc +++ b/TrkReco/src/SelectReflections_module.cc @@ -129,7 +129,7 @@ namespace mu2e { if(matches.size() >0 ){ ibest = 0; if(matches.size()>1){ - if(debug_ > 1) std::cout << "Selecting best reflection pair from " << matches.size() << " candidates " << std::endl; + if(debug_ > 0) std::cout << "Selecting best reflection pair from " << matches.size() << " candidates " << std::endl; double value = std::numeric_limits::max(); for (size_t imatch = 0; imatch < matches.size(); ++imatch) { auto const& match = matches[imatch]; From 787018383baa92611e8fcfecc0024396acc9e19e Mon Sep 17 00:00:00 2001 From: David Nathan Brown Date: Thu, 2 Jul 2026 16:55:15 -0700 Subject: [PATCH 3/3] Add diagnostics. Update config to save everything needed for reflection constraint --- Mu2eKinKal/fcl/prolog.fcl | 12 ++++++-- TrkReco/fcl/prolog.fcl | 7 +++-- TrkReco/src/SelectReflections_module.cc | 39 +++++++++++++++++-------- 3 files changed, 40 insertions(+), 18 deletions(-) diff --git a/Mu2eKinKal/fcl/prolog.fcl b/Mu2eKinKal/fcl/prolog.fcl index 12ae1391c6..0d64d20281 100644 --- a/Mu2eKinKal/fcl/prolog.fcl +++ b/Mu2eKinKal/fcl/prolog.fcl @@ -587,9 +587,12 @@ Mu2eKinKal : { SampleSurfaces : ["ST_Outer","ST_Front","ST_Back"] AddHits : false SaveDomains : true - SaveTrajectory : Detector + SaveTrajectory : Full + } + ExtrapolationSettings :{ + @table::Mu2eKinKal.LHDRIFTXTRAP + BackToTracker: true } - ExtrapolationSettings : @local::Mu2eKinKal.LHDRIFTXTRAP Extend: false } @@ -636,5 +639,8 @@ Mu2eKinKal.producers.KKUe.ModuleSettings.FitParticle : @local::Particle.eminus Mu2eKinKal.producers.KKUe.FitDirection : @local::FitDir.upstream Mu2eKinKal.producers.KKUmu.ModuleSettings.FitParticle : @local::Particle.muminus Mu2eKinKal.producers.KKUmu.FitDirection : @local::FitDir.upstream -# +# propagate upstream fits back to the tracker +Mu2eKinKal.producers.KKUe.ExtrapolationSettings.BackToTracker : true +Mu2eKinKal.producers.KKUmu.ExtrapolationSettings.BackToTracker : true + END_PROLOG diff --git a/TrkReco/fcl/prolog.fcl b/TrkReco/fcl/prolog.fcl index 49aea5119f..6f6f1f8c64 100644 --- a/TrkReco/fcl/prolog.fcl +++ b/TrkReco/fcl/prolog.fcl @@ -138,8 +138,9 @@ TrkReco: { @table::TrkReco # standard config for reflection selection SelectReflections: { module_type : SelectReflections - MaxDeltaT : 20 # ns - MaxDeltaP : 10 # MeV + Surface : TT_Front + MaxDeltaT : 15 # ns + MaxDeltaP : 15 # MeV # optionally only merge selected candidates # a simple selector is used here, but any selector implemented as a tool can be used Selector : { @@ -150,7 +151,7 @@ TrkReco: { @table::TrkReco MinDeltaNHitFraction : 0.0 MinActiveHits : 0 } - SelectBestPair : 0 # select highest-momentum match + SelectBestPair : 3 # select pair with highest # of active hits debugLevel : 0 } } diff --git a/TrkReco/src/SelectReflections_module.cc b/TrkReco/src/SelectReflections_module.cc index f88190bed9..5612bd784c 100644 --- a/TrkReco/src/SelectReflections_module.cc +++ b/TrkReco/src/SelectReflections_module.cc @@ -34,31 +34,36 @@ namespace mu2e { class SelectReflections : public art::EDFilter { public: - enum bestpair {mom=0, deltat,deltap}; // selection options for defining the best pair + enum bestpair {mom=0, deltat,deltap, nactive}; // selection options for defining the best pair using Name=fhicl::Name; using Comment=fhicl::Comment; struct Config { fhicl::Atom debug{ Name("debugLevel"), Comment("Debug Level"), 0}; fhicl::Atom upstreamTag {Name("UpstreamKalSeedCollection"), Comment("Upstream KalSeed collection") }; fhicl::Atom downstreamTag {Name("DownstreamKalSeedCollection"), Comment("Downstream KalSeed collection") }; - fhicl::Atom maxdt{ Name("MaxDeltaT"), Comment("Maximum time difference at tracker entrance")}; - fhicl::Atom maxdp{ Name("MaxDeltaP"), Comment("Maximum scalar momentum difference at tracker entrance")}; + fhicl::Atom surface{ Name("Surface"), Comment("Surface to compare fits at")}; + fhicl::Atom maxdt{ Name("MaxDeltaT"), Comment("Maximum time difference at comparison surface")}; + fhicl::Atom maxdp{ Name("MaxDeltaP"), Comment("Maximum scalar momentum difference at comparison surface")}; fhicl::DelegatedParameter selector{ Name("Selector"), Comment("Selector parameters")}; fhicl::Atom selectbest{ Name("SelectBestPair"), Comment("Method to select the best pair if multiple pairs are found ")}; }; using Parameters = art::EDFilter::Table; explicit SelectReflections(const Parameters& config); bool filter(art::Event& evt) override; + void endJob() override; private: int debug_; + SurfaceId sid_; double maxdt_, maxdp_; art::ProductToken upstreamtoken_, downstreamtoken_; std::unique_ptr selector_; bestpair selbest_; + unsigned nevts_=0, nref_=0, nmultref_=0; }; SelectReflections::SelectReflections(const Parameters& config) : art::EDFilter{config}, debug_(config().debug()), + sid_(config().surface()), maxdt_(config().maxdt()), maxdp_(config().maxdp()), upstreamtoken_ { consumes (config().upstreamTag()) }, @@ -69,7 +74,7 @@ namespace mu2e { } bool SelectReflections::filter(art::Event& event) { - static const SurfaceId trkfront (SurfaceIdDetail::TT_Front); + ++nevts_; // create output std::unique_ptr mkseeds(new KalSeedPtrCollection); auto const& upksch = event.getValidHandle(upstreamtoken_); @@ -79,7 +84,7 @@ namespace mu2e { if(debug_ > 1) std::cout << "Upstream " << upksc.size() << " , Downstream " << downksc.size() << " KalSeeds" << std::endl; bool keep(false); if(upksc.size() > 0 && downksc.size() > 0){ - std::vector> matches; // matched up and downstream track, with (downstream) mom, dt and dmom + std::vector> matches; // matched up and downstream track, with (downstream) mom, dt and dmom for(size_t iup = 0; iup select(upks)){ @@ -88,7 +93,7 @@ namespace mu2e { auto uptrkiinter = upks.intersections().end(); for(auto upiinter = upks.intersections().begin(); upiinter != upks.intersections().end(); ++upiinter){ auto const& upinter = *upiinter; - if(upinter.surfaceId() == trkfront && upinter.momentum3().Z() > 0.0){ // correct surface and direction + if(upinter.surfaceId() == sid_ && upinter.momentum3().Z() > 0.0){ // correct surface and direction if(debug_ > 1) std::cout << "Found upstream intersection mom " << upinter.momentum3() << " time " << upinter.time() << std::endl; uptrkiinter = upiinter; break; @@ -105,7 +110,7 @@ namespace mu2e { auto downtrkiinter = downks.intersections().end(); for(auto downiinter = downks.intersections().begin(); downiinter != downks.intersections().end(); ++downiinter){ auto const& downinter = *downiinter; - if(downinter.surfaceId() == trkfront && downinter.momentum3().Z() > 0.0){ // correct surface and direction + if(downinter.surfaceId() == sid_ && downinter.momentum3().Z() > 0.0){ // correct surface and direction if(debug_ > 1) std::cout << "Found downstream intersection mom " << downinter.momentum3() << " time " << downinter.time() << std::endl; downtrkiinter = downiinter; break; @@ -116,8 +121,8 @@ namespace mu2e { double dt = fabs(uptrkiinter->time() - downtrkiinter->time()); double dmom = fabs(uptrkiinter->mom() - downtrkiinter->mom()); if( dt < maxdt_ && dmom < maxdp_){ - if(debug_ > 1) std::cout << "Found matching track pair " << std::endl; - matches.emplace_back(iup,idown,downtrkiinter->mom(),dt,dmom); + if(debug_ > 1) std::cout << "Found matching track pair, dt " << dt << " dmom " << dmom << std::endl; + matches.emplace_back(iup,idown,downtrkiinter->mom(),dt,dmom,upks.nHits()+downks.nHits()); } } } @@ -127,21 +132,27 @@ namespace mu2e { // if there are >1 matches select the best int ibest = -1; if(matches.size() >0 ){ + ++nref_; ibest = 0; if(matches.size()>1){ - if(debug_ > 0) std::cout << "Selecting best reflection pair from " << matches.size() << " candidates " << std::endl; + ++nmultref_; + if(debug_ > 0) std::cout << "Selecting best reflection pair from " << matches.size() << " candidates " << " according to " << selbest_ << " event " << event.id() << std::endl; double value = std::numeric_limits::max(); for (size_t imatch = 0; imatch < matches.size(); ++imatch) { auto const& match = matches[imatch]; - if(selbest_ == mom && -std::get<2>(match) < value){ + if(debug_ > 1)std::cout << "Match " << imatch << " has downstream momentum " << std::get<2>(match) << " dt " << std::get<3>(match) << " dp " << std::get<4>(match) << " nactive " << std::get<5>(match) << std::endl; + if(selbest_ == mom && std::get<2>(match) > value){ ibest = imatch; - value = -std::get<2>(match); // sign to make highest momentum best + value = std::get<2>(match); } else if(selbest_ == deltat && std::get<3>(match) < value){ ibest = imatch; value = std::get<3>(match); } else if(selbest_ == deltap && std::get<4>(match) < value){ ibest = imatch; value = std::get<4>(match); + } else if(selbest_ == nactive && std::get<5>(match) > value){ + ibest = imatch; + value = std::get<5>(match); } } } @@ -158,5 +169,9 @@ namespace mu2e { event.put(std::move(mkseeds)); return keep; } + + void SelectReflections::endJob() { + std::cout << moduleDescription().moduleLabel() << " processed " << nevts_ << " events and found " << nref_ << " reflections with " << nmultref_ << " multiple reflections" << std::endl; + } } DEFINE_ART_MODULE(mu2e::SelectReflections)