ZMassConstraint::ConstraintFit Class Reference

#include <ConstraintFit.h>

Inheritance diagram for ZMassConstraint::ConstraintFit:

Collaboration diagram for ZMassConstraint::ConstraintFit:

Public Member Functions
	ConstraintFit (const std::string &name)
	Create a proper constructor for Athena. More...
	~ConstraintFit (void)
StatusCode	initialize ()
	Initialize constraint fit. More...
StatusCode	doMassFit (const ConstraintFitInput &input, ConstraintFitOutput &output)
	Perform the constrained mass fit. More...
double	getMassError (const ConstraintFitInput &firstInput, const ConstraintFitInput &secondInput=ConstraintFitInput())
	Calculate the mass error without fit - use just the inputs. More...
double	getMassError (const ConstraintFitOutput &fitOutput, const ConstraintFitInput &extraInput=ConstraintFitInput())
	Calculate the mass error with the result of fit and possible additional input (4-vecs/cov mat) that has NOT been fit. More...
double	getMassError (const ConstraintFitOutput &fitOutput, const ConstraintFitOutput &secondFitOutput)
	Calculate the mass error with the result of two mass fits More...
void	addParticle (const xAOD::Muon &part, ConstraintFitInput &input, MassConstraintMuonType muonType=isCombMCMT)
	Add muon to input, must provide the resolution Scale Factor. More...
void	addParticle (const xAOD::Electron &part, float elEnergyRes, ConstraintFitInput &input)
	Add electron to input, must provide the electron energy resolution. More...
void	addFSRParticle (const xAOD::IParticle &part, const TLorentzVector &fsr4Vec, ConstraintFitInput &input)
	Add in FSR photon to input, (energy resolution is obtain in method) More...
virtual void	print () const =0
	Print the state of the tool. More...
virtual void	print () const
	Print the state of the tool. 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

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
enum	TrackParameters {   d0 = 0, z0 = 1, phi0 = 2, theta = 3,   qOverP = 4, qP = 4, x = 0, y = 1,   z = 2 }
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
int	massFitInterface (const ConstraintFitInput &theInput)
void	massFitRun (ConstraintFitOutput &output, double zresol=-1.)
double	getMassError (const std::vector< TLorentzVector > &particles, const std::vector< AmgMatrix(3, 3)> &covariances)
double	likelihoodMass2 (void)
double	likelihoodMass (double)
void	setIgnore (bool val)
bool	doSanityChecksOnCovariance (const TLorentzVector &vector, const AmgMatrix(5, 5)&covar) const
void	convertCovd0z0PhiThetaPToXYZ (const TLorentzVector &fourVec, const AmgMatrix(5, 5)&covard0z0PhiThetaP, AmgMatrix(5, 5)&covarXYZ) const
void	convertCovXYZTod0z0PhiThetaP (const std::vector< TLorentzVector > &particleList, const Amg::MatrixX &covarXYZ, Amg::MatrixX &covard0z0PhiThetaP) const
double	massFitCalculation (const Amg::MatrixX &var, double mass, Amg::MatrixX &p0)
double	massFit (const Amg::MatrixX &, const Amg::MatrixX &var, double mass, Amg::MatrixX &pOut, Amg::MatrixX &)
float	retrieve_eta_calo (const xAOD::IParticle &part) const
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
bool	m_conHasWidth
double	m_conMass
double	m_conWidth
double	m_resolution
bool	m_ignoreInputChecks
ToolHandle< CP::IEgammaCalibrationAndSmearingTool >	m_energyRescaler
ToolHandle< CP::IMuonCalibrationAndSmearingTool >	m_mu_resolSFTool
ConstraintFitInput	m_theInput
unsigned int	m_parameters
unsigned int	m_nobj
std::vector< double >	m_objmass
Amg::MatrixX	m_parametersInit
Amg::MatrixX	m_covarianceInit
Amg::MatrixX	m_parametersFinal
Amg::MatrixX	m_covarianceFinal
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 30 of file ConstraintFit.h.

Member Typedef Documentation

StoreGateSvc_t

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

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Member Enumeration Documentation

TrackParameters

enum ZMassConstraint::ConstraintFit::TrackParameters

private

Enumerator
d0
z0
phi0
theta
qOverP
qP
x
y
z

Definition at line 82 of file ConstraintFit.h.

                             {
            d0     = 0,
            z0     = 1,
            phi0   = 2,
            theta  = 3,
            qOverP = 4,
            qP     = 4,
            x      = 0,
            y      = 1,
            z      = 2
        };

Constructor & Destructor Documentation

ConstraintFit()

ZMassConstraint::ConstraintFit::ConstraintFit

(

const std::string &

name

)

Create a proper constructor for Athena.

Create a constructor for standalone usage

Definition at line 13 of file ConstraintFit.cxx.

                                                       :
        asg::AsgTool( name ),
        m_conHasWidth(true),
        m_conMass(91187.6),
        m_conWidth(2495.2),
        m_resolution(0.01),
        m_ignoreInputChecks(false),
        m_parameters(3)
    {
        declareProperty( "Z_pdg_mass",                     m_conMass);
        declareProperty( "Z_pdg_width",                    m_conWidth);
        declareProperty( "hasWidth",                       m_conHasWidth);
        declareProperty( "ignoreInputChecks",              m_ignoreInputChecks);
        declareProperty( "EgammaCalibAndSmearingTool",     m_energyRescaler = ToolHandle<CP::IEgammaCalibrationAndSmearingTool>("CP::EgammaCalibrationAndSmearingTool"));
        declareProperty( "MuonCalibrationAndSmearingTool", m_mu_resolSFTool = ToolHandle<CP::IMuonCalibrationAndSmearingTool>("CP::MuonCalibrationAndSmearingTool"));
    }

~ConstraintFit()

ZMassConstraint::ConstraintFit::~ConstraintFit

(

void

)

Definition at line 30 of file ConstraintFit.cxx.

    {}

Member Function Documentation

addFSRParticle()

void ZMassConstraint::ConstraintFit::addFSRParticle ( const xAOD::IParticle &  part, const TLorentzVector &  fsr4Vec, ConstraintFitInput &  input  )

virtual

Add in FSR photon to input, (energy resolution is obtain in method)

Implements ZMassConstraint::IConstraintFit.

Definition at line 83 of file ConstraintFit.cxx.

    {
 
        ATH_MSG_DEBUG ("addFSRParticle: *** 4vecFsr: " << fsr4Vec.Pt() << "/" << fsr4Vec.Eta() << "/" << fsr4Vec.Phi());
        // Get energy resolution - done differently for photons and electrons
        float e_res = 0;
        PATCore::ParticleType::Type phType     = PATCore::ParticleType::Electron;
        float                       cl_etaCalo = retrieve_eta_calo(part);
        if (!m_energyRescaler.empty()) {
            if (xAOD::Type::Photon ==  part.type()) {
                const xAOD::Photon* fsrPf = dynamic_cast<const xAOD::Photon*>(&part);
                phType = (xAOD::EgammaHelpers::isConvertedPhoton(fsrPf)) ?
                    PATCore::ParticleType::ConvertedPhoton : 
                    PATCore::ParticleType::UnconvertedPhoton;
            }
            e_res = m_energyRescaler->resolution( fsr4Vec.E(), fsr4Vec.Eta(), cl_etaCalo, phType)*fsr4Vec.E();
        }
        else {
            ATH_MSG_ERROR( "addFSRParticle:  - cannot find EgammaCalibrationAndSmearingTool. ");
            ATH_MSG_ERROR( "addFSRParticle:  - please set property 'EgammaCalibAndSmearingToolName' with the name of your EgammaCalibrationAndSmearingTool - default is 'EgammaCalibrationAndSmearingTool'. ");
        }
    
        AmgMatrix(5,5) photonCovarianceMatrix;
        photonCovarianceMatrix.setZero();
        photonCovarianceMatrix(d0,d0)         = 0.000001;
        photonCovarianceMatrix(z0,z0)         = 0.000001;
        photonCovarianceMatrix(phi0,phi0)     = 0.000001;
        photonCovarianceMatrix(theta,theta)   = 0.000001;
        photonCovarianceMatrix(qOverP,qOverP) = e_res*e_res;
        photonCovarianceMatrix(d0,z0)         = 0;
        photonCovarianceMatrix(d0,phi0)       = 0;
        photonCovarianceMatrix(d0,theta)      = 0;
        photonCovarianceMatrix(d0,qOverP)     = 0;
        photonCovarianceMatrix(z0,phi0)       = 0;
        photonCovarianceMatrix(z0,theta)      = 0;
        photonCovarianceMatrix(z0,qOverP)     = 0;
        photonCovarianceMatrix(phi0,theta)    = 0;
        photonCovarianceMatrix(phi0,qOverP)   = 0;
        photonCovarianceMatrix(theta,qOverP)  = 0;
 
        ATH_MSG_DEBUG ("addFSRParticle: FSR calib type,e, e_res  " << phType << " " << fsr4Vec.E() << " " << e_res);
 
        // symmetrize it
        for (int ii = 0; ii < 5; ii++)
            for (int jj = ii + 1; jj < 5; jj++)
                photonCovarianceMatrix(jj, ii) = photonCovarianceMatrix(ii, jj);
 
        bool isOK = doSanityChecksOnCovariance(fsr4Vec, photonCovarianceMatrix);
        // Get covariance in cartesian CS
        AmgMatrix(5,5) covarXYZ;
        covarXYZ.setZero();
        convertCovd0z0PhiThetaPToXYZ(fsr4Vec, photonCovarianceMatrix, covarXYZ);
        // Save in input
        input.addConstituent_FourVector_d0z0PhiThetaP(fsr4Vec, photonCovarianceMatrix, covarXYZ, isOK);
 
        ATH_MSG_DEBUG( "addFSRParticle:  - added fsr");
 
    }

