diff --git a/DDDigi/include/DDDigi/DigiData.h b/DDDigi/include/DDDigi/DigiData.h index 5212f621d56293796109816519829acf56cd6c07..2ba807fdf77f51483ff24aab9b46eb62eef9aa87 100644 --- a/DDDigi/include/DDDigi/DigiData.h +++ b/DDDigi/include/DDDigi/DigiData.h @@ -501,21 +501,23 @@ namespace dd4hep { }; /// Hit position - Position position { }; + Position position { }; /// Hit direction - Direction momentum { }; + Direction momentum { }; /// Length of the track segment contributing to this hit - double length { 0 }; + double length { 0e0 }; /// Total energy deposit - double deposit { 0 }; + double deposit { 0e0 }; + /// Total energy deposit + double depositError { -1e0 }; /// Proper creation time of the deposit with rescpect to beam crossing - double time { 0 }; + double time { 0e0 }; /// Optional flag for user masks - uint64_t flag { 0UL }; + uint64_t flag { 0UL }; /// Source mask of this deposit - Key::mask_type mask { 0 }; + Key::mask_type mask { 0 }; /// Deposit history - History history { }; + History history { }; public: /// Default constructor diff --git a/DDDigi/io/Digi2edm4hep.cpp b/DDDigi/io/Digi2edm4hep.cpp index a32a64d78f139f3883611367b94698157f6e6768..564f582ae754e54fffa25714b843acfbdd3cfe2f 100644 --- a/DDDigi/io/Digi2edm4hep.cpp +++ b/DDDigi/io/Digi2edm4hep.cpp @@ -91,11 +91,10 @@ namespace dd4hep { virtual void execute(context_t& context) const; }; - /// Actor to select energy deposits according to the supplied segmentation - /** Actor to select energy deposits according to the supplied segmentation + /// Actor to save individual data containers to edm4hep + /** Actor to save individual data containers to edm4hep * - * The selected deposits are placed in the output container - * supplied by the arguments. + * This is a typical worker action of the Digi2edm4hepAction * * \author M.Frank * \version 1.0 @@ -105,9 +104,14 @@ namespace dd4hep { friend class Digi2edm4hepAction; protected: + /// Reference to the edm4hep engine std::shared_ptr<Digi2edm4hepAction::internals_t> internals; - float pointResoRPhi = 0.004; - float pointResoZ = 0.004; + /// Property: RPhi resolution + float m_pointResoutionRPhi = 0.004; + /// Property: Z resolution + float m_pointResoutionZ = 0.004; + /// Hit type for hit processor + int m_hit_type = 0; public: /// Standard constructor @@ -173,14 +177,24 @@ namespace dd4hep { /// Namespace for the Digitization part of the AIDA detector description toolkit namespace digi { + /// Helper class to create output in edm4hep format + /** Helper class to create output in edm4hep format + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ class Digi2edm4hepAction::internals_t { public: - Digi2edm4hepAction* m_parent { nullptr }; - std::unique_ptr<podio::EventStore> m_store { }; - std::unique_ptr<podio::ROOTWriter> m_file { }; - edm4hep::EventHeaderCollection* m_header { nullptr }; + Digi2edm4hepAction* m_parent { nullptr }; + /// Reference to podio store + std::unique_ptr<podio::EventStore> m_store { }; + /// Reference to podio writer + std::unique_ptr<podio::ROOTWriter> m_file { }; + /// edm4hep event header collection + edm4hep::EventHeaderCollection* m_header { nullptr }; /// MC particle collection - edm4hep::MCParticleCollection* m_particles { nullptr }; + edm4hep::MCParticleCollection* m_particles { nullptr }; /// Collection of all edm4hep object collections std::map<std::string, podio::CollectionBase*> m_collections; @@ -382,6 +396,9 @@ namespace dd4hep { Digi2edm4hepProcessor::Digi2edm4hepProcessor(const DigiKernel& krnl, const std::string& nam) : DigiContainerProcessor(krnl, nam) { + declareProperty("point_resolution_RPhi", m_pointResoutionRPhi); + declareProperty("point_resolution_Z", m_pointResoutionZ); + declareProperty("hit_type", m_hit_type = 0); } void Digi2edm4hepProcessor::convert_particles(DigiContext& ctxt, @@ -401,10 +418,12 @@ namespace dd4hep { const predicate_t& predicate, edm4hep::TrackerHitCollection* collection) const { - std::array<float,6> covMat = {0., 0., pointResoRPhi*pointResoRPhi, 0., 0., pointResoZ*pointResoZ}; + std::array<float,6> covMat = {0., 0., m_pointResoutionRPhi*m_pointResoutionRPhi, + 0., 0., m_pointResoutionZ*m_pointResoutionZ + }; for ( const auto& depo : cont ) { if ( predicate(depo) ) { - digi_io::_to_edm4hep(depo, covMat, collection); + digi_io::_to_edm4hep(depo, covMat, m_hit_type /* edm4hep::SIMTRACKERHIT */, collection); } } } @@ -416,7 +435,7 @@ namespace dd4hep { { for ( const auto& depo : cont ) { if ( predicate(depo) ) { - digi_io::_to_edm4hep(depo, collection); + digi_io::_to_edm4hep(depo, m_hit_type /* edm4hep::SIMCALORIMETERHIT */, collection); } } } diff --git a/DDDigi/io/DigiEdm4hepInput.cpp b/DDDigi/io/DigiEdm4hepInput.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4e39c45c3e1aa7863fcaab3694c3c6194935339b --- /dev/null +++ b/DDDigi/io/DigiEdm4hepInput.cpp @@ -0,0 +1,126 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include <DDDigi/DigiROOTInput.h> +#include <DDDigi/DigiFactories.h> +#include "DigiIO.h" + +#include <DDG4/Geant4Data.h> + +// ROOT include files +#include <TBranch.h> +#include <TClass.h> + +// C/C++ include files + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + using namespace std::placeholders; + + class DigiEdm4hepROOT : public DigiROOTInput { + public: + static constexpr double epsilon = std::numeric_limits<double>::epsilon(); + + /// Property to generate extra history records + bool m_keep_raw { true }; + mutable TClass* m_trackerHitClass { nullptr }; + mutable TClass* m_caloHitClass { nullptr }; + mutable TClass* m_particlesClass { nullptr }; + + TClass* get_class(TClass* c) const { + if ( c == m_particlesClass || c == m_trackerHitClass || c == m_caloHitClass ) + return c; + const char* nam = c->GetName(); + if ( strstr(nam, "vector<dd4hep::sim::Geant4Particle*") ) + return m_particlesClass = c; + else if ( strstr(nam, "vector<dd4hep::sim::Geant4Tracker::Hit*") ) + return m_trackerHitClass = c; + else if ( strstr(nam, "vector<dd4hep::sim::Geant4Calorimeter::Hit*") ) + return m_caloHitClass = c; + except("Unknown data class: %s Cannot be handled", nam); + return c; + } + + public: + /// Initializing constructor + DigiEdm4hepROOT(const DigiKernel& krnl, const std::string& nam) + : DigiROOTInput(krnl, nam) + { + declareProperty("keep_raw", m_keep_raw); + } + + template <typename T> + void conv_hits(DigiContext& context, DataSegment& segment, + const std::string& tag, + Key::mask_type mask, + const char* nam, + void* ptr) const + { + DepositVector out(nam, mask, SegmentEntry::UNKNOWN); + std::map<CellID, std::shared_ptr<T> > hits; + std::size_t len = 0; + if ( ptr ) { + input_data<T> data(ptr); + const DepositPredicate<EnergyCut> predicate ({ this->epsilon }); + len = data.size(); + data_io<ddg4_input>()._to_digi_if(data.get(), hits, predicate); + data_io<ddg4_input>()._to_digi(Key(nam, segment.id, mask), hits, out); + data.clear(); + } + info("%s+++ %-24s Converted %6ld Edm4hep %-14s hits to %6ld cell deposits", + context.event->id(), nam, len, tag.c_str(), out.size()); + put_data(segment, Key(out.name, mask), out); + if ( m_keep_raw ) { + put_data(segment, Key(std::string(nam)+".ddg4", mask, segment.id), hits); + } + } + + void conv_particles(DigiContext& context, DataSegment& segment, + int mask, const std::string& nam, void* ptr) const + { + ParticleMapping particles(nam, mask); + if ( ptr ) { + input_data<sim::Geant4Particle> data(ptr); + data_io<ddg4_input>()._to_digi(Key(nam, segment.id, mask), data.get(), particles); + data.clear(); + } + debug("%s+++ Converted %ld Edm4hep particles", context.event->id(), particles.size()); + put_data(segment, Key(nam, mask), particles); + } + + /// Callback to handle single branch + virtual void operator()(DigiContext& context, work_t& work) const override { + TBranch& br = work.branch; + int msk = work.input_key.mask(); + auto& seg = work.input_segment; + void** add = (void**)br.GetAddress(); + TClass* cls = get_class(&work.cl); + const char* nam = br.GetName(); + if ( cls == m_caloHitClass ) + conv_hits<sim::Geant4Calorimeter::Hit>(context, seg, "calorimeter", msk, nam, *add); + else if ( cls == m_trackerHitClass ) + conv_hits<sim::Geant4Tracker::Hit>(context, seg, "tracker", msk, nam, *add); + else if ( cls == m_particlesClass ) + conv_particles(context, seg, msk, nam, *add); + else + except("Unknown data type encountered in branch: %s", nam); + } + }; + } // End namespace digi +} // End namespace dd4hep +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiEdm4hepROOT) diff --git a/DDDigi/io/DigiIO.cpp b/DDDigi/io/DigiIO.cpp index b54490689995ad439756cfe92e457f7b908ab5f5..73023c497e3b6bc33836e33fc528f70a0aa58729 100644 --- a/DDDigi/io/DigiIO.cpp +++ b/DDDigi/io/DigiIO.cpp @@ -16,6 +16,9 @@ #include <DD4hep/Printout.h> #include "DigiIO.h" +/// C/C++ include files +#include <limits> + // ========================================================================= // EDM4HEP specific stuff // ========================================================================= @@ -354,8 +357,8 @@ namespace dd4hep { cnv_to_digi(key, p, out); } - } // End namespace digi -} // End namespace dd4hep + } // End namespace digi +} // End namespace dd4hep #endif // DD4HEP_USE_DDG4 /// Namespace for the AIDA detector description toolkit @@ -398,33 +401,43 @@ namespace dd4hep { template <> void digi_io::_to_edm4hep(const std::pair<const CellID, EnergyDeposit>& dep, const std::array<float, 6>& covMat, + int hit_type, edm4hep::TrackerHitCollection* collection) { const EnergyDeposit& de = dep.second; auto hit = collection->create(); - hit.setType(0 /* edm4hep::SIMTRACKERHIT */); - hit.setTime(de.time); - hit.setCovMatrix(covMat); - hit.setCellID(dep.first); - hit.setEDep(de.deposit); - hit.setEDepError(de.deposit/1e0); - hit.setEdx(de.deposit/de.length); + double dep_error = de.depositError; + if ( dep_error < -std::numeric_limits<double>::epsilon() ) { + dep_error = 0e0; + } + hit.setType( hit_type ); + hit.setTime( de.time ); + hit.setCovMatrix( covMat ); + hit.setCellID( dep.first ); + hit.setEDep( de.deposit ); + hit.setEDepError( dep_error ); + hit.setEdx( de.deposit/de.length ); hit.setPosition( _toVectorD(de.position) ); } template <> void digi_io::_to_edm4hep(const std::pair<const CellID, EnergyDeposit>& dep, + int hit_type, edm4hep::CalorimeterHitCollection* collection) { const EnergyDeposit& de = dep.second; auto hit = collection->create(); - hit.setType(0 /* edm4hep::SIMCALORIMETERHIT */); - hit.setTime(de.time); - hit.setCellID(dep.first); - hit.setEnergy(de.deposit); - hit.setEnergyError(de.deposit/10e0); + double dep_error = de.depositError; + if ( dep_error < -std::numeric_limits<double>::epsilon() ) { + dep_error = 0e0; + } + hit.setType( hit_type ); + hit.setTime( de.time ); + hit.setCellID( dep.first ); + hit.setEnergy( de.deposit ); + hit.setEnergyError( dep_error ); hit.setPosition( _toVectorF(de.position) ); } - } // End namespace digi -} // End namespace dd4hep + } // End namespace digi +} // End namespace dd4hep #endif // DD4HEP_USE_EDM4HEP diff --git a/DDDigi/io/DigiIO.h b/DDDigi/io/DigiIO.h index 84cd054ccfaa9734a332605e335d4ec91886a2e2..b429c62b952dc0804bac08e3e2a464d6c0c68cd4 100644 --- a/DDDigi/io/DigiIO.h +++ b/DDDigi/io/DigiIO.h @@ -74,8 +74,11 @@ namespace dd4hep { template <typename FIRST_TYPE, typename OUTPUT_TYPE> static void _to_edm4hep(const FIRST_TYPE& first, OUTPUT_TYPE output); + template <typename FIRST_TYPE, typename OUTPUT_TYPE> static + void _to_edm4hep(const FIRST_TYPE& first, int hit_type, OUTPUT_TYPE output); + template <typename FIRST_TYPE, typename SECOND_TYPE, typename OUTPUT_TYPE> static - void _to_edm4hep(const FIRST_TYPE& first, const SECOND_TYPE& second, OUTPUT_TYPE output); + void _to_edm4hep(const FIRST_TYPE& first, const SECOND_TYPE& second, int hit_type, OUTPUT_TYPE output); }; struct ddg4_input { @@ -118,6 +121,9 @@ namespace dd4hep { template <typename FIRST, typename SECOND> static void _to_edm4hep(const FIRST& cont, SECOND coll); + template <typename FIRST, typename SECOND> static + void _to_edm4hep(const FIRST& cont, SECOND coll, int hit_type); + template <typename FIRST, typename SECOND, typename THIRD> static void _to_digi(FIRST first, const SECOND& second, THIRD& third); diff --git a/DDDigi/plugins/DigiDepositSmearEnergy.cpp b/DDDigi/plugins/DigiDepositSmearEnergy.cpp index 25e8e6b4290d6b2294b54194ec49738059f1dc09..66a781707e7e974ca7857a1f562facc99bef17f4 100644 --- a/DDDigi/plugins/DigiDepositSmearEnergy.cpp +++ b/DDDigi/plugins/DigiDepositSmearEnergy.cpp @@ -57,6 +57,13 @@ namespace dd4hep { * * (See all: https://handwiki.org/wiki/Physics:Energy_resolution_in_calorimeters) * + * Properties: + * intrinsic_fluctuation: Intrinsic energy resolution constant (gaussian ~ std::sqrt(energy)) + * systematic_resolution: Systematic energy resolution constant (gaussian ~deposit energy) + * instrumentation_resolution: Instrumentation energy resolution constant (gaussian(CONSTANT)) + * pair_ionization_energy: Ionization constant to create electron ION pair in absorber + * ionization_fluctuation: Flag to use ionization fluctuation smearing (poisson ~ # e-ion-pairs) + * modify_energy: Flag to modify energy deposit according to error - otherwise only compute error * * \author M.Frank * \version 1.0 @@ -74,6 +81,8 @@ namespace dd4hep { double m_pair_ionization_energy { 0e0 }; /// Property: Flag to use ionization fluctuation smearing (poisson ~ # e-ion-pairs) bool m_ionization_fluctuation { false }; + /// Property: Flag to modify energy deposit according to error - otherwise only compute error + bool m_modify_energy { true }; public: /// Standard constructor @@ -85,6 +94,7 @@ namespace dd4hep { declareProperty("systematic_resolution", m_systematic_resolution = 0e0); declareProperty("pair_ionisation_energy", m_pair_ionization_energy = 3.6*dd4hep::eV); declareProperty("ionization_fluctuation", m_ionization_fluctuation = false); + declareProperty("modify_energy", m_modify_energy = true); DEPOSIT_PROCESSOR_BIND_HANDLERS(DigiDepositSmearEnergy::smear) } @@ -127,9 +137,14 @@ namespace dd4hep { } /// delta_E is in GeV delta_E *= dd4hep::GeV; - if ( m_monitor ) m_monitor->energy_shift(dep, delta_E); - depo.deposit = deposit + (delta_E); - depo.flag |= EnergyDeposit::ENERGY_SMEARED; + depo.depositError = delta_E; + if ( m_monitor ) { + m_monitor->energy_shift(dep, delta_E); + } + if ( m_modify_energy ) { + depo.deposit = deposit + delta_E; + depo.flag |= EnergyDeposit::ENERGY_SMEARED; + } ++updated; } } @@ -137,8 +152,26 @@ namespace dd4hep { context.event->id(), cont.name.c_str(), updated, cont.size(), cont.key.mask()); } }; + + /// Actor to only set energy error (as above, but with preset option + /** + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiDepositSetEnergyError : public DigiDepositSmearEnergy { + public: + /// Standard constructor + DigiDepositSetEnergyError(const DigiKernel& krnl, const std::string& nam) + : DigiDepositSmearEnergy(krnl, nam) + { + m_modify_energy = false; + } + }; } // End namespace digi } // End namespace dd4hep #include <DDDigi/DigiFactories.h> DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearEnergy) +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSetEnergyError) diff --git a/DDDigi/python/dddigi.py b/DDDigi/python/dddigi.py index 4d35ef7b553f5136e29f1e8094bd8595666d2ab3..80922be74b995c1ef375b468ea31592472cccd78 100644 --- a/DDDigi/python/dddigi.py +++ b/DDDigi/python/dddigi.py @@ -335,7 +335,6 @@ _props('DigiSequentialActionSequence', adopt_action=_adopt_sequence_action) _props('DigiContainerSequenceAction', adopt_container_processor=_adopt_container_processor) _props('DigiMultiContainerProcessor', adopt_processor=_adopt_processor) _props('DigiSegmentSplitter', adopt_segment_processor=_adopt_segment_processor) - # --------------------------------------------------------------------------- diff --git a/DDDigi/src/DigiData.cpp b/DDDigi/src/DigiData.cpp index f89f3560d61cf11d6b14197fe1958318ace6a589..fac5b73cf3e64d00806e0552269eaf22795f7e29 100644 --- a/DDDigi/src/DigiData.cpp +++ b/DDDigi/src/DigiData.cpp @@ -157,23 +157,27 @@ Direction History::average_particle_momentum(const DigiEvent& event) const { /// Update the deposit using deposit weighting void EnergyDeposit::update_deposit_weighted(const EnergyDeposit& upda) { double sum = deposit + upda.deposit; + double err = depositError + upda.depositError; Position pos = ((deposit / sum) * position) + ((upda.deposit / sum) * upda.position); Direction mom = ((deposit / sum) * momentum) + ((upda.deposit / sum) * upda.momentum); position = pos; momentum = mom; deposit = sum; + depositError = err; history.update(upda.history); } /// Update the deposit using deposit weighting void EnergyDeposit::update_deposit_weighted(EnergyDeposit&& upda) { double sum = deposit + upda.deposit; + double err = depositError + upda.depositError; Position pos = ((deposit / sum) * position) + ((upda.deposit / sum) * upda.position); Direction mom = ((deposit / sum) * momentum) + ((upda.deposit / sum) * upda.momentum); position = pos; momentum = mom; - deposit = sum; - history.update(upda.history); + deposit = sum; + depositError = err; + history.update(std::move(upda.history)); } /// Merge new deposit map onto existing map diff --git a/examples/DDDigi/scripts/DigiTest.py b/examples/DDDigi/scripts/DigiTest.py index 34db86b5161fecfa195f88380e0c61a268c04738..651b834c89316cd41f56f81398e685b712bdc46c 100644 --- a/examples/DDDigi/scripts/DigiTest.py +++ b/examples/DDDigi/scripts/DigiTest.py @@ -120,12 +120,13 @@ class Test(dddigi.Digitize): evt_done = self.run(num_events=num_events, num_threads=num_threads, parallel=parallel) if evt_done == num_events: result = "PASSED" - self.always('%s Test finished after processing %d events. [%d parallel threads, %d parallel events]' \ + self.always('%s Test finished after processing %d events. [%d parallel threads, %d parallel events]' % (result, evt_done, num_threads, parallel, )) self.always('Test done. Exiting') self.kernel().terminate() return evt_done + # ========================================================================================================== def test_setup_1(digi, print_level=WARNING, parallel=True): """ diff --git a/examples/DDDigi/scripts/TestSpillover.py b/examples/DDDigi/scripts/TestSpillover.py index 6986b143313ce11a33f5268e5d4bae39b4f03054..b60267bd3126d9d27908c72c6b0b7361996eecb9 100644 --- a/examples/DDDigi/scripts/TestSpillover.py +++ b/examples/DDDigi/scripts/TestSpillover.py @@ -31,7 +31,7 @@ def run(): # ======================================================================================================== input_action.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) digi.info('Created input.signal') @@ -41,7 +41,7 @@ def run(): spillover = input_action.adopt_action('DigiSequentialActionSequence/Spillover-25') evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-25ns', mask=0x1, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-25ns', @@ -58,7 +58,7 @@ def run(): spillover = input_action.adopt_action('DigiSequentialActionSequence/Spillover-50') evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-50ns', mask=0x2, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-50ns', @@ -73,7 +73,7 @@ def run(): spillover = input_action.adopt_action('DigiSequentialActionSequence/Spillover-75') evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-75ns', mask=0x3, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-75ns', @@ -105,7 +105,7 @@ def run(): spillover = input_action.adopt_action('DigiSequentialActionSequence/Spillover+50') evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader+50ns', mask=0x5, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att+50ns', @@ -120,7 +120,7 @@ def run(): spillover = input_action.adopt_action('DigiSequentialActionSequence/Spillover+75') evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader+75ns', mask=0x6, - keep_raw = keep_raw, + keep_raw=keep_raw, input=[digi.next_input()], OutputLevel=rdr_output) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att+75ns',