#include <CscBipolarStripFitter.h>

Inheritance diagram for CscBipolarStripFitter:

Collaboration diagram for CscBipolarStripFitter:

Public Types
typedef std::vector< float >	ChargeList

Public Member Functions
	CscBipolarStripFitter (const std::string &, const std::string &, const IInterface *)
	~CscBipolarStripFitter ()=default
StatusCode	initialize () override
Result	fit (const ChargeList &charges, double samplingTime, Identifier &stripId) const
virtual Result	fit (const ChargeList &ChargeList, double samplingTime, bool samplingPhase, Identifier &sid) const
virtual Result	fit (const Muon::CscStripPrepData &strip) const
ServiceHandle< StoreGateSvc > &	evtStore ()
	The standard `StoreGateSvc` (event store) Returns (kind of) a pointer to the `StoreGateSvc`.
const ServiceHandle< StoreGateSvc > &	detStore () const
	The standard `StoreGateSvc/DetectorStore` Returns (kind of) a pointer to the `StoreGateSvc`.
virtual StatusCode	sysInitialize () override
	Perform system initialization for an algorithm.
virtual StatusCode	sysStart () override
	Handle START transition.
virtual std::vector< Gaudi::DataHandle * >	inputHandles () const override
	Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * >	outputHandles () const override
	Return this algorithm's output handles.
Gaudi::Details::PropertyBase &	declareProperty (Gaudi::Property< T, V, H > &t)
void	updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream &	msg () 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
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.

Private Types
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
double	FindPow (double z) const
void	Derivative (double A[][3], double fp[][1], double p0[][1], int imeas, const int *meas) const
int	TheFitter (double x, const double ex, const double initValues, int imeas, int meas, int ipar, int par, double chi2, double result) const
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey>

Static Private Member Functions
static double	FindInitValues (double x, double initValues, int *maxsample)
static void	InvertMatrix (double matrix[][3], const int dim, const int *)
static void	InvertSymmetric4x4 (double W[][4])

Private Attributes
const CscIdHelper *	m_phelper
ToolHandle< ICscCalibTool >	m_cscCalibTool
double	m_qerr
double	m_terr
double	m_qerr_fail
double	m_terr_fail
double	m_qerrprop
double	m_n
double	m_n2
double	m_zmax
double	m_bipolarNormalization
double	m_tsampling
StoreGateSvc_t	m_evtStore
	Pointer to StoreGate (event store by default)
StoreGateSvc_t	m_detStore
	Pointer to StoreGate (detector store by default)
std::vector< SG::VarHandleKeyArray * >	m_vhka
bool	m_varHandleArraysDeclared

Detailed Description

Definition at line 27 of file CscBipolarStripFitter.h.

Member Typedef Documentation

◆ ChargeList

typedef std::vector<float> ICscStripFitter::ChargeList

inherited

Definition at line 55 of file ICscStripFitter.h.

◆ StoreGateSvc_t

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

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ CscBipolarStripFitter()

CscBipolarStripFitter::CscBipolarStripFitter	(	const std::string &	type,
		const std::string &	aname,
		const IInterface *	parent )

Definition at line 29 of file CscBipolarStripFitter.cxx.

                                                                                                                    :
    AthAlgTool(type, aname, parent), m_phelper(nullptr), m_n(0), m_n2(0), m_zmax(0), m_bipolarNormalization(0), m_tsampling(0) {
    declareInterface<ICscStripFitter>(this);
    declareProperty("chargeError", m_qerr = 5500.0);
    declareProperty("timeError", m_terr = 5.0);
    declareProperty("failChargeError", m_qerr_fail = 5000.0);
    declareProperty("failTimeError", m_terr_fail = 50.0);
    declareProperty("chargeCalibrationError", m_qerrprop = 0.002);
}

◆ ~CscBipolarStripFitter()

CscBipolarStripFitter::~CscBipolarStripFitter ( )

default

Member Function Documentation

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty	(	Gaudi::Property< T, V, H > &	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());
 
  }

◆ declareProperty()

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

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

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

◆ Derivative()

void CscBipolarStripFitter::Derivative	(	double	A[][3],
		double	fp[][1],
		double	p0[][1],
		int	imeas,
		const int *	meas ) const

private

