diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt
index a10513d9f89a01e5155cdc43753252d87fdabd1b..4dc59bec52a2862bd5b72822e44ab5fd365fa892 100644
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/CMakeLists.txt
@@ -5,21 +5,6 @@
# Declare the package name:
atlas_subdir( MdtSegmentT0Fitter )
-# Declare the package's dependencies:
-atlas_depends_on_subdirs( PUBLIC
- Control/AthenaBaseComps
- GaudiKernel
- MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibSvc
- MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerInterfaces
- Tracking/TrkUtilityPackages/TrkDriftCircleMath
- PRIVATE
- MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibData
- MuonSpectrometer/MuonCalib/MuonCalibTools
- MuonSpectrometer/MuonDetDescr/MuonReadoutGeometry
- MuonSpectrometer/MuonIdHelpers
- MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonPrepRawData
- MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonRIO_OnTrack )
-
# External dependencies:
find_package( ROOT COMPONENTS Minuit Core Tree MathCore Hist RIO pthread MathMore Minuit2 Matrix Physics HistPainter Rint )
@@ -32,4 +17,3 @@ atlas_add_component( MdtSegmentT0Fitter
# Install files from the package:
atlas_install_headers( MdtSegmentT0Fitter )
-
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/MdtSegmentT0Fitter/MdtSegmentT0Fitter.h b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/MdtSegmentT0Fitter/MdtSegmentT0Fitter.h
index a2813ee3795f7cca67216b47868be8e25f8d92b3..76b0e27f529047c92f0a010f090e5b7a4a9624b9 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/MdtSegmentT0Fitter/MdtSegmentT0Fitter.h
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/MdtSegmentT0Fitter/MdtSegmentT0Fitter.h
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
// Fitting for in situ calibration of segments
@@ -19,9 +19,7 @@
#include <vector>
#include <memory>
-class IIdToFixedIdTool;
-class MdtCalibrationDbTool;
-class TMinuit;
+#include "TMinuit.h"
namespace TrkDriftCircleMath {
@@ -29,13 +27,12 @@ namespace TrkDriftCircleMath {
public:
MdtSegmentT0Fitter(const std::string&,const std::string&,const IInterface*);
- virtual ~MdtSegmentT0Fitter ();
+ virtual ~MdtSegmentT0Fitter()=default;
virtual StatusCode initialize() override;
- virtual StatusCode finalize () override;
+ virtual StatusCode finalize() override;
virtual bool fit( Segment& result, const Line& line, const DCOnTrackVec& dcs, double t0Seed ) const override;
virtual bool fit( Segment& result, const Line& line, const DCOnTrackVec& dcs, const HitSelection& selection, double t0Seed ) const override;
-
virtual const DCSLFitter* getFitter() const override { return this; }
@@ -60,23 +57,21 @@ namespace TrkDriftCircleMath {
};
private:
- bool m_trace; // debug - traces operation
- bool m_dumpToFile; // debug - dumps some performance info
- bool m_dumpNoFit; // debug - print hit info where fit doesn't run
-
- bool m_useInternalRT; // whether to use an internal RT function or the one from Calibration Service
- ToolHandle<MdtCalibrationDbTool> m_calibrationDbTool;
-
- int m_minHits; // minimum number of selected hits for t0 fit. Otherwise use default
-
- bool m_constrainShifts; // whether to constrain t0 shifts to a 50 ns window
- double m_constrainT0Error; // t0 error that is used in the constraint
- bool m_rejectWeakTopologies; // Reject topolgies that do not have at least one +- combination in one Multilayer
- bool m_scaleErrors; //!< rescale errors in fit
-
- bool m_propagateErrors; //!< propagate errors
-
- TMinuit* m_minuit;
+ ToolHandle<MdtCalibrationDbTool> m_calibrationDbTool{this,"CalibrationDbTool","MdtCalibrationDbTool"};
+
+ Gaudi::Property<bool> m_trace{this,"TraceOperation",false,"debug - traces operation"};
+ Gaudi::Property<bool> m_dumpToFile{this,"DumpToFile",false,"debug - dumps some performance info"};
+ Gaudi::Property<bool> m_dumpNoFit{this,"DumpNoFit",false,"debug - print hit info where fit does not run"};
+ Gaudi::Property<bool> m_useInternalRT{this,"UseInternalRT",false,"whether to use an internal RT function or the one from Calibration Service"};
+ Gaudi::Property<bool> m_constrainShifts{this,"ConstrainShifts",false,"whether to constrain t0 shifts to a 50 ns window"};
+ Gaudi::Property<bool> m_rejectWeakTopologies{this,"RejectWeakTopologies",true,"reject topolgies that do not have at least one +- combination in one multilayer"};
+ Gaudi::Property<bool> m_scaleErrors{this,"RescaleErrors",true,"rescale errors in fit"};
+ Gaudi::Property<bool> m_propagateErrors{this,"PropagateErrors",true,"propagate errors"};
+ Gaudi::Property<int> m_minHits{this,"MinimumHits",4,"minimum number of selected hits for t0 fit. Otherwise use default"};
+ Gaudi::Property<double> m_constrainT0Error{this,"ConstrainT0Error",10,"t0 error that is used in the constraint"};
+ Gaudi::Property<float> m_dRTol{this,"dRTolerance",0.1};
+
+ std::unique_ptr<TMinuit> m_minuit;
// counters
mutable std::atomic_uint m_ntotalCalls;
@@ -86,7 +81,6 @@ namespace TrkDriftCircleMath {
mutable std::atomic_uint m_npassedMinHits;
mutable std::atomic_uint m_npassedMinuitFit;
- std::unique_ptr<MdtSegmentT0FcnData> m_fcnData; //!< Struct to hold data to pass to/from TMinuit fit function
};
inline bool MdtSegmentT0Fitter::fit( Segment& result, const Line& line, const DCOnTrackVec& dcs, double t0Seed ) const {
@@ -95,6 +89,4 @@ namespace TrkDriftCircleMath {
}
}
-
-
#endif
diff --git a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx
index aebf4ebb2e41f02934b3b8751e75a1f8d43163b7..60f43aac7fe2133b773a84ddf3e0b88304571d51 100755
--- a/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx
+++ b/MuonSpectrometer/MuonReconstruction/MuonSegmentMakers/MuonSegmentMakerTools/MdtSegmentT0Fitter/src/MdtSegmentT0Fitter.cxx
@@ -8,8 +8,6 @@
#include "MuonReadoutGeometry/MdtReadoutElement.h"
#include "MuonReadoutGeometry/MuonDetectorManager.h"
-#include "MdtCalibSvc/MdtCalibrationDbTool.h"
-#include "MuonCalibTools/IdToFixedIdTool.h"
#include "MdtCalibData/IRtRelation.h"
#include "MdtCalibData/IRtResolution.h"
#include "MdtCalibData/MdtRtRelation.h"
@@ -22,10 +20,6 @@
#include <atomic>
#include <mutex>
-#include "TMinuit.h"
-// tdc count bin size -- seems to be unused -> commenting for the moment as it would be wrong
-// #define TDCBINSIZE 0.78125 //1tdc=0.78125ns; 1tdc=0.1953125ns for BMG!
-
// number of fit parameters
#define NUMPAR 3
@@ -52,42 +46,19 @@ namespace TrkDriftCircleMath {
MdtSegmentT0Fitter::MdtSegmentT0Fitter(const std::string& ty,const std::string& na,const IInterface* pa)
: AthAlgTool(ty,na,pa),
DCSLFitter(),
- m_calibrationDbTool("MdtCalibrationDbTool",this),
m_ntotalCalls(0),
m_npassedNHits(0),
m_npassedSelectionConsistency(0),
m_npassedNSelectedHits(0),
m_npassedMinHits(0),
- m_npassedMinuitFit(0),
- m_fcnData(std::make_unique<MdtSegmentT0FcnData>())
- {
+ m_npassedMinuitFit(0) {
declareInterface <IDCSLFitProvider> (this);
- declareProperty("MinimumHits", m_minHits = 4);
- declareProperty("DumpToFile", m_dumpToFile = false);
- declareProperty("UseInternalRT", m_useInternalRT = false);
- declareProperty("TraceOperation", m_trace = false);
- declareProperty("DumpNoFit", m_dumpNoFit = false);
- declareProperty("ConstrainShifts", m_constrainShifts = false);
- declareProperty("ConstrainT0Error", m_constrainT0Error = 10.);
- declareProperty("RejectWeakTopologies", m_rejectWeakTopologies = true);
- declareProperty("RescaleErrors",m_scaleErrors = true );
- declareProperty("PropagateErrors",m_propagateErrors = true );
- declareProperty("CalibrationDbTool",m_calibrationDbTool);
}
-
- MdtSegmentT0Fitter::~MdtSegmentT0Fitter(){}
-
StatusCode MdtSegmentT0Fitter::initialize() {
- ATH_CHECK ( AthAlgTool::initialize() );
-
- m_fcnData->use_hardcoded = m_useInternalRT;
- m_fcnData->use_shift_constraint = m_constrainShifts;
- m_fcnData->constrainT0Error = m_constrainT0Error;
- //count = 0;
TMinuit* oldMinuit = gMinuit;
- m_minuit = new TMinuit(3);
+ m_minuit = std::make_unique<TMinuit>(3);
m_minuit->SetPrintLevel(-1); // -1: no output, 1: std output
if( msgLvl(MSG::VERBOSE) ) m_minuit->SetPrintLevel(1);
gMinuit = oldMinuit;
@@ -106,9 +77,8 @@ namespace TrkDriftCircleMath {
<< " sel. hits > 2 " << std::setw(10) << m_npassedNSelectedHits << " " << scaleFactor*m_npassedNSelectedHits << "\n"
<< " Hits > min hits " << std::setw(10) << m_npassedMinHits << " " << scaleFactor*m_npassedMinHits << "\n"
<< " Passed Fit " << std::setw(10) << m_npassedMinuitFit << " " << scaleFactor*m_npassedMinuitFit );
- if( gMinuit == m_minuit ) gMinuit = 0;
- delete m_minuit;
- return AthAlgTool::finalize();
+ if(gMinuit == m_minuit.get()) gMinuit = nullptr;
+ return StatusCode::SUCCESS;
}
@@ -119,7 +89,7 @@ namespace TrkDriftCircleMath {
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(double tin) {
+ double rcorr(const double tin) {
double rc;
if(tin < 16.) {
rc = -50.;
@@ -136,7 +106,7 @@ namespace TrkDriftCircleMath {
return rc;
}
- double t2r(double tin) {
+ double t2r(const double tin) {
if(tin < 0.) return 0.;
if(tin > MAX_DRIFT) return 20.;
@@ -152,7 +122,7 @@ namespace TrkDriftCircleMath {
}
/// derivatives of RT function, use to get errors
- double rcorrprime(double tin) {
+ double rcorrprime(const double tin) {
double rc;
if(tin < 16.) {
rc = 0.;
@@ -169,7 +139,7 @@ namespace TrkDriftCircleMath {
return rc;
}
- double t2rprime(double tin) {
+ double t2rprime(const double tin) {
if(tin < 0.) return 0.;
if(tin > MAX_DRIFT) return 20.;
@@ -187,7 +157,7 @@ namespace TrkDriftCircleMath {
/// 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(double r) {
+ double r2t(const double r) {
double ta = 0;
double tb = MAX_DRIFT;
if(r<t2r(ta) ) {
@@ -200,7 +170,7 @@ namespace TrkDriftCircleMath {
while (ta <= tb) {
double tm = (ta + tb) / 2; // compute mid point.
double rtm = t2r(tm);
- if(fabs(rtm - r) < 0.001 ) {
+ if(std::abs(rtm - r) < 0.001 ) {
return tm;
}
else if (r > rtm) {
@@ -216,11 +186,8 @@ namespace TrkDriftCircleMath {
return -1; // failed to find key
}
- double r2t_ext( std::vector<const MuonCalib::IRtRelation*> *rtpointers, double r, int i) {
+ double r2t_ext(std::vector<const MuonCalib::IRtRelation*> *rtpointers, double r, int i) {
const MuonCalib::IRtRelation* rtrel = rtpointers->at(i);
- // double ta = 0;
- // double tb =
-
double ta = rtrel->tLower();
double tb = rtrel->tUpper();
if(r<rtrel->radius(ta) ) {
@@ -233,7 +200,7 @@ namespace TrkDriftCircleMath {
while (ta <= tb) {
double tm = (ta + tb) / 2; // compute mid point.
double rtm = rtrel->radius(tm);
- if(fabs(rtm - r) < 0.001 ) {
+ if(std::abs(rtm - r) < 0.001 ) {
return tm;
}
else if (r > rtm) {
@@ -259,14 +226,8 @@ namespace TrkDriftCircleMath {
double b = par[1];
double t0 = par[2];
- double cosin = cos(ang);
- double sinus = sin(ang);
-// if ( sinus < 0.0 ) {
-// sinus = -sinus;
-// cosin = -cosin;
-// } else if ( sinus == 0.0 && cosin < 0.0 ) {
-// cosin = -cosin;
-// }
+ double cosin = std::cos(ang);
+ double sinus = std::sin(ang);
fval = 0.;
// Add t0 constraint
@@ -281,7 +242,7 @@ namespace TrkDriftCircleMath {
z = g_fcnData->data[i].z;
y = g_fcnData->data[i].y;
w = g_fcnData->data[i].w;
- dist = fabs(b*cosin + z*sinus - y*cosin); // same thing as fabs(a*z - y + b)/sqrt(1. + a*a);
+ dist = std::abs(b*cosin + z*sinus - y*cosin); // same thing as fabs(a*z - y + b)/sqrt(1. + a*a);
double uppercut = g_fcnData->use_hardcoded ? TUBE_TIME : g_fcnData->data[i].rt->tUpper();
double lowercut = g_fcnData->use_hardcoded ? 0 : g_fcnData->data[i].rt->tLower();
// Penalty for t<lowercut and t >uppercut
@@ -315,12 +276,19 @@ namespace TrkDriftCircleMath {
bool MdtSegmentT0Fitter::fit( Segment& result, const Line& line, const DCOnTrackVec& dcs, const HitSelection& selection, double t0Seed ) const {
++m_ntotalCalls;
- if(m_trace) ATH_MSG_DEBUG( "New seg: " );
+ if(m_trace) ATH_MSG_DEBUG("New seg: ");
const DCOnTrackVec& dcs_keep = dcs;
unsigned int N = dcs_keep.size();
- m_fcnData->used=0;
+
+ std::unique_ptr<MdtSegmentT0FcnData> fcnData = std::make_unique<MdtSegmentT0FcnData>();
+ fcnData->use_hardcoded = m_useInternalRT;
+ fcnData->use_shift_constraint = m_constrainShifts;
+ fcnData->constrainT0Error = m_constrainT0Error;
+
+
+ fcnData->used=0;
result.setT0Shift(-99999,-99999);
if(N<2) {
@@ -328,21 +296,20 @@ namespace TrkDriftCircleMath {
}
++m_npassedNHits;
if( selection.size() != N ) {
- ATH_MSG_ERROR( "MdtSegmentT0Fitter.cxx:fit with t0 <bad HitSelection>" );
+ ATH_MSG_ERROR("MdtSegmentT0Fitter.cxx:fit with t0 <bad HitSelection>");
return false;
}
++m_npassedSelectionConsistency;
for(unsigned int i=0;i<N;++i){
- if( selection[i] == 0 ) ++(m_fcnData->used);
- // if(m_trace) *m_log << MSG::DEBUG << " selection flag " << selection[i] << endmsg;
+ if( selection[i] == 0 ) ++(fcnData->used);
}
- if(m_fcnData->used < 2){
- if(m_trace) ATH_MSG_DEBUG( "TOO FEW HITS SELECTED" );
+ if(fcnData->used < 2){
+ if(m_trace) ATH_MSG_DEBUG("TOO FEW HITS SELECTED");
return false;
}
++m_npassedNSelectedHits;
- if(m_fcnData->used < m_minHits) {
- if(m_trace) ATH_MSG_DEBUG( "FEWER THAN Minimum HITS N " << m_minHits << " total hits " <<N<<" used " << m_fcnData->used );
+ if(fcnData->used < m_minHits) {
+ if(m_trace) ATH_MSG_DEBUG("FEWER THAN Minimum HITS N " << m_minHits << " total hits " <<N<<" used " << fcnData->used);
//
// Copy driftcircles and reset the drift radii as they might have been overwritten
@@ -357,7 +324,7 @@ namespace TrkDriftCircleMath {
double chi2p = 0.;
for(int i=0; it!=it_end; ++it, ++i ){
const DriftCircle* ds = & dcs[i];
- if(fabs(ds->r()-ds->rot()->driftRadius())>0.1&&m_trace) std::cout << " Different radii on dc " << ds->r() << " rot " << ds->rot()->driftRadius() << std::endl;
+ 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());
@@ -371,58 +338,60 @@ namespace TrkDriftCircleMath {
if(t>tUp) chi2p += (t-tUp)*(t-tUp)*0.1;
}
}
- if(m_trace&&chi2p>0) std::cout << " NO Minuit Fit TOO few hits Chi2 penalty " << chi2p << std::endl;
+ 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 );
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() );
+ result.set(chi2p, result.ndof(), result.dtheta(), result.dy0());
int iok = 0;
if(oldrefit) iok = 1;
- if(m_trace) std::cout << " chi2 total " << result.chi2() << " angle " << result.line().phi() << " y0 " << result.line().y0() << " nhits "<< selection.size() << " refit ok " << iok << std::endl;
+ if(m_trace) ATH_MSG_DEBUG(" chi2 total " << result.chi2() << " angle " << result.line().phi() << " y0 " << result.line().y0() << " nhits "<< selection.size() << " refit ok " << iok);
return oldrefit;
} else {
- if(m_trace) ATH_MSG_DEBUG( "FITTING FOR T0 N "<<N<<" used " << m_fcnData->used );
+ if(m_trace) ATH_MSG_DEBUG("FITTING FOR T0 N "<<N<<" used " << fcnData->used);
}
++m_npassedMinHits;
-
-
- if (m_trace) ATH_MSG_DEBUG(" in MdtSegmentT0Fitter::fit with N dcs "<< N << " hit selection size " << selection.size() );
- if(m_trace) ATH_MSG_DEBUG("in fit "<<result.hasT0Shift()<< " " <<result.t0Shift() );
+ 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, Yc=0, S=0, Sz=0, Sy=0;
+ double Zc(0);
+ double Yc(0);
+ double S(0);
+ double Sz(0);
+ double Sy(0);
std::vector<double> y(N);
std::vector<double> z(N);
std::vector<double> w(N);
std::vector<double> r(N);
std::vector<double> dr(N);
std::vector<double> t(N);
- // std::vector<double> driftt(N);
std::vector<const MuonCalib::IRtRelation*> rtpointers(N);
// allocate memory for data
- // allocate memory for data
- if( m_fcnData->data.capacity() != 50 ) m_fcnData->data.reserve(50);
- m_fcnData->data.resize(m_fcnData->used);
+ if(fcnData->data.capacity() != 50) fcnData->data.reserve(50);
+ fcnData->data.resize(fcnData->used);
{
DCOnTrackVec::const_iterator it = dcs_keep.begin();
DCOnTrackVec::const_iterator it_end = dcs_keep.end();
for(int ii=0 ;it!=it_end; ++it, ++ii ){
const Muon::MdtDriftCircleOnTrack *roto = it->rot();
if (!roto) {
- ATH_MSG_ERROR(" MdtSegmentT0Fitter: NO DC ROT pointer found " );
+ ATH_MSG_ERROR("MdtSegmentT0Fitter: NO DC ROT pointer found");
return false;
}
- ATH_MSG_DEBUG( "hit # "<<ii );
+ ATH_MSG_DEBUG("hit # "<<ii );
y[ii] = it->y();
z[ii] = it->x();
- r[ii] = std::abs( roto->driftRadius() );
+ r[ii] = std::abs(roto->driftRadius());
dr[ii] = it->dr();
const Muon::MdtPrepData *peerd;
- peerd = dynamic_cast<const Muon::MdtPrepData*>( roto->prepRawData() );
+ peerd = dynamic_cast<const Muon::MdtPrepData*>(roto->prepRawData());
if(!peerd) {
- ATH_MSG_ERROR("MdtSegmentT0Fitter: Can't convert to MdtPrepData* !! Not fitting for t0" );
+ ATH_MSG_ERROR("MdtSegmentT0Fitter: Can't convert to MdtPrepData* !! Not fitting for t0");
return false;
}
Identifier id = roto->identify();
@@ -437,17 +406,17 @@ namespace TrkDriftCircleMath {
w[ii]*=w[ii];
if(r[ii]<0){
r[ii] = 0.;
- if(m_trace) ATH_MSG_DEBUG( "MdtSegmentT0Fitter (using times) ERROR: <Negative r> " << r[ii] );
+ if(m_trace) ATH_MSG_DEBUG("MdtSegmentT0Fitter (using times) ERROR: <Negative r> " << r[ii]);
}
- 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(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];
Sy+= w[ii]*y[ii];
if(r[ii] > MAX_RAD ) {
- if(m_trace) ATH_MSG_DEBUG("Radius is too big" );
+ if(m_trace) ATH_MSG_DEBUG("Radius is too big");
}
}
}
@@ -456,7 +425,7 @@ namespace TrkDriftCircleMath {
Zc = Sz*inv_S;
Yc = Sy*inv_S;
- if(m_trace) ATH_MSG_DEBUG( "Yc " << Yc << " Zc " << Zc );
+ 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) {
@@ -468,8 +437,8 @@ namespace TrkDriftCircleMath {
DCOnTrackVec::const_iterator it = dcs_keep.begin();
DCOnTrackVec::const_iterator it_end = dcs_keep.end();
- m_fcnData->t_lo = 1e10;
- m_fcnData->t_hi = -1e10;
+ fcnData->t_lo = 1e10;
+ fcnData->t_hi = -1e10;
// replicate for the case where the external rt is used...
// each hit has an rt function with some range...we want to fit such that
// tlower_i < ti - t0 < tupper_i
@@ -484,13 +453,13 @@ namespace TrkDriftCircleMath {
for(int ii=0 ;it!=it_end; ++it, ++ii ){
if( selection[ii] == 0 ) {
- m_fcnData->data[selcount].z = z[ii];
- m_fcnData->data[selcount].y = y[ii];
- m_fcnData->data[selcount].r = r[ii];
- m_fcnData->data[selcount].w = w[ii];
- m_fcnData->data[selcount].rt = rtpointers[ii];
+ fcnData->data[selcount].z = z[ii];
+ fcnData->data[selcount].y = y[ii];
+ fcnData->data[selcount].r = r[ii];
+ fcnData->data[selcount].w = w[ii];
+ fcnData->data[selcount].rt = rtpointers[ii];
double r2tval;
- if(m_fcnData->use_hardcoded) {
+ if(fcnData->use_hardcoded) {
r2tval = r2t(r[ii]);
} else {
r2tval = r2t_ext(&rtpointers, r[ii], ii) ;
@@ -507,7 +476,7 @@ namespace TrkDriftCircleMath {
<<" r "<<r[ii]
<<" t "<<t[ii]
<<" t0 "<<tee0;
- if(!m_fcnData->use_hardcoded) {
+ if(!fcnData->use_hardcoded) {
msg() << MSG::DEBUG <<" tLower "<<tl;
msg() << MSG::DEBUG <<" tUpper "<<th;
}
@@ -515,26 +484,26 @@ namespace TrkDriftCircleMath {
}
t0seed += tee0;
st0 += tee0*tee0;
- if(tee0 < min_t0 && fabs(r2tval) < R2TSPURIOUS) min_t0 = tee0;
+ if(tee0 < min_t0 && std::abs(r2tval) < R2TSPURIOUS) min_t0 = tee0;
- m_fcnData->data[selcount].t = t[ii];
- if(t[ii]< m_fcnData->t_lo) m_fcnData->t_lo = t[ii];
- if(t[ii] > m_fcnData->t_hi) m_fcnData->t_hi = t[ii];
+ fcnData->data[selcount].t = t[ii];
+ if(t[ii]< fcnData->t_lo) fcnData->t_lo = t[ii];
+ if(t[ii] > fcnData->t_hi) fcnData->t_hi = t[ii];
selcount++;
}
}
t0seed /= selcount;
st0 = st0/selcount - t0seed*t0seed;
- st0 = st0 > 0. ? sqrt(st0) : 0.;
- m_fcnData->t_hi -= MAX_DRIFT;
+ st0 = st0 > 0. ? std::sqrt(st0) : 0.;
+ fcnData->t_hi -= MAX_DRIFT;
- if(!m_fcnData->use_hardcoded) {
- m_fcnData->t_hi = max_tupper;
- m_fcnData->t_lo = min_tlower;
+ if(!fcnData->use_hardcoded) {
+ fcnData->t_hi = max_tupper;
+ fcnData->t_lo = min_tlower;
}
- if(m_trace) ATH_MSG_DEBUG("t_hi "<<m_fcnData->t_hi<<" t_lo "<<m_fcnData->t_lo<<" t0seed "<<t0seed<<" sigma "<<st0<< " min_t0 "<<min_t0 );
- if(m_fcnData->t_hi > m_fcnData->t_lo ) {
+ if(m_trace) ATH_MSG_DEBUG("t_hi "<<fcnData->t_hi<<" t_lo "<<fcnData->t_lo<<" t0seed "<<t0seed<<" sigma "<<st0<< " min_t0 "<<min_t0);
+ if(fcnData->t_hi > fcnData->t_lo ) {
if(m_dumpNoFit) {
std::ofstream gg;
gg.open("fitnotdone.txt", std::ios::out | std::ios::app);
@@ -566,8 +535,8 @@ namespace TrkDriftCircleMath {
// ************************* seed the parameters
double theta = line.phi();
- double cosin = cos(theta);
- double sinus = sin(theta);
+ double cosin = std::cos(theta);
+ double sinus = std::sin(theta);
if ( sinus < 0.0 ) {
sinus = -sinus;
@@ -575,15 +544,15 @@ namespace TrkDriftCircleMath {
} else if ( sinus == 0.0 && cosin < 0.0 ) {
cosin = -cosin;
}
- if(m_trace) ATH_MSG_DEBUG("before fit theta "<<theta<<" sinus "<<sinus<< " 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 );
+ 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
@@ -596,8 +565,7 @@ namespace TrkDriftCircleMath {
it = dcs_keep.begin();
for(int ii=0 ;it!=it_end; ++it, ++ii ){
if( selection[ii] != 0 ) continue;
- sdist = d*cosin + z[ii]*sinus - y[ii]*cosin; // same thing as fabs(a*z - y + b)/sqrt(1. + a*a);
-// if (r[ii] < 1) continue;
+ sdist = d*cosin + z[ii]*sinus - y[ii]*cosin; // same thing as |a*z - y + b|/sqrt(1. + a*a);
if(it->id().ml()==0&&sdist > 0) nml1p++;
if(it->id().ml()==0&&sdist < 0) nml1n++;
if(it->id().ml()==1&&sdist > 0) nml2p++;
@@ -605,12 +573,12 @@ namespace TrkDriftCircleMath {
}
// Define t0 constraint in Minuit
- m_fcnData->t0Error = STRONG_TOPO_T0ERROR;
- if (nml1p+nml2p < 2 || nml1n+nml2n < 2) m_fcnData->t0Error = WEAK_TOPO_T0ERROR;
+ fcnData->t0Error = STRONG_TOPO_T0ERROR;
+ if (nml1p+nml2p < 2 || nml1n+nml2n < 2) fcnData->t0Error = WEAK_TOPO_T0ERROR;
// Reject topologies where in one of the Multilayers no +- combination is present
if((nml1p<1||nml1n<1)&&(nml2p<1||nml2n<1)&&m_rejectWeakTopologies) {
- if(m_trace) ATH_MSG_DEBUG(" Combination rejected for positive radii ML1 " << nml1p << " ML2 " << nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << m_fcnData->used << " t0 Error " << m_fcnData->t0Error );
+ if(m_trace) ATH_MSG_DEBUG("Combination rejected for positive radii ML1 " << nml1p << " ML2 " << nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << fcnData->used << " t0 Error " << fcnData->t0Error);
it = dcs.begin();
it_end = dcs.end();
double chi2p = 0.;
@@ -618,7 +586,7 @@ namespace TrkDriftCircleMath {
dcs_new.reserve(dcs.size());
for(int i=0; it!=it_end; ++it, ++i ){
const DriftCircle* ds = & dcs[i];
- if(fabs(ds->r()-ds->rot()->driftRadius())>0.1&&m_trace) std::cout << " Different radii on dc " << ds->r() << " rot " << ds->rot()->driftRadius() << std::endl;
+ 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());
@@ -632,7 +600,7 @@ namespace TrkDriftCircleMath {
if(t>tUp) chi2p += (t-tUp)*(t-tUp)*0.1;
}
}
- if(m_trace&&chi2p>0) ATH_MSG_DEBUG( " Rejected weak topology Chi2 penalty " << chi2p );
+ if(m_trace&&chi2p>0) ATH_MSG_DEBUG(" Rejected weak topology Chi2 penalty " << chi2p);
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
@@ -640,34 +608,7 @@ namespace TrkDriftCircleMath {
return oldrefit;
} // 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 " << m_fcnData->used << " t0 Error " << m_fcnData->t0Error );
-
-// m_minuit->BuildArrays(3);
- m_minuit->SetFCN(mdtSegmentT0Fcn);
- Double_t arglist[3];
-
- // Set printout level
- arglist[0] = -1;
- // Clear previous fit
- m_minuit->mncler();
- arglist[0] = -1;
-
- Double_t vstart[3];
- vstart[0] = theta;
- // x' = x - xc, y' = y - yc => y' = m x' + b + m xc - yc
- // and b = yc - m xc
- vstart[1] = d ;
- vstart[2] = 0 ;
-
- // if t0Seed value from outside use this
- if(t0Seed > -999.) vstart[2] = t0Seed;
-
- Double_t step[3] = {0.01 , 0.01 , 0.1 };
-
- if(m_trace) {
- ATH_MSG_DEBUG("\n____________INITIAL VALUES________________" );
- ATH_MSG_DEBUG("Theta " << theta << "("<<tan(theta)<<") d " << vstart[1] );
- }
+ if(m_trace) ATH_MSG_DEBUG("positive radii ML1 " << nml1p << " ML2 " << nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << fcnData->used << " t0 Error " << fcnData->t0Error);
int errFlag = 0;
int errFlag1 = 0;
@@ -675,21 +616,49 @@ namespace TrkDriftCircleMath {
int errFlag3 = 0;
int errFlag4 = 0;
int errFlag5 = 0;
- m_minuit->mnparm(0, "a", vstart[0], step[0], -0.5, 3.5,errFlag);
- m_minuit->mnparm(1, "b", vstart[1], step[1], 0., 0.,errFlag);
- // m_minuit->SetParameter(2, "t0", vstart[2], step[2], t_hi, t_lo,errFlag);
- m_minuit->mnparm(2, "t0", vstart[2], step[2], 0., 0.,errFlag);
-
- if (errFlag !=0 &&m_trace) {
- ATH_MSG_DEBUG( " ALARM with steering of Minuit variable " << errFlag );
- }
-
- m_minuit->FixParameter(0);
- m_minuit->FixParameter(1);
{
// The part where the fcn gets used, requires locking to protect fcnData from concurrent access
std::scoped_lock lock(g_fcnDataMutex);
- g_fcnData = m_fcnData.get();
+ g_fcnData = fcnData.get();
+
+ m_minuit->SetFCN(mdtSegmentT0Fcn);
+ Double_t arglist[3];
+
+ // Set printout level
+ arglist[0] = -1;
+ // Clear previous fit
+ m_minuit->mncler();
+ arglist[0] = -1;
+
+ Double_t vstart[3];
+ vstart[0] = theta;
+ // x' = x - xc, y' = y - yc => y' = m x' + b + m xc - yc
+ // and b = yc - m xc
+ vstart[1] = d ;
+ vstart[2] = 0 ;
+
+ // if t0Seed value from outside use this
+ if(t0Seed > -999.) vstart[2] = t0Seed;
+
+ Double_t step[3] = {0.01 , 0.01 , 0.1 };
+
+ if(m_trace) {
+ ATH_MSG_DEBUG("\n____________INITIAL VALUES________________" );
+ ATH_MSG_DEBUG("Theta " << theta << "("<<std::tan(theta)<<") d " << vstart[1]);
+ }
+
+
+ m_minuit->mnparm(0, "a", vstart[0], step[0], -0.5, 3.5,errFlag);
+ m_minuit->mnparm(1, "b", vstart[1], step[1], 0., 0.,errFlag);
+ m_minuit->mnparm(2, "t0", vstart[2], step[2], 0., 0.,errFlag);
+
+ if (errFlag !=0 &&m_trace) {
+ ATH_MSG_DEBUG("ALARM with steering of Minuit variable " << errFlag);
+ }
+
+ m_minuit->FixParameter(0);
+ m_minuit->FixParameter(1);
+
m_minuit->mnexcm("MINIMIZE", arglist, 0,errFlag1);
m_minuit->Release(0);
@@ -703,21 +672,23 @@ namespace TrkDriftCircleMath {
if(errFlag5!=0&&errFlag4==0) {
m_minuit->mnexcm("MINIMIZE", arglist, 0,errFlag4);
if (errFlag4 == 0 &&m_trace) {
- ATH_MSG_DEBUG( " ALARM Fall back to MINIMIZE " << errFlag4 );
+ ATH_MSG_DEBUG(" ALARM Fall back to MINIMIZE " << errFlag4);
}
}
g_fcnData = nullptr;
- }
+ } // end of scoped_lock
// Get the chisquare and errors
- double chisq;
- double edm;
- double errdef;
- int npari, nparx, icstat;
+ double chisq(0);
+ double edm(0);
+ double errdef(0);
+ int npari(0);
+ int nparx(0);
+ int icstat(0);
m_minuit->mnstat(chisq, edm, errdef, npari, nparx, icstat);
int fitretval = errFlag5;
- if (npari<0 || nparx < 0 || chisq < 0) ATH_MSG_ERROR( " MdtSegmentT0Fitter MINUIT problem " << " chisq "<<chisq<<" npari "<<npari<<" nparx "<< nparx <<" fitretval "<< fitretval );
+ if (npari<0 || nparx < 0 || chisq < 0) ATH_MSG_ERROR("MdtSegmentT0Fitter MINUIT problem " << " chisq "<<chisq<<" npari "<<npari<<" nparx "<< nparx <<" fitretval "<< fitretval);
if (fitretval !=0 &&m_trace) {
ATH_MSG_DEBUG( " ALARM return value " << fitretval );
@@ -725,19 +696,17 @@ namespace TrkDriftCircleMath {
if(m_dumpToFile) {
std::ofstream ff;
ff.open("fitres.txt", std::ios::out | std::ios::app);
-// ff << fitretval<<" 0 "<<chisq<<" "<<t_hi<<" "<<t_lo<<" "<<t0seed<<" "<<st0<<" "<<min_t0<<" "<<m_minuit->GetParameter(2)<<" "<<m_minuit->GetParError(2)<< endmsg;
-// TMinut * tm = (TMinuit*) minuit->GetObjectFit();
ff<<npari<<" "<<nparx<<" "<<fitretval<<" "<<chisq<< " "<< std::endl;
ff.close();
}
- if(m_trace) ATH_MSG_DEBUG("chisq "<<chisq<<" npari "<<npari<<" nparx "<<nparx<<" fitretval "<<fitretval );
+ if(m_trace) ATH_MSG_DEBUG("chisq "<<chisq<<" npari "<<npari<<" nparx "<<nparx<<" fitretval "<<fitretval);
// Do standard DCSL fit if minuit is bad
// Do standard DCSL fit if all fits from minuit are bad
- if (errFlag5!=0&&errFlag4!=0&&errFlag3!=0&&errFlag2!=0&&errFlag1!=0) {
- if(m_trace)ATH_MSG_DEBUG( " ALARM Minimize fix " << errFlag1 << " ALARM minimize release " << errFlag2
+ if (errFlag5!=0&&errFlag4!=0&&errFlag3!=0&&errFlag2!=0&&errFlag1!=0) {
+ if(m_trace)ATH_MSG_DEBUG(" ALARM Minimize fix " << errFlag1 << " ALARM minimize release " << errFlag2
<< " ALARM migrad fix 1 " << errFlag3 << " ALARM minimize all free" << errFlag4 << " ALARM migrad all free "
- << errFlag5 );
+ << errFlag5);
DCOnTrackVec dcs_new;
dcs_new.reserve(dcs.size());
@@ -745,64 +714,53 @@ namespace TrkDriftCircleMath {
DCOnTrackVec::const_iterator it_end = dcs.end();
for(int i=0; it!=it_end; ++it, ++i ){
- const DriftCircle* ds = & dcs[i];
- 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());
- dcs_new.push_back( dc_new );
+ const DriftCircle* ds = & dcs[i];
+ 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());
+ dcs_new.push_back( dc_new );
}
-
bool oldrefit = DCSLFitter::fit( result, line , dcs_new, selection );
- //int iok = 0;
- //if(oldrefit) iok = 1;
- return oldrefit;
- }
+ return oldrefit;
+ }
++m_npassedMinuitFit;
// Get the fit values
- double aret(0),aErr(0);
+ double aret(0);
+ double aErr(0);
m_minuit->GetParameter(0,aret,aErr); // theta returned
double dtheta = aErr;
-
- double tana = tan(aret); // tangent of angle
-
+ double tana = std::tan(aret); // tangent of angle
double ang = aret; // between zero and pi
- cosin = cos(ang);
- sinus = sin(ang);
+ cosin = std::cos(ang);
+ sinus = std::sin(ang);
if ( sinus < 0.0 ) {
sinus = -sinus;
cosin = -cosin;
} else if ( sinus == 0.0 && cosin < 0.0 ) {
cosin = -cosin;
}
- ang = atan2(sinus, cosin);
+ ang = std::atan2(sinus, cosin);
- double b(0),bErr(0);
+ double b(0);
+ double bErr(0);
m_minuit->GetParameter(1,b,bErr); // intercept
-
-
- double t0(0),t0Err(0);
+ double t0(0);
+ double t0Err(0);
m_minuit->GetParameter(2,t0,t0Err);
-
-
- // double y0 = b * cosin ; // parallel distance
double dy0 = cosin * bErr - b * sinus * aErr;
double del_t;
- if(m_fcnData->use_hardcoded) del_t = fabs(t2r((t0+t0Err)) - t2r(t0));
- else del_t = fabs(rtpointers[0]->radius((t0+t0Err)) - rtpointers[0]->radius(t0)) ;
+ if(fcnData->use_hardcoded) del_t = std::abs(t2r((t0+t0Err)) - t2r(t0));
+ else del_t = std::abs(rtpointers[0]->radius((t0+t0Err)) - rtpointers[0]->radius(t0)) ;
if(m_trace) {
ATH_MSG_DEBUG("____________FINAL VALUES________________" );
- ATH_MSG_DEBUG("Values: a "<<tana<<" d "<<b * cosin <<" t0 "<<t0 );
- ATH_MSG_DEBUG("Errors: a "<<aErr<<" b "<<dy0 <<" t0 "<<t0Err );
+ ATH_MSG_DEBUG("Values: a "<<tana<<" d "<<b * cosin <<" t0 "<<t0);
+ ATH_MSG_DEBUG("Errors: a "<<aErr<<" b "<<dy0 <<" t0 "<<t0Err);
}
-
-
d = b * cosin;
-
-
double covar[3][3];
m_minuit->mnemat(&covar[0][0],3); // 3x3
if(m_trace) {
@@ -819,7 +777,7 @@ namespace TrkDriftCircleMath {
result.clusters().clear();
result.emptyTubes().clear();
- if(m_trace) ATH_MSG_DEBUG("after fit theta "<<ang<<" sinus "<<sinus<< " cosin "<< cosin );
+ if(m_trace) ATH_MSG_DEBUG("after fit theta "<<ang<<" sinus "<<sinus<< " cosin "<< cosin);
double chi2 = 0;
unsigned int nhits(0);
@@ -829,14 +787,14 @@ namespace TrkDriftCircleMath {
// calculate predicted hit positions from track parameters
it = dcs_keep.begin();
it_end = dcs_keep.end();
- ATH_MSG_DEBUG("------NEW HITS------" );
+ ATH_MSG_DEBUG("------NEW HITS------");
for(int i=0; it!=it_end; ++it, ++i ){
double rad, drad;
- double uppercut = m_fcnData->use_hardcoded ? TUBE_TIME : rtpointers[i]->tUpper();
- double lowercut = m_fcnData->use_hardcoded ? 0 : rtpointers[i]->tLower();
- if(m_fcnData->use_hardcoded) {
+ double uppercut = fcnData->use_hardcoded ? TUBE_TIME : rtpointers[i]->tUpper();
+ double lowercut = fcnData->use_hardcoded ? 0 : rtpointers[i]->tLower();
+ if(fcnData->use_hardcoded) {
rad = t2r(t[i]-t0);
if(t[i]-t0<lowercut) rad = t2r(lowercut);
if(t[i]-t0>uppercut) rad = t2r(uppercut);
@@ -845,22 +803,24 @@ namespace TrkDriftCircleMath {
if(t[i]-t0<lowercut) rad = rtpointers[i]->radius(lowercut);
if(t[i]-t0>uppercut) rad = rtpointers[i]->radius(uppercut);
}
-
- drad = 1.0/sqrt(w[i]) ;
+ if (w[i]==0) {
+ ATH_MSG_WARNING("w[i]==0, continuing");
+ continue;
+ }
+ 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 = fabs(yl)/sqrt((1+tana*tana)) - rad;
+ 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<<" " );
+ ATH_MSG_DEBUG(" dth "<<dth<<" dy0 "<<dy0<<" del_t "<<del_t);
}
- double errorResiduals = sqrt( dth*dth + dy0*dy0 + del_t*del_t);
+ double errorResiduals = std::hypot(dth, dy0, del_t);
// derivatives of the residual 'R'
double deriv[3];
@@ -871,7 +831,7 @@ namespace TrkDriftCircleMath {
deriv[1] = sign(dd) * cosin ;
// del R / del t0
- if(m_fcnData->use_hardcoded) deriv[2] = -1* t2rprime(t[i]-t0);
+ if(fcnData->use_hardcoded) deriv[2] = -1* t2rprime(t[i]-t0);
else deriv[2] = -1* rtpointers[i]->driftvelocity(t[i]-t0);
double covsq=0;
@@ -881,35 +841,33 @@ namespace TrkDriftCircleMath {
}
}
if(m_trace) {
- ATH_MSG_DEBUG( " covsquared " << covsq );
+ ATH_MSG_DEBUG(" covsquared " << covsq);
if( covsq < 0. ){
for(int rr=0; rr<3; rr++) {
for(int cc=0; cc<3; cc++) {
double dot = deriv[rr]*covar[rr][cc]* deriv[cc];
- ATH_MSG_DEBUG( " adding term " << dot << " dev1 " << deriv[rr] << " cov " << covar[rr][cc]
- << " dev2 " << deriv[cc] );
+ ATH_MSG_DEBUG(" adding term " << dot << " dev1 " << deriv[rr] << " cov " << covar[rr][cc] << " dev2 " << deriv[cc]);
}
}
}
}
- covsq = covsq > 0. ? sqrt(covsq) : 0.;
+ covsq = covsq > 0. ? std::sqrt(covsq) : 0.;
const DriftCircle* ds = & dcs_keep[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(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;
-// chi2 += residuals*residuals*hitwt;
if (!m_propagateErrors) {
chi2 += residuals*residuals*w[i];
} else {
chi2 += residuals*residuals/(drad*drad);
}
- if(m_trace) std::cout <<" T0 Segment hit res "<<residuals<<" eres "<<errorResiduals<<" covsq "<<covsq<<" ri " << r[i]<<" radius after t0 "<<rad<<" radius error "<< drad << " original error " << dr[i] << std::endl;
+ 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
chi2 += (t[i]-t0-uppercut)* (t[i]-t0-uppercut)*0.1;
@@ -922,7 +880,7 @@ namespace TrkDriftCircleMath {
double oldshift;
oldshift = result.t0Shift();
- if(m_trace) ATH_MSG_DEBUG("end fit old "<<oldshift<< " new " <<t0 );
+ if(m_trace) ATH_MSG_DEBUG("end fit old "<<oldshift<< " new " <<t0);
// Final Minuit Fit result
if(nhits==NUMPAR) {
nhits++;
@@ -932,32 +890,18 @@ namespace TrkDriftCircleMath {
result.line().set( LocPos( Zc - sinus*d, Yc + cosin*d ), ang );
if(t0==0.) t0=0.00001;
result.setT0Shift(t0,t0Err);
- if(m_trace) std::cout << "Minuit Fit complete: Chi2 " << chi2 << " angle " << result.line().phi() << " nhits "<< nhits << " t0result " << t0 << std::endl;
-
if(m_trace) {
- ATH_MSG_DEBUG( "Minuit Fit complete: Chi2 " << chi2 << " angle " << result.line().phi() << " nhits "<<nhits<<" numpar "<<NUMPAR << " per dof " << chi2/(nhits-NUMPAR) );
+ ATH_MSG_DEBUG("Minuit Fit complete: Chi2 " << chi2 << " angle " << result.line().phi() << " nhits "<< nhits << " t0result " << t0);
+ ATH_MSG_DEBUG("Minuit Fit complete: Chi2 " << chi2 << " angle " << result.line().phi() << " nhits "<<nhits<<" numpar "<<NUMPAR << " per dof " << chi2/(nhits-NUMPAR));
- ATH_MSG_DEBUG( "Fit complete: Chi2 " << chi2 <<" nhits "<<nhits<<" numpar "<<NUMPAR << " per dof " << chi2/(nhits-NUMPAR) );
+ ATH_MSG_DEBUG("Fit complete: Chi2 " << chi2 <<" nhits "<<nhits<<" numpar "<<NUMPAR << " per dof " << chi2/(nhits-NUMPAR));
if(chi2/(nhits-NUMPAR) > 5) {
- ATH_MSG_DEBUG( "_______NOT GOOD " );
+ ATH_MSG_DEBUG("_______NOT GOOD ");
}
- ATH_MSG_DEBUG("chi2 "<<chi2<<" per dof "<<chi2/(nhits-NUMPAR)<<" root chisq "<<chisq );
+ ATH_MSG_DEBUG("chi2 "<<chi2<<" per dof "<<chi2/(nhits-NUMPAR)<<" root chisq "<<chisq);
}
-
- // Do standard DCSL fit using the t0 updated results
- // if minuit MIGRAD failed to converge (and t0 fit converged)
-// if (errFlag5!=0&&errFlag4!=0) {
-// bool oldrefit = DCSLFitter::fit( result, result.line() , result.dcs(), selection );
-// int iok = 0;
-// if(oldrefit) iok = 1;
-// if(m_trace)*m_log << MSG::DEBUG << " ALARM Minimize fix " << errFlag1 << " ALARM minimize release " << errFlag2 << " ALARM migrad fix 1 " << errFlag3 << " ALARM minimize all free " << errFlag4 << " ALARM migrad all free " << errFlag5 << endmsg;
-// if(m_trace) std::cout << " ALARM Migrad and Minimize failed DCSL Fit Chi2 " << result.chi2() << " angle " << result.line().phi() << " y0 " << result.line().y0() << " nhits "<< selection.size() << " fitretval " << fitretval << " refit ok " << iok << std::endl;
-// return oldrefit;
-// }
-
-
return true;
}
-
}
+