addParticle() [1/2]

void ZMassConstraint::ConstraintFit::addParticle ( const xAOD::Electron &  part, float  elEnergyRes, ConstraintFitInput &  input  )

virtual

Add electron to input, must provide the electron energy resolution.

Implements ZMassConstraint::IConstraintFit.

Definition at line 310 of file ConstraintFit.cxx.

    {
        float el_E_res = elEnergyRes;
        // if (!m_energyRescaler.empty()) {
        //     // Use el.e() which is assumed to be equal to the corrected cluster energy
        //     el_E_res = m_energyRescaler->resolution(el.caloCluster()->e(), el.caloCluster()->eta())*
        //         el.caloCluster()->e();
        // }
        // else ATH_MSG_ERROR( "addParticle: - energyRescaler is NOT AVAILABLE. ");
 
        const xAOD::TrackParticle* track = el.trackParticle();
 
        if (!track) {
            ATH_MSG_ERROR( "addParticle: - ERROR could not get track particle for electron");
            return;
        }
 
        ATH_MSG_DEBUG ("addParticle: *** 4vecEl: " << el.pt() << "/" << el.eta() << "/" << el.phi());
 
        AmgMatrix(5,5) covmatrix;
        covmatrix.setZero();
        covmatrix(d0,d0)         = track->definingParametersCovMatrix()(d0,d0);
        covmatrix(z0,z0)         = track->definingParametersCovMatrix()(z0,z0);
        covmatrix(phi0,phi0)     = track->definingParametersCovMatrix()(phi0,phi0);
        covmatrix(theta,theta)   = track->definingParametersCovMatrix()(theta,theta);
        covmatrix(qOverP,qOverP) = el_E_res*el_E_res;
        covmatrix(d0,z0)         = track->definingParametersCovMatrix()(d0,z0);
        covmatrix(d0,phi0)       = track->definingParametersCovMatrix()(d0,phi0);
        covmatrix(d0,theta)      = track->definingParametersCovMatrix()(d0,theta);
        covmatrix(d0,qOverP)     = 0;
        covmatrix(z0,phi0)       = track->definingParametersCovMatrix()(z0,phi0);
        covmatrix(z0,theta)      = track->definingParametersCovMatrix()(z0,theta);
        covmatrix(z0,qOverP)     = 0;
        covmatrix(phi0,theta)    = track->definingParametersCovMatrix()(phi0,theta);
        covmatrix(phi0,qOverP)   = 0;
        covmatrix(theta,qOverP)  = 0;
 
        ATH_MSG_DEBUG( "addParticle: el p,cov" 
                       << " " << el_E_res << " " << el.e()
                       << " " << track->definingParametersCovMatrix()(d0,    d0)
                       << " " << track->definingParametersCovMatrix()(z0,    z0)
                       << " " << track->definingParametersCovMatrix()(phi0,  phi0)
                       << " " << track->definingParametersCovMatrix()(theta, theta)
                       << " " << track->definingParametersCovMatrix()(d0,    z0)
                       << " " << track->definingParametersCovMatrix()(d0,    phi0)
                       << " " << track->definingParametersCovMatrix()(d0,    theta)
                       << " " << track->definingParametersCovMatrix()(z0,    phi0)
                       << " " << track->definingParametersCovMatrix()(z0,    theta)
                       << " " << track->definingParametersCovMatrix()(phi0,  theta));
        
 
 
        for(int ii=0;ii<5;ii++)
            for(int jj=ii+1;jj<5;jj++)
                covmatrix(jj,ii) = covmatrix(ii,jj);
 
        bool isOK = doSanityChecksOnCovariance(el.p4(), covmatrix);
        // Get covariance in cartesian CS
        AmgMatrix(5,5) covarXYZ;
        covarXYZ.setZero();
        convertCovd0z0PhiThetaPToXYZ(el.p4(), covmatrix, covarXYZ);
        // Save in input
        input.addConstituent_FourVector_d0z0PhiThetaP(el.p4(), covmatrix, covarXYZ, isOK);
 
        ATH_MSG_DEBUG( "addParticle: - el \n" << covmatrix);
 
        ATH_MSG_VERBOSE( "addParticle: - el\n" << input.getCovarianceCartesian(input.getNConstituents()-1));
 
        ATH_MSG_VERBOSE( "addParticle: - added electron");
    }

addParticle() [2/2]

void ZMassConstraint::ConstraintFit::addParticle ( const xAOD::Muon &  part, ConstraintFitInput &  input, MassConstraintMuonType  muonType = isCombMCMT  )

virtual

Add muon to input, must provide the resolution Scale Factor.

Implements ZMassConstraint::IConstraintFit.

Definition at line 145 of file ConstraintFit.cxx.

    {
 
        // For muons, we allow the mass constraint to be done in one of three ways:
        //  1)   standard way - use muon momentum and primary track for cov matrix, or
        //  2,3) ID or MS     - use the ID or MS tracks for both the momentum and cov matrix
 
        TLorentzVector mu4vec;  // to be filled according to muon type
 
        // setup accessors for muon decorations - set in MuonCalibrationAndSmearingTool
        static const SG::AuxElement::ConstAccessor<float> muonSpectrometerPt ("MuonSpectrometerPt");
        static const SG::AuxElement::ConstAccessor<float> innerDetectorPt    ("InnerDetectorPt");
 
        // For the momentum, 
        // Get the track particle according to the requested muon type for the covariance matrix
        const xAOD::TrackParticle* track = mu.primaryTrackParticle();
        bool set4vec = false;
        if (((isMS_MCMT == muonType) || (isID_MCMT == muonType)) && xAOD::Muon::Combined == mu.muonType()) {
            if (isMS_MCMT  == muonType) { 
                track = mu.trackParticle(xAOD::Muon::MuonSpectrometerTrackParticle);
                if (!muonSpectrometerPt.isAvailable(mu))
                    ATH_MSG_ERROR( "addParticle:  - could not get muonSpectrometerPt from muon. Please applyCorrection with the MuonCalibAndSmearTool");
                if (!track) {
                    ATH_MSG_ERROR( "addParticle:  - Combined muon is missing MS track particle. Using combined track particle");
                    track = mu.primaryTrackParticle();
                }
                mu4vec.SetPtEtaPhiM(muonSpectrometerPt(mu), track->eta(), track->phi(), mu.m());
                set4vec = true;
                ATH_MSG_DEBUG ("addParticle: set muon track to MS ");
            }
            else {
                track = mu.trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
                if (!innerDetectorPt.isAvailable(mu))
                    ATH_MSG_ERROR( "addParticle:  - could not get innerDetectorPt from muon. Please applyCorrection with the MuonCalibAndSmearTool");
                mu4vec.SetPtEtaPhiM(innerDetectorPt(mu), track->eta(), track->phi(), mu.m());
                set4vec = true;
                ATH_MSG_DEBUG ("addParticle: set muon track to ID ");
            }
        }
 
        if (!set4vec) mu4vec.SetPtEtaPhiM(mu.pt(), mu.eta(), mu.phi(), mu.m());
        
        if (!track) {
            ATH_MSG_ERROR( "addParticle:  - could not get track particle for muon");
            return;
        }
 
        ATH_MSG_DEBUG ("addParticle: *** 4vecMu: " << mu4vec.Pt() << "/" << mu4vec.Eta() << "/" << mu4vec.Phi());
 
 
        float muSF = 1;
        int   type = -1;      
 
        // Comment out use of MuonCalibrationAndSmearingTool until the SF is implemented by Giacomo
        
        // if (m_mu_resolSFTool) {
        //     // Get scale factor for the muon covariance qoverp resolution.
        //     // For muonType == isCombMCMT, we get a SF differently according to the type of muon: Comb/ST/CT/SA
        //     //     muonType == isID_MCMT,  we use the type 2 for ID tracks
        //     //     muonType == isMS_MCMT,  we use the type 3 for MS tracks
 
        //     if (muonType == isCombMCMT) {
        //         if (xAOD::Muon::Combined            == mu.muonType()) type = 1;
        //         else if (xAOD::Muon::SegmentTagged  == mu.muonType()) type = 2;
        //         else if (xAOD::Muon::CaloTagged     == mu.muonType()) type = 2;
        //         else if (xAOD::Muon::MuonStandAlone == mu.muonType()) type = 3;
        //         else ATH_MSG_ERROR("addParticle: - **** MUON type is not CB/ST/SA/Calo **** ");
        //     }
        //     else if (muonType == isID_MCMT) {
        //         type = 2;
        //         if (xAOD::Muon::MuonStandAlone == mu.muonType()) type = 3;
        //     }
        //     else if (muonType == isMS_MCMT) {
        //         type = 3;
        //         if (xAOD::Muon::SegmentTagged  == mu.muonType() ||
        //             xAOD::Muon::CaloTagged     == mu.muonType()) type = 2;
        //     }
 
        //     if (muonType < isCombMCMT) {
        //         ATH_MSG_ERROR( "addParticle: - Invalid muon type requested, please use .");
        //     }
        //     else {
        //         muSF = m_mu_resolSF->getResolutionScaleFactor(mu4vec, type);
        //     }
 
        //     ATH_MSG_DEBUG ("addParticle: 4v,type,errSF  " << mu4vec.Eta() << " " << type << " " << muSF);
            
            
        // }
        // else {
        //     ATH_MSG_ERROR( "addParticle: - MuonResSFTool is NOT set. Please call setMuonResSFTool with MuonResolutionAndMomentumScaleFactors tool.");
        // }
 
        
        AmgMatrix(5,5) covmatrix;
        covmatrix.setZero();
        covmatrix(d0,d0)         = track->definingParametersCovMatrix()(d0,d0);
        covmatrix(z0,z0)         = track->definingParametersCovMatrix()(z0,z0);
        covmatrix(phi0,phi0)     = track->definingParametersCovMatrix()(phi0,phi0);
        covmatrix(theta,theta)   = track->definingParametersCovMatrix()(theta,theta);
        covmatrix(qOverP,qOverP) = track->definingParametersCovMatrix()(qOverP,qOverP)*muSF*muSF;
        covmatrix(d0,z0)         = track->definingParametersCovMatrix()(d0,z0);
        covmatrix(d0,phi0)       = track->definingParametersCovMatrix()(d0,phi0);
        covmatrix(d0,theta)      = track->definingParametersCovMatrix()(d0,theta);
        covmatrix(d0,qOverP)     = track->definingParametersCovMatrix()(d0,qOverP)*muSF;
        covmatrix(z0,phi0)       = track->definingParametersCovMatrix()(z0,phi0);
        covmatrix(z0,theta)      = track->definingParametersCovMatrix()(z0,theta);
        covmatrix(z0,qOverP)     = track->definingParametersCovMatrix()(z0,qOverP)*muSF;
        covmatrix(phi0,theta)    = track->definingParametersCovMatrix()(phi0,theta);
        covmatrix(phi0,qOverP)   = track->definingParametersCovMatrix()(phi0,qOverP)*muSF;
        covmatrix(theta,qOverP)  = track->definingParametersCovMatrix()(theta,qOverP)*muSF;
 
        ATH_MSG_DEBUG ("addParticle: type,errSF,mu p,cov  " << type << " " << muSF
                       << " " << mu4vec.E()
                       << " " << track->definingParametersCovMatrix()(d0,    d0)
                       << " " << track->definingParametersCovMatrix()(z0,    z0)
                       << " " << track->definingParametersCovMatrix()(phi0,  phi0)
                       << " " << track->definingParametersCovMatrix()(theta, theta)
                       << " " << track->definingParametersCovMatrix()(qOverP,qOverP)*muSF*muSF
                       << " " << track->definingParametersCovMatrix()(d0,    z0)
                       << " " << track->definingParametersCovMatrix()(d0,    phi0)
                       << " " << track->definingParametersCovMatrix()(d0,    theta)
                       << " " << track->definingParametersCovMatrix()(d0,qOverP)*muSF
                       << " " << track->definingParametersCovMatrix()(z0,    phi0)
                       << " " << track->definingParametersCovMatrix()(z0,    theta)
                       << " " << track->definingParametersCovMatrix()(z0,qOverP)*muSF
                       << " " << track->definingParametersCovMatrix()(phi0,  theta)
                       << " " << track->definingParametersCovMatrix()(phi0,qOverP)*muSF
                       << " " << track->definingParametersCovMatrix()(theta,qOverP)*muSF);
 
        for (int ii = 0; ii < 5; ii++)
            for (int jj = ii+1; jj < 5; jj++)
                covmatrix(jj,ii) = covmatrix(ii,jj);
 
        ATH_MSG_VERBOSE( "addParticle: cov matrix \n" << covmatrix );
 
        // Transform cov of q/p to cov p in covariance matrix
        double P     = mu4vec.P();
        AmgMatrix(5,5) jacobian0;
        jacobian0.setZero();
        jacobian0(0,0) = 1.;
        jacobian0(1,1) = 1.;
        jacobian0(2,2) = 1.;
        jacobian0(3,3) = 1.;
        jacobian0(4,4) = -P*P;
 
        AmgMatrix(5,5) newmatrix = jacobian0.transpose() * covmatrix * jacobian0;
 
        ATH_MSG_DEBUG( "addParticle: new matrix \n" << newmatrix );
 
        bool isOK = doSanityChecksOnCovariance(mu4vec, newmatrix);
        // Get covariance in cartesian CS
        AmgMatrix(5,5) covarXYZ;
        covarXYZ.setZero();
        convertCovd0z0PhiThetaPToXYZ(mu4vec, newmatrix, covarXYZ);
        // Save in input
        input.addConstituent_FourVector_d0z0PhiThetaP(mu4vec, newmatrix, covarXYZ, isOK);
 
        ATH_MSG_VERBOSE( "addParticle: - added muon");
 
    }

convertCovd0z0PhiThetaPToXYZ()

void ZMassConstraint::ConstraintFit::convertCovd0z0PhiThetaPToXYZ ( const TLorentzVector &  fourVec, const AmgMatrix(5, 5)&  covard0z0PhiThetaP, AmgMatrix(5, 5)&  covarXYZ  ) const

private

Definition at line 1127 of file ConstraintFit.cxx.

    {
        //going from d0,z0,phi,theta,P --> d0,z0,px,py,pz
        double phi   = fourVec.Phi();
        double theta = fourVec.Theta();
        double P     = fourVec.P();
        //std::cout << ":::::L \t"<<P<<"\t"<< phi<<"\t"<<theta<<std::endl;
        AmgMatrix(5,5) jacobian;
        jacobian.setZero();
        jacobian(0,0) = 1.;
        jacobian(1,1) = 1.;
        jacobian(2,2) = -P * TMath::Sin(theta) * TMath::Sin(phi);
        jacobian(2,3) =  P * TMath::Sin(theta) * TMath::Cos(phi);
        jacobian(3,2) =  P * TMath::Cos(theta) * TMath::Cos(phi);
        jacobian(3,3) =  P * TMath::Cos(theta) * TMath::Sin(phi);
        jacobian(3,4) = -P * TMath::Sin(theta);
        jacobian(4,2) =      TMath::Sin(theta) * TMath::Cos(phi);
        jacobian(4,3) =      TMath::Sin(theta) * TMath::Sin(phi);
        jacobian(4,4) =      TMath::Cos(theta);
    
        covarXYZ = jacobian.transpose() * covard0z0PhiThetaP * jacobian;
 
        // std::cout << "initial\n"<< covar << "final \n" << newcovariance << "jac\n"
        //         << jacobian << std::endl;
 
    }

convertCovXYZTod0z0PhiThetaP()

void ZMassConstraint::ConstraintFit::convertCovXYZTod0z0PhiThetaP ( const std::vector< TLorentzVector > &  particleList, const Amg::MatrixX &  covarXYZ, Amg::MatrixX &  covard0z0PhiThetaP  ) const

private

build the Jacobian of the phi,theta,P --> px,py,pz

