diff --git a/Calorimeter/CaloRecAlgs/CMakeLists.txt b/Calorimeter/CaloRecAlgs/CMakeLists.txt index 219db92a3a8927ae9cc025bea46985387e072ce0..75d9d09a23eb857555e55848ee16b1d2b98b89f3 100644 --- a/Calorimeter/CaloRecAlgs/CMakeLists.txt +++ b/Calorimeter/CaloRecAlgs/CMakeLists.txt @@ -9,7 +9,7 @@ atlas_subdir( CaloRecAlgs ) atlas_add_component( CaloRecAlgs src/*.cxx src/*.h src/components/*.cxx - LINK_LIBRARIES AthenaBaseComps StoreGateLib xAODFaserCalorimeter xAODFaserWaveform) + LINK_LIBRARIES AthenaKernel GaudiKernel AthenaBaseComps AthenaPoolUtilities StoreGateLib xAODFaserCalorimeter xAODFaserWaveform CaloRecToolsLib) atlas_install_python_modules( python/*.py ) diff --git a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py index 034bf08b5b36a4620567f6716eca5d31338a1c04..83839b9cba1a4c6a50eef52852ed8ee2127dfaa7 100644 --- a/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py +++ b/Calorimeter/CaloRecAlgs/python/CaloRecAlgsConfig.py @@ -1,30 +1,31 @@ -""" Define methods used to instantiate configured Waveform reconstruction tools and algorithms +""" Define methods used to instantiate configured Calorimeter Calibration reconstruction tools and algorithms -Copyright (C) 2020 CERN for the benefit of the FASER collaboration +Copyright (C) 2022 CERN for the benefit of the FASER collaboration """ from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator from AthenaConfiguration.ComponentFactory import CompFactory - +from IOVDbSvc.IOVDbSvcConfig import addFolders from OutputStreamAthenaPool.OutputStreamConfig import OutputStreamCfg -#CalorimeterReconstructionTool = CompFactory.CalorimeterReconstructionTool +from CaloRecTools.CaloRecToolConfig import CaloRecToolCfg # 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) + kwargs.setdefault("PreshowerHitContainerKey", "PreshowerHits") recoAlg = CompFactory.CaloRecAlg("CaloRecAlg", **kwargs) - #recoAlg.CalorimeterReconstructionTool = tool acc.addEventAlgo(recoAlg) + dbInstance = kwargs.get("dbInstance", "TRIGGER_OFL") # what should this be set to??? + + acc.merge(CaloRecToolCfg(flags)) + return acc def CalorimeterReconstructionOutputCfg(flags, **kwargs): diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx index 4b4d6f158f29478a60443381c5292d254d1c665f..12ddca3a44a5ba4691ff3feffa2adae37d19e2d5 100644 --- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx +++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.cxx @@ -1,4 +1,5 @@ #include "CaloRecAlg.h" +#include <math.h> CaloRecAlg::CaloRecAlg(const std::string& name, ISvcLocator* pSvcLocator) @@ -6,13 +7,9 @@ CaloRecAlg::CaloRecAlg(const std::string& name, } -StatusCode -CaloRecAlg::initialize() { +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() ); @@ -21,19 +18,21 @@ CaloRecAlg::initialize() { // Set key to write container ATH_CHECK( m_caloHitContainerKey.initialize() ); + ATH_CHECK( m_preshowerHitContainerKey.initialize() ); + + // Initalize tools + ATH_CHECK( m_recoCalibTool.retrieve() ); return StatusCode::SUCCESS; } - -StatusCode -CaloRecAlg::finalize() { +//---------------------------------------------------------------------- +StatusCode CaloRecAlg::finalize() { ATH_MSG_INFO(name() << "::finalize()"); return StatusCode::SUCCESS; } - -StatusCode -CaloRecAlg::execute(const EventContext& ctx) const { +//---------------------------------------------------------------------- +StatusCode CaloRecAlg::execute(const EventContext& ctx) const { ATH_MSG_DEBUG("Executing"); ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() @@ -63,61 +62,115 @@ CaloRecAlg::execute(const EventContext& ctx) const { ATH_CHECK( caloHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(), std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) ); + SG::WriteHandle<xAOD::CalorimeterHitContainer> preshowerHitContainerHandle(m_preshowerHitContainerKey, ctx); + ATH_CHECK( preshowerHitContainerHandle.record( std::make_unique<xAOD::CalorimeterHitContainer>(), + std::make_unique<xAOD::CalorimeterHitAuxContainer>() ) ); + ATH_MSG_DEBUG("WaveformsHitContainer '" << caloHitContainerHandle.name() << "' initialized"); + ATH_MSG_DEBUG("WaveformsHitContainer '" << preshowerHitContainerHandle.name() << "' initialized"); + + // Loop over calo hits and calibrate each primary hit + for( const auto& hit : *caloWaveHitHandle ) { + if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue; + + // Create a new calo hit + xAOD::CalorimeterHit* calo_hit = new xAOD::CalorimeterHit(); + caloHitContainerHandle->push_back(calo_hit); + + ATH_MSG_DEBUG("calo_hit in channel " << hit->channel() ); - // 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; + float MIPcharge_ref = m_recoCalibTool->getMIPcharge_ref(hit->channel()); // get reference MIP charge from database + + float hit_Edep = 0; + if (m_isMC) { + hit_Edep = hit->integral() / MIPcharge_ref; // MC MIPcharge_ref from database is really digitization scaling factors that take simulated Edep in MeV to mV*ns + ATH_MSG_DEBUG("isMC == True: calo hit Edep = " << hit_Edep << " MeV"); + } else { + float hit_charge = hit->integral()/50.0; // divide by 50 ohms to get charge + ATH_MSG_DEBUG("calo_hit filled has charge of " << hit_charge << " pC"); + + float HVgain_extrap = extrapolateHVgain(hit->channel()); + ATH_MSG_DEBUG("HV gain = " << HVgain_extrap ); + + float hit_MIP = (hit_charge * HVgain_extrap) / MIPcharge_ref; + ATH_MSG_DEBUG("Edep in MIPs = " << hit_MIP ); + + hit_Edep = hit_MIP * m_MIP_sim_Edep_calo; + ATH_MSG_DEBUG("Edep in MeV = " << hit_Edep ); } - } - // 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); + calo_hit->set_Edep(hit_Edep); // set calibrated Edep value + + float hit_raw_Edep = 0.0; + if (hit->integral() != 0.0) { // avoid possibility of division by zero error + hit_raw_Edep = hit_Edep * (hit->raw_integral() / hit->integral()); + } + calo_hit->set_raw_Edep(hit_raw_Edep); // set calibrated raw_Edep value + + calo_hit->set_channel(hit->channel()); // set channel number + + calo_hit->clearWaveformLinks(); + calo_hit->addHit(caloWaveHitHandle.get(), hit); // create link to calo waveform hit } ATH_MSG_DEBUG("CaloHitContainer '" << caloHitContainerHandle.name() << "' filled with "<< caloHitContainerHandle->size() <<" items"); - return StatusCode::SUCCESS; -} + for( const auto& hit : *preshowerWaveHitHandle ) { + if (hit->status_bit(xAOD::WaveformStatus::SECONDARY)) continue; + + // Create a new preshower hit + xAOD::CalorimeterHit* preshower_hit = new xAOD::CalorimeterHit(); + preshowerHitContainerHandle->push_back(preshower_hit); + + ATH_MSG_DEBUG("preshower_hit in channel " << hit->channel() ); -const xAOD::WaveformHit* -CaloRecAlg::findPeakHit(const xAOD::WaveformHitContainer& hitContainer) const { + float MIPcharge_ref = m_recoCalibTool->getMIPcharge_ref(hit->channel()); // get reference MIP charge from database - const xAOD::WaveformHit* peakHit = NULL; - for( const auto& hit : hitContainer ) { - if (peakHit == NULL) { - peakHit = hit; + float hit_Edep = 0; + if (m_isMC) { + hit_Edep = hit->integral() / MIPcharge_ref; // MC MIPcharge_ref from database is really digitization scaling factors that take simulated Edep in MeV to mV*ns + ATH_MSG_DEBUG("isMC == True: preshower hit Edep = " << hit_Edep << " MeV"); } else { - if ( hit->peak() > peakHit->peak() ) peakHit = hit; + float hit_charge = hit->integral()/50.0; // divide by 50 ohms to get charge + ATH_MSG_DEBUG("preshower_hit filled has charge of " << hit_charge << " pC"); + + float hit_MIP = hit_charge / MIPcharge_ref; + ATH_MSG_DEBUG("Edep in MIPs = " << hit_MIP ); + + hit_Edep = hit_MIP * m_MIP_sim_Edep_preshower; + ATH_MSG_DEBUG("Edep in GeV = " << hit_Edep ); + } + + preshower_hit->set_Edep(hit_Edep); // set calibrated Edep value + + float hit_raw_Edep = 0.0; + if (hit->integral() != 0.0) { // avoid possibility of division by zero error + hit_raw_Edep = hit_Edep * (hit->raw_integral() / hit->integral()); } + preshower_hit->set_raw_Edep(hit_raw_Edep); // set calibrated rsw_Edep value + + preshower_hit->set_channel(hit->channel()); // set channel number + + preshower_hit->clearWaveformLinks(); + preshower_hit->addHit(preshowerWaveHitHandle.get(), hit); // create link to preshower waveform hit } - // Didn't find anything? - if (peakHit == NULL) return NULL; - if (peakHit->status_bit(xAOD::WaveformStatus::THRESHOLD_FAILED)) return NULL; - return peakHit; + ATH_MSG_DEBUG("PreshowerHitContainer '" << preshowerHitContainerHandle.name() << "' filled with "<< preshowerHitContainerHandle->size() <<" items"); + + + return StatusCode::SUCCESS; } + +//---------------------------------------------------------------------- +float CaloRecAlg::extrapolateHVgain(int channel) const { + float PMT_hv = m_recoCalibTool->getHV(channel); + float PMT_hv_ref = m_recoCalibTool->getHV_ref(channel); + TF1 gaincurve = m_recoCalibTool->get_PMT_HV_curve(channel); + + float gaincurve_atHV = gaincurve.Eval(PMT_hv); + float gaincurve_atHVref = gaincurve.Eval(PMT_hv_ref); + + return ( gaincurve_atHVref / gaincurve_atHV ) * pow( PMT_hv_ref / PMT_hv , 6.6); +} +//---------------------------------------------------------------------- + diff --git a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h index bd98e9afd6b693b8cc4518ebd87ddd72207fddef..c9ff4c5146d9e59c94a279345f6bb5977519c997 100644 --- a/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h +++ b/Calorimeter/CaloRecAlgs/src/CaloRecAlg.h @@ -12,7 +12,11 @@ #include "xAODFaserCalorimeter/CalorimeterHitAuxContainer.h" // Tool classes -//#include "CaloRecTools/ICalorimeterReconstructionTool.h" +#include "CaloRecTools/ICaloRecTool.h" + +// Include Athena stuff for Conditions db reading +#include "AthenaPoolUtilities/CondAttrListCollection.h" +#include "StoreGate/ReadCondHandleKey.h" // Handles #include "StoreGate/ReadHandleKey.h" @@ -21,9 +25,16 @@ // Gaudi #include "GaudiKernel/ServiceHandle.h" #include "GaudiKernel/ToolHandle.h" +#include "GaudiKernel/EventContext.h" +#include "GaudiKernel/ICondSvc.h" +#include "Gaudi/Property.h" + +// ROOT +#include "TF1.h" // STL #include <string> +#include <vector> class CaloRecAlg : public AthReentrantAlgorithm { @@ -51,32 +62,33 @@ class CaloRecAlg : public AthReentrantAlgorithm { /** * @name Reconstruction tool */ - //ToolHandle<ICalorimeterReconstructionTool> m_recoTool - //{this, "CalorimeterReconstructionTool", "CalorimeterReconstructionTool"}; + ToolHandle<ICaloRecTool> m_recoCalibTool {this, "CaloRecTool", "CaloRecTool"}; /** * @name Input raw waveform data using SG::ReadHandleKey */ //@{ - SG::ReadHandleKey<xAOD::WaveformHitContainer> m_caloWaveHitContainerKey - {this, "CaloWaveHitContainerKey", ""}; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_caloWaveHitContainerKey {this, "CaloWaveHitContainerKey", "CaloWaveformHits"}; //@} //@{ - SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerWaveHitContainerKey - {this, "PreshowerWaveHitContainerKey", ""}; + SG::ReadHandleKey<xAOD::WaveformHitContainer> m_preshowerWaveHitContainerKey {this, "PreshowerWaveHitContainerKey", "PreshowerWaveformHits"}; //@} /** * @name Output data using SG::WriteHandleKey */ //@{ - SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_caloHitContainerKey - {this, "CaloHitContainerKey", ""}; + SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_caloHitContainerKey {this, "CaloHitContainerKey", "CaloHits"}; + SG::WriteHandleKey<xAOD::CalorimeterHitContainer> m_preshowerHitContainerKey {this, "PreshowerHitContainerKey", "PreshowerHits"}; //@} - const xAOD::WaveformHit* - findPeakHit(const xAOD::WaveformHitContainer& hitContainer) const; + float extrapolateHVgain(int channel) const; + + FloatProperty m_MIP_sim_Edep_calo {this, "MIP_sim_Edep_calo", 58.5}; // MIP deposits 5.85 MeV of energy in calo + FloatProperty m_MIP_sim_Edep_preshower {this, "MIP_sim_Edep_preshower", 4.894}; // MIP deposits 4.894 MeV of energy in a preshower layer + + Gaudi::Property<bool> m_isMC {this, "isMC", false}; }; diff --git a/Calorimeter/CaloRecTools/CMakeLists.txt b/Calorimeter/CaloRecTools/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..f0e9f07b53f5426b360c2639d4f7953a44d57792 --- /dev/null +++ b/Calorimeter/CaloRecTools/CMakeLists.txt @@ -0,0 +1,30 @@ +################################################################################ +# Package: CaloRecTools +################################################################################ + +# Declare the package name: +atlas_subdir( CaloRecTools ) + +# External dependencies: +find_package( ROOT ) + +# Component(s) in the package: +atlas_add_library( CaloRecToolsLib + CaloRecTools/*.h src/*.cxx src/*.h + PUBLIC_HEADERS CaloRecTools + PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES AthenaBaseComps AthenaPoolUtilities AthenaKernel GaudiKernel + PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} + ) + +atlas_add_component( CaloRecTools + src/components/*.cxx + INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} + LINK_LIBRARIES ${ROOT_LIBRARIES} + AthenaBaseComps AthenaPoolUtilities AthenaKernel GaudiKernel + CaloRecToolsLib) + +# Install files from the package: +atlas_install_python_modules( python/*.py ) + + diff --git a/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h b/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h new file mode 100644 index 0000000000000000000000000000000000000000..feb141d915d6214a681f437d94b8118e14c33820 --- /dev/null +++ b/Calorimeter/CaloRecTools/CaloRecTools/ICaloRecTool.h @@ -0,0 +1,45 @@ +/* + Copyright (C) 2022 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file ICaloRecTool.h + * Header file for the ICaloRecTool class + * @author Deion Fellers, 2022 + */ + + +#ifndef CALORECTOOL_ICALORECTOOL_H +#define CALORECTOOL_ICALORECTOOL_H + +// Base class +#include "GaudiKernel/IAlgTool.h" +#include "GaudiKernel/ToolHandle.h" +#include "GaudiKernel/EventContext.h" + +#include "TF1.h" + +///Interface for Calo reco algorithms +class ICaloRecTool : virtual public IAlgTool +{ + public: + + // InterfaceID + DeclareInterfaceID(ICaloRecTool, 1, 0); + + virtual ~ICaloRecTool() = default; //!< Destructor + + // methods to return calibration database data + virtual float getHV(const EventContext& ctx, int channel) const = 0; + virtual float getHV(int channel) const = 0; + + virtual float getHV_ref(const EventContext& ctx, int channel) const = 0; + virtual float getHV_ref(int channel) const = 0; + + virtual float getMIPcharge_ref(const EventContext& ctx, int channel) const = 0; + virtual float getMIPcharge_ref(int channel) const = 0; + + virtual TF1 get_PMT_HV_curve(int channel) const = 0; +}; + +#endif // CALORECTOOL_ICALORECTOOL_H diff --git a/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py b/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py new file mode 100644 index 0000000000000000000000000000000000000000..e87ac3658ec8f92caf00bc19e061d3d9d95c8e5d --- /dev/null +++ b/Calorimeter/CaloRecTools/python/CaloRecToolConfig.py @@ -0,0 +1,28 @@ +""" Define methods to configure CaloRecTool + +Copyright (C) 2022 CERN for the benefit of the FASER collaboration +""" +from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator +from AthenaConfiguration.ComponentFactory import CompFactory +from IOVDbSvc.IOVDbSvcConfig import addFolders +CaloRecTool=CompFactory.CaloRecTool + + +def CaloRecToolCfg(flags, name="CaloRecTool", **kwargs): + """ Return configured ComponentAccumulator and tool for Calo Calibration + + CaloRecTool may be provided in kwargs + """ + + acc = ComponentAccumulator() + # tool = kwargs.get("CaloRecTool", CaloRecTool(flags)) + # Probably need to figure this out! + dbInstance = kwargs.get("dbInstance", "TRIGGER_OFL") + + acc.merge(addFolders(flags, "/WAVE/Calibration/HV", dbInstance, className="CondAttrListCollection")) + acc.merge(addFolders(flags, "/WAVE/Calibration/HV_ref", dbInstance, className="CondAttrListCollection")) + acc.merge(addFolders(flags, "/WAVE/Calibration/MIP_charge_ref", dbInstance, className="CondAttrListCollection")) + + return acc + + diff --git a/Calorimeter/CaloRecTools/src/CaloRecTool.cxx b/Calorimeter/CaloRecTools/src/CaloRecTool.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a6bb42c5140d10c1d2efdc486201c2bab81c0365 --- /dev/null +++ b/Calorimeter/CaloRecTools/src/CaloRecTool.cxx @@ -0,0 +1,195 @@ +/* + Copyright (C) 2022 CERN for the benefit of the FASER collaboration +*/ + +/** + * @file CaloRecTool.cxx + * Implementation file for the CaloRecTool.cxx + * @ author Deion Fellers, 2022 + **/ + +#include "CaloRecTool.h" + +// Constructor +CaloRecTool::CaloRecTool(const std::string& type, const std::string& name, const IInterface* parent) : + base_class(type, name, parent) +{ +} + +//-------------------------------------------------------------- +StatusCode +CaloRecTool::initialize() { + ATH_MSG_INFO( name() << "::initalize()" ); + + // Set keys to read calibratiion variables from data base + ATH_CHECK( m_PMT_HV_ReadKey.initialize() ); + ATH_CHECK( m_PMT_HV_ref_ReadKey.initialize() ); + ATH_CHECK( m_MIPcharge_ref_ReadKey.initialize() ); + + // access and store calo PMT HV gain curves + HVgaincurves_rootFile = TFile::Open(m_PMT_HV_Gain_Curve_file.value().c_str(),"read"); // open root file that has TF1 gain curves stored + + chan0_HVgaincurve_pntr = HVgaincurves_rootFile->Get<TF1>(m_chan0_HVgaincurve_name.value().c_str()); // make pointers to the gain curves + chan1_HVgaincurve_pntr = HVgaincurves_rootFile->Get<TF1>(m_chan1_HVgaincurve_name.value().c_str()); + chan2_HVgaincurve_pntr = HVgaincurves_rootFile->Get<TF1>(m_chan2_HVgaincurve_name.value().c_str()); + chan3_HVgaincurve_pntr = HVgaincurves_rootFile->Get<TF1>(m_chan3_HVgaincurve_name.value().c_str()); + + m_HVgainCurveMap[0] = *chan0_HVgaincurve_pntr; // store TF1 objects in an array mapped to the calo channel numbers + m_HVgainCurveMap[1] = *chan1_HVgaincurve_pntr; + m_HVgainCurveMap[2] = *chan2_HVgaincurve_pntr; + m_HVgainCurveMap[3] = *chan3_HVgaincurve_pntr; + + HVgaincurves_rootFile->Close(); // close the root file + + + return StatusCode::SUCCESS; +} + +//-------------------------------------------------------------- +StatusCode +CaloRecTool::finalize() { + // Print where you are + return StatusCode::SUCCESS; +} + +//-------------------------------------------------------------- +TF1 CaloRecTool::get_PMT_HV_curve(int channel) const { + if (channel <= 4) { + return m_HVgainCurveMap[channel]; + } else { + ATH_MSG_WARNING("channel "<<channel<<" is not <= 4 and thus not a calorimeter channel, so no HV gain-curve exists!"); + } + return TF1("default", "1.0", 0.0, 2000.0); +} + +//-------------------------------------------------------------- +float CaloRecTool::getHV(const EventContext& ctx, int channel) const { + + ATH_MSG_DEBUG("in getHV("<<channel<<")"); + + float HV=0.; + + // Read Cond Handle + SG::ReadCondHandle<CondAttrListCollection> readHandle{m_PMT_HV_ReadKey, ctx}; + const CondAttrListCollection* readCdo{*readHandle}; + if (readCdo==nullptr) { + ATH_MSG_FATAL("Null pointer to the read conditions object"); + return HV; + } + // Get the validitiy range + EventIDRange rangeW; + if (not readHandle.range(rangeW)) { + ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key()); + return HV; + } + ATH_MSG_DEBUG("Size of CondAttrListCollection " << readHandle.fullKey() << " readCdo->size()= " << readCdo->size()); + ATH_MSG_DEBUG("Range of input is " << rangeW); + + // Read offset for specific channel + const CondAttrListCollection::AttributeList& payload{readCdo->attributeList(channel)}; + + if (payload.exists("HV") and not payload["HV"].isNull()) { + HV = payload["HV"].data<float>(); + ATH_MSG_DEBUG("Found digitizer channel " << channel << ", HV as " << HV); + } else { + ATH_MSG_WARNING("No valid HV found for channel "<<channel<<"!"); + } + + return HV; + +} + +float CaloRecTool::getHV(int channel) const { + const EventContext& ctx(Gaudi::Hive::currentContext()); + return CaloRecTool::getHV(ctx, channel); +} + +//---------------------------------------------------------------------- +float CaloRecTool::getHV_ref(const EventContext& ctx, int channel) const { + + ATH_MSG_DEBUG("in getHV_ref("<<channel<<")"); + + float HV_ref=0.; + + // Read Cond Handle + SG::ReadCondHandle<CondAttrListCollection> readHandle{m_PMT_HV_ref_ReadKey, ctx}; + const CondAttrListCollection* readCdo{*readHandle}; + if (readCdo==nullptr) { + ATH_MSG_FATAL("Null pointer to the read conditions object"); + return HV_ref; + } + // Get the validitiy range + EventIDRange rangeW; + if (not readHandle.range(rangeW)) { + ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key()); + return HV_ref; + } + ATH_MSG_DEBUG("Size of CondAttrListCollection " << readHandle.fullKey() << " readCdo->size()= " << readCdo->size()); + ATH_MSG_DEBUG("Range of input is " << rangeW); + + // Read offset for specific channel + const CondAttrListCollection::AttributeList& payload{readCdo->attributeList(channel)}; + + if (payload.exists("HV_ref") and not payload["HV_ref"].isNull()) { + HV_ref = payload["HV_ref"].data<float>(); + ATH_MSG_DEBUG("Found digitizer channel " << channel << ", HV_ref as " << HV_ref); + } else { + ATH_MSG_WARNING("No valid HV_ref found for channel "<<channel<<"!"); + } + + return HV_ref; + +} + +float CaloRecTool::getHV_ref(int channel) const { + const EventContext& ctx(Gaudi::Hive::currentContext()); + return CaloRecTool::getHV_ref(ctx, channel); +} + +//---------------------------------------------------------------------- +float CaloRecTool::getMIPcharge_ref(const EventContext& ctx, int channel) const { + + ATH_MSG_DEBUG("in getMIPcharge_ref("<<channel<<")"); + + float MIP_charge_ref =0.; + + // Read Cond Handle + SG::ReadCondHandle<CondAttrListCollection> readHandle{m_MIPcharge_ref_ReadKey, ctx}; + const CondAttrListCollection* readCdo{*readHandle}; + if (readCdo==nullptr) { + ATH_MSG_FATAL("Null pointer to the read conditions object"); + return MIP_charge_ref; + } + // Get the validitiy range + EventIDRange rangeW; + if (not readHandle.range(rangeW)) { + ATH_MSG_FATAL("Failed to retrieve validity range for " << readHandle.key()); + return MIP_charge_ref; + } + ATH_MSG_DEBUG("Size of CondAttrListCollection " << readHandle.fullKey() << " readCdo->size()= " << readCdo->size()); + ATH_MSG_DEBUG("Range of input is " << rangeW); + + // Read offset for specific channel + const CondAttrListCollection::AttributeList& payload{readCdo->attributeList(channel)}; + + if (payload.exists("MIP_charge_ref") and not payload["MIP_charge_ref"].isNull()) { + MIP_charge_ref = payload["MIP_charge_ref"].data<float>(); + ATH_MSG_DEBUG("Found digitizer channel " << channel << ", MIP_charge_ref as " << MIP_charge_ref); + } else { + ATH_MSG_WARNING("No valid MIP_charge_ref found for channel "<<channel<<"!"); + } + + return MIP_charge_ref; + +} + +float CaloRecTool::getMIPcharge_ref(int channel) const { + const EventContext& ctx(Gaudi::Hive::currentContext()); + return CaloRecTool::getMIPcharge_ref(ctx, channel); +} + +//-------------------------------------------------------------- + + + + diff --git a/Calorimeter/CaloRecTools/src/CaloRecTool.h b/Calorimeter/CaloRecTools/src/CaloRecTool.h new file mode 100644 index 0000000000000000000000000000000000000000..c6dc951818cb59472fe6a1ff077e97e6753ec509 --- /dev/null +++ b/Calorimeter/CaloRecTools/src/CaloRecTool.h @@ -0,0 +1,86 @@ +/* + Copyright (C) 2021 CERN for the benefit of the FASER collaboration +*/ + +/** @file CaloRecTool.h + * Header file for CaloRecTool.h + * + */ +#ifndef CALORECTOOLS_CALORECTOOL_H +#define CALORECTOOLS_CALORECTOOL_H + +// Include interface class +#include "AthenaBaseComps/AthAlgTool.h" +#include "CaloRecTools/ICaloRecTool.h" + +// Include Athena stuff +#include "AthenaPoolUtilities/CondAttrListCollection.h" +#include "StoreGate/ReadCondHandleKey.h" + +#include "GaudiKernel/ICondSvc.h" +#include "Gaudi/Property.h" + +// Include Gaudi classes +#include "GaudiKernel/ICondSvc.h" +#include "Gaudi/Property.h" +#include "GaudiKernel/EventContext.h" + +// Include ROOT classes +#include "TF1.h" +#include "TFile.h" + +#include <string> + +class CaloRecTool: public extends<AthAlgTool, ICaloRecTool> { + public: + + /// Normal constructor for an AlgTool; 'properties' are also declared here + CaloRecTool(const std::string& type, const std::string& name, const IInterface* parent); + + // Standard Gaudi functions + virtual StatusCode initialize() override; //!< Gaudi initialiser + virtual StatusCode finalize() override; //!< Gaudi finaliser + + // methods for returning calibration data + // Channels indexed by digitizer channel number + // HV is in Volts and MIPcharge is in pC + // + virtual float getHV(const EventContext& ctx, int channel) const override; + virtual float getHV(int channel) const override; + + virtual float getHV_ref(const EventContext& ctx, int channel) const override; + virtual float getHV_ref(int channel) const override; + + virtual float getMIPcharge_ref(const EventContext& ctx, int channel) const override; + virtual float getMIPcharge_ref(int channel) const override; + + // method for returning PMT HV calibration curves from root file + virtual TF1 get_PMT_HV_curve(int channel) const override; + + TFile* HVgaincurves_rootFile; + + TF1* chan0_HVgaincurve_pntr; + TF1* chan1_HVgaincurve_pntr; + TF1* chan2_HVgaincurve_pntr; + TF1* chan3_HVgaincurve_pntr; + + TF1 m_HVgainCurveMap[4]; + + private: + // Propert that points to the location of the root file which contains the HV gain curves for the calorimeter PMTs + StringProperty m_PMT_HV_Gain_Curve_file{this, "PMT_HV_Gain_Curve_file", "/cvmfs/faser.cern.ch/repo/sw/database/DBRelease/current/pmtgain/CaloGainCurves.root"}; + + // properties that map channel name to PMT HV gain curves. PMt names found at https://twiki.cern.ch/twiki/bin/viewauth/FASER/CaloScintillator + StringProperty m_chan0_HVgaincurve_name{this, "Chan0HVgaincurve", "pol5_HV_Gain_Curve_PMT_LB8770"}; + StringProperty m_chan1_HVgaincurve_name{this, "Chan1HVgaincurve", "pol5_HV_Gain_Curve_PMT_LB8733"}; + StringProperty m_chan2_HVgaincurve_name{this, "Chan2HVgaincurve", "pol5_HV_Gain_Curve_PMT_LB8786"}; + StringProperty m_chan3_HVgaincurve_name{this, "Chan3HVgaincurve", "pol5_HV_Gain_Curve_PMT_LB8732"}; + + // Read Cond Handle + SG::ReadCondHandleKey<CondAttrListCollection> m_PMT_HV_ReadKey{this, "PMT_HV_ReadKey", "/WAVE/Calibration/HV", "Key of folder for PMT HV reading"}; + SG::ReadCondHandleKey<CondAttrListCollection> m_PMT_HV_ref_ReadKey{this, "PMT_HV_ref_ReadKey", "/WAVE/Calibration/HV_ref", "Key of folder for PMT HV at time of MIP calibration measurement"}; + SG::ReadCondHandleKey<CondAttrListCollection> m_MIPcharge_ref_ReadKey{this, "MIPcharge_ref_ReadKey", "/WAVE/Calibration/MIP_charge_ref", "Key of folder for MIPcharge calibration measurment"}; + +}; + +#endif // CALORECTOOLS_CALORECTOOL_H diff --git a/Calorimeter/CaloRecTools/src/components/CaloRecTools_entries.cxx b/Calorimeter/CaloRecTools/src/components/CaloRecTools_entries.cxx new file mode 100644 index 0000000000000000000000000000000000000000..c9c5438eb029be7532820e35dee8e41710f5e781 --- /dev/null +++ b/Calorimeter/CaloRecTools/src/components/CaloRecTools_entries.cxx @@ -0,0 +1,3 @@ +#include "../CaloRecTool.h" + +DECLARE_COMPONENT( CaloRecTool ) diff --git a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py index 4e233bf4fcb8a014393bf8858c55ef25d7690519..8a0b5a3cd15dc76526ae4b02d8efb4c2f3fb7adf 100755 --- a/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py +++ b/Control/CalypsoExample/Reconstruction/scripts/faser_reco.py @@ -180,9 +180,13 @@ 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)) +# Calorimeter reconstruction +if args.isMC: + # Not ready for MC quite yet + pass +else: + from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionCfg + acc.merge(CalorimeterReconstructionCfg(ConfigFlags)) # Tracker clusters from TrackerPrepRawDataFormation.TrackerPrepRawDataFormationConfig import FaserSCT_ClusterizationCfg @@ -252,8 +256,8 @@ from WaveRecAlgs.WaveRecAlgsConfig import WaveformReconstructionOutputCfg acc.merge(WaveformReconstructionOutputCfg(ConfigFlags)) # Calorimeter reconstruction output -# from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg -# acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags)) +from CaloRecAlgs.CaloRecAlgsConfig import CalorimeterReconstructionOutputCfg +acc.merge(CalorimeterReconstructionOutputCfg(ConfigFlags)) # Check what we have print( "Writing out xAOD objects:" ) diff --git a/xAOD/xAODFaserCalorimeter/Root/CalorimeterHitAuxContainer_v1.cxx b/xAOD/xAODFaserCalorimeter/Root/CalorimeterHitAuxContainer_v1.cxx index 6f6edf3e4c520dbf32bb816addb74f4ad3eae23f..f73a83dd127cc0d24d999c0e8f140e913a329e76 100644 --- a/xAOD/xAODFaserCalorimeter/Root/CalorimeterHitAuxContainer_v1.cxx +++ b/xAOD/xAODFaserCalorimeter/Root/CalorimeterHitAuxContainer_v1.cxx @@ -10,11 +10,11 @@ namespace xAOD { CalorimeterHitAuxContainer_v1::CalorimeterHitAuxContainer_v1() : AuxContainerBase() { - AUX_VARIABLE(localtime); - AUX_VARIABLE(bcid_time); - AUX_VARIABLE(raw_energy); + AUX_VARIABLE(channel); + AUX_VARIABLE(Edep); + AUX_VARIABLE(raw_Edep); - AUX_VARIABLE(caloLinks); + AUX_VARIABLE(WaveformLink); } } // namespace xAOD diff --git a/xAOD/xAODFaserCalorimeter/Root/CalorimeterHit_v1.cxx b/xAOD/xAODFaserCalorimeter/Root/CalorimeterHit_v1.cxx index 5e7a69ecc024dc3ecb70df95f1e623c4bc4e09fc..f085cb0c92acca759173260ba56c26813fd72935 100644 --- a/xAOD/xAODFaserCalorimeter/Root/CalorimeterHit_v1.cxx +++ b/xAOD/xAODFaserCalorimeter/Root/CalorimeterHit_v1.cxx @@ -13,41 +13,41 @@ 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, unsigned int, channel, set_channel ) - AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, bcidtime, set_bcidtime ) + AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, Edep, set_Edep ) - AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, raw_energy, set_raw_energy ) + AUXSTORE_PRIMITIVE_SETTER_AND_GETTER( CalorimeterHit_v1, float, raw_Edep, set_raw_Edep ) // setters and getters for the Calo WaveformHit links AUXSTORE_OBJECT_SETTER_AND_GETTER( CalorimeterHit_v1, CalorimeterHit_v1::WaveformHitLinks_t, - caloWaveformLinks, - setCaloWaveformLinks ) + WaveformLinks, + setWaveformLinks ) - static const SG::AuxElement::Accessor< CalorimeterHit_v1::WaveformHitLinks_t > caloHitAcc( "caloWaveformLinks" ); + static const SG::AuxElement::Accessor< CalorimeterHit_v1::WaveformHitLinks_t > HitAcc( "WaveformLinks" ); - const WaveformHit* CalorimeterHit_v1::caloHit( size_t i ) const { - return ( *caloHitAcc( *this )[ i ] ); + const WaveformHit* CalorimeterHit_v1::Hit( size_t i ) const { + return ( *HitAcc( *this )[ i ] ); } - size_t CalorimeterHit_v1::nCaloHits() const { - return caloHitAcc( *this ).size(); + size_t CalorimeterHit_v1::nHits() const { + return HitAcc( *this ).size(); } - void CalorimeterHit_v1::addCaloHit( const xAOD::WaveformHitContainer* pWaveformHitContainer, + void CalorimeterHit_v1::addHit( const xAOD::WaveformHitContainer* pWaveformHitContainer, const xAOD::WaveformHit* pWaveformHit) { ElementLink< xAOD::WaveformHitContainer > linkToWaveformHit; linkToWaveformHit.toContainedElement(*pWaveformHitContainer, pWaveformHit); - caloHitAcc( *this ).push_back( linkToWaveformHit ); + HitAcc( *this ).push_back( linkToWaveformHit ); return; } - void CalorimeterHit_v1::clearCaloWaveformLinks() { - caloHitAcc( *this ).clear(); + void CalorimeterHit_v1::clearWaveformLinks() { + HitAcc( *this ).clear(); return; } @@ -57,8 +57,8 @@ namespace xAOD { std::ostream& operator<<(std::ostream& s, const xAOD::CalorimeterHit_v1& hit) { s << "xAODCalorimeterHit:" - << " local time=" << hit.localtime() - << " raw_energy=" << hit.raw_energy() + << " channel = " << hit.channel() + << ", Edep = " << hit.Edep() << std::endl; return s; diff --git a/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h b/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h index 7e61261eb4bd10f106d03b2a85f9918b82ce7944..41a89a0325e9dd5ce78dab40ad6ae08458eb37c3 100644 --- a/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h +++ b/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHitAuxContainer_v1.h @@ -29,13 +29,14 @@ namespace xAOD { private: /// @name Basic variables ///@ { - std::vector<float> localtime; - std::vector<float> bcid_time; - std::vector<float> raw_energy; + std::vector<unsigned int> channel; + std::vector<float> Edep; + std::vector<float> raw_Edep; typedef std::vector< ElementLink< WaveformHitContainer > > WaveformHitLink_t; - std::vector< WaveformHitLink_t > caloLinks; - std::vector< WaveformHitLink_t > preshowerLinks; + std::vector< WaveformHitLink_t > WaveformLink; + //std::vector< WaveformHitLink_t > caloLinks; + //std::vector< WaveformHitLink_t > preshowerLinks; ///@} diff --git a/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHit_v1.h b/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHit_v1.h index d9bb35047e6a75d6dac12dc0ba0eb47e3bceb684..fefc0d5ae3dac013ee734e34ed82c600f8fbde5d 100644 --- a/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHit_v1.h +++ b/xAOD/xAODFaserCalorimeter/xAODFaserCalorimeter/versions/CalorimeterHit_v1.h @@ -35,29 +35,29 @@ namespace xAOD { /// @{ /// Best results - float localtime() const; - void set_localtime(float value); + unsigned int channel() const; + void set_channel(unsigned int value); - float bcidtime() const; - void set_bcidtime(float value); + float Edep() const; + void set_Edep(float value); - float raw_energy() const; - void set_raw_energy(float value); + float raw_Edep() const; + void set_raw_Edep(float value); // Waveform Hits typedef std::vector< ElementLink< xAOD::WaveformHitContainer > > WaveformHitLinks_t; // Contributing Calorimeter Waveform Hits - const WaveformHitLinks_t& caloWaveformLinks() const; - void setCaloWaveformLinks( const WaveformHitLinks_t& caloWaveforms ); + const WaveformHitLinks_t& WaveformLinks() const; + void setWaveformLinks( const WaveformHitLinks_t& Waveforms ); // Remove all waveform hits - void clearCaloWaveformLinks(); + void clearWaveformLinks(); // Get the pointer to a given waveform hit - const WaveformHit* caloHit( size_t i ) const; + const WaveformHit* Hit( size_t i ) const; // Get the number of waveform hits - size_t nCaloHits() const; + size_t nHits() const; // Add a waveform hit - void addCaloHit( const xAOD::WaveformHitContainer*, const xAOD::WaveformHit*); + void addHit( const xAOD::WaveformHitContainer*, const xAOD::WaveformHit*); // // Contributing Preshower Waveform Hits // const WaveformHitLinks_t& preshowerWaveformLinks() const;