TrkDriftCircleMath::MdtSegmentT0Fitter Class Reference

#include <MdtSegmentT0Fitter.h>

Inheritance diagram for TrkDriftCircleMath::MdtSegmentT0Fitter:

Collaboration diagram for TrkDriftCircleMath::MdtSegmentT0Fitter:

Public Member Functions
	MdtSegmentT0Fitter (const std::string &, const std::string &, const IInterface *)
virtual	~MdtSegmentT0Fitter ()=default
virtual StatusCode	initialize () 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
	fitter factory More...
ServiceHandle< StoreGateSvc > &	evtStore ()
	The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc. More...
const ServiceHandle< StoreGateSvc > &	evtStore () const
	The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc. More...
const ServiceHandle< StoreGateSvc > &	detStore () const
	The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc. More...
virtual StatusCode	sysInitialize () override
	Perform system initialization for an algorithm. More...
virtual StatusCode	sysStart () override
	Handle START transition. More...
virtual std::vector< Gaudi::DataHandle * >	inputHandles () const override
	Return this algorithm's input handles. More...
virtual std::vector< Gaudi::DataHandle * >	outputHandles () const override
	Return this algorithm's output handles. More...
Gaudi::Details::PropertyBase &	declareProperty (Gaudi::Property< T > &t)
Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleKey &hndl, const std::string &doc, const SG::VarHandleKeyType &)
	Declare a new Gaudi property. More...
Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleBase &hndl, const std::string &doc, const SG::VarHandleType &)
	Declare a new Gaudi property. More...
Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, SG::VarHandleKeyArray &hndArr, const std::string &doc, const SG::VarHandleKeyArrayType &)
Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, T &property, const std::string &doc, const SG::NotHandleType &)
	Declare a new Gaudi property. More...
Gaudi::Details::PropertyBase *	declareProperty (const std::string &name, T &property, const std::string &doc="none")
	Declare a new Gaudi property. More...
void	updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream &	msg () const
MsgStream &	msg (const MSG::Level lvl) const
bool	msgLvl (const MSG::Level lvl) const

Static Public Member Functions
static const InterfaceID &	interfaceID ()

Protected Member Functions
void	renounceArray (SG::VarHandleKeyArray &handlesArray)
	remove all handles from I/O resolution More...
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void >	renounce (T &h)
void	extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
	Add StoreName to extra input/output deps as needed. More...

Private Types
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey> More...
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyArrayType &)
	specialization for handling Gaudi::Property<SG::VarHandleKeyArray> More...
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleType &)
	specialization for handling Gaudi::Property<SG::VarHandleBase> More...
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T > &t, const SG::NotHandleType &)
	specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray> More...

Private Attributes
SG::ReadCondHandleKey< MuonCalib::MdtCalibDataContainer >	m_calibDbKey
ServiceHandle< Muon::IMuonIdHelperSvc >	m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}
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< float >	m_dRTol {this,"dRTolerance",0.1}
Gaudi::Property< bool >	m_floatDir
std::atomic_uint	m_ntotalCalls {0}
std::atomic_uint	m_npassedNHits {0}
std::atomic_uint	m_npassedSelectionConsistency {0}
std::atomic_uint	m_npassedNSelectedHits {0}
std::atomic_uint	m_npassedMinHits {0}
std::atomic_uint	m_npassedMinuitFit {0}
StoreGateSvc_t	m_evtStore
	Pointer to StoreGate (event store by default) More...
StoreGateSvc_t	m_detStore
	Pointer to StoreGate (detector store by default) More...
std::vector< SG::VarHandleKeyArray * >	m_vhka
bool	m_varHandleArraysDeclared

Detailed Description

Definition at line 25 of file MdtSegmentT0Fitter.h.

Member Typedef Documentation

StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

MdtSegmentT0Fitter()

TrkDriftCircleMath::MdtSegmentT0Fitter::MdtSegmentT0Fitter	(	const std::string &	ty,
		const std::string &	na,
		const IInterface *	pa
	)

Definition at line 168 of file MdtSegmentT0Fitter.cxx.

  : AthAlgTool(ty,na,pa),
    DCSLFitter() {
    declareInterface <IDCSLFitProvider> (this);
  }

~MdtSegmentT0Fitter()

virtual TrkDriftCircleMath::MdtSegmentT0Fitter::~MdtSegmentT0Fitter ( )

virtualdefault

Member Function Documentation

declareGaudiProperty() [1/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T > &  hndl, const SG::VarHandleKeyArrayType &    )

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKeyArray>

Definition at line 170 of file AthCommonDataStore.h.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
 
  }

declareGaudiProperty() [2/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T > &  hndl, const SG::VarHandleKeyType &    )

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
 
  }

declareGaudiProperty() [3/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T > &  hndl, const SG::VarHandleType &    )

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleBase>

Definition at line 184 of file AthCommonDataStore.h.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
  }

declareGaudiProperty() [4/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T > &  t, const SG::NotHandleType &    )

inlineprivateinherited

specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray>

Definition at line 199 of file AthCommonDataStore.h.

  {
    return PBASE::declareProperty(t);
  }

declareProperty() [1/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name, SG::VarHandleBase &  hndl, const std::string &  doc, const SG::VarHandleType &    )

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
hndl	Object holding the property value.
doc	Documentation string for the property.

This is the version for types that derive from SG::VarHandleBase. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 245 of file AthCommonDataStore.h.

  {
    this->declare(hndl.vhKey());
    hndl.vhKey().setOwner(this);
 
    return PBASE::declareProperty(name,hndl,doc);
  }

declareProperty() [2/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name, SG::VarHandleKey &  hndl, const std::string &  doc, const SG::VarHandleKeyType &    )

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
hndl	Object holding the property value.
doc	Documentation string for the property.

This is the version for types that derive from SG::VarHandleKey. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 221 of file AthCommonDataStore.h.

  {
    this->declare(hndl);
    hndl.setOwner(this);
 
    return PBASE::declareProperty(name,hndl,doc);
  }

declareProperty() [3/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name, SG::VarHandleKeyArray &  hndArr, const std::string &  doc, const SG::VarHandleKeyArrayType &    )

inlineinherited

Definition at line 259 of file AthCommonDataStore.h.

  {
 
    // std::ostringstream ost;
    // ost << Algorithm::name() << " VHKA declareProp: " << name 
    //     << " size: " << hndArr.keys().size() 
    //     << " mode: " << hndArr.mode() 
    //     << "  vhka size: " << m_vhka.size()
    //     << "\n";
    // debug() << ost.str() << endmsg;
 
    hndArr.setOwner(this);
    m_vhka.push_back(&hndArr);
 
    Gaudi::Details::PropertyBase* p =  PBASE::declareProperty(name, hndArr, doc);
    if (p != 0) {
      p->declareUpdateHandler(&AthCommonDataStore<PBASE>::updateVHKA, this);
    } else {
      ATH_MSG_ERROR("unable to call declareProperty on VarHandleKeyArray " 
                    << name);
    }
 
    return p;
 
  }

declareProperty() [4/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name, T &  property, const std::string &  doc, const SG::NotHandleType &    )

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
property	Object holding the property value.
doc	Documentation string for the property.

This is the generic version, for types that do not derive from SG::VarHandleKey. It just forwards to the base class version of declareProperty.

Definition at line 333 of file AthCommonDataStore.h.

  {
    return PBASE::declareProperty(name, property, doc);
  }

declareProperty() [5/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name, T &  property, const std::string &  doc = "none"  )

inlineinherited

Declare a new Gaudi property.

Parameters

name	Name of the property.
property	Object holding the property value.
doc	Documentation string for the property.

This dispatches to either the generic declareProperty or the one for VarHandle/Key/KeyArray.

Definition at line 352 of file AthCommonDataStore.h.

  {
    typedef typename SG::HandleClassifier<T>::type htype;
    return declareProperty (name, property, doc, htype());
  }

declareProperty() [6/6]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T > &  t )

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

                                                                 {
    typedef typename SG::HandleClassifier<T>::type htype;
    return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
  }

detStore()

const ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const

inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

{ return m_detStore; }

evtStore() [1/2]

ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

{ return m_evtStore; }

evtStore() [2/2]

