Belle II Software  release-06-02-00
BeamBkgMixerModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 // Own include
10 #include <background/modules/BeamBkgMixer/BeamBkgMixerModule.h>
11 
12 
13 
14 // framework - DataStore
15 #include <framework/datastore/DataStore.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/datastore/StoreObjPtr.h>
18 #include <framework/datastore/RelationArray.h>
19 
20 // framework aux
21 #include <framework/core/ModuleParam.templateDetails.h>
22 #include <framework/logging/Logger.h>
23 
24 // SimHits
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>
34 
35 // MetaData
36 #include <framework/dataobjects/EventMetaData.h>
37 #include <framework/dataobjects/BackgroundInfo.h>
38 
39 // Root
40 #include <TFile.h>
41 #include <TRandom3.h>
42 
43 //std::find
44 #include <algorithm>
45 
46 using namespace std;
47 
48 namespace Belle2 {
54  //-----------------------------------------------------------------
55  // Register module
56  //-----------------------------------------------------------------
57 
58  REG_MODULE(BeamBkgMixer)
59 
60  //-----------------------------------------------------------------
61  // Implementation
62  //-----------------------------------------------------------------
63 
65  {
66  // set module description (e.g. insert text)
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.");
71 
72  // Add parameters
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(),
85  m_scaleFactors);
86  addParam("components", m_components,
87  "Detector components to be included, empty list means all components",
88  m_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);
99 
100  addParam("beambackhits", m_BeamBackHits,
101  "If true also add the BeamBackHits collection for the selected "
102  "subdetectors to the output file", false);
103 
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);
107 
108  addParam("cacheSize", m_cacheSize,
109  "file cache size in Mbytes. If negative, use root default", 0);
110  }
111 
112  BeamBkgMixerModule::~BeamBkgMixerModule()
113  {
114  }
115 
116  void BeamBkgMixerModule::initialize()
117  {
118 
119  // included components
120 
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");
129 
130  // ignore these ones
131 
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");
144 
145  if (!components.empty()) {
146  std::string str;
147  for (unsigned i = 0; i < components.size(); ++i) str = str + " " + components[i];
148  B2WARNING("Unknown components:" << str);
149  }
150 
151  // check files and append them to sample container
152 
153  for (auto file : m_backgroundFiles) {
154 
155  // wildcarding is not allowed anymore
156  if (TString(file.c_str()).Contains("*")) {
157  B2ERROR(file << ": wildcards are not allowed");
158  continue;
159  }
160 
161  // check the file existance
162  TFile* f = TFile::Open(file.c_str(), "READ");
163  if (!f) {
164  B2ERROR(file << ": file not found");
165  continue;
166  }
167  if (!f->IsOpen()) {
168  B2ERROR(file << ": can't open file");
169  continue;
170  }
171  f->Close();
172 
173  TChain persistent("persistent");
174  int nFiles = persistent.Add(file.c_str());
175  if (nFiles == 0) {
176  B2ERROR(file << ": no such files");
177  continue;
178  }
179  if (persistent.GetEntries() == 0) {
180  B2ERROR(file << ": tree 'persistent' has no entries");
181  continue;
182  }
183 
184  TObject* bkgMetaData = 0; // Note: allocation left to root
185  TBranch* branchBMD = persistent.GetBranch("BackgroundMetaData");
186  if (!branchBMD) {
187  B2ERROR(file << ": branch 'BackgroundMetaData' not found");
188  continue;
189  }
190  branchBMD->SetAddress(&bkgMetaData);
191 
192  std::vector<BackgroundMetaData::BG_TAG> tags;
193  std::vector<std::string> types;
194  std::vector<BackgroundMetaData::EFileType> fileTypes;
195  double realTime = 0;
196  for (unsigned k = 0; k < persistent.GetEntries(); k++) {
197  persistent.GetEntry(k);
198  BackgroundMetaData* bgMD = static_cast<BackgroundMetaData*>(bkgMetaData);
199  tags.push_back(bgMD->getBackgroundTag());
200  types.push_back(bgMD->getBackgroundType());
201  fileTypes.push_back(bgMD->getFileType());
202  realTime += bgMD->getRealTime();
203  }
204  if (realTime <= 0) {
205  B2ERROR(file << ": invalid realTime: " << realTime);
206  continue;
207  }
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");
211  continue;
212  }
213  }
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");
217  continue;
218  }
219  }
220 
221  appendSample(tags[0], types[0], file, realTime, fileTypes[0]);
222 
223  }
224 
225 
226  // set scale factors
227 
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);
236  }
237  }
238 
239  // open files for reading SimHits
240 
241  for (auto& bkg : m_backgrounds) {
242 
243  // define TChain for reading SimHits
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());
247  }
248 
249  bkg.numEvents = bkg.tree->GetEntries();
250  bkg.rate = bkg.numEvents / bkg.realTime * bkg.scaleFactor;
251 
252  if (m_cacheSize >= 0) bkg.tree->SetCacheSize(m_cacheSize * 1024 * 1024);
253 
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);
270 
271  if (m_BeamBackHits and bkg.tree->GetBranch("BeamBackHits"))
272  bkg.tree->SetBranchAddress("BeamBackHits", &m_simHits.BeamBackHits);
273 
274  // print INFO
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";}
280 
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");
287  }
288 
289 
290  // SimHits registration
291 
292  StoreArray<PXDSimHit> pxdSimHits;
293  if (m_PXD) pxdSimHits.registerInDataStore();
294 
295  StoreArray<SVDSimHit> svdSimHits;
296  if (m_SVD) svdSimHits.registerInDataStore();
297 
298  StoreArray<CDCSimHit> cdcSimHits;
299  if (m_CDC) cdcSimHits.registerInDataStore();
300 
301  StoreArray<TOPSimHit> topSimHits;
302  if (m_TOP) topSimHits.registerInDataStore();
303 
304  StoreArray<ARICHSimHit> arichSimHits;
305  if (m_ARICH) arichSimHits.registerInDataStore();
306 
307  StoreArray<ECLHit> eclHits;
308  if (m_ECL) eclHits.registerInDataStore();
309 
310  StoreArray<BKLMSimHit> bklmSimHits;
311  if (m_KLM) bklmSimHits.registerInDataStore();
312 
313  StoreArray<EKLMSimHit> eklmSimHits;
314  if (m_KLM) eklmSimHits.registerInDataStore();
315 
316  StoreArray<BeamBackHit> beamBackHits;
317  if (m_BeamBackHits) beamBackHits.registerInDataStore();
318 
319 
320  // add BackgroundInfo to persistent tree
321 
322  StoreObjPtr<BackgroundInfo> bkgInfo("", DataStore::c_Persistent);
323  bkgInfo.registerInDataStore();
324  bkgInfo.create();
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) {
337  descr.tag = bkg.tag;
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;
345  descr.reused = 0;
346  bkgInfo->appendBackgroundDescr(descr);
347  }
348 
349  }
350 
351  void BeamBkgMixerModule::beginRun()
352  {
353  }
354 
355  void BeamBkgMixerModule::event()
356  {
357  StoreArray<PXDSimHit> pxdSimHits;
358  StoreArray<SVDSimHit> svdSimHits;
359  StoreArray<CDCSimHit> cdcSimHits;
360  StoreArray<TOPSimHit> topSimHits;
361  StoreArray<ARICHSimHit> arichSimHits;
362  StoreArray<ECLHit> eclHits;
363  StoreArray<BKLMSimHit> bklmSimHits;
364  StoreArray<EKLMSimHit> eklmSimHits;
365  StoreArray<BeamBackHit> beamBackHits;
366  StoreObjPtr<BackgroundInfo> bkgInfo("", DataStore::c_Persistent);
367 
368  for (auto& bkg : m_backgrounds) {
369 
370  if (bkg.fileType != BackgroundMetaData::c_Usual) continue;
371 
372  double mean = bkg.rate * (m_maxTime - m_minTime);
373  int nev = gRandom->Poisson(mean);
374 
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);
378 
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);
390  } else {
391  iev--;
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;
395  m_rejectedCount++;
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)");
401  }
402  }
403 
404  bkg.eventCount++;
405  if (bkg.eventCount >= bkg.numEvents) {
406  bkg.eventCount = 0;
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);
411  }
412  }
413  }
414 
415 
416  for (auto& bkg : m_backgrounds) {
417 
418  if (bkg.fileType != BackgroundMetaData::c_ECL) continue;
419 
420  double mean = bkg.rate * (m_maxTimeECL - m_minTimeECL);
421  int nev = gRandom->Poisson(mean);
422 
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);
427 
428  if (acceptEvent(m_simHits.ECL)) {
429  double minTime = m_minTimeECL;
430  double maxTime = m_maxTimeECL;
431  if (timeShift <= m_minTime) {
432  maxTime = m_minTime;
433  } else {
434  minTime = m_maxTime;
435  }
436  addSimHits(eclHits, m_simHits.ECL, timeShift, minTime, maxTime);
437  } else {
438  iev--;
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;
442  m_rejectedCount++;
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)");
448  }
449  }
450 
451  bkg.eventCount++;
452  if (bkg.eventCount >= bkg.numEvents) {
453  bkg.eventCount = 0;
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);
458  }
459  }
460 
461  }
462 
463 
464  for (auto& bkg : m_backgrounds) {
465 
466  if (bkg.fileType != BackgroundMetaData::c_PXD) continue;
467 
468  double mean = bkg.rate * (m_maxTimePXD - m_minTimePXD);
469  int nev = gRandom->Poisson(mean);
470 
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);
475 
476  double minTime = m_minTimePXD;
477  double maxTime = m_maxTimePXD;
478  if (timeShift <= m_minTime) {
479  maxTime = m_minTime;
480  } else {
481  minTime = m_maxTime;
482  }
483  addSimHits(pxdSimHits, m_simHits.PXD, timeShift, minTime, maxTime);
484 
485  bkg.eventCount++;
486  if (bkg.eventCount >= bkg.numEvents) {
487  bkg.eventCount = 0;
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);
492  }
493  }
494 
495  }
496 
497  }
498 
499 
500  void BeamBkgMixerModule::endRun()
501  {
502  }
503 
504  void BeamBkgMixerModule::terminate()
505  {
506 
507  B2INFO("BeamBkgMixer - reused samples:");
508  for (const auto& message : m_reused) {
509  B2INFO(message.first << "(occured " << message.second << " times)");
510  }
511  B2INFO("BeamBkgMixer - rejected events:");
512  for (const auto& message : m_rejected) {
513  B2INFO(message.first << "(occured " << message.second << " times)");
514  }
515 
516  for (auto& bkg : m_backgrounds) {
517  bkg.tree.reset();
518  }
519 
520  }
521 
522 
523  bool BeamBkgMixerModule::isComponentIncluded(std::vector<std::string>& components,
524  const std::string& component)
525  {
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);
530  return true;
531  }
532  return false;
533  }
534 
535 
536  void BeamBkgMixerModule::appendSample(BackgroundMetaData::BG_TAG tag,
537  const std::string& type,
538  const std::string& fileName,
539  double realTime,
541  {
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;
546  return;
547  }
548  }
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));
555  }
556 
557 
558  bool BeamBkgMixerModule::acceptEvent(TClonesArray* cloneArrayECL)
559  {
560  if (!cloneArrayECL) return true;
561  if (m_maxEdepECL == 0) return true;
562 
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;
567  }
568  return true;
569  }
570 
571 
573 } // end Belle2 namespace
574 
Metadata information about the beam background file.
EFileType
Enum for BG file types.
const std::string & getBackgroundType() const
Returns the type of background.
float getRealTime() const
Returns real time that corresponds to this background sample.
EFileType getFileType() const
Returns file type.
BG_TAG
Enum for background tags.
BG_TAG getBackgroundTag() const
Returns background tag value.
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...
Definition: ECLHit.h:25
double getEnergyDep() const
Get deposit energy.
Definition: ECLHit.h:70
Base class for Modules.
Definition: Module.h:72
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.
Definition: StoreArray.h:113
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
Abstract base class for different kinds of events.
Structure for background description.
structure to hold samples of a particular background type