diff --git a/CalPatRec/inc/CalHelixFinderData.hh b/CalPatRec/inc/CalHelixFinderData.hh index 858093d042..874bffc5bd 100644 --- a/CalPatRec/inc/CalHelixFinderData.hh +++ b/CalPatRec/inc/CalHelixFinderData.hh @@ -132,33 +132,27 @@ namespace mu2e { //----------------------------------------------------------------------------- // circle parameters; the z center is ignored. //----------------------------------------------------------------------------- - ::LsqSums4 _sxy; - ::LsqSums4 _szphi; + ::LsqSums4 _sxy; // circle fitter + ::LsqSums4 _szphi; // dphi / dz line fitter - XYZVectorF _center; - float _radius; + XYZVectorF _center; // circle fit results + float _radius; + float _circle_chisq_dof; // circle fit quality - // float _chi2; -//----------------------------------------------------------------------------- -// 2015-02-06 P.Murat: fit with non-equal weights - XY-only -//----------------------------------------------------------------------------- - // ::LsqSums4 _sxyw; - // XYZVectorF _cw; - // float _rw; - // float _chi2w; //----------------------------------------------------------------------------- // Z parameters; dfdz is the slope of phi vs z (=-sign(1.0,qBzdir)/(R*tandip)), // fz0 is the phi value of the particle where it goes through z=0 // note that dfdz has a physical ambiguity in q*zdir. //----------------------------------------------------------------------------- - float _dfdz; - float _fz0; + float _dfdz; // dphi / dz fit results + float _fz0; + float _dfdz_chisq_dof; // dphi / dz fit quality //----------------------------------------------------------------------------- // diagnostics, histogramming //----------------------------------------------------------------------------- Diag_t _diag; //----------------------------------------------------------------------------- -// structure used to organize thei strawHits for the pattern recognition +// structure used to organize the strawHits for the pattern recognition //----------------------------------------------------------------------------- // PanelZ_t _oTracker[kNTotalPanels]; // std::array _hitsUsed; diff --git a/CalPatRec/inc/CalHelixFinder_module.hh b/CalPatRec/inc/CalHelixFinder_module.hh index d9d04fd100..b2036e591d 100644 --- a/CalPatRec/inc/CalHelixFinder_module.hh +++ b/CalPatRec/inc/CalHelixFinder_module.hh @@ -75,7 +75,6 @@ namespace mu2e { // event object labels //----------------------------------------------------------------------------- std::string _shLabel ; // MakeStrawHit label (makeSH) - // std::string _shpLabel; std::string _timeclLabel; int _minNHitsTimeCluster; //min nhits within a TimeCluster after check of Delta-ray hits @@ -96,7 +95,6 @@ namespace mu2e { const ComboHitCollection* _chcol; const TimeClusterCollection* _timeclcol; - HelixTraj* _helTraj; CalHelixFinderAlg _hfinder; CalHelixFinderData _hfResult; std::vector _hels; // helicity values to fit @@ -161,6 +159,10 @@ namespace mu2e { int goodHitsTimeCluster(const TimeCluster* TimeCluster); void pickBestHelix(std::vector& HelVec, int &Index_best); + + void fillDiagnosticInfo(const std::vector& helix_seed_vec, + const std::vector& nHitsRatio_vec, + int index_best, const int nGoodTClusterHits); }; } #endif diff --git a/CalPatRec/inc/CalHelixFinder_types.hh b/CalPatRec/inc/CalHelixFinder_types.hh index 1338af1436..59141946d2 100644 --- a/CalPatRec/inc/CalHelixFinder_types.hh +++ b/CalPatRec/inc/CalHelixFinder_types.hh @@ -39,7 +39,6 @@ namespace mu2e { struct Data_t { const art::Event* event; - std::string shLabel; enum { kMaxSeeds = 100, kMaxHits = 200 }; diff --git a/CalPatRec/src/CalHelixFinderAlg.cc b/CalPatRec/src/CalHelixFinderAlg.cc index edaa872271..b0c8d3400e 100644 --- a/CalPatRec/src/CalHelixFinderAlg.cc +++ b/CalPatRec/src/CalHelixFinderAlg.cc @@ -318,11 +318,6 @@ namespace mu2e { //----------------------------------------------------------------------------- if (_diag > 0) saveResults(Helix, 0); - // if (_filter) { - // filterDist(Helix); - // if (_diag > 0) saveResults(_xyzp,Helix,1); - // } - doPatternRecognition(Helix); //--------------------------------------------------------------------------- // 2014-11-11 gianipez changed the following if() statement to test the @@ -331,7 +326,7 @@ namespace mu2e { //--------------------------------------------------------------------------- if (_debug != 0) { printf("[CalHelixFinderAlg::findHelix] Helix._nXYSh = %i goodPointsTrkCandidate = %i\n", - Helix._nXYSh, Helix._nStrawHits);//_goodPointsTrkCandidate); + Helix._nXYSh, Helix._nStrawHits); } if (Helix._nStrawHits < _minNHits ) { @@ -340,10 +335,12 @@ namespace mu2e { else if ((Helix._radius < _rmin) || (Helix._radius > _rmax)) { Helix._fit = TrkErrCode(TrkErrCode::fail,2); // initialization failure } - else if ((Helix._nXYSh < _minNHits) || (Helix._sxy.chi2DofCircle() > _chi2xyMax)) { + else if ((Helix._nXYSh < _minNHits) || (Helix._circle_chisq_dof > _chi2xyMax)) { Helix._fit = TrkErrCode(TrkErrCode::fail,3); // xy reconstruction failure } - else if ((Helix._nZPhiSh < _minNHits) || (Helix._szphi.chi2DofLine() > _chi2zphiMax) || + else if ((Helix._nZPhiSh < _minNHits) || + (Helix._dfdz_chisq_dof < 0.f) || // no valid line fit was found + (Helix._dfdz_chisq_dof > _chi2zphiMax) || (fabs(Helix._dfdz) < _minDfDz) || (fabs(Helix._dfdz) > _maxDfDz)) { Helix._fit = TrkErrCode(TrkErrCode::fail,4); // phi-z reconstruction failure } @@ -534,7 +531,7 @@ namespace mu2e { //----------------------------------------------------------------------------- for (int n=nmin; n<=nmax; n++) { const float x = dphidz + n*dx; - const int bin = std::min(nbinsX-1, int((x-minX)/stepX)); + const int bin = std::max(0, std::min(nbinsX-1, int((x-minX)/stepX))); hist[bin] += weight; } } @@ -1005,6 +1002,8 @@ namespace mu2e { } CHECK_RESIDUALS:; + constexpr bool endAtThreshold(false); // Whether or not to keep iterating down to the threshold or stop once it's acceptable + if (endAtThreshold && Helix._szphi.chi2DofLine() <= _chi2zphiMax) goto NEXT_ITERATION_END; dphi_max = _maxXDPhi; //reset the coordinates of the worst hit iworst.face = -1; @@ -1046,6 +1045,7 @@ namespace mu2e { goto NEXT_ITERATION; } } + NEXT_ITERATION_END:; } } //----------------------------------------------------------------------------- @@ -1063,6 +1063,7 @@ namespace mu2e { else if (success) { // update helix results Helix._fz0 = Helix._szphi.phi0(); Helix._dfdz = Helix._szphi.dfdz(); + Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); } if ((SeedIndex.face == 0 ) && (SeedIndex.panel == 0) && (SeedIndex.panelHitIndex == 0) && (_diag > 0)) { @@ -1151,6 +1152,9 @@ namespace mu2e { if (!success) { Helix._hitsUsed = hitsUsed; + // _szphi was cleared and partially rebuilt during the failed fit; + // invalidate the cached chi2 field so it cannot disagree with _szphi + Helix._dfdz_chisq_dof = -1.f; } return success; @@ -1461,6 +1465,7 @@ namespace mu2e { //----------------------------------------------------------------------------- Helix._center.SetXYZ(x0, y0, 0.0); Helix._radius = radius; + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); Helix._nComboHits += rescuedPoints; Helix._nStrawHits += rescuedStrawHits; @@ -1518,13 +1523,14 @@ namespace mu2e { void CalHelixFinderAlg::filterUsingPatternRecognition(CalHelixFinderData& Helix) { if (Helix._seedIndex.panel < 0) return; + if (Helix._dfdz == 0.) return; // undefined helix results int nActive(0), nActive_hel(0); int nSh = Helix._nFiltStrawHits; float straw_mean_radius(0), chi2_global_helix(0), total_weight(0); float x_center(Helix._center.x()), y_center(Helix._center.y()), radius(Helix._radius); - float fz0(Helix._fz0), lambda(Helix._dfdz != 0. ? 1./Helix._dfdz : 0.); + float fz0(Helix._fz0), dfdz(Helix._dfdz); XYZVectorF hel_pred(0., 0., 0.); PanelZ_t* panelz(0); @@ -1563,7 +1569,7 @@ namespace mu2e { float x = hit->pos().x(); float y = hit->pos().y(); float z = Helix._zFace[f]; - float phi_pred = fz0 + z/lambda; + float phi_pred = fz0 + z*dfdz; float x_pred = x_center + radius*cos(phi_pred); float y_pred = y_center + radius*sin(phi_pred); hel_pred.SetX(x_pred); @@ -1767,6 +1773,9 @@ namespace mu2e { rescueHits(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), usePhiResid); if ((Helix._nXYSh - 1) != Helix._nZPhiSh) rc = doLinearFitPhiZ(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), useIntelligentWeight);//the factor "-1" takes into account that the XY fit includes the target center + // rescueHits may have added points to _szphi directly (UsePhiResiduals==1); + // if doLinearFitPhiZ was skipped above, sync the field from the live fitter + else if(usePhiResid) Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); if (_debug != 0) printInfo(Helix); //-------------------------------------------------------------------------------------------------------------- @@ -1781,6 +1790,7 @@ namespace mu2e { if (rs == 1) { // update Helix Z-phi part Helix._dfdz = _hdfdz; Helix._fz0 = _hphi0; + Helix._dfdz_chisq_dof = -1.f; } } @@ -1790,6 +1800,9 @@ namespace mu2e { usePhiResid = 1; rescueHits(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), usePhiResid); if ((Helix._nXYSh - 1) != Helix._nZPhiSh) rc = doLinearFitPhiZ(Helix, HitInfo_t(StrawId::_ntotalfaces-1,0,-1), useIntelligentWeight); //the factor "-1" takes into account that the XY fit includes the target center + // rescueHits may have added points to _szphi directly (UsePhiResiduals==1); + // if doLinearFitPhiZ was skipped above, sync the field from the live fitter + else if(usePhiResid) Helix._dfdz_chisq_dof = Helix._szphi.chi2DofLine(); if (_debug != 0) printInfo(Helix); strcpy(banner,"refineHelixParameters-after-doLinearFitPhiZ"); @@ -1921,9 +1934,12 @@ namespace mu2e { Helix._nXYSh = 0; Helix._nComboHits = 0; - Helix._sxy.addPoint(fCaloX,fCaloY,1./100.); - Helix._nXYSh += 1; - Helix._nComboHits += 1; + // Only add the calo cluster position if one is associated with this time cluster + if(fCaloTime > 0.) { + Helix._sxy.addPoint(fCaloX,fCaloY,1./100.); + Helix._nXYSh += 1; + Helix._nComboHits += 1; + } //------------------------------------------------------------------------------- // add stopping target center with a position error of 100 mm/sqrt(12) ~ 30mm => wt = 1/900 //------------------------------------------------------------------------------- @@ -1987,6 +2003,7 @@ namespace mu2e { Helix._center.SetX(Helix._sxy.x0()); Helix._center.SetY(Helix._sxy.y0()); Helix._radius = Helix._sxy.radius(); + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); if (_debug > 5) { @@ -2103,7 +2120,7 @@ namespace mu2e { chi2 = sxy.chi2DofCircle(); - if ((chi2 < chi2_min) || ( (i == SeedIndex.panelHitIndex) && (p == SeedIndex.panel) && (f == SeedIndex.face)) ) { + if (chi2 < chi2_min) { chi2_min = chi2; IWorst.face = f; IWorst.panel = p; @@ -2287,11 +2304,10 @@ namespace mu2e { //----------------------------------------------------------------------------- // update circle parameters //----------------------------------------------------------------------------- - // Trk._sxyw = sxyw; Trk._center.SetX(Trk._sxy.x0()); Trk._center.SetY(Trk._sxy.y0()); Trk._radius = Trk._sxy.radius(); - // Trk._chi2 = Trk._sxy.chi2DofCircle(); + Trk._circle_chisq_dof = Trk._sxy.chi2DofCircle(); }else { Trk._hitsUsed = hitsUsed; //restore the info of the used-hits that was originally passed to the procedure @@ -2590,7 +2606,7 @@ namespace mu2e { Helix._center.SetX(Helix._sxy.x0()); Helix._center.SetY(Helix._sxy.y0()); Helix._radius = Helix._sxy.radius(); - // Helix._chi2 = Helix._sxy.chi2DofCircle(); + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); F_END:; if (_debug > 5 ) { @@ -2945,10 +2961,12 @@ namespace mu2e { Helix._nXYSh = NPoints; Helix._radius = sxy.radius(); Helix._center.SetXYZ(sxy.x0(), sxy.y0(), 0.0); + Helix._circle_chisq_dof = sxy.chi2DofCircle(); Helix._nStrawHits = NPoints; Helix._nComboHits = NComboHits; Helix._dfdz = dfdz; Helix._fz0 = phi0 - dfdz*z_phi0; // *FLOAT_CHECK* + Helix._dfdz_chisq_dof = -1.f; radius_end = Helix._radius; //breakpoint -- radius before refine $ @@ -2958,11 +2976,8 @@ namespace mu2e { //----------------------------------------------------------------------------- if (rc < 0) return; //breakpoint -- radius after refine $ - // Helix._center.SetXYZ(Helix._cw.x(), Helix._cw.y(), 0.0); - // Helix._radius = Helix._rw; radius_end = Helix._radius; - // Helix._sxy = Helix._sxyw; // doWeightedCircleFit still adds the ST and the cluster Chi2 = Helix._sxy.chi2DofCircle(); NPoints = Helix._nStrawHits; // *FIXME* in principle, the fit can remove ST as well as the cluster @@ -2972,9 +2987,10 @@ namespace mu2e { // 2015-01-22 G. Pezzullo and P. Murat; update the dfdz value using all hits //----------------------------------------------------------------------------- int rs = findDfDz(Helix, SeedIndex); - if (rs ==1 ) { + if (rs ==1 ) { // 1 = success Helix._dfdz = _hdfdz; Helix._fz0 = _hphi0; + Helix._dfdz_chisq_dof = -1.f; // fill diag vector dfdzRes[1] = _hdfdz; dphi0Res[1] = _hphi0; @@ -2982,10 +2998,18 @@ namespace mu2e { //----------------------------------------------------------------------------- // 2015-01-23 G. Pezzu and P. Murat: when it fails, doLinearFitPhiZ returns negative value // in this case, use the previous value for dfdz and phi0 +// NOTE: _hphi0 is initialised to -9999. in resetTrackParamters() and is only updated by +// findDfDz() when >=2 stations have hits. If findDfDz() also failed, the fallback +// sets phi0_end = -9999., which is a sentinel value and NOT a valid phi0. +// The helix candidate is still stored in this state and compared against other +// candidates in searchBestTriplet; the final quality cuts in fitHelix will reject it +// only if the sentinel survives as _fz0 through to that check. +// A future improvement would be to return early from findTrack when rcPhiZ is false +// AND _hphi0 is still the sentinel value. //----------------------------------------------------------------------------- bool rcPhiZ = doLinearFitPhiZ(Helix, SeedIndex); - if (rcPhiZ) { + if (rcPhiZ) { // fit was successful dfdz_end = Helix._dfdz; phi0_end = Helix._fz0; // fill diagnostic vector @@ -2998,27 +3022,6 @@ namespace mu2e { //FIXME! implement a function countUsedHits(Helix, SeedIndex, NComboHits, NPoints); - // for (int f=SeedIndex.face; fpanelZs[p]; - // int nhitsPerPanel = panelz->nChHits(); - // int seedPanelIndex(0); - // if (nhitsPerPanel == 0) continue; - // if ( (f==SeedIndex.face) && (p==SeedIndex.panel) && (SeedIndex.panelHitIndex >=0)) seedPanelIndex = SeedIndex.panelHitIndex - panelz->idChBegin; - - // for (int i=seedPanelIndex; iidChBegin + i; - // hit = &Helix._chHitsToProcess[index]; - // if (Helix._hitsUsed[index] > 0 ) { - // ++NComboHits; - // NPoints += hit->nStrawHits(); - // } - // } - // }//endl panels loop - // } } else { dfdz_end = _hdfdz; @@ -3075,8 +3078,12 @@ namespace mu2e { //----------------------------------------------------------------------------- Helix._center.SetXYZ(Helix._sxy.x0(),Helix._sxy.y0(), 0.0); Helix._radius = radius_end; + Helix._circle_chisq_dof = Helix._sxy.chi2DofCircle(); Helix._fz0 = phi0_end; Helix._dfdz = dfdz_end; + // _dfdz_chisq_dof is already set by doLinearFitPhiZ when rcPhiZ==true; + // only reset it to 0 when the fit failed (histogram estimate, no valid chi2) + if (!rcPhiZ) Helix._dfdz_chisq_dof = -1.f; if (_diag > 0){ Helix._diag.loopId_4 = _findTrackLoopIndex; @@ -3106,7 +3113,7 @@ namespace mu2e { int firstPanel(0); if (f == SeedIndex.face) firstPanel = SeedIndex.panel; for (int p=firstPanel; ppanelZs[p];//Helix._oTracker[p]; + panelz = &facez->panelZs[p]; int nhitsPerPanel = panelz->nChHits(); int seedPanelIndex(0); if (nhitsPerPanel == 0) continue; @@ -3116,17 +3123,8 @@ namespace mu2e { int index = panelz->idChBegin + i; hit = &Helix._chHitsToProcess[index]; - // int index = facez->evalUniqueHitIndex(f,p,i);//p*FaceZ_t::kNMaxHitsPerPanel + i; if (Helix._hitsUsed[index] != 1) continue; - // if (j < Helix.maxIndex()) { - // Helix._diag.dist[j] = hit->_drFromPred; - // Helix._diag.dz [j] = hit->_dzFromSeed; - // ++j; - // } - // else { - // printf("ERROR in CalHelixFinderAlg::findTrack : index out limits. IGNORE; \n"); - // } }//end loop over the hits within a panel }//end panel loop }//end face loop diff --git a/CalPatRec/src/CalHelixFinderData.cc b/CalPatRec/src/CalHelixFinderData.cc index 7b91201e47..d2126ffa35 100644 --- a/CalPatRec/src/CalHelixFinderData.cc +++ b/CalPatRec/src/CalHelixFinderData.cc @@ -113,9 +113,11 @@ namespace mu2e { _szphi.clear(); _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nFiltPoints = 0; _nFiltStrawHits = 0; @@ -185,9 +187,11 @@ void CalHelixFinderData::clearHelixInfo() { _szphi.clear(); _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nXYSh = 0; _nZPhiSh = 0; @@ -216,9 +220,11 @@ void CalHelixFinderData::clearHelixInfo() { _szphi.clear(); // _chi2 = -1.; _radius = -1.; + _circle_chisq_dof = 1e10; _dfdz = -1.e6; _fz0 = -1.e6; + _dfdz_chisq_dof = 1e10; _nXYSh = 0; diff --git a/CalPatRec/src/CalHelixFinder_module.cc b/CalPatRec/src/CalHelixFinder_module.cc index 697e5f5de6..ecad489323 100644 --- a/CalPatRec/src/CalHelixFinder_module.cc +++ b/CalPatRec/src/CalHelixFinder_module.cc @@ -20,7 +20,6 @@ #include "Offline/GeometryService/inc/DetectorSystem.hh" #include "Offline/CalorimeterGeom/inc/DiskCalorimeter.hh" #include "Offline/ConfigTools/inc/ConfigFileLookupPolicy.hh" -// #include "CalPatRec/inc/KalFitResult.hh" #include "Offline/RecoDataProducts/inc/StrawHitIndex.hh" #include "Offline/RecoDataProducts/inc/HelixHit.hh" @@ -41,6 +40,8 @@ #include "TSystem.h" #include "TInterpreter.h" +#include + using namespace boost::accumulators; using CLHEP::HepVector; using CLHEP::Hep3Vector; @@ -65,11 +66,14 @@ namespace mu2e { consumes(_shLabel); consumes(_timeclLabel); + // Initialize the list of helicities to consider std::vector helvals = config().Helicities(); for(auto hv : helvals) { Helicity hel(hv); _hels.push_back(hel); } + + // In single output, put all helices by all helicities into one helix collection, otherwise produce a collection per helicity if (_doSingleOutput){ produces(); } else { @@ -80,11 +84,8 @@ namespace mu2e { _tpart = (TrkParticle::type)_fitparticle; //----------------------------------------------------------------------------- -// provide for interactive disanostics +// provide for interactive diagnostics //----------------------------------------------------------------------------- - _helTraj = 0; - - _data.shLabel = _shLabel; if (_debugLevel != 0) _printfreq = 1; @@ -96,13 +97,14 @@ namespace mu2e { // destructor //----------------------------------------------------------------------------- CalHelixFinder::~CalHelixFinder() { - if (_helTraj) delete _helTraj; } //----------------------------------------------------------------------------- void CalHelixFinder::beginJob(){ - art::ServiceHandle tfs; - _hmanager->bookHistograms(tfs); + if(_diagLevel != 0) { // only book histograms if running diagnostics + art::ServiceHandle tfs; + _hmanager->bookHistograms(tfs); + } } //----------------------------------------------------------------------------- @@ -182,16 +184,6 @@ namespace mu2e { _shLabel.data()); } - // art::Handle shposH; - // if (evt.getByLabel(_shpLabel,shposH)) { - // _shpcol = shposH.product(); - // } - // else { - // _shpcol = 0; - // printf(" >>> ERROR in CalHelixFinder::findData: StrawHitPositionCollection with label=%s not found.\n", - // _shpLabel.data()); - // } - if (evt.getByLabel(_timeclLabel, _timeclcolH)) { _timeclcol = _timeclcolH.product(); } @@ -203,8 +195,7 @@ namespace mu2e { //----------------------------------------------------------------------------- // done //----------------------------------------------------------------------------- - return (_chcol != nullptr) //&& (_shpcol != nullptr) - && (_timeclcol != nullptr); + return (_chcol != nullptr) && (_timeclcol != nullptr); } //----------------------------------------------------------------------------- @@ -228,11 +219,11 @@ namespace mu2e { _data.nseeds [counter] = 0; ++counter; } - }else { + } else { helcols[0] = std::unique_ptr(new HelixSeedCollection()); _data.nseeds [counter] = 0; } - // unique_ptr outseeds(new HelixSeedCollection); + //----------------------------------------------------------------------------- // find the data //----------------------------------------------------------------------------- @@ -246,8 +237,6 @@ namespace mu2e { _hfResult._tpart = _tpart; _hfResult._fdir = _fdir; _hfResult._chcol = _chcol; - // _hfResult._shpos = _shpcol; - //_hfResult._shfcol = _shfcol; _data.nTimePeaks = _timeclcol->size(); for (int ipeak=0; ipeak<_data.nTimePeaks; ipeak++) { @@ -262,7 +251,7 @@ namespace mu2e { // create track definitions for the helix fit from this initial information // track fitting objects for this peak //----------------------------------------------------------------------------- - _hfResult.clearTempVariables();//clearTimeClusterInfo(); + _hfResult.clearTempVariables(); _hfResult._timeCluster = tc; _hfResult._timeClusterPtr = art::Ptr(_timeclcolH,ipeak); @@ -335,102 +324,16 @@ namespace mu2e { } } - // helix_seed_vec[index_best]._status.merge(TrkFitFlag::helixOK); - // outseeds->push_back(helix_seed_vec[index_best]); if (_diagLevel > 0) { -//-------------------------------------------------------------------------------- -// fill diagnostic information -//-------------------------------------------------------------------------------- - int nhitsMin(15); - double mm2MeV = (3/10.)*_bz0; - - int loc = _data.nseeds[0]; - if (loc < _data.maxSeeds()) { - if (index_best == 2){ - index_best = 0; - } - int nhits = helix_seed_vec[index_best]._hhits.size(); - _data.ntclhits[loc]= nGoodTClusterHits; - _data.nhits[loc] = nhits; - _data.radius[loc] = helix_seed_vec[index_best].helix().radius(); - _data.pT[loc] = mm2MeV*_data.radius[loc]; - _data.p[loc] = _data.pT[loc]/std::cos( std::atan(helix_seed_vec[index_best].helix().lambda()/_data.radius[loc])); - - _data.chi2XY[loc] = _hfResult._sxy.chi2DofCircle(); - _data.chi2ZPhi[loc] = _hfResult._szphi.chi2DofLine(); - - _data.nHitsRatio[loc] = nHitsRatio_vec[index_best]; - _data.eDepAvg[loc] = helix_seed_vec[index_best]._eDepAvg; - - _data.nseeds[0]++; - _data.good[loc] = 0; - if (nhits >= nhitsMin) { - _data.nseeds[1]++; - _data.good[loc] = 1; - } - _data.nStationPairs[loc] = _hfResult._diag.nStationPairs; - - _data.dr [loc] = _hfResult._diag.dr; - _data.shmeanr [loc] = _hfResult._diag.straw_mean_radius; - _data.chi2d_helix [loc] = _hfResult._diag.chi2d_helix; - if (_hfResult._diag.chi2d_helix>3) printf("[%s] : chi2Helix = %10.3f event number %8i\n", oname,_hfResult._diag.chi2d_helix,_iev); -//----------------------------------------------------------------------------- -// info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition -//----------------------------------------------------------------------------- - _data.loopId [loc] = _hfResult._diag.loopId_4; - if (_hfResult._diag.loopId_4 == 1) { - _data.chi2d_loop0 [loc] = _hfResult._diag.chi2_dof_circle_12; - _data.chi2d_line_loop0 [loc] = _hfResult._diag.chi2_dof_line_13; - _data.npoints_loop0 [loc] = _hfResult._diag.n_active_11; - - } - if (_hfResult._diag.loopId_4 == 2){ - _data.chi2d_loop1 [loc] = _hfResult._diag.chi2_dof_circle_12; - _data.chi2d_line_loop1 [loc] = _hfResult._diag.chi2_dof_line_13; - _data.npoints_loop1 [loc] = _hfResult._diag.n_active_11; - } - -//-------------------------------------------------------------------------------- -// info of the track candidate during the CAlHelixFinderAlg::findTrack loop -//-------------------------------------------------------------------------------- - int counter(0); - for (unsigned i=0; i<_hfResult._hitsUsed.size(); ++i){ - if (_hfResult._hitsUsed[i] != 1) continue; - ++counter; - } - // for (int f=0; fpanelZs[p];//&_hfResult._oTracker[p]; - // int nhits = panelz->fNHits; - // if (nhits == 0) continue; - - // for (int i=0; i_chHitsToProcess.at(i); - // int index = facez->evalUniqueHitIndex(f,p,i);//p*CalHelixFinderData::kNMaxHitsPerPanel + i; - // if (_hfResult._hitsUsed[index] != 1) continue; - - // // double dzFromSeed = hit->_dzFromSeed; //distance form the hit used to seed the 3D-search - // // double drFromPred = hit->_drFromPred; //distance from prediction - // // _data.hitDzSeed[loc][counter] = dzFromSeed; - // // _data.hitDrPred[loc][counter] = drFromPred; - // ++counter; - // }//end loop over the hits within a panel - // }//end panels loop - // }//end faces loop - } - else { - printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); - } + fillDiagnosticInfo(helix_seed_vec, nHitsRatio_vec, + index_best, nGoodTClusterHits); } - } //-------------------------------------------------------------------------------- // fill histograms //-------------------------------------------------------------------------------- - if (_diagLevel > 0) { - _hmanager->fillHistograms(&_data); - } + if (_diagLevel > 0) _hmanager->fillHistograms(&_data); + //----------------------------------------------------------------------------- // put reconstructed tracks into the event record //----------------------------------------------------------------------------- @@ -451,8 +354,6 @@ namespace mu2e { // //----------------------------------------------------------------------------- void CalHelixFinder::endJob() { - // does this cause the file to close? - art::ServiceHandle tfs; } //-------------------------------------------------------------------------------- @@ -465,8 +366,6 @@ namespace mu2e { double helixRadius = 1./fabs(hel->omega()); double impactParam = hel->d0(); double phi0 = hel->phi0(); - // double x0 = -(helixRadius + impactParam)*sin(phi0)*sig; - // double y0 = (helixRadius + impactParam)*cos(phi0)*sig; double x0 = -(1/hel->omega()+impactParam)*sin(phi0); double y0 = (1/hel->omega()+impactParam)*cos(phi0); @@ -504,11 +403,7 @@ namespace mu2e { // cluster hits assigned to the reconsturcted Helix int nhits = HfResult.nGoodHits(); - // printf("[CalHelixFinder::initHelixSeed] radius = %2.3f x0 = %2.3f y0 = %2.3f dfdz = %2.3e nhits = %i chi2XY = %2.3f chi2PHIZ = %2.3f\n", - // helixRadius, center.x(), center.y(), dfdz, nhits, HfResult._sxyw.chi2DofCircle(), HfResult._srphi.chi2DofLine()); - // printf("[CalHelixFinder::initHelixSeed] Index X Y Z PHI\n"); - // double z_start(0); HelSeed._hhits.setParent(_chcol->parent()); for (int i=0; inhits(); int ngoodhits(0); - // double minT(500.), maxT(2000.); for (int i=0; ihits().at(i); - // StrawHitFlag flag = _shfcol->at(index); ComboHit sh = _chcol ->at(index); StrawHitFlag flag = sh.flag(); int bkg_hit = flag.hasAnyProperty(StrawHitFlag::bkg); if (bkg_hit) continue; - // if ( (sh.time() < minT) || (sh.time() > maxT) ) continue; ngoodhits += sh.nStrawHits(); } @@ -616,86 +504,133 @@ namespace mu2e { return; } +//----------------------------------------------------------------------------- +// at Mu2e, 2 helices with different helicity could be duplicates of each other +//----------------------------------------------------------------------------- + const HelixSeed *h1, *h2; const ComboHitCollection *tlist, *clist; int nh1, nh2, natc(0); - const mu2e::HelixHit *hitt, *hitc; + constexpr float overlap_threshold(0.5f); // 50% overlap to match helices + + // Helix 1 h1 = &HelVec[0]; -//------------------------------------------------------------------------------ -// check if an AlgorithmID collection has been created by the process -//----------------------------------------------------------------------------- tlist = &h1->hits(); nh1 = tlist->size(); + // Helix 2 h2 = &HelVec[1]; -//----------------------------------------------------------------------------- -// at Mu2e, 2 helices with different helicity could be duplicates of each other -//----------------------------------------------------------------------------- clist = &h2->hits(); nh2 = clist->size(); //----------------------------------------------------------------------------- // check the number of common hits //----------------------------------------------------------------------------- - for (int k=0; kat(k); - for (int l=0; lat(l); - if (hitt->index() == hitc->index()) { - natc += 1; - break; - } - } + std::unordered_set hits_1; // unordered set for O(1) lookups + hits_1.reserve(nh1); + + // Fill the set + for(int k = 0; k < nh1; ++k) { + hits_1.insert(tlist->at(k).index()); + } + + // Count overlapping hits + for(int l = 0; l < nh2; l++) { + if(hits_1.contains(clist->at(l).index())) ++natc; } - if ((natc > nh1/2.) || (natc > nh2/2.)) { + // Check if at least one shares half of its hits with the other helix + if ((natc > nh1*overlap_threshold) || (natc > nh2*overlap_threshold)) { //----------------------------------------------------------------------------- // pick the helix with the largest number of hits //----------------------------------------------------------------------------- - if (nh2 > nh1) { -//----------------------------------------------------------------------------- -// h2 is a winner, no need to save h1 -//----------------------------------------------------------------------------- - Index_best = 1; - return; - } - else if (nh1 > nh2){ -//----------------------------------------------------------------------------- -// h1 is a winner, mark h2 in hope that it will be OK, continue looping -//----------------------------------------------------------------------------- - Index_best = 0; - return; - } -//----------------------------------------------------------------------------- -// in case they have the exact amount of hits, pick the one with better chi2dZphi -//----------------------------------------------------------------------------- - if (nh1 == nh2) { - float chi2dZphi_h1 = h1->helix().chi2dZPhi(); - float chi2dZphi_h2 = h2->helix().chi2dZPhi(); - if (chi2dZphi_h1 < chi2dZphi_h2){ - Index_best = 0; - return; - }else { - Index_best = 1; - return; - } - } - }else { + if (nh2 > nh1) Index_best = 1; // helix two has more hits + else if (nh1 > nh2) Index_best = 0; // helix one has more hits + else Index_best = (h1->helix().chi2dZPhi() < h2->helix().chi2dZPhi()) ? 0 : 1; // equal hits --> pick the best chi^2 + } else { //----------------------------------------------------------------------------- // this is the case where we consider the two helices independent, so we want // to store both //----------------------------------------------------------------------------- Index_best = 2; - return; } } +//----------------------------------------------------------------------------- +// fill diagnostic information +//----------------------------------------------------------------------------- + void CalHelixFinder::fillDiagnosticInfo(const std::vector& helix_seed_vec, + const std::vector& nHitsRatio_vec, + int index_best, const int nGoodTClusterHits) { + const char* oname = "CalHelixFinder::fillDiagnosticsInfo"; + int nhitsMin(15); + double mm2MeV = (3/10.)*_bz0; + + int loc = _data.nseeds[0]; + if (loc < _data.maxSeeds()) { + if (index_best == 2){ + index_best = 0; + } + int nhits = helix_seed_vec[index_best]._hhits.size(); + _data.ntclhits[loc]= nGoodTClusterHits; + _data.nhits[loc] = nhits; + _data.radius[loc] = helix_seed_vec[index_best].helix().radius(); + _data.pT[loc] = mm2MeV*_data.radius[loc]; + _data.p[loc] = _data.pT[loc]/std::cos( std::atan(helix_seed_vec[index_best].helix().lambda()/_data.radius[loc])); + + _data.chi2XY[loc] = _hfResult._sxy.chi2DofCircle(); + _data.chi2ZPhi[loc] = _hfResult._szphi.chi2DofLine(); + + _data.nHitsRatio[loc] = nHitsRatio_vec[index_best]; + _data.eDepAvg[loc] = helix_seed_vec[index_best]._eDepAvg; + + _data.nseeds[0]++; + _data.good[loc] = 0; + if (nhits >= nhitsMin) { + _data.nseeds[1]++; + _data.good[loc] = 1; + } + _data.nStationPairs[loc] = _hfResult._diag.nStationPairs; + + _data.dr [loc] = _hfResult._diag.dr; + _data.shmeanr [loc] = _hfResult._diag.straw_mean_radius; + _data.chi2d_helix [loc] = _hfResult._diag.chi2d_helix; + if (_hfResult._diag.chi2d_helix>3) printf("[%s] : chi2Helix = %10.3f event number %8i\n", oname,_hfResult._diag.chi2d_helix,_iev); + //----------------------------------------------------------------------------- + // info of the track candidate after the first loop with findtrack on CalHelixFinderAlg::doPatternRecognition + //----------------------------------------------------------------------------- + _data.loopId [loc] = _hfResult._diag.loopId_4; + if (_hfResult._diag.loopId_4 == 1) { + _data.chi2d_loop0 [loc] = _hfResult._diag.chi2_dof_circle_12; + _data.chi2d_line_loop0 [loc] = _hfResult._diag.chi2_dof_line_13; + _data.npoints_loop0 [loc] = _hfResult._diag.n_active_11; + + } + if (_hfResult._diag.loopId_4 == 2){ + _data.chi2d_loop1 [loc] = _hfResult._diag.chi2_dof_circle_12; + _data.chi2d_line_loop1 [loc] = _hfResult._diag.chi2_dof_line_13; + _data.npoints_loop1 [loc] = _hfResult._diag.n_active_11; + } + + //-------------------------------------------------------------------------------- + // info of the track candidate during the CAlHelixFinderAlg::findTrack loop + //-------------------------------------------------------------------------------- + int counter(0); + for (unsigned i=0; i<_hfResult._hitsUsed.size(); ++i){ + if (_hfResult._hitsUsed[i] != 1) continue; + ++counter; + } + } + else { + printf(" N(seeds) > %i, IGNORE SEED\n",_data.maxSeeds()); + } + } } using mu2e::CalHelixFinder;