const ServiceHandle<StoreGateSvc>& AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( ) const

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 90 of file AthCommonDataStore.h.

{ return m_evtStore; }

extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase &  ExtraDeps )

protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

finalize()

StatusCode TrkDriftCircleMath::MdtSegmentT0Fitter::finalize ( )

overridevirtual

Definition at line 179 of file MdtSegmentT0Fitter.cxx.

                                          {
 
    double scaleFactor = m_ntotalCalls != 0 ? 1./(double)m_ntotalCalls : 1.;
 
    ATH_MSG_INFO( "Summarizing fitter statistics " << "\n"
                  << " Total fits       " << std::setw(10) << m_ntotalCalls << "   " << scaleFactor*m_ntotalCalls << "\n"
                  << " hits > 2         " << std::setw(10) << m_npassedNHits << "   " << scaleFactor*m_npassedNHits << "\n"
                  << " hit consis.      " << std::setw(10) << m_npassedSelectionConsistency << "   " << scaleFactor*m_npassedSelectionConsistency << "\n"
                  << " 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  );
    return StatusCode::SUCCESS;
  }

fit() [1/2]

bool TrkDriftCircleMath::MdtSegmentT0Fitter::fit ( Segment &  result, const Line &  line, const DCOnTrackVec &  dcs, const HitSelection &  selection, double  t0Seed  ) const

overridevirtual

Normalize the mean positions

go to coordinates centered at the average of the hits

Reimplemented from TrkDriftCircleMath::DCSLFitter.

Definition at line 193 of file MdtSegmentT0Fitter.cxx.

                                                                                                                                               {
    ++m_ntotalCalls;
    const MdtIdHelper& id_helper{m_idHelperSvc->mdtIdHelper()};
    ATH_MSG_DEBUG("New seg: ");
    const EventContext& ctx{Gaudi::Hive::currentContext()};
    SG::ReadCondHandle<MuonCalib::MdtCalibDataContainer> calibData{m_calibDbKey, ctx};
    if (!calibData.isValid()) {
      ATH_MSG_FATAL("Failed to retrieve Mdt calibration object "<<m_calibDbKey.fullKey());
      return false;
    }
    std::array<const MuonCalib::MdtRtRelation*, 2> rtRelations{};
    {
      unsigned int nRel{0};
      for (unsigned int i = 0; i < dcs.size() ; ++i) {        
        const Identifier id = dcs[i].rot()->identify();
        const int mlIdx = id_helper.multilayer(id) -1;
        if (rtRelations[mlIdx]) continue;
        rtRelations[mlIdx] = calibData->getCalibData(id, msgStream())->rtRelation.get();
        ++nRel;
        if (nRel == 2) break;
      }
    }
    const DCOnTrackVec& dcs_keep = dcs;
 
    unsigned int N = dcs_keep.size();
 
    result.setT0Shift(-99999,-99999);
 
    if(N<2) {
      return false;
    }
    ++m_npassedNHits;
    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){
      if( selection[i] == 0 ) ++used;
    }
    if(used < 2){
      ATH_MSG_DEBUG("TOO FEW HITS SELECTED");
      return false;
    }
    ++m_npassedNSelectedHits;
    
    if(used < m_minHits) {
      ATH_MSG_DEBUG("FEWER THAN Minimum HITS N " << m_minHits << " total hits " <<N<<" used " << used);
 
      //
      //     Copy driftcircles and reset the drift radii as they might have been overwritten
      //     after a succesfull t0 fit
      //
 
      DCOnTrackVec dcs_new;
      dcs_new.reserve(dcs.size());
      
      double chi2p = 0.;    
      int n_elements = dcs.size();
      for(int i=0; i< n_elements; ++i ){
          const DriftCircle* ds  = & dcs[i];
          if(std::abs(ds->r()-ds->rot()->driftRadius())>m_dRTol) 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->driftState(), ds->id(), ds->rot(), ds->index());
          DCOnTrack dc_new(dc_keep, 0., 0.);
          
          dc_new.state(dcs[i].state());
          dcs_new.push_back( dc_new );
          if( selection[i] == 0 ){
            double t = ds->rot()->driftTime();
            const unsigned int mlIdx = id_helper.multilayer(ds->rot()->identify()) - 1;
            const MuonCalib::MdtRtRelation *rtInfo = rtRelations[mlIdx];
            
            double tUp = rtInfo->rt()->tUpper();
            double tLow = rtInfo->rt()->tLower();
            
            if(t<tLow) chi2p += sq(t-tLow)*0.1;
            else if(t>tUp) chi2p += sq(t-tUp)*0.1;
            }
      }
      
      if(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());
      int iok = 0;
      if(oldrefit) iok = 1;
      ATH_MSG_DEBUG(" chi2 total " << result.chi2() << " angle " << result.line().phi() << " y0 " << result.line().y0()  << " nhits "<< selection.size() << " refit ok " << iok);
      return oldrefit;
    
    }
    
    ATH_MSG_DEBUG("FITTING FOR T0 N "<<N<<" used " << used);
    
 
    ++m_npassedMinHits;
 
   
    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};
 
    std::vector<HitCoords> hits{};
    hits.reserve(N);
    FunctionToMinimize minFunct(used);
 
    {
      unsigned int ii{0};
      for(const DCOnTrack& keep_me : dcs_keep ){
        const Muon::MdtDriftCircleOnTrack *roto = keep_me.rot();
        const unsigned int mlIdx = id_helper.multilayer(roto->identify()) - 1;
        const MuonCalib::MdtRtRelation *rtInfo = rtRelations[mlIdx];
 
        const double newerror = m_scaleErrors ? keep_me.drPrecise() : keep_me.dr();
        const double w = newerror >0. ? sq(1./newerror) : 0.;
        hits.emplace_back(keep_me.x(), roto->driftTime(), keep_me.y(), w, std::abs(roto->driftRadius()), rtInfo->rt());
        HitCoords& coords = hits.back();
        coords.dr = keep_me.dr();
        coords.rejected = selection[ii];
        ATH_MSG_DEBUG("DC:  (" << coords.y << "," << coords.z << ") R = " << coords.r << " W " << coords.w 
                               <<" t " <<coords.t<< " id: "<<keep_me.id()<<" sel " <<coords.rejected);
        if (!coords.rejected) {          
           S += coords.w;
           Sz+= coords.z* coords.w;
           Sy+= coords.y * coords.w;
        }
        ++ii;
      }
    }
    
    const double inv_S = 1. / S;
    Zc = Sz*inv_S;
    Yc = Sy*inv_S;
 
    ATH_MSG_DEBUG("Yc " << Yc << " Zc " << Zc);
 
    for(HitCoords& coords : hits) {
      coords.y -= Yc;
      coords.z -= Zc;
    }
 
    int selcount{0};
  
    // 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
    double min_tlower{std::numeric_limits<float>::max()}, max_tupper{-std::numeric_limits<float>::max()};
 
    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
  
    for(HitCoords& coords : hits) {
      if(coords.rejected) continue;
      
      double r2tval = r2t_ext(coords.rt, coords.r) ;
      const double tl = coords.rt->tLower();
      const double th = coords.rt->tUpper();
      const double tee0 = coords.t - r2tval;
 
      min_tlower = std::min(min_tlower, coords.t - tl);
      max_tupper = std::max(max_tupper, coords.t - th);
      
      
      ATH_MSG_DEBUG(" z "<<coords.z <<" y "<<coords.y<<" r "<<coords.r
                  <<" t "<<coords.t<<" t0 "<<tee0<<" tLower "<<tl<<" tUpper "<<th);
      t0seed += tee0;
      st0 += sq(tee0);
      if(tee0 < min_t0 && std::abs(r2tval) < R2TSPURIOUS) min_t0 = tee0;
 
      minFunct.addCoords(coords);
 
      selcount++;
      
    }
    t0seed /= selcount;
    st0 = st0/selcount - t0seed*t0seed;
    st0 = st0 > 0. ? std::sqrt(st0) : 0.;
 
    ATH_MSG_DEBUG(" t0seed "<<t0seed<<" sigma "<<st0<< " min_t0 "<<min_t0);
 
    // ************************* seed the parameters
    const double theta = line.phi();
    double cosin = std::cos(theta);
    double sinus = std::sin(theta);
 
    if ( sinus < 0.0 ) {
      sinus = -sinus;
      cosin = -cosin;
    } else if ( sinus == 0.0 && cosin < 0.0 ) {
      cosin = -cosin;
    }
    
    ATH_MSG_DEBUG("before fit theta "<<theta<<" sinus "<<sinus<< " cosin "<< cosin);
 
    double d = line.y0() + 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
 
    int nml1p{0}, nml2p{0}, nml1n{0}, nml2n{0};
    int ii{-1};
    for(const DCOnTrack& keep_me : dcs_keep){
      ++ii;
      const HitCoords& coords = hits[ii];
      if(coords.rejected) continue;
      const double sdist = d*cosin + coords.z*sinus - coords.y*cosin; // same thing as |a*z - y + b|/sqrt(1. + a*a);
      nml1p+=(keep_me.id().ml()==0&&sdist > 0);
      nml1n+=(keep_me.id().ml()==0&&sdist < 0);
      nml2p+=(keep_me.id().ml()==1&&sdist > 0);
      nml2n+=(keep_me.id().ml()==1&&sdist < 0);
    }
 
// Define t0 constraint in Minuit
    int t0Error = STRONG_TOPO_T0ERROR;
    if (nml1p+nml2p < 2 || nml1n+nml2n < 2) t0Error = WEAK_TOPO_T0ERROR;
 
    minFunct.setT0Error(t0Error);
 
// Reject topologies where in one of the Multilayers no +- combination is present
    if((nml1p<1||nml1n<1)&&(nml2p<1||nml2n<1)&&m_rejectWeakTopologies) {
       ATH_MSG_DEBUG("Combination rejected for positive radii ML1 " <<  nml1p << " ML2 " <<  nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << used << " t0 Error " << t0Error);
      DCOnTrackVec::const_iterator it = dcs.begin();
      DCOnTrackVec::const_iterator it_end = dcs.end();
      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) 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->driftState(), ds->id(), ds->rot(), ds->index() );
        DCOnTrack dc_new(dc_keep, 0., 0.);
        dc_new.state(dcs[i].state());
        dcs_new.push_back( std::move(dc_new) );
        if( selection[i] == 0 ){
          double t = ds->rot()->driftTime();
          const unsigned int mlIdx = id_helper.multilayer(ds->rot()->identify()) - 1;
          const MuonCalib::MdtRtRelation *rtInfo = rtRelations[mlIdx];
          double tUp = rtInfo->rt()->tUpper();
          double tLow = rtInfo->rt()->tLower();
          if(t<tLow) chi2p += sq(t-tLow)*0.1;
          if(t>tUp) chi2p += sq(t-tUp)*0.1;
        }
      }
      if(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
      result.set( chi2p, result.ndof(),  result.dtheta(),  result.dy0() );
      return oldrefit;
    }  // end rejection of weak topologies
 
    ATH_MSG_DEBUG("positive radii ML1 " <<  nml1p << " ML2 " <<  nml2p << " negative radii ML1 " << nml1n << " ML " << nml2n << " used hits " << used << " t0 Error " << t0Error);
 
    constexpr std::array<Double_t,3> step{0.01 , 0.01 , 0.1 };
    // starting point
    std::array<Double_t,3> variable{theta,d,0};
    // if t0Seed value from outside use this
    if(t0Seed > -999.) variable[2] = t0Seed;
 
    ROOT::Minuit2::Minuit2Minimizer minimum("algoName");
    minimum.SetMaxFunctionCalls(10000);
    minimum.SetTolerance(0.001);
    minimum.SetPrintLevel(-1);
    if(msgLvl(MSG::VERBOSE)) minimum.SetPrintLevel(1);
 
    if (m_floatDir){
      minimum.SetVariable(0,"a",variable[0], step[0]);
      minimum.SetVariable(1,"b",variable[1], step[1]);
    } else {
      minimum.SetFixedVariable(0,"a", variable[0]);
      minimum.SetFixedVariable(1,"b", variable[1]);
    }
    
    minimum.SetVariable(2,"t0",variable[2], step[2]);
    
    minimum.SetFunction(minFunct);
 
    // do the minimization
    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());
 
    ++m_npassedMinuitFit;
 
    // Get the fit values
    double aret=results[0];
    double aErr=errors[0];
    double dtheta = aErr;
    double tana = std::tan(aret); // tangent of angle
    double ang = aret;  // between zero and pi
    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 = std::atan2(sinus, cosin);
    double b=results[1];
    double bErr=errors[1];
    double t0=results[2];
    double t0Err=errors[2];
    double dy0 = cosin * bErr - b * sinus * aErr;
 
    const double del_t = std::abs(hits[0].rt->radius((t0+t0Err)) - hits[0].rt->radius(t0)) ;
 
    
    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);
    
    d = b * cosin;
    if(msg().level() <=MSG::DEBUG) {
      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 << endmsg;
      }
    }
 
    result.dcs().clear();
    result.clusters().clear();
    result.emptyTubes().clear();
 
     ATH_MSG_DEBUG("after fit theta "<<ang<<" sinus "<<sinus<< " cosin "<< cosin);
 
    double chi2 = 0;
    unsigned int nhits(0);
    
    // calculate predicted hit positions from track parameters
  
    ATH_MSG_DEBUG("------NEW HITS------");
    int i{-1};
    for(const HitCoords& coords : hits){
      ++i;
      const DCOnTrack& keep_me{dcs_keep[i]};
      const double uppercut = coords.rt->tUpper();
      const double lowercut = coords.rt->tLower();
      
      double rad = coords.rt->radius(coords.t-t0);
      if(coords.t-t0<lowercut) rad = coords.rt->radius(lowercut);
      if(coords.t-t0>uppercut) rad = coords.rt->radius(uppercut);
      if (coords.w==0) {
        ATH_MSG_WARNING("coords.w==0, continuing");
        continue;
      }
      double drad = 1.0/std::sqrt(coords.w) ;
 
      double yl = (coords.y -  tana*coords.z - b);
      
      ATH_MSG_DEBUG("i "<<i<<" ");
      
 
      double dth = -(sinus*coords.y + cosin*coords.z)*dtheta;
      double residuals = std::abs(yl)/std::sqrt(1+tana*tana) - rad;
      
      ATH_MSG_DEBUG(" dth "<<dth<<" dy0 "<<dy0<<" del_t "<<del_t);
      
 
      double errorResiduals = std::hypot(dth, dy0, del_t);
 
      // derivatives of the residual 'R'
      std::array<double,3> deriv{};
      // del R/del theta
      double dd = coords.z * sinus + b *cosin - coords.y * cosin;
      deriv[0] = sign(dd) * (coords.z * cosin - b * sinus + coords.y * sinus);
      // del R / del b
      deriv[1] = sign(dd) * cosin ;
      // del R / del t0
 
      deriv[2] = -1* coords.rt->driftvelocity(coords.t-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];
        }
      }
      ATH_MSG_DEBUG(" covsquared " << covsq);
      if( covsq < 0. && msgLvl(MSG::DEBUG)){
        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]);
            }
        }
      }
      
      covsq = covsq > 0. ? std::sqrt(covsq) : 0.;
      const DriftCircle* ds  = & keep_me;
      if (m_propagateErrors) drad = coords.dr;
      
      DriftCircle dc_newrad(keep_me.position(), rad, drad, ds->driftState(), keep_me.id(), ds->rot(), keep_me.index() );
      DCOnTrack dc_new(dc_newrad, residuals, covsq);
      dc_new.state(keep_me.state());
 
      ATH_MSG_DEBUG("T0 Segment hit res "<<residuals<<" eres "<<errorResiduals<<" covsq "<<covsq<<" ri " << coords.r<<" ro "<<rad<<" drad "<<drad << " sel "<<selection[i]<< " inv error " << coords.w);
 
      if(!coords.rejected) {
        ++nhits;
        if (!m_propagateErrors) {
          chi2 += sq(residuals)*coords.w;
        } else {
          chi2 += sq(residuals)/sq(drad);
        }
        ATH_MSG_DEBUG("T0 Segment hit res "<<residuals<<" eres "<<errorResiduals<<" covsq "<<covsq<<" ri " << coords.r<<" radius after t0 "<<rad<<" radius error "<< drad <<  " original error " << coords.dr);
// Put chi2 penalty for drift times outside window
        if (coords.t-t0> uppercut ) { // too large
              chi2  += sq(coords.t-t0-uppercut)*0.1;
        }else if (coords.t-t0 < lowercut ) {// too small
              chi2 += sq(coords.t-t0-lowercut)*0.1;
        }
      }
      result.dcs().push_back( dc_new );
    }
 
    double oldshift = result.t0Shift();
    ATH_MSG_DEBUG("end fit old "<<oldshift<< " new " <<t0);
    // Final Minuit Fit result
    if(nhits==NUMPAR) {
      nhits++;
      chi2 += 1.;
    }
    result.set( chi2, nhits-NUMPAR, dtheta, -1.*dy0 );
    result.line().set( LocVec2D( Zc - sinus*d, Yc + cosin*d ), ang );
    if(t0==0.) t0=0.00001;
    result.setT0Shift(t0,t0Err);
   
    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)<<(chi2/(nhits-NUMPAR) > 5 ? " NOT ":" ")<< "GOOD");
    ATH_MSG_DEBUG("chi2 "<<chi2<<" per dof "<<chi2/(nhits-NUMPAR));
    
    return true;
  }

fit() [2/2]

bool TrkDriftCircleMath::MdtSegmentT0Fitter::fit ( Segment &  result, const Line &  line, const DCOnTrackVec &  dcs, double  t0Seed  ) const

inlineoverridevirtual

Reimplemented from TrkDriftCircleMath::DCSLFitter.

Definition at line 63 of file MdtSegmentT0Fitter.h.

                                                                                                                       { 
    HitSelection selection(dcs.size(),0);
    return fit( result, line, dcs, selection, t0Seed ); 
  }

getFitter()

virtual const DCSLFitter* TrkDriftCircleMath::MdtSegmentT0Fitter::getFitter ( ) const

inlineoverridevirtual

fitter factory

Returns

provides pointer to fitter, ownsership not passed to client

Implements Muon::IDCSLFitProvider.

Definition at line 36 of file MdtSegmentT0Fitter.h.

{ return this; }

initialize()

StatusCode TrkDriftCircleMath::MdtSegmentT0Fitter::initialize ( )

overridevirtual

Definition at line 174 of file MdtSegmentT0Fitter.cxx.

                                            {
    ATH_CHECK(m_calibDbKey.initialize());
    return StatusCode::SUCCESS;
  }

inputHandles()

virtual std::vector<Gaudi::DataHandle*> AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const

overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

interfaceID()

static const InterfaceID& Muon::IDCSLFitProvider::interfaceID ( )

inlinestaticinherited

Definition at line 21 of file IDCSLFitProvider.h.

                                             {
          static const InterfaceID IID_IDCSLFitProvider("Muon::IDCSLFitProvider",1,0);
         return IID_IDCSLFitProvider;
      }

msg() [1/2]

MsgStream& AthCommonMsg< AlgTool >::msg ( ) const

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

msg() [2/2]

MsgStream& AthCommonMsg< AlgTool >::msg ( const MSG::Level  lvl ) const

inlineinherited

Definition at line 27 of file AthCommonMsg.h.

                                                  {
    return this->msgStream(lvl);
  }

msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level  lvl ) const

inlineinherited

Definition at line 30 of file AthCommonMsg.h.

                                               {
    return this->msgLevel(lvl);
  }

outputHandles()

virtual std::vector<Gaudi::DataHandle*> AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const

overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

renounce()

std::enable_if_t<std::is_void_v<std::result_of_t<decltype(&T::renounce)(T)> > && !std::is_base_of_v<SG::VarHandleKeyArray, T> && std::is_base_of_v<Gaudi::DataHandle, T>, void> AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T &  h )

inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

  {
    h.renounce();
    PBASE::renounce (h);
  }

renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray &  handlesArray )

inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

                                                          {
    handlesArray.renounce();
  }

sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )

overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in DerivationFramework::CfAthAlgTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and asg::AsgMetadataTool.

sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )

overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase &  )

inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

                                               {
    // debug() << "updateVHKA for property " << p.name() << " " << p.toString() 
    //         << "  size: " << m_vhka.size() << endmsg;
    for (auto &a : m_vhka) {
      std::vector<SG::VarHandleKey*> keys = a->keys();
      for (auto k : keys) {
        k->setOwner(this);
      }
    }
  }

Member Data Documentation

m_calibDbKey

SG::ReadCondHandleKey<MuonCalib::MdtCalibDataContainer> TrkDriftCircleMath::MdtSegmentT0Fitter::m_calibDbKey

private

Initial value:

{this, "CalibDataKey", "MdtCalibConstants",
                                                                       "Conditions object containing the calibrations"}

Definition at line 40 of file MdtSegmentT0Fitter.h.

m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore

privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

m_dRTol

Gaudi::Property<float> TrkDriftCircleMath::MdtSegmentT0Fitter::m_dRTol {this,"dRTolerance",0.1}

private

Definition at line 49 of file MdtSegmentT0Fitter.h.

m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore

privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

m_floatDir

Gaudi::Property<bool> TrkDriftCircleMath::MdtSegmentT0Fitter::m_floatDir

private

Initial value:

{this, "FloatSegDirection", false, 
                                      "If set to true, the line of the segment is simultaenously fitted with t0"}

Definition at line 51 of file MdtSegmentT0Fitter.h.

m_idHelperSvc

ServiceHandle<Muon::IMuonIdHelperSvc> TrkDriftCircleMath::MdtSegmentT0Fitter::m_idHelperSvc {this, "MuonIdHelperSvc", "Muon::MuonIdHelperSvc/MuonIdHelperSvc"}

private

Definition at line 43 of file MdtSegmentT0Fitter.h.

m_minHits

Gaudi::Property<int> TrkDriftCircleMath::MdtSegmentT0Fitter::m_minHits {this,"MinimumHits",4,"minimum number of selected hits for t0 fit. Otherwise use default"}

private

Definition at line 48 of file MdtSegmentT0Fitter.h.

m_npassedMinHits

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_npassedMinHits {0}

mutableprivate

Definition at line 58 of file MdtSegmentT0Fitter.h.

m_npassedMinuitFit

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_npassedMinuitFit {0}

mutableprivate

Definition at line 59 of file MdtSegmentT0Fitter.h.

m_npassedNHits

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_npassedNHits {0}

mutableprivate

Definition at line 55 of file MdtSegmentT0Fitter.h.

m_npassedNSelectedHits

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_npassedNSelectedHits {0}

mutableprivate

Definition at line 57 of file MdtSegmentT0Fitter.h.

m_npassedSelectionConsistency

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_npassedSelectionConsistency {0}

mutableprivate

Definition at line 56 of file MdtSegmentT0Fitter.h.

m_ntotalCalls

std::atomic_uint TrkDriftCircleMath::MdtSegmentT0Fitter::m_ntotalCalls {0}

mutableprivate

Definition at line 54 of file MdtSegmentT0Fitter.h.

m_propagateErrors

Gaudi::Property<bool> TrkDriftCircleMath::MdtSegmentT0Fitter::m_propagateErrors {this,"PropagateErrors",true,"propagate errors"}

private

Definition at line 47 of file MdtSegmentT0Fitter.h.

m_rejectWeakTopologies

Gaudi::Property<bool> TrkDriftCircleMath::MdtSegmentT0Fitter::m_rejectWeakTopologies {this,"RejectWeakTopologies",true,"reject topolgies that do not have at least one +- combination in one multilayer"}

private

Definition at line 45 of file MdtSegmentT0Fitter.h.

m_scaleErrors

Gaudi::Property<bool> TrkDriftCircleMath::MdtSegmentT0Fitter::m_scaleErrors {this,"RescaleErrors",true,"rescale errors in fit"}

private

Definition at line 46 of file MdtSegmentT0Fitter.h.

m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared

privateinherited

Definition at line 399 of file AthCommonDataStore.h.

m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka

privateinherited

Definition at line 398 of file AthCommonDataStore.h.

---

The documentation for this class was generated from the following files:

• MdtSegmentT0Fitter.h
• MdtSegmentT0Fitter.cxx