diff --git a/Simulation/DetSimMixing/CMakeLists.txt b/Simulation/DetSimMixing/CMakeLists.txt
index bfaea3d388836406822388202c23122e72b439c2..2718d65aa24b3876cba39ba0483311ff7354bf67 100755
--- a/Simulation/DetSimMixing/CMakeLists.txt
+++ b/Simulation/DetSimMixing/CMakeLists.txt
@@ -4,6 +4,7 @@ gaudi_add_module(DetSimMixing
EDM4HEP::edm4hep EDM4HEP::edm4hepDict
podio::podioRootIO
Gaudi::GaudiKernel
+ k4FWCore::k4FWCore
)
install(TARGETS DetSimMixing
diff --git a/Simulation/DetSimMixing/src/BackgroundBatch.hh b/Simulation/DetSimMixing/src/BackgroundBatch.hh
index b8f334b2afafc2c06ad5da0c01d089c3ccf689b1..bdd16ba14e18dbca6bd5a1aafaf467a62ccad9ef 100755
--- a/Simulation/DetSimMixing/src/BackgroundBatch.hh
+++ b/Simulation/DetSimMixing/src/BackgroundBatch.hh
@@ -7,8 +7,8 @@
class BackgroundBatch {
public:
- TTimeStamp start_time; // the start time of the batch
- double duration; // the duration of the batch (in ns)
+ 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
// below are the metadata of each type of backgrouds in the batch
diff --git a/Simulation/DetSimMixing/src/BackgroundLoader.hh b/Simulation/DetSimMixing/src/BackgroundLoader.hh
index 792ac1e7195f8fb843d1911f64a8609ba1bb575b..7ab42bf9c7128ed71d9378ce847bc98ab2e3b112 100755
--- a/Simulation/DetSimMixing/src/BackgroundLoader.hh
+++ b/Simulation/DetSimMixing/src/BackgroundLoader.hh
@@ -35,37 +35,109 @@ public:
// std::cout << "available collection: " << col_name << std::endl;
evt.collection_index[col_name] = idx;
++idx;
+
+ // find the corresponding subdetector type.
+ for (size_t subdet = 0; subdet < BackgroundEvent::kNSubDetType; ++subdet) {
+ for (auto key: evt.subdet2colnames[subdet]) {
+ if (col_name.find(key) != std::string::npos) {
+ evt.collection_subdet[col_name] = subdet;
+ break;
+ }
+ }
+ }
}
}
for (auto [name, colidx]: evt.collection_index) {
auto col_ = frame.get(name);
+ // check whether the collection has a valid subdetector type.
+ if (evt.collection_subdet.find(name) == evt.collection_subdet.end()) {
+ continue;
+ }
+
+ auto time_window = evt.subdet2twindow[evt.collection_subdet[name]];
+
+ // if the current time is out of the time window, then skip this collection.
+ if (current_time_in_ns < -time_window || current_time_in_ns > time_window) {
+ continue;
+ }
+
+ std::cout << "Collection: " << name
+ << ", index: " << colidx
+ << ", time window: " << time_window << " ns"
+ << ", current time: " << current_time_in_ns;
+
+ // debug only. don't create any hits.
+ // conclusion: no memory leakage in the above code.
+ // continue;
+
+ int counter = 0;
+
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);
+ auto t = oldhit.getTime() + current_time_in_ns;
+ if (t < -time_window || t > time_window) {
+ // if the hit is not in the time window, skip.
+ continue;
+ }
+
+ auto newhit = trk_col.create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEDep(oldhit.getEDep());
+ newhit.setTime(t); // new
+ newhit.setPathLength(oldhit.getPathLength());
+ newhit.setQuality(oldhit.getQuality());
+ newhit.setPosition(oldhit.getPosition());
+ newhit.setMomentum(oldhit.getMomentum());
+
+ // extra
newhit.setOverlay(true);
- trk_col.push_back(newhit);
+
+ ++counter;
}
} else if (auto col = dynamic_cast<const edm4hep::SimCalorimeterHitCollection*>(col_)) {
auto& calo_col = evt.calorimeter_hits[colidx];
+ auto& calo_contrib_col = evt.calo_contribs[colidx];
for (auto oldhit: *col) {
- auto newhit = oldhit.clone();
+ // check whether the hit is in the time window.
+ bool is_in_window = false;
+ for (auto contrib: oldhit.getContributions()) {
+ auto t = contrib.getTime() + current_time_in_ns;
+ if (t < -time_window || t > time_window) {
+ continue;
+ }
+ is_in_window = true;
+ break;
+ }
+ if (not is_in_window) {
+ continue;
+ }
+
+ auto newhit = calo_col.create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEnergy(oldhit.getEnergy());
+ newhit.setPosition(oldhit.getPosition());
+ // todo: fixme
// loop all the contributions and add the time.
for (auto contrib: oldhit.getContributions()) {
- auto newcontrib = contrib.clone();
+ auto newcontrib = calo_contrib_col.create();
+ newcontrib.setPDG(contrib.getPDG());
+ newcontrib.setEnergy(contrib.getEnergy());
+ newcontrib.setStepPosition(contrib.getStepPosition());
newcontrib.setTime(contrib.getTime() + current_time_in_ns);
newhit.addToContributions(newcontrib);
}
- calo_col.push_back(newhit);
+ ++counter;
}
} else {
continue;
}
+ std::cout << ", counter: " << counter << std::endl;
+
}
++m_current_evt;
diff --git a/Simulation/DetSimMixing/src/DetSimMixingAlg.cc b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc
index f054f3b059631c8accba8b3d655d83d1f957625c..23c9bcb56c4df945ef938afa23a5c3904891b076 100755
--- a/Simulation/DetSimMixing/src/DetSimMixingAlg.cc
+++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.cc
@@ -38,20 +38,58 @@ StatusCode DetSimMixingAlg::initialize() {
// prepare the loaders
for (size_t i = 0; i < m_event_types.size(); ++i) {
- m_total_rates += m_event_rates[i];
+ // only the positive rates are considered into total rates
+ if (m_event_rates[i] > 0) {
+ 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;
+ if (m_total_rates < 0) {
+ error() << "The total rate of the background events should be >= 0." << endmsg;
return StatusCode::FAILURE;
+ } else if (m_total_rates == 0) {
+ info() << "All background events will be loaded in fixed time window mode." << endmsg;
+ } else {
+ m_total_tau = 1.0 / m_total_rates;
}
-
- 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;
+ if (m_event_rates[i] > 0) {
+ info() << " Event type: " << m_event_types[i] << ", rate: " << m_event_rates[i] << " Hz" << endmsg;
+ } else {
+ info() << " Event type: " << m_event_types[i] << ", time window: " << fabs(m_event_rates[i]) << " ns" << endmsg;
+ }
+ }
+
+ // prepare for the output
+ // for simplicity, the collection names are with prefix. the key in the map is without prefix, so it is same to the original input.
+ std::string prefix = "Mixed";
+ // SimTrackerHitCollections
+ for (const auto& name : m_trackerColNames) {
+ std::string name_col = name + "Collection";
+ auto col = new DataHandle<edm4hep::SimTrackerHitCollection>(prefix + name_col, Gaudi::DataHandle::Writer, this);
+ m_trackerColMap[name_col] = col;
+
+ auto sig_col = new DataHandle<edm4hep::SimTrackerHitCollection>(name_col, Gaudi::DataHandle::Reader, this);
+ m_sig_trackerColMap[name_col] = col;
+
+ }
+
+ // SimCalorimeterHitCollections
+ for (const auto& name : m_calorimeterColNames) {
+ std::string name_col = name + "Collection";
+ std::string name_contrib_col = name + "ContributionCollection";
+ auto col = new DataHandle<edm4hep::SimCalorimeterHitCollection>(prefix + name_col, Gaudi::DataHandle::Writer, this);
+ m_calorimeterColMap[name_col] = col;
+ auto col_contrib = new DataHandle<edm4hep::CaloHitContributionCollection>(prefix + name_contrib_col, Gaudi::DataHandle::Writer, this);
+ m_caloContribColMap[name_col] = col_contrib;
+
+ auto sig_col = new DataHandle<edm4hep::SimCalorimeterHitCollection>(name_col, Gaudi::DataHandle::Reader, this);
+ m_sig_calorimeterColMap[name_col] = sig_col;
+ auto sig_col_contrib = new DataHandle<edm4hep::CaloHitContributionCollection>(name_contrib_col, Gaudi::DataHandle::Reader, this);
+ m_sig_caloContribColMap[name_col] = sig_col_contrib;
+
}
return sc;
@@ -62,18 +100,22 @@ StatusCode DetSimMixingAlg::execute() {
info() << "Execute DetSimMixingAlg... " << endmsg;
+ // ========================================================================
+ // Reset
+ // ========================================================================
+ double start_time{0.0}; // ns, the start time of the event, always align to the signal
+
// ========================================================================
// Prepare the batches
// ========================================================================
std::vector<BackgroundBatch> batches;
- int nbatches = 5;
- double duration = 1e3; // 1 us = 1,000 ns of each batch
- TTimeStamp start_time;
+ 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
- for (int i = 0; i < nbatches; i++) {
+ for (int i = -nbatches; i < nbatches; ++i) {
BackgroundBatch batch;
- batch.start_time = start_time + TTimeStamp(0, static_cast<int>(i*duration));
+ batch.start_time = start_time + i*duration;
batch.duration = duration;
// todo: need to sample according to the rate
@@ -84,12 +126,40 @@ StatusCode DetSimMixingAlg::execute() {
double current_time = 0;
- while (true) {
+ // =======================================================================
+ // fixed time window mode
+ // =======================================================================
+ // 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) {
+ continue;
+ }
+ double time_window_event = fabs(m_event_rates[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) {
+ double this_current_time = j * time_window_event; // ns
+ ++batch.num_events;
+ batch.event_types.push_back(evttype);
+ batch.event_times.push_back(this_current_time);
+ batch.event_numbers_groupby_type[evttype] += 1;
+ }
+ }
+
+ // =======================================================================
+ // sampled mode
+ // =======================================================================
+ while (m_total_tau>0) {
// 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) {
+ // skip the fixed time mode type
+ if (m_event_rates[evttype] <= 0) {
+ continue;
+ }
accumulated += m_event_rates[evttype];
if (r < accumulated) {
selected_evttype = evttype;
@@ -121,13 +191,15 @@ StatusCode DetSimMixingAlg::execute() {
// of BackgroundEvent.
// ========================================================================
BackgroundEvent bkg_evt;
-
+ info() << "Creating a BackgroundEvent..." << endmsg;
for (size_t i = 0; i < batches.size(); ++i) {
+ info() << " Batch " << i << ": " << endmsg;
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];
+ double evttime = batch.start_time + batch.event_times[j];
+ info() << " Event type: " << m_event_types[evttype] << ", time: " << evttime << " ns" << endmsg;
IBackgroundLoader* loader = m_event_loaders[evttype];
@@ -138,10 +210,11 @@ StatusCode DetSimMixingAlg::execute() {
}
}
-
+ info() << "BackgroundEvent created." << endmsg;
+ info() << "Summary: " << endmsg;
// dump the information
for (auto [name, idx]: bkg_evt.collection_index) {
- info() << "Collection: " << name << ", index: " << idx;
+ 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)) {
@@ -154,6 +227,99 @@ StatusCode DetSimMixingAlg::execute() {
// ========================================================================
// Put the BackgroundEvent into the event store
// ========================================================================
+ // add prefix 'Mixed'.
+ for (auto [col_name, colidx]: bkg_evt.collection_index) {
+
+ if (bkg_evt.tracker_hits.count(colidx)) {
+ auto& col = bkg_evt.tracker_hits[colidx];
+
+ auto newcol = m_trackerColMap[col_name]->createAndPut();
+
+ // background
+ for (auto oldhit: col) {
+ auto newhit = newcol->create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEDep(oldhit.getEDep());
+ newhit.setTime(oldhit.getTime()); // new
+ newhit.setPathLength(oldhit.getPathLength());
+ newhit.setQuality(oldhit.getQuality());
+ newhit.setPosition(oldhit.getPosition());
+ newhit.setMomentum(oldhit.getMomentum());
+ }
+
+ // signal
+ auto sig_col = m_sig_trackerColMap[col_name]->get();
+ if (!sig_col) {
+ continue;
+ }
+ for (auto oldhit: *sig_col) {
+ auto newhit = newcol->create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEDep(oldhit.getEDep());
+ newhit.setTime(oldhit.getTime()); // new
+ newhit.setPathLength(oldhit.getPathLength());
+ newhit.setQuality(oldhit.getQuality());
+ newhit.setPosition(oldhit.getPosition());
+ newhit.setMomentum(oldhit.getMomentum());
+ }
+
+ } else if (bkg_evt.calorimeter_hits.count(colidx)) {
+ auto& col = bkg_evt.calorimeter_hits[colidx];
+ auto& col_contrib = bkg_evt.calo_contribs[colidx];
+
+ auto newcol = m_calorimeterColMap[col_name]->createAndPut();
+ auto newcol_contrib = m_caloContribColMap[col_name]->createAndPut();
+
+ for (auto oldhit: col) {
+ auto newhit = newcol->create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEnergy(oldhit.getEnergy());
+ newhit.setPosition(oldhit.getPosition());
+ // todo: fixme
+ // loop all the contributions and add the time.
+ for (auto contrib: oldhit.getContributions()) {
+ auto newcontrib = newcol_contrib->create();
+ newcontrib.setPDG(contrib.getPDG());
+ newcontrib.setEnergy(contrib.getEnergy());
+ newcontrib.setStepPosition(contrib.getStepPosition());
+ newcontrib.setTime(contrib.getTime());
+ newhit.addToContributions(newcontrib);
+ }
+ }
+
+ // signal
+ auto sig_col = m_sig_calorimeterColMap[col_name]->get();
+ auto sig_col_contrib = m_sig_caloContribColMap[col_name]->get();
+
+ if (!sig_col || !sig_col_contrib) {
+ continue;
+ }
+
+ for (auto oldhit: *sig_col) {
+ auto newhit = newcol->create();
+ newhit.setCellID(oldhit.getCellID());
+ newhit.setEnergy(oldhit.getEnergy());
+ newhit.setPosition(oldhit.getPosition());
+ // todo: fixme
+ // loop all the contributions and add the time.
+ for (auto contrib: oldhit.getContributions()) {
+ auto newcontrib = newcol_contrib->create();
+ newcontrib.setPDG(contrib.getPDG());
+ newcontrib.setEnergy(contrib.getEnergy());
+ newcontrib.setStepPosition(contrib.getStepPosition());
+ newcontrib.setTime(contrib.getTime());
+ newhit.addToContributions(newcontrib);
+ }
+ }
+
+ } else {
+ error() << "The collection " << col_name
+ << " with idx" << colidx
+ << " is not found." << endmsg;
+ continue;
+ }
+ }
+
return sc;
}
@@ -164,4 +330,4 @@ StatusCode DetSimMixingAlg::finalize() {
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
index d678dca31c93e7c989d8b32ea3129fa5fc5f100b..55aa56a7def63ac36ddd554ca49d161a81b41171 100755
--- a/Simulation/DetSimMixing/src/DetSimMixingAlg.hh
+++ b/Simulation/DetSimMixing/src/DetSimMixingAlg.hh
@@ -59,7 +59,7 @@ public:
private:
double m_total_rates{0}; // Hz.
- double m_total_tau{0}; // s. the total time of the background events
+ double m_total_tau{0}; // s.
private:
std::vector<std::string> m_event_types; // the event types of the background events
@@ -67,9 +67,38 @@ private:
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
+ // following maps are used to load signal collections from memory
+ std::map<std::string, DataHandle<edm4hep::SimTrackerHitCollection>*> m_sig_trackerColMap;
+ std::map<std::string, DataHandle<edm4hep::SimCalorimeterHitCollection>*> m_sig_calorimeterColMap;
+ std::map<std::string, DataHandle<edm4hep::CaloHitContributionCollection>*> m_sig_caloContribColMap;
+
+private:
+ // following maps are used to put collections into memory
+ std::map<std::string, DataHandle<edm4hep::SimTrackerHitCollection>*> m_trackerColMap;
+ std::map<std::string, DataHandle<edm4hep::SimCalorimeterHitCollection>*> m_calorimeterColMap;
+ std::map<std::string, DataHandle<edm4hep::CaloHitContributionCollection>*> m_caloContribColMap;
+
+ // the name here is without suffix "Collection"
+ Gaudi::Property<std::vector<std::string>> m_trackerColNames{this,
+ "TrackerCollections",
+ {"VXD", "ITKBarrel", "ITKEndcap", "TPC", "TPCLowPt", "TPCSpacePoint",
+ "OTKBarrel", "OTKEndcap", "MuonBarrel", "MuonEndcap"},
+ "Names of the Tracker collections (without suffix Collection)"};
+ Gaudi::Property<std::vector<std::string>> m_calorimeterColNames{this,
+ "CalorimeterCollections",
+ {"EcalBarrel", "EcalEndcaps", "EcalEndcapRing",
+ "HcalBarrel", "HcalEndcaps", "HcalEndcapRing"},
+ "Names of the Calorimeter collections (without suffix Collection)"};
+
private:
- // properties for user side
- Gaudi::Property<std::map<std::string, double>> m_background_rates{this, "BackgroundRates", {}, "The rates of the background events"};
+ // Following properties are used to configure different types of background events.
+ // * Key: label of one type
+ // * Value:
+ // * rate:
+ // * 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, std::vector<std::string>>> m_background_filelists{this, "BackgroundFileLists", {}, "The input file lists for the background events"};
};
diff --git a/Simulation/DetSimMixing/src/IBackgroundLoader.hh b/Simulation/DetSimMixing/src/IBackgroundLoader.hh
index 588d928b0d8c3fb3910fbde501e588ce5f7b96ac..273bb01a9264c1eebc78cf0b43a54c54f9d9da72 100755
--- a/Simulation/DetSimMixing/src/IBackgroundLoader.hh
+++ b/Simulation/DetSimMixing/src/IBackgroundLoader.hh
@@ -24,6 +24,32 @@ struct BackgroundEvent {
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"}
+ };
};
class IBackgroundLoader {