Skip to content
Snippets Groups Projects
Commit e7730896 authored by Eric Torrence's avatar Eric Torrence
Browse files

Add calorimeter reconstruction code

Code framework added to take waveforms and reconstruct total calorimeter.
Placeholder for the future, not ready for primetime just yet.

Also some minor fixes to the waveform hit reco when there isn't a
real signal present.
parent 50735cc2
Branches
No related tags found
No related merge requests found
Pipeline #3227756 passed
Showing
with 576 additions and 3 deletions
################################################################################
# Package: CaloRecAlgs
################################################################################
# Declare the package name:
atlas_subdir( CaloRecAlgs )
# Component(s) in the package:
atlas_add_component( CaloRecAlgs
src/*.cxx src/*.h
src/components/*.cxx
LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserCalorimeter xAODFaserWaveform)
atlas_install_python_modules( python/*.py )
""" Define methods used to instantiate configured Waveform reconstruction tools and algorithms
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
"""
from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
from AthenaConfiguration.ComponentFactory import CompFactory
from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg
#CalorimeterReconstructionTool = CompFactory.CalorimeterReconstructionTool
# One stop shopping for normal FASER data
def CalorimeterReconstructionCfg(flags, **kwargs):
""" Return all algorithms and tools for Waveform reconstruction """
acc = ComponentAccumulator()
#tool = CalorimeterReconstructionTool(name="CaloRecTool", **kwargs)
kwargs.setdefault("CaloWaveHitContainerKey", "CaloWaveformHits")
kwargs.setdefault("PreshowerWaveHitContainerKey", "PreshowerWaveformHits")
kwargs.setdefault("CaloHitContainerKey", "CaloHits")
#kwargs.setdefault("CalorimeterReconstructionTool", tool)
recoAlg = CompFactory.CaloRecAlg("CaloRecAlg", **kwargs)
#recoAlg.CalorimeterReconstructionTool = tool
acc.addEventAlgo(recoAlg)
return acc
def CalorimeterReconstructionOutputCfg(flags, **kwargs):
""" Return ComponentAccumulator with output for Calorimeter Reco"""
acc = ComponentAccumulator()
ItemList = [
"xAOD::CalorimeterHitContainer#*"
, "xAOD::CalorimeterHitAuxContainer#*"
]
acc.merge(OutputStreamCfg(flags, "xAOD", ItemList))
# ostream = acc.getEventAlgo("OutputStreamRDO")
# ostream.TakeItemsFromInput = True # Don't know what this does
return acc
#include "CaloRecAlg.h"
CaloRecAlg::CaloRecAlg(const std::string& name,
ISvcLocator* pSvcLocator)
: AthReentrantAlgorithm(name, pSvcLocator) {
}
StatusCode
CaloRecAlg::initialize() {
ATH_MSG_INFO(name() << "::initalize()" );
// Initalize tools
//ATH_CHECK( m_recoTool.retrieve() );
// Set key to read calo hits from
ATH_CHECK( m_caloWaveHitContainerKey.initialize() );
// Set key to read preshower hits from
ATH_CHECK( m_preshowerWaveHitContainerKey.initialize() );
// Set key to write container
ATH_CHECK( m_caloHitContainerKey.initialize() );
return StatusCode::SUCCESS;
}
StatusCode
CaloRecAlg::finalize() {
ATH_MSG_INFO(name() << "::finalize()");
return StatusCode::SUCCESS;
}
StatusCode
CaloRecAlg::execute(const EventContext& ctx) const {
ATH_MSG_DEBUG("Executing");
ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number()
<< " Event: " << ctx.eventID().event_number());
// Find the input waveform hit containers
SG::ReadHandle<xAOD::WaveformHitContainer> caloWaveHitHandle(m_caloWaveHitContainerKey, ctx);
ATH_CHECK( caloWaveHitHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for WaveformHitContainer " << m_caloWaveHitContainerKey);
if (caloWaveHitHandle->size() == 0) {
ATH_MSG_DEBUG("Calorimeter Waveform Hit container found with zero length!");
}
SG::ReadHandle<xAOD::WaveformHitContainer> preshowerWaveHitHandle(m_preshowerWaveHitContainerKey, ctx);
ATH_CHECK( preshowerWaveHitHandle.isValid() );
ATH_MSG_DEBUG("Found ReadHandle for WaveformHitContainer " << m_preshowerWaveHitContainerKey);
if (preshowerWaveHitHandle->size() == 0) {
ATH_MSG_DEBUG("Preshower Waveform Hit container found with zero length!");
}
// Find the output waveform container
SG::WriteHandle<xAOD::CalorimeterHitContainer> caloHitContainerHandle(m_caloHitContainerKey, ctx);
ATH_CHECK( caloHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(),
std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) );
ATH_MSG_DEBUG("WaveformsHitContainer '" << caloHitContainerHandle.name() << "' initialized");
// Reconstruct all waveforms
//CHECK( m_recoTool->reconstructAll(*waveformHandle, clockptr, hitContainerHandle.ptr()) );
// Find peak time (most significant hit)
const xAOD::WaveformHit* peakHit = findPeakHit(*caloWaveHitHandle);
if (peakHit == NULL) return StatusCode::SUCCESS;
// Create a new calo hit
xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit();
caloHitContainerHandle->push_back(calo_hit);
calo_hit->set_raw_energy(-1.); // Dummy value
// Find closest hits in time per channel
std::map<int, const xAOD::WaveformHit*> hitMap;
for ( const auto& hit : *caloWaveHitHandle ) {
int channel = hit->channel();
if (hitMap.count(channel) == 0)
hitMap[channel] = hit;
else {
if (abs(hitMap[channel]->localtime() - peakHit->localtime()) >
abs(hit->localtime() - peakHit->localtime()))
hitMap[channel] = hit;
}
}
// For each hit found, insert these into the caloHit
// Clear before association
calo_hit->clearCaloWaveformLinks();
for ( const auto& [chan, hit] : hitMap ) {
ATH_MSG_VERBOSE("Found hit " << *hit);
calo_hit->addCaloHit(caloWaveHitHandle.get(), hit);
}
ATH_MSG_DEBUG("CaloHitContainer '" << caloHitContainerHandle.name() << "' filled with "<< caloHitContainerHandle->size() <<" items");
return StatusCode::SUCCESS;
}
const xAOD::WaveformHit*
CaloRecAlg::findPeakHit(const xAOD::WaveformHitContainer& hitContainer) const {
const xAOD::WaveformHit* peakHit = NULL;
for( const auto& hit : hitContainer ) {
if (peakHit == NULL) {
peakHit = hit;
} else {
if ( hit->peak() > peakHit->peak() ) peakHit = hit;
}
}
// Didn't find anything?
if (peakHit == NULL) return NULL;
if (peakHit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) return NULL;
return peakHit;
}
#ifndef CALORECALGS_CALORECALG_H
#define CALORECALGS_CALORECALG_H
// Base class
#include "AthenaBaseComps/AthReentrantAlgorithm.h"
// Data classes
#include "xAODFaserWaveform/WaveformHitContainer.h"
#include "xAODFaserCalorimeter/CalorimeterHit.h"
#include "xAODFaserCalorimeter/CalorimeterHitContainer.h"
#include "xAODFaserCalorimeter/CalorimeterHitAuxContainer.h"
// Tool classes
//#include "CaloRecTools/ICalorimeterReconstructionTool.h"
// Handles
#include "StoreGate/ReadHandleKey.h"
#include "StoreGate/WriteHandleKey.h"
// Gaudi
#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
// STL
#include <string>
class CaloRecAlg : public AthReentrantAlgorithm {
public:
// Constructor
CaloRecAlg(const std::string& name, ISvcLocator* pSvcLocator);
virtual ~CaloRecAlg() = default;
/** @name Usual algorithm methods */
//@{
virtual StatusCode initialize() override;
virtual StatusCode execute(const EventContext& ctx) const override;
virtual StatusCode finalize() override;
//@}
private:
/** @name Disallow default instantiation, copy, assignment */
//@{
CaloRecAlg() = delete;
CaloRecAlg(const CaloRecAlg&) = delete;
CaloRecAlg &operator=(const CaloRecAlg&) = delete;
//@}
/**
* @name Reconstruction tool
*/
//ToolHandle<ICalorimeterReconstructionTool> m_recoTool
//{this, "CalorimeterReconstructionTool", "CalorimeterReconstructionTool"};
/**
* @name Input raw waveform data using SG::ReadHandleKey
*/
//@{
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_caloWaveHitContainerKey
{this, "CaloWaveHitContainerKey", ""};
//@}
//@{
SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerWaveHitContainerKey
{this, "PreshowerWaveHitContainerKey", ""};
//@}
/**
* @name Output data using SG::WriteHandleKey
*/
//@{
SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_caloHitContainerKey
{this, "CaloHitContainerKey", ""};
//@}
const xAOD::WaveformHit*
findPeakHit(const xAOD::WaveformHitContainer& hitContainer) const;
};
#endif // CALORECALGS_CALORECALG_H
#include "../CaloRecAlg.h"
DECLARE_COMPONENT( CaloRecAlg )
......@@ -123,6 +123,10 @@ acc.merge(FaserGeometryCfg(ConfigFlags))
from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionCfg
acc.merge(WaveformReconstructionCfg(ConfigFlags))
# Not ready for primetime
# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg
# acc.merge(CalorimeterReconstructionCfg(ConfigFlags))
# Tracker clusters
from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg
acc.merge(FaserSCT_ClusterizationCfg(ConfigFlags))
......@@ -154,6 +158,10 @@ acc.merge(OutputStreamCfg(ConfigFlags, "xAOD", itemList))
from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg
acc.merge(WaveformReconstructionOutputCfg(ConfigFlags))
# Calorimeter reconstruction output
# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg
# acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags))
# Check what we have
print( "Writing out xAOD objects:" )
print( acc.getEventAlgo("OutputStreamxAOD").ItemList )
......
......
......@@ -658,7 +658,12 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons
rfit.integral = 2*tot; // Factor of 2 because at 500 MHz, dt = 2 ns
rfit.time = rfit.mean;
ATH_MSG_DEBUG( "Initial Mean: " << rfit.mean << " RMS: " << rfit.sigma << " Peak: " << rfit.peak << " Integral: " << rfit.integral);
ATH_MSG_DEBUG(
"Initial Mean: " << rfit.mean
<< " RMS: " << rfit.sigma
<< " Peak: " << rfit.peak
<< " Integral: " << rfit.integral
);
// Fix bad values
if (isnan(rfit.sigma)) {
......@@ -666,6 +671,12 @@ WaveformReconstructionTool::findRawHitValues(const std::vector<float> time, cons
rfit.sigma = 2.;
}
if (rfit.mean < time.front() || rfit.mean > time.back()) {
rfit.mean = (3*time.front() + time.back())/4.; // Set time to 1/4 of way through window
ATH_MSG_DEBUG("Found mean: " << rfit.time << " out of range " << time.front() << "-" << time.back() << " Replace with " << rfit.mean);
rfit.time = rfit.mean;
}
return rfit;
}
......
......
......@@ -84,8 +84,8 @@ class WaveformReconstructionTool: public extends<AthAlgTool, IWaveformReconstruc
//
// Window to define fitting range, in samples (2ns/sample)
IntegerProperty m_windowStart{this, "FitWindowStart", -10};
IntegerProperty m_windowWidth{this, "FitWindowWidth", 50};
IntegerProperty m_windowStart{this, "FitWindowStart", -15};
IntegerProperty m_windowWidth{this, "FitWindowWidth", 60};
//
// Remove overflow values from CB fit
......
......
# Copyright (C) 2020 CERN for the benefit of the FASER collaboration
# Declare the package name.
atlas_subdir( xAODFaserCalorimeter )
# External dependencies.
find_package( xAODUtilities )
# Component(s) in the package.
atlas_add_library( xAODFaserCalorimeter
xAODFaserCalorimeter/*.h xAODFaserCalorimeter/versions/*.h Root/*.cxx
PUBLIC_HEADERS xAODFaserCalorimeter
LINK_LIBRARIES Identifier xAODCore xAODFaserWaveform )
atlas_add_xaod_smart_pointer_dicts(
INPUT xAODFaserCalorimeter/selection.xml
OUTPUT _selectionFile
CONTAINERS "xAOD::CalorimeterHitContainer_v1")
atlas_add_dictionary( xAODFaserCalorimeterDict
xAODFaserCalorimeter/xAODFaserCalorimeterDict.h
${_selectionFile}
LINK_LIBRARIES Identifier xAODCore xAODFaserCalorimeter
EXTRA_FILES Root/dict/*.cxx )
ElementLinks added following DiTauJet example here:
https://gitlab.cern.ch/atlas/athena/-/tree/master/Event/xAOD/xAODTau
\ No newline at end of file
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
// Local include(s):
#include "xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h"
namespace xAOD {
CalorimeterHitAuxContainer_v1::CalorimeterHitAuxContainer_v1()
: AuxContainerBase() {
AUX_VARIABLE(localtime);
AUX_VARIABLE(bcid_time);
AUX_VARIABLE(raw_energy);
AUX_VARIABLE(caloLinks);
}
} // namespace xAOD
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
// EDM include(s):
#include "xAODCore/AuxStoreAccessorMacros.h"
// Local include(s):
#include "xAODFaserCalorimeter/versions/CalorimeterHit_v1.h"
namespace xAOD {
CalorimeterHit_v1::CalorimeterHit_v1() : SG::AuxElement() {
}
AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, localtime, set_localtime )
AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, bcidtime, set_bcidtime )
AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, raw_energy, set_raw_energy )
// setters and getters for the Calo WaveformHit links
AUXSTORE_OBJECT_SETTER_AND_GETTER( CalorimeterHit_v1,
CalorimeterHit_v1::WaveformHitLinks_t,
caloWaveformLinks,
setCaloWaveformLinks )
static const SG::AuxElement::Accessor< CalorimeterHit_v1::WaveformHitLinks_t > caloHitAcc( "caloWaveformLinks" );
const WaveformHit* CalorimeterHit_v1::caloHit( size_t i ) const {
return ( *caloHitAcc( *this )[ i ] );
}
size_t CalorimeterHit_v1::nCaloHits() const {
return caloHitAcc( *this ).size();
}
void CalorimeterHit_v1::addCaloHit( const xAOD::WaveformHitContainer* pWaveformHitContainer,
const xAOD::WaveformHit* pWaveformHit) {
ElementLink< xAOD::WaveformHitContainer > linkToWaveformHit;
linkToWaveformHit.toContainedElement(*pWaveformHitContainer, pWaveformHit);
caloHitAcc( *this ).push_back( linkToWaveformHit );
return;
}
void CalorimeterHit_v1::clearCaloWaveformLinks() {
caloHitAcc( *this ).clear();
return;
}
} // namespace xAOD
namespace xAOD {
std::ostream& operator<<(std::ostream& s, const xAOD::CalorimeterHit_v1& hit) {
s << "xAODCalorimeterHit:"
<< " local time=" << hit.localtime()
<< " raw_energy=" << hit.raw_energy()
<< std::endl;
return s;
}
} // namespace xAOD
// EDM include(s):
#include "xAODCore/AddDVProxy.h"
// Local include(s):
#include "xAODFaserCalorimeter/CalorimeterHitContainer.h"
// Set up the collection proxies:
ADD_NS_DV_PROXY( xAOD, CalorimeterHitContainer );
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
//simple includes to force the CLASS_DEF etc to be encountered during compile
#include "xAODFaserCalorimeter/CalorimeterHitContainer.h"
#include "xAODFaserCalorimeter/CalorimeterHitAuxContainer.h"
// Dear emacs, this is -*- c++ -*-
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
#ifndef XAODFASERCALORIMETER_CALORIMETERHIT_H
#define XAODFASERCALORIMETER_CALORIMETERHIT_H
// Local include(s):
#include "xAODFaserCalorimeter/versions/CalorimeterHit_v1.h"
namespace xAOD {
/// Declare the latest version of the class
typedef CalorimeterHit_v1 CalorimeterHit;
}
// Set up a CLID for the container:
#include "xAODCore/CLASS_DEF.h"
CLASS_DEF( xAOD::CalorimeterHit, 195445226, 1 )
#endif // XAODFASERCALORIMETER_CALORIMETERHIT_H
// Dear emacs, this is -*- c++ -*-
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
#ifndef XAODFASERCALORIMETER_CALORIMETERHITAUXCONTAINER_H
#define XAODFASERCALORIMETER_CALORIMETERHITAUXCONTAINER_H
// Local include(s):
#include "xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h"
namespace xAOD {
/// Declare the latest version of the class
typedef CalorimeterHitAuxContainer_v1 CalorimeterHitAuxContainer;
}
// Set up a CLID for the container:
#include "xAODCore/CLASS_DEF.h"
CLASS_DEF( xAOD::CalorimeterHitAuxContainer, 1074912337, 1 )
#endif // XAODFASERCALORIMETER_CALORIMETERHITAUXCONTAINER_H
// Dear emacs, this is -*- c++ -*-
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
#ifndef XAODFASERCALORIMETER_CALORIMETERHITCONTAINER_H
#define XAODFASERCALORIMETER_CALORIMETERHITCONTAINER_H
// Local include(s):
#include "xAODFaserCalorimeter/versions/CalorimeterHitContainer_v1.h"
namespace xAOD {
/// Declare the latest version of the class
typedef CalorimeterHitContainer_v1 CalorimeterHitContainer;
}
// Set up a CLID for the container:
#include "xAODCore/CLASS_DEF.h"
CLASS_DEF( xAOD::CalorimeterHitContainer, 1147607550, 1 )
#endif // XAODFASERCALORIMETER_CALORIMETERHITCONTAINER_H
<!-- Copyright (C) 2020 CERN for the benefit of the FASER collaboration -->
<lcgdict>
<class name="xAOD::CalorimeterHit_v1" />
<typedef name="xAOD::CalorimeterHit" />
<class name="xAOD::CalorimeterHitContainer_v1"
id="1172f311-f09a-4e3b-a946-68991c2c0646" />
<typedef name="xAOD::CalorimeterHitContainer" />
<class name="xAOD::CalorimeterHitAuxContainer_v1"
id="755f71cf-2148-4736-b26e-3c007d627861" />
<typedef name="xAOD::CalorimeterHitAuxContainer" />
</lcgdict>
// Dear emacs, this is -*- c++ -*-
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
#ifndef XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITAUXCONTAINER_V1_H
#define XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITAUXCONTAINER_V1_H
// STL include(s):
#include <vector>
// EDM include(s):
#include "xAODCore/AuxContainerBase.h"
#include "xAODFaserWaveform/WaveformHitContainer.h"
namespace xAOD {
/// Auxiliary container for CalorimeterHit containers
class CalorimeterHitAuxContainer_v1 : public AuxContainerBase {
public:
/// Default constructor
CalorimeterHitAuxContainer_v1();
/// Destructor
~CalorimeterHitAuxContainer_v1() {}
private:
/// @name Basic variables
///@ {
std::vector<float> localtime;
std::vector<float> bcid_time;
std::vector<float> raw_energy;
typedef std::vector< ElementLink< WaveformHitContainer > > WaveformHitLink_t;
std::vector< WaveformHitLink_t > caloLinks;
std::vector< WaveformHitLink_t > preshowerLinks;
///@}
}; // class CalorimeterHitAuxContainer_v1
} // namespace xAOD
// Set up a CLID and StoreGate inheritance for the class:
#include "xAODCore/BaseInfo.h"
SG_BASE( xAOD::CalorimeterHitAuxContainer_v1, xAOD::AuxContainerBase );
#endif // XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITAUXCONTAINER_V1_H
// Dear emacs, this is -*- c++ -*-
/*
Copyright (C) 2020 CERN for the benefit of the FASER collaboration
*/
#ifndef XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITCONTAINER_V1_H
#define XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITCONTAINER_V1_H
// System include(s):
extern "C" {
# include "stdint.h"
}
// EDM include(s):
#include "AthContainers/DataVector.h"
// Local includes:
#include "xAODFaserCalorimeter/versions/CalorimeterHit_v1.h"
namespace xAOD {
// Define the container as a simple DataVector
typedef DataVector<CalorimeterHit_v1> CalorimeterHitContainer_v1;
}
#endif // XAODFASERCALORIMETER_VERSIONS_CALORIMETERHITCONTAINER_V1_H
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment