diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt index 4dc59bec52a2862bd5b72822e44ab5fd365fa892..e985bf243958d3130cf854322d2dd43c43730d63 100644 --- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt +++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt @@ -1,19 +1,15 @@ -################################################################################ -# Package: MdtSegmentT0Fitter -################################################################################ +# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration # Declare the package name: atlas_subdir( MdtSegmentT0Fitter ) # External dependencies: -find_package( ROOT COMPONENTS Minuit Core Tree MathCore Hist RIO pthread MathMore Minuit2 Matrix Physics HistPainter Rint ) +find_package( ROOT COMPONENTS Core MathCore Minuit2 ) # Component(s) in the package: atlas_add_component( MdtSegmentT0Fitter - src/*.cxx - src/components/*.cxx - INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} - LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel MdtCalibSvcLib TrkDriftCircleMath MdtCalibData MuonCalibToolsLib MuonReadoutGeometry MuonIdHelpersLib MuonPrepRawData MuonRIO_OnTrack MuonSegmentMakerInterfacesLib ) - -# Install files from the package: -atlas_install_headers( MdtSegmentT0Fitter ) + MdtSegmentT0Fitter/*.h src/*.cxx src/components/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} AthenaBaseComps GaudiKernel MdtCalibSvcLib + TrkDriftCircleMath MdtCalibData MuonReadoutGeometry MuonPrepRawData + MuonRIO_OnTrack MuonSegmentMakerInterfacesLib ) diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx index 5c9ab78eaacfed57406362fdb1933ddc59b855a5..e5a9311aec54c486c9c34a6d8b1c591d3c0f4cf6 100755 --- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx +++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx @@ -57,15 +57,15 @@ namespace TrkDriftCircleMath { FunctionToMinimize(const int used) : m_data(),m_used(used),m_t0Error() {} ~FunctionToMinimize() {m_data.clear();} double DoEval(const double* xx) const { - double ang = xx[0]; + double ang = xx[0]; double b = xx[1]; double t0 = xx[2]; - + double cosin = std::cos(ang); double sinus = std::sin(ang); - + double fval = 0.; - // Add t0 constraint + // Add t0 constraint if (m_t0Error == WEAK_TOPO_T0ERROR ) { fval += xx[2]*xx[2]/(1.0 *m_t0Error*m_t0Error); } @@ -79,7 +79,7 @@ namespace TrkDriftCircleMath { double uppercut = m_data[i].rt->tUpper(); double lowercut = m_data[i].rt->tLower(); // Penalty for t<lowercut and t >uppercut - if (t> uppercut ) { // too large + if (t> uppercut ) { // too large fval += (t-uppercut)* (t-uppercut)*0.1; } else if (t < 0 ) {// too small fval += (t-lowercut)*(t-lowercut)*0.1; @@ -112,7 +112,7 @@ namespace TrkDriftCircleMath { MdtSegmentT0Fitter::MdtSegmentT0Fitter(const std::string& ty,const std::string& na,const IInterface* pa) : AthAlgTool(ty,na,pa), - DCSLFitter(), + DCSLFitter(), m_ntotalCalls(0), m_npassedNHits(0), m_npassedSelectionConsistency(0), @@ -121,14 +121,14 @@ namespace TrkDriftCircleMath { m_npassedMinuitFit(0) { declareInterface <IDCSLFitProvider> (this); } - + StatusCode MdtSegmentT0Fitter::initialize() { ATH_CHECK(m_calibrationDbTool.retrieve()); return StatusCode::SUCCESS; } - + StatusCode MdtSegmentT0Fitter::finalize() { - + double scaleFactor = m_ntotalCalls != 0 ? 1./(double)m_ntotalCalls : 1.; ATH_MSG_INFO( "Summarizing fitter statistics " << "\n" @@ -140,15 +140,15 @@ namespace TrkDriftCircleMath { << " Passed Fit " << std::setw(10) << m_npassedMinuitFit << " " << scaleFactor*m_npassedMinuitFit ); return StatusCode::SUCCESS; } - - + + /***********************************************************************************/ /// RT function from Craig Blocker /// ok for now, possibly replace with actual RT function used to calibrate run - + constexpr double T2R_A[] = {1.184169e-1, 3.32382e-2, 4.179808e-4, -5.012896e-6, 2.61497e-8, -7.800677e-11, 1.407393e-13, -1.516193e-16, 8.967997e-20, -2.238627e-23}; constexpr double RCORR_A[] = {234.3413, -5.803375, 5.061677e-2, -1.994959e-4, 4.017433e-7, -3.975037e-10, 1.522393e-13}; - + double rcorr(const double tin) { double rc; if(tin < 16.) { @@ -162,25 +162,25 @@ namespace TrkDriftCircleMath { } } rc = rc*1e-3; - + return rc; } - + double t2r(const double tin) { if(tin < 0.) return 0.; if(tin > MAX_DRIFT) return 20.; - + double tpow = 1.; double rc = 0.; for(int j = 0; j < 10; j++) { rc += T2R_A[j]*tpow; tpow *= tin; } - + rc -= rcorr(tin); return rc; } - + /// derivatives of RT function, use to get errors double rcorrprime(const double tin) { double rc; @@ -195,26 +195,26 @@ namespace TrkDriftCircleMath { } } rc = rc*1e-3; - + return rc; } - + double t2rprime(const double tin) { if(tin < 0.) return 0.; if(tin > MAX_DRIFT) return 20.; - + double tpow = 1.; double rc = 0.; for(int j = 1; j < 10; j++) { rc += j*T2R_A[j]*tpow; tpow *= tin; } - + rc -= rcorrprime(tin); return rc; } - - + + /// use a binary search to get rt-inverse from rt /// assumes the function is monotonic, obviously not true for these polynomial parametrizations for all t double r2t(const double r) { @@ -225,7 +225,7 @@ namespace TrkDriftCircleMath { } else if(r>t2r(tb)) { return R2TSPURIOUS; } - + int itr = 0; while (ta <= tb) { double tm = (ta + tb) / 2; // compute mid point. @@ -239,13 +239,13 @@ namespace TrkDriftCircleMath { else if (r < rtm ) { tb = tm; // repeat search in bottom half. } - + itr++; if(itr>50) return -1; } return -1; // failed to find key } - + double r2t_ext(std::vector<const MuonCalib::IRtRelation*> *rtpointers, double r, int i) { const MuonCalib::IRtRelation* rtrel = rtpointers->at(i); double ta = rtrel->tLower(); @@ -255,7 +255,7 @@ namespace TrkDriftCircleMath { } else if(r>rtrel->radius(tb)) { return R2TSPURIOUS; } - + int itr = 0; while (ta <= tb) { double tm = (ta + tb) / 2; // compute mid point. @@ -269,7 +269,7 @@ namespace TrkDriftCircleMath { else if (r < rtm ) { tb = tm; // repeat search in bottom half. } - + itr++; if(itr>50) return -1; } @@ -285,18 +285,12 @@ namespace TrkDriftCircleMath { if(m_trace) ATH_MSG_DEBUG("New seg: "); - ROOT::Math::Minimizer* minimum = new ROOT::Minuit2::Minuit2Minimizer("algoName"); - minimum->SetMaxFunctionCalls(1000000); - minimum->SetTolerance(0.001); - minimum->SetPrintLevel(-1); - if(msgLvl(MSG::VERBOSE)) minimum->SetPrintLevel(1); - const DCOnTrackVec& dcs_keep = dcs; unsigned int N = dcs_keep.size(); result.setT0Shift(-99999,-99999); - + if(N<2) { return false; } @@ -304,7 +298,7 @@ namespace TrkDriftCircleMath { if( selection.size() != N ) { ATH_MSG_ERROR("MdtSegmentT0Fitter.cxx:fit with t0 <bad HitSelection>"); return false; - } + } ++m_npassedSelectionConsistency; int used=0; for(unsigned int i=0;i<N;++i){ @@ -320,12 +314,12 @@ namespace TrkDriftCircleMath { // // Copy driftcircles and reset the drift radii as they might have been overwritten - // after a succesfull t0 fit - // + // after a succesfull t0 fit + // DCOnTrackVec::const_iterator it = dcs.begin(); DCOnTrackVec::const_iterator it_end = dcs.end(); - + DCOnTrackVec dcs_new; dcs_new.reserve(dcs.size()); double chi2p = 0.; @@ -346,7 +340,7 @@ namespace TrkDriftCircleMath { } } if(m_trace&&chi2p>0) ATH_MSG_DEBUG("NO Minuit Fit TOO few hits Chi2 penalty " << chi2p); - bool oldrefit = DCSLFitter::fit( result, line, dcs_new, selection ); + bool oldrefit = DCSLFitter::fit( result, line, dcs_new, selection ); chi2p += result.chi2(); // add chi2 penalty for too large or too small driftTimes t < 0 or t> t upper result.set(chi2p, result.ndof(), result.dtheta(), result.dy0()); @@ -356,15 +350,15 @@ namespace TrkDriftCircleMath { return oldrefit; } else { if(m_trace) ATH_MSG_DEBUG("FITTING FOR T0 N "<<N<<" used " << used); - } + } ++m_npassedMinHits; - + if (m_trace) { ATH_MSG_DEBUG(" in MdtSegmentT0Fitter::fit with N dcs "<< N << " hit selection size " << selection.size()); ATH_MSG_DEBUG("in fit "<<result.hasT0Shift()<< " " <<result.t0Shift()); } - + double Zc(0); double Yc(0); double S(0); @@ -372,7 +366,7 @@ namespace TrkDriftCircleMath { double Sy(0); std::vector<double> y(N); std::vector<double> z(N); - std::vector<double> w(N); + std::vector<double> w(N); std::vector<double> r(N); std::vector<double> dr(N); std::vector<double> t(N); @@ -393,7 +387,7 @@ namespace TrkDriftCircleMath { y[ii] = it->y(); z[ii] = it->x(); r[ii] = std::abs(roto->driftRadius()); - dr[ii] = it->dr(); + dr[ii] = it->dr(); const Muon::MdtPrepData *peerd; peerd = dynamic_cast<const Muon::MdtPrepData*>(roto->prepRawData()); if(!peerd) { @@ -405,7 +399,7 @@ namespace TrkDriftCircleMath { rtpointers[ii] = rtInfo->rt(); t[ii] = roto->driftTime(); - double newerror = m_scaleErrors ? it->drPrecise() : it->dr(); + double newerror = m_scaleErrors ? it->drPrecise() : it->dr(); if( newerror > 0.) w[ii] = 1./(newerror); else w[ii] = 0.; @@ -416,7 +410,7 @@ namespace TrkDriftCircleMath { } if(m_trace) ATH_MSG_DEBUG("DC: (" << y[ii] << "," << z[ii] << ") R = " << r[ii] << " W " << w[ii] <<" t " <<t[ii]<< " id: "<<it->id()<<" sel " << selection[ii]); - + if( selection[ii] == 0 ){ S+=w[ii]; Sz+= w[ii]*z[ii]; @@ -430,15 +424,15 @@ namespace TrkDriftCircleMath { const double inv_S = 1. / S; Zc = Sz*inv_S; Yc = Sy*inv_S; - + if(m_trace) ATH_MSG_DEBUG("Yc " << Yc << " Zc " << Zc); - + /// go to coordinates centered at the average of the hits for(unsigned int i=0;i<N;++i) { y[i] -= Yc; z[i] -= Zc; } - + int selcount(0); DCOnTrackVec::const_iterator it = dcs_keep.begin(); DCOnTrackVec::const_iterator it_end = dcs_keep.end(); @@ -448,13 +442,13 @@ namespace TrkDriftCircleMath { // tlower_i < ti - t0 < tupper_i double min_tlower=1e10; double max_tupper=-1e10; - + double t0seed=0; // the average t0 of the hit double st0 = 0; // the std deviation of the hit t0s double min_t0 = 1e10; // the smallest t0 seen double tee0, tl, th; - - + + for(int ii=0 ;it!=it_end; ++it, ++ii ){ if( selection[ii] == 0 ) { double r2tval = r2t_ext(&rtpointers, r[ii], ii) ; @@ -481,7 +475,7 @@ namespace TrkDriftCircleMath { minFunct.addCoords(z[ii], t[ii], y[ii], w[ii], r[ii], rtpointers[ii]); selcount++; - } + } } t0seed /= selcount; st0 = st0/selcount - t0seed*t0seed; @@ -493,7 +487,7 @@ namespace TrkDriftCircleMath { double theta = line.phi(); double cosin = std::cos(theta); double sinus = std::sin(theta); - + if ( sinus < 0.0 ) { sinus = -sinus; cosin = -cosin; @@ -501,17 +495,17 @@ namespace TrkDriftCircleMath { cosin = -cosin; } if(m_trace) ATH_MSG_DEBUG("before fit theta "<<theta<<" sinus "<<sinus<< " cosin "<< cosin); - + double d = line.y0() + Zc*sinus-Yc*cosin; - + if(m_trace) { ATH_MSG_DEBUG(" line x y "<<line.position().x()<<" "<<line.position().y()); ATH_MSG_DEBUG(" Zc Yc "<< Zc <<" "<<Yc); ATH_MSG_DEBUG(" line x0 y0 "<<line.x0()<<" "<<line.y0()); ATH_MSG_DEBUG(" hit shift " << -Zc*sinus+Yc*cosin); - } + } -// Calculate signed radii +// Calculate signed radii int nml1p = 0; int nml2p = 0; @@ -539,12 +533,12 @@ namespace TrkDriftCircleMath { if(m_trace) ATH_MSG_DEBUG("Combination rejected for positive radii ML1 " << nml1p << " ML2 " << nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << used << " t0 Error " << t0Error); it = dcs.begin(); it_end = dcs.end(); - double chi2p = 0.; + double chi2p = 0.; DCOnTrackVec dcs_new; dcs_new.reserve(dcs.size()); for(int i=0; it!=it_end; ++it, ++i ){ const DriftCircle* ds = & dcs[i]; - if(std::abs(ds->r()-ds->rot()->driftRadius())>m_dRTol && m_trace) ATH_MSG_DEBUG("Different radii on dc " << ds->r() << " rot " << ds->rot()->driftRadius()); + if(std::abs(ds->r()-ds->rot()->driftRadius())>m_dRTol && m_trace) ATH_MSG_DEBUG("Different radii on dc " << ds->r() << " rot " << ds->rot()->driftRadius()); DriftCircle dc_keep(ds->position(), ds->rot()->driftRadius(), ds->dr(), ds->drPrecise(), ds->state(), ds->id(), ds->index(),ds->rot() ); DCOnTrack dc_new(dc_keep, 0., 0.); dc_new.state(dcs[i].state()); @@ -556,15 +550,15 @@ namespace TrkDriftCircleMath { double tLow = rtInfo->rt()->tLower(); if(t<tLow) chi2p += (t-tLow)*(t-tLow)*0.1; if(t>tUp) chi2p += (t-tUp)*(t-tUp)*0.1; - } + } } if(m_trace&&chi2p>0) ATH_MSG_DEBUG(" Rejected weak topology Chi2 penalty " << chi2p); - bool oldrefit = DCSLFitter::fit( result, line, dcs_new, selection ); + bool oldrefit = DCSLFitter::fit( result, line, dcs_new, selection ); chi2p += result.chi2(); // add chi2 penalty for too large or too small driftTimes t < 0 or t> t upper result.set( chi2p, result.ndof(), result.dtheta(), result.dy0() ); return oldrefit; - } // end rejection of weak topologies + } // end rejection of weak topologies if(m_trace) ATH_MSG_DEBUG("positive radii ML1 " << nml1p << " ML2 " << nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << used << " t0 Error " << t0Error); @@ -572,20 +566,26 @@ namespace TrkDriftCircleMath { // starting point Double_t variable[3] = {theta,d,0}; // if t0Seed value from outside use this - if(t0Seed > -999.) variable[2] = t0Seed; + if(t0Seed > -999.) variable[2] = t0Seed; + + ROOT::Minuit2::Minuit2Minimizer minimum("algoName"); + minimum.SetMaxFunctionCalls(1000000); + minimum.SetTolerance(0.001); + minimum.SetPrintLevel(-1); + if(msgLvl(MSG::VERBOSE)) minimum.SetPrintLevel(1); - minimum->SetVariable(0,"a",variable[0], step[0]); - minimum->SetVariable(1,"b",variable[1], step[1]); - minimum->SetVariable(2,"t0",variable[2], step[2]); + minimum.SetVariable(0,"a",variable[0], step[0]); + minimum.SetVariable(1,"b",variable[1], step[1]); + minimum.SetVariable(2,"t0",variable[2], step[2]); - minimum->SetFunction(minFunct); + minimum.SetFunction(minFunct); // do the minimization - minimum->Minimize(); + minimum.Minimize(); - const double *results = minimum->X(); - const double *errors = minimum->Errors(); - ATH_MSG_DEBUG("Minimum: f(" << results[0] << "+-" << errors[0] << "," << results[1]<< "+-" << errors[1]<< "," << results[2] << "+-" << errors[2]<< "): " << minimum->MinValue()); + const double *results = minimum.X(); + const double *errors = minimum.Errors(); + ATH_MSG_DEBUG("Minimum: f(" << results[0] << "+-" << errors[0] << "," << results[1]<< "+-" << errors[1]<< "," << results[2] << "+-" << errors[2]<< "): " << minimum.MinValue()); ++m_npassedMinuitFit; @@ -609,7 +609,7 @@ namespace TrkDriftCircleMath { double t0=results[2]; double t0Err=errors[2]; double dy0 = cosin * bErr - b * sinus * aErr; - + double del_t; del_t = std::abs(rtpointers[0]->radius((t0+t0Err)) - rtpointers[0]->radius(t0)) ; @@ -623,7 +623,7 @@ namespace TrkDriftCircleMath { msg() << MSG::DEBUG <<"COVAR "; for(int it1=0; it1<3; it1++) { for(int it2=0; it2<3; it2++) { - msg() << MSG::DEBUG <<minimum->CovMatrix(it1,it2)<<" "; + msg() << MSG::DEBUG <<minimum.CovMatrix(it1,it2)<<" "; } msg() << MSG::DEBUG << endmsg; } @@ -643,34 +643,34 @@ namespace TrkDriftCircleMath { it = dcs_keep.begin(); it_end = dcs_keep.end(); ATH_MSG_DEBUG("------NEW HITS------"); - + for(int i=0; it!=it_end; ++it, ++i ){ double rad, drad; - + double uppercut = rtpointers[i]->tUpper(); double lowercut = rtpointers[i]->tLower(); rad = rtpointers[i]->radius(t[i]-t0); - if(t[i]-t0<lowercut) rad = rtpointers[i]->radius(lowercut); - if(t[i]-t0>uppercut) rad = rtpointers[i]->radius(uppercut); + if(t[i]-t0<lowercut) rad = rtpointers[i]->radius(lowercut); + if(t[i]-t0>uppercut) rad = rtpointers[i]->radius(uppercut); if (w[i]==0) { ATH_MSG_WARNING("w[i]==0, continuing"); continue; } - drad = 1.0/std::sqrt(w[i]) ; - + drad = 1.0/std::sqrt(w[i]) ; + yl = (y[i] - tana*z[i] - b); if(m_trace) { ATH_MSG_DEBUG("i "<<i<<" "); } - + double dth = -(sinus*y[i] + cosin*z[i])*dtheta; double residuals = std::abs(yl)/std::sqrt(1+tana*tana) - rad; if(m_trace) { ATH_MSG_DEBUG(" dth "<<dth<<" dy0 "<<dy0<<" del_t "<<del_t); } - + double errorResiduals = std::hypot(dth, dy0, del_t); - + // derivatives of the residual 'R' double deriv[3]; // del R/del theta @@ -679,13 +679,13 @@ namespace TrkDriftCircleMath { // del R / del b deriv[1] = sign(dd) * cosin ; // del R / del t0 - + deriv[2] = -1* rtpointers[i]->driftvelocity(t[i]-t0); - + double covsq=0; for(int rr=0; rr<3; rr++) { for(int cc=0; cc<3; cc++) { - covsq += deriv[rr]*minimum->CovMatrix(rr,cc)* deriv[cc]; + covsq += deriv[rr]*minimum.CovMatrix(rr,cc)* deriv[cc]; } } if(m_trace) { @@ -693,21 +693,21 @@ namespace TrkDriftCircleMath { if( covsq < 0. ){ for(int rr=0; rr<3; rr++) { for(int cc=0; cc<3; cc++) { - double dot = deriv[rr]*minimum->CovMatrix(rr,cc)* deriv[cc]; - ATH_MSG_DEBUG(" adding term " << dot << " dev1 " << deriv[rr] << " cov " << minimum->CovMatrix(rr,cc) << " dev2 " << deriv[cc]); + double dot = deriv[rr]*minimum.CovMatrix(rr,cc)* deriv[cc]; + ATH_MSG_DEBUG(" adding term " << dot << " dev1 " << deriv[rr] << " cov " << minimum.CovMatrix(rr,cc) << " dev2 " << deriv[cc]); } } } } covsq = covsq > 0. ? std::sqrt(covsq) : 0.; const DriftCircle* ds = & dcs_keep[i]; - if (m_propagateErrors) drad = dr[i]; + if (m_propagateErrors) drad = dr[i]; DriftCircle dc_newrad(dcs_keep[i].position(), rad, drad, ds->state(), dcs_keep[i].id(), dcs_keep[i].index(),ds->rot() ); DCOnTrack dc_new(dc_newrad, residuals, covsq); dc_new.state(dcs_keep[i].state()); - + if(m_trace) ATH_MSG_DEBUG("T0 Segment hit res "<<residuals<<" eres "<<errorResiduals<<" covsq "<<covsq<<" ri " << r[i]<<" ro "<<rad<<" drad "<<drad << " sel "<<selection[i]<< " inv error " << w[i]); - + if( selection[i] == 0 ) { ++nhits; if (!m_propagateErrors) { @@ -717,15 +717,15 @@ namespace TrkDriftCircleMath { } if(m_trace) ATH_MSG_DEBUG("T0 Segment hit res "<<residuals<<" eres "<<errorResiduals<<" covsq "<<covsq<<" ri " << r[i]<<" radius after t0 "<<rad<<" radius error "<< drad << " original error " << dr[i]); // Put chi2 penalty for drift times outside window - if (t[i]-t0> uppercut ) { // too large + if (t[i]-t0> uppercut ) { // too large chi2 += (t[i]-t0-uppercut)* (t[i]-t0-uppercut)*0.1; }else if (t[i]-t0 < lowercut ) {// too small chi2 += ((t[i]-t0-lowercut)*(t[i]-t0-lowercut))*0.1; } - } + } result.dcs().push_back( dc_new ); } - + double oldshift; oldshift = result.t0Shift(); if(m_trace) ATH_MSG_DEBUG("end fit old "<<oldshift<< " new " <<t0); @@ -750,6 +750,5 @@ namespace TrkDriftCircleMath { } return true; } - -} +}