Definition at line 267 of file CscBipolarStripFitter.cxx.

                                                                                                                      {
    // calculate the derivatives and the 0th order approximation
    // around the ADC samplings
    double norm = p0[0][0];
    // double parm[3]={1.,p0[1][0],p0[2][0]};
    for (int i = 0; i < imeas; i++) {
        int ii = meas[i];
        double z = (ii - p0[1][0]) * m_tsampling / p0[2][0];
        double repquant = 0.;
        double dFdzNormalized = 0.;
        if (z > 0.) {
            repquant = FindPow(z) * std::exp(-z) / m_bipolarNormalization;
            dFdzNormalized = repquant * (m_n / z + z / (m_n2 + 1) - (m_n + 1.) / (m_n2 + 1.) - 1.);  // repquant*(m_n/z+z/(m_n+1)-2.);
        }
 
        A[ii][0] = repquant * (1. - z / (m_n2 + 1.));  // repquant*(1.-z/(m_n+1.));
        // A[ii][0] = bipolar(&z,parm);
        fp[ii][0] = norm * A[ii][0];
 
        // double normOverZmax = norm/m_bipolarNormalization;
        double commonpart = norm * dFdzNormalized;  //(z,parm);
        A[ii][1] = commonpart * (-m_tsampling / p0[2][0]);
        A[ii][2] = commonpart * (-z / p0[2][0]);
    }
    // end of derivative/zeroth order calculations
}

◆ 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.

95{ return m_detStore; }

AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore

StoreGateSvc_t m_detStore

Definition AthCommonDataStore.h:393

◆ evtStore()

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.

85{ return m_evtStore; }

AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore

StoreGateSvc_t m_evtStore

Definition AthCommonDataStore.h:390

◆ 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

◆ FindInitValues()

double CscBipolarStripFitter::FindInitValues	(	double *	x,
		double *	initValues,
		int *	maxsample )

staticprivate

Definition at line 129 of file CscBipolarStripFitter.cxx.

                                                                                          {
    // find maximum sample imax:
    double peakingTime = -99.;  // interpolated peaking time in samples
    double amplitude = -99.;    // interpolated amplitude
    double amin, amax;
    int imax = 0;
    const int numSample = 4;
    amax = x[0];
    amin = x[0];
    for (int j = 1; j < numSample; j++) {
        if (amax < x[j]) {
            amax = x[j];
            imax = j;
        }
        if (amin > x[j]) amin = x[j];
    }
 
    // calculate peak and amplitude:
    (*maxsample) = imax;
    if ((imax == 0) || (imax == numSample - 1))  // no interpolation possible
    {
        peakingTime = imax;
        amplitude = amax;
    } else  // do parabolic interpolation
    {
        double a, b, c;  // coeffients of parabola y=a*x*x+b*x+c
        a = 0.5 * (x[imax + 1] + x[imax - 1] - 2.0 * x[imax]);
        b = 0.5 * (x[imax + 1] - x[imax - 1]);
        c = x[imax];
 
        peakingTime = -0.5 * b / a;
        amplitude = a * peakingTime * peakingTime + b * peakingTime + c;
        peakingTime += imax;  // it was relative to imax
    }
 
    initValues[0] = amplitude;
    initValues[1] = peakingTime;
    // - m_zmax*initValues[2]/m_tsampling;
    return x[imax];
}

◆ FindPow()

double CscBipolarStripFitter::FindPow ( double z ) const

private

Definition at line 170 of file CscBipolarStripFitter.cxx.

                                                    {
    double zpower = std::exp(m_n * std::log(z));
    return zpower;
}

◆ fit() [1/3]

Result ICscStripFitter::fit	(	const ChargeList &	ChargeList,
		double	samplingTime,
		bool	samplingPhase,
		Identifier &	sid ) const

virtual

Reimplemented from ICscStripFitter.

Definition at line 72 of file ICscStripFitter.cxx.

                                                                                      {
    return {};
}

◆ fit() [2/3]

Result CscBipolarStripFitter::fit	(	const ChargeList &	charges,
		double	samplingTime,
		Identifier &	stripId ) const

Definition at line 74 of file CscBipolarStripFitter.cxx.

                                                                                                  {
    ATH_MSG_DEBUG("Fitting with sample time " << period);
    ATH_MSG_DEBUG(" Number of charges =" << chgs.size() << "  Charges: " << chgs);
    Result res;
 
    IdentifierHash stripHash;
    if (m_phelper->get_channel_hash(stripId, stripHash)) {
        ATH_MSG_WARNING("Unable to get CSC striphash id "
                        << " the identifier is ");
        stripId.show();
    }
    ATH_MSG_DEBUG("CalibCscStripFitter:: " << stripId << "  " << (unsigned int)stripHash);
    double initValues[3];
    for (int i = 0; i < 3; ++i) { initValues[i] = 0.; }
    //  int imax = -1;
    initValues[2] = 8.;  // m_width
 
    double noise = m_cscCalibTool->stripNoise(stripHash);
    res.status = 0;
    if (m_qerr > 0)
        res.dcharge = m_qerr;
    else
        res.dcharge = noise;
 
    res.dtime = m_terr;
    res.charge = initValues[0];
    res.time = initValues[1];
 
    initValues[1] = res.time - 2.5;
    //  int imeas =4;
    //  int meas[4] = {true,true,true,true};
    //  int ipar = 3;
    //  int par[3] = {true,true,false};
    //  double m_chi2;
    //double result[3];
 
    // Add an error proportional to the charge.
    double dqprop = m_qerrprop * res.charge;
    res.dcharge = sqrt(res.dcharge * res.dcharge + dqprop * dqprop);
    // Display result.
    ATH_MSG_DEBUG("  Status: " << res.status);
    ATH_MSG_DEBUG("  Charge: " << res.charge);
    ATH_MSG_DEBUG("  Charge err: " << res.dcharge);
    ATH_MSG_DEBUG("    Time: " << res.time);
    //ATH_MSG_DEBUG("   fit  : " << result[0] << " " << result[1] << " " << result[2]);
    return res;
}

◆ fit() [3/3]

Result ICscStripFitter::fit ( const Muon::CscStripPrepData & strip ) const

virtual

Reimplemented from ICscStripFitter.

Definition at line 79 of file ICscStripFitter.cxx.

                                                               {
    Identifier sid = strip.identify();
    //  IdentifierHash coll_hash = strip.collectionHash();
 
    Result res = fit(strip.sampleCharges(), strip.samplingTime(), strip.samplingPhase(), sid);
    res.strip = &strip;
    if (res.status) return res;
    res.time += strip.timeOfFirstSample();
    // Do we also need a phase correction here?
    return res;
}

◆ initialize()

StatusCode CscBipolarStripFitter::initialize ( )

override

CSC calibratin tool for the Condtiions Data base access

Definition at line 41 of file CscBipolarStripFitter.cxx.

                                             {
    ATH_MSG_DEBUG("Initializing " << name());
 
    ATH_MSG_DEBUG("Properties for " << name() << ":");
    ATH_MSG_DEBUG("              Charge error: " << m_qerr);
    ATH_MSG_DEBUG("                Time error: " << m_terr);
    ATH_MSG_DEBUG("         Fail charge error: " << m_qerr_fail);
    ATH_MSG_DEBUG("           Fail time error: " << m_terr_fail);
    ATH_MSG_DEBUG("  Charge calibration error: " << m_qerrprop);
    ATH_CHECK_RECOVERABLE(m_cscCalibTool.retrieve());
 
    const MuonGM::MuonDetectorManager *muDetMgr = nullptr;
    ATH_CHECK_RECOVERABLE(detStore()->retrieve(muDetMgr));
    ATH_MSG_DEBUG("Retrieved geometry.");
 
    m_phelper = muDetMgr->cscIdHelper();
 
    m_n = m_cscCalibTool->getNumberOfIntegration();    // 12.;
    m_n2 = m_cscCalibTool->getNumberOfIntegration2();  // 11.66;
    m_zmax = m_n + m_n2 + 2. - sqrt((m_n + 2. + m_n2) * (m_n + 2. + m_n2) - 4. * m_n * (m_n2 + 1));
    m_zmax *= 0.5;
    // m_n+1 - sqrt(m_n+1.0);
    m_bipolarNormalization = FindPow(m_zmax) * (1 - m_zmax / (m_n2 + 1)) * std::exp(-m_zmax);
    // FindPow(m_zmax)*(1-m_zmax/(m_n+1))*std::exp(-m_zmax);
    m_tsampling = m_cscCalibTool->getSamplingTime() / 2.;  //   50/2 =   25.;
 
    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()

const InterfaceID & ICscStripFitter::interfaceID ( )

inlinestaticinherited

Must declare this, with name of interface

Definition at line 59 of file ICscStripFitter.h.

                                            { 
        static const InterfaceID IID_ICscStripFitter("ICscStripFitter", 1, 0);
        return IID_ICscStripFitter;
    }

◆ InvertMatrix()

void CscBipolarStripFitter::InvertMatrix	(	double	matrix[][3],
		const int	dim,
		const int *	correspdim )

staticprivate

Definition at line 175 of file CscBipolarStripFitter.cxx.

                                                                                                 {
    // invert 2x2 or 3x3 symmetric matrix
    if (dim == 1) {
        int ii = correspdim[0];
        matrix[ii][ii] = 1. / matrix[ii][ii];
 
    } else if (dim == 2) {
        int ii = correspdim[0];
        int jj = correspdim[1];
        double determinant = -matrix[jj][ii] * matrix[ii][jj] + matrix[ii][ii] * matrix[jj][jj];
        // if(std::abs(determinant)<1.e-13)
        // std::cout<<" zero determinant "<<std::endl;
        double i00 = matrix[ii][ii];
        matrix[ii][ii] = matrix[jj][jj] / determinant;
        matrix[jj][jj] = i00 / determinant;
        matrix[ii][jj] = -matrix[ii][jj] / determinant;
        matrix[jj][ii] = matrix[ii][jj];
 
    } else if (dim == 3) {
        double sm12 = matrix[0][1] * matrix[0][1];
        double sm13 = matrix[0][2] * matrix[0][2];
        double sm23 = matrix[1][2] * matrix[1][2];
        double determinant = matrix[0][0] * matrix[1][1] * matrix[2][2] - matrix[0][0] * sm23 - sm12 * matrix[2][2] +
                             2. * matrix[0][1] * matrix[0][2] * matrix[1][2] - matrix[1][1] * sm13;
 
        // if(std::abs(determinant)<1.e-13)
        // std::cout << "zero determinant"<<std::endl;
        double i00 = matrix[1][1] * matrix[2][2] - sm23;
        double i11 = matrix[0][0] * matrix[2][2] - sm13;
        double i22 = matrix[0][0] * matrix[1][1] - sm12;
 
        matrix[1][0] = (matrix[0][2] * matrix[1][2] - matrix[2][2] * matrix[0][1]) / determinant;
        matrix[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1]) / determinant;
        matrix[2][1] = (matrix[0][1] * matrix[0][2] - matrix[0][0] * matrix[1][2]) / determinant;
        matrix[0][1] = matrix[1][0];
        matrix[0][2] = matrix[2][0];
        matrix[1][2] = matrix[2][1];
        matrix[0][0] = i00 / determinant;
        matrix[1][1] = i11 / determinant;
        matrix[2][2] = i22 / determinant;
 
    } else {
        // int ierr;
        // matrix.invert(ierr);
        printf("this is not a 1x1 or 2x2 or 3x3 matrix\n");
    }
}

◆ InvertSymmetric4x4()

void CscBipolarStripFitter::InvertSymmetric4x4 ( double W[][4] )

staticprivate

Definition at line 223 of file CscBipolarStripFitter.cxx.

                                                            {
    double Determinant =
        W[0][3] * W[0][3] * W[1][2] * W[1][2] - 2. * W[0][2] * W[0][3] * W[1][2] * W[1][3] + W[0][2] * W[0][2] * W[1][3] * W[1][3] -
        W[0][3] * W[0][3] * W[1][1] * W[2][2] + 2. * W[0][1] * W[0][3] * W[1][3] * W[2][2] - W[0][0] * W[1][3] * W[1][3] * W[2][2] +
        2. * W[0][2] * W[0][3] * W[1][1] * W[2][3] - 2. * W[0][1] * W[0][3] * W[1][2] * W[2][3] -
        2. * W[0][1] * W[0][2] * W[1][3] * W[2][3] + 2. * W[0][0] * W[1][2] * W[1][3] * W[2][3] + W[0][1] * W[0][1] * W[2][3] * W[2][3] -
        W[0][0] * W[1][1] * W[2][3] * W[2][3] - W[0][2] * W[0][2] * W[1][1] * W[3][3] + 2. * W[0][1] * W[0][2] * W[1][2] * W[3][3] -
        W[0][0] * W[1][2] * W[1][2] * W[3][3] - W[0][1] * W[0][1] * W[2][2] * W[3][3] + W[0][0] * W[1][1] * W[2][2] * W[3][3];
 
    W[0][1] = W[3][0] * W[3][1] * W[2][2] - W[3][0] * W[2][1] * W[3][2] - W[2][0] * W[3][1] * W[3][2] + W[1][0] * W[3][2] * W[3][2] +
              W[2][0] * W[2][1] * W[3][3] - W[1][0] * W[2][2] * W[3][3];
    W[0][2] = -W[3][0] * W[2][1] * W[3][1] + W[2][0] * W[3][1] * W[3][1] + W[3][0] * W[1][1] * W[3][2] - W[1][0] * W[3][1] * W[3][2] -
              W[2][0] * W[1][1] * W[3][3] + W[1][0] * W[2][1] * W[3][3];
    W[0][3] = W[3][0] * W[2][1] * W[2][1] - W[2][0] * W[2][1] * W[3][1] - W[3][0] * W[1][1] * W[2][2] + W[1][0] * W[3][1] * W[2][2] +
              W[2][0] * W[1][1] * W[3][2] - W[1][0] * W[2][1] * W[3][2];
    W[1][2] = W[3][0] * W[3][0] * W[2][1] - W[2][0] * W[3][0] * W[3][1] - W[1][0] * W[3][0] * W[3][2] + W[0][0] * W[3][1] * W[3][2] +
              W[1][0] * W[2][0] * W[3][3] - W[0][0] * W[2][1] * W[3][3];
 
    W[1][3] = -W[2][0] * W[3][0] * W[2][1] + W[2][0] * W[2][0] * W[3][1] + W[1][0] * W[3][0] * W[2][2] - W[0][0] * W[3][1] * W[2][2] -
              W[1][0] * W[2][0] * W[3][2] + W[0][0] * W[2][1] * W[3][2];
    W[2][3] = W[2][0] * W[3][0] * W[1][1] - W[1][0] * W[3][0] * W[2][1] - W[1][0] * W[2][0] * W[3][1] + W[0][0] * W[2][1] * W[3][1] +
              W[1][0] * W[1][0] * W[3][2] - W[0][0] * W[1][1] * W[3][2];
 
    double W00 = -W[3][1] * W[3][1] * W[2][2] + 2. * W[2][1] * W[3][1] * W[3][2] - W[1][1] * W[3][2] * W[3][2] -
                 W[2][1] * W[2][1] * W[3][3] + W[1][1] * W[2][2] * W[3][3];
    double W11 = -W[3][0] * W[3][0] * W[2][2] + 2. * W[2][0] * W[3][0] * W[3][2] - W[0][0] * W[3][2] * W[3][2] -
                 W[2][0] * W[2][0] * W[3][3] + W[0][0] * W[2][2] * W[3][3];
    double W22 = -W[3][0] * W[3][0] * W[1][1] + 2. * W[1][0] * W[3][0] * W[3][1] - W[0][0] * W[3][1] * W[3][1] -
                 W[1][0] * W[1][0] * W[3][3] + W[0][0] * W[1][1] * W[3][3];
    double W33 = -W[2][0] * W[2][0] * W[1][1] + 2. * W[1][0] * W[2][0] * W[2][1] - W[0][0] * W[2][1] * W[2][1] -
                 W[1][0] * W[1][0] * W[2][2] + W[0][0] * W[1][1] * W[2][2];
 
    for (int i = 0; i < 3; i++) {
        for (int j = 1; j < 4; j++) {
            if (i >= j) continue;
            W[i][j] = W[i][j] / Determinant;
            W[j][i] = W[i][j];
        }
    }
    W[0][0] = W00 / Determinant;
    W[1][1] = W11 / Determinant;
    W[2][2] = W22 / Determinant;
    W[3][3] = W33 / Determinant;
}

◆ msg()

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

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

◆ 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 asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and DerivationFramework::CfAthAlgTool.

◆ 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.

◆ TheFitter()

int CscBipolarStripFitter::TheFitter	(	double *	x,
		const double	ex,
		const double *	initValues,
		int	imeas,
		int *	meas,
		int	ipar,
		int *	par,
		double *	chi2,
		double *	result ) const

private

Definition at line 294 of file CscBipolarStripFitter.cxx.

                                                           {
    // maximum iterations
    const int maxIter = 7;
    // tolerances
    double fitTolerance0 = 0.1;
    double fitTolerance1 = 0.01;
    // CLHEP::HepMatrix p0(3,1,0); // the matrix of the initial fit parameters
    double p0[3][1];
    for (int j = 0; j < 3; j++) p0[j][0] = initValues[j];
 
    // CLHEP::HepMatrix m(4,1,0); // the matrix of ADC measurements (samples: 0,1,2,3)
    double m[4][1];
    // CLHEP::HepMatrix W(4,4,0); // the error matrix of the ADC measurements
    double W[4][4];
    for (int i = 0; i < 4; i++) {
        m[i][0] = x[i];
        W[i][i] = ex * ex;
    }
    // covariances
    W[0][1] = 0.03 * ex * ex;
    W[0][2] = -0.411 * ex * ex;
    W[0][3] = -0.188 * ex * ex;
    W[1][2] = 0.0275 * ex * ex;
    W[1][3] = -0.4303 * ex * ex;
    W[2][3] = 0. * ex * ex;
    W[1][0] = W[0][1];
    W[2][0] = W[0][2];
    W[3][0] = W[0][3];
    W[2][1] = W[1][2];
    W[3][1] = W[1][3];
    W[3][2] = W[2][3];
 
    // WW.invert(ierr);
    InvertSymmetric4x4(W);
 
    // Taylor std::expansion of the bipolar pulse model around the
    // samplings : F(x) = F(p0) + A *(p-p0) + higher.order
    // CLHEP::HepMatrix fp(4,1,0); // the matrix of 0th order approximation
    double fp[4][1];
    // CLHEP::HepMatrix A(4,3,0);  // the matrix of derivatives
    double A[4][3];
    for (int i = 0; i < 4; ++i) {
        fp[i][0] = 0.;
        for (int j = 0; j < 3; ++j) { A[i][j] = 0.; }
    }
    // remarks :
    // if the pulse peaks in the last sampling fit with a constant shaping time
    // if the pulse peaks in the first sampling fit without using the last sampling
    // (too large contribution to the chi2
    int counter = 0;
    bool converged = false;
    double amplitudeChangeOld = 0.;
    bool diverganceCandidate = false;
    // CLHEP::HepMatrix weight(3,3,1); // weight matrix allocated once
    // the non-fitted parts are taken care appropriately
    // at least if the fitting parameters or measurements
    // don't change during the fitting procedure
    double weight[3][3];
 
    // CLHEP::HepMatrix residFactor(3,4,0); // residFactor allocated once
    double residFactor[3][4];
 
    for (int i = 0; i < 3; ++i) {
        for (int j = 0; j < 4; ++j) {
            if (j < 3) { weight[i][j] = 0.; }
            residFactor[i][j] = 0.;
        }
    }
    weight[0][0] = 1.;
    weight[1][1] = 1.;
    weight[2][2] = 1.;
 
    while (!converged && counter < maxIter)  // fit loop
    {
        Derivative(A, fp, p0, imeas, meas);  // calculate the matrix of derivatives and 0th order approximation
        // matrix multiplication
        // the weight matrix is symmetric
        // weight= A.T()*W*A;//.assign(A.T()*W*A);
 
        double helpmatrix[4][3];
        for (int i = 0; i < 4; ++i) {
            for (int j = 0; j < 3; ++j) { helpmatrix[i][j] = 0.; }
        }
 
        for (int i = 0; i < imeas; i++) {
            int ii = meas[i];
            for (int j = 0; j < ipar; j++) {
                int jj = par[j];
                for (int k = 0; k < imeas; k++) {
                    int kk = meas[k];
                    helpmatrix[ii][jj] += W[ii][kk] * A[kk][jj];
                }
            }
        }
        for (int i = 0; i < ipar; i++) {
            int ii = par[i];
            for (int j = i; j < ipar; j++) {
                int jj = par[j];
                weight[ii][jj] = 0.;
                for (int k = 0; k < imeas; k++) {
                    int kk = meas[k];
                    weight[ii][jj] += A[kk][ii] * helpmatrix[kk][jj];  // A[kk][ii]*A[kk][jj];
                }
                // weight[ii][jj]*=W[0][0];
                weight[jj][ii] = weight[ii][jj];
            }
        }
        // weight.invert(ierr); // inversion of weight matrix
        // hand-made inversion of 2x2 or 3x3 symmetric matrix
        InvertMatrix(weight, ipar, par);
 
        // calculate W*(A.T()*W)
        // residFactor = weight*(A.T()*W);
        double helpmatrix2[3][4];
        for (int i = 0; i < 3; ++i) {
            for (int j = 0; j < 4; ++j) { helpmatrix2[i][j] = 0.; }
        }
 
        for (int i = 0; i < ipar; i++) {
            int ii = par[i];
            for (int j = 0; j < imeas; j++) {
                int jj = meas[j];
                for (int k = 0; k < imeas; k++) {
                    int kk = meas[k];
                    helpmatrix2[ii][jj] += A[kk][ii] * W[kk][jj];
                }
            }
        }
 
        for (int i = 0; i < ipar; i++) {
            int ii = par[i];
            for (int j = 0; j < imeas; j++) {
                int jj = meas[j];
                residFactor[ii][jj] = 0.;
                for (int k = 0; k < ipar; k++) {
                    int kk = par[k];
                    residFactor[ii][jj] += weight[ii][kk] * helpmatrix2[kk][jj];
                }
                // residFactor[ii][jj]*=W[0][0];
            }
        }
 
        double paramDiff[3];
        for (int i = 0; i < 3; ++i) { paramDiff[i] = 0.; }
 
        for (int i = 0; i < ipar; i++) {
            int ii = par[i];
            // estimation of new parameters
            // paramDiff[i][0] += (weight*(A.T()*W)*(m-fp))[i][0];
            for (int j = 0; j < imeas; j++) {
                int jj = meas[j];
                paramDiff[ii] += residFactor[ii][jj] * (m[jj][0] - fp[jj][0]);
            }
            p0[ii][0] += paramDiff[ii];
        }
        std::cout << "##### " << p0[0][0] << " " << p0[1][0] << " " << p0[2][0] << std::endl;
        // if the parameters are not physical, keep them sensible
        // if peaking time less than -0.5
        double peakingTime = p0[1][0] + m_zmax * p0[2][0] / m_tsampling;
        if (peakingTime < -0.5 || peakingTime > 3.) p0[1][0] = initValues[1];
 
        if (p0[0][0] < 0.) p0[0][0] = initValues[0];
        double amplitudeChangeNew = std::abs(paramDiff[0]);
        if (std::abs(paramDiff[0]) < fitTolerance0 && std::abs(paramDiff[1]) < fitTolerance1) {
            converged = true;
            // calculate chi2
            // (m-fp).T()*W*(m-fp)
            double residual[4];
            for (int i = 0; i < 4; ++i) { residual[i] = 0.; }
 
            for (int i = 0; i < imeas; i++) {
                int ii = meas[i];
                residual[i] = m[ii][0] - fp[ii][0];
            }
 
            double tmpChi2 = 0.;
            double helpmatrixchi2[4][1];
            for (int i = 0; i < 4; ++i) { helpmatrixchi2[i][0] = 0.; }
            for (int i = 0; i < imeas; i++) {
                int ii = meas[i];
                for (int k = 0; k < imeas; k++) {
                    int kk = meas[k];
                    helpmatrixchi2[ii][0] += W[ii][kk] * residual[kk];
                }
            }
            for (int k = 0; k < imeas; k++) {
                int kk = meas[k];
                tmpChi2 += residual[kk] * helpmatrixchi2[kk][0];
                // std::cout << residual[kk] << " ";//*W[kk][kk]*residual[kk]<< " ";
            }
            //    std::cout<<std::endl;
            (*chi2) = tmpChi2;
        } else if (counter > 4 && (amplitudeChangeNew > 4. * amplitudeChangeOld)) {
            if (diverganceCandidate) {
                // diverging fit
                // return parabola interpolation
                printf("%3.2f %3.2f %3.2f %3.2f\n ", x[0], x[1], x[2], x[3]);
                printf("Diverging fit\n");
                return 4;
            } else
                diverganceCandidate = true;
        }
        if (p0[0][0] < 0.) {
            // negative amplitude
            // fit diverged
            // return parabola
            return 4;
        }
        // if after a couple of iterations the amplitude is low
        // reduce the tolerances (or the maximum iterations)
        // low amplitude pulses tend to oscillate and exhaust all iterations
        if (p0[0][0] < 20.) {
            fitTolerance0 = 0.1;
            fitTolerance1 = 0.05;
        }
        amplitudeChangeOld = amplitudeChangeNew;
        counter++;
    }
    /*
    std::cout << " Error matrix "<<std::endl;
    for(int i=0;i<3;i++)
      {
        for(int j=0;j<3;j++)
          std::cout << weight[i][j]<<"\t";
        std::cout<<std::endl;
      }
    */
 
    result[0] = p0[0][0];
    result[1] = m_zmax * p0[2][0] / m_tsampling + p0[1][0];
    result[2] = p0[2][0];
 
    if (counter == maxIter) return 3;
    return 0;
}

◆ 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_bipolarNormalization

double CscBipolarStripFitter::m_bipolarNormalization

private

Definition at line 74 of file CscBipolarStripFitter.h.

◆ m_cscCalibTool

ToolHandle<ICscCalibTool> CscBipolarStripFitter::m_cscCalibTool

private

Initial value:

{
        this,
        "cscCalibTool",
        "",
    }

Definition at line 51 of file CscBipolarStripFitter.h.

                                            {
        this,
        "cscCalibTool",
        "",
    };

◆ 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_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_n

double CscBipolarStripFitter::m_n

private

Definition at line 71 of file CscBipolarStripFitter.h.

◆ m_n2

double CscBipolarStripFitter::m_n2

private

Definition at line 72 of file CscBipolarStripFitter.h.

◆ m_phelper

const CscIdHelper* CscBipolarStripFitter::m_phelper

private

Definition at line 48 of file CscBipolarStripFitter.h.

◆ m_qerr

double CscBipolarStripFitter::m_qerr

private

Definition at line 66 of file CscBipolarStripFitter.h.

◆ m_qerr_fail

double CscBipolarStripFitter::m_qerr_fail

private

Definition at line 68 of file CscBipolarStripFitter.h.

◆ m_qerrprop

double CscBipolarStripFitter::m_qerrprop

private

Definition at line 70 of file CscBipolarStripFitter.h.

◆ m_terr

double CscBipolarStripFitter::m_terr

private

Definition at line 67 of file CscBipolarStripFitter.h.

◆ m_terr_fail

double CscBipolarStripFitter::m_terr_fail

private

Definition at line 69 of file CscBipolarStripFitter.h.

◆ m_tsampling

double CscBipolarStripFitter::m_tsampling

private

Definition at line 75 of file CscBipolarStripFitter.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.

◆ m_zmax

double CscBipolarStripFitter::m_zmax

private

Definition at line 73 of file CscBipolarStripFitter.h.

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

Public Types

Public Member Functions

Static Public Member Functions

Protected Member Functions

Private Types

Private Member Functions

Static Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

◆ ChargeList

◆ StoreGateSvc_t

Constructor & Destructor Documentation

◆ CscBipolarStripFitter()

◆ ~CscBipolarStripFitter()

Member Function Documentation

◆ declareGaudiProperty()

◆ declareProperty()

◆ Derivative()

◆ detStore()

◆ evtStore()

◆ extraDeps_update_handler()

◆ FindInitValues()

◆ FindPow()

◆ fit() [1/3]

◆ fit() [2/3]

◆ fit() [3/3]

◆ initialize()

◆ inputHandles()

◆ interfaceID()

◆ InvertMatrix()

◆ InvertSymmetric4x4()

◆ msg()

◆ msgLvl()

◆ outputHandles()

◆ renounce()

◆ renounceArray()

◆ sysInitialize()

◆ sysStart()

◆ TheFitter()

◆ updateVHKA()

Member Data Documentation

◆ m_bipolarNormalization

◆ m_cscCalibTool

◆ m_detStore

◆ m_evtStore

◆ m_n

◆ m_n2

◆ m_phelper

◆ m_qerr

◆ m_qerr_fail

◆ m_qerrprop

◆ m_terr

◆ m_terr_fail

◆ m_tsampling

◆ m_varHandleArraysDeclared

◆ m_vhka

◆ m_zmax