10 #include <background/modules/BeamBkgMixer/BeamBkgMixerModule.h>
15 #include <framework/datastore/DataStore.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/datastore/StoreObjPtr.h>
18 #include <framework/datastore/RelationArray.h>
21 #include <framework/core/ModuleParam.templateDetails.h>
22 #include <framework/logging/Logger.h>
25 #include <pxd/dataobjects/PXDSimHit.h>
26 #include <svd/dataobjects/SVDSimHit.h>
27 #include <cdc/dataobjects/CDCSimHit.h>
28 #include <top/dataobjects/TOPSimHit.h>
29 #include <arich/dataobjects/ARICHSimHit.h>
30 #include <ecl/dataobjects/ECLHit.h>
31 #include <klm/dataobjects/bklm/BKLMSimHit.h>
32 #include <klm/dataobjects/eklm/EKLMSimHit.h>
33 #include <simulation/dataobjects/BeamBackHit.h>
36 #include <framework/dataobjects/EventMetaData.h>
37 #include <framework/dataobjects/BackgroundInfo.h>
67 setDescription(
"Beam background mixer at SimHit level that uses beam background"
68 " simulation output directly (collision files) and not ROF files. "
69 "Each background event is shifted in time randomly within "
70 "a time window specified with minTime and maxTime.");
73 addParam(
"backgroundFiles", m_backgroundFiles,
74 "List of background (collision) files (wildcards not allowed - "
75 "use python glob.glob() to expand to list of files)");
76 addParam(
"minTime", m_minTime,
77 "Time window lower edge in nano seconds", -1000.0);
78 addParam(
"maxTime", m_maxTime,
79 "Time window upper edge in nano seconds", 800.0);
80 addParam(
"overallScaleFactor", m_overallScaleFactor,
81 "Overall factor to scale rates of backgrounds", 1.0);
82 addParam(
"scaleFactors", m_scaleFactors,
83 "Factors to scale rates of backgrounds. "
84 "Possible tag names: " + m_bgTypes.getBGTypes(),
86 addParam(
"components", m_components,
87 "Detector components to be included, empty list means all components",
89 addParam(
"wrapAround", m_wrapAround,
90 "if true wrap around events passing time window upper edge",
true);
91 addParam(
"minTimeECL", m_minTimeECL,
92 "Time window lower edge for ECL in nano seconds", -17600.0);
93 addParam(
"maxTimeECL", m_maxTimeECL,
94 "Time window upper edge for ECL in nano seconds", 8500.0);
95 addParam(
"minTimePXD", m_minTimePXD,
96 "Time window lower edge for PXD in nano seconds", -10000.0);
97 addParam(
"maxTimePXD", m_maxTimePXD,
98 "Time window upper edge for PXD in nano seconds", 10000.0);
100 addParam(
"beambackhits", m_BeamBackHits,
101 "If true also add the BeamBackHits collection for the selected "
102 "subdetectors to the output file",
false);
104 addParam(
"maxEdepECL", m_maxEdepECL,
105 "maximal deposited energy of ECLHit to accept BG event for mixing"
106 "(0 means accept all events)", 1.0);
108 addParam(
"cacheSize", m_cacheSize,
109 "file cache size in Mbytes. If negative, use root default", 0);
112 BeamBkgMixerModule::~BeamBkgMixerModule()
116 void BeamBkgMixerModule::initialize()
121 std::vector<std::string> components = m_components;
122 m_PXD = isComponentIncluded(components,
"PXD");
123 m_SVD = isComponentIncluded(components,
"SVD");
124 m_CDC = isComponentIncluded(components,
"CDC");
125 m_TOP = isComponentIncluded(components,
"TOP");
126 m_ARICH = isComponentIncluded(components,
"ARICH");
127 m_ECL = isComponentIncluded(components,
"ECL");
128 m_KLM = isComponentIncluded(components,
"KLM");
132 isComponentIncluded(components,
"MagneticField2d");
133 isComponentIncluded(components,
"MagneticField3d");
134 isComponentIncluded(components,
"MagneticField");
135 isComponentIncluded(components,
"MagneticFieldConstant4LimitedRCDC");
136 isComponentIncluded(components,
"MagneticFieldConstant4LimitedRSVD");
137 isComponentIncluded(components,
"BeamPipe");
138 isComponentIncluded(components,
"Cryostat");
139 isComponentIncluded(components,
"FarBeamLine");
140 isComponentIncluded(components,
"HeavyMetalShield");
141 isComponentIncluded(components,
"COIL");
142 isComponentIncluded(components,
"STR");
143 isComponentIncluded(components,
"VXDService");
145 if (!components.empty()) {
147 for (
unsigned i = 0; i < components.size(); ++i) str = str +
" " + components[i];
148 B2WARNING(
"Unknown components:" << str);
153 for (
auto file : m_backgroundFiles) {
156 if (TString(file.c_str()).Contains(
"*")) {
157 B2ERROR(file <<
": wildcards are not allowed");
162 TFile* f = TFile::Open(file.c_str(),
"READ");
164 B2ERROR(file <<
": file not found");
168 B2ERROR(file <<
": can't open file");
173 TChain persistent(
"persistent");
174 int nFiles = persistent.Add(file.c_str());
176 B2ERROR(file <<
": no such files");
179 if (persistent.GetEntries() == 0) {
180 B2ERROR(file <<
": tree 'persistent' has no entries");
184 TObject* bkgMetaData = 0;
185 TBranch* branchBMD = persistent.GetBranch(
"BackgroundMetaData");
187 B2ERROR(file <<
": branch 'BackgroundMetaData' not found");
190 branchBMD->SetAddress(&bkgMetaData);
192 std::vector<BackgroundMetaData::BG_TAG> tags;
193 std::vector<std::string> types;
194 std::vector<BackgroundMetaData::EFileType> fileTypes;
196 for (
unsigned k = 0; k < persistent.GetEntries(); k++) {
197 persistent.GetEntry(k);
205 B2ERROR(file <<
": invalid realTime: " << realTime);
208 for (
unsigned i = 1; i < tags.size(); ++i) {
209 if (tags[i] != tags[0]) {
210 B2ERROR(file <<
": files with mixed background types not supported");
214 for (
unsigned i = 1; i < fileTypes.size(); ++i) {
215 if (fileTypes[i] != fileTypes[0]) {
216 B2ERROR(file <<
": files with mixed file types not supported");
221 appendSample(tags[0], types[0], file, realTime, fileTypes[0]);
228 for (
auto scaleFactor : m_scaleFactors) {
229 std::string type = std::get<0>(scaleFactor);
230 if (m_bgTypes.getTag(type) == 0)
231 B2ERROR(
"Unknown beam background type found in 'scaleFactors': " << type <<
"\n"
232 "Possible are: " + m_bgTypes.getBGTypes());
233 for (
auto& bkg : m_backgrounds) {
234 if (bkg.type.find(type) != std::string::npos)
235 bkg.scaleFactor *= std::get<1>(scaleFactor);
241 for (
auto& bkg : m_backgrounds) {
244 bkg.tree.reset(
new TChain(
"tree"));
245 for (
unsigned i = 0; i < bkg.fileNames.size(); ++i) {
246 bkg.numFiles += bkg.tree->Add(bkg.fileNames[i].c_str());
249 bkg.numEvents = bkg.tree->GetEntries();
250 bkg.rate = bkg.numEvents / bkg.realTime * bkg.scaleFactor;
252 if (m_cacheSize >= 0) bkg.tree->SetCacheSize(m_cacheSize * 1024 * 1024);
254 if (m_PXD and bkg.tree->GetBranch(
"PXDSimHits"))
255 bkg.tree->SetBranchAddress(
"PXDSimHits", &m_simHits.PXD);
256 if (m_SVD and bkg.tree->GetBranch(
"SVDSimHits"))
257 bkg.tree->SetBranchAddress(
"SVDSimHits", &m_simHits.SVD);
258 if (m_CDC and bkg.tree->GetBranch(
"CDCSimHits"))
259 bkg.tree->SetBranchAddress(
"CDCSimHits", &m_simHits.CDC);
260 if (m_TOP and bkg.tree->GetBranch(
"TOPSimHits"))
261 bkg.tree->SetBranchAddress(
"TOPSimHits", &m_simHits.TOP);
262 if (m_ARICH and bkg.tree->GetBranch(
"ARICHSimHits"))
263 bkg.tree->SetBranchAddress(
"ARICHSimHits", &m_simHits.ARICH);
264 if (m_ECL and bkg.tree->GetBranch(
"ECLHits"))
265 bkg.tree->SetBranchAddress(
"ECLHits", &m_simHits.ECL);
266 if (m_KLM and bkg.tree->GetBranch(
"BKLMSimHits"))
267 bkg.tree->SetBranchAddress(
"BKLMSimHits", &m_simHits.BKLM);
268 if (m_KLM and bkg.tree->GetBranch(
"EKLMSimHits"))
269 bkg.tree->SetBranchAddress(
"EKLMSimHits", &m_simHits.EKLM);
271 if (m_BeamBackHits and bkg.tree->GetBranch(
"BeamBackHits"))
272 bkg.tree->SetBranchAddress(
"BeamBackHits", &m_simHits.BeamBackHits);
275 std::string unit(
" ns");
276 double realTime = bkg.realTime;
277 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" us";}
278 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" ms";}
279 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" s";}
281 B2INFO(
"BeamBkgMixer: " << bkg.type <<
282 " files=" << bkg.numFiles <<
283 " events=" << bkg.numEvents <<
284 " realTime=" << realTime << unit <<
285 " scaleFactor=" << bkg.scaleFactor <<
286 " rate=" << bkg.rate * 1000 <<
" MHz");
325 bkgInfo->setMethod(BackgroundInfo::c_Mixing);
326 bkgInfo->setComponents(m_components);
327 bkgInfo->setMinTime(m_minTime);
328 bkgInfo->setMaxTime(m_maxTime);
329 bkgInfo->setMinTimeECL(m_minTimeECL);
330 bkgInfo->setMaxTimeECL(m_maxTimeECL);
331 bkgInfo->setMinTimePXD(m_minTimePXD);
332 bkgInfo->setMaxTimePXD(m_maxTimePXD);
333 bkgInfo->setWrapAround(m_wrapAround);
334 bkgInfo->setMaxEdepECL(m_maxEdepECL);
335 for (
auto& bkg : m_backgrounds) {
338 descr.type = bkg.type;
339 descr.fileType = bkg.fileType;
340 descr.fileNames = bkg.fileNames;
341 descr.realTime = bkg.realTime;
342 descr.numEvents = bkg.numEvents;
343 descr.scaleFactor = bkg.scaleFactor;
344 descr.rate = bkg.rate;
346 bkgInfo->appendBackgroundDescr(descr);
351 void BeamBkgMixerModule::beginRun()
355 void BeamBkgMixerModule::event()
368 for (
auto& bkg : m_backgrounds) {
370 if (bkg.fileType != BackgroundMetaData::c_Usual)
continue;
372 double mean = bkg.rate * (m_maxTime - m_minTime);
373 int nev = gRandom->Poisson(mean);
375 for (
int iev = 0; iev < nev; iev++) {
376 double timeShift = gRandom->Rndm() * (m_maxTime - m_minTime) + m_minTime;
377 bkg.tree->GetEntry(bkg.eventCount);
379 if (acceptEvent(m_simHits.ECL)) {
380 addSimHits(pxdSimHits, m_simHits.PXD, timeShift, m_minTime, m_maxTime);
381 addSimHits(svdSimHits, m_simHits.SVD, timeShift, m_minTime, m_maxTime);
382 addSimHits(cdcSimHits, m_simHits.CDC, timeShift, m_minTime, m_maxTime);
383 addSimHits(topSimHits, m_simHits.TOP, timeShift, m_minTime, m_maxTime);
384 addSimHits(arichSimHits, m_simHits.ARICH, timeShift, m_minTime, m_maxTime);
385 addSimHits(eclHits, m_simHits.ECL, timeShift, m_minTime, m_maxTime);
386 addSimHits(bklmSimHits, m_simHits.BKLM, timeShift, m_minTime, m_maxTime);
387 addSimHits(eklmSimHits, m_simHits.EKLM, timeShift, m_minTime, m_maxTime);
388 addBeamBackHits(beamBackHits, m_simHits.BeamBackHits, timeShift,
389 m_minTime, m_maxTime);
392 std::string message =
"BeamBkgMixer: event " + to_string(bkg.eventCount)
393 +
" of " + bkg.type +
" rejected due to large energy deposit in ECL";
394 m_rejected[message] += 1;
396 if (m_rejectedCount < 10) {
397 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL");
398 }
else if (m_rejectedCount == 10) {
399 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL "
400 <<
"(message will be suppressed now)");
405 if (bkg.eventCount >= bkg.numEvents) {
407 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
408 m_reused[message] += 1;
409 if (m_reused[message] == 1) B2INFO(message);
410 bkgInfo->incrementReusedCounter(bkg.index);
416 for (
auto& bkg : m_backgrounds) {
418 if (bkg.fileType != BackgroundMetaData::c_ECL)
continue;
420 double mean = bkg.rate * (m_maxTimeECL - m_minTimeECL);
421 int nev = gRandom->Poisson(mean);
423 for (
int iev = 0; iev < nev; iev++) {
424 double timeShift = gRandom->Rndm() * (m_maxTimeECL - m_minTimeECL) + m_minTimeECL;
425 if (timeShift > m_minTime and timeShift < m_maxTime)
continue;
426 bkg.tree->GetEntry(bkg.eventCount);
428 if (acceptEvent(m_simHits.ECL)) {
429 double minTime = m_minTimeECL;
430 double maxTime = m_maxTimeECL;
431 if (timeShift <= m_minTime) {
436 addSimHits(eclHits, m_simHits.ECL, timeShift, minTime, maxTime);
439 std::string message =
"BeamBkgMixer: event " + to_string(bkg.eventCount)
440 +
" of " + bkg.type +
" rejected due to large energy deposit in ECL";
441 m_rejected[message] += 1;
443 if (m_rejectedCount < 10) {
444 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL");
445 }
else if (m_rejectedCount == 10) {
446 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL "
447 <<
"(message will be suppressed now)");
452 if (bkg.eventCount >= bkg.numEvents) {
454 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
455 m_reused[message] += 1;
456 if (m_reused[message] == 1) B2INFO(message);
457 bkgInfo->incrementReusedCounter(bkg.index);
464 for (
auto& bkg : m_backgrounds) {
466 if (bkg.fileType != BackgroundMetaData::c_PXD)
continue;
468 double mean = bkg.rate * (m_maxTimePXD - m_minTimePXD);
469 int nev = gRandom->Poisson(mean);
471 for (
int iev = 0; iev < nev; iev++) {
472 double timeShift = gRandom->Rndm() * (m_maxTimePXD - m_minTimePXD) + m_minTimePXD;
473 if (timeShift > m_minTime and timeShift < m_maxTime)
continue;
474 bkg.tree->GetEntry(bkg.eventCount);
476 double minTime = m_minTimePXD;
477 double maxTime = m_maxTimePXD;
478 if (timeShift <= m_minTime) {
483 addSimHits(pxdSimHits, m_simHits.PXD, timeShift, minTime, maxTime);
486 if (bkg.eventCount >= bkg.numEvents) {
488 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
489 m_reused[message] += 1;
490 if (m_reused[message] == 1) B2INFO(message);
491 bkgInfo->incrementReusedCounter(bkg.index);
500 void BeamBkgMixerModule::endRun()
504 void BeamBkgMixerModule::terminate()
507 B2INFO(
"BeamBkgMixer - reused samples:");
508 for (
const auto& message : m_reused) {
509 B2INFO(message.first <<
"(occured " << message.second <<
" times)");
511 B2INFO(
"BeamBkgMixer - rejected events:");
512 for (
const auto& message : m_rejected) {
513 B2INFO(message.first <<
"(occured " << message.second <<
" times)");
516 for (
auto& bkg : m_backgrounds) {
523 bool BeamBkgMixerModule::isComponentIncluded(std::vector<std::string>& components,
524 const std::string& component)
526 if (m_components.empty())
return true;
527 auto iterator = std::find(components.begin(), components.end(), component);
528 if (iterator != components.end()) {
529 components.erase(iterator);
537 const std::string& type,
538 const std::string& fileName,
542 for (
auto& bkg : m_backgrounds) {
543 if (tag == bkg.tag and fileType == bkg.fileType) {
544 bkg.fileNames.push_back(fileName);
545 bkg.realTime += realTime;
549 std::string ftype = type;
550 if (fileType == BackgroundMetaData::c_ECL) ftype +=
"(ECL)";
551 if (fileType == BackgroundMetaData::c_PXD) ftype +=
"(PXD)";
552 unsigned index = m_backgrounds.size();
553 m_backgrounds.push_back(
BkgFiles(tag, ftype, fileName, realTime,
554 m_overallScaleFactor, fileType, index));
558 bool BeamBkgMixerModule::acceptEvent(TClonesArray* cloneArrayECL)
560 if (!cloneArrayECL)
return true;
561 if (m_maxEdepECL == 0)
return true;
563 int numEntries = cloneArrayECL->GetEntriesFast();
564 for (
int i = 0; i < numEntries; i++) {
565 ECLHit* simHit =
static_cast<ECLHit*
>(cloneArrayECL->AddrAt(i));
566 if (simHit->
getEnergyDep() > m_maxEdepECL)
return false;
New beam background mixer; this one doesn't need ROF files.
Class to store simulated hits which equate to average of ECLSImHit on crystals input for digitization...
double getEnergyDep() const
Get deposit energy.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
bool create(bool replace=false)
Create a default object in the data store.
Accessor to arrays stored in the data store.
Type-safe access to single objects in the data store.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Abstract base class for different kinds of events.
Structure for background description.
structure to hold samples of a particular background type