12 #include <background/modules/BeamBkgMixer/BeamBkgMixerModule.h>
17 #include <framework/datastore/DataStore.h>
18 #include <framework/datastore/StoreArray.h>
19 #include <framework/datastore/StoreObjPtr.h>
20 #include <framework/datastore/RelationArray.h>
23 #include <framework/core/ModuleParam.templateDetails.h>
24 #include <framework/logging/Logger.h>
27 #include <pxd/dataobjects/PXDSimHit.h>
28 #include <svd/dataobjects/SVDSimHit.h>
29 #include <cdc/dataobjects/CDCSimHit.h>
30 #include <top/dataobjects/TOPSimHit.h>
31 #include <arich/dataobjects/ARICHSimHit.h>
32 #include <ecl/dataobjects/ECLHit.h>
33 #include <klm/dataobjects/bklm/BKLMSimHit.h>
34 #include <klm/dataobjects/eklm/EKLMSimHit.h>
35 #include <simulation/dataobjects/BeamBackHit.h>
38 #include <framework/dataobjects/EventMetaData.h>
39 #include <framework/dataobjects/BackgroundInfo.h>
69 setDescription(
"Beam background mixer at SimHit level that uses beam background"
70 " simulation output directly (collision files) and not ROF files. "
71 "Each background event is shifted in time randomly within "
72 "a time window specified with minTime and maxTime.");
75 addParam(
"backgroundFiles", m_backgroundFiles,
76 "List of background (collision) files (wildcards not allowed - "
77 "use python glob.glob() to expand to list of files)");
78 addParam(
"minTime", m_minTime,
79 "Time window lower edge in nano seconds", -1000.0);
80 addParam(
"maxTime", m_maxTime,
81 "Time window upper edge in nano seconds", 800.0);
82 addParam(
"overallScaleFactor", m_overallScaleFactor,
83 "Overall factor to scale rates of backgrounds", 1.0);
84 addParam(
"scaleFactors", m_scaleFactors,
85 "Factors to scale rates of backgrounds. "
86 "Possible tag names: " + m_bgTypes.getBGTypes(),
88 addParam(
"components", m_components,
89 "Detector components to be included, empty list means all components",
91 addParam(
"wrapAround", m_wrapAround,
92 "if true wrap around events passing time window upper edge",
true);
93 addParam(
"minTimeECL", m_minTimeECL,
94 "Time window lower edge for ECL in nano seconds", -17600.0);
95 addParam(
"maxTimeECL", m_maxTimeECL,
96 "Time window upper edge for ECL in nano seconds", 8500.0);
97 addParam(
"minTimePXD", m_minTimePXD,
98 "Time window lower edge for PXD in nano seconds", -10000.0);
99 addParam(
"maxTimePXD", m_maxTimePXD,
100 "Time window upper edge for PXD in nano seconds", 10000.0);
102 addParam(
"beambackhits", m_BeamBackHits,
103 "If true also add the BeamBackHits collection for the selected "
104 "subdetectors to the output file",
false);
106 addParam(
"maxEdepECL", m_maxEdepECL,
107 "maximal deposited energy of ECLHit to accept BG event for mixing"
108 "(0 means accept all events)", 1.0);
110 addParam(
"cacheSize", m_cacheSize,
111 "file cache size in Mbytes. If negative, use root default", 0);
114 BeamBkgMixerModule::~BeamBkgMixerModule()
118 void BeamBkgMixerModule::initialize()
123 std::vector<std::string> components = m_components;
124 m_PXD = isComponentIncluded(components,
"PXD");
125 m_SVD = isComponentIncluded(components,
"SVD");
126 m_CDC = isComponentIncluded(components,
"CDC");
127 m_TOP = isComponentIncluded(components,
"TOP");
128 m_ARICH = isComponentIncluded(components,
"ARICH");
129 m_ECL = isComponentIncluded(components,
"ECL");
130 m_KLM = isComponentIncluded(components,
"KLM");
134 isComponentIncluded(components,
"MagneticField2d");
135 isComponentIncluded(components,
"MagneticField3d");
136 isComponentIncluded(components,
"MagneticField");
137 isComponentIncluded(components,
"MagneticFieldConstant4LimitedRCDC");
138 isComponentIncluded(components,
"MagneticFieldConstant4LimitedRSVD");
139 isComponentIncluded(components,
"BeamPipe");
140 isComponentIncluded(components,
"Cryostat");
141 isComponentIncluded(components,
"FarBeamLine");
142 isComponentIncluded(components,
"HeavyMetalShield");
143 isComponentIncluded(components,
"COIL");
144 isComponentIncluded(components,
"STR");
145 isComponentIncluded(components,
"VXDService");
147 if (!components.empty()) {
149 for (
unsigned i = 0; i < components.size(); ++i) str = str +
" " + components[i];
150 B2WARNING(
"Unknown components:" << str);
155 for (
auto file : m_backgroundFiles) {
158 if (TString(file.c_str()).Contains(
"*")) {
159 B2ERROR(file <<
": wildcards are not allowed");
164 TFile* f = TFile::Open(file.c_str(),
"READ");
166 B2ERROR(file <<
": file not found");
170 B2ERROR(file <<
": can't open file");
175 TChain persistent(
"persistent");
176 int nFiles = persistent.Add(file.c_str());
178 B2ERROR(file <<
": no such files");
181 if (persistent.GetEntries() == 0) {
182 B2ERROR(file <<
": tree 'persistent' has no entries");
186 TObject* bkgMetaData = 0;
187 TBranch* branchBMD = persistent.GetBranch(
"BackgroundMetaData");
189 B2ERROR(file <<
": branch 'BackgroundMetaData' not found");
192 branchBMD->SetAddress(&bkgMetaData);
194 std::vector<BackgroundMetaData::BG_TAG> tags;
195 std::vector<std::string> types;
196 std::vector<BackgroundMetaData::EFileType> fileTypes;
198 for (
unsigned k = 0; k < persistent.GetEntries(); k++) {
199 persistent.GetEntry(k);
207 B2ERROR(file <<
": invalid realTime: " << realTime);
210 for (
unsigned i = 1; i < tags.size(); ++i) {
211 if (tags[i] != tags[0]) {
212 B2ERROR(file <<
": files with mixed background types not supported");
216 for (
unsigned i = 1; i < fileTypes.size(); ++i) {
217 if (fileTypes[i] != fileTypes[0]) {
218 B2ERROR(file <<
": files with mixed file types not supported");
223 appendSample(tags[0], types[0], file, realTime, fileTypes[0]);
230 for (
auto scaleFactor : m_scaleFactors) {
231 std::string type = std::get<0>(scaleFactor);
232 if (m_bgTypes.getTag(type) == 0)
233 B2ERROR(
"Unknown beam background type found in 'scaleFactors': " << type <<
"\n"
234 "Possible are: " + m_bgTypes.getBGTypes());
235 for (
auto& bkg : m_backgrounds) {
236 if (bkg.type.find(type) != std::string::npos)
237 bkg.scaleFactor *= std::get<1>(scaleFactor);
243 for (
auto& bkg : m_backgrounds) {
246 bkg.tree.reset(
new TChain(
"tree"));
247 for (
unsigned i = 0; i < bkg.fileNames.size(); ++i) {
248 bkg.numFiles += bkg.tree->Add(bkg.fileNames[i].c_str());
251 bkg.numEvents = bkg.tree->GetEntries();
252 bkg.rate = bkg.numEvents / bkg.realTime * bkg.scaleFactor;
254 if (m_cacheSize >= 0) bkg.tree->SetCacheSize(m_cacheSize * 1024 * 1024);
256 if (m_PXD and bkg.tree->GetBranch(
"PXDSimHits"))
257 bkg.tree->SetBranchAddress(
"PXDSimHits", &m_simHits.PXD);
258 if (m_SVD and bkg.tree->GetBranch(
"SVDSimHits"))
259 bkg.tree->SetBranchAddress(
"SVDSimHits", &m_simHits.SVD);
260 if (m_CDC and bkg.tree->GetBranch(
"CDCSimHits"))
261 bkg.tree->SetBranchAddress(
"CDCSimHits", &m_simHits.CDC);
262 if (m_TOP and bkg.tree->GetBranch(
"TOPSimHits"))
263 bkg.tree->SetBranchAddress(
"TOPSimHits", &m_simHits.TOP);
264 if (m_ARICH and bkg.tree->GetBranch(
"ARICHSimHits"))
265 bkg.tree->SetBranchAddress(
"ARICHSimHits", &m_simHits.ARICH);
266 if (m_ECL and bkg.tree->GetBranch(
"ECLHits"))
267 bkg.tree->SetBranchAddress(
"ECLHits", &m_simHits.ECL);
268 if (m_KLM and bkg.tree->GetBranch(
"BKLMSimHits"))
269 bkg.tree->SetBranchAddress(
"BKLMSimHits", &m_simHits.BKLM);
270 if (m_KLM and bkg.tree->GetBranch(
"EKLMSimHits"))
271 bkg.tree->SetBranchAddress(
"EKLMSimHits", &m_simHits.EKLM);
273 if (m_BeamBackHits and bkg.tree->GetBranch(
"BeamBackHits"))
274 bkg.tree->SetBranchAddress(
"BeamBackHits", &m_simHits.BeamBackHits);
277 std::string unit(
" ns");
278 double realTime = bkg.realTime;
279 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" us";}
280 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" ms";}
281 if (realTime >= 1000.0) {realTime /= 1000.0; unit =
" s";}
283 B2INFO(
"BeamBkgMixer: " << bkg.type <<
284 " files=" << bkg.numFiles <<
285 " events=" << bkg.numEvents <<
286 " realTime=" << realTime << unit <<
287 " scaleFactor=" << bkg.scaleFactor <<
288 " rate=" << bkg.rate * 1000 <<
" MHz");
295 if (m_PXD) pxdSimHits.registerInDataStore();
298 if (m_SVD) svdSimHits.registerInDataStore();
301 if (m_CDC) cdcSimHits.registerInDataStore();
304 if (m_TOP) topSimHits.registerInDataStore();
307 if (m_ARICH) arichSimHits.registerInDataStore();
310 if (m_ECL) eclHits.registerInDataStore();
313 if (m_KLM) bklmSimHits.registerInDataStore();
316 if (m_KLM) eklmSimHits.registerInDataStore();
319 if (m_BeamBackHits) beamBackHits.registerInDataStore();
325 bkgInfo.registerInDataStore();
327 bkgInfo->setMethod(BackgroundInfo::c_Mixing);
328 bkgInfo->setComponents(m_components);
329 bkgInfo->setMinTime(m_minTime);
330 bkgInfo->setMaxTime(m_maxTime);
331 bkgInfo->setMinTimeECL(m_minTimeECL);
332 bkgInfo->setMaxTimeECL(m_maxTimeECL);
333 bkgInfo->setMinTimePXD(m_minTimePXD);
334 bkgInfo->setMaxTimePXD(m_maxTimePXD);
335 bkgInfo->setWrapAround(m_wrapAround);
336 bkgInfo->setMaxEdepECL(m_maxEdepECL);
337 for (
auto& bkg : m_backgrounds) {
340 descr.type = bkg.type;
341 descr.fileType = bkg.fileType;
342 descr.fileNames = bkg.fileNames;
343 descr.realTime = bkg.realTime;
344 descr.numEvents = bkg.numEvents;
345 descr.scaleFactor = bkg.scaleFactor;
346 descr.rate = bkg.rate;
348 bkgInfo->appendBackgroundDescr(descr);
353 void BeamBkgMixerModule::beginRun()
357 void BeamBkgMixerModule::event()
370 for (
auto& bkg : m_backgrounds) {
372 if (bkg.fileType != BackgroundMetaData::c_Usual)
continue;
374 double mean = bkg.rate * (m_maxTime - m_minTime);
375 int nev = gRandom->Poisson(mean);
377 for (
int iev = 0; iev < nev; iev++) {
378 double timeShift = gRandom->Rndm() * (m_maxTime - m_minTime) + m_minTime;
379 bkg.tree->GetEntry(bkg.eventCount);
381 if (acceptEvent(m_simHits.ECL)) {
382 addSimHits(pxdSimHits, m_simHits.PXD, timeShift, m_minTime, m_maxTime);
383 addSimHits(svdSimHits, m_simHits.SVD, timeShift, m_minTime, m_maxTime);
384 addSimHits(cdcSimHits, m_simHits.CDC, timeShift, m_minTime, m_maxTime);
385 addSimHits(topSimHits, m_simHits.TOP, timeShift, m_minTime, m_maxTime);
386 addSimHits(arichSimHits, m_simHits.ARICH, timeShift, m_minTime, m_maxTime);
387 addSimHits(eclHits, m_simHits.ECL, timeShift, m_minTime, m_maxTime);
388 addSimHits(bklmSimHits, m_simHits.BKLM, timeShift, m_minTime, m_maxTime);
389 addSimHits(eklmSimHits, m_simHits.EKLM, timeShift, m_minTime, m_maxTime);
390 addBeamBackHits(beamBackHits, m_simHits.BeamBackHits, timeShift,
391 m_minTime, m_maxTime);
394 std::string message =
"BeamBkgMixer: event " + to_string(bkg.eventCount)
395 +
" of " + bkg.type +
" rejected due to large energy deposit in ECL";
396 m_rejected[message] += 1;
398 if (m_rejectedCount < 10) {
399 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL");
400 }
else if (m_rejectedCount == 10) {
401 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL "
402 <<
"(message will be suppressed now)");
407 if (bkg.eventCount >= bkg.numEvents) {
409 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
410 m_reused[message] += 1;
411 if (m_reused[message] == 1) B2INFO(message);
412 bkgInfo->incrementReusedCounter(bkg.index);
418 for (
auto& bkg : m_backgrounds) {
420 if (bkg.fileType != BackgroundMetaData::c_ECL)
continue;
422 double mean = bkg.rate * (m_maxTimeECL - m_minTimeECL);
423 int nev = gRandom->Poisson(mean);
425 for (
int iev = 0; iev < nev; iev++) {
426 double timeShift = gRandom->Rndm() * (m_maxTimeECL - m_minTimeECL) + m_minTimeECL;
427 if (timeShift > m_minTime and timeShift < m_maxTime)
continue;
428 bkg.tree->GetEntry(bkg.eventCount);
430 if (acceptEvent(m_simHits.ECL)) {
431 double minTime = m_minTimeECL;
432 double maxTime = m_maxTimeECL;
433 if (timeShift <= m_minTime) {
438 addSimHits(eclHits, m_simHits.ECL, timeShift, minTime, maxTime);
441 std::string message =
"BeamBkgMixer: event " + to_string(bkg.eventCount)
442 +
" of " + bkg.type +
" rejected due to large energy deposit in ECL";
443 m_rejected[message] += 1;
445 if (m_rejectedCount < 10) {
446 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL");
447 }
else if (m_rejectedCount == 10) {
448 B2INFO(
"BeamBkgMixer: event rejected due to large energy deposit in ECL "
449 <<
"(message will be suppressed now)");
454 if (bkg.eventCount >= bkg.numEvents) {
456 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
457 m_reused[message] += 1;
458 if (m_reused[message] == 1) B2INFO(message);
459 bkgInfo->incrementReusedCounter(bkg.index);
466 for (
auto& bkg : m_backgrounds) {
468 if (bkg.fileType != BackgroundMetaData::c_PXD)
continue;
470 double mean = bkg.rate * (m_maxTimePXD - m_minTimePXD);
471 int nev = gRandom->Poisson(mean);
473 for (
int iev = 0; iev < nev; iev++) {
474 double timeShift = gRandom->Rndm() * (m_maxTimePXD - m_minTimePXD) + m_minTimePXD;
475 if (timeShift > m_minTime and timeShift < m_maxTime)
continue;
476 bkg.tree->GetEntry(bkg.eventCount);
478 double minTime = m_minTimePXD;
479 double maxTime = m_maxTimePXD;
480 if (timeShift <= m_minTime) {
485 addSimHits(pxdSimHits, m_simHits.PXD, timeShift, minTime, maxTime);
488 if (bkg.eventCount >= bkg.numEvents) {
490 std::string message =
"BeamBkgMixer: events of " + bkg.type +
" will be re-used";
491 m_reused[message] += 1;
492 if (m_reused[message] == 1) B2INFO(message);
493 bkgInfo->incrementReusedCounter(bkg.index);
502 void BeamBkgMixerModule::endRun()
506 void BeamBkgMixerModule::terminate()
509 B2INFO(
"BeamBkgMixer - reused samples:");
510 for (
const auto& message : m_reused) {
511 B2INFO(message.first <<
"(occured " << message.second <<
" times)");
513 B2INFO(
"BeamBkgMixer - rejected events:");
514 for (
const auto& message : m_rejected) {
515 B2INFO(message.first <<
"(occured " << message.second <<
" times)");
518 for (
auto& bkg : m_backgrounds) {
525 bool BeamBkgMixerModule::isComponentIncluded(std::vector<std::string>& components,
526 const std::string& component)
528 if (m_components.empty())
return true;
529 auto iterator = std::find(components.begin(), components.end(), component);
530 if (iterator != components.end()) {
531 components.erase(iterator);
539 const std::string& type,
540 const std::string& fileName,
544 for (
auto& bkg : m_backgrounds) {
545 if (tag == bkg.tag and fileType == bkg.fileType) {
546 bkg.fileNames.push_back(fileName);
547 bkg.realTime += realTime;
551 std::string ftype = type;
552 if (fileType == BackgroundMetaData::c_ECL) ftype +=
"(ECL)";
553 if (fileType == BackgroundMetaData::c_PXD) ftype +=
"(PXD)";
554 unsigned index = m_backgrounds.size();
555 m_backgrounds.push_back(
BkgFiles(tag, ftype, fileName, realTime,
556 m_overallScaleFactor, fileType, index));
560 bool BeamBkgMixerModule::acceptEvent(TClonesArray* cloneArrayECL)
562 if (!cloneArrayECL)
return true;
563 if (m_maxEdepECL == 0)
return true;
565 int numEntries = cloneArrayECL->GetEntriesFast();
566 for (
int i = 0; i < numEntries; i++) {
567 ECLHit* simHit =
static_cast<ECLHit*
>(cloneArrayECL->AddrAt(i));
568 if (simHit->
getEnergyDep() > m_maxEdepECL)
return false;