diff --git a/Analyses/CMakeLists.txt b/Analyses/CMakeLists.txt index 1520f73fb9..4aedb8b320 100644 --- a/Analyses/CMakeLists.txt +++ b/Analyses/CMakeLists.txt @@ -1,3 +1,16 @@ +cet_build_plugin(AprHelixTruthAna art::module + REG_SOURCE src/AprHelixTruthAna_module.cc + LIBRARIES REG + art_root_io::TFileService_service + + Offline::BFieldGeom + Offline::GeometryService + Offline::MCDataProducts + Offline::Mu2eUtilities + Offline::RecoDataProducts + ROOT::Tree +) + cet_build_plugin(BFieldPlotter art::module REG_SOURCE src/BFieldPlotter_module.cc LIBRARIES REG diff --git a/Analyses/src/AprHelixTruthAna_module.cc b/Analyses/src/AprHelixTruthAna_module.cc new file mode 100644 index 0000000000..a713f860aa --- /dev/null +++ b/Analyses/src/AprHelixTruthAna_module.cc @@ -0,0 +1,332 @@ +/////////////////////////////////////////////////////////////////////////////// +// AprHelixTruthAna: per-event TTree comparing helix seeds (e.g. the trigger APR +// helices) with the MC truth of the main signal particle, to study helix finding +// efficiency and wrong-lambda / wrong-circle helices. Optionally records the +// decisions of trigger paths of the current process (triggerPaths), so trigger +// efficiencies can be split by the truth of the signal particle. +// +// The signal particle is the highest-momentum e+- or mu+- with >= minHits straw hits +// (via the StrawDigiMCs). Its truth at the earliest StrawGasStep gives p, pT, pz, the +// circle center and radius, and lambda = pz/(c B). Its combo hits in the +// helix finder's input collection give the z span and the largest z gap. +// For each helix: lambda, radius, center, N(straw hits), N(signal straw hits), +// z span and largest z gap of its hits. +/////////////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileService.h" +#include "fhiclcpp/types/Atom.h" + +#include "Offline/BFieldGeom/inc/BFieldManager.hh" +#include "Offline/DataProducts/inc/PDGCode.hh" +#include "Offline/GeometryService/inc/DetectorSystem.hh" +#include "Offline/GeometryService/inc/GeomHandle.hh" +#include "Offline/MCDataProducts/inc/StrawDigiMC.hh" +#include "Offline/RecoDataProducts/inc/ComboHit.hh" +#include "Offline/RecoDataProducts/inc/HelixSeed.hh" +#include "Offline/Mu2eUtilities/inc/TriggerResultsNavigator.hh" +#include "canvas/Persistency/Common/TriggerResults.h" +#include "fhiclcpp/types/Sequence.h" + +#include "CLHEP/Units/PhysicalConstants.h" +#include "TTree.h" + +#include +#include +#include +#include + +namespace mu2e { + + class AprHelixTruthAna : public art::EDAnalyzer { + public: + struct Config { + using Name = fhicl::Name; + using Comment = fhicl::Comment; + fhicl::Atom hsTag {Name("helixSeedTag"), Comment("helix seed collection")}; + fhicl::Atom chTag {Name("comboHitTag"), Comment("helix finder input combo hits")}; + fhicl::Atom sdmcTag {Name("strawDigiMCTag"), Comment("StrawDigiMC collection")}; + fhicl::Atom minHits {Name("minHits"), Comment("min straw hits of the signal particle"), 5}; + fhicl::Sequence paths {Name("triggerPaths"), Comment("trigger paths of this process to record (pass = 1)"), std::vector{}}; + }; + using Parameters = art::EDAnalyzer::Table; + explicit AprHelixTruthAna(const Parameters& config); + void beginJob() override; + void beginRun(const art::Run&) override; + void analyze(const art::Event& event) override; + + private: + // z span, largest gap and N(distinct z planes) of a set of z values + static void zInfo(std::vector z, float& zmin, float& zmax, float& gap, int& nz); + + art::InputTag _hsTag, _chTag, _sdmcTag; + int _minHits; + std::vector _paths; + std::vector _pass; // per entry of _paths + float _bzConv = 0.; // MeV/c per mm + TTree* _tree = nullptr; + + // event + int _run, _subrun, _event; + bool _hasHelices; // the helix collection exists (its producer ran) + // truth of the signal particle + int _pdg, _nsh; // pdg, N(straw hits) + float _p, _pt, _pz, _rtrue, _lambdatrue, _xctrue, _yctrue, _d0true; + int _nchsig, _nzsig; // signal combo hits in the input collection, distinct z + float _zminsig, _zmaxsig, _gapsig; + // helices + std::vector _hlambda, _hradius, _hxc, _hyc, _hd0, _hzmin, _hzmax, _hgap, _hchi2xy, _hchi2zphi; + std::vector _hnsh, _hnshsig, _hnz; + // circle chi2/ndof of the helix's hits for the reco and the true circle (finder error model), + // and N(signal combo hits in the input collection) within 3 sigma of each circle + std::vector _hchi2creco, _hchi2ctrue; + std::vector _hnsig3reco, _hnsig3true; + int _nsig3true; + // signal combo hits: position, N(straw hits), bit mask of the helices containing it + std::vector _sx, _sy, _sz; + std::vector _snsh, _shel; + static float circleChi2(const ComboHit& hit, float xc, float yc, float r); + }; + + AprHelixTruthAna::AprHelixTruthAna(const Parameters& config) : + art::EDAnalyzer{config}, _hsTag(config().hsTag()), _chTag(config().chTag()), _sdmcTag(config().sdmcTag()), _minHits(config().minHits()), _paths(config().paths()), _pass(_paths.size(), 0) {} + + void AprHelixTruthAna::beginJob() { + art::ServiceHandle tfs; + _tree = tfs->make("helices", "helix seeds vs signal truth"); + _tree->Branch("run", &_run); + _tree->Branch("subrun", &_subrun); + _tree->Branch("event", &_event); + _tree->Branch("hasHelices", &_hasHelices); + _tree->Branch("pdg", &_pdg); + _tree->Branch("nsh", &_nsh); + _tree->Branch("p", &_p); + _tree->Branch("pt", &_pt); + _tree->Branch("pz", &_pz); + _tree->Branch("rtrue", &_rtrue); + _tree->Branch("lambdatrue", &_lambdatrue); + _tree->Branch("xctrue", &_xctrue); + _tree->Branch("yctrue", &_yctrue); + _tree->Branch("d0true", &_d0true); + _tree->Branch("nchsig", &_nchsig); + _tree->Branch("nzsig", &_nzsig); + _tree->Branch("zminsig", &_zminsig); + _tree->Branch("zmaxsig", &_zmaxsig); + _tree->Branch("gapsig", &_gapsig); + _tree->Branch("hlambda", &_hlambda); + _tree->Branch("hradius", &_hradius); + _tree->Branch("hxc", &_hxc); + _tree->Branch("hyc", &_hyc); + _tree->Branch("hd0", &_hd0); + _tree->Branch("hzmin", &_hzmin); + _tree->Branch("hzmax", &_hzmax); + _tree->Branch("hgap", &_hgap); + _tree->Branch("hchi2xy", &_hchi2xy); + _tree->Branch("hchi2zphi", &_hchi2zphi); + _tree->Branch("hnsh", &_hnsh); + _tree->Branch("hnshsig", &_hnshsig); + _tree->Branch("hnz", &_hnz); + _tree->Branch("hchi2creco", &_hchi2creco); + _tree->Branch("hchi2ctrue", &_hchi2ctrue); + _tree->Branch("hnsig3reco", &_hnsig3reco); + _tree->Branch("hnsig3true", &_hnsig3true); + _tree->Branch("nsig3true", &_nsig3true); + _tree->Branch("sx", &_sx); + _tree->Branch("sy", &_sy); + _tree->Branch("sz", &_sz); + _tree->Branch("snsh", &_snsh); + _tree->Branch("shel", &_shel); + for (size_t i = 0; i < _paths.size(); ++i) _tree->Branch(("trig_" + _paths[i]).c_str(), &_pass[i]); + } + + void AprHelixTruthAna::beginRun(const art::Run&) { + GeomHandle bfmgr; + GeomHandle det; + const double bz = bfmgr->getBField(det->toMu2e(CLHEP::Hep3Vector(0., 0., 0.))).z(); + _bzConv = bz * CLHEP::c_light / 1000.; // MeV/c per mm + } + + // squared circle residual in units of the AgnosticHelixFinder circle error + float AprHelixTruthAna::circleChi2(const ComboHit& hit, float xc, float yc, float r) { + const float dx = hit.pos().x() - xc, dy = hit.pos().y() - yc; + const float d2 = dx*dx + dy*dy; + const float dxn = dx * hit.vDir().x() + dy * hit.vDir().y(); + const float costh2 = dxn * dxn / d2; + const float err2 = hit.wireVar() * (1.f - costh2) + hit.transVar() * costh2; + const float res = r - std::sqrt(d2); + return res * res / err2; + } + + void AprHelixTruthAna::zInfo(std::vector z, float& zmin, float& zmax, float& gap, int& nz) { + zmin = zmax = gap = 0.; + nz = 0; + if (z.empty()) return; + std::sort(z.begin(), z.end()); + zmin = z.front(); + zmax = z.back(); + nz = 1; + for (size_t i = 1; i < z.size(); ++i) { + const float dz = z[i] - z[i-1]; + gap = std::max(gap, dz); + if (dz > 10.) ++nz; // a new plane face + } + } + + void AprHelixTruthAna::analyze(const art::Event& event) { + _run = event.run(); _subrun = event.subRun(); _event = event.event(); + if (!_paths.empty()) { + art::Handle trH; + event.getByLabel(art::InputTag("TriggerResults", "", moduleDescription().processName()), trH); + for (size_t i = 0; i < _paths.size(); ++i) { + _pass[i] = 0; + if (trH.isValid()) { + TriggerResultsNavigator nav(trH.product()); + if (nav.validPath(_paths[i])) _pass[i] = nav.accepted(_paths[i]); + } + } + } + const auto& sdmcs = *event.getValidHandle(_sdmcTag); + art::Handle chH; + event.getByLabel(_chTag, chH); + static const ComboHitCollection noHits; + const auto& chcol = chH.isValid() ? *chH : noHits; + // the helix finder may not have run (an earlier filter in its path failed) + art::Handle hsH; + event.getByLabel(_hsTag, hsH); + _hasHelices = hsH.isValid(); + static const HelixSeedCollection empty; + const auto& hscol = _hasHelices ? *hsH : empty; + + // signal particle: the highest-momentum e+- or mu+- (at its earliest straw hit) with >= minHits straw hits + // (the conversion electron in CE samples; pileup loopers have more hits but low momentum) + std::map, std::pair>> info; // N(hits), (earliest time, p) + for (const auto& mc : sdmcs) { + const auto& sgs = mc.earlyStrawGasStep(); + if (sgs.isNull()) continue; + const auto& simp = sgs->simParticle(); + if (std::abs(simp->pdgId()) != PDGCode::e_minus && std::abs(simp->pdgId()) != PDGCode::mu_minus) continue; // e+- and mu+- + auto& in = info[simp]; + if (in.first == 0 || sgs->time() < in.second.first) in.second = {float(sgs->time()), float(sgs->momvec().mag())}; + ++in.first; + } + art::Ptr sig; + _nsh = 0; + float pbest = 0.; + for (const auto& [simp, in] : info) { + if (in.first < _minHits || in.second.second <= pbest) continue; + pbest = in.second.second; + sig = simp; + _nsh = in.first; + } + + // truth at the earliest step of the signal particle + _pdg = 0; _p = _pt = _pz = _rtrue = _lambdatrue = _xctrue = _yctrue = _d0true = 0.; + float tmin = 1.e12; + for (const auto& mc : sdmcs) { + const auto& sgs = mc.earlyStrawGasStep(); + if (sgs.isNull() || sgs->simParticle() != sig || sgs->time() >= tmin) continue; + tmin = sgs->time(); + const auto mom = sgs->momvec(); + const auto pos = sgs->position(); + _pdg = sig->pdgId(); + _p = mom.mag(); _pt = mom.perp(); _pz = mom.z(); + const float q = (sig->pdgId() == PDGCode::e_minus || sig->pdgId() == PDGCode::mu_minus) ? -1. : 1.; // only e+- and mu+- reach here + _rtrue = _pt / std::fabs(_bzConv); + _lambdatrue = _pz / std::fabs(_bzConv); + // circle center: pos + R * q * sign(B) * (py, -px)/pT + const float s = q * (_bzConv > 0 ? 1. : -1.); + _xctrue = pos.x() + _rtrue * s * mom.y() / _pt; + _yctrue = pos.y() - _rtrue * s * mom.x() / _pt; + _d0true = std::sqrt(_xctrue*_xctrue + _yctrue*_yctrue) - _rtrue; + } + + // signal combo hits in the helix finder input collection: hits with a majority of signal straws + std::vector zsig; + std::vector sigHits; // signal combo hits in the input collection + std::vector shids; + for (size_t ich = 0; ich < chcol.size(); ++ich) { + shids.clear(); + chcol.fillStrawHitIndices(ich, shids); + int nsig = 0; + // the straw hit index is the StrawDigiMC index: true while the hit makers don't filter (FilterHits false) + for (auto sh : shids) { + const auto& sgs = sdmcs.at(sh).earlyStrawGasStep(); + if (sgs.isNonnull() && sgs->simParticle() == sig) ++nsig; + } + if (2 * nsig > (int)shids.size()) { zsig.push_back(chcol[ich].pos().z()); sigHits.push_back(ich); } + } + _nchsig = zsig.size(); + _sx.clear(); _sy.clear(); _sz.clear(); _snsh.clear(); _shel.clear(); + for (auto ich : sigHits) { + const auto& hit = chcol[ich]; + _sx.push_back(hit.pos().x()); _sy.push_back(hit.pos().y()); _sz.push_back(hit.pos().z()); + _snsh.push_back(hit.nStrawHits()); + int mask = 0; + for (size_t ih = 0; ih < hscol.size() && ih < 30; ++ih) { + for (const auto& hh : hscol[ih].hits()) { + if ((hh.pos() - hit.pos()).Mag2() < 1.e-6) { mask |= (1 << ih); break; } + } + } + _shel.push_back(mask); + } + _nsig3true = 0; + for (auto ich : sigHits) if (circleChi2(chcol[ich], _xctrue, _yctrue, _rtrue) < 9.f) ++_nsig3true; + zInfo(zsig, _zminsig, _zmaxsig, _gapsig, _nzsig); + + // helices + _hlambda.clear(); _hradius.clear(); _hxc.clear(); _hyc.clear(); _hd0.clear(); _hzmin.clear(); _hzmax.clear(); + _hgap.clear(); _hchi2xy.clear(); _hchi2zphi.clear(); _hnsh.clear(); _hnshsig.clear(); _hnz.clear(); + _hchi2creco.clear(); _hchi2ctrue.clear(); _hnsig3reco.clear(); _hnsig3true.clear(); + for (const auto& hs : hscol) { + const auto& h = hs.helix(); + _hlambda.push_back(h.lambda()); + _hradius.push_back(h.radius()); + _hxc.push_back(h.centerx()); + _hyc.push_back(h.centery()); + _hd0.push_back(h.rcent() - h.radius()); + _hchi2xy.push_back(h.chi2dXY()); + _hchi2zphi.push_back(h.chi2dZPhi()); + const auto& hits = hs.hits(); + std::vector z; + int nsh = 0, nshsig = 0; + for (size_t i = 0; i < hits.size(); ++i) { + z.push_back(hits[i].pos().z()); + shids.clear(); + hits.fillStrawHitIndices(i, shids); + for (auto sh : shids) { + ++nsh; + const auto& sgs = sdmcs.at(sh).earlyStrawGasStep(); + if (sgs.isNonnull() && sgs->simParticle() == sig) ++nshsig; + } + } + float zmin, zmax, gap; + int nz; + zInfo(z, zmin, zmax, gap, nz); + _hzmin.push_back(zmin); _hzmax.push_back(zmax); _hgap.push_back(gap); _hnz.push_back(nz); + _hnsh.push_back(nsh); _hnshsig.push_back(nshsig); + float c2r = 0., c2t = 0.; + for (size_t i = 0; i < hits.size(); ++i) { + c2r += circleChi2(hits[i], h.centerx(), h.centery(), h.radius()); + c2t += circleChi2(hits[i], _xctrue, _yctrue, _rtrue); + } + const float ndof = std::max(1.f, hits.size() - 3.f); + _hchi2creco.push_back(c2r / ndof); + _hchi2ctrue.push_back(c2t / ndof); + int n3r = 0, n3t = 0; + for (auto ich : sigHits) { + if (circleChi2(chcol[ich], h.centerx(), h.centery(), h.radius()) < 9.f) ++n3r; + if (circleChi2(chcol[ich], _xctrue, _yctrue, _rtrue) < 9.f) ++n3t; + } + _hnsig3reco.push_back(n3r); + _hnsig3true.push_back(n3t); + } + _tree->Fill(); + } + +} // namespace mu2e + +DEFINE_ART_MODULE(mu2e::AprHelixTruthAna) diff --git a/CalPatRec/fcl/prolog.fcl b/CalPatRec/fcl/prolog.fcl index 4c21c55218..851b984ce5 100644 --- a/CalPatRec/fcl/prolog.fcl +++ b/CalPatRec/fcl/prolog.fcl @@ -546,6 +546,12 @@ CalPatRec : { @table::CalPatRec caloClusterSigma : 30.0 minNHelixComboHits : 10 minNHelixStrawHits : 12 + maxHelixCandidates : 1 # save the first viable candidate + acceptHelixStrawHits : 0 # (maxHelixCandidates > 1) save a candidate with this many straw hits at once; 0: never + minAliasGap : 0. # (mm) 2pi phi-z alias test across the largest z gap; 0: off + targetRefitSigma : 0. # (mm) target-constrained refit of the best candidate; 0: off + targetRefitPhiWindow : 0.5 # (rad) target refit dphi/dz scan window + targetRefitDphiDzBins : 100 # target refit dphi/dz scan bins per helicity minHelixPerpMomentum : 20.0 maxHelixPerpMomentum : 140.0 minHelixMomentum : 30.0 diff --git a/CalPatRec/src/AgnosticHelixFinder_module.cc b/CalPatRec/src/AgnosticHelixFinder_module.cc index 38a6c5dbca..6b82bf3e91 100644 --- a/CalPatRec/src/AgnosticHelixFinder_module.cc +++ b/CalPatRec/src/AgnosticHelixFinder_module.cc @@ -94,6 +94,12 @@ namespace mu2e { fhicl::Atom caloClusterSigma {Name("caloClusterSigma" ), Comment("error assigned to calo clust") }; fhicl::Atom minNHelixStrawHits {Name("minNHelixStrawHits" ), Comment("straw hit save threshold" ) }; fhicl::Atom minNHelixComboHits {Name("minNHelixComboHits" ), Comment("combo hit save threshold" ) }; + fhicl::Atom maxHelixCandidates {Name("maxHelixCandidates" ), Comment("save the best (most straw hits) of up to this many candidates per search; 1: save the first")}; + fhicl::Atom acceptHelixStrawHits {Name("acceptHelixStrawHits" ), Comment("save a candidate with at least this many straw hits without searching further; <= 0: never")}; + fhicl::Atom minAliasGap {Name("minAliasGap" ), Comment("test the 2pi phi-z alias across the largest z gap of a helix's hits if it is at least this (mm); <= 0: off")}; + fhicl::Atom targetRefitSigma {Name("targetRefitSigma" ), Comment("if the best of the maxHelixCandidates (> 1) candidates has |d0| > this (mm), refit its circle through the stopping target (error = this) with a dphi/dz scan; the refit replaces it if it has at least as many straw hits; <= 0: off")}; + fhicl::Atom targetRefitPhiWindow {Name("targetRefitPhiWindow" ), Comment("phi window (rad) for hits to count on a line in the target refit's dphi/dz scan")}; + fhicl::Atom targetRefitDphiDzBins {Name("targetRefitDphiDzBins"), Comment("dphi/dz bins per helicity in the target refit's dphi/dz scan")}; fhicl::Atom minHelixPerpMomentum {Name("minHelixPerpMomentum" ), Comment("min pt of helix" ) }; fhicl::Atom maxHelixPerpMomentum {Name("maxHelixPerpMomentum" ), Comment("max pt of helix" ) }; fhicl::Atom minHelixMomentum {Name("minHelixMomentum" ), Comment("min momentum of helix" ) }; @@ -174,6 +180,12 @@ namespace mu2e { float _caloClusterSigmaSq; int _minNHelixStrawHits; int _minNHelixComboHits; + int _maxHelixCandidates; + int _acceptHelixStrawHits; + float _minAliasGap; + float _targetRefitSigma; + float _targetRefitPhiWindow; + int _targetRefitDphiDzBins; float _minHelixPerpMomentum; float _maxHelixPerpMomentum; float _minHelixMomentum; @@ -211,6 +223,17 @@ namespace mu2e { bool _intenseEvent; bool _intenseCluster; + // best helix candidate of the current search (maxHelixCandidates > 1) + struct helixCandidate { + std::vector tcHits; + ::LsqSums4 circleFitter; + ::LsqSums2 lineFitter; + int nStrawHits = 0, nComboHits = 0, nStrawHitsInTC = 0, nComboHitsInTC = 0; + float chi2 = 0.f; + float xc = 0.f, yc = 0.f, r = 0.f; // circle + }; + std::vector _candidates; + //----------------------------------------------------------------------------- // stuff for tool //----------------------------------------------------------------------------- @@ -264,8 +287,14 @@ namespace mu2e { bool refinePhiLine (size_t lineIndex); void initFinalSeed (LoopCondition& outcome); bool recoverPoints (); + bool resolveGapAlias (); + bool scanDphiDz (float window); + bool targetRefit (const helixCandidate& best, helixCandidate& refit); void checkHelixViability (LoopCondition& outcome); void saveHelix (size_t tc, HelixSeedCollection& HSColl); + bool saveAndContinue (size_t tc, HelixSeedCollection& HSColl, int nStrawHitsInHelix, int nComboHitsInHelix, + int nStrawHitsInTimeCluster, int nComboHitsInTimeCluster); + bool saveCandidates (size_t tc, HelixSeedCollection& HSColl); bool validHelixDirection (TrkFitDirection::FitDirection direction); bool goodTimeCluster (const TimeCluster& tc); @@ -317,6 +346,12 @@ namespace mu2e { _caloClusterSigma (config().caloClusterSigma() ), _minNHelixStrawHits (config().minNHelixStrawHits() ), _minNHelixComboHits (config().minNHelixComboHits() ), + _maxHelixCandidates (config().maxHelixCandidates() ), + _acceptHelixStrawHits (config().acceptHelixStrawHits() ), + _minAliasGap (config().minAliasGap() ), + _targetRefitSigma (config().targetRefitSigma() ), + _targetRefitPhiWindow (config().targetRefitPhiWindow() ), + _targetRefitDphiDzBins (config().targetRefitDphiDzBins() ), _minHelixPerpMomentum (config().minHelixPerpMomentum() ), _maxHelixPerpMomentum (config().maxHelixPerpMomentum() ), _minHelixMomentum (config().minHelixMomentum() ), @@ -824,13 +859,15 @@ namespace mu2e { // now we loop over triplets LoopCondition loopCondition; triplet tripletInfo; + _candidates.clear(); // viable helix candidates of this search (maxHelixCandidates > 1) for (size_t i = 0; i < _tcHits.size() - 2; i++) { if(_doTiming > 1) _watch->Increment("triplet-i"); if(_doTiming > 2) _watch->StopTime ("triplet-j"); if(shouldExitForTimeout(__func__)) { if(_doTiming > 0) _watch->StopTime(__func__); if(_doTiming > 1) _watch->StopTime("triplet-i"); - return false; + // keep the best candidate found before the timeout (maxHelixCandidates > 1) + return _candidates.empty() ? false : saveCandidates(tc, HSColl); } bool uselessSeed = true; setTripletI(i, tripletInfo, loopCondition); @@ -868,6 +905,7 @@ namespace mu2e { initFinalSeed(loopCondition); if (loopCondition == CONTINUE) { continue; } while (recoverPoints()); // while points are recovered, keep trying to recover + if (_minAliasGap > 0.f && resolveGapAlias()) { while (recoverPoints()); } checkHelixViability(loopCondition); if (loopCondition == CONTINUE) { continue; } // before saving helix we make sure it has enough hits @@ -887,15 +925,30 @@ namespace mu2e { nStrawHitsInHelix += nStrawHitsInHit; ++nComboHitsInHelix; } - if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits) { - saveHelix(tc, HSColl); - if (_diagLevel > 0) { _diagInfo.nHelices++; } - // we only want to search for another helix if we have enough remaining hits after saving helix - const int remainingStrawHits = nStrawHitsInTimeCluster - nStrawHitsInHelix; - const int remainingComboHits = nComboHitsInTimeCluster - nComboHitsInHelix; - const bool findAnotherHelix = remainingStrawHits >= _minNHelixStrawHits && remainingComboHits >= _minNHelixComboHits; - if(findAnotherHelix && _doAverageFlag) resetFlags(); // reset flags if needed - + if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits && + _maxHelixCandidates > 1 && !(_acceptHelixStrawHits > 0 && nStrawHitsInHelix >= _acceptHelixStrawHits)) { + // keep the best candidate (most straw hits, then the best fit) and continue the search + helixCandidate c; + c.tcHits = _tcHits; + c.circleFitter = _circleFitter; + c.lineFitter = _lineFitter; + c.nStrawHits = nStrawHitsInHelix; + c.nComboHits = nComboHitsInHelix; + c.nStrawHitsInTC = nStrawHitsInTimeCluster; + c.nComboHitsInTC = nComboHitsInTimeCluster; + c.chi2 = _circleFitter.chi2DofCircle() + _lineFitter.chi2Dof(); + c.xc = _circleFitter.x0(); c.yc = _circleFitter.y0(); c.r = _circleFitter.radius(); + _candidates.push_back(std::move(c)); + if ((int)_candidates.size() < _maxHelixCandidates) continue; + // enough candidates: save the best one + if(_doTiming > 0) _watch->StopTime(__func__); + if(_doTiming > 1) _watch->StopTime("triplet-i"); + if(_doTiming > 2) _watch->StopTime("triplet-j"); + if(_doTiming > 3) _watch->StopTime("triplet-k"); + return saveCandidates(tc, HSColl); + } else if (nStrawHitsInHelix >= _minNHelixStrawHits && nComboHitsInHelix >= _minNHelixComboHits) { + const bool findAnotherHelix = saveAndContinue(tc, HSColl, nStrawHitsInHelix, nComboHitsInHelix, + nStrawHitsInTimeCluster, nComboHitsInTimeCluster); // Helix found, so we can exit this search if(_doTiming > 0) _watch->StopTime(__func__); if(_doTiming > 1) _watch->StopTime("triplet-i"); @@ -915,10 +968,103 @@ namespace mu2e { if(_doTiming > 2) _watch->StopTime("triplet-j"); if(_doTiming > 3) _watch->StopTime("triplet-k"); + // search finished with fewer than maxHelixCandidates candidates: save the best one + if (!_candidates.empty()) return saveCandidates(tc, HSColl); + // No helix found, no need to continue searching in this time cluster return false; } + //----------------------------------------------------------------------------- + // target-constrained refit of a candidate with |d0| > targetRefitSigma: the circle is refit + // with its hits plus the stopping target (0,0) with error targetRefitSigma, then the phi-z line + // is found by a dphi/dz scan on that circle, the hits recovered and the helix checked as usual. + // Short arcs leave the radius and the 2pi phi-z count ambiguous (half radius / half lambda + // solutions fit the hits too); this tries the solution through the target. + //----------------------------------------------------------------------------- + bool AgnosticHelixFinder::targetRefit(const helixCandidate& best, helixCandidate& refit) { + const float d0 = std::sqrt(best.xc * best.xc + best.yc * best.yc) - best.r; + if (std::fabs(d0) <= _targetRefitSigma) return false; + _tcHits = best.tcHits; + const float wT = 1.f / (_targetRefitSigma * _targetRefitSigma); + // circle: the candidate's hits plus the target, errors updated once + for (int iter = 0; iter < 2; ++iter) { + const float xC = iter == 0 ? best.xc : _circleFitter.x0(); + const float yC = iter == 0 ? best.yc : _circleFitter.y0(); + _circleFitter.clear(); + _circleFitter.addPoint(0.f, 0.f, wT); + for (size_t i = 0; i < _tcHits.size(); i++) { + const auto& h = _tcHits[i]; + if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET) continue; // a used calo cluster is kept, as in the other fits + computeCircleError2(i, xC, yC); + _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / _tcHits[i].circleError2); + } + } + const float pt = computeHelixPerpMomentum(_circleFitter.radius()); + if (pt < _minHelixPerpMomentum || pt > _maxHelixPerpMomentum) return false; + if (!scanDphiDz(_targetRefitPhiWindow)) return false; + while (recoverPoints()); + LoopCondition outcome; + checkHelixViability(outcome); + if (outcome == CONTINUE) return false; + int nsh = 0, nch = 0; + for (const auto& h : _tcHits) { + if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue; + nsh += h.hit->nStrawHits(); + ++nch; + } + if (nsh < _minNHelixStrawHits || nch < _minNHelixComboHits) return false; + refit.tcHits = _tcHits; + refit.circleFitter = _circleFitter; + refit.lineFitter = _lineFitter; + refit.nStrawHits = nsh; + refit.nComboHits = nch; + refit.nStrawHitsInTC = best.nStrawHitsInTC; + refit.nComboHitsInTC = best.nComboHitsInTC; + refit.chi2 = _circleFitter.chi2DofCircle() + _lineFitter.chi2Dof(); + refit.xc = _circleFitter.x0(); refit.yc = _circleFitter.y0(); refit.r = _circleFitter.radius(); + return true; + } + + //----------------------------------------------------------------------------- + // save the best candidate (most straw hits, then the best fit), or its target-constrained + // refit (targetRefitSigma > 0) if that has at least as many straw hits. Returns whether to + // search for another helix. + //----------------------------------------------------------------------------- + bool AgnosticHelixFinder::saveCandidates(size_t tc, HelixSeedCollection& HSColl) { + size_t ib = 0; + for (size_t i = 1; i < _candidates.size(); ++i) { + const auto& c = _candidates[i]; + const auto& b = _candidates[ib]; + if (c.nStrawHits > b.nStrawHits || (c.nStrawHits == b.nStrawHits && c.chi2 < b.chi2)) ib = i; + } + if (_targetRefitSigma > 0.f) { + helixCandidate refit; + if (targetRefit(_candidates[ib], refit) && refit.nStrawHits >= _candidates[ib].nStrawHits) + _candidates[ib] = std::move(refit); + } + const auto& best = _candidates[ib]; + _tcHits = best.tcHits; + _circleFitter = best.circleFitter; + _lineFitter = best.lineFitter; + return saveAndContinue(tc, HSColl, best.nStrawHits, best.nComboHits, best.nStrawHitsInTC, best.nComboHitsInTC); + } + + //----------------------------------------------------------------------------- + // save the current helix (fitters and used hits) and decide whether to search for another one + //----------------------------------------------------------------------------- + bool AgnosticHelixFinder::saveAndContinue(size_t tc, HelixSeedCollection& HSColl, int nStrawHitsInHelix, int nComboHitsInHelix, + int nStrawHitsInTimeCluster, int nComboHitsInTimeCluster) { + saveHelix(tc, HSColl); + if (_diagLevel > 0) { _diagInfo.nHelices++; } + // we only want to search for another helix if we have enough remaining hits after saving helix + const int remainingStrawHits = nStrawHitsInTimeCluster - nStrawHitsInHelix; + const int remainingComboHits = nComboHitsInTimeCluster - nComboHitsInHelix; + const bool findAnotherHelix = remainingStrawHits >= _minNHelixStrawHits && remainingComboHits >= _minNHelixComboHits; + if(findAnotherHelix && _doAverageFlag) resetFlags(); // reset flags if needed + return findAnotherHelix; + } + //----------------------------------------------------------------------------- // check flags to see if point is good for triplet-ing with //----------------------------------------------------------------------------- @@ -1567,6 +1713,158 @@ namespace mu2e { return recoveries; } + //----------------------------------------------------------------------------- + // test the 2pi alias of the phi-z line across the largest z gap between the used hits: the + // phi of the hits on one side of the gap is shifted by +-2pi and the line refit. The + // alternative with the lowest chi2/dof, if lower than the current line's and within the + // momentum range, replaces it. Returns true if the line changed. + //----------------------------------------------------------------------------- + bool AgnosticHelixFinder::resolveGapAlias() { + if(_doTiming > 4) _watch->SetTime(__func__); + // used tracker hits, in decreasing z (the _tcHits order); a used calo cluster is in the line + // fits, as in recoverPoints, on the high-z side of any tracker gap, but not in the gap search + std::vector idx, calo; + for (size_t i = 0; i < _tcHits.size(); i++) { + const auto& h = _tcHits[i]; + if (!h.used || h.inHelix || h.hitIndice == HitType::STOPPINGTARGET) continue; + if (h.hitIndice == HitType::CALOCLUSTER) calo.push_back(i); + else idx.push_back(i); + } + if (idx.size() < 4) { if(_doTiming > 4) _watch->StopTime(__func__); return false; } + size_t igap = 0; // hits idx[0..igap] are on the high-z side + float gap = 0.f; + for (size_t k = 0; k + 1 < idx.size(); ++k) { + const float dz = getPos(idx[k]).z() - getPos(idx[k+1]).z(); + if (dz > gap) { gap = dz; igap = k; } + } + if (gap < _minAliasGap) { if(_doTiming > 4) _watch->StopTime(__func__); return false; } + + // the current line (shift 0) is refit over the same points as the alternatives + const float radius = _circleFitter.radius(); + float bestChi2 = 0.f; + int bestShift = 0; + ::LsqSums2 bestFitter; + for (int shift : {0, -1, 1}) { + ::LsqSums2 fitter; + for (size_t k = 0; k < idx.size(); ++k) { + const auto& h = _tcHits[idx[k]]; + const int n = h.helixPhiCorrection + (k > igap ? shift : 0); + fitter.addPoint(getPos(idx[k]).z(), h.helixPhi + _twopi * n, 1.f / h.helixPhiError2); + } + for (auto i : calo) { + const auto& h = _tcHits[i]; + fitter.addPoint(getPos(i).z(), h.helixPhi + _twopi * h.helixPhiCorrection, 1.f / h.helixPhiError2); + } + if (shift == 0) { bestChi2 = fitter.chi2Dof(); continue; } + const float mom2 = computeHelixMomentum2(radius, fitter.dydx()); + if (mom2 < _minHelixMomentumSq || mom2 > _maxHelixMomentumSq) continue; + if (fitter.chi2Dof() < bestChi2) { + bestChi2 = fitter.chi2Dof(); + bestShift = shift; + bestFitter = fitter; + } + } + if (bestShift != 0) { + for (size_t k = igap + 1; k < idx.size(); ++k) _tcHits[idx[k]].helixPhiCorrection += bestShift; + _lineFitter = bestFitter; + } + if(_doTiming > 4) _watch->StopTime(__func__); + return bestShift != 0; + } + + //----------------------------------------------------------------------------- + // scan the phi-z slope on the current circle (target refit): for each dphi/dz allowed by the + // momentum range (both helicities, targetRefitDphiDzBins bins each), count the time-cluster + // hits on the circle (within maxCircleRecoverSigma) whose phi lies on a common line, modulo + // 2pi, within window. The hits on the best line replace the current ones and the line and + // circle are refit. Returns false if no line with >= 3 hits is found. + //----------------------------------------------------------------------------- + bool AgnosticHelixFinder::scanDphiDz(float window) { + if(_doTiming > 4) _watch->SetTime(__func__); + const float xC = _circleFitter.x0(), yC = _circleFitter.y0(), rC = _circleFitter.radius(); + const float maxRes2 = _maxCircleRecoverSigma * _maxCircleRecoverSigma; + std::vector cand; + for (size_t i = 0; i < _tcHits.size(); i++) { + auto& h = _tcHits[i]; + if (h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue; + computeCircleError2(i, xC, yC); + if (computeCircleResidual2(i, xC, yC, rC) > maxRes2) continue; + cand.push_back(i); + } + auto done = [&](bool res) { if(_doTiming > 4) _watch->StopTime(__func__); return res; }; + if (cand.size() < 3) return done(false); + std::vector z(cand.size()), phi(cand.size()); + for (size_t k = 0; k < cand.size(); ++k) { + computeHelixPhi(cand[k], xC, yC); + z[k] = getPos(cand[k]).z(); + phi[k] = _tcHits[cand[k]].helixPhi; + } + // slope range from the momentum limits: p^2 = (R^2 + lambda^2) (cB)^2, |dphi/dz| <= maxDphiDz + const float lam2Max = _maxHelixMomentumSq / _bz0ConvSq - rC * rC; + const float lam2Min = _minHelixMomentumSq / _bz0ConvSq - rC * rC; + const float sLow = (lam2Max > 0.f) ? 1.f / std::sqrt(lam2Max) : 0.f; + float sHigh = _maxDphiDz; + if (lam2Min > 0.f) sHigh = std::min(sHigh, 1.f / std::sqrt(lam2Min)); + if (sHigh <= sLow) return done(false); + + int bestN = 0; + float bestS = 0.f, bestPhi0 = 0.f; + std::vector ph(cand.size()), ph2(2 * cand.size()); + const int nb = std::max(1, _targetRefitDphiDzBins); + for (int sign : {-1, 1}) { + for (int b = 0; b <= nb; ++b) { + const float slope = sign * (sLow + (sHigh - sLow) * b / nb); + for (size_t k = 0; k < cand.size(); ++k) { + float r = std::fmod(phi[k] - slope * z[k], float(_twopi)); + if (r < 0.f) r += _twopi; + ph[k] = r; + } + std::sort(ph.begin(), ph.end()); + for (size_t k = 0; k < ph.size(); ++k) { ph2[k] = ph[k]; ph2[k + ph.size()] = ph[k] + _twopi; } + // largest number of phases within a window (circular) + size_t j = 0; + for (size_t k = 0; k < ph.size(); ++k) { + if (j < k) j = k; + while (j + 1 < ph2.size() && ph2[j + 1] - ph2[k] <= window) ++j; + const int n = j - k + 1; + if (n > bestN) { bestN = n; bestS = slope; bestPhi0 = 0.5f * (ph2[k] + ph2[j]); } + } + } + } + if (bestN < 3) return done(false); + + // adopt the best line: hits within the window are used, with their 2pi corrections + for (size_t i = 0; i < _tcHits.size(); i++) { + auto& h = _tcHits[i]; + if (h.inHelix || h.hitIndice == HitType::STOPPINGTARGET || h.hitIndice == HitType::CALOCLUSTER) continue; + h.used = false; + } + _lineFitter.clear(); + _circleFitter.clear(); + for (size_t k = 0; k < cand.size(); ++k) { + auto& h = _tcHits[cand[k]]; + const float pred = bestS * z[k] + bestPhi0; + const int n = std::round((pred - phi[k]) / _twopi); + if (std::fabs(phi[k] + _twopi * n - pred) > 0.5f * window + 1.e-4f) continue; + h.used = true; + h.helixPhiCorrection = n; + _lineFitter.addPoint(z[k], phi[k] + _twopi * n, 1.f / h.helixPhiError2); + _circleFitter.addPoint(getPos(cand[k]).x(), getPos(cand[k]).y(), 1.f / h.circleError2); + } + // keep a used calo cluster in the circle and line fits, as in the other fits + for (size_t i = 0; i < _tcHits.size(); i++) { + auto& h = _tcHits[i]; + if (h.hitIndice != HitType::CALOCLUSTER || !h.used || h.inHelix) continue; + computeCircleError2(i, xC, yC); + computeHelixPhi(i, xC, yC); + const float zC = getPos(i).z(); + h.helixPhiCorrection = std::round((bestS * zC + bestPhi0 - h.helixPhi) / _twopi); + _circleFitter.addPoint(getPos(i).x(), getPos(i).y(), 1.f / h.circleError2); + _lineFitter.addPoint(zC, h.helixPhi + _twopi * h.helixPhiCorrection, 1.f / h.helixPhiError2); + } + return done(true); + } + //----------------------------------------------------------------------------- // check if helix is savable based on fcl defined cuts //----------------------------------------------------------------------------- diff --git a/Mu2eUtilities/inc/HelixTool.hh b/Mu2eUtilities/inc/HelixTool.hh index 7c4c19e567..370c5e78ca 100644 --- a/Mu2eUtilities/inc/HelixTool.hh +++ b/Mu2eUtilities/inc/HelixTool.hh @@ -4,6 +4,8 @@ // Original author G. Pezzullo // +#include + #include "Offline/DataProducts/inc/StrawId.hh" #include "Math/GenVector/VectorUtil.h" #include "Math/VectorUtil.h" @@ -13,6 +15,7 @@ #include "Offline/TrackerGeom/inc/Tracker.hh" namespace mu2e { class Tracker; } +namespace mu2e { class Plane; } namespace mu2e { struct HelixSeed; } using namespace ROOT::Math::VectorUtil; @@ -22,35 +25,44 @@ namespace mu2e { class HelixTool{ public: - HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker); + // returns true for a straw that can't produce a hit (e.g. TrackerStatus noSignal or suppress). + // An empty function means all straws are live. A function is used instead of TrackerStatus + // to avoid a library dependency of Mu2eUtilities on TrackerConditions. + using StrawNoSignal = std::function; + + HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker, StrawNoSignal NoSignal = StrawNoSignal()); // Accept compiler supplied d'tor, copy c'tor and assignment operator. // Accessors float nLoops () const { return _nLoops; } - int nMinHitsLoop () const { return _nMinHitsLoop; } float meanHitRadialDist() const { return _meanHitRadialDist; } float d0 () const { return _d0; } - float nstrawhits () const { return _nStrawHits; } + int nstrawhits () const { return _nStrawHits; } - //function that evaluates the ratio between the measured ComboHits and the - //expected intesections of the helix with the tracker planes. This function - //models the tracker as a perfect cylinder - float hitRatio () const { return _hitRatio; } + //function that evaluates the ratio between the measured straw hits and the + //expected intersections of the helix with the tracker plane faces inside the tracker + //active radius. If a StrawNoSignal function is given, each face crossing is weighted by + //the fraction of live straws (nearest straw in each layer) at the crossing point, so that + //dead panels and straws don't lower the ratio. A crossing with no straw of the face at that + //point (outside the panel's straws) keeps the weight 1. Evaluated on the first call. + float hitRatio () const; void dirOfProp(float& slope, float& slopeErr, float& chi2ndof); private: + float liveFraction(const Plane& Pln, uint16_t Face, float X, float Y) const; + const HelixSeed* _hel; - int _nMinHitsLoop; float _nLoops; int _nStrawHits; float _meanHitRadialDist; float _d0; - float _hitRatio; - float _trackerRIn, _trackerROut, _trackerLength; + mutable float _hitRatio = -1.; // < 0: not evaluated yet + float _trackerRIn, _trackerROut; const mu2e::Tracker*_tracker; + StrawNoSignal _noSignal; }; } // namespace mu2e diff --git a/Mu2eUtilities/src/HelixTool.cc b/Mu2eUtilities/src/HelixTool.cc index dcfb88070c..afb9eb712b 100644 --- a/Mu2eUtilities/src/HelixTool.cc +++ b/Mu2eUtilities/src/HelixTool.cc @@ -1,6 +1,8 @@ #include #include +#include #include +#include #include #include "CLHEP/Vector/ThreeVector.h" @@ -22,9 +24,10 @@ namespace mu2e { - HelixTool::HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker) { + HelixTool::HelixTool(const HelixSeed *Helix, const mu2e::Tracker*MyTracker, StrawNoSignal NoSignal) { _hel = Helix; _tracker = MyTracker; + _noSignal = std::move(NoSignal); _trackerRIn = _tracker->g4Tracker()->getInnerTrackerEnvelopeParams().innerRadius(); _trackerROut = _tracker->g4Tracker()->getInnerTrackerEnvelopeParams().outerRadius(); @@ -66,38 +69,106 @@ namespace mu2e { //here we evaluate once the impact parameter _d0 = robustHel->rcent () - robustHel->radius (); + } + + float HelixTool::hitRatio() const { + if (_hitRatio >= 0.) return _hitRatio; - // we now estiamte the ratio of the number of measured hits to the number of the expected ones - // we make a few assumptions and appriximations: + // ratio of the number of measured hits to the number of expected ones: count the plane + // faces the helix crosses inside the tracker active radius, weighted by the fraction of + // live straws at the crossing + const mu2e::RobustHelix& robustHel = _hel->helix(); float expected_faces(1e-6); for (size_t planeId=0; planeId<_tracker->nPlanes(); planeId++) { - const Plane* pln = &_tracker->getPlane(planeId); - int nPanels = pln->nPanels(); - if (nPanels == 0 ) continue; - std::array idPanels = {0, (nPanels-1)}; + const Plane& pln = _tracker->getPlane(planeId); + if (pln.nPanels() == 0 ) continue; - for (size_t ipn=0; ipngetPanel(ipn); - float z = (panel->getStraw(0).getMidPoint().z()+panel->getStraw(1).getMidPoint().z())/2.; - // if (z < z_first_hit ) continue; - // if (z > z_last_hit ) continue; + for (uint16_t face=0; faceposition(pos); + robustHel.position(pos); //now check that we are in the active area of tracker float hitR = sqrtf(pos.x()*pos.x() + pos.y()*pos.y()); if (hitR < _trackerRIn ) continue; if (hitR > _trackerROut) continue; - ++expected_faces; - }//end loop over the two panels + expected_faces += _noSignal ? liveFraction(pln, face, pos.x(), pos.y()) : 1.; + }//end loop over the two faces }//end loop over the planes - //each plane has two layers of panels. We are assuming 2 ComboHits per face _hitRatio = _nStrawHits/expected_faces; + return _hitRatio; + } + + + float HelixTool::liveFraction(const Plane& Pln, uint16_t Face, float X, float Y) const { + // find the panel of this face covering the crossing point: the one whose outward direction, + // perpendicular to its straws in the transverse plane, is closest to the crossing direction + const Panel* best(nullptr); + float bestRho(0), nx(0), ny(0); + for (int ipn=Face; ipn bestRho) { + best = &panel; + bestRho = rho; + nx = px; + ny = py; + } + } + if (best == nullptr) return 1.; + + // straws are ordered radially (straw 0 innermost), alternating between the two layers: + // estimate the straw at the crossing and take the nearest straw of each layer around it + auto strawRho = [&](int s) { + const auto& mid = best->getStraw(s).getMidPoint(); + return float(mid.x()*nx + mid.y()*ny); + }; + const int nStraws = best->nStraws(); + float rho0 = strawRho(0); + float rhoMax = strawRho(nStraws-1); + // no straw of this face at the crossing (radially outside the panel, or beyond the straw + // ends, e.g. between the panels of the face): keep the weight without the status + if (bestRho < rho0 || bestRho > rhoMax) return 1.; + int sEst = (rhoMax > rho0) ? std::lround((bestRho-rho0)/(rhoMax-rho0)*(nStraws-1)) : 0; + sEst = std::clamp(sEst, 0, nStraws-1); + { + const Straw& straw = best->getStraw(sEst); + const auto& mid = straw.getMidPoint(); + const auto& dir = straw.getDirection(); + if (std::fabs((X-mid.x())*dir.x() + (Y-mid.y())*dir.y()) > straw.halfLength()) return 1.; + } + + float nLive(0); + for (int layer=0; layergetStraw(s).id().layer() != layer) continue; + float d = std::fabs(strawRho(s) - bestRho); + if (nearest < 0 || d < dmin) { + nearest = s; + dmin = d; + } + } + if (nearest >= 0 && !_noSignal(best->getStraw(nearest).id())) ++nLive; + } + return nLive/StrawId::_nlayers; } diff --git a/Trigger/CMakeLists.txt b/Trigger/CMakeLists.txt index 1c8df531b2..28ba03d424 100644 --- a/Trigger/CMakeLists.txt +++ b/Trigger/CMakeLists.txt @@ -32,7 +32,9 @@ cet_build_plugin(ReadTriggerInfo art::module Offline::GlobalConstantsService Offline::MCDataProducts Offline::Mu2eUtilities + Offline::ProditionsService Offline::RecoDataProducts + Offline::TrackerConditions Offline::TrackerGeom KinKal::Trajectory ) diff --git a/Trigger/src/ReadTriggerInfo_module.cc b/Trigger/src/ReadTriggerInfo_module.cc index 6b61db3c49..9b673a47d3 100644 --- a/Trigger/src/ReadTriggerInfo_module.cc +++ b/Trigger/src/ReadTriggerInfo_module.cc @@ -32,6 +32,8 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/TrackerConditions/inc/TrackerStatus.hh" //Data products #include "Offline/RecoDataProducts/inc/CaloCluster.hh" @@ -192,6 +194,7 @@ namespace mu2e { fhicl::Atom genCountTag {Name("genCount" ), Comment("GenEventCount label" ), "genCounter" }; fhicl::Atom PBITag {Name("PBITag" ), Comment("ProtonBunchIntensity label" ), "PBISim" }; fhicl::Atom processName {Name("processName" ), Comment("globalTrigger" ), "" }; + fhicl::Atom useTrackerStatus {Name("useTrackerStatus" ), Comment("don't count straws without signal (TrackerStatus) as expected hits in the helix hit ratio; set as in the HelixFilter"), true }; }; explicit ReadTriggerInfo(const art::EDAnalyzer::Table& config); @@ -260,6 +263,9 @@ namespace mu2e { occupancyHist_ _occupancyHist; const mu2e::Tracker* _tracker = nullptr; + ProditionsHandle _trackerStatus_h; + bool _useTrackerStatus = true; + HelixTool::StrawNoSignal _noSignal; // straws not counted as expected hits in the helix hit ratio; empty: all straws live float _minPOT, _maxPOT; const int _genOccIndex; @@ -284,6 +290,7 @@ namespace mu2e { _maxPOT (3.5e8 ), _genOccIndex (2*_nTrackTrig + _nCaloTrig ) { + _useTrackerStatus = config().useTrackerStatus(); _trigAll. resize(_nMaxTrig); _trigFinal. resize(_nMaxTrig); @@ -765,6 +772,11 @@ namespace mu2e { //--------------------------------------------------------------------------------// void ReadTriggerInfo::analyze(const art::Event& event) { + // same definition of the helix hit ratio as in HelixFilter + if (_useTrackerStatus) { + const TrackerStatus* status = &_trackerStatus_h.get(event.id()); + _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); }; + } if(!_useNGen) _nProcess += 1; @@ -1011,7 +1023,7 @@ namespace mu2e { //--------------------------------------------------------------------------------// void ReadTriggerInfo::fillHelixTrigInfo(int HelTrigIndex, const HelixSeed* HSeed, helixInfoHist_& Hist) { if(!Hist._hHelInfo[HelTrigIndex][0]) throw cet::exception("BADCONFIG") << __func__ << ": Helix histogram index " << HelTrigIndex << " not initialized\n"; - HelixTool helTool(HSeed, _tracker); + HelixTool helTool(HSeed, _tracker, _noSignal); int nch = (int)HSeed->hits().size(); int nsh(0); diff --git a/Trigger/src/SConscript b/Trigger/src/SConscript index 65b2a77290..f197921806 100644 --- a/Trigger/src/SConscript +++ b/Trigger/src/SConscript @@ -23,6 +23,7 @@ helper.make_plugins( [ mainlib, 'mu2e_MCDataProducts', 'mu2e_RecoDataProducts', 'mu2e_TrackerConditions', + 'mu2e_DbTables', 'mu2e_GeometryService', 'mu2e_BFieldGeom', 'mu2e_CalorimeterGeom', diff --git a/TrkFilters/CMakeLists.txt b/TrkFilters/CMakeLists.txt index d7dbca9c1c..90886f149c 100644 --- a/TrkFilters/CMakeLists.txt +++ b/TrkFilters/CMakeLists.txt @@ -6,7 +6,9 @@ cet_build_plugin(HelixFilter art::module Offline::GeneralUtilities Offline::GeometryService Offline::Mu2eUtilities + Offline::ProditionsService Offline::RecoDataProducts + Offline::TrackerConditions Offline::TrackerGeom ) @@ -24,7 +26,9 @@ cet_build_plugin(MultiHelixFilter art::module Offline::GeneralUtilities Offline::GeometryService Offline::Mu2eUtilities + Offline::ProditionsService Offline::RecoDataProducts + Offline::TrackerConditions Offline::TrackerGeom ) diff --git a/TrkFilters/src/HelixFilter_module.cc b/TrkFilters/src/HelixFilter_module.cc index 4da6f163dd..6606338814 100644 --- a/TrkFilters/src/HelixFilter_module.cc +++ b/TrkFilters/src/HelixFilter_module.cc @@ -17,6 +17,8 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/TrackerConditions/inc/TrackerStatus.hh" // data #include "Offline/RecoDataProducts/inc/HelixSeed.hh" #include "Offline/DataProducts/inc/Helicity.hh" @@ -104,6 +106,8 @@ namespace mu2e HelixCutsTool():_configured(false){}; void setTrackerGeomHandle(const Tracker* TrackerGeom) { _myTracker = TrackerGeom; } + // straws without signal are not counted as expected hits in the hit ratio; nullptr: all straws live + void setTrackerStatus(const TrackerStatus* Status) { _trackerStatus = Status; } int evalIPAPresc(const float &phi0){ int val= (_prescalerPar._amplitude - (_prescalerPar._amplitude-1)*sin(_prescalerPar._frequency*phi0 + _prescalerPar._phase)); @@ -116,7 +120,11 @@ namespace mu2e //check the helicity if (_doHelicityCheck && !(Helix.helix().helicity() == Helicity(_hel))) return false; float mm2MeV = 3./10.;//FIXME! - HelixTool helTool(&Helix, _myTracker); + HelixTool::StrawNoSignal noSignal; + if (_trackerStatus) { + noSignal = [status = _trackerStatus](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); }; + } + HelixTool helTool(&Helix, _myTracker, noSignal); // compute the helix momentum. Note this is in units of mm!!! float hmom = Helix.helix().momentum()*mm2MeV; int nstrawhits = helTool.nstrawhits(); @@ -126,7 +134,6 @@ namespace mu2e float d0 = Helix.helix().rcent() - Helix.helix().radius(); float lambda = std::fabs(Helix.helix().lambda()); float nLoops = helTool.nLoops(); - float hRatio = helTool.hitRatio(); const float slope = Helix.recoDir().slope(); const float slopeErr = std::fabs(Helix.recoDir().slopeErr()); const float slopeSignedSig = (slopeErr > 0.f) ? slope/slopeErr : 0.f; //signifance from 0 signed by the slope direction @@ -135,7 +142,7 @@ namespace mu2e std::cout << "[HelixFilter] : status = " << Helix.status() << " nhits = " << nstrawhits << " mom = " << hmom << std::endl << "[HelixFilter] : chi2XY = " << chi2XY << " chi2ZPHI = " << chi2PhiZ << " d0 = " << d0 << " lambda = "<< lambda << " p = " << hmom << " pT = " << hpT - << " nLoops = " << nLoops << " hRatio = "<< hRatio << " slopeSignedSig = " << slopeSignedSig << std::endl; + << " nLoops = " << nLoops << " hRatio = "<< helTool.hitRatio() << " slopeSignedSig = " << slopeSignedSig << std::endl; } if( Helix.status().hasAllProperties(_goodh) && (!_hascc || Helix.caloCluster().isNonnull()) && @@ -153,7 +160,7 @@ namespace mu2e hmom <= _maxmom && (!_useSlopeSigMin || slopeSignedSig > _slopeSigMin) && (!_useSlopeSigMax || slopeSignedSig < _slopeSigMax) && - hRatio >= _minHitRatio ) { + (_minHitRatio <= 0. || helTool.hitRatio() >= _minHitRatio) ) { // hit ratio last: it is the most expensive //now check if we want to prescake or not if (_prescaleUsingD0Phi) { float phiAtD0 = Helix.helix().fcent(); @@ -193,6 +200,7 @@ namespace mu2e bool _prescaleUsingD0Phi = false; PhiPrescalingParams _prescalerPar; const Tracker*_myTracker = nullptr; + const TrackerStatus* _trackerStatus = nullptr; }; struct Config{ @@ -205,6 +213,7 @@ namespace mu2e fhicl::Atom noFilter { Name("noFilter"), Comment("don't apply the filter decision"), false}; fhicl::Atom minNHelices { Name("minNHelices"), Comment("minimum number of helices passing the cuts"), 1}; fhicl::Atom maxDt0 { Name("maxDt0"), Comment("maximum difference in time between helices"), -1}; + fhicl::Atom useTrackerStatus { Name("useTrackerStatus"), Comment("don't count straws without signal (TrackerStatus) as expected hits in the hit ratio"), true}; }; @@ -229,6 +238,8 @@ namespace mu2e bool _noFilter; unsigned _minNHelices; float _maxDt0; + bool _useTrackerStatus; + ProditionsHandle _trackerStatus_h; bool checkHelixFromHelicity(const HelixSeed&helix); }; @@ -243,7 +254,8 @@ namespace mu2e _npass(0), _noFilter (config().noFilter()), _minNHelices (config().minNHelices()), - _maxDt0 (config().maxDt0()) + _maxDt0 (config().maxDt0()), + _useTrackerStatus (config().useTrackerStatus()) { produces(); if(_minNHelices < 2 && _maxDt0 >= 0.) throw cet::exception("BADCONFIG") << "Requested a timing difference cut of " << _maxDt0 << " ns between helices but only " << _minNHelices << " helices required"; @@ -291,6 +303,11 @@ namespace mu2e // create output std::unique_ptr triginfo(new TriggerInfo); ++_nevt; + if (_useTrackerStatus) { + const TrackerStatus& trackerStatus = _trackerStatus_h.get(evt.id()); + _posHelCuts.setTrackerStatus(&trackerStatus); + _negHelCuts.setTrackerStatus(&trackerStatus); + } // find the collection auto hsH = evt.getValidHandle(_hsTag); const HelixSeedCollection* hscol = hsH.product(); diff --git a/TrkFilters/src/MultiHelixFilter_module.cc b/TrkFilters/src/MultiHelixFilter_module.cc index d313cc0cdb..11ba3ab223 100644 --- a/TrkFilters/src/MultiHelixFilter_module.cc +++ b/TrkFilters/src/MultiHelixFilter_module.cc @@ -17,6 +17,8 @@ #include "Offline/GeometryService/inc/GeomHandle.hh" #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/TrackerGeom/inc/Tracker.hh" +#include "Offline/ProditionsService/inc/ProditionsHandle.hh" +#include "Offline/TrackerConditions/inc/TrackerStatus.hh" // data #include "Offline/RecoDataProducts/inc/CaloCluster.hh" @@ -104,6 +106,7 @@ namespace mu2e } void setB(const float bz0) { _bz0 = bz0; } void setTracker(const Tracker* tracker) { _tracker = tracker; } + void setNoSignal(const HelixTool::StrawNoSignal& noSignal) { _noSignal = noSignal; } bool applyCuts(const HelixSeed& helix) { if(_debugLevel > 2) std::cout << "Testing helix:\n"; @@ -121,7 +124,7 @@ namespace mu2e if(_minPt > 0. && pt < _minPt) return false; // initialize the helix tool - HelixTool hTool(&helix, _tracker); + HelixTool hTool(&helix, _tracker, _noSignal); // geometric selections const Helicity hel = helix.helix().helicity(); @@ -148,16 +151,15 @@ namespace mu2e const float chi2XY = helix.helix().chi2dXY(); const float chi2PhiZ = helix.helix().chi2dZPhi(); const int nhits = hTool.nstrawhits(); - const float hRatio = hTool.hitRatio(); const bool hasCluster = helix.caloCluster().isNonnull(); if(_debugLevel > 2) std::cout << " chi2XY = " << chi2XY << " chi2PhiZ = " << chi2PhiZ - << " nhits = " << nhits << " hitRatio = " << hRatio << " hasCC = " << hasCluster << std::endl; + << " nhits = " << nhits << " hitRatio = " << hTool.hitRatio() << " hasCC = " << hasCluster << std::endl; if(chi2XY > _maxChi2XY) return false; if(chi2PhiZ > _maxChi2PhiZ) return false; if(nhits < _minNStrawHits) return false; - if(hRatio < _minHitRatio) return false; if(_requireCaloCluster && !hasCluster) return false; + if(_minHitRatio > 0. && hTool.hitRatio() < _minHitRatio) return false; // last: the most expensive if(_debugLevel > 2) std::cout << " --> passed all cuts\n"; // passes all cuts @@ -186,6 +188,7 @@ namespace mu2e double _bz0 = 1. ; const Tracker* _tracker = nullptr; + HelixTool::StrawNoSignal _noSignal; // empty: all straws live int _debugLevel = 0 ; }; @@ -202,6 +205,7 @@ namespace mu2e fhicl::Atom minArcGap { Name("minArcGap"), Comment("Minimum distance between helix arcs of same curve direction"), -1.}; fhicl::Sequence> HelixCuts { Name("HelixCuts"), Comment("Helix cut sets")}; fhicl::Atom doHistograms { Name("doHistograms"), Comment("Make debug histograms"), false}; + fhicl::Atom useTrackerStatus { Name("useTrackerStatus"), Comment("don't count straws without signal (TrackerStatus) as expected hits in the hit ratio"), true}; }; @@ -258,9 +262,12 @@ namespace mu2e // cut sets std::vector _helixCuts; bool _doHistograms; + bool _useTrackerStatus; double _bz0; const Tracker* _tracker; + ProditionsHandle _trackerStatus_h; + HelixTool::StrawNoSignal _noSignal; // empty: all straws live // counters unsigned _nevt, _npass; @@ -280,6 +287,7 @@ namespace mu2e , _maxOverlap (config().maxOverlap()) , _minArcGap (config().minArcGap()) , _doHistograms (config().doHistograms()) + , _useTrackerStatus (config().useTrackerStatus()) , _nevt (0) , _npass (0) { @@ -353,7 +361,7 @@ namespace mu2e float helSum = 0.; for(auto index : indices) { const HelixSeed& helix = hcol->at(index); - HelixTool hTool(&helix, _tracker); + HelixTool hTool(&helix, _tracker, _noSignal); constexpr float mmToMeV = 3./10.; Hist->p ->Fill(_bz0 * mmToMeV * helix.helix().momentum()); @@ -516,6 +524,13 @@ namespace mu2e bool MultiHelixFilter::filter(art::Event& event){ ++_nevt; + // straws without signal are not counted as expected hits in the helix hit ratio + if(_useTrackerStatus) { + const TrackerStatus* status = &_trackerStatus_h.get(event.id()); + _noSignal = [status](StrawId const& sid) { return status->noSignal(sid) || status->suppress(sid); }; + for(auto& cuts : _helixCuts) cuts.setNoSignal(_noSignal); + } + // create the output std::unique_ptr triginfo(new TriggerInfo); diff --git a/TrkFilters/src/SConscript b/TrkFilters/src/SConscript index 2cf1887502..0b961dd325 100644 --- a/TrkFilters/src/SConscript +++ b/TrkFilters/src/SConscript @@ -16,6 +16,8 @@ mainlib = helper.make_mainlib ( [] ) # Fixme: split into link lists for each module. helper.make_plugins( [ mainlib, 'mu2e_Mu2eUtilities', + 'mu2e_TrackerConditions', + 'mu2e_DbTables', 'mu2e_SeedService', 'mu2e_GeometryService', 'mu2e_BFieldGeom',