diff --git a/Simulation/CMakeLists.txt b/Simulation/CMakeLists.txt index 3ce6540c7301bfc59af6ade09e7cc41de7802f23..673b31237aa8327da637571e08fea3d08f1e94a5 100644 --- a/Simulation/CMakeLists.txt +++ b/Simulation/CMakeLists.txt @@ -6,3 +6,4 @@ add_subdirectory(DetSimFastModel) add_subdirectory(DetSimGeom) add_subdirectory(DetSimInterface) add_subdirectory(DetSimSD) +add_subdirectory(DetSimMixing) \ No newline at end of file diff --git a/Simulation/DetSimMixing/CMakeLists.txt b/Simulation/DetSimMixing/CMakeLists.txt new file mode 100755 index 0000000000000000000000000000000000000000..bfaea3d388836406822388202c23122e72b439c2 --- /dev/null +++ b/Simulation/DetSimMixing/CMakeLists.txt @@ -0,0 +1,13 @@ +gaudi_add_module(DetSimMixing + SOURCES src/DetSimMixingAlg.cc + LINK DetSimInterface + EDM4HEP::edm4hep EDM4HEP::edm4hepDict + podio::podioRootIO + Gaudi::GaudiKernel +) + +install(TARGETS DetSimMixing + EXPORT CEPCSWTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) diff --git a/Simulation/DetSimMixing/src/BackgroundBatch.hh b/Simulation/DetSimMixing/src/BackgroundBatch.hh new file mode 100755 index 0000000000000000000000000000000000000000..b8f334b2afafc2c06ad5da0c01d089c3ccf689b1 --- /dev/null +++ b/Simulation/DetSimMixing/src/BackgroundBatch.hh @@ -0,0 +1,22 @@ +#ifndef BackgroundBatch_hh +#define BackgroundBatch_hh + +#include <cstdint> +#include <vector> +#include "TTimeStamp.h" + +class BackgroundBatch { +public: + TTimeStamp start_time; // the start time of the batch + double duration; // the duration of the batch (in ns) + int num_events; // the number of events in the batch + + // below are the metadata of each type of backgrouds in the batch + std::vector<int> event_numbers_groupby_type; // the number of events for each type + + // below are the metadata of each event in the batch + std::vector<size_t> event_types; // the event types in the batch + std::vector<double> event_times; // the relative event times in the batch (in ns) +}; + +#endif diff --git a/Simulation/DetSimMixing/src/BackgroundLoader.hh b/Simulation/DetSimMixing/src/BackgroundLoader.hh new file mode 100755 index 0000000000000000000000000000000000000000..792ac1e7195f8fb843d1911f64a8609ba1bb575b --- /dev/null +++ b/Simulation/DetSimMixing/src/BackgroundLoader.hh @@ -0,0 +1,84 @@ +#ifndef BackgroundLoader_hh +#define BackgroundLoader_hh + +#include "IBackgroundLoader.hh" +#include <iostream> +#include <vector> +#include <string> + +#include <podio/Frame.h> +#include <podio/ROOTFrameReader.h> + +class BackgroundLoader: public IBackgroundLoader { +public: + + BackgroundLoader(const std::vector<std::string>& filenames) { + m_reader.openFiles(filenames); + m_max_evt = m_reader.getEntries(podio::Category::Event); + } + virtual ~BackgroundLoader() = default; + + bool fill(BackgroundEvent& evt, const double current_time_in_ns) override { + // before load the next event, check whether end or not. + // if reach the end, then go back to the first event. + if (m_current_evt >= m_max_evt) { + m_current_evt = 0; + } + + // get the frame of the current event + auto frame = podio::Frame(m_reader.readEntry(podio::Category::Event, m_current_evt)); + + // before unpack the hits, check the collections name to index map. + if (evt.collection_index.empty()) { + size_t idx = 0; + for (auto col_name: frame.getAvailableCollections()) { + // std::cout << "available collection: " << col_name << std::endl; + evt.collection_index[col_name] = idx; + ++idx; + } + } + + for (auto [name, colidx]: evt.collection_index) { + auto col_ = frame.get(name); + + if (auto col = dynamic_cast<const edm4hep::SimTrackerHitCollection*>(col_)) { + auto& trk_col = evt.tracker_hits[colidx]; + for (auto oldhit: *col) { + auto newhit = oldhit.clone(); + newhit.setTime(oldhit.getTime() + current_time_in_ns); + newhit.setOverlay(true); + trk_col.push_back(newhit); + } + + } else if (auto col = dynamic_cast<const edm4hep::SimCalorimeterHitCollection*>(col_)) { + auto& calo_col = evt.calorimeter_hits[colidx]; + for (auto oldhit: *col) { + auto newhit = oldhit.clone(); + // loop all the contributions and add the time. + for (auto contrib: oldhit.getContributions()) { + auto newcontrib = contrib.clone(); + newcontrib.setTime(contrib.getTime() + current_time_in_ns); + newhit.addToContributions(newcontrib); + } + calo_col.push_back(newhit); + } + } else { + continue; + } + + } + + ++m_current_evt; + return true; + } + +private: + podio::ROOTFrameReader m_reader; + + unsigned int m_current_evt{0}; + unsigned int m_max_evt{0}; + +}; + + +#endif \ No newline at end of file diff --git a/Simulation/DetSimMixing/src/DetSimMixingAlg.cc b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc new file mode 100755 index 0000000000000000000000000000000000000000..f054f3b059631c8accba8b3d655d83d1f957625c --- /dev/null +++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc @@ -0,0 +1,167 @@ +#include "DetSimMixingAlg.hh" +#include "GaudiKernel/IEventProcessor.h" +#include "GaudiKernel/IAppMgrUI.h" +#include "GaudiKernel/GaudiException.h" +#include "GaudiKernel/IRndmEngine.h" + +#include <CLHEP/Random/RandExponential.h> +#include <CLHEP/Random/RandFlat.h> +#include <CLHEP/Random/RandPoisson.h> + +#include "BackgroundLoader.hh" + +DECLARE_COMPONENT(DetSimMixingAlg) + +DetSimMixingAlg::DetSimMixingAlg(const std::string& name, ISvcLocator* pSvcLocator) +: Algorithm(name, pSvcLocator) { +} + +StatusCode DetSimMixingAlg::initialize() { + StatusCode sc; + + info() << "Initialize DetSimMixingAlg... " << endmsg; + // preparation according to user properties + if (m_background_rates.value().size() != m_background_filelists.value().size()) { + error() << "The size of the background rates and filelists should be the same." << endmsg; + return StatusCode::FAILURE; + } + + for (auto [type, rate]: m_background_rates.value()) { + if (m_background_filelists.value().find(type) == m_background_filelists.value().end()) { + error() << "The input file lists for the background type " << type << " is not provided." << endmsg; + return StatusCode::FAILURE; + } + m_event_types.push_back(type); + m_event_rates.push_back(rate); + m_input_lists.push_back(m_background_filelists.value()[type]); + } + + // prepare the loaders + for (size_t i = 0; i < m_event_types.size(); ++i) { + m_total_rates += m_event_rates[i]; + m_event_loaders.push_back(new BackgroundLoader(m_input_lists[i])); + } + + if (m_total_rates <= 0) { + error() << "The total rate of the background events should be positive." << endmsg; + return StatusCode::FAILURE; + } + + m_total_tau = 1.0 / m_total_rates; + + info() << "Summary of the background events: " << endmsg; + for (size_t i = 0; i < m_event_types.size(); ++i) { + info() << " Event type: " << m_event_types[i] << ", rate: " << m_event_rates[i] << " Hz" << endmsg; + } + + return sc; +} + +StatusCode DetSimMixingAlg::execute() { + StatusCode sc; + + info() << "Execute DetSimMixingAlg... " << endmsg; + + // ======================================================================== + // Prepare the batches + // ======================================================================== + std::vector<BackgroundBatch> batches; + + int nbatches = 5; + double duration = 1e3; // 1 us = 1,000 ns of each batch + TTimeStamp start_time; + + for (int i = 0; i < nbatches; i++) { + BackgroundBatch batch; + batch.start_time = start_time + TTimeStamp(0, static_cast<int>(i*duration)); + batch.duration = duration; + + // todo: need to sample according to the rate + batch.num_events = 0; + for (size_t evttype = 0; evttype < m_event_types.size(); ++evttype) { + batch.event_numbers_groupby_type.push_back(0); + } + + double current_time = 0; + + while (true) { + // sampling and get an event + size_t selected_evttype = 0; + double r = CLHEP::RandFlat::shoot(m_total_rates); + double accumulated = 0; + for (size_t evttype = 0; evttype < m_event_types.size(); ++evttype) { + accumulated += m_event_rates[evttype]; + if (r < accumulated) { + selected_evttype = evttype; + break; + } + } + + // sampling the time + current_time += CLHEP::RandExponential::shoot(m_total_tau)*1e9; // ns + + if (current_time > duration) { + break; + } + + // add the event to Batch. + ++batch.num_events; + batch.event_types.push_back(selected_evttype); + batch.event_times.push_back(current_time); + + batch.event_numbers_groupby_type[selected_evttype] += 1; + } + + batches.push_back(batch); + } + + + // ======================================================================== + // Unpack the hits from Loader and fill the existing hits into collections + // of BackgroundEvent. + // ======================================================================== + BackgroundEvent bkg_evt; + + for (size_t i = 0; i < batches.size(); ++i) { + const BackgroundBatch& batch = batches[i]; + + for (size_t j = 0; j < batch.num_events; ++j) { + size_t evttype = batch.event_types[j]; + double evttime = batch.event_times[j]; + + IBackgroundLoader* loader = m_event_loaders[evttype]; + + if (not loader->fill(bkg_evt, evttime)) { + error() << "Failed to load the next event." << endmsg; + return StatusCode::FAILURE; + } + + } + } + + // dump the information + for (auto [name, idx]: bkg_evt.collection_index) { + info() << "Collection: " << name << ", index: " << idx; + if (bkg_evt.tracker_hits.count(idx)) { + info() << " Tracker hits: " << bkg_evt.tracker_hits[idx].size() << endmsg; + } else if (bkg_evt.calorimeter_hits.count(idx)) { + info() << " Calorimeter hits: " << bkg_evt.calorimeter_hits[idx].size() << endmsg; + } else { + info() << endmsg; + } + } + + // ======================================================================== + // Put the BackgroundEvent into the event store + // ======================================================================== + + return sc; +} + +StatusCode DetSimMixingAlg::finalize() { + StatusCode sc; + + info() << "Finalize DetSimMixingAlg... " << endmsg; + + return sc; +} \ No newline at end of file diff --git a/Simulation/DetSimMixing/src/DetSimMixingAlg.hh b/Simulation/DetSimMixing/src/DetSimMixingAlg.hh new file mode 100755 index 0000000000000000000000000000000000000000..d678dca31c93e7c989d8b32ea3129fa5fc5f100b --- /dev/null +++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.hh @@ -0,0 +1,76 @@ +#ifndef DetSimMixingAlg_hh +#define DetSimMixingAlg_hh +/* + * Description: + * This algorithm is used to generate a mixing event. + * Assume there will be one physics event and multiple background events. + * + * The main algorithm first creates enough BackgroundBatch instances. + * All the metadata will be prepared into the batch. + * Then, the main algorithm will start the mixing and load the real hits. + * The reason to unpack the hits in main algorithm is to reduce the data copy. + * + * About the BackgroundBatch: + * + * -----------------------------------------------> time + * | o x o o x | + * | | | + * | | \-> group by event type 'x', num of events = 2 + * | -> group by event type 'o', num of events = 3 + * \-> start time + * + * For each event, it consists of several collections. For one batch: + * + * | EventType | Collection 1 | Collection 2 | ... | Collection N | + * | o | | | ... | | + * | x | | | ... | | + * | o | | | ... | | + * | o | | | ... | | + * | x | | | ... | | + * + * The main algorithm needs to loop BackgroundLoader and unpack the collections. + * For each iteration in the loop of the above table, a current event could be accessed. + * According to the current event, the main algorithm could unpack the collections. + * + * Authors: + * Tao Lin <lintao AT ihep.ac.cn> + */ + +#include <string> +#include <vector> +#include <map> + +#include <GaudiKernel/Algorithm.h> +#include <Gaudi/Property.h> +#include <GaudiKernel/ToolHandle.h> + +#include "k4FWCore/DataHandle.h" + +#include "BackgroundBatch.hh" +#include "IBackgroundLoader.hh" + +class DetSimMixingAlg: public Algorithm { +public: + DetSimMixingAlg(const std::string& name, ISvcLocator* pSvcLocator); + + StatusCode initialize() override; + StatusCode execute() override; + StatusCode finalize() override; + +private: + double m_total_rates{0}; // Hz. + double m_total_tau{0}; // s. the total time of the background events + +private: + std::vector<std::string> m_event_types; // the event types of the background events + std::vector<double> m_event_rates; // the rates of the background events + std::vector<IBackgroundLoader*> m_event_loaders; // the loaders of the background events + std::vector<std::vector<std::string>> m_input_lists; // input files for each backgroud + +private: + // properties for user side + Gaudi::Property<std::map<std::string, double>> m_background_rates{this, "BackgroundRates", {}, "The rates of the background events"}; + Gaudi::Property<std::map<std::string, std::vector<std::string>>> m_background_filelists{this, "BackgroundFileLists", {}, "The input file lists for the background events"}; +}; + +#endif diff --git a/Simulation/DetSimMixing/src/IBackgroundLoader.hh b/Simulation/DetSimMixing/src/IBackgroundLoader.hh new file mode 100755 index 0000000000000000000000000000000000000000..588d928b0d8c3fb3910fbde501e588ce5f7b96ac --- /dev/null +++ b/Simulation/DetSimMixing/src/IBackgroundLoader.hh @@ -0,0 +1,38 @@ +#ifndef IBackgroundLoader_hh +#define IBackgroundLoader_hh + +#include "BackgroundBatch.hh" +#include <map> + +#include "edm4hep/SimTrackerHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/CaloHitContributionCollection.h" + +/* + * Description: + * As there are multiple collections for one event, this BackgrounEvent + * is used to organize these collections. + * + * The key is the index of collection. + * + * BackgroundLoader is resposible to fill its hit collections into the + * BackgroundEvent. + */ +struct BackgroundEvent { + std::map<size_t, edm4hep::SimTrackerHitCollection> tracker_hits; + std::map<size_t, edm4hep::SimCalorimeterHitCollection> calorimeter_hits; + std::map<size_t, edm4hep::CaloHitContributionCollection> calo_contribs; + + std::map<std::string, size_t> collection_index; // key is collection name, value is index. +}; + +class IBackgroundLoader { +public: + virtual ~IBackgroundLoader() = default; + + // the loader is responsible to unpack the hits into background event. + // if failed, then return false. + virtual bool fill(BackgroundEvent& bkg_event, const double current_time_in_ns) = 0; +}; + +#endif \ No newline at end of file