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..c17a9a51ea 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") }; }; 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; + 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 054c99a9eb..c1997bd38f 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" @@ -114,6 +116,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)); @@ -126,7 +130,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(); @@ -136,7 +144,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 @@ -151,7 +158,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 + << " nLoops = " << nLoops << " hRatio = "<< helTool.hitRatio() << " slopeSignedSig = " << slopeSignedSig << " eDepAvg = " << edep << " mass = " << mass << std::endl; } if( Helix.status().hasAllProperties(_goodh) && @@ -174,7 +181,7 @@ namespace mu2e (!_useMaxMass || mass <= _maxmass) && (!_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(); @@ -220,6 +227,7 @@ namespace mu2e bool _prescaleUsingD0Phi = false; PhiPrescalingParams _prescalerPar; const Tracker*_myTracker = nullptr; + const TrackerStatus* _trackerStatus = nullptr; }; struct Config{ @@ -232,6 +240,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")}; }; @@ -256,6 +265,8 @@ namespace mu2e bool _noFilter; unsigned _minNHelices; float _maxDt0; + bool _useTrackerStatus; + ProditionsHandle _trackerStatus_h; bool checkHelixFromHelicity(const HelixSeed&helix); }; @@ -270,7 +281,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"; @@ -318,6 +330,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..98a92c7bc0 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")}; }; @@ -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',