Rec::MuonCaloEnergyTool Class Reference

#include <MuonCaloEnergyTool.h>

Inheritance diagram for Rec::MuonCaloEnergyTool:

Collaboration diagram for Rec::MuonCaloEnergyTool:

Public Member Functions
	MuonCaloEnergyTool (const std::string &, const std::string &, const IInterface *)
virtual	~MuonCaloEnergyTool ()=default
virtual StatusCode	initialize () override
void	calculateMuonEnergies (const Trk::Track *trk, double deltaE, double meanIoni, double sigmaIoni, double &E, double &sigma, double &E_FSR, double &E_expected, double &E_em_meas, double &E_em_exp, double &E_tile_meas, double &E_tile_exp, double &E_HEC_meas, double &E_HEC_exp, double &E_dead_exp, std::vector< Identifier > *crossedCells=0, std::vector< double > *sigmaNoise_cell=0, std::vector< double > *E_exp_cell=0) const override
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	thresholdCorrection (double E_observed, double E_expected, double sigma_Noise) 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	etaCorr (double eta)

Private Attributes
ToolHandle< Trk::IParticleCaloExtensionTool >	m_caloExtensionTool {this, "ParticleCaloExtensionTool", "", "Tool to make the step-wise extrapolation"}
ToolHandle< Rec::IParticleCaloCellAssociationTool >	m_caloCellAssociationTool {this, "ParticleCaloCellAssociationTool", "", "Tool to make the cell association"}
ToolHandle< Trk::ITrackParticleCreatorTool >	m_particleCreator {this,"TrackParticleCreator", "", "The CB Particle Creator Tool"}
SG::ReadCondHandleKey< CaloNoise >	m_caloNoiseCDOKey {this, "CaloNoiseKey", "totalNoise", "SG Key of CaloNoise data object"}
Gaudi::Property< double >	m_sigmasAboveNoise {this, "SigmasAboveNoise", 4.0 }
Gaudi::Property< double >	m_emEtCut {this, "EmEtCut", 2.5* Gaudi::Units::GeV }
Gaudi::Property< double >	m_emF1Cut {this, "EmF1Cut", 0.15 }
SG::ReadHandleKey< xAOD::TrackParticleContainer >	m_indetTrackParticleLocation
SG::ReadHandleKey< xAOD::TrackParticleContainer >	m_muonTrackParticleLocation
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

Static Private Attributes
static constexpr double	m_emipEM {0.42}
static constexpr double	m_emipTile {0.86}
static constexpr double	m_emipHEC {0.65}

Detailed Description

Definition at line 36 of file MuonCaloEnergyTool.h.

Member Typedef Documentation

StoreGateSvc_t

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

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

MuonCaloEnergyTool()

Rec::MuonCaloEnergyTool::MuonCaloEnergyTool	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p )

Definition at line 23 of file MuonCaloEnergyTool.cxx.

                                                                                                      :
        AthAlgTool(t, n, p)
    {
        declareInterface<IMuonCaloEnergyTool>(this);
    }

~MuonCaloEnergyTool()

virtual Rec::MuonCaloEnergyTool::~MuonCaloEnergyTool ( )

virtualdefault

Member Function Documentation

calculateMuonEnergies()

void Rec::MuonCaloEnergyTool::calculateMuonEnergies ( const Trk::Track * trk, double deltaE, double meanIoni, double sigmaIoni, double & E, double & sigma, double & E_FSR, double & E_expected, double & E_em_meas, double & E_em_exp, double & E_tile_meas, double & E_tile_exp, double & E_HEC_meas, double & E_HEC_exp, double & E_dead_exp, std::vector< Identifier > * crossedCells = 0, std::vector< double > * sigmaNoise_cell = 0, std::vector< double > * E_exp_cell = 0 ) const

overridevirtual

Implements Rec::IMuonCaloEnergyTool.

Definition at line 42 of file MuonCaloEnergyTool.cxx.

                                                                                                                            {
        const EventContext& ctx = Gaudi::Hive::currentContext();
        //
        // Input parameters trk:        (muon) track pointer
        //                  deltaE:     Mean Energy loss in Calorimeter
        //                  meanIoni:   Mean ionization Energy loss in Calorimeter
        //                  sigmaIoni:  sigma of the ionization Energy loss in Calorimeter
        //
        // Ouput parameters E:          Calorimeter energy measured (corrected for dead material etc.)
        //                  sigma:      Error on measured energy
        //                  E_FSR:      Energy of Final State Radiation: e.m. = photon energy (if above ET and F1 cut)
        //                              if EFSR>0 the corrected hadronic energy is used for the measured Energy E
        //                  E_expected: expected meanIonization Energyloss from the TrackingGeometry
 
        E = 0.;
        sigma = 0.;
        E_FSR = 0.;
        E_expected = 0.;
 
        E_em_meas = 0.;
        E_em_exp = 0.;
        E_tile_meas = 0.;
        E_tile_exp = 0.;
        E_HEC_meas = 0.;
        E_HEC_exp = 0.;
        E_dead_exp = 0.;
 
        bool storeCells = false;
        if (crossedCells != nullptr && sigmaNoise_cell != nullptr && E_exp_cell != nullptr) {
            storeCells = true;
            crossedCells->clear();
            sigmaNoise_cell->clear();
            E_exp_cell->clear();
        }
        //
        // For the expected Eloss in the dead or not instrumented material we will use the meanIonization loss
        //     this meanIoni is stored on the extended Eloss object
        //     the sigmaIoni corresponds to the Landau width
        //     a factor 3.59524 gives the sigma of the upper tail
        //
        // The Eloss derived from the momentum differences on the calo extension correponds to
        //     deltaE on the Eloss object = mean Eloss including ionization and radiation
        //
        // To go to the  meanIonization loss one needs a scale factor scale_Ionization
        //
        if (!trk) return;
        if (!trk->perigeeParameters()) return;
        if (trk->perigeeParameters()->parameters()[Trk::qOverP] == 0.) return;
 
        //    int isign = 1;
        //    if(meanIoni<0) isign = -1;
        //    double mopIoni =  meanIoni - isign*3.59524*std::abs(sigmaIoni);
        double mopIoni = meanIoni;
 
        double scale_Ionization = 0.9;
        if (std::abs(deltaE) > 0 && std::abs(meanIoni) > 0) scale_Ionization = mopIoni / deltaE;
        double error_ratio = 0.3;
        if (std::abs(meanIoni) > 0 && sigmaIoni > 0) error_ratio = 3.59524 * std::abs(sigmaIoni / meanIoni);
 
        ATH_MSG_DEBUG(" Inputs deltaE " << deltaE << " meanIoni " << meanIoni << " sigmaIoni " << sigmaIoni << " scale_Ionization "
                                        << scale_Ionization << " error_ratio " << error_ratio);
        if (deltaE == 0 || meanIoni == 0 || sigmaIoni < 0)
            ATH_MSG_WARNING(" Strange Inputs deltaE " << deltaE << " meanIoni " << meanIoni << " sigmaIoni " << sigmaIoni);
 
        // associate muon to calorimeter
 
        const xAOD::TrackParticle* tp = nullptr;
 
        SG::ReadHandle<xAOD::TrackParticleContainer> indetTrackParticles(m_indetTrackParticleLocation,ctx);
        if (indetTrackParticles.isValid()) {
            // check ID trackparticles
            for (const auto *it : *indetTrackParticles) {
                if ((*it).track() == trk) {
                    tp = &(*it);
                    break;
                }
            }
            if (tp) ATH_MSG_DEBUG(" ID xAOD::TrackParticle found " << tp);
        } else {
            ATH_MSG_WARNING("Failed to retrieve " << m_indetTrackParticleLocation.key());
        }
 
        // look for Muon trackparticles
 
        SG::ReadHandle<xAOD::TrackParticleContainer> muonTrackParticles(m_muonTrackParticleLocation,ctx);
        if (!tp && muonTrackParticles.isValid()) {
            for (const auto *it : *muonTrackParticles) {
                if ((*it).track() == trk) {
                    tp = &(*it);
                    break;
                }
            }
            if (tp) ATH_MSG_DEBUG(" Muon xAOD::TrackParticle found " << tp);
        } else if (!tp) {
            ATH_MSG_WARNING("Failed to retrieve " << m_muonTrackParticleLocation.key());
        }
 
        std::unique_ptr<const xAOD::TrackParticle> tpholder;
        if (!tp) {
            tpholder = std::unique_ptr<const xAOD::TrackParticle>(m_particleCreator->createParticle(*trk, nullptr, nullptr, xAOD::muon));
 
            tp = tpholder.get();
            if (tp) ATH_MSG_DEBUG(" xAOD::TrackParticle created from scratch " << tp);
        }
 
        if (!tp) return;
 
        std::unique_ptr<const Rec::ParticleCellAssociation> association =
            m_caloCellAssociationTool->particleCellAssociation(*tp, 0.2, nullptr);
 
        if (!association) { return; }
        ATH_MSG_DEBUG(" particleCellAssociation done  " << association.get());
 
        std::vector<std::pair<const CaloCell*, Rec::ParticleCellIntersection*> > cellIntersections = association->cellIntersections();
 
        const Trk::CaloExtension& caloExtension = association->caloExtension();
 
        if (!caloExtension.caloEntryLayerIntersection()) {
            ATH_MSG_DEBUG("calculateMuonEnergies() - No caloEntryLayerIntersection found");
            return;
        }
 
        ATH_MSG_DEBUG(" nr of cell intersections " << cellIntersections.size());
        if (cellIntersections.size() < 3) ATH_MSG_DEBUG(" Strange nr of calorimeter cell intersections " << cellIntersections.size());
 
        double theta = caloExtension.caloEntryLayerIntersection()->position().theta();
 
        double Eloss = 0.;
        double scaleTG = 1.0;
        if (!caloExtension.muonEntryLayerIntersection()) {
            if (caloExtension.caloEntryLayerIntersection()->momentum().mag() - std::abs(deltaE) > 1000.) {
                ATH_MSG_WARNING(
                    " No muonEntryLayerIntersection found and therefore the expected Eloss is not calculated properly for momentum "
                    << caloExtension.caloEntryLayerIntersection()->momentum().mag());
            }
            ATH_MSG_DEBUG(" No muonEntryLayerIntersection found and therefore the expected Eloss is not calculated properly ");
        } else {
            Eloss =
                caloExtension.caloEntryLayerIntersection()->momentum().mag() - caloExtension.muonEntryLayerIntersection()->momentum().mag();
            ATH_MSG_DEBUG(" Energy loss from CaloExtension "
                          << Eloss << " R muon Entry " << caloExtension.muonEntryLayerIntersection()->position().perp() << " Z muon Entry "
                          << caloExtension.muonEntryLayerIntersection()->position().z() << " R calo entry "
                          << caloExtension.caloEntryLayerIntersection()->position().perp() << " Z calo entry "
                          << caloExtension.caloEntryLayerIntersection()->position().z());
            if (Eloss > 25000. && Eloss > 0.1 * caloExtension.muonEntryLayerIntersection()->momentum().mag()) {
                ATH_MSG_WARNING(" Crazy Energy loss from CaloExtension "
                                << Eloss << " p at CaloEntry " << caloExtension.caloEntryLayerIntersection()->momentum().mag()
                                << " p at Muon Entry " << caloExtension.muonEntryLayerIntersection()->momentum().mag());
                scaleTG = mopIoni / Eloss;
            }
        }
 
        CaloExtensionHelpers::EntryExitLayerMap entryExitLayerMap;
        CaloExtensionHelpers::entryExitLayerMap(caloExtension, entryExitLayerMap);
        ATH_MSG_DEBUG("EntryExitLayerMap " << entryExitLayerMap.size());
 
        CaloExtensionHelpers::ScalarLayerMap eLossLayerMap;
        CaloExtensionHelpers::eLossLayerMap(caloExtension, eLossLayerMap);
        ATH_MSG_DEBUG("eLossLayerMap " << eLossLayerMap.size());
 
        double scale_em_expected = 0.97;    // 0.73
        double scale_tile_expected = 1.17;  // 0.89
        double scale_HEC_expected = 1.11;   // 0.84
                                            //
                                            // TG expectations
                                            //
        double EE_emB = 0.;
        double EE_emE = 0.;
        double EE_HEC = 0.;
        double EE_tile = 0.;
        //
        //    const char* sampname[24] = {
        //      "PreSamplerB", "EMB1", "EMB2", "EMB3",
        //      "PreSamplerE", "EME1", "EME2", "EME3",
        //      "HEC0", "HEC1","HEC2", "HEC3",
        //      "TileBar0", "TileBar1", "TileBar2",
        //      "TileExt0", "TileExt1", "TileExt2",
        //      "TileGap1", "TileGap2", "TileGap3",
        //      "FCAL0", "FCAL1", "FCAL2"};
        //
        for (int i = CaloSampling::PreSamplerB; i < CaloSampling::CaloSample::Unknown; i++) {
            CaloSampling::CaloSample sample = static_cast<CaloSampling::CaloSample>(i);
            auto pos = entryExitLayerMap.find(sample);
            if (pos == entryExitLayerMap.end()) continue;
            auto eLossPair = eLossLayerMap.find(sample);
            double eLoss = 0.;
            if (eLossPair != eLossLayerMap.end()) {
                eLoss = eLossPair->second;
                ATH_MSG_DEBUG(" sample " << i << " eLoss " << scale_Ionization * eLoss);
                if (i <= 3) {
                    EE_emB += scale_em_expected * scale_Ionization * eLoss;
                } else if (i <= 7) {
                    EE_emE += scale_em_expected * scale_Ionization * eLoss;
                } else if (i <= 11) {
                    EE_HEC += scale_HEC_expected * scale_Ionization * eLoss;
                } else if (i <= 20) {
                    EE_tile += scale_tile_expected * scale_Ionization * eLoss;
                }
            }
        }
        if (caloExtension.caloEntryLayerIntersection() && caloExtension.muonEntryLayerIntersection()) {
            double Eloss =
                caloExtension.caloEntryLayerIntersection()->momentum().mag() - caloExtension.muonEntryLayerIntersection()->momentum().mag();
            ATH_MSG_DEBUG(" eta " << caloExtension.caloEntryLayerIntersection()->momentum().eta() << " expected energies from TG EE_emB "
                                  << EE_emB << " EE_emE " << EE_emE << " EE_tile " << EE_tile << " EE_HEC " << EE_HEC << " total Eloss "
                                  << Eloss);
        }
        //
        // do not continue with crazy values for the Eloss from TG
        //
        if (scaleTG != 1.0) return;
 
        // measured and expected energies
 
        // Get Calo-Noise CDO:
        SG::ReadCondHandle<CaloNoise> caloNoiseHdl{m_caloNoiseCDOKey};
        const CaloNoise* caloNoise = *caloNoiseHdl;
 
        double E_em = 0.;
        double E_em_expected = 0.;
        double E_em_exptot = 0.;
        double E_em_expall = 0.;
        double E_em_threshold = 0.;
        int nlay_em = 0;
        double sigma_Noise_em = 0.;
 
        double E_HEC = 0.;
        double E_HEC_expected = 0.;
        double E_HEC_threshold = 0.;
        double E_HEC_exptot = 0.;
        double E_HEC_expall = 0.;
        int nlay_HEC = 0;
        double sigma_Noise_HEC = 0.;
 
        double E_tile = 0.;
        double E_tile_expected = 0.;
        double E_tile_threshold = 0.;
        double E_tile_exptot = 0.;
        double E_tile_expall = 0.;
        int nlay_tile = 0;
        double sigma_Noise_tile = 0.;
 
        E_expected = 0.;
        const Amg::Vector3D&  extension_pos = caloExtension.caloEntryLayerIntersection()->position();
        const double thetaPos = caloExtension.caloEntryLayerIntersection()->position().theta();
 
        for (auto it : cellIntersections) {
            const CaloCell* curr_cell = it.first;
            const Amg::Vector3D cell_vec(curr_cell->x(),curr_cell->y(), curr_cell->z());
            const double cell_perp = cell_vec.perp();
            int cellSampling = curr_cell->caloDDE()->getSampling();
            bool badCell = curr_cell->badcell();
            double cellEn = curr_cell->energy();
            Identifier id = curr_cell->ID();
 
            double f_exp = (it.second)->pathLength();
            double E_exp = (it.second)->expectedEnergyLoss();
 
            if (f_exp > 1.) f_exp = 1.;
 
            //      if(f_exp<0.1) f_exp = 0.1;
 
            double sigma_Noise = caloNoise->getEffectiveSigma(id, curr_cell->gain(), cellEn);
 
           
            double dtheta = cell_vec.theta() - thetaPos;
            double dphi = cell_vec.deltaPhi(extension_pos);
            double celldr = curr_cell->caloDDE()->dr();
            double celldz = curr_cell->caloDDE()->dz();
            double celldtheta = celldr / (cell_perp ? cell_perp : 1. );
            double celldphi = curr_cell->caloDDE()->dphi();
 
            ATH_MSG_DEBUG(" cell sampling " << cellSampling << " energy " << cellEn << " position R "
                                            << cell_perp << " z "
                                            << curr_cell->z() << " f_exp " << f_exp << " E_exp " << E_exp << " dtheta " << dtheta
                                            << " dphi " << dphi << " cell dtheta " << celldtheta << " cell dr " << celldr << " cell dz "
                                            << celldz << " cell dphi " << celldphi);
            //      ATH_MSG_DEBUG( " sigma_Noise totalNoiseRMS from Alan " << sigma_NoiseA << " sigma_Noise getEffectiveSigma " <<
            //      sigma_Noise);
 
            //   Use expected energy if measured energy below noise level
            // - this will bias the energy measurement towards too high values
            // - the bias is corrected in the thresholdCorrection function below
            //
            // sum measured, expected energies for crossed cells after noise cuts
            //
            if (cellSampling == CaloSampling::PreSamplerB || cellSampling == CaloSampling::PreSamplerE) {
                if (f_exp > 0 && cellEn > m_sigmasAboveNoise * sigma_Noise && !badCell) {
                    if (storeCells) {
                        crossedCells->push_back(id);
                        sigmaNoise_cell->push_back(sigma_Noise);
                        E_exp_cell->push_back(scale_Ionization * f_exp * E_exp *
                                              scale_em_expected);  // I don't want sum, but value for each cell
                    }
                }
            }
            if (curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::LAREM) {
                E_em_exptot += scale_Ionization * f_exp * E_exp * scale_em_expected;
                if (f_exp > 0) E_em_expall += scale_Ionization * E_exp * scale_em_expected;
                if (f_exp > 0 && cellEn > m_sigmasAboveNoise * sigma_Noise && !badCell) {
                    //        if(f_exp>0&&cellEn>m_sigmasAboveNoise*sigma_Noise&&!badCell) {
                    E_em += cellEn;
                    E_em_expected += scale_Ionization * f_exp * E_exp * scale_em_expected;
                    ATH_MSG_VERBOSE(" EM cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
                    if (storeCells) {
                        crossedCells->push_back(id);
                        sigmaNoise_cell->push_back(sigma_Noise);
                        E_exp_cell->push_back(scale_Ionization * f_exp * E_exp * scale_em_expected);
                    }
                }
                if (f_exp > 0 && !badCell) nlay_em++;
                if (f_exp > 0 && !badCell) sigma_Noise_em += sigma_Noise;
            } else if (curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::TILE) {
                E_tile_exptot += scale_Ionization * f_exp * E_exp * scale_tile_expected;
                if (f_exp > 0) E_tile_expall += scale_Ionization * E_exp * scale_tile_expected;
                if (f_exp > 0 && cellEn > m_sigmasAboveNoise * sigma_Noise && !badCell) {
                    //                       &&f_exp*E_exp>m_sigmasAboveNoise*sigma_Noise/4) {
                    E_tile += cellEn;
                    E_tile_expected += scale_Ionization * f_exp * E_exp * scale_tile_expected;
                    ATH_MSG_VERBOSE(" Tile cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
                    if (storeCells) {
                        crossedCells->push_back(id);
                        sigmaNoise_cell->push_back(sigma_Noise);
                        E_exp_cell->push_back(scale_Ionization * f_exp * E_exp * scale_tile_expected);
                    }
                }
                if (f_exp > 0 && !badCell) nlay_tile++;
                if (f_exp > 0 && !badCell) sigma_Noise_tile += sigma_Noise;
            } else if (curr_cell->caloDDE()->getSubCalo() == CaloCell_ID::LARHEC) {
                E_HEC_exptot += scale_Ionization * f_exp * E_exp * scale_HEC_expected;
                if (f_exp > 0) E_HEC_expall += scale_Ionization * E_exp * scale_HEC_expected;
                // lower threshold for HEC
                if (f_exp > 0 && cellEn > m_sigmasAboveNoise * sigma_Noise / 2. && !badCell) {
                    //                        &&f_exp*E_exp>m_sigmasAboveNoise*sigma_Noise/4) {
                    E_HEC += cellEn;
                    E_HEC_expected += scale_Ionization * f_exp * E_exp * scale_HEC_expected;
                    ATH_MSG_VERBOSE(" HEC cell " << cellEn << " sigma_Noise " << sigma_Noise << " f_exp " << f_exp << " E_exp " << E_exp);
                    if (storeCells) {
                        crossedCells->push_back(id);
                        sigmaNoise_cell->push_back(sigma_Noise);
                        E_exp_cell->push_back(scale_Ionization * f_exp * E_exp * scale_HEC_expected);
                    }
                }
                if (f_exp > 0 && !badCell) nlay_HEC++;
                if (f_exp > 0 && !badCell) sigma_Noise_HEC += sigma_Noise;
            }
            E_expected += scale_Ionization * f_exp * E_exp;
        }
        ATH_MSG_DEBUG(" After cuts measured Energies em " << E_em << " tile " << E_tile << " HEC " << E_HEC);
        ATH_MSG_DEBUG(" After cuts expected Energies em " << E_em_expected << " tile " << E_tile_expected << " HEC " << E_HEC_expected);
        ATH_MSG_DEBUG(" No cuts with length expected Energies em " << E_em_exptot << " tile " << E_tile_exptot << " HEC " << E_HEC_exptot);
        ATH_MSG_DEBUG(" No cuts NO length expected Energies em " << E_em_expall << " tile " << E_tile_expall << " HEC " << E_HEC_expall);
 
        // E resolution of calorimeters
        // Resolution of Hadronic at low energies is not 100% but more like 50%/sqrt(E)
 
        double scaleError_tile = 0.20 + 0.80 * E_tile / (10000. + E_tile);
        double scaleError_HEC = 0.20 + 0.80 * E_HEC / (10000. + E_HEC);
        double sigma_em = std::sqrt(500. * E_em);
        double sigma_tile = scaleError_tile * std::sqrt(1000. * E_tile);
        double sigma_HEC = scaleError_HEC * std::sqrt(1000. * E_HEC);
 
        // go from e.m. scale to Muon Energy scale
 
        E_em *= m_emipEM;
        E_tile *= m_emipTile;
        E_HEC *= m_emipHEC;
        ATH_MSG_DEBUG(" e.m. to Muon scale measured Energies em " << E_em << " tile " << E_tile << " HEC " << E_HEC);
 
        // also for errors
 
        sigma_em *= m_emipEM;
        sigma_tile *= m_emipTile;
        sigma_HEC *= m_emipHEC;
 
        // Average Noise per layer
 
        if (nlay_em > 0) sigma_Noise_em /= nlay_em;
        if (nlay_tile > 0) sigma_Noise_tile /= nlay_tile;
        if (nlay_HEC > 0) sigma_Noise_HEC /= nlay_HEC;
        ATH_MSG_DEBUG(" Average Noise em " << sigma_Noise_em << " nlay " << nlay_em << "  tile " << sigma_Noise_tile << " nlay "
                                           << nlay_tile << "  HEC " << sigma_Noise_HEC << " nlay " << nlay_HEC);
 
        //  apply energy correction using the expected Eloss for noise cut
 
        E_em_threshold += thresholdCorrection(E_em, E_em_expected, m_sigmasAboveNoise * sigma_Noise_em / m_emipEM / 4);
        E_tile_threshold += thresholdCorrection(E_tile, E_tile_expected, m_sigmasAboveNoise * sigma_Noise_tile / m_emipTile / 4);
        E_HEC_threshold += thresholdCorrection(E_HEC, E_HEC_expected, m_sigmasAboveNoise * sigma_Noise_HEC / m_emipHEC / 8);
 
        ATH_MSG_DEBUG(" Threshold correction to Energies em " << E_em_threshold << " tile " << E_tile_threshold << " HEC "
                                                              << E_HEC_threshold);
 
        ATH_MSG_DEBUG(" total Energies em " << E_em + E_em_threshold << " tile " << E_tile + E_tile_threshold << " HEC "
                                            << E_HEC + E_HEC_threshold);
 
        //  Eloss in dead material (so everything not accounted for in the measurement)
 
        //     if(scale_Ionization*Eloss>E_expected) E_expected = scale_Ionization*Eloss;
        E_expected = scale_Ionization * Eloss;
 
        // Corrections for dead material 0.15*E_tile_expected + 0.20*E_HEC_expected;
 
        double E_dead = E_expected - E_em_expected - E_tile_expected - E_HEC_expected + 0.12 * E_tile_expected + 0.27 * E_HEC_expected;
        
 
        //  treatment of FSR
 
        E_FSR = 0.;
        double E_measured = 0.;
        double E_measured_expected = E_em_expected + E_tile_expected + E_HEC_expected;
        if (E_em * std::sin(theta) > m_emEtCut) {
            // large e.m. deposit starting in first e.m. layer
            E_FSR = E_em;
            // do not use measured e.m. energy for muons and use expected (tile and HEC are fine)
            E = E_em_expected + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold + E_dead;
            double sigma_expdead = error_ratio * (E_em_expected + E_dead);
            sigma = std::sqrt(sigma_tile * sigma_tile + sigma_HEC * sigma_HEC + sigma_expdead * sigma_expdead);
            E_measured = E_tile + E_tile_threshold + E_HEC + E_HEC_threshold;
        } else {
            // no FSR
            E = E_em + E_em_threshold + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold + E_dead;
            double sigma_threshold = error_ratio * E_dead;
            sigma = std::sqrt(sigma_em * sigma_em + sigma_tile * sigma_tile + sigma_HEC * sigma_HEC + sigma_threshold * sigma_threshold);
            E_measured = E_em + E_em_threshold + E_tile + E_tile_threshold + E_HEC + E_HEC_threshold;
        }
 
        // eta dependent correction
 
        double etaPos = caloExtension.caloEntryLayerIntersection()->position().eta();
        E += etaCorr(etaPos) * E_expected;
 
        ATH_MSG_DEBUG(" Total energy " << E << " sigma " << sigma << " E calo measured in cells " << E_measured
                                       << " E calo expected in cells " << E_measured_expected << " E_expected meanIoni from TG "
                                       << E_expected);
 
        if (E_em + E_tile + E_HEC < 0.1 * E && E_em + E_tile + E_HEC > 1000.) {
            ATH_MSG_DEBUG(" Real Measured Calorimeter energy " << E_em + E_tile + E_HEC << " is much lower than total " << E
                                                                 << " expected energy too large " << E_expected
                                                                 << " put measured energy to zero ");
            E = 0.;
            sigma = 0.;
        }
        //
        // add for validation (temporarily)
        //
        //   E_em_meas,E_em_exp,E_tile_meas,E_tile_exp,E_HEC_meas,E_HEC_exp,E_dead_exp
        //
 
        //     E_em_meas   = E_em + E_em_threshold;
        //     E_em_exp    = E_em_expected;
        //     E_tile_meas = E_tile + E_tile_threshold;
        //     E_tile_exp  = E_tile_expected;
        //     E_HEC_meas  = E_HEC + E_HEC_threshold;
        //     E_HEC_exp   = E_HEC_expected;
        //     E_dead_exp  = E_dead;
 
        //     double EE_dead = E_expected - EE_emB - EE_emE - EE_tile - EE_HEC + 0.15*EE_tile + 0.20*EE_HEC;
        double EE_dead = E_expected - EE_emB - EE_emE - EE_tile - EE_HEC + 0.12 * EE_tile + 0.27 * EE_HEC;
        //
        //   move expected (to match the TG expected Eloss) and measured energies
        //
        E_em_meas = E_em + E_em_threshold + EE_emB + EE_emE - E_em_expected;
        E_em_exp = EE_emB + EE_emE;
        E_tile_meas = E_tile + E_tile_threshold + EE_tile - E_tile_expected;
        E_tile_exp = EE_tile;
        E_HEC_meas = E_HEC + E_HEC_threshold + EE_HEC - E_HEC_expected;
        E_HEC_exp = EE_HEC;
        E_dead_exp = EE_dead;
 
    }  // calculateMuonEnergies

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

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

etaCorr()

double Rec::MuonCaloEnergyTool::etaCorr ( double eta )

staticprivate

Definition at line 517 of file MuonCaloEnergyTool.cxx.

                                                 {
        // measured energy* = measured energy + etaCorr(eta) * expected
 
        int eta_index = int(std::abs(eta) * (60. / 3.));
        if (eta_index > 59) return 0;
 
       static constexpr double corr[60] = {0.00368146, 0.0419526,  0.0419322,  0.0392922,  0.030304,   0.0262424,  0.0156346,  0.00590235, 0.00772249,
                           -0.0141775, -0.0152247, -0.0174432, -0.0319056, -0.0670813, -0.128678,  -0.0982326, -0.0256406, -0.200244,
                           -0.178975,  -0.120156,  -0.124606,  -0.0961311, -0.0163201, 0.00357829, 0.0292199,  0.0110466,  -0.0375598,
                           0.0417912,  0.0386369,  -0.0689454, -0.21496,   -0.126157,  -0.210573,  -0.215757,  -0.202019,  -0.164546,
                           -0.168607,  -0.128602,  -0.124629,  -0.0882631, -0.100869,  -0.0460344, -0.0709039, -0.0163041, -0.0521138,
                           -0.0125259, -0.0378681, 0.00396062, -0.0636308, -0.032199,  -0.0588335, -0.0470752, -0.0450315, -0.0301302,
                           -0.087378,  -0.0115615, -0.0152037, 0,          0,          0};  // corrections
 
        // additional correction release 21
 
       static constexpr double cor21[20] = {0.999793, 1.00017, 0.990946, 0.995358, 1.01377, 1.02676, 1.03111, 1.01483, 0.995585, 1.00465,
                            1.05224,  1.05238, 1.03208,  1.02373,  1.02305, 1.03975, 1.06547, 1.0364,  1.0361,   1.0361};
 
        int i21 = std::abs(eta * 20. / 3.);
        if (i21 > 19) i21 = 19;
 
        return (cor21[i21] * corr[eta_index]);
    }  // etaCorr

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.

{ 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

initialize()

StatusCode Rec::MuonCaloEnergyTool::initialize ( )

overridevirtual

Definition at line 29 of file MuonCaloEnergyTool.cxx.

                                              {
        // RETRIEVE TOOLS
        ATH_CHECK(m_caloExtensionTool.retrieve());
        ATH_CHECK(m_caloCellAssociationTool.retrieve());
        ATH_CHECK(m_particleCreator.retrieve());
 
        ATH_CHECK(m_caloNoiseCDOKey.initialize());
        ATH_CHECK(m_indetTrackParticleLocation.initialize());
        ATH_CHECK(m_muonTrackParticleLocation.initialize());
 
        return StatusCode::SUCCESS;
    }

inputHandles()

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

overridevirtualinherited

Return this algorithm's input handles.

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

interfaceID()

const InterfaceID & Rec::IMuonCaloEnergyTool::interfaceID ( )

inlinestaticinherited

Definition at line 41 of file IMuonCaloEnergyTool.h.

                                                             {
    return IID_IMuonCaloEnergyTool; 
  }

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.

thresholdCorrection()

double Rec::MuonCaloEnergyTool::thresholdCorrection ( double E_observed, double E_expected, double sigma_Noise ) const

private

Definition at line 542 of file MuonCaloEnergyTool.cxx.

                                                                                                                 {
        //
        // Total energy observed and expected in cells of LAr, Tile, or HEC after 4*sigma_Noise cut
        //                       for a resolution of dE/E of 20%
        //
        // returns a correction to the energy based on the expected energy
        //
        //   one should so use: Etotal = E_observed + E_dead + ThresholdCorrection(E_observed,E_expected,sigma_Noise)
 
        static constexpr double par[6] = {1.01431e+00, -2.26317e+01, 1.80901e+02, -2.33105e+01, 1.82930e+02, 1.11356e-02};
 
        double E = E_expected;
 
        if (E == 0) return 0.;
 
        double x = sigma_Noise / E;
        if (x > 1.) x = 1.;
 
        // formula is for a 4 sigma Noise cut and obtained by a small simulation programm
 
        double fraction = (par[0] + par[1] * x + par[2] * x * x) / (1. + par[3] * x + par[4] * x * x) + par[5] * x;
 
        // for low x values (low thresholds) the fraction is bigger than 1
        // (1.33) because the observed energy overestimates
 
        ATH_MSG_DEBUG(" ThresholdCorrection E " << E << " E_observed not used " << E_observed << " sigma_Noise " << sigma_Noise
                                                << " x = sigma_Noise/E  " << x << " fraction " << fraction << " E correction "
                                                << (1 - fraction) * E);
 
        return (1. - fraction) * E;
    }  // thresholdCorrection

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_caloCellAssociationTool

ToolHandle<Rec::IParticleCaloCellAssociationTool> Rec::MuonCaloEnergyTool::m_caloCellAssociationTool {this, "ParticleCaloCellAssociationTool", "", "Tool to make the cell association"}

private

Definition at line 52 of file MuonCaloEnergyTool.h.

{this, "ParticleCaloCellAssociationTool", "", "Tool to make the cell association"}; 

m_caloExtensionTool

ToolHandle<Trk::IParticleCaloExtensionTool> Rec::MuonCaloEnergyTool::m_caloExtensionTool {this, "ParticleCaloExtensionTool", "", "Tool to make the step-wise extrapolation"}

private

Definition at line 51 of file MuonCaloEnergyTool.h.

{this, "ParticleCaloExtensionTool", "", "Tool to make the step-wise extrapolation"};

m_caloNoiseCDOKey

SG::ReadCondHandleKey<CaloNoise> Rec::MuonCaloEnergyTool::m_caloNoiseCDOKey {this, "CaloNoiseKey", "totalNoise", "SG Key of CaloNoise data object"}

private

Definition at line 55 of file MuonCaloEnergyTool.h.

{this, "CaloNoiseKey", "totalNoise", "SG Key of CaloNoise data object"};

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_emEtCut

Gaudi::Property<double> Rec::MuonCaloEnergyTool::m_emEtCut {this, "EmEtCut", 2.5* Gaudi::Units::GeV }

private

Definition at line 59 of file MuonCaloEnergyTool.h.

{this, "EmEtCut", 2.5* Gaudi::Units::GeV };

m_emF1Cut

Gaudi::Property<double> Rec::MuonCaloEnergyTool::m_emF1Cut {this, "EmF1Cut", 0.15 }

private

Definition at line 60 of file MuonCaloEnergyTool.h.

{this, "EmF1Cut", 0.15 }; 

m_emipEM

double Rec::MuonCaloEnergyTool::m_emipEM {0.42}

staticconstexprprivate

Definition at line 62 of file MuonCaloEnergyTool.h.

{0.42};      

m_emipHEC

double Rec::MuonCaloEnergyTool::m_emipHEC {0.65}

staticconstexprprivate

Definition at line 64 of file MuonCaloEnergyTool.h.

{0.65};     

m_emipTile

double Rec::MuonCaloEnergyTool::m_emipTile {0.86}

staticconstexprprivate

Definition at line 63 of file MuonCaloEnergyTool.h.

{0.86};    

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_indetTrackParticleLocation

SG::ReadHandleKey<xAOD::TrackParticleContainer> Rec::MuonCaloEnergyTool::m_indetTrackParticleLocation

private

Initial value:

{this, "TrackParticleLocation", "InDetTrackParticles",
                                                                                     "ID track particles"}

Definition at line 66 of file MuonCaloEnergyTool.h.

                                                                                {this, "TrackParticleLocation", "InDetTrackParticles",
                                                                                     "ID track particles"};

m_muonTrackParticleLocation

SG::ReadHandleKey<xAOD::TrackParticleContainer> Rec::MuonCaloEnergyTool::m_muonTrackParticleLocation

private

Initial value:

{this, "MuonSpectrometerTrackPArticleLocation",
                                                                                    "MuonSpectrometerTrackParticles", "MS track particles"}

Definition at line 68 of file MuonCaloEnergyTool.h.

                                                                               {this, "MuonSpectrometerTrackPArticleLocation",
                                                                                    "MuonSpectrometerTrackParticles", "MS track particles"};

m_particleCreator

ToolHandle<Trk::ITrackParticleCreatorTool> Rec::MuonCaloEnergyTool::m_particleCreator {this,"TrackParticleCreator", "", "The CB Particle Creator Tool"}

private

Definition at line 53 of file MuonCaloEnergyTool.h.

{this,"TrackParticleCreator",  "", "The CB Particle Creator Tool"};

m_sigmasAboveNoise

Gaudi::Property<double> Rec::MuonCaloEnergyTool::m_sigmasAboveNoise {this, "SigmasAboveNoise", 4.0 }

private

Definition at line 58 of file MuonCaloEnergyTool.h.

{this, "SigmasAboveNoise", 4.0 }; 

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:

• MuonCaloEnergyTool.h
• MuonCaloEnergyTool.cxx