Definition at line 1158 of file ConstraintFit.cxx.

    {
        unsigned int matrixSize = 5 * particleList.size();
        Amg::MatrixX Jacobian(matrixSize, matrixSize);
        Jacobian.setZero();
        for(unsigned int ii = 0; ii < particleList.size(); ii++) {
            double phi   = particleList.at(ii).Phi();
            double theta = particleList.at(ii).Theta();
            double P     = particleList.at(ii).P();
            Jacobian(    5*ii,     5*ii) = 1.;
            Jacobian(1 + 5*ii, 1 + 5*ii) = 1.;
 
            Jacobian(2 + 5*ii, 2 + 5*ii) = -P * TMath::Sin(theta) * TMath::Sin(phi);
            Jacobian(2 + 5*ii, 3 + 5*ii) =  P * TMath::Sin(theta) * TMath::Cos(phi);
            Jacobian(3 + 5*ii, 2 + 5*ii) =  P * TMath::Cos(theta) * TMath::Cos(phi);
            Jacobian(3 + 5*ii, 3 + 5*ii) =  P * TMath::Cos(theta) * TMath::Sin(phi);
            Jacobian(3 + 5*ii, 4 + 5*ii) = -P * TMath::Sin(theta);
            Jacobian(4 + 5*ii, 2 + 5*ii) =      TMath::Sin(theta) * TMath::Cos(phi);
            Jacobian(4 + 5*ii, 3 + 5*ii) =      TMath::Sin(theta) * TMath::Sin(phi);
            Jacobian(4 + 5*ii, 4 + 5*ii) =      TMath::Cos(theta);
        }
 
        Amg::MatrixX Jacobianinverse(matrixSize, matrixSize);
        Jacobianinverse=Jacobian.inverse();
        covard0z0PhiThetaP = Jacobianinverse.transpose() * covarXYZ * Jacobianinverse;
    }

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; }

doMassFit()

StatusCode ZMassConstraint::ConstraintFit::doMassFit ( const ConstraintFitInput &  input, ConstraintFitOutput &  output  )

virtual

Perform the constrained mass fit.

Implements ZMassConstraint::IConstraintFit.

Definition at line 54 of file ConstraintFit.cxx.

    {
 
        // For the moment, call the 'old' methods
        if (massFitInterface(input)) {
            ATH_MSG_ERROR( "Could not add in input" );
            return StatusCode::FAILURE;
        }
    
        massFitRun(output);
    
        // Return gracefully:
        return StatusCode::SUCCESS;
    }

doSanityChecksOnCovariance()

bool ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance ( const TLorentzVector &  vector, const AmgMatrix(5, 5)&  covar  ) const

private

Definition at line 1090 of file ConstraintFit.cxx.

    {
        bool isOK = true;
        if ( sqrt(covar(qP,qP))/vector.P()>1.0 ) { 
            isOK = false;
            ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Fractional uncertainty on P = "
                            << sqrt(covar(qP,qP))/vector.P() << " is > 1 ...this is not ok...");
            
        }
 
        if( covar(theta,theta) > 1.e-2) {
            isOK = false;
            ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Uncertainty on Theta = "
                            << covar(theta,theta) << " is > 10^-2 ...this is not ok");
        }
        if( covar(phi0,phi0) > 1.e-2) {                                 
            isOK = false;
            ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Uncertainty on Phi = "
                            << covar(phi0,phi0) << " is > 10^-2 ...this is not ok");
        }
 
        for(int i=0;i<5;i++) {
            for(int j=0;j<5;j++) {
                if(i==j) continue;
                if(fabs(covar(i,j))/sqrt(covar(i,i))/sqrt(covar(j,j)) > 1.) {
                    ATH_MSG_WARNING( "ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Off-diagonal term "
                                     << i << " " << j << " is > 1 - doesn't look ok");
                    isOK = false;
                }
            }
        }
 
        return isOK;
    }

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

getKey()

SG::sgkey_t asg::AsgTool::getKey ( const void *  ptr ) const

inherited

Get the (hashed) key of an object that is in the event store.

This is a bit of a special one. StoreGateSvc and xAOD::TEvent both provide ways for getting the SG::sgkey_t key for an object that is in the store, based on a bare pointer. But they provide different interfaces for doing so.

In order to allow tools to efficiently perform this operation, they can use this helper function.

See also

asg::AsgTool::getName

Parameters

ptr	The bare pointer to the object that the event store should know about

Returns

The hashed key of the object in the store. If not found, an invalid (zero) key.

Definition at line 119 of file AsgTool.cxx.

                                                    {
 
#ifdef XAOD_STANDALONE
      // In case we use @c xAOD::TEvent, we have a direct function call
      // for this.
      return evtStore()->event()->getKey( ptr );
#else
      const SG::DataProxy* proxy = evtStore()->proxy( ptr );
      return ( proxy == nullptr ? 0 : proxy->sgkey() );
#endif // XAOD_STANDALONE
   }

getMassError() [1/4]

double ZMassConstraint::ConstraintFit::getMassError ( const ConstraintFitInput &  firstInput, const ConstraintFitInput &  secondInput = ConstraintFitInput()  )

virtual

Calculate the mass error without fit - use just the inputs.

Implements ZMassConstraint::IConstraintFit.

Definition at line 905 of file ConstraintFit.cxx.

    {
        unsigned int nPartFirst  = firstInput.getNConstituents();
        unsigned int nPartSecond = secondInput.getNConstituents();
        unsigned int nPart = nPartFirst + nPartSecond;
        std::vector<TLorentzVector> particles(nPart);
        std::vector<AmgMatrix(3,3)> covariances(nPart);
        for ( unsigned int iobj = 0; iobj < firstInput.getNConstituents(); ++iobj) {
            covariances[iobj].setZero();
            particles[iobj] = firstInput.getConstituentFourVector(iobj);
            firstInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
        }
        for ( unsigned int iobj = 0; iobj < secondInput.getNConstituents(); ++iobj) {
            covariances[iobj + nPartFirst].setZero();
            particles[iobj + nPartFirst] = secondInput.getConstituentFourVector(iobj);
            secondInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
        }
        return getMassError( particles, covariances );
    }

getMassError() [2/4]

double ZMassConstraint::ConstraintFit::getMassError ( const ConstraintFitOutput &  fitOutput, const ConstraintFitInput &  extraInput = ConstraintFitInput()  )

virtual

Calculate the mass error with the result of fit and possible additional input (4-vecs/cov mat) that has NOT been fit.

Implements ZMassConstraint::IConstraintFit.

Definition at line 929 of file ConstraintFit.cxx.

    {
        unsigned int nPartFirst  = fitOutput.getNConstituents();
        unsigned int nPartSecond = extraInput.getNConstituents();
        unsigned int nPart = nPartFirst + nPartSecond;
        std::vector<TLorentzVector> particles(nPart);
        std::vector<AmgMatrix(3,3)> covariances(nPart);
        for ( unsigned int iobj = 0; iobj < fitOutput.getNConstituents(); ++iobj) {
            covariances[iobj].setZero();
            particles[iobj] = fitOutput.getConstituentFourVector(iobj);
            fitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
        }
        for ( unsigned int iobj = 0; iobj < extraInput.getNConstituents(); ++iobj) {
            covariances[iobj + nPartFirst].setZero();
            particles[iobj + nPartFirst] = extraInput.getConstituentFourVector(iobj);
            extraInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
        }
        return getMassError( particles, covariances );
    }

getMassError() [3/4]

double ZMassConstraint::ConstraintFit::getMassError ( const ConstraintFitOutput &  fitOutput, const ConstraintFitOutput &  secondFitOutput  )

virtual

Calculate the mass error with the result of two mass fits

Implements ZMassConstraint::IConstraintFit.

Definition at line 954 of file ConstraintFit.cxx.

    {
        unsigned int nPartFirst  = fitOutput.getNConstituents();
        unsigned int nPartSecond = secondFitOutput.getNConstituents();
        unsigned int nPart = nPartFirst + nPartSecond;
        std::vector<TLorentzVector> particles(nPart);
        std::vector<AmgMatrix(3,3)> covariances(nPart);
        for ( unsigned int iobj = 0; iobj < fitOutput.getNConstituents(); ++iobj) {
            covariances[iobj].setZero();
            particles[iobj] = fitOutput.getConstituentFourVector(iobj);
            fitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
        }
        for ( unsigned int iobj = 0; iobj < secondFitOutput.getNConstituents(); ++iobj) {
            covariances[iobj + nPartFirst].setZero();
            particles[iobj + nPartFirst] = secondFitOutput.getConstituentFourVector(iobj);
            secondFitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
        }
        return getMassError( particles, covariances );
    }

getMassError() [4/4]

double ZMassConstraint::ConstraintFit::getMassError ( const std::vector< TLorentzVector > &  particles, const std::vector< AmgMatrix(3, 3)> &  covariances  )

private

Definition at line 976 of file ConstraintFit.cxx.

    {
        double masserror=0;
 
        if (particles.size() != covariances.size()) {
            ATH_MSG_ERROR( "getMassError: Number of particles is not equal to the number of covariance matrices");
            return 1E11;
        }
 
        ATH_MSG_DEBUG ("getMassError: 1 " << particles.size() << "/" << covariances.size());
        
        //composite lorentz vector and its invariant mass
        TLorentzVector combv;
        for ( auto lv : particles ) combv += lv;
        double invmass = combv.M();
 
        ATH_MSG_VERBOSE ("getMassError: 1.1 invmass " << invmass);
 
        //Calculation of 1x4 jacobian - used in derivation of var(M)
        Amg::MatrixX jacobianNPtoM(1,4);
        jacobianNPtoM.setZero();
        jacobianNPtoM(0,0) = -1.*combv.Px()/invmass;  // dM/dpxN
        jacobianNPtoM(0,1) = -1.*combv.Py()/invmass;  // dM/dpyN
        jacobianNPtoM(0,2) = -1.*combv.Pz()/invmass;  // dM/dpzN
        jacobianNPtoM(0,3) =     combv.E()/invmass;   // dM/dEN
 
        ATH_MSG_VERBOSE ("getMassError: 2\n" << jacobianNPtoM);
 
        int iobj = -1;
        for ( auto lv : particles ) {
            ++iobj;
 
            ATH_MSG_VERBOSE ("getMassError: 1.2 iobj " << iobj);
            
            // Constituent 3*3 covariance matrix (Phi, Theta, P) 
            const AmgMatrix(3,3)& covPhiThetaP = covariances[iobj];
 
            ATH_MSG_VERBOSE ("getMassError: 1.3 covPhiThetaP\n " << covPhiThetaP);
 
            double theta = lv.Theta();
            double phi   = lv.Phi();
            double p     = lv.P();
            double e     = lv.E();
            double m     = lv.M(); 
 
            ATH_MSG_VERBOSE ("getMassError: 2.1 iobj " << iobj << "covPhiThetaP\n" << covPhiThetaP);
 
            // Convert (Phi, Theta, P) to (Phi, eta, P), as covPhiEtaP
            AmgMatrix(3,3) jacobian1;
            jacobian1.setZero();
            jacobian1(0,0) =  1.;                   // 
            jacobian1(1,1) = -1*(1./sin(theta));    // deta/dtheta
            jacobian1(2,2) =  1.;                   // 
            AmgMatrix(3,3) covPhiEtaP = jacobian1 * covPhiThetaP * jacobian1.transpose();
 
            ATH_MSG_VERBOSE ("getMassError: 2.2 jacobian1\n" << jacobian1);
            ATH_MSG_VERBOSE ("getMassError: 2.4 covPhiEtaP\n" << covPhiEtaP);
 
            
            // Rearrange things to get (E, eta, phi, M),  as covEEtaPhiM
            AmgMatrix(4,4) covEEtaPhiM;
            covEEtaPhiM.setZero();
            for(int i = 0; i < 3; i++) {
                for(int j = 0; j < 3; j++) {
                    covEEtaPhiM(i,j) = covEEtaPhiM(j,i) = covPhiEtaP(2-i,2-j);
                }
            } 
    
            ATH_MSG_VERBOSE ("getMassError: 3 covEEtaPhiM\n" << covEEtaPhiM);
 
            //convert (E, eta, phi, M) to (px,py,pz,E) representation, as covPxPyPzE
            AmgMatrix(4,4) jacobian2;
            jacobian2.setZero();
            jacobian2(0,0) =  e/p * sin(theta)*cos(phi); // dpx/dE
            jacobian2(1,0) =  e/p * sin(theta)*sin(phi); // dpy/dE
            jacobian2(2,0) =  e/p * cos(theta);          // dpy/dE 
            jacobian2(3,0) =  1;                     // dE/dE
    
            jacobian2(0,1) =  -p * cos(phi) * cos(theta)*sin(theta); // dpx/deta
            jacobian2(1,1) =  -p * sin(phi) * cos(theta)*sin(theta); // dpy/deta
            jacobian2(2,1) =   p * sin(theta)*sin(theta);            // dpz/deta;
 
            jacobian2(0,2) =  -p * sin(theta) * sin(phi);  // dpx/dphi
            jacobian2(1,2) =   p * sin(theta) * cos(phi);  // dpy/dphi
    
            jacobian2(0,3) =  -m/p * sin(theta)*cos(phi); // dpx/dM
            jacobian2(1,3) =  -m/p * sin(theta)*sin(phi); // dpy/dM
            jacobian2(2,3) =  -m/p * cos(theta);          // dpz/dM
 
            AmgMatrix(4,4) covPxPyPzE = jacobian2 * covEEtaPhiM * jacobian2.transpose();
 
            ATH_MSG_VERBOSE ("getMassError: 4 jacobian2\n" << jacobian2);
            ATH_MSG_VERBOSE ("getMassError: 5 covPxPyPzE\n" << covPxPyPzE);
 
            //Get mass variance
            // (Similarity transform: C = A*B*A^T (from egammaFourMomentumError/GeneralUtils.cxx))
            Amg::MatrixX  em = jacobianNPtoM * covPxPyPzE * jacobianNPtoM.transpose(); 
            masserror += em(0,0);
 
            ATH_MSG_VERBOSE ("getMassError: 6 em\n" << em);
 
        }
 
        //Return square root of variance.
        if (masserror < 0.) {
            ATH_MSG_WARNING( "getMassError: Mass covariance element less than zero! Returning 1E11 ...");
            masserror = 1E11;
        }
        return sqrt(masserror);
 
    }

getName()

const std::string & asg::AsgTool::getName ( const void *  ptr ) const

inherited

Get the name of an object that is / should be in the event store.

This is a bit of a special one. StoreGateSvc and xAOD::TEvent both provide ways for getting the std::string name for an object that is in the store, based on a bare pointer. But they provide different interfaces for doing so.

In order to allow tools to efficiently perform this operation, they can use this helper function.

See also

asg::AsgTool::getKey

Parameters

ptr	The bare pointer to the object that the event store should know about

Returns

The string name of the object in the store. If not found, an empty string.

Definition at line 106 of file AsgTool.cxx.

                                                            {
 
#ifdef XAOD_STANDALONE
      // In case we use @c xAOD::TEvent, we have a direct function call
      // for this.
      return evtStore()->event()->getName( ptr );
#else
      const SG::DataProxy* proxy = evtStore()->proxy( ptr );
      static const std::string dummy = "";
      return ( proxy == nullptr ? dummy : proxy->name() );
#endif // XAOD_STANDALONE
   }

getProperty()

template<class T >

const T* asg::AsgTool::getProperty ( const std::string &  name ) const

inherited

Get one of the tool's properties.

initialize()

StatusCode ZMassConstraint::ConstraintFit::initialize ( )

virtual

Initialize constraint fit.

Reimplemented from asg::AsgTool.

Definition at line 34 of file ConstraintFit.cxx.

    {
        ATH_MSG_INFO ("initialize: retrieve tools" << m_energyRescaler.name() << " and "
                       << m_mu_resolSFTool.name());
 
        if (!m_energyRescaler.retrieve().isSuccess()) {
            ATH_MSG_ERROR ("initialize: unable to retrieve EgammaCalibrationAndSmearingTool");
            return StatusCode::FAILURE;
        }
 
        if (!m_mu_resolSFTool.retrieve().isSuccess()) {
            ATH_MSG_ERROR ("initialize: unable to retrieve MuonCalibrationAndSmearingTool");
            return StatusCode::FAILURE;
        }
        
        // Return gracefully:
        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.

likelihoodMass()

double ZMassConstraint::ConstraintFit::likelihoodMass ( double  MassResol )

private

Definition at line 434 of file ConstraintFit.cxx.

    {
        const int rows =4;
        const int columns = m_nobj;
        double **p = 0;
        p = new double *[rows] ;
        for( int i = 0 ; i < rows ; i++ )
            p[i] = new double[columns];
 
        for(unsigned int iobj=0; iobj<m_nobj; iobj++)
            {
                p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
                for(int i=0;i<3;i++)
                    {
                        p[i][iobj] = m_parametersInit(i + 3*iobj, 0);
                        p[3][iobj] += p[i][iobj]*p[i][iobj];
                    }
                p[3][iobj] = sqrt(p[3][iobj]);
            }
        double Etot = 0.;
        double Xtot = 0.;
        double Ytot = 0.;
        double Ztot = 0.;
        for(unsigned int ij=0;ij<m_nobj;ij++)
            {
                Etot += p[3][ij];
                Xtot += p[0][ij];
                Ytot += p[1][ij];
                Ztot += p[2][ij];
            }
        for( int i = 0 ; i < rows ; i++ )
            delete [] p[i];
        delete [] p;
 
        double initMass =  Etot * Etot - Xtot * Xtot - Ytot * Ytot - Ztot * Ztot;
        initMass = sqrt(initMass);
 
        double sig = MassResol;
        double xLeft = initMass;
        double xRight = m_conMass;
        if (m_conMass < initMass) {
            xLeft = m_conMass;
            xRight = initMass;
        }
        double dLinitL = (initMass - xLeft) / sig / sig - 2.*(xLeft - m_conMass) / ((xLeft - m_conMass) * (xLeft - m_conMass) + m_conWidth * m_conWidth);
        double dLinitR = (initMass - xRight) / sig / sig - 2.*(xRight - m_conMass) / ((xRight - m_conMass) * (xRight - m_conMass) + m_conWidth * m_conWidth);
        if (dLinitL * dLinitR < 0.) {
            while (xRight - xLeft > 1.) { //1 MeV
                double xM = (xRight + xLeft) / 2.;
                double dL = (initMass - xM) / sig / sig - 2.*(xM - m_conMass) / ((xM - m_conMass) * (xM - m_conMass) + m_conWidth * m_conWidth);
                if (dL * dLinitL < 0.) {
                    xRight = xM;
                    dLinitR = dL;
                } else {
                    xLeft = xM;
                    dLinitL = dL;
                }
            }
            return (xLeft + xRight) / 2.;
        } else {
            if (dLinitL > dLinitR)
                return xLeft;
            else
                return xRight;
        }
    }

likelihoodMass2()

double ZMassConstraint::ConstraintFit::likelihoodMass2 ( void  )

private

Definition at line 502 of file ConstraintFit.cxx.

    {
        ATH_MSG_VERBOSE( "likelihoodMass2: entering" );
        const int rows =4;
        const int columns = m_nobj;
        double **p = 0;
        p = new double *[rows] ;
        for( int i = 0 ; i < rows ; i++ )
            p[i] = new double[columns];
 
        for(unsigned int iobj=0; iobj<m_nobj; iobj++) {
            p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
            for(int i=0;i<3;i++)                 {
                p[i][iobj] = m_parametersInit(i + 3*iobj, 0);
                p[3][iobj] += p[i][iobj]*p[i][iobj];
            }
            p[3][iobj] = sqrt(p[3][iobj]);
        }
        double etot = 0.;
        double xtot = 0.;
        double ytot = 0.;
        double ztot = 0.;
        for(unsigned int ij=0;ij<m_nobj;ij++) {
            etot += p[3][ij];
            xtot += p[0][ij];
            ytot += p[1][ij];
            ztot += p[2][ij];
        }
    
        double initMass =  etot * etot - xtot * xtot - ytot * ytot - ztot * ztot;
        initMass = sqrt(initMass);
        //what is the uncertainty on initMass?
        Amg::MatrixX jacobianMass(3*m_nobj, 1);
        jacobianMass.setZero();
        for(unsigned int ii=0;ii<m_nobj;ii++) {
            jacobianMass(3*ii + 0, 0) = (2.*(etot) * p[0][ii] / p[3][ii] - 2.* xtot);
            jacobianMass(3*ii + 1, 0) = (2.*(etot) * p[1][ii] / p[3][ii] - 2.* ytot);
            jacobianMass(3*ii + 2, 0) = (2.*(etot) * p[2][ii] / p[3][ii] - 2.* ztot);
        }
        ATH_MSG_VERBOSE( "likelihoodMass2: the jacobian mass/matrix \n" << jacobianMass << "\n" << m_covarianceInit );
        
 
        Amg::MatrixX sigm = jacobianMass.transpose() * m_covarianceInit * jacobianMass;
 
        ATH_MSG_VERBOSE( "likelihoodMass2: sigm \n" << sigm );
 
        double sig = sigm(0, 0);
        sig = sqrt(sig);
        m_resolution = sig;
        sig /= 2 * initMass;
        double MassResol = sig;
        double maxmass = initMass;
 
        ATH_MSG_VERBOSE( "likelihoodMass2: 2 " );
 
        double max = -(maxmass - initMass) * (maxmass - initMass) / 2. / MassResol / MassResol 
            - log((maxmass * maxmass - m_conMass * m_conMass) * (maxmass * maxmass - m_conMass * m_conMass) + m_conMass * m_conMass * m_conWidth * m_conWidth);
 
        for (int i = 1; i < 401; i++) {
            double ytest = initMass + (m_conMass - initMass) / 400 * i;
 
            double val = -(ytest - initMass) * (ytest - initMass) / 2. / MassResol / MassResol 
                - log((ytest * ytest - m_conMass * m_conMass) * (ytest * ytest - m_conMass * m_conMass) + m_conMass * m_conMass * m_conWidth * m_conWidth);
 
            if (val > max) {
                max = val;
                maxmass = ytest;
            }
        }
  
        ATH_MSG_VERBOSE( "likelihoodMass2: 3 " );
  
        for( int i = 0 ; i < rows ; i++ )
            delete [] p[i];
        delete [] p;
 
        ATH_MSG_VERBOSE( "likelihoodMass2: maxmass " << maxmass );
 
        return maxmass;
    }

massFit()

double ZMassConstraint::ConstraintFit::massFit ( const Amg::MatrixX &  , const Amg::MatrixX &  var, double  mass, Amg::MatrixX &  pOut, Amg::MatrixX &    )

private

Definition at line 856 of file ConstraintFit.cxx.

    {
        Amg::MatrixX ivar = var.inverse();
 
        bool doIter = true;
        int  maxIterations = 20;
        int  iIter  = 0;
        double constraintValue = -1.e10;
        ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose() 
                       << " Iterations " << iIter
                       << " parameters \n" << m_parametersFinal
                       << " cov matrix \n" << m_covarianceFinal
                       << " Iteration  " << iIter);
        while (doIter) {
            Amg::MatrixX p_old(pOut);
            constraintValue = massFitCalculation(var, mass, pOut); //call the fit
 
            // Convergence criteria
            // 1. the parameters should not change too much (<1.e-6 relative change)
            // 2. the constraint should be satisfied very well (<1.e-6 absolute)
            double maxDiff = 5.e15;
            Amg::MatrixX diff = (pOut - p_old);
            for (unsigned int i = 0; i < m_parameters * m_nobj; i++) {
                diff(i, 0) = diff(i, 0) / p_old(i, 0);
                if (maxDiff > diff(i, 0))
                    maxDiff = diff(i, 0);
            }
            if ((maxDiff < 1.e-4 && fabs(constraintValue) < 5.e-5) || maxIterations <= iIter)
                doIter = false;
            ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose() 
                           << " Iterations " << iIter
                           << " parameters \n" << m_parametersFinal
                           << " cov matrix \n" << m_covarianceFinal
                           << " \n Iteration  " << iIter << " doIter " << doIter);
            iIter++;
 
        }
        ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose() 
                       << " Iterations " << iIter
                       << " parameters \n" << m_parametersFinal
                       << " cov matrix \n" << m_covarianceFinal
                       << " Iteration  " << iIter);
        // double chi2 = 0.;//calculateChi2(p0, var, mass);
        return 0;
    }

massFitCalculation()

double ZMassConstraint::ConstraintFit::massFitCalculation ( const Amg::MatrixX &  var, double  mass, Amg::MatrixX &  p0  )

private

Definition at line 786 of file ConstraintFit.cxx.

    {
        const int rows =4;
        const int columns = m_nobj;
 
        double **p = 0;
        p = new double *[rows] ;
        for( int i = 0 ; i < rows ; i++ )
            p[i] = new double[columns];
 
        for(unsigned int iobj=0; iobj<m_nobj; iobj++) {
            p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
            for(int i=0;i<3;i++)
                {
                    p[i][iobj] = p0(i+3*iobj, 0);
                    p[3][iobj] += p[i][iobj]*p[i][iobj];
                }
            p[3][iobj] = sqrt(p[3][iobj]);
        }
 
        double Etot = 0.;
        double Xtot = 0.;
        double Ytot = 0.;
        double Ztot = 0.;
        for(unsigned int ij=0;ij<m_nobj;ij++)
            {
                Etot += p[3][ij];
                Xtot += p[0][ij];
                Ytot += p[1][ij];
                Ztot += p[2][ij];
            }
        double XYZtot[3] = {Xtot, Ytot, Ztot};
        std::vector<double> constraintD;
        for(unsigned int ij=0;ij<m_nobj;ij++)
            for(int ik=0;ik<3;ik++)
                constraintD.push_back((2.*(Etot) * p[ik][ij] / p[3][ij] - 2.* XYZtot[ik]));
 
        double constraintd =  Etot * Etot - Xtot * Xtot - Ytot * Ytot - Ztot * Ztot;
        constraintd = constraintd - mass * mass;
        //std::cout << "constraint = " <<constraintd << " mass value " <<mass<<std::endl;
        Amg::MatrixX D(1,3*m_nobj);
        for(unsigned int ij=0;ij<3*m_nobj;ij++)
            D(0, ij) = constraintD.at(ij);
        AmgMatrix(1,1) d;
        d(0, 0) = constraintd;
 
        Amg::MatrixX DVD(D * var * D.transpose());
 
        Amg::MatrixX DVDinverse(1,1);
        DVDinverse=DVD.inverse();
 
        Amg::MatrixX VD(DVDinverse);
        Amg::MatrixX lambda(VD * d);
        Amg::MatrixX test(p0 - var * D.transpose() * lambda);
        Amg::MatrixX Vp(var - var * D.transpose() * VD * D * var);
        for(unsigned int ij=0;ij<3*m_nobj;ij++) {
            p0(ij, 0) = test(ij, 0);
            m_parametersFinal(ij,0) = p0(ij,0);
            for(unsigned int jk = 0; jk < 3 * m_nobj; jk++)
                m_covarianceFinal(ij,jk) = Vp(ij,jk);
        }
        double constraintValue = d(0, 0);
 
        for( int i = 0 ; i < rows ; i++ )
            delete [] p[i];
        delete [] p;
        return constraintValue;
    }

massFitInterface()

int ZMassConstraint::ConstraintFit::massFitInterface ( const ConstraintFitInput &  theInput )

private

Definition at line 384 of file ConstraintFit.cxx.

    {
        m_objmass.clear();
        m_theInput = theInput;
 
        // the m_parameters first parameters of pin are the fitted parameters
        // the next one is the mass of the particle
        if(theInput.getNConstituents() <2 ) {
            ATH_MSG_ERROR( "massFitInterface: Found number of input particles less than 2");
            return 1; 
        }
        else {
            ATH_MSG_DEBUG( "massFitInterface: nparticles " << theInput.getNConstituents());
        }
        
        m_nobj            = theInput.getNConstituents();
        unsigned int dimension    = m_parameters * m_nobj;
        m_parametersInit = Amg::MatrixX(dimension, 1);
        m_parametersInit.setZero();
        m_covarianceInit = Amg::MatrixX(dimension, dimension);
        m_covarianceInit.setZero();
        m_parametersFinal= Amg::MatrixX(dimension, 1);
        m_parametersFinal.setZero();
        m_covarianceFinal= Amg::MatrixX(dimension, dimension);
        m_covarianceFinal.setZero();
 
        for (unsigned int iobj = 0; iobj < m_nobj; iobj++)
            {
                m_objmass.push_back(theInput.getConstituentFourVector(iobj).M());
                m_parametersInit(0 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Px();
                m_parametersInit(1 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Py();
                m_parametersInit(2 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Pz();
 
                for (unsigned int i = 0; i < m_parameters; i++)
                    for (unsigned int j = 0; j < m_parameters; j++)
                        m_covarianceInit(i + 3*iobj, j + 3*iobj) = (theInput.getCovarianceCartesian(iobj))(i + 2, j + 2); //keep only the lower right 3x3 part
            }
  
        for(unsigned int i=0;i<dimension;i++)
            {
                m_parametersFinal(i,0) = m_parametersInit(i,0);
                for(unsigned int j=0;j<dimension;j++)
                    m_covarianceFinal(i,j) = m_covarianceInit(i,j);
            }
        ATH_MSG_VERBOSE( "massFitInterface: parameters \n" << m_parametersInit );
        ATH_MSG_VERBOSE( "massFitInterface: cov matrix \n" << m_covarianceInit );
        return 0;
    }

massFitRun()

void ZMassConstraint::ConstraintFit::massFitRun ( ConstraintFitOutput &  output, double  zresol = -1.  )

private

Definition at line 584 of file ConstraintFit.cxx.

    {
        ATH_MSG_VERBOSE( "massFitRun: entering ");
 
        if(!m_ignoreInputChecks && !m_theInput.isOK()) {
            ATH_MSG_WARNING("ConstraintFit::massFitRun() the input is not OK.");
            ATH_MSG_WARNING("ConstraintFit::massFitRun()(cont) Following H->4l group decision (Feb 1, 2013), this event will not be fit");
            ATH_MSG_WARNING("ConstraintFit::massFitRun()(cont) To bypass this behaviour set the ingore flag");
 
            std::vector<TLorentzVector> particleList;
            for(unsigned int i=0;i<m_nobj;i++)
                {
                    TLorentzVector particle;
                    particle.SetXYZM(m_parametersInit(0 + i*m_parameters,0),
                                     m_parametersInit(1 + i*m_parameters,0),
                                     m_parametersInit(2 + i*m_parameters,0),
                                     m_objmass[i]);
                    particleList.push_back(particle);
      
                }
            Amg::MatrixX covOut(5*m_nobj,5*m_nobj);
            covOut.setZero();
            for(unsigned int iobj=0;iobj<m_nobj;iobj++)
                for(int ii=0;ii<5;ii++)
                    for(int jj=0;jj<5;jj++)
                        covOut(ii+5*iobj,jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii,jj);
 
            ATH_MSG_VERBOSE( "massFitRun: cov 1 ");
 
            // convert covariance matrix from cartesian to d0z0PhiThetaP
            Amg::MatrixX covard0z0PhiThetaP;
            covard0z0PhiThetaP.setZero();
            convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
            // save output
            output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
 
            // ATH_MSG_VERBOSE( "massFitRun: cov 1 " << covOut);
 
            return;
        }
 
        ATH_MSG_VERBOSE( "massFitRun: 1 ");
 
        Amg::MatrixX parametersFit = m_parametersInit;
        double mass = -1.;
        if (!m_conHasWidth) 
            mass = m_conMass;
        else if(zresol<0.)
            mass = likelihoodMass2();
        else
            mass = likelihoodMass(zresol);
 
        ATH_MSG_VERBOSE( "massFitRun: resolution/mass " << m_resolution << "/" << mass);
 
        // do not use chi2:
        // double chi2 = massFit(&m_parametersInit, &m_covarianceInit, mass, &parametersFit, &m_covarianceFinal);
        massFit(m_parametersInit, m_covarianceInit, mass, parametersFit, m_covarianceFinal);
 
        ATH_MSG_VERBOSE( "massFitRun: 2");
 
        std::vector<TLorentzVector> particleList;
        for(unsigned int i=0;i<m_nobj;i++) {
            TLorentzVector particle;
            particle.SetXYZM(m_parametersFinal(0 + i*m_parameters,0),
                             m_parametersFinal(1 + i*m_parameters,0),
                             m_parametersFinal(2 + i*m_parameters,0),
                             m_objmass[i]);
            particleList.push_back(particle);
        }
 
        ATH_MSG_VERBOSE( "massFitRun: 3 ");
 
        Amg::MatrixX covOut(5*m_nobj,5*m_nobj);
        covOut.setZero();
        for(unsigned int iobj=0;iobj<m_nobj;iobj++)
            for(int ii=0;ii<2;ii++)
                for(int jj=0;jj<2;jj++)
                    covOut(ii+5*iobj,jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii,jj);
        for(unsigned int iobj=0;iobj<m_nobj;iobj++)
            for(int ii=2;ii<5;ii++)
                for(int jj=2;jj<5;jj++)
                    covOut(ii+5*iobj,jj+5*iobj) = m_covarianceFinal(ii-2+iobj*3,jj-2+iobj*3);
 
        ATH_MSG_VERBOSE( "massFitRun: 4 ");
 
        // empirical corrections for constraint width
        if(!m_conHasWidth) {
            // correlation terms
            for(unsigned int iobj1=0;iobj1<m_nobj;iobj1++) {
                for(unsigned int iobj2=0;iobj2<m_nobj;iobj2++) {
                    if(iobj1 == iobj2) continue;
                    for(int ii=2;ii<5;ii++)
                        for(int jj=2;jj<5;jj++)
                            covOut(ii + 5*iobj1, jj+5*iobj2) = m_covarianceFinal(ii-2+iobj1*3,jj-2+iobj2*3);
                }
            }
        }
 
        // ATH_MSG_DEBUG( "massFitRun: cov 2 ");
 
        // convert covariance matrix from cartesian to d0z0PhiThetaP
        Amg::MatrixX covard0z0PhiThetaP;
        covard0z0PhiThetaP.setZero();
        convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
        // save output
        output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
 
        // ATH_MSG_DEBUG( "massFitRun: cov 2 " << covOut);
 
 
        ATH_MSG_VERBOSE( "massFitRun: 5 conHasWidth " << m_conHasWidth);
 
        // ATH_MSG_DEBUG( "massFitRun: 5.5 conHasWidth " << m_conHasWidth);
 
        if(m_conHasWidth) {
            TLorentzVector aa;
            for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
                aa += m_theInput.getConstituentFourVector(iobj);
            }
            double massInit = aa.M();
 
            // ATH_MSG_DEBUG( "massFitRun: 5.6 massInit " << massInit);
 
            if(massInit<83.e3) {
 
                ATH_MSG_VERBOSE( "massFitRun: 6 \n" << covOut);
 
                for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
                    for(int ii = 0; ii < 5; ii++) {
                        for(int jj = 0; jj < 5; jj++) {
                            covOut(ii+5*iobj, jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii, jj);
                        }
                    }
                }
 
                ATH_MSG_VERBOSE( "massFitRun: 6.1 \n" << covOut);
            }
            else {
 
                // ATH_MSG_DEBUG( "massFitRun: 7 ");
 
                std::vector<double> resolInit;
                std::vector<double> resolFinal;
                std::vector<double> resolFinalnew;
                resolInit.reserve(m_nobj);
                resolFinal.reserve(m_nobj);
                resolFinalnew.reserve(m_nobj);
                double sumResolInit=0.;
                const double constraintWidthSquare = /*m_resolution*m_resolution + */m_conWidth*m_conWidth/2.35/2.35;
                for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
 
                    // ATH_MSG_DEBUG( "massFitRun: 7 iobj " << iobj << " " << m_theInput.getNConstituents());
 
                    resolInit.push_back((m_theInput.getConstituentCovariance(iobj))(4,4));
 
                    // ATH_MSG_DEBUG( "massFitRun: 7 1 ");
 
                    AmgMatrix(5,5) tmpMtx;
                    tmpMtx.setZero();
                    output.getConstituentCovarianced0z0PhiThetaP(iobj, tmpMtx);
                    resolFinal.push_back(tmpMtx(4,4));
                    // resolFinal.push_back((output.getConstituentCovarianced0z0PhiThetaP(iobj))(4,4));
 
                    // ATH_MSG_DEBUG( "massFitRun: 7 2 ");
 
                    sumResolInit += resolInit[iobj];
                    //std::cout << " ### "<< iobj <<"\t"<< resolInit[iobj] <<"\t"<<resolFinal[iobj]<<std::endl;
                }
      
 
                ATH_MSG_VERBOSE( "massFitRun: 8 \n" << covOut);
 
                for(unsigned int iobj=0;iobj<m_nobj;iobj++) {
                    resolFinalnew.push_back( resolFinal[iobj] + constraintWidthSquare * 
                                             resolInit[iobj]*resolInit[iobj]
                                             /sumResolInit/sumResolInit );
                    covOut(4+5*iobj,4+5*iobj) =  covOut(4+5*iobj,4+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
                    covOut(3+5*iobj,3+5*iobj) =  covOut(3+5*iobj,3+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
                    covOut(2+5*iobj,2+5*iobj) =  covOut(2+5*iobj,2+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
                    //std::cout << " ### "<< iobj <<"\t"<< resolInit[iobj] <<"\t"<<resolFinal[iobj]<<"\t"<< resolFinalnew[iobj]<<std::endl;
                }
                ATH_MSG_VERBOSE( "massFitRun: 8.1 \n" << covOut);
            }
 
            ATH_MSG_VERBOSE( "massFitRun: 9 ");
 
            // ATH_MSG_DEBUG( "massFitRun: cov 3 ");
 
            // convert covariance matrix from cartesian to d0z0PhiThetaP
            Amg::MatrixX covard0z0PhiThetaP;
            convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
            // save output
            output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
 
            ATH_MSG_VERBOSE( "massFitRun: cov 3 \n" << covOut);
 
        }
 
        ATH_MSG_VERBOSE( "massFitRun: end ");
    }

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);
  }

msg_level_name()

const std::string & asg::AsgTool::msg_level_name ( ) const

inherited

A deprecated function for getting the message level's name.

Instead of using this, weirdly named function, user code should get the string name of the current minimum message level (in case they really need it...), with:

MSG::name( msg().level() )

This function's name doesn't follow the ATLAS coding rules, and as such will be removed in the not too distant future.

Returns

The string name of the current minimum message level that's printed

Definition at line 101 of file AsgTool.cxx.

                                                  {
 
      return MSG::name( msg().level() );
   }

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.

print() [1/2]

void asg::AsgTool::print ( ) const

virtualinherited

Print the state of the tool.

Implements asg::IAsgTool.

Reimplemented in JetRecTool, JetFinder, JetModifiedMassDrop, JetFromPseudojet, JetReclusterer, JetReclusteringTool, JetTruthLabelingTool, JetPileupLabelingTool, HI::HIPileupTool, top::TopObjectSelection, JetDumper, JetBottomUpSoftDrop, JetRecursiveSoftDrop, JetSoftDrop, JetConstituentsRetriever, JetSubStructureMomentToolsBase, JetSplitter, JetToolRunner, JetPruner, JetPseudojetRetriever, JetTrimmer, AsgHelloTool, and KtDeltaRTool.

Definition at line 131 of file AsgTool.cxx.

                             {
 
      ATH_MSG_INFO( "AsgTool " << name() << " @ " << this );
      return;
   }

print() [2/2]

virtual void asg::IAsgTool::print ( ) const

pure virtualinherited

Print the state of the tool.

Implemented in JetRecTool, JetFinder, JetModifiedMassDrop, JetFromPseudojet, JetReclusterer, JetReclusteringTool, JetTruthLabelingTool, JetPileupLabelingTool, HI::HIPileupTool, top::TopObjectSelection, asg::AsgTool, JetDumper, JetBottomUpSoftDrop, JetRecursiveSoftDrop, JetSoftDrop, JetConstituentsRetriever, JetSubStructureMomentToolsBase, JetSplitter, JetToolRunner, JetPruner, JetPseudojetRetriever, JetTrimmer, AsgHelloTool, and KtDeltaRTool.

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();
  }

retrieve_eta_calo()

float ZMassConstraint::ConstraintFit::retrieve_eta_calo ( const xAOD::IParticle &  part ) const

private

Definition at line 1189 of file ConstraintFit.cxx.

    {
        double eta_calo = 0;
        const xAOD::Egamma* eg = dynamic_cast<const xAOD::Egamma*> (&part);
        if (!eg) {
            ATH_MSG_ERROR("retrieve_eta_calo - unable to cast to Egamma");
            return eta_calo;
        }
        const xAOD::CaloCluster& cluster = *eg->caloCluster();
        if (cluster.retrieveMoment(xAOD::CaloCluster::ETACALOFRAME,
                                   eta_calo)) { }
        else if (cluster.isAvailable<float>("etaCalo")) {
            eta_calo = cluster.auxdata<float>("etaCalo");
        }
        else {
            ATH_MSG_ERROR("retrieve_eta_calo - etaCalo not available as auxilliary variable");
            ATH_MSG_WARNING("retrieve_eta_calo - using eta as etaCalo");
            eta_calo = cluster.eta();
        }
        return eta_calo;
    }

setIgnore()

void ZMassConstraint::ConstraintFit::setIgnore ( bool  val )

inlineprivate

Definition at line 100 of file ConstraintFit.h.

{ m_ignoreInputChecks=val; }

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_conHasWidth

bool ZMassConstraint::ConstraintFit::m_conHasWidth

private

Definition at line 109 of file ConstraintFit.h.

m_conMass

double ZMassConstraint::ConstraintFit::m_conMass

private

Definition at line 110 of file ConstraintFit.h.

m_conWidth

double ZMassConstraint::ConstraintFit::m_conWidth

private

Definition at line 111 of file ConstraintFit.h.

m_covarianceFinal

Amg::MatrixX ZMassConstraint::ConstraintFit::m_covarianceFinal

private

Definition at line 126 of file ConstraintFit.h.

m_covarianceInit

Amg::MatrixX ZMassConstraint::ConstraintFit::m_covarianceInit

private

Definition at line 124 of file ConstraintFit.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_energyRescaler

ToolHandle<CP::IEgammaCalibrationAndSmearingTool> ZMassConstraint::ConstraintFit::m_energyRescaler

private

Definition at line 116 of file ConstraintFit.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_ignoreInputChecks

bool ZMassConstraint::ConstraintFit::m_ignoreInputChecks

private

Definition at line 113 of file ConstraintFit.h.

m_mu_resolSFTool

ToolHandle<CP::IMuonCalibrationAndSmearingTool> ZMassConstraint::ConstraintFit::m_mu_resolSFTool

private

Definition at line 117 of file ConstraintFit.h.

m_nobj

unsigned int ZMassConstraint::ConstraintFit::m_nobj

private

Definition at line 121 of file ConstraintFit.h.

m_objmass

std::vector<double> ZMassConstraint::ConstraintFit::m_objmass

private

Definition at line 122 of file ConstraintFit.h.

m_parameters

unsigned int ZMassConstraint::ConstraintFit::m_parameters

private

Definition at line 120 of file ConstraintFit.h.

m_parametersFinal

Amg::MatrixX ZMassConstraint::ConstraintFit::m_parametersFinal

private

Definition at line 125 of file ConstraintFit.h.

m_parametersInit

Amg::MatrixX ZMassConstraint::ConstraintFit::m_parametersInit

private

Definition at line 123 of file ConstraintFit.h.

m_resolution

double ZMassConstraint::ConstraintFit::m_resolution

private

Definition at line 112 of file ConstraintFit.h.

m_theInput

ConstraintFitInput ZMassConstraint::ConstraintFit::m_theInput

private

Definition at line 119 of file ConstraintFit.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:

• ConstraintFit.h
• ConstraintFit.cxx