diff --git a/Simulation/DetSimMixing/CMakeLists.txt b/Simulation/DetSimMixing/CMakeLists.txt old mode 100755 new mode 100644 diff --git a/Simulation/DetSimMixing/src/BackgroundBatch.hh b/Simulation/DetSimMixing/src/BackgroundBatch.hh old mode 100755 new mode 100644 index bdd16ba14e18dbca6bd5a1aafaf467a62ccad9ef..9672e6a5408e237a224db4b6bde5afac36ef83a1 --- a/Simulation/DetSimMixing/src/BackgroundBatch.hh +++ b/Simulation/DetSimMixing/src/BackgroundBatch.hh @@ -5,8 +5,7 @@ #include <vector> #include "TTimeStamp.h" -class BackgroundBatch { -public: +struct BackgroundBatch { double start_time; // ns, the start time of the batch, relative to signal event double duration; // the duration of the batch int num_events; // the number of events in the batch diff --git a/Simulation/DetSimMixing/src/BackgroundEvent.hh b/Simulation/DetSimMixing/src/BackgroundEvent.hh new file mode 100644 index 0000000000000000000000000000000000000000..c6a8bea9ed6b55b14749aaf960f447b06742829f --- /dev/null +++ b/Simulation/DetSimMixing/src/BackgroundEvent.hh @@ -0,0 +1,55 @@ +#ifndef BackgroundEvent_hh +#define BackgroundEvent_hh + +#include <map> +#include <vector> +#include <string> +#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. + std::map<std::string, size_t> collection_subdet; // key is collection name, value is subdetector type. + + // in order to reduce the memory usage, need additional filter setup here. + // for example, for TPC, we need to load all the hits. + // however, for the tracker/calo, we only need to load the hits which are + // in an interested event time window. + enum SubDetType { + kVXD = 0, // 200 ns + kITK = 1, // 200 ns + kTPC = 2, // 34 us + kOTK = 3, // 1 us + kECAL = 4, // 150 ns + kHCAL = 5, // 1 us + kMUON = 6, // 100 ns // from Xiaolong Wang + kNSubDetType + }; + std::vector<double> subdet2twindow = {200, 200, 34000, 1000, 150, 1000, 100}; // key is subdet, value is time window. + std::vector<std::vector<std::string>> subdet2colnames = { + {"VXD"}, + {"ITK"}, + {"TPC"}, + {"OTK"}, + {"Ecal"}, + {"Hcal"}, + {"Muon"} + }; +}; + +#endif diff --git a/Simulation/DetSimMixing/src/BackgroundLoader.hh b/Simulation/DetSimMixing/src/BackgroundLoader.hh old mode 100755 new mode 100644 index 7ab42bf9c7128ed71d9378ce847bc98ab2e3b112..6839ec53215523ae09148156ca6a5ec65e108034 --- a/Simulation/DetSimMixing/src/BackgroundLoader.hh +++ b/Simulation/DetSimMixing/src/BackgroundLoader.hh @@ -2,6 +2,7 @@ #define BackgroundLoader_hh #include "IBackgroundLoader.hh" +#include "BackgroundEvent.hh" #include <iostream> #include <vector> #include <string> diff --git a/Simulation/DetSimMixing/src/DetSimMixingAlg.cc b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc old mode 100755 new mode 100644 index 23c9bcb56c4df945ef938afa23a5c3904891b076..cd0e030fbed29ee5fff773ce49475a7b4072ee77 --- a/Simulation/DetSimMixing/src/DetSimMixingAlg.cc +++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc @@ -9,6 +9,7 @@ #include <CLHEP/Random/RandPoisson.h> #include "BackgroundLoader.hh" +#include "BackgroundEvent.hh" DECLARE_COMPONENT(DetSimMixingAlg) @@ -21,26 +22,26 @@ StatusCode DetSimMixingAlg::initialize() { info() << "Initialize DetSimMixingAlg... " << endmsg; // preparation according to user properties - if (m_background_rates.value().size() != m_background_filelists.value().size()) { + if (m_background_timings.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()) { + for (auto [type, timing]: m_background_timings.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_event_timings.push_back(timing); m_input_lists.push_back(m_background_filelists.value()[type]); } // prepare the loaders for (size_t i = 0; i < m_event_types.size(); ++i) { // only the positive rates are considered into total rates - if (m_event_rates[i] > 0) { - m_total_rates += m_event_rates[i]; + if (m_event_timings[i] > 0) { + m_total_rates += m_event_timings[i]; } m_event_loaders.push_back(new BackgroundLoader(m_input_lists[i])); } @@ -55,10 +56,10 @@ StatusCode DetSimMixingAlg::initialize() { } info() << "Summary of the background events: " << endmsg; for (size_t i = 0; i < m_event_types.size(); ++i) { - if (m_event_rates[i] > 0) { - info() << " Event type: " << m_event_types[i] << ", rate: " << m_event_rates[i] << " Hz" << endmsg; + if (m_event_timings[i] > 0) { + info() << " Event type: " << m_event_types[i] << ", rate: " << m_event_timings[i] << " Hz" << endmsg; } else { - info() << " Event type: " << m_event_types[i] << ", time window: " << fabs(m_event_rates[i]) << " ns" << endmsg; + info() << " Event type: " << m_event_types[i] << ", time window: " << fabs(m_event_timings[i]) << " ns" << endmsg; } } @@ -110,8 +111,9 @@ StatusCode DetSimMixingAlg::execute() { // ======================================================================== std::vector<BackgroundBatch> batches; - int nbatches = 10; // the total number of batches will be 2*nbatches, [-nbatches, nbatches) - double duration = 3460; // duration of each batch (ns): Nbunch x TBunchSpacing = 10 x 346 ns = 3460 ns + int nbatches = m_nbatches.value(); // the total number of batches will be 2*nbatches, [-nbatches, nbatches) + int nbunches = m_nbunches_per_batch.value(); // the total number of bunches per batch + double duration = nbunches * m_bunch_crossing_spacing.value(); // duration of each batch (ns): Nbunch x TBunchSpacing = 10 x 277 ns = 2770 ns for (int i = -nbatches; i < nbatches; ++i) { BackgroundBatch batch; @@ -132,10 +134,10 @@ StatusCode DetSimMixingAlg::execute() { // insert at beginning, align to the begin of batch for (size_t evttype = 0; evttype < m_event_types.size(); ++evttype) { // skip the sampled mode type - if (m_event_rates[evttype] > 0) { + if (m_event_timings[evttype] > 0) { continue; } - double time_window_event = fabs(m_event_rates[evttype]); + double time_window_event = fabs(m_event_timings[evttype]); int n_events = std::max(1, static_cast<int>(batch.duration/time_window_event)); // ns for (size_t j = 0; j < n_events; ++j) { @@ -157,10 +159,10 @@ StatusCode DetSimMixingAlg::execute() { double accumulated = 0; for (size_t evttype = 0; evttype < m_event_types.size(); ++evttype) { // skip the fixed time mode type - if (m_event_rates[evttype] <= 0) { + if (m_event_timings[evttype] <= 0) { continue; } - accumulated += m_event_rates[evttype]; + accumulated += m_event_timings[evttype]; if (r < accumulated) { selected_evttype = evttype; break; @@ -191,6 +193,15 @@ StatusCode DetSimMixingAlg::execute() { // of BackgroundEvent. // ======================================================================== BackgroundEvent bkg_evt; + // setup the time window + bkg_evt.subdet2twindow[BackgroundEvent::kVXD] = m_vxd_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kITK] = m_itk_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kTPC] = m_tpc_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kOTK] = m_otk_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kECAL] = m_ecal_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kHCAL] = m_hcal_time_window.value(); + bkg_evt.subdet2twindow[BackgroundEvent::kMUON] = m_muon_time_window.value(); + info() << "Creating a BackgroundEvent..." << endmsg; for (size_t i = 0; i < batches.size(); ++i) { info() << " Batch " << i << ": " << endmsg; @@ -228,8 +239,34 @@ StatusCode DetSimMixingAlg::execute() { // Put the BackgroundEvent into the event store // ======================================================================== // add prefix 'Mixed'. - for (auto [col_name, colidx]: bkg_evt.collection_index) { + + // Need to prepare a complete list of collection names. + std::map<std::string, size_t> col_names; + for (auto [name, idx]: bkg_evt.collection_index) { + ++col_names[name]; + } + for (auto [name, idx]: m_sig_trackerColMap) { + if (col_names.find(name) == col_names.end()) { + warning() << "Collection " << name << " is not in the background event." << endmsg; + } + ++col_names[name]; + } + for (auto [name, idx]: m_sig_calorimeterColMap) { + if (col_names.find(name) == col_names.end()) { + warning() << "Collection " << name << " is not in the background event." << endmsg; + } + ++col_names[name]; + } + + for (auto [col_name, cnt]: col_names) { + + if (bkg_evt.collection_index.find(col_name) == bkg_evt.collection_index.end()) { + debug() << "Collection " << col_name << " is not in the background event." << endmsg; + continue; + } + auto colidx = bkg_evt.collection_index[col_name]; + if (bkg_evt.tracker_hits.count(colidx)) { auto& col = bkg_evt.tracker_hits[colidx]; @@ -252,7 +289,17 @@ StatusCode DetSimMixingAlg::execute() { if (!sig_col) { continue; } + + // keep the same time windows for signal and backgrounds. + auto time_window = bkg_evt.subdet2twindow[bkg_evt.collection_subdet[col_name]]; + for (auto oldhit: *sig_col) { + auto t = oldhit.getTime(); + if (t < -time_window || t > time_window) { + // if the hit is not in the time window, skip. + continue; + } + auto newhit = newcol->create(); newhit.setCellID(oldhit.getCellID()); newhit.setEDep(oldhit.getEDep()); @@ -295,7 +342,25 @@ StatusCode DetSimMixingAlg::execute() { continue; } + // keep the same time windows for signal and backgrounds. + auto time_window = bkg_evt.subdet2twindow[bkg_evt.collection_subdet[col_name]]; + + for (auto oldhit: *sig_col) { + // check whether the hit is in the time window. + bool is_in_window = false; + for (auto contrib: oldhit.getContributions()) { + auto t = contrib.getTime(); + if (t < -time_window || t > time_window) { + continue; + } + is_in_window = true; + break; + } + if (not is_in_window) { + continue; + } + auto newhit = newcol->create(); newhit.setCellID(oldhit.getCellID()); newhit.setEnergy(oldhit.getEnergy()); @@ -313,7 +378,7 @@ StatusCode DetSimMixingAlg::execute() { } } else { - error() << "The collection " << col_name + debug() << "The collection " << col_name << " with idx" << colidx << " is not found." << endmsg; continue; diff --git a/Simulation/DetSimMixing/src/DetSimMixingAlg.hh b/Simulation/DetSimMixing/src/DetSimMixingAlg.hh old mode 100755 new mode 100644 index 55aa56a7def63ac36ddd554ca49d161a81b41171..d50c9b3f17332f1992a0a47e27a0b92c33117d6f --- a/Simulation/DetSimMixing/src/DetSimMixingAlg.hh +++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.hh @@ -46,8 +46,9 @@ #include "k4FWCore/DataHandle.h" -#include "BackgroundBatch.hh" #include "IBackgroundLoader.hh" +#include "BackgroundBatch.hh" +#include "BackgroundEvent.hh" class DetSimMixingAlg: public Algorithm { public: @@ -63,7 +64,7 @@ private: 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<double> m_event_timings; // the rates or durations 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 @@ -94,12 +95,29 @@ private: // Following properties are used to configure different types of background events. // * Key: label of one type // * Value: - // * rate: + // * timing: // * FixedTimeWindow: if the rate is less than 0, then always mix the background events. fabs(rate) is the time window of the background events. // * Sampled: if the rate is larger than 0, then the background events will be sampled according to the rate. // * filelist: the input file list for this type of background events. - Gaudi::Property<std::map<std::string, double>> m_background_rates{this, "BackgroundRates", {}, "The rates (positive, Hz) or time windows (negative, ns) of the background events"}; + Gaudi::Property<std::map<std::string, double>> m_background_timings{this, "BackgroundTimings", {}, "The rates (positive, Hz) or time windows (negative, ns) 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"}; + + + // Bunch Crossing information + Gaudi::Property<double> m_bunch_crossing_spacing{this, "BunchCrossingSpacing", 277.0, "The bunch crossing spacing in ns"}; + Gaudi::Property<int> m_nbunches_per_batch{this, "NbunchesPerBatch", 10, "The number of bunches per batch"}; + Gaudi::Property<int> m_nbatches{this, "Nbatches", 13, "The number of batches to be mixed"}; // 34us/(277ns*10) = ~13 + +private: + // Time window for VXD, ITK, TPC, OTK, ECAL, HCAL, MUON + // Only the hits between [-T, T] ns are loaded + Gaudi::Property<double> m_vxd_time_window{this, "VXDTimeWindow", 200.0, "The time window for VXD in ns"}; + Gaudi::Property<double> m_itk_time_window{this, "ITKTimeWindow", 200.0, "The time window for ITK in ns"}; + Gaudi::Property<double> m_tpc_time_window{this, "TPCTimeWindow", 34000.0, "The time window for TPC in ns"}; + Gaudi::Property<double> m_otk_time_window{this, "OTKTimeWindow", 1000.0, "The time window for OTK in ns"}; + Gaudi::Property<double> m_ecal_time_window{this, "EcalTimeWindow", 150.0, "The time window for ECAL in ns"}; + Gaudi::Property<double> m_hcal_time_window{this, "HcalTimeWindow", 1000.0, "The time window for HCAL in ns"}; + Gaudi::Property<double> m_muon_time_window{this, "MuonTimeWindow", 100.0, "The time window for MUON in ns"}; }; #endif diff --git a/Simulation/DetSimMixing/src/IBackgroundLoader.hh b/Simulation/DetSimMixing/src/IBackgroundLoader.hh old mode 100755 new mode 100644 index 273bb01a9264c1eebc78cf0b43a54c54f9d9da72..9142a9404ffe21f71937c67963cd375176014a2f --- a/Simulation/DetSimMixing/src/IBackgroundLoader.hh +++ b/Simulation/DetSimMixing/src/IBackgroundLoader.hh @@ -1,56 +1,7 @@ #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. - std::map<std::string, size_t> collection_subdet; // key is collection name, value is subdetector type. - - // in order to reduce the memory usage, need additional filter setup here. - // for example, for TPC, we need to load all the hits. - // however, for the tracker/calo, we only need to load the hits which are - // in an interested event time window. - enum SubDetType { - kVXD = 0, // 200 ns - kITK = 1, // 200 ns - kTPC = 2, // 34 us - kOTK = 3, // 1 us - kECAL = 4, // 150 ns - kHCAL = 5, // 1 us - kMUON = 6, // 1 us // todo - kNSubDetType - }; - std::vector<double> subdet2twindow = {200, 200, 34000, 1000, 150, 1000, 1000}; // key is subdet, value is time window. - std::vector<std::vector<std::string>> subdet2colnames = { - {"VXD"}, - {"ITK"}, - {"TPC"}, - {"OTK"}, - {"Ecal"}, - {"Hcal"}, - {"Muon"} - }; -}; +struct BackgroundEvent; class IBackgroundLoader { public: