diff --git a/DDDigi/ddg4/DigiDDG4Input.cpp b/DDDigi/ddg4/DigiDDG4Input.cpp index e8c65782dfbdd51f0d6da30c67e33ed643a099b2..9c9c5101fe92ada613847a3905352d23f8251824 100644 --- a/DDDigi/ddg4/DigiDDG4Input.cpp +++ b/DDDigi/ddg4/DigiDDG4Input.cpp @@ -12,55 +12,81 @@ //========================================================================== // Framework include files -#include <DD4hep/Printout.h> -#include <DD4hep/Factories.h> #include <DD4hep/DD4hepUnits.h> + #include <DDG4/Geant4Data.h> #include <DDG4/Geant4Particle.h> + #include <DDDigi/DigiData.h> -#include <CLHEP/Units/SystemOfUnits.h> +#include <DDDigi/DigiROOTInput.h> +#include <DDDigi/DigiFactories.h> // ROOT include files +#include <TBranch.h> // C/C++ include files #include <stdexcept> #include <any> -using namespace dd4hep; -using namespace dd4hep::digi; - -typedef std::function<void(DataSegment& segment, int, const char*, void*)> func_t; - -namespace { - - double epsilon = std::numeric_limits<double>::epsilon(); +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + class DigiDDG4ROOT : public DigiROOTInput { + public: + static constexpr double epsilon = std::numeric_limits<double>::epsilon(); + /// Property to generate extra history records + bool m_seperate_history { false }; + mutable TClass* m_trackerHitClass { nullptr }; + mutable TClass* m_caloHitClass { nullptr }; + mutable TClass* m_particlesClass { nullptr }; + + template <typename T> union input_data { + const void* raw; + std::vector<T*>* items; + input_data(const void* p) { this->raw = p; } + }; + + TClass* get_class_pointers(TClass* c , const char* nam) const { + if ( strstr(nam, "vector<dd4hep::sim::Geant4Particle*") ) + return m_particlesClass = c; + if ( strstr(nam, "vector<dd4hep::sim::Geant4Tracker::Hit*") ) + return m_trackerHitClass = c; + if ( strstr(nam, "vector<dd4hep::sim::Geant4Calorimeter::Hit*") ) + return m_caloHitClass = c; + except("Unknown data class: %s Cannot be handled", nam); + return c; + } - template <typename T> union input_data { - const void* raw; - std::vector<T*>* items; - input_data(const void* p) { this->raw = p; } - }; + void add_particle_history(const sim::Geant4Calorimeter::Hit* hit, Key key, History& hist) const { + for( const auto& truth : hit->truth ) { + key.set_item(truth.trackID); + hist.particles.emplace_back(key, truth.deposit); + } + } + void add_particle_history(const sim::Geant4Tracker::Hit* hit, Key key, History& hist) const { + key.set_item(hit->truth.trackID); + hist.particles.emplace_back(key, hit->truth.deposit); + } - template <typename T> void add_particle_history(const T* hit, Key key, History& hist); + public: + /// Initializing constructor + DigiDDG4ROOT(const DigiKernel& krnl, const std::string& nam) + : DigiROOTInput(krnl, nam) + { + declareProperty("seperate_history", m_seperate_history); + } - template <> void add_particle_history(const sim::Geant4Calorimeter::Hit* hit, Key key, History& hist) { - for( const auto& truth : hit->truth ) { - key.set_item(truth.trackID); - hist.particles.emplace_back(key, truth.deposit); - } - } - template <> void add_particle_history(const sim::Geant4Tracker::Hit* hit, Key key, History& hist) { - key.set_item(hit->truth.trackID); - hist.particles.emplace_back(key, hit->truth.deposit); - } - template <typename T> void* ddg4_hit_convert_function(const char* desc) { - std::string tag = desc; - func_t* cnv = new func_t([tag] (DataSegment& segment, int mask, const char* name, void* ptr) { + template <typename T> + void conv_hits(DataSegment& segment, const std::string& tag, int mask, const char* name, void* ptr) const { using wrap_t = std::shared_ptr<sim::Geant4HitData>; std::size_t len = 0; std::string nam = name; std::vector<wrap_t> geant4_hits; - DepositVector out(name, mask); + DepositVector out(nam, mask); + DetectorHistory hist(nam+".hist", mask); if ( ptr ) { input_data<T> data(ptr); for(size_t i=0; i < data.items->size(); ++i) { @@ -68,79 +94,94 @@ namespace { if ( p->energyDeposit > epsilon ) { Key history_key; EnergyDeposit dep { }; + Position pos = p->position; + pos *= 1./dd4hep::mm; + dep.flag = p->flag; dep.deposit = p->energyDeposit; - dep.position = p->position; + dep.position = pos; history_key.set_mask(Key::mask_type(mask)); history_key.set_item(geant4_hits.size()); history_key.set_segment(segment.id); - dep.history.hits.emplace_back(history_key, p->energyDeposit); + dep.history.hits.emplace_back(history_key, dep.deposit); add_particle_history(p, history_key, dep.history); out.emplace(p->cellID, std::move(dep)); + if ( m_seperate_history ) { + hist.insert(p->cellID, dep.history); + } geant4_hits.emplace_back(wrap_t(p)); } } len = data.items->size(); data.items->clear(); } - printout(DEBUG,"DDG4Converter","++ Converted %ld %s to %ld cell deposits", - len, tag.c_str(), out.size()); - Key depo_key(nam, mask); + debug("++ Converted %ld %s to %ld cell deposits", len, tag.c_str(), out.size()); + Key depo_key(out.name, mask); segment.emplace(depo_key, std::make_any<DepositVector>(std::move(out))); Key src_key(nam+".ddg4", mask); segment.emplace(src_key, std::make_any<std::vector<wrap_t>>(std::move(geant4_hits))); - }); - return cnv; - } -} - -static void* convert_sim_geant4calorimeter_hits() { - return ddg4_hit_convert_function<sim::Geant4Calorimeter::Hit>("DDG4 calorimeter hits"); -} -DECLARE_CREATE(DD4hep_DDDigiConverter_vector_dd4hep_sim_Geant4Calorimeter_Hit_,convert_sim_geant4calorimeter_hits) - -static void* convert_sim_geant4tracker_hits() { - return ddg4_hit_convert_function<sim::Geant4Tracker::Hit>("DDG4 tracker hits"); -} -DECLARE_CREATE(DD4hep_DDDigiConverter_vector_dd4hep_sim_Geant4Tracker_Hit_,convert_sim_geant4tracker_hits) - -static void* convert_sim_geant4particles() { - func_t* cnv = new func_t([] (DataSegment& segment, int mask, const char* name, void* ptr) { - using wrap_t = std::shared_ptr<sim::Geant4Particle>; - ParticleMapping out(name, mask); - std::map<Key, wrap_t> source; - Key part_key(name, mask); - if ( ptr ) { - auto* items = (std::vector<sim::Geant4Particle*>*)ptr; - Key key; - key.set_segment(segment.id); - key.set_mask(mask); - for( auto* p : *items ) { - Particle part; - p->mask = mask; - part.start_position = Position(p->vsx, p->vsy, p->vsz); - part.end_position = Position(p->vex, p->vey, p->vez); - part.momentum = Direction(p->psx, p->psy, p->psz); - part.pdgID = p->pdgID; - part.charge = p->charge; - part.mass = p->mass; - - part.history = key.set_item(part_key.item()); - key.set_item(out.size()); - out.push(key, std::move(part)); - source[key] = wrap_t(p); + if ( m_seperate_history ) { + Key hist_key(hist.name, mask); + segment.emplace(hist_key, std::make_any<DetectorHistory>(std::move(hist))); + } + } + + void conv_particles(DataSegment& segment, int mask, const char* name, void* ptr) const { + using wrap_t = std::shared_ptr<sim::Geant4Particle>; + ParticleMapping out(name, mask); + std::map<Key, wrap_t> source; + Key part_key(name, mask); + if ( ptr ) { + auto* items = (std::vector<sim::Geant4Particle*>*)ptr; + Key key; + key.set_segment(segment.id); + key.set_mask(mask); + for( auto* p : *items ) { + Particle part; + p->mask = mask; + part.start_position = Position(p->vsx, p->vsy, p->vsz); + part.end_position = Position(p->vex, p->vey, p->vez); + part.momentum = Direction(p->psx, p->psy, p->psz); + part.pdgID = p->pdgID; + part.charge = p->charge; + part.mass = p->mass; + + part.history = key.set_item(part_key.item()); + key.set_item(out.size()); + out.push(key, std::move(part)); + source[key] = wrap_t(p); + } + items->clear(); } - items->clear(); + debug("++ Converted %ld DDG4 particles", out.size()); + std::string nam = name; + segment.emplace(part_key, std::make_any<ParticleMapping>(std::move(out))); + Key src_key(nam+".ddg4", mask); + segment.emplace(src_key, std::make_any<std::map<Key, wrap_t> >(std::move(source))); + } + + /// Callback to handle single branch + virtual void operator()(DigiContext& /* context */, work_t& work) const override { + TBranch& br = work.branch; + TClass& cl = work.cl; + TClass* c = get_class_pointers(&cl, br.GetClassName()); + void** addr = (void**)br.GetAddress(); + const char* nam = br.GetName(); + int mask = work.input_key.mask(); + auto& seg = work.input_segment; + if ( c == m_caloHitClass ) + conv_hits<sim::Geant4Calorimeter::Hit>(seg, "DDG4 calorimeter hits", mask, nam, *addr); + else if ( c == m_trackerHitClass ) + conv_hits<sim::Geant4Tracker::Hit>(seg, "DDG4 tracker hits", mask, nam, *addr); + else if ( c == m_particlesClass ) + conv_particles(seg, mask, nam, *addr); + else + except("Unknown data type encountered in branch: %s", nam); } - printout(DEBUG,"DDG4Converter","++ Converted %ld DDG4 particles", out.size()); - std::string nam = name; - segment.emplace(part_key, std::make_any<ParticleMapping>(std::move(out))); - Key src_key(nam+".ddg4", mask); - segment.emplace(src_key, std::make_any<std::map<Key, wrap_t> >(std::move(source))); - }); - return cnv; -} -DECLARE_CREATE(DD4hep_DDDigiConverter_vector_dd4hep_sim_Geant4Particle_,convert_sim_geant4particles) + }; + } // End namespace digi +} // End namespace dd4hep +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDDG4ROOT) diff --git a/DDDigi/include/DDDigi/DigiData.h b/DDDigi/include/DDDigi/DigiData.h index f31bf506e3f171fea62c37419e898706b09d9900..f31d7715e44d5563931664b2752507f1508bde42 100644 --- a/DDDigi/include/DDDigi/DigiData.h +++ b/DDDigi/include/DDDigi/DigiData.h @@ -19,6 +19,7 @@ /// C/C++ include files #include <cstdint> #include <memory> +#include <limits> #include <mutex> #include <map> #include <any> @@ -26,6 +27,10 @@ /// Namespace for the AIDA detector description toolkit namespace dd4hep { + namespace detail { + static constexpr double numeric_epsilon = 10.0 * std::numeric_limits<double>::epsilon(); + } + /// Namespace for the Digitization part of the AIDA detector description toolkit namespace digi { @@ -410,8 +415,8 @@ namespace dd4hep { hist_entry_t& operator=(hist_entry_t&& copy) = default; hist_entry_t& operator=(const hist_entry_t& copy) = default; ~hist_entry_t() = default; - const Particle& get_particle(DigiEvent& event) const; - const EnergyDeposit& get_deposit(DigiEvent& event, Key::itemkey_type container_item) const; + const Particle& get_particle(const DigiEvent& event) const; + const EnergyDeposit& get_deposit(const DigiEvent& event, Key::itemkey_type container_item) const; }; /// Sources contributing to the deposit indexed by the cell identifier std::vector<hist_entry_t> hits; @@ -435,6 +440,17 @@ namespace dd4hep { void update(const History& upda); /// Drop history information std::pair<std::size_t,std::size_t> drop(); + + /// Number of hits contributing to history entry + std::size_t num_hits() const { + return hits.size(); + } + /// Number of particles contributing to history entry + std::size_t num_particles() const { + return particles.size(); + } + /// Retrieve the weighted momentum of all contributing particles + Direction average_particle_momentum(const DigiEvent& event) const; }; inline History::hist_entry_t::hist_entry_t(Key s, double w) @@ -450,6 +466,12 @@ namespace dd4hep { */ class EnergyDeposit { public: + enum { + KILLED = 1 << 0, + ENERGY_SMEARED = 1 << 1, + POSITION_SMEARED = 1 << 2 + }; + /// Hit position Position position { }; /// Hit direction @@ -549,6 +571,8 @@ namespace dd4hep { const_iterator begin() const { return this->data.begin(); } /// End iteration (CONST) const_iterator end() const { return this->data.end(); } + /// Remove entry + void remove(iterator position); }; /// Initializing constructor @@ -602,6 +626,8 @@ namespace dd4hep { std::size_t merge(DepositVector&& updates); /// Merge new deposit map onto existing map (not thread safe!) std::size_t insert(const DepositVector& updates); + /// Emplace entry + void emplace(CellID cell, EnergyDeposit&& deposit); /// Access container size std::size_t size() const { return this->data.size(); } @@ -619,6 +645,8 @@ namespace dd4hep { const_iterator begin() const { return this->data.begin(); } /// End iteration (CONST) const_iterator end() const { return this->data.end(); } + /// Remove entry + void remove(iterator position); }; /// Initializing constructor @@ -942,15 +970,19 @@ namespace dd4hep { const char* id() const { return this->m_id.c_str(); } /// Retrieve data segment from the event structure by name DataSegment& get_segment(const std::string& name); + /// Retrieve data segment from the event structure by name (CONST) + const DataSegment& get_segment(const std::string& name) const; /// Retrieve data segment from the event structure by identifier DataSegment& get_segment(Key::segment_type id); + /// Retrieve data segment from the event structure by identifier (CONST) + const DataSegment& get_segment(Key::segment_type id) const; }; /// Static global functions std::string digiTypeName(const std::type_info& info); std::string digiTypeName(const std::any& data); - const Particle& get_history_particle(DigiEvent& event, Key history_key); - const EnergyDeposit& get_history_deposit(DigiEvent& event, Key::itemkey_type container_item, Key history_key); + const Particle& get_history_particle(const DigiEvent& event, Key history_key); + const EnergyDeposit& get_history_deposit(const DigiEvent& event, Key::itemkey_type container_item, Key history_key); } // End namespace digi } // End namespace dd4hep diff --git a/DDDigi/include/DDDigi/DigiInputAction.h b/DDDigi/include/DDDigi/DigiInputAction.h index 87fc6496692b75dfaad23fd5b16cb50c53833311..584775b0768599a7be62ea29cf0b7c7c57bbc6a3 100644 --- a/DDDigi/include/DDDigi/DigiInputAction.h +++ b/DDDigi/include/DDDigi/DigiInputAction.h @@ -39,9 +39,13 @@ namespace dd4hep { protected: /// Property: Input data specification - std::vector<std::string> m_inputs { }; + std::vector<std::string> m_input_sources { }; + /// Property: Input data segment name + std::string m_input_segment { "inputs" }; /// Property: Mask to flag input source items - int m_mask { NO_MASK }; + int m_input_mask { NO_MASK }; + /// Property: Loop on inputs and restart at end + bool m_input_rescan { true }; protected: /// Define standard assignments and constructors @@ -53,6 +57,18 @@ namespace dd4hep { /// Default destructor virtual ~DigiInputAction(); + /// Access to container names containing data + const std::vector<std::string>& inputs() const { + return m_input_sources; + } + /// Access to container names containing data + const std::string& input_segment() const { + return m_input_segment; + } + /// Access to container names containing data + int input_mask() const { + return m_input_mask; + } /// Callback to read event input virtual void execute(DigiContext& context) const override; }; diff --git a/DDDigi/include/DDDigi/DigiROOTInput.h b/DDDigi/include/DDDigi/DigiROOTInput.h index 6107693fd9789a0a360c1cac913ef34adc943488..ee5d6fcede70681a784fe3078b89b43871cb025e 100644 --- a/DDDigi/include/DDDigi/DigiROOTInput.h +++ b/DDDigi/include/DDDigi/DigiROOTInput.h @@ -19,6 +19,9 @@ // C/C++ include files #include <memory> +/// Forward declarations +class TBranch; + /// Namespace for the AIDA detector description toolkit namespace dd4hep { @@ -37,18 +40,24 @@ namespace dd4hep { */ class DigiROOTInput : public DigiInputAction { - protected: + public: /// Helper classes class internals_t; + class inputsource_t; + + class work_t { + public: + DataSegment& input_segment; + Key input_key; + TBranch& branch; + TClass& cl; + }; + + protected: /// Property: Name of the tree to connect to std::string m_tree_name { }; /// Property: Container names to be loaded std::vector<std::string> m_containers { }; - /// Property: Segment name to place input data default: inputs - std::string m_location { }; - - /// Current input id - mutable int m_curr_input { 0 }; /// Connection parameters to the "current" input source mutable std::unique_ptr<internals_t> imp; @@ -64,8 +73,20 @@ namespace dd4hep { DigiROOTInput(const DigiKernel& kernel, const std::string& nam); /// Default destructor virtual ~DigiROOTInput(); + + /// Access to tree name containing data + const std::string& tree_name() const { + return m_tree_name; + } + /// Access to container names containing data + const std::vector<std::string>& container_names() const { + return m_containers; + } + /// Callback to read event input virtual void execute(DigiContext& context) const override; + /// Callback to handle single branch + virtual void operator()(DigiContext& context, work_t& work) const = 0; }; } // End namespace digi diff --git a/DDDigi/plugins/Components.cpp b/DDDigi/plugins/Components.cpp index b9ef72ec906249a0c1a86537b408798849161761..3d7c772e1c274794e67cee9a8de1ba5b831cf77b 100644 --- a/DDDigi/plugins/Components.cpp +++ b/DDDigi/plugins/Components.cpp @@ -24,9 +24,6 @@ #include <DDDigi/DigiInputAction.h> DECLARE_DIGIACTION_NS(dd4hep::digi,DigiInputAction) -#include <DDDigi/DigiROOTInput.h> -DECLARE_DIGIACTION_NS(dd4hep::digi,DigiROOTInput) - #include <DDDigi/DigiOutputAction.h> DECLARE_DIGIACTION_NS(dd4hep::digi,DigiOutputAction) diff --git a/DDDigi/plugins/DigiDepositSmearEnergy.cpp b/DDDigi/plugins/DigiDepositSmearEnergy.cpp new file mode 100644 index 0000000000000000000000000000000000000000..b76f7453041d2c1f24bda5e9e2362c8ca374d751 --- /dev/null +++ b/DDDigi/plugins/DigiDepositSmearEnergy.cpp @@ -0,0 +1,215 @@ +//========================================================================== +// 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/DigiContainerProcessor.h> +#include <DD4hep/InstanceCount.h> +#include <DD4hep/DD4hepUnits.h> + +/// C/C++ include files +#include <limits> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + /// Actor to smear geographically energy deposits + /** + * Intrinsic fluctuation: sigma(E)/E = const(fluctuation) / sqrt(E) + * sigma(E) = const(fluctuation) * sqrt(E) + * Sampling fluctuation: sigma(E)/E = const(sampling) * sqrt(deltaE/E) + * sigma(E) = const(sampling) * sqrt(deltaE*E) + * deltaE = energy loss to create single electron/ion + * pair in sampling layer + * Instrumetation: sigma(E)/E = const(instrumentaion) / E + * sigma(E) = const(instrumentaion) + * Systematics: sigma(E)/E = const(systematics) + * sigma(E) = const(systematics) * E + * + * electron conversion + * statictical fluctuation: poissonian change of converted e-ion pairs + * + * Electromagnetic showers: + * Shower counters: + * Fluctuation resolution: sigma(E)/E ~ 0.005 / sqrt(E) + * + * Sampling calorimeters: + * Sampling resolution: sigma(E)/E ~ 0.04 * sqrt(1000*dE/E) + * dE = energy loss of a single charged particle in one sampling layer + * + * Hadronic showers: + * Shower counters: + * Fluctuation resolution: sigma(E)/E ~ 0.45 / sqrt(E) + * With compensation for nuclear effects: sigma(E)/E ~ 0.25 / sqrt(E) + * Sampling resolution: sigma(E)/E ~ 0.09 * sqrt(1000*dE/E) + * dE = energy loss of a single charged particle in one sampling layer + * + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiDepositSmearEnergy : public DigiContainerProcessor { + protected: + /// Property: Energy cutoff. No hits will be merged with a deposit smaller + double m_deposit_cutoff { std::numeric_limits<double>::epsilon() }; + /// Property: Intrinsic energy resolution constant (gaussian ~ std::sqrt(energy)) + double m_intrinsic_fluctuation { 0e0 }; + /// Property: Systematic energy resolution constant (gaussian ~deposit energy) + double m_systematic_resolution { 0e0 }; + /// Property: Instrumentation energy resolution constant (gaussian(CONSTANT)) + double m_instrumentation_resolution { 0e0 }; + /// Property: Ionization constant to create electron ION pair in absorber + 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 drop killed entries from the deposits + bool m_drop_killed { false }; + /// Property: Flag to update existing container in-place or create a new container + bool m_update_in_place { true }; + + public: + /// Standard constructor + DigiDepositSmearEnergy(const DigiKernel& krnl, const std::string& nam) + : DigiContainerProcessor(krnl, nam) + { + declareProperty("drop_killed", m_drop_killed = false); + declareProperty("deposit_cutoff", m_deposit_cutoff); + declareProperty("intrinsic_fluctuation", m_intrinsic_fluctuation = 0e0); + declareProperty("instrumentation_resolution", m_instrumentation_resolution = 0e0); + 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("update_in_place", m_update_in_place = true); + InstanceCount::increment(this); + } + + /// Default destructor + virtual ~DigiDepositSmearEnergy() { + InstanceCount::decrement(this); + } + + template <typename T> std::size_t drop_killed_entries(T& cont) const { + //CellID last_cell = 0UL; + std::size_t killed = 0; + for( auto iter = cont.begin(); iter != cont.end(); ++iter ) { + auto& dep = iter->second; + if ( dep.flag&EnergyDeposit::KILLED ) { + cont.remove(iter); + iter = cont.begin(); // (last_cell != 0UL) ? cont.find(last_cell) : cont.begin(); + ++killed; + continue; + } + //last_cell = iter->first; + } + return killed; + } + + /// Create deposit mapping with updates on same cellIDs + template <typename T> void smear_deposit_energy(DigiContext& context, T& cont, work_t& work) const { + auto& random = context.randomGenerator(); + T output_cont(cont.name, work.output.mask); + std::size_t killed = 0UL; + std::size_t updated = 0UL; + std::size_t created = 0UL; + + if ( dd4hep::isActivePrintLevel(outputLevel()) ) { + print("%s+++ Smearing container: %s %6ld entries", + context.event->id(), cont.name.c_str(), cont.size()); + } + for( auto& dep : cont ) { + CellID cell = dep.first; + const EnergyDeposit& depo = dep.second; + double deposit = depo.deposit; + if ( deposit >= m_deposit_cutoff ) { + double delta_E = 0e0; + double energy = deposit / dd4hep::GeV; // E in units of GeV + double sigma_E_systematic = m_systematic_resolution * energy; + double sigma_E_intrin_fluct = m_intrinsic_fluctuation * std::sqrt(energy); + double sigma_E_instrument = m_instrumentation_resolution / dd4hep::GeV; + double delta_ion = 0e0, num_pairs = 0e0; + constexpr static double eps = std::numeric_limits<double>::epsilon(); + if ( sigma_E_systematic > eps ) { + delta_E += sigma_E_systematic * random.gaussian(0e0, 1e0); + } + if ( sigma_E_intrin_fluct > eps ) { + delta_E += sigma_E_intrin_fluct * random.gaussian(0e0, 1e0); + } + if ( sigma_E_instrument > eps ) { + delta_E += sigma_E_instrument * random.gaussian(0e0, 1e0); + } + if ( m_ionization_fluctuation ) { + num_pairs = energy / (m_pair_ionization_energy/dd4hep::GeV); + delta_ion = energy * (random.poisson(num_pairs)/num_pairs); + delta_E += delta_ion; + } + if ( dd4hep::isActivePrintLevel(outputLevel()) ) { + print("+++ %016lX [GeV] E:%9.2e [%9.2e %9.2e] intrin_fluct:%9.2e systematic:%9.2e instrument:%9.2e ioni:%9.2e/%.0f %s", + cell, energy, deposit/dd4hep::GeV, delta_E, + sigma_E_intrin_fluct, sigma_E_systematic, sigma_E_instrument, delta_ion, num_pairs, + deposit < m_deposit_cutoff ? "KILLED" : ""); + } + /// delta_E is in GeV + if ( deposit > m_deposit_cutoff ) { + if ( m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.deposit = deposit + (delta_E * dd4hep::GeV); + d.flag |= EnergyDeposit::ENERGY_SMEARED; + ++updated; + } + else { + EnergyDeposit d(depo); + d.deposit = deposit + (delta_E * dd4hep::GeV); + d.flag |= EnergyDeposit::ENERGY_SMEARED; + output_cont.emplace(cell, EnergyDeposit(depo)); + ++created; + } + continue; + } + } + else if ( !m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.flag |= EnergyDeposit::KILLED; + ++killed; + } + } + + if ( !m_update_in_place ) { + work.output.data.put(output_cont.key, std::move(output_cont)); + } + + killed = (m_drop_killed && m_update_in_place && killed > 0) + ? drop_killed_entries(cont) + : killed; + + info("%s+++ %-32s Smearing: created %6ld updated %6ld killed %6ld entries from mask: %04X", + context.event->id(), cont.name.c_str(), created, updated, killed, cont.key.mask()); + } + + /// Main functional callback + virtual void execute(DigiContext& context, work_t& work) const override final { + if ( auto* v = work.get_input<DepositVector>() ) + smear_deposit_energy(context, *v, work); + else if ( auto* m = work.get_input<DepositMapping>() ) + smear_deposit_energy(context, *m, work); + else + except("Request to handle unknown data type: %s", work.input_type_name().c_str()); + } + }; + } // End namespace digi +} // End namespace dd4hep + +#include <DDDigi/DigiFactories.h> +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearEnergy) diff --git a/DDDigi/plugins/DigiDepositSmearPositionResolution.cpp b/DDDigi/plugins/DigiDepositSmearPositionResolution.cpp new file mode 100644 index 0000000000000000000000000000000000000000..960d65c268af41edf056f43f2627a1394b47530c --- /dev/null +++ b/DDDigi/plugins/DigiDepositSmearPositionResolution.cpp @@ -0,0 +1,141 @@ +//========================================================================== +// 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/DigiContainerProcessor.h> +#include <DDDigi/DigiKernel.h> + +#include <DD4hep/InstanceCount.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/Detector.h> +#include <DD4hep/VolumeManager.h> + +/// C/C++ include files +#include <limits> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + /// Actor to smear geographically energy deposits + /** + * + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiDepositSmearPositionResolution : public DigiContainerProcessor { + protected: + /// Property: Energy cutoff. No hits will be merged with a deposit smaller + double m_deposit_cutoff { detail::numeric_epsilon }; + /// Property: Resolution in loacl sensor's U coordinates along (1, 0, 0) + double m_resolution_u { 0e0 * dd4hep::mm }; + /// Property: Resolution in loacl sensor's V coordinates along (0, 1, 0) + double m_resolution_v { 0e0 * dd4hep::mm }; + /// Property: Flag to update existing container in-place or create a new container + bool m_update_in_place { true }; + + public: + /// Standard constructor + DigiDepositSmearPositionResolution(const DigiKernel& krnl, const std::string& nam) + : DigiContainerProcessor(krnl, nam) + { + declareProperty("deposit_cutoff", m_deposit_cutoff); + declareProperty("update_in_place", m_update_in_place = true); + declareProperty("resolution_u", m_resolution_u); + declareProperty("resolution_v", m_resolution_v); + InstanceCount::increment(this); + } + + /// Default destructor + virtual ~DigiDepositSmearPositionResolution() { + InstanceCount::decrement(this); + } + + /// Create deposit mapping with updates on same cellIDs + template <typename T> void smear_deposit_position(DigiContext& context, T& cont, work_t& work) const { + T output_cont(cont.name, work.output.mask); + auto& random = context.randomGenerator(); + std::size_t killed = 0UL; + std::size_t updated = 0UL; + std::size_t created = 0UL; + + if ( dd4hep::isActivePrintLevel(outputLevel()) ) { + print("%s+++ Smearing container: %s %6ld entries", + context.event->id(), cont.name.c_str(), cont.size()); + } + VolumeManager volMgr = m_kernel.detectorDescription().volumeManager(); + for( auto& dep : cont ) { + CellID cell = dep.first; + const EnergyDeposit& depo = dep.second; + double deposit = depo.deposit; + if ( deposit >= m_deposit_cutoff ) { + auto* ctxt = volMgr.lookupContext(cell); + Position local_pos = ctxt->worldToLocal(depo.position); + double delta_u = m_resolution_u * random.gaussian(); + double delta_v = m_resolution_v * random.gaussian(); + Position delta_pos(delta_u, delta_v, 0e0); + Position oldpos = depo.position; + Position newpos = ctxt->localToWorld(local_pos + delta_pos); + + delta_pos = newpos - oldpos; + info("+++ %016lX Smeared position [%9.2e %9.2e %9.2e] by [%9.2e %9.2e %9.2e] to [%9.2e %9.2e %9.2e] [mm]", + oldpos.X(), oldpos.Y(), oldpos.Z(), delta_pos.X(), delta_pos.Y(), delta_pos.Z(), + newpos.X(), newpos.Y(), newpos.Z()); + + /// Update / create deposit with smeared position + if ( m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.position = newpos; + d.flag |= EnergyDeposit::POSITION_SMEARED; + ++updated; + } + else { + EnergyDeposit d(depo); + d.position = newpos; + d.flag |= EnergyDeposit::POSITION_SMEARED; + output_cont.emplace(cell, std::move(d)); + ++created; + } + } + else if ( !m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.flag |= EnergyDeposit::KILLED; + ++killed; + } + } + + if ( !m_update_in_place ) { + work.output.data.put(output_cont.key, std::move(output_cont)); + } + info("%s+++ %-32s Smearing: created %6ld updated %6ld killed %6ld entries from mask: %04X", + context.event->id(), cont.name.c_str(), created, updated, killed, cont.key.mask()); + } + + /// Main functional callback + virtual void execute(DigiContext& context, work_t& work) const override final { + if ( auto* v = work.get_input<DepositVector>() ) + smear_deposit_position(context, *v, work); + else if ( auto* m = work.get_input<DepositMapping>() ) + smear_deposit_position(context, *m, work); + else + except("Request to handle unknown data type: %s", work.input_type_name().c_str()); + } + }; + } // End namespace digi +} // End namespace dd4hep + +#include <DDDigi/DigiFactories.h> +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearPositionResolution) diff --git a/DDDigi/plugins/DigiDepositSmearPositionTrack.cpp b/DDDigi/plugins/DigiDepositSmearPositionTrack.cpp new file mode 100644 index 0000000000000000000000000000000000000000..652ae28da8a61ec2c6bf2b1bb76a412ab34e0c51 --- /dev/null +++ b/DDDigi/plugins/DigiDepositSmearPositionTrack.cpp @@ -0,0 +1,153 @@ +//========================================================================== +// 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/DigiContainerProcessor.h> +#include <DDDigi/DigiKernel.h> + +#include <DD4hep/InstanceCount.h> +#include <DD4hep/DD4hepUnits.h> +#include <DD4hep/Detector.h> +#include <DD4hep/VolumeManager.h> + +/// C/C++ include files +#include <limits> + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + /// Namespace for the Digitization part of the AIDA detector description toolkit + namespace digi { + + /// Actor to smear geographically energy deposits + /** + * + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ + class DigiDepositSmearPositionTrack : public DigiContainerProcessor { + protected: + /// Property: Energy cutoff. No hits will be merged with a deposit smaller + double m_deposit_cutoff { detail::numeric_epsilon }; + /// Property: Resolution in loacl sensor's U coordinates along (1, 0, 0) + double m_resolution_u { 0e0 * dd4hep::mm }; + /// Property: Resolution in loacl sensor's V coordinates along (0, 1, 0) + double m_resolution_v { 0e0 * dd4hep::mm }; + /// Property: Flag to update existing container in-place or create a new container + bool m_update_in_place { true }; + + public: + /// Standard constructor + DigiDepositSmearPositionTrack(const DigiKernel& krnl, const std::string& nam) + : DigiContainerProcessor(krnl, nam) + { + declareProperty("deposit_cutoff", m_deposit_cutoff); + declareProperty("update_in_place", m_update_in_place = true); + declareProperty("resolution_u", m_resolution_u); + declareProperty("resolution_v", m_resolution_v); + InstanceCount::increment(this); + } + + /// Default destructor + virtual ~DigiDepositSmearPositionTrack() { + InstanceCount::decrement(this); + } + + /// Create deposit mapping with updates on same cellIDs + template <typename T> void smear_deposit_position(DigiContext& context, T& cont, work_t& work) const { + constexpr double eps = detail::numeric_epsilon; + T output_cont(cont.name, work.output.mask); + auto& random = context.randomGenerator(); + const auto& ev = *(context.event); + std::size_t killed = 0UL; + std::size_t updated = 0UL; + std::size_t created = 0UL; + + if ( dd4hep::isActivePrintLevel(outputLevel()) ) { + print("%s+++ Smearing container: %s %6ld entries", + context.event->id(), cont.name.c_str(), cont.size()); + } + VolumeManager volMgr = m_kernel.detectorDescription().volumeManager(); + for( auto& dep : cont ) { + CellID cell = dep.first; + const EnergyDeposit& depo = dep.second; + double deposit = depo.deposit; + if ( deposit >= m_deposit_cutoff ) { + auto* ctxt = volMgr.lookupContext(cell); + Direction part_momentum = depo.history.average_particle_momentum(ev); + Position local_pos = ctxt->worldToLocal(depo.position); + Position local_dir = ctxt->worldToLocal(part_momentum).unit(); + double cos_u = local_dir.Dot(Position(1,0,0)); + double sin_u = std::sqrt(1e0 - cos_u*cos_u); + double tan_u = sin_u/(std::abs(cos_u)>eps ? cos_u : eps); + double delta_u = tan_u * m_resolution_u * random.gaussian(); + + double cos_v = local_dir.Dot(Position(0,1,0)); + double sin_v = std::sqrt(1e0 - cos_v*cos_v); + double tan_v = sin_v/(std::abs(cos_v)>eps ? cos_v : eps); + double delta_v = tan_v * m_resolution_v * random.gaussian(); + + Position delta_pos(delta_u, delta_v, 0e0); + Position oldpos = depo.position; + Position newpos = ctxt->localToWorld(local_pos + delta_pos); + + delta_pos = newpos - oldpos; + info("+++ %016lX Smeared position [%9.2e %9.2e %9.2e] by [%9.2e %9.2e %9.2e] to [%9.2e %9.2e %9.2e] [mm]", + oldpos.X(), oldpos.Y(), oldpos.Z(), delta_pos.X(), delta_pos.Y(), delta_pos.Z(), + newpos.X(), newpos.Y(), newpos.Z()); + + /// Update / create deposit with smeared position + if ( m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.position = newpos; + d.flag |= EnergyDeposit::POSITION_SMEARED; + ++updated; + } + else { + EnergyDeposit d(depo); + d.position = newpos; + d.flag |= EnergyDeposit::POSITION_SMEARED; + output_cont.emplace(cell, std::move(d)); + ++created; + } + } + else if ( !m_update_in_place ) { + EnergyDeposit& d = dep.second; + d.flag |= EnergyDeposit::KILLED; + ++killed; + } + } + + if ( !m_update_in_place ) { + work.output.data.put(output_cont.key, std::move(output_cont)); + } + info("%s+++ %-32s Smearing: created %6ld updated %6ld killed %6ld entries from mask: %04X", + context.event->id(), cont.name.c_str(), created, updated, killed, cont.key.mask()); + } + + /// Main functional callback + virtual void execute(DigiContext& context, work_t& work) const override final { + if ( auto* v = work.get_input<DepositVector>() ) + smear_deposit_position(context, *v, work); + else if ( auto* m = work.get_input<DepositMapping>() ) + smear_deposit_position(context, *m, work); + else + except("Request to handle unknown data type: %s", work.input_type_name().c_str()); + } + }; + } // End namespace digi +} // End namespace dd4hep + +#include <DDDigi/DigiFactories.h> +DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearPositionTrack) diff --git a/DDDigi/plugins/DigiDepositSmearing.cpp b/DDDigi/plugins/DigiDepositSmearing.cpp deleted file mode 100644 index da28a26234e2de48b53edd7a0a39627ff6ae3a79..0000000000000000000000000000000000000000 --- a/DDDigi/plugins/DigiDepositSmearing.cpp +++ /dev/null @@ -1,87 +0,0 @@ -//========================================================================== -// 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/DigiContainerProcessor.h> - -/// C/C++ include files -#include <limits> - -/// Namespace for the AIDA detector description toolkit -namespace dd4hep { - /// Namespace for the Digitization part of the AIDA detector description toolkit - namespace digi { - - /// Actor to smear geographically energy deposits - /** - * - * \author M.Frank - * \version 1.0 - * \ingroup DD4HEP_DIGITIZATION - */ - class DigiDepositSmearing : public DigiContainerProcessor { - protected: - /// Property: Energy cutoff. No hits will be merged with a deposit smaller - double m_cutoff { std::numeric_limits<double>::epsilon() }; - /// Property: register empty containers - bool m_register_empty_containers { true }; - - public: - /// Standard constructor - DigiDepositSmearing(const DigiKernel& krnl, const std::string& nam) - : DigiContainerProcessor(krnl, nam) - { - declareProperty("deposit_cutoff", m_cutoff); - declareProperty("m_register_empty_containers", m_register_empty_containers); - } - - /// Create deposit mapping with updates on same cellIDs - template <typename T> void create_deposits(const T& cont, work_t& work) const { - Key key(cont.name, work.output.mask); - DepositMapping m(cont.name, work.output.mask); - std::size_t dropped = 0UL, updated = 0UL, added = 0UL; - for( const auto& dep : cont ) { - const EnergyDeposit& depo = dep.second; - if ( depo.deposit >= m_cutoff ) { - CellID cell = dep.first; - auto iter = m.data.find(cell); - if ( iter == m.data.end() ) - m.data.emplace(dep.first, depo), ++added; - else - iter->second.update_deposit_weighted(depo), ++updated; - continue; - } - ++dropped; - } - if ( m_register_empty_containers ) { - work.output.data.put(m.key, std::move(m)); - } - info("+++ %-32s added %6ld updated %6ld dropped %6ld entries (now: %6ld) from mask: %04X to mask: %04X", - cont.name.c_str(), added, updated, dropped, m.size(), cont.key.mask(), m.key.mask()); - } - - /// Main functional callback - virtual void execute(DigiContext&, work_t& work) const override final { - if ( const auto* v = work.get_input<DepositVector>() ) - create_deposits(*v, work); - else if ( const auto* m = work.get_input<DepositMapping>() ) - create_deposits(*m, work); - else - except("Request to handle unknown data type: %s", work.input_type_name().c_str()); - } - }; - } // End namespace digi -} // End namespace dd4hep - -#include <DDDigi/DigiFactories.h> -DECLARE_DIGIACTION_NS(dd4hep::digi,DigiDepositSmearing) diff --git a/DDDigi/plugins/DigiDepositWeightedPosition.cpp b/DDDigi/plugins/DigiDepositWeightedPosition.cpp index 0c2a2f786281d18eda0f3789cde9a29199c76f2c..7de2ee4e716028c082667cf8db72157d8929ef91 100644 --- a/DDDigi/plugins/DigiDepositWeightedPosition.cpp +++ b/DDDigi/plugins/DigiDepositWeightedPosition.cpp @@ -13,6 +13,7 @@ // Framework include files #include <DDDigi/DigiContainerProcessor.h> +#include <DD4hep/InstanceCount.h> /// C/C++ include files #include <limits> @@ -37,8 +38,6 @@ namespace dd4hep { protected: /// Property: Energy cutoff. No hits will be merged with a deposit smaller double m_cutoff { -std::numeric_limits<double>::epsilon() }; - /// Property: register empty containers - bool m_register_empty_containers { true }; public: /// Standard constructor @@ -46,7 +45,12 @@ namespace dd4hep { : DigiContainerProcessor(krnl, nam) { declareProperty("deposit_cutoff", m_cutoff); - declareProperty("m_register_empty_containers", m_register_empty_containers); + InstanceCount::increment(this); + } + + /// Default destructor + virtual ~DigiDepositWeightedPosition() { + InstanceCount::decrement(this); } /// Create deposit mapping with updates on same cellIDs @@ -67,11 +71,9 @@ namespace dd4hep { } ++dropped; } - if ( m_register_empty_containers ) { - work.output.data.put(m.key, std::move(m)); - } info("%s+++ %-32s added %6ld updated %6ld dropped %6ld entries (now: %6ld) from mask: %04X to mask: %04X", tag, cont.name.c_str(), added, updated, dropped, m.size(), cont.key.mask(), m.key.mask()); + work.output.data.put(m.key, std::move(m)); } /// Main functional callback diff --git a/DDDigi/plugins/DigiSegmentDepositPrint.cpp b/DDDigi/plugins/DigiSegmentDepositPrint.cpp index a76fa9f06c81d482ee60dcda043e13469300d509..a97281da960e58fdd975978c241aed403b31ec0a 100644 --- a/DDDigi/plugins/DigiSegmentDepositPrint.cpp +++ b/DDDigi/plugins/DigiSegmentDepositPrint.cpp @@ -33,11 +33,16 @@ namespace dd4hep { DigiSegmentDepositPrint(const DigiKernel& kernel, const std::string& nam) : DigiSegmentProcessor(kernel, nam) {} - void print_deposit(const char* format, CellID cell, const EnergyDeposit& depo) const { - info(format, segment.split_id(cell), cell, - depo.history.hits.size(), - depo.history.particles.size(), - depo.deposit); + /// Single container printout + template <typename T> void print_deposits(const char* fmt, const T& cont) const { + for(const auto& dep : cont) { + if( this->segment.matches(dep.first) ) { + info(fmt, segment.split_id(dep.first), dep.first, + dep.second.history.hits.size(), + dep.second.history.particles.size(), + dep.second.deposit); + } + } } /// Main functional callback virtual void execute(DigiContext& context, work_t& work) const override final { @@ -45,14 +50,10 @@ namespace dd4hep { ::snprintf(format, sizeof(format), "%s[%s] %s-id: %%d [processor:%d] Cell: %%016lX mask: %016lX hist:%%4ld hits %%4ld parts. entries deposit: %%f", context.event->id(), segment.idspec.name(), segment.cname(), segment.id, segment.split_mask); - auto call = [this, format](const std::pair<CellID,EnergyDeposit>& d) { - if( this->segment.matches(d.first) ) - this->print_deposit(format, d.first, d.second); - }; if ( const auto* m = work.get_input<DepositMapping>() ) - std::for_each(m->begin(), m->end(), call); + print_deposits(format, *m); else if ( const auto* v = work.get_input<DepositVector>() ) - std::for_each(v->begin(), v->end(), call); + print_deposits(format, *v); else error("+++ Request to dump an invalid container %s", Key::key_name(work.input.key).c_str()); } diff --git a/DDDigi/plugins/DigiSimpleADCResponse.cpp b/DDDigi/plugins/DigiSimpleADCResponse.cpp index b540f1d9a0c7479c713341e99f60840fc9b2b969..540901261a7b54192a8e02f325c667e9e4afc99f 100644 --- a/DDDigi/plugins/DigiSimpleADCResponse.cpp +++ b/DDDigi/plugins/DigiSimpleADCResponse.cpp @@ -29,6 +29,7 @@ namespace dd4hep { double m_use_segmentation { false }; double m_signal_cutoff { std::numeric_limits<double>::epsilon() }; double m_signal_saturation { std::numeric_limits<double>::max() }; + double m_adc_offset { 0e0 }; std::size_t m_adc_resolution { 1024 }; public: @@ -47,15 +48,17 @@ namespace dd4hep { template <typename T, typename P> void create_adc_counts(const char* tag, const T& input, work_t& work, const P& predicate) const { std::string postfix = m_use_segmentation ? "."+segment.identifier() : std::string(); - std::string response_name = input.name + postfix + m_response_postfix; std::string history_name = input.name + postfix + m_history_postfix; - DetectorResponse response(response_name, work.output.mask); + std::string response_name = input.name + postfix + m_response_postfix; DetectorHistory history (history_name, work.output.mask); + DetectorResponse response(response_name, work.output.mask); for( const auto& dep : input ) { if ( predicate(dep) ) { CellID cell = dep.first; const auto& depo = dep.second; - ADCValue::value_t adc_count = std::round((depo.deposit * m_adc_resolution) / m_signal_saturation); + double offset_ene = depo.deposit-m_adc_offset; + ADCValue::value_t adc_count = std::round(((offset_ene) * m_adc_resolution) / m_signal_saturation); + adc_count = std::min(adc_count, ADCValue::value_t(m_adc_resolution)); response.emplace(cell, {adc_count, ADCValue::address_t(cell)}); history.insert(cell, depo.history); } diff --git a/DDDigi/src/DigiContainerCombine.cpp b/DDDigi/src/DigiContainerCombine.cpp index ce2aec85d7f96a1901810f72035d030522baa9ec..e6dc17af81469f4f5bf573b701db90d3dcce5649 100644 --- a/DDDigi/src/DigiContainerCombine.cpp +++ b/DDDigi/src/DigiContainerCombine.cpp @@ -84,14 +84,17 @@ public: /// Specialized deposit merger: implicitly assume identical item types are mapped sequentially template<typename OUT, typename IN> void merge_depos(OUT& output, IN& input, int thr) { std::size_t cnt = 0; + std::string nam = input.name; + auto mask = input.key.mask(); if ( combine->m_erase_combined ) cnt = output.merge(std::move(input)); else cnt = output.insert(input); - combine->debug(this->format, thr, input.name.c_str(), input.key.mask(), cnt, "deposits"); + combine->debug(this->format, thr, nam.c_str(), mask, cnt, "deposits"); this->cnt_depos += cnt; this->cnt_conts++; } + /// Merge history records: implicitly assume identical item types are mapped sequentially void merge_hist(const std::string& nam, size_t start, int thr) { std::size_t cnt; @@ -100,8 +103,9 @@ public: for( std::size_t j=start; j < keys.size(); ++j ) { if ( keys[j].item() == key.item() ) { DetectorHistory* next = std::any_cast<DetectorHistory>(work[j]); + std::string next_name = next->name; cnt = (combine->m_erase_combined) ? out.merge(std::move(*next)) : out.insert(*next); - combine->debug(format, thr, next->name.c_str(), keys[j].mask(), cnt, "histories"); + combine->debug(format, thr, next_name.c_str(), keys[j].mask(), cnt, "histories"); used_keys_insert(keys[j]); cnt_hist += cnt; cnt_conts++; @@ -118,8 +122,9 @@ public: for( std::size_t j=start; j < keys.size(); ++j ) { if ( keys[j].item() == key.item() ) { DetectorResponse* next = std::any_cast<DetectorResponse>(work[j]); + std::string next_name = next->name; cnt = (combine->m_erase_combined) ? out.merge(std::move(*next)) : out.insert(*next); - combine->debug(format, thr, next->name.c_str(), keys[j].mask(), cnt, "histories"); + combine->debug(format, thr, next_name.c_str(), keys[j].mask(), cnt, "histories"); used_keys_insert(keys[j]); cnt_response += cnt; cnt_conts++; @@ -136,8 +141,9 @@ public: for( std::size_t j=start; j < keys.size(); ++j ) { if ( keys[j].item() == key.item() ) { ParticleMapping* next = std::any_cast<ParticleMapping>(work[j]); + std::string next_name = next->name; cnt = (combine->m_erase_combined) ? out.merge(std::move(*next)) : out.insert(*next); - combine->debug(format, thr, next->name.c_str(), keys[j].mask(), cnt, "particles"); + combine->debug(format, thr, next_name.c_str(), keys[j].mask(), cnt, "particles"); used_keys_insert(keys[j]); cnt_parts += cnt; cnt_conts++; diff --git a/DDDigi/src/DigiContainerProcessor.cpp b/DDDigi/src/DigiContainerProcessor.cpp index 758f75162a8246d43572294c78d691488810ac60..df186c91bd84893aef2a624f61a754a41f016f90 100644 --- a/DDDigi/src/DigiContainerProcessor.cpp +++ b/DDDigi/src/DigiContainerProcessor.cpp @@ -159,6 +159,7 @@ void dd4hep::digi::DigiContainerSequenceAction::initialize() { worker_t* w = new worker_t(ent.second, m_registered_workers.size()); m_registered_workers.emplace(ent.first, w); m_workers.insert(w); + ent.second->release(); } } @@ -294,7 +295,7 @@ void DigiMultiContainerProcessor::initialize() { m_workers.insert(w); action_workers[action] = w; } - for( auto key : worker_keys ) + for( const auto& key : worker_keys ) m_worker_map[key.item()].push_back(w); } /// Add some printout about the configuration diff --git a/DDDigi/src/DigiData.cpp b/DDDigi/src/DigiData.cpp index 1eb0245ace5eddb12fb27544c6039ba806c769a9..9f40c8773052702cc4d847c02103fded09895d9a 100644 --- a/DDDigi/src/DigiData.cpp +++ b/DDDigi/src/DigiData.cpp @@ -82,7 +82,7 @@ std::string Key::key_name(const Key& k) { return "UNKNOWN"; } -const Particle& dd4hep::digi::get_history_particle(DigiEvent& event, Key history_key) { +const Particle& dd4hep::digi::get_history_particle(const DigiEvent& event, Key history_key) { static Key item_key("MCParticles",0x0); Key key; const auto& segment = event.get_segment(history_key.segment()); @@ -93,7 +93,7 @@ const Particle& dd4hep::digi::get_history_particle(DigiEvent& event, Key history return particles.get(history_key); } -const EnergyDeposit& dd4hep::digi::get_history_deposit(DigiEvent& event, Key::itemkey_type container_item, Key history_key) { +const EnergyDeposit& dd4hep::digi::get_history_deposit(const DigiEvent& event, Key::itemkey_type container_item, Key history_key) { Key key; const auto& segment = event.get_segment(history_key.segment()); key.set_segment(history_key.segment()); @@ -117,7 +117,7 @@ std::pair<std::size_t,std::size_t> History::drop() { return ret; } -const Particle& History::hist_entry_t::get_particle(DigiEvent& event) const { +const Particle& History::hist_entry_t::get_particle(const DigiEvent& event) const { static Key item_key("MCParticles",0x0); Key key(this->source); const auto& segment = event.get_segment(this->source.segment()); @@ -125,13 +125,28 @@ const Particle& History::hist_entry_t::get_particle(DigiEvent& event) const { return particles.get(this->source); } -const EnergyDeposit& History::hist_entry_t::get_deposit(DigiEvent& event, Key::itemkey_type container_item) const { +const EnergyDeposit& History::hist_entry_t::get_deposit(const DigiEvent& event, Key::itemkey_type container_item) const { Key key(this->source); const auto& segment = event.get_segment(this->source.segment()); const auto& hits = segment.get<DepositVector>(key.set_item(container_item)); return hits.at(this->source.item()); } +/// Retrieve the weighted momentum of all contributing particles +Direction History::average_particle_momentum(const DigiEvent& event) const { + Direction momentum; + if ( !particles.empty() ) { + double weight = 0e0; + for( const auto& p : particles ) { + const Particle& par = p.get_particle(event); + momentum += par.momentum * p.weight; + weight += p.weight; + } + momentum /= std::max(weight, detail::numeric_epsilon); + } + return momentum; +} + /// Update the deposit using deposit weighting void EnergyDeposit::update_deposit_weighted(const EnergyDeposit& upda) { double sum = deposit + upda.deposit; @@ -157,7 +172,8 @@ void EnergyDeposit::update_deposit_weighted(EnergyDeposit&& upda) { /// Merge new deposit map onto existing map std::size_t DepositVector::merge(DepositVector&& updates) { std::size_t update_size = updates.size(); - data.reserve(data.size()+updates.size()); + std::size_t newlen = std::max(2*data.size(), data.size()+updates.size()); + data.reserve(newlen); for( auto& c : updates ) { data.emplace_back(c); } @@ -167,7 +183,8 @@ std::size_t DepositVector::merge(DepositVector&& updates) { /// Merge new deposit map onto existing map (keep inputs) std::size_t DepositVector::merge(DepositMapping&& updates) { std::size_t update_size = updates.size(); - data.reserve(data.size()+updates.size()); + std::size_t newlen = std::max(2*data.size(), data.size()+updates.size()); + data.reserve(newlen); for( const auto& c : updates ) { data.emplace_back(c); } @@ -177,7 +194,8 @@ std::size_t DepositVector::merge(DepositMapping&& updates) { /// Merge new deposit map onto existing map (keep inputs) std::size_t DepositVector::insert(const DepositVector& updates) { std::size_t update_size = updates.size(); - data.reserve(data.size()+updates.size()); + std::size_t newlen = std::max(2*data.size(), data.size()+updates.size()); + data.reserve(newlen); for( const auto& c : updates ) { data.emplace_back(c); } @@ -187,7 +205,8 @@ std::size_t DepositVector::insert(const DepositVector& updates) { /// Merge new deposit map onto existing map (keep inputs) std::size_t DepositVector::insert(const DepositMapping& updates) { std::size_t update_size = updates.size(); - data.reserve(data.size()+updates.size()); + std::size_t newlen = std::max(2*data.size(), data.size()+updates.size()); + data.reserve(newlen); for( const auto& c : updates ) { data.emplace_back(c); } @@ -210,6 +229,11 @@ const EnergyDeposit& DepositVector::at(std::size_t cell) const { return data.at(cell).second; } +/// Remove entry +void DepositVector::remove(iterator position) { + //data.erase(position); +} + /// Merge new deposit map onto existing map std::size_t DepositMapping::merge(DepositVector&& updates) { std::size_t update_size = updates.size(); @@ -258,6 +282,11 @@ std::size_t DepositMapping::insert(const DepositMapping& updates) { return update_size; } +/// Emplace entry +void DepositMapping::emplace(CellID cell, EnergyDeposit&& deposit) { + data.emplace(cell, std::move(deposit)); +} + /// Access energy deposit by key const EnergyDeposit& DepositMapping::get(CellID cell) const { auto iter = data.find(cell); @@ -267,6 +296,11 @@ const EnergyDeposit& DepositMapping::get(CellID cell) const { throw std::runtime_error("Failed to access deposit by CellID"); } +/// Remove entry +void DepositMapping::remove(iterator position) { + data.erase(position); +} + /// Merge new deposit map onto existing map std::size_t ParticleMapping::insert(const ParticleMapping& updates) { std::size_t update_size = updates.size(); @@ -353,10 +387,10 @@ DataSegment::DataSegment(std::mutex& l, Key::segment_type i) /// Remove data item from segment bool DataSegment::emplace(Key key, std::any&& item) { bool has_value = item.has_value(); -#if 0 - printout(INFO, "DataSegment", "PUT Key No.%4d: %-32s %016lX -> %04X %04X %08Xld %s", +#if DD4HEP_DDDIGI_DEBUG + printout(INFO, "DataSegment", "PUT Key No.%4d: %-32s %016lX -> %04X %04X %08Xld Value:%s %s", data.size(), Key::key_name(key).c_str(), key.value(), key.segment(), key.mask(), key.item(), - digiTypeName(item.type()).c_str()); + yes_no(has_value), digiTypeName(item.type()).c_str()); #endif std::lock_guard<std::mutex> l(lock); bool ret = data.emplace(key, std::move(item)).second; @@ -394,7 +428,7 @@ bool DataSegment::erase(Key key) { std::size_t DataSegment::erase(const std::vector<Key>& keys) { std::size_t count = 0; std::lock_guard<std::mutex> l(lock); - for(auto key : keys) { + for(const auto& key : keys) { auto iter = data.find(key); if ( iter != data.end() ) { data.erase(iter); @@ -510,6 +544,30 @@ DataSegment& DigiEvent::get_segment(const std::string& name) { throw std::runtime_error("Invalid segment name"); } +/// Retrieve data segment from the event structure +const DataSegment& DigiEvent::get_segment(const std::string& name) const { + switch(::toupper(name[0])) { + case 'I': + return *m_inputs; + case 'C': + return *m_counts; + case 'D': + switch(::toupper(name[1])) { + case 'E': + return *m_deposits; + default: + return *m_data; + } + break; + case 'O': + return *m_outputs; + default: + break; + } + except("DigiEvent", "Invalid segment name: %s", name.c_str()); + throw std::runtime_error("Invalid segment name"); +} + /// Retrieve data segment from the event structure by identifier DataSegment& DigiEvent::get_segment(Key::segment_type id) { switch(id) { @@ -530,3 +588,23 @@ DataSegment& DigiEvent::get_segment(Key::segment_type id) { throw std::runtime_error("Invalid segment name"); } +/// Retrieve data segment from the event structure by identifier +const DataSegment& DigiEvent::get_segment(Key::segment_type id) const { + switch(id) { + case 1: + return *m_inputs; + case 2: + return *m_counts; + case 3: + return *m_deposits; + case 4: + return *m_outputs; + case 0xAB: + return *m_data; + default: + break; + } + except("DigiEvent", "Invalid segment identifier: %d", id); + throw std::runtime_error("Invalid segment name"); +} + diff --git a/DDDigi/src/DigiInputAction.cpp b/DDDigi/src/DigiInputAction.cpp index 8b141e319fdb8f1f2d56fa3b6fb0d9b653ad1655..160803e90bce6c7371af1478d71f7e28812e6321 100644 --- a/DDDigi/src/DigiInputAction.cpp +++ b/DDDigi/src/DigiInputAction.cpp @@ -26,8 +26,10 @@ using namespace dd4hep::digi; DigiInputAction::DigiInputAction(const DigiKernel& kernel, const string& nam) : DigiEventAction(kernel, nam) { - declareProperty("input", m_inputs); - declareProperty("mask", m_mask); + declareProperty("input", m_input_sources); + declareProperty("segment", m_input_segment); + declareProperty("mask", m_input_mask); + declareProperty("rescan_inputs", m_input_rescan); InstanceCount::increment(this); } diff --git a/DDDigi/src/DigiKernel.cpp b/DDDigi/src/DigiKernel.cpp index 8afe11f7c2e4d131d00f3c9748485d1d8293820d..920dbd630f5015d2527ed85b41c1aa9b04b7cfe8 100644 --- a/DDDigi/src/DigiKernel.cpp +++ b/DDDigi/src/DigiKernel.cpp @@ -435,7 +435,7 @@ int DigiKernel::run() { internals->events_todo = internals->numEvents; printout(INFO,"DigiKernel","+++ Total number of events: %d",internals->numEvents); #ifdef DD4HEP_USE_TBB - if ( !internals->tbb_init && internals->num_threads>=0 ) { + if ( !internals->tbb_init && internals->num_threads > 0 ) { using ctrl_t = tbb::global_control; if ( 0 == internals->num_threads ) { internals->num_threads = ctrl_t::max_allowed_parallelism; @@ -443,22 +443,25 @@ int DigiKernel::run() { printout(INFO, "DigiKernel", "+++ Number of TBB threads to: %d",internals->num_threads); printout(INFO, "DigiKernel", "+++ Number of parallel events: %d",internals->maxEventsParallel); internals->tbb_init = std::make_unique<ctrl_t>(ctrl_t::max_allowed_parallelism,internals->num_threads+1); - if ( internals->maxEventsParallel > 1 ) { + if ( internals->maxEventsParallel >= 0 ) { int todo_evt = internals->events_todo; int num_proc = std::min(todo_evt,internals->maxEventsParallel); tbb::task_group main_group; for(int i=0; i < num_proc; ++i) main_group.run(Processor(*this)); main_group.wait(); - printout(DEBUG, "DigiKernel", "+++ All event processing threads Synchronized --- Done!"); } + printout(DEBUG, "DigiKernel", "+++ All event processing threads Synchronized --- Done!"); } + else #endif - while ( internals->events_todo > 0 && !internals->stop ) { - Processor proc(*this); - proc(); - ++internals->events_submitted; - } + { + while ( internals->events_todo > 0 && !internals->stop ) { + Processor proc(*this); + proc(); + ++internals->events_submitted; + } + } std::chrono::duration<double> duration = std::chrono::system_clock::now() - start; double sec = std::chrono::duration_cast<std::chrono::seconds>(duration).count(); printout(DEBUG, "DigiKernel", "+++ %d Events out of %d processed. " diff --git a/DDDigi/src/DigiROOTInput.cpp b/DDDigi/src/DigiROOTInput.cpp index 8d353559422ea3ed73516656de0584625f58ff7c..5a8382ee4398e641fcecb64b8f3b8eb67bb2e572 100644 --- a/DDDigi/src/DigiROOTInput.cpp +++ b/DDDigi/src/DigiROOTInput.cpp @@ -28,183 +28,190 @@ using namespace dd4hep::digi; -/// Helper class to hide connection parameters -/** - * - * \author M.Frank - * \version 1.0 - * \ingroup DD4HEP_DIGITIZATION - */ -class DigiROOTInput::internals_t { +class DigiROOTInput::inputsource_t { public: - typedef std::function<void(DataSegment&, int, const char*, void*)> func_t; - class converter_t { public: TBranch* branch { nullptr }; TClass* cls { nullptr }; - std::unique_ptr<func_t> call { }; }; - /// Mutex to allow re-use of a single source for multiple input streams - std::mutex m_lock { }; - +public: /// Branches present in the current file - std::map<unsigned long, converter_t> branches; + std::map<Key, converter_t> branches { }; /// Reference to the current ROOT file to be read - TFile* file { }; + TFile* file { nullptr }; /// Reference to the ROOT tree to be read TTree* tree { nullptr }; - /// Pointer to current input source - int input { INPUT_START }; /// Current entry inside the current input source Long64_t entry { -1 }; public: - /// Default constructor - internals_t () = default; - /// Default destructor - ~internals_t () { + inputsource_t() = default; + ~inputsource_t() { if ( file ) { file->Close(); delete file; file = nullptr; } } - /// Coverter function from input data to DDDige data - std::unique_ptr<func_t> input_converter(TClass* cl) const { - using detail::str_replace; - std::string cln = str_replace(cl->GetName(),"::","_"); - cln = str_replace(str_replace(str_replace(cln,'<','_'),'>','_'),'*',""); - std::string factory = "DD4hep_DDDigiConverter_"+cln; - func_t *fptr = (func_t*)PluginService::Create<void*>(factory.c_str()); - if ( 0 == fptr ) { - dd4hep::except("DigiROOTInput","Failed to access function pointer to read ROOT datatype: %s",cl->GetName()); - } - return std::unique_ptr<func_t>(fptr); + inputsource_t& next() { + ++entry; + tree->LoadTree( entry ); + return *this; } }; +/// Helper class to hide internal ROOT stuff +/** + * + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DIGITIZATION + */ +class DigiROOTInput::internals_t { +public: + using handle_t = std::unique_ptr<inputsource_t>; + /// Pointer to current input source + int m_curr_input { INPUT_START }; + /// Handle to input source + handle_t m_input_handle; + /// Reference to parent action + DigiROOTInput* m_input_action { nullptr }; + +public: + /// Default constructor + internals_t (DigiROOTInput* p); + /// Default destructor + ~internals_t () = default; + /// Access the next valid event entry + inputsource_t& next(); + /// Open the next input source + std::unique_ptr<inputsource_t> open_next_data_source(); +}; -/// Standard constructor -DigiROOTInput::DigiROOTInput(const DigiKernel& kernel, const std::string& nam) - : DigiInputAction(kernel, nam) +/// Default constructor +DigiROOTInput::internals_t::internals_t (DigiROOTInput* p) + : m_input_action(p) { - declareProperty("tree", m_tree_name = "EVENT"); - declareProperty("location", m_location = "inputs"); - declareProperty("containers", m_containers); - InstanceCount::increment(this); } -/// Default destructor -DigiROOTInput::~DigiROOTInput() { - InstanceCount::decrement(this); +DigiROOTInput::inputsource_t& DigiROOTInput::internals_t::next() { + if ( !m_input_handle ) { + m_input_handle = open_next_data_source(); + } + Int_t total = m_input_handle->tree->GetEntries(); + if ( (1+m_input_handle->entry) >= total ) { + m_input_handle.reset(); + m_input_handle = open_next_data_source(); + } + return m_input_handle->next(); } -/// Open new input file -void DigiROOTInput::open_new_file() const { - int len = this->m_inputs.size(); - imp = std::make_unique<internals_t>(); - if ( this->m_inputs.empty() ) imp->input = 0; - while ( (imp->input+1) < len ) { - const auto& fname = m_inputs[++imp->input]; - imp->file = TFile::Open(fname.c_str()); - if ( imp->file && imp->file->IsZombie() ) { - delete imp->file; - imp->file = nullptr; - error("OpenInput ++ Failed to open input source %s", fname.c_str()); +std::unique_ptr<DigiROOTInput::inputsource_t> +DigiROOTInput::internals_t::open_next_data_source() { + const auto& inputs = m_input_action->inputs(); + const auto& tree_name = m_input_action->tree_name(); + + int len = inputs.size(); + if ( inputs.empty() ) m_curr_input = 0; + while ( (m_curr_input+1) < len ) { + const auto& fname = inputs[++m_curr_input]; + auto* file = TFile::Open(fname.c_str()); + if ( file && file->IsZombie() ) { + delete file; + file = nullptr; + m_input_action->error("OpenInput ++ Failed to open input source %s", fname.c_str()); } - else if ( imp->file ) { - imp->tree = (TTree*)imp->file->Get(this->m_tree_name.c_str()); - if ( !imp->tree ) { - error("OpenInput ++ Failed to access tree: %s in input: %s", - this->m_tree_name.c_str(), fname.c_str()); + else if ( file ) { + auto* tree = (TTree*)file->Get(tree_name.c_str()); + if ( !tree ) { + m_input_action->error("OpenInput ++ Failed to access tree: %s in input: %s", + tree_name.c_str(), fname.c_str()); + continue; + } + Int_t total = tree->GetEntries(); + if ( total <= 0 ) { + m_input_action->error("OpenInput ++ TTree %s exists, but has no data. input: %s", + tree_name.c_str(), fname.c_str()); continue; } - imp->branches.clear(); - auto* branches = imp->tree->GetListOfBranches(); + auto source = std::make_unique<inputsource_t>(); + source->branches.clear(); + source->file = file; + source->tree = tree; + source->entry = -1; + + auto* branches = tree->GetListOfBranches(); + const auto& containers = m_input_action->container_names(); + int mask = m_input_action->input_mask(); TObjArrayIter it(branches); - for(Int_t i=0; i<branches->GetEntriesFast(); ++i) { + for(Int_t i=0; i < branches->GetEntriesFast(); ++i) { TBranch* b = (TBranch*)branches->At(i); TClass* cls = gROOT->GetClass( b->GetClassName(), kTRUE ); /// If there are no required branches, we convert everything - if ( this->m_containers.empty() ) { - Key key(b->GetName(), this->m_mask); - b->SetAutoDelete( kFALSE ); - imp->branches[key.value()] = {b, cls, imp->input_converter(cls)}; + if ( containers.empty() ) { + Key key(b->GetName(), mask); + source->branches[key] = {b, cls}; continue; } /// Otherwise only the entities asked for - for( const auto& bname : this->m_containers ) { + for( const auto& bname : containers ) { if ( bname == b->GetName() ) { - Key key(b->GetName(), this->m_mask); - b->SetAutoDelete( kFALSE ); - imp->branches[key.value()] = {b, cls, imp->input_converter(cls)}; + Key key(b->GetName(), mask); + source->branches[key] = {b, cls}; break; } } } - break; - } - else { - error("OpenInput ++ Failed to open input source %s", fname.c_str()); + if ( source->branches.empty() ) { + m_input_action->except("+++ No branches to be loaded. Configuration error!"); + } + return source; } - imp->file = nullptr; - } - if ( imp->input >= len ) { - imp.reset(); - except("OpenInput ++ Failed to open NEW input source"); } + m_input_action->except("+++ No open file present. Configuration error?"); + throw std::runtime_error("+++ No open file present"); +} + +/// Standard constructor +DigiROOTInput::DigiROOTInput(const DigiKernel& kernel, const std::string& nam) + : DigiInputAction(kernel, nam) +{ + imp = std::make_unique<internals_t>(this); + declareProperty("tree", m_tree_name = "EVENT"); + declareProperty("containers", m_containers); + InstanceCount::increment(this); +} + +/// Default destructor +DigiROOTInput::~DigiROOTInput() { + imp.reset(); + InstanceCount::decrement(this); } /// Pre-track action callback void DigiROOTInput::execute(DigiContext& context) const { - if ( imp ) { - std::lock_guard<std::mutex> lock(context.global_io_lock()); - Int_t total = imp->tree->GetEntries(); - if ( (1+imp->entry) >= total ) { - imp.reset(); - } - } - /// If necessary open a new file - if ( !imp ) { - std::lock_guard<std::mutex> lock(context.global_io_lock()); - open_new_file(); - } - - if ( nullptr == imp->file ) { - except("+++ No open file present. Configuration error?"); - } - if ( imp->branches.empty() ) { - except("+++ No branches to be loaded. Configuration error!"); - } - std::size_t input_len = 0; // // We have to lock all ROOT based actions. Consequences are SEGV otherwise. // - auto& event = context.event; std::lock_guard<std::mutex> lock(context.global_io_lock()); - ++imp->entry; - imp->tree->LoadTree( imp->entry ); + auto& event = context.event; + auto& source = imp->next(); + std::size_t input_len = 0; /// We only get here with a valid input - for( const auto& br : imp->branches ) { - void* obj = br.second.cls->New(); - br.second.branch->SetAddress(&obj); - } - DataSegment& segment = event->get_segment(this->m_location); - for( const auto& br : imp->branches ) { - TBranch* b = br.second.branch; - int nb = b->GetEntry( imp->entry ); - if ( nb < 0 ) { // This is definitely an error...ROOT says if reads fail, -1 is issued. - continue; + DataSegment& segment = event->get_segment(this->m_input_segment); + for( auto& b : source.branches ) { + auto& ent = b.second; + Long64_t bytes = ent.branch->GetEntry( source.entry ); + if ( bytes > 0 ) { + work_t work { segment, b.first, *ent.branch, *ent.cls }; + (*this)(context, work); + input_len += bytes; } - debug("%s+++ Loaded %8d bytes from branch %s", event->id(), nb, b->GetName()); - input_len += nb; - const auto& func = *br.second.call; - void** addr = (void**)b->GetAddress(); - func(segment, this->m_mask, b->GetName(), *addr); + debug("%s+++ Loaded %8ld bytes from branch %s", event->id(), bytes, ent.branch->GetName()); } info("%s+++ Read event %6ld [%ld bytes] from tree %s file: %s", - event->id(), imp->entry, input_len, imp->tree->GetName(), imp->file->GetName()); + event->id(), source.entry, input_len, source.tree->GetName(), source.file->GetName()); } diff --git a/DDDigi/src/DigiSegmentationTool.cpp b/DDDigi/src/DigiSegmentationTool.cpp index c6052d5e1f4a6c07cf659f4692df96efb199c951..a41d92856ce13e7ff9f402c7cddabe34837a31f7 100644 --- a/DDDigi/src/DigiSegmentationTool.cpp +++ b/DDDigi/src/DigiSegmentationTool.cpp @@ -192,7 +192,7 @@ DigiSegmentationTool::split_segmentation(const string& split_by) const "%-24s has %ld parallel entries when splitting by \"%s\"", det, segmentation_splits.size(), split_by.c_str()); stringstream str; - for( auto id : segmentation_splits ) + for( const auto& id : segmentation_splits ) str << setw(16) << hex << setfill('0') << id.first << " "; printout(INFO,"DigiSegmentationTool","%-24s --> Parallel Entries: %s", det, str.str().c_str()); diff --git a/DDDigi/src/DigiStoreDump.cpp b/DDDigi/src/DigiStoreDump.cpp index 5cdc4290f87c9d6c3bf19b9ac2bce1b770ce2940..a07881ee06573f955cb1195ad83a0ed0bb0117e0 100644 --- a/DDDigi/src/DigiStoreDump.cpp +++ b/DDDigi/src/DigiStoreDump.cpp @@ -84,9 +84,9 @@ DigiStoreDump::dump_history(DigiContext& context, Key history_key = item.history; std::vector<std::string> records; History::hist_entry_t hist { particle_key, 1.0 }; - const Particle& p = hist.get_particle(ev); - const Direction& mom = p.momentum; - const Position& vtx = p.start_position; + const Particle& par = hist.get_particle(ev); + const Direction& mom = par.momentum; + const Position& vtx = par.start_position; str << Key::key_name(key) << "[" << seq_no << "]:"; std::string line = @@ -96,7 +96,7 @@ DigiStoreDump::dump_history(DigiContext& context, long(&data.second)); records.emplace_back(line); line = format("| PDG:%6d Charge:%-2d Mass:%7.2f v:%7.5g %7.5g %7.5g p:%10g %10g %10g %016lX", - p.pdgID, int(p.charge), p.mass, vtx.X(), vtx.Y(), vtx.Z(), mom.X(), mom.Y(), mom.Z(), long(&p)); + par.pdgID, int(par.charge), par.mass, vtx.X(), vtx.Y(), vtx.Z(), mom.X(), mom.Y(), mom.Z(), long(&par)); records.emplace_back(line); return records; } @@ -120,30 +120,30 @@ DigiStoreDump::dump_history(DigiContext& context, cell, item.history.hits.size(), item.history.particles.size()); records.emplace_back(line); for( std::size_t i=0; i<item.history.hits.size(); ++i ) { - const auto& e = item.history.hits[i]; - const EnergyDeposit& dep = e.get_deposit(ev, container_key.item()); - const Position& p = dep.position; - const Position& m = dep.momentum; - Key k = e.source; + const auto& entry = item.history.hits[i]; + const EnergyDeposit& dep = entry.get_deposit(ev, container_key.item()); + const Position& pos = dep.position; + const Position& mom = dep.momentum; + Key k = entry.source; str.str(""); str << "| Hit-history[" << i << "]:"; line = format("%-30s Segment:%04X Mask:%04X Cell:%08X %.8g", - str.str().c_str(), k.segment(), k.mask(), k.item(), e.weight); + str.str().c_str(), k.segment(), k.mask(), k.item(), entry.weight); records.emplace_back(line); line = format("| pos: %7.3f %7.3f %7.3f p: %7.3f %7.3f %7.3f deposit: %7.3f", - p.X(), p.Y(), p.Z(), m.X(), m.Y(), m.Z(), dep.deposit); + pos.X(), pos.Y(), pos.Z(), mom.X(), mom.Y(), mom.Z(), dep.deposit); records.emplace_back(line); } for( std::size_t i=0; i<item.history.particles.size(); ++i ) { - const auto& e = item.history.particles[i]; - const Particle& par = e.get_particle(ev); + const auto& ent = item.history.particles[i]; + const Particle& par = ent.get_particle(ev); const Direction& mom = par.momentum; const Position& vtx = par.start_position; - Key k = e.source; + Key key = ent.source; str.str(""); str << "| Part-history[" << i << "]:"; line = format("%-30s Segment:%04X Mask:%04X Key: %08X %.8g", - str.str().c_str(), k.segment(), k.mask(), k.item(), e.weight); + str.str().c_str(), key.segment(), key.mask(), key.item(), ent.weight); records.emplace_back(line); line = format("| PDG:%6d Charge:%-2d Mass:%7.3f v:%7.3f %7.3f %7.3f p:%7.3f %7.3f %7.3f", par.pdgID, int(par.charge), par.mass, vtx.X(), vtx.Y(), vtx.Z(), mom.X(), mom.Y(), mom.Z()); @@ -212,9 +212,9 @@ void DigiStoreDump::dump_history(DigiContext& context, std::lock_guard<std::mutex> lock(segment.lock); records.emplace_back(format("+--- %-12s segment: %ld entries", tag.c_str(), segment.size())); - for ( const auto& e : segment ) { - Key key {e.first}; - const std::any& data = e.second; + for ( const auto& entry : segment ) { + Key key {entry.first}; + const std::any& data = entry.second; bool use = m_containers.empty() || std::find(m_container_items.begin(), m_container_items.end(), key.item()) != m_container_items.end(); use &= m_masks.empty() || @@ -257,9 +257,9 @@ void DigiStoreDump::dump_headers(const std::string& tag, std::vector<std::string> records; info("%s+--- %-12s segment: %ld entries", event.id(), tag.c_str(), segment.size()); std::lock_guard<std::mutex> lock(segment.lock); - for ( const auto& e : segment ) { - Key key {e.first}; - const std::any& data = e.second; + for ( const auto& entry : segment ) { + Key key {entry.first}; + const std::any& data = entry.second; if ( const auto* mapping = std::any_cast<DepositMapping>(&data) ) str = data_header(key, "deposits", *mapping); else if ( const auto* vector = std::any_cast<DepositVector>(&data) ) diff --git a/examples/DDDigi/scripts/DigiTest.py b/examples/DDDigi/scripts/DigiTest.py index cb3df5887d099c8a73982fef8382e80313d0d23f..a3ecf0c81e204bf5909c5bde4a0e3523c645c9e8 100644 --- a/examples/DDDigi/scripts/DigiTest.py +++ b/examples/DDDigi/scripts/DigiTest.py @@ -31,8 +31,10 @@ attenuation = {'SiVertexEndcapHits': 50 * ns, 'MuonBarrelHits': 50 * ns, 'BeamCalHits': 50 * ns, 'LumiCalHits': 50 * ns, - } +} + +# ========================================================================================================== class Test(dddigi.Digitize): def __init__(self, geometry=None): @@ -45,29 +47,38 @@ class Test(dddigi.Digitize): self.input = None self.main_sequencer() self.attenuation = attenuation - self.inputs = [ - 'CLICSiD_2022-11-01_15-54.root', - 'CLICSiD_2022-11-01_15-10.root', - 'CLICSiD_2022-10-31_17-26.root', - 'CLICSiD_2022-10-31_17-55.root', - 'CLICSiD_2022-10-31_17-55.root', - 'CLICSiD_2022-10-31_18-20.root', - 'CLICSiD_2022-10-31_18-40.root', - 'CLICSiD_2022-10-31_18-59.root', - 'CLICSiD_2022-10-31_19-18.root', - 'CLICSiD_2022-10-31_19-36.root', - 'CLICSiD_2022-10-31_19-53.root', - 'CLICSiD_2022-10-31_20-11.root'] self.used_inputs = [] + self.inputs = ['CLICSiD_2022-11-01_15-54.root', + 'CLICSiD_2022-11-01_15-10.root', + 'CLICSiD_2022-10-31_17-26.root', + 'CLICSiD_2022-10-31_17-55.root', + 'CLICSiD_2022-10-31_17-55.root', + 'CLICSiD_2022-10-31_18-20.root', + 'CLICSiD_2022-10-31_18-40.root', + 'CLICSiD_2022-10-31_18-59.root', + 'CLICSiD_2022-10-31_19-18.root', + 'CLICSiD_2022-10-31_19-36.root', + 'CLICSiD_2022-10-31_19-53.root', + 'CLICSiD_2022-10-31_20-11.root'] def segment_action(self, nam, **options): obj = dddigi.Interface.createSegmentAction(self.kernel(), str(nam)) return obj - def load_geo(self): + def load_geo(self, volume_manager=None): fname = "file:" + os.environ['DD4hepINSTALL'] + "/DDDetectors/compact/SiD.xml" self.kernel().loadGeometry(str(fname)) self.printDetectors() + if volume_manager: + vm = self.description.volumeManager() + if not vm.isValid(): + self.description.processXMLString(str("""<plugins> + <plugin name="DD4hep_VolumeManager"/> + </plugins>""")) + self.volumeManager = self.description.volumeManager() + if self.volumeManager.isValid(): + self.info('+++ Successfully created DD4hep VolumeManager') + return self def data_containers(self): return list(self.attenuation.keys()) @@ -112,3 +123,46 @@ class Test(dddigi.Digitize): result = "PASSED" self.always('%s Test finished after processing %d events.' % (result, evt_done,)) self.always('Test done. Exiting') + + +# ========================================================================================================== +def test_setup_1(digi): + """ + Create default setup for tests. Simply too bad to repeat the same lines over and over again. + + \author M.Frank + \version 1.0 + """ + # ======================================================================================================== + digi.info('Created SIGNAL input') + input = digi.input_action('DigiParallelActionSequence/READER') + input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0xCBAA, input=[digi.next_input()]) + # ======================================================================================================== + digi.info('Creating collision overlay....') + # ======================================================================================================== + overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-1') + overlay.adopt_action('DigiDDG4ROOT/Read-1', mask=0xCBEE, input=[digi.next_input()]) + digi.info('Created input.overlay-1') + # ======================================================================================================== + event = digi.event_action('DigiSequentialActionSequence/EventAction') + combine = event.adopt_action('DigiContainerCombine/Combine', + parallel=False, + input_masks=[0xCBAA, 0xCBEE], + output_mask=0xAAA0, + output_segment='deposits') + combine.erase_combined = False + proc = event.adopt_action('DigiContainerSequenceAction/HitP2', + parallel=False, + input_mask=0xAAA0, + input_segment='deposits', + output_mask=0xEEE5, + output_segment='deposits') + combine = digi.create_action('DigiDepositWeightedPosition/DepoCombine') + proc.adopt_container_processor(combine, digi.containers()) + conts = [c for c in digi.containers()] + event.adopt_action('DigiContainerDrop/Drop', + containers=conts, + input_segment='deposits', + input_masks=[0xAAA0]) + + return event diff --git a/examples/DDDigi/scripts/TestDepositCount.py b/examples/DDDigi/scripts/TestDepositCount.py index 3a7c490c787c55a081e4947e6b28e41b227c0ff7..a6be8433a4781d047255d09adcef38c70540e53e 100644 --- a/examples/DDDigi/scripts/TestDepositCount.py +++ b/examples/DDDigi/scripts/TestDepositCount.py @@ -19,7 +19,7 @@ def run(): # ======================================================================================================== digi.info('Created SIGNAL input') signal = input.adopt_action('DigiSequentialActionSequence/Signal') - reader = signal.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + reader = signal.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) sequence = signal.adopt_action('DigiContainerSequenceAction/Counter', parallel=True, input_mask=0x0, input_segment='inputs') count = digi.create_action('DigiCellMultiplicityCounter/CellCounter') diff --git a/examples/DDDigi/scripts/TestDepositSmearEnergy.py b/examples/DDDigi/scripts/TestDepositSmearEnergy.py new file mode 100644 index 0000000000000000000000000000000000000000..ff925982e86e0efd3bd9fb5cc62acf95324259ae --- /dev/null +++ b/examples/DDDigi/scripts/TestDepositSmearEnergy.py @@ -0,0 +1,47 @@ +# ========================================================================== +# 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. +# +# ========================================================================== +from __future__ import absolute_import + + +def run(): + import math + import DigiTest + from dd4hep import units + digi = DigiTest.Test(geometry=None) + + event = DigiTest.test_setup_1(digi) + proc = event.adopt_action('DigiContainerSequenceAction/Smearing', + parallel=False, + input_mask=0xEEE5, + input_segment='deposits', + output_mask=0xFFF0, + output_segment='outputs') + smear = digi.create_action('DigiDepositSmearing/Smear') + smear.deposit_cutoff = 1e-55 + # sigma(E)[GeV]/E[GeV] = 0.02 (systematic) + 0.005 / sqrt(E[GeV]) (intrinsic) + 1keV/E[GeV] (instrumentation) + # sigma(E)[GeV] = 0.02*E[GeV] (systematic) + 0.005 * sqrt(E[GeV]) (intrinsic) + 1keV (instrumentation) + smear.intrinsic_fluctuation = 0.005 / math.sqrt(units.GeV) + smear.systematic_resolution = 0.02 / units.GeV + smear.instrumentation_resolution = 1 * units.keV + smear.pair_ionisation_energy = 10 * units.eV + smear.ionization_fluctuation = True + smear.update_in_place = False + # proc.adopt_container_processor(smear, conts) + proc.adopt_container_processor(smear, ['EcalBarrelHits','EcalEndcapHits','HcalBarrelHits']) + + event.adopt_action('DigiStoreDump/HeaderDump') + # ======================================================================================================== + digi.info('Starting digitization core') + digi.run_checked(num_events=3, num_threads=-1, parallel=5) + + +if __name__ == '__main__': + run() diff --git a/examples/DDDigi/scripts/TestDepositWeighted.py b/examples/DDDigi/scripts/TestDepositWeighted.py index feb01bbfc4715f56908511fd16b6efc146bc32dd..067eab2140dffcc6dca9799f9ea000899a988971 100644 --- a/examples/DDDigi/scripts/TestDepositWeighted.py +++ b/examples/DDDigi/scripts/TestDepositWeighted.py @@ -18,18 +18,18 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================================================== digi.info('Created SIGNAL input') - signal = input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + signal = input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) digi.check_creation([signal]) # ======================================================================================================== digi.info('Creating collision overlays....') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-1') - evtreader = overlay.adopt_action('DigiROOTInput/Read-1', mask=0x1, input=[digi.next_input()]) + evtreader = overlay.adopt_action('DigiDDG4ROOT/Read-1', mask=0x1, input=[digi.next_input()]) digi.check_creation([overlay, evtreader]) digi.info('Created input.overlay-1') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-2') - evtreader = overlay.adopt_action('DigiROOTInput/Read-2', mask=0x2, input=[digi.next_input()]) + evtreader = overlay.adopt_action('DigiDDG4ROOT/Read-2', mask=0x2, input=[digi.next_input()]) digi.check_creation([overlay, evtreader]) digi.info('Created input.overlay-2') # ======================================================================================================== diff --git a/examples/DDDigi/scripts/TestHitHistory.py b/examples/DDDigi/scripts/TestHitHistory.py index 46c4d8ed0c6043f60a6b567295831df93fe4b944..8026f28b57402c3d314a13c5ebb8e0fc4feb0f58 100644 --- a/examples/DDDigi/scripts/TestHitHistory.py +++ b/examples/DDDigi/scripts/TestHitHistory.py @@ -18,18 +18,18 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================================================== digi.info('Created SIGNAL input') - input.adopt_action('DigiROOTInput/SignalReader', mask=0xCBAA, input=[digi.next_input()]) + input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0xCBAA, input=[digi.next_input()]) # ======================================================================================================== digi.info('Creating collision overlay....') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-1') - overlay.adopt_action('DigiROOTInput/Read-1', mask=0xCBEE, input=[digi.next_input()]) + overlay.adopt_action('DigiDDG4ROOT/Read-1', mask=0xCBEE, input=[digi.next_input()]) digi.info('Created input.overlay-1') # ======================================================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') combine = event.adopt_action('DigiContainerCombine/Combine', parallel=False, - input_masks=[0xCBAA,0xCBEE], + input_masks=[0xCBAA, 0xCBEE], output_mask=0xAAA0, output_segment='deposits') combine.erase_combined = False @@ -51,7 +51,7 @@ def run(): dump_history=True, containers=['SiTrackerBarrelHits', 'MCParticles'], segments=['deposits'], - masks=[0xAAA0,0xEEE5]) + masks=[0xAAA0, 0xEEE5]) # ======================================================================================================== digi.info('Starting digitization core') digi.run_checked(num_events=3, num_threads=1, parallel=5) diff --git a/examples/DDDigi/scripts/TestHitProcessing.py b/examples/DDDigi/scripts/TestHitProcessing.py index 0634c2094f775b7e073e9c1eebfb6225a618473e..755c8de6b2c019ca00d8d26e835708cc555ee085 100644 --- a/examples/DDDigi/scripts/TestHitProcessing.py +++ b/examples/DDDigi/scripts/TestHitProcessing.py @@ -18,16 +18,16 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================================================== digi.info('Created SIGNAL input') - input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) # ======================================================================================================== digi.info('Creating collision overlays....') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-1') - overlay.adopt_action('DigiROOTInput/Read-1', mask=0x1, input=[digi.next_input()]) + overlay.adopt_action('DigiDDG4ROOT/Read-1', mask=0x1, input=[digi.next_input()]) digi.info('Created input.overlay-1') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-2') - overlay.adopt_action('DigiROOTInput/Read-2', mask=0x2, input=[digi.next_input()]) + overlay.adopt_action('DigiDDG4ROOT/Read-2', mask=0x2, input=[digi.next_input()]) digi.info('Created input.overlay-2') # ======================================================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') diff --git a/examples/DDDigi/scripts/TestIPMove.py b/examples/DDDigi/scripts/TestIPMove.py index 0e84ef875d8b1905922b863f783ad5f782aa5870..cdfcf49bb78dd5b16affe306c9200439f65f785a 100644 --- a/examples/DDDigi/scripts/TestIPMove.py +++ b/examples/DDDigi/scripts/TestIPMove.py @@ -19,7 +19,7 @@ def run(): # ======================================================================================================== digi.info('Created SIGNAL input') signal = input.adopt_action('DigiSequentialActionSequence/Signal') - signal.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + signal.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) set_ip = signal.adopt_action('DigiIPCreate/SignalIP') set_ip.offset_ip = [1, 2, 3] set_ip.sigma_ip = [.5, .5, 3.0] diff --git a/examples/DDDigi/scripts/TestInput.py b/examples/DDDigi/scripts/TestInput.py index a35b18f93453a19c97ebb0d2a84eb632b3df02f4..c07424504b4ee4aa4f9bead3c20b489eb29f80f7 100644 --- a/examples/DDDigi/scripts/TestInput.py +++ b/examples/DDDigi/scripts/TestInput.py @@ -14,8 +14,8 @@ from __future__ import absolute_import def run(): import DigiTest digi = DigiTest.Test(geometry=None) - read = digi.input_action('DigiROOTInput/SignalReader', mask=0x0, input=['CLICSiD_2022-10-05_13-21.root']) - dump = digi.event_action('DigiStoreDump/StoreDump') + read = digi.input_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) + dump = digi.event_action('DigiStoreDump/StoreDump', parallel=False) digi.check_creation([read, dump]) digi.run_checked(num_events=5, num_threads=5, parallel=3) diff --git a/examples/DDDigi/scripts/TestMultiContainerParallel.py b/examples/DDDigi/scripts/TestMultiContainerParallel.py index 72fff4e167862e2b79091dccf0cb697068f73783..dc4a33e7eebefea9885955f0494426f49b801cd6 100644 --- a/examples/DDDigi/scripts/TestMultiContainerParallel.py +++ b/examples/DDDigi/scripts/TestMultiContainerParallel.py @@ -19,7 +19,7 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================== digi.info('Created SIGNAL input') - signal = input.adopt_action('DigiROOTInput/SignalReader', mask=0xFEED, input=[digi.next_input()]) + signal = input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0xFEED, input=[digi.next_input()]) digi.check_creation([signal]) # ======================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') diff --git a/examples/DDDigi/scripts/TestMultiInteractions.py b/examples/DDDigi/scripts/TestMultiInteractions.py index b20a45553f06c12b9d93830499a55a45a9c0b87d..cf4fe3152ee8edbe0ede260341905236e707468e 100644 --- a/examples/DDDigi/scripts/TestMultiInteractions.py +++ b/examples/DDDigi/scripts/TestMultiInteractions.py @@ -18,25 +18,25 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================================================== digi.info('Created SIGNAL input') - signal = input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + signal = input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) digi.check_creation([signal]) # ======================================================================================================== digi.info('Creating collision overlays....') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-1') - evtreader = overlay.adopt_action('DigiROOTInput/Reader-1', mask=0x1, input=[digi.next_input()]) + evtreader = overlay.adopt_action('DigiDDG4ROOT/Reader-1', mask=0x1, input=[digi.next_input()]) hist_drop = overlay.adopt_action('DigiHitHistoryDrop/Drop-1', masks=[evtreader.mask]) digi.check_creation([overlay, evtreader, hist_drop]) digi.info('Created input.overlay25') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-2') - evtreader = overlay.adopt_action('DigiROOTInput/Reader-2', mask=0x2, input=[digi.next_input()]) + evtreader = overlay.adopt_action('DigiDDG4ROOT/Reader-2', mask=0x2, input=[digi.next_input()]) hist_drop = overlay.adopt_action('DigiHitHistoryDrop/Drop-2', masks=[evtreader.mask]) digi.check_creation([overlay, evtreader, hist_drop]) digi.info('Created input.overlay50') # ======================================================================================================== overlay = input.adopt_action('DigiSequentialActionSequence/Overlay-3') - evtreader = overlay.adopt_action('DigiROOTInput/Reader-3', mask=0x3, input=[digi.next_input()]) + evtreader = overlay.adopt_action('DigiDDG4ROOT/Reader-3', mask=0x3, input=[digi.next_input()]) hist_drop = overlay.adopt_action('DigiHitHistoryDrop/Drop-3', masks=[evtreader.mask]) digi.check_creation([overlay, evtreader, hist_drop]) digi.info('Created input.overlay75') diff --git a/examples/DDDigi/scripts/TestPositionSmearResolution.py b/examples/DDDigi/scripts/TestPositionSmearResolution.py new file mode 100644 index 0000000000000000000000000000000000000000..63df6d078c3df775dc1b85ae250c656ad3268ce2 --- /dev/null +++ b/examples/DDDigi/scripts/TestPositionSmearResolution.py @@ -0,0 +1,42 @@ +# ========================================================================== +# 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. +# +# ========================================================================== +from __future__ import absolute_import + + +def run(): + import math + import DigiTest + from dd4hep import units + digi = DigiTest.Test(geometry=None) + digi.load_geo(volume_manager=True) + + event = DigiTest.test_setup_1(digi) + # ======================================================================================================== + proc = event.adopt_action('DigiContainerSequenceAction/Smearing', + parallel=False, + input_mask=0xEEE5, + input_segment='deposits') + smear = digi.create_action('DigiDepositSmearPositionResolution/Smear', + deposit_cutoff = 1e-55 * units.MeV, + resolution_u=1e-2 * units.mm, + resolution_v=1e-2 * units.mm) + proc.adopt_container_processor(smear, ['SiVertexEndcapHits','SiVertexBarrelHits', + 'SiTrackerEndcapHits','SiTrackerBarrelHits', + 'SiTrackerForwardHits']) + + event.adopt_action('DigiStoreDump/HeaderDump') + # ======================================================================================================== + digi.info('Starting digitization core') + digi.run_checked(num_events=3, num_threads=-1, parallel=5) + + +if __name__ == '__main__': + run() diff --git a/examples/DDDigi/scripts/TestPositionSmearTrack.py b/examples/DDDigi/scripts/TestPositionSmearTrack.py new file mode 100644 index 0000000000000000000000000000000000000000..7054e0056c93a459ff658b144e2596e5495194c9 --- /dev/null +++ b/examples/DDDigi/scripts/TestPositionSmearTrack.py @@ -0,0 +1,42 @@ +# ========================================================================== +# 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. +# +# ========================================================================== +from __future__ import absolute_import + + +def run(): + import math + import DigiTest + from dd4hep import units + digi = DigiTest.Test(geometry=None) + digi.load_geo(volume_manager=True) + + event = DigiTest.test_setup_1(digi) + # ======================================================================================================== + proc = event.adopt_action('DigiContainerSequenceAction/Smearing', + parallel=False, + input_mask=0xEEE5, + input_segment='deposits') + smear = digi.create_action('DigiDepositSmearPositionTrack/Smear', + deposit_cutoff = 1e-55 * units.MeV, + resolution_u=1e-2 * units.mm, + resolution_v=1e-2 * units.mm) + proc.adopt_container_processor(smear, ['SiVertexEndcapHits','SiVertexBarrelHits', + 'SiTrackerEndcapHits','SiTrackerBarrelHits', + 'SiTrackerForwardHits']) + + event.adopt_action('DigiStoreDump/HeaderDump') + # ======================================================================================================== + digi.info('Starting digitization core') + digi.run_checked(num_events=3, num_threads=-1, parallel=5) + + +if __name__ == '__main__': + run() diff --git a/examples/DDDigi/scripts/TestResegmentation.py b/examples/DDDigi/scripts/TestResegmentation.py index bbac3cf75760ed39b78c66913724d41117e969ec..4bba47447085e959528862d9602fe16c2d262e48 100644 --- a/examples/DDDigi/scripts/TestResegmentation.py +++ b/examples/DDDigi/scripts/TestResegmentation.py @@ -18,7 +18,7 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================== digi.info('Created SIGNAL input') - signal = input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + signal = input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) digi.check_creation([signal]) # ======================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') diff --git a/examples/DDDigi/scripts/TestSegmentationSplit.py b/examples/DDDigi/scripts/TestSegmentationSplit.py index 9ffb740a63fa3b17b70c0299ec31ee0d58c633fb..6ff313bad24285e4c942e1a9699b7c806c48ef02 100644 --- a/examples/DDDigi/scripts/TestSegmentationSplit.py +++ b/examples/DDDigi/scripts/TestSegmentationSplit.py @@ -18,7 +18,7 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================== digi.info('Created SIGNAL input') - signal = input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + signal = input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) digi.check_creation([signal]) # ======================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') diff --git a/examples/DDDigi/scripts/TestSimpleADCResponse.py b/examples/DDDigi/scripts/TestSimpleADCResponse.py index fb1529cc5817a86da00d612b7f43d11db3dac400..aa61c7d985d27c076aab96911df588ffe7214e09 100644 --- a/examples/DDDigi/scripts/TestSimpleADCResponse.py +++ b/examples/DDDigi/scripts/TestSimpleADCResponse.py @@ -17,7 +17,7 @@ def run(): # ======================================================================================================== input = digi.input_action('DigiSequentialActionSequence/READER') - input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) # ======================================================================================================== event = digi.event_action('DigiSequentialActionSequence/EventAction') event.adopt_action('DigiStoreDump/DumpInput') diff --git a/examples/DDDigi/scripts/TestSpillover.py b/examples/DDDigi/scripts/TestSpillover.py index e62e21c4bddda182b4addc7071dde87d0e7940f9..dd3a9a4b3c1c0845ecb95bd07229e731484cf2bf 100644 --- a/examples/DDDigi/scripts/TestSpillover.py +++ b/examples/DDDigi/scripts/TestSpillover.py @@ -19,13 +19,13 @@ def run(): input = digi.input_action('DigiParallelActionSequence/READER') # ======================================================================================================== - input.adopt_action('DigiROOTInput/SignalReader', mask=0x0, input=[digi.next_input()]) + input.adopt_action('DigiDDG4ROOT/SignalReader', mask=0x0, input=[digi.next_input()]) digi.info('Created input.signal') # ======================================================================================================== digi.info('Creating spillover sequence for EARLIER bunch crossings.....') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover-25') - evtreader = spillover.adopt_action('DigiROOTInput/Reader-25ns', mask=0x1, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-25ns', mask=0x1, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-25ns', t0=-25 * ns, signal_decay='exponential', @@ -38,7 +38,7 @@ def run(): digi.info('Created input.spillover-25') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover-50') - evtreader = spillover.adopt_action('DigiROOTInput/Reader-50ns', mask=0x2, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-50ns', mask=0x2, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-50ns', t0=-50 * ns, input_mask=evtreader.mask, containers=attenuation) hist_drop = spillover.adopt_action('DigiHitHistoryDrop/Drop-50ns', masks=[evtreader.mask]) @@ -46,7 +46,7 @@ def run(): digi.info('Created input.spillover-50') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover-75') - evtreader = spillover.adopt_action('DigiROOTInput/Reader-75ns', mask=0x3, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader-75ns', mask=0x3, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att-75ns', t0=-75 * ns, input_mask=evtreader.mask, containers=attenuation) hist_drop = spillover.adopt_action('DigiHitHistoryDrop/Drop-75ns', masks=[evtreader.mask]) @@ -56,7 +56,7 @@ def run(): digi.info('Creating spillover sequence for LATER bunch crossings.....') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover+25') - evtreader = spillover.adopt_action('DigiROOTInput/Reader+25ns', mask=0x4, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader+25ns', mask=0x4, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att+25ns', t0=25 * ns, input_mask=evtreader.mask, containers=attenuation) hist_drop = spillover.adopt_action('DigiHitHistoryDrop/Drop+25ns', masks=[evtreader.mask]) @@ -64,7 +64,7 @@ def run(): digi.info('Created input.spillover+25') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover+50') - evtreader = spillover.adopt_action('DigiROOTInput/Reader+50ns', mask=0x5, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader+50ns', mask=0x5, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att+50ns', t0=50 * ns, input_mask=evtreader.mask, containers=attenuation) hist_drop = spillover.adopt_action('DigiHitHistoryDrop/Drop_50ns', masks=[evtreader.mask]) @@ -72,7 +72,7 @@ def run(): digi.info('Created input.spillover+50') # ======================================================================================================== spillover = input.adopt_action('DigiSequentialActionSequence/Spillover+75') - evtreader = spillover.adopt_action('DigiROOTInput/Reader+75ns', mask=0x6, input=[digi.next_input()]) + evtreader = spillover.adopt_action('DigiDDG4ROOT/Reader+75ns', mask=0x6, input=[digi.next_input()]) attenuate = spillover.adopt_action('DigiAttenuatorSequence/Att+75ns', t0=75 * ns, input_mask=evtreader.mask, containers=attenuation) hist_drop = spillover.adopt_action('DigiHitHistoryDrop/Drop+75ns', masks=[evtreader.mask])