Belle II Software  release-05-01-25
BeamBkgGeneratorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Staric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // Own include
12 #include <background/modules/BeamBkgGenerator/BeamBkgGeneratorModule.h>
13 
14 // framework - DataStore
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/datastore/StoreObjPtr.h>
17 
18 // framework aux
19 #include <framework/gearbox/Unit.h>
20 #include <framework/gearbox/Const.h>
21 #include <framework/logging/Logger.h>
22 #include <framework/gearbox/GearDir.h>
23 
24 // data objects
25 #include <framework/dataobjects/EventMetaData.h>
26 #include <mdst/dataobjects/MCParticle.h>
27 #include <generators/SAD/dataobjects/SADMetaHit.h>
28 
29 // random generator
30 #include <TRandom.h>
31 
32 
33 using namespace std;
34 
35 namespace Belle2 {
41  //-----------------------------------------------------------------
42  // Register module
43  //-----------------------------------------------------------------
44 
45  REG_MODULE(BeamBkgGenerator)
46 
47  //-----------------------------------------------------------------
48  // Implementation
49  //-----------------------------------------------------------------
50 
52 
53  {
54  // set module description
55  setDescription("Beam background generator based on SAD files. "
56  "The generator picks up particles from the SAD file randomly "
57  "according to their rates. "
58  "Number of events is determined from 'realTime' and overall rate, and "
59  "the generator terminates the execution when this number is reached.");
60 
61  // Add parameters
62  addParam("fileName", m_fileName, "name of the SAD file converted to root");
63  addParam("treeName", m_treeName, "name of the TTree in the SAD file", string("sad"));
64  addParam("ringName", m_ringName, "name of the superKEKB ring (LER or HER)");
65  addParam("realTime", m_realTime,
66  "equivalent superKEKB running time to generate sample [ns].");
67 
68  }
69 
70  BeamBkgGeneratorModule::~BeamBkgGeneratorModule()
71  {
72  }
73 
74  void BeamBkgGeneratorModule::initialize()
75  {
76  // register MCParticles
77 
78  StoreArray<MCParticle> mcParticles;
79  mcParticles.registerInDataStore();
80 
81  StoreArray<SADMetaHit> sadHits;
82  sadHits.registerInDataStore();
83 
84 
85  // check steering parameters
86 
87  if (m_ringName != "LER" and m_ringName != "HER") {
88  B2ERROR("ringName can only be LER or HER");
89  }
90  if (m_realTime <= 0) {
91  B2ERROR("realTime must be positive");
92  }
93 
94  // open SAD file and get TTree
95 
96  m_file = TFile::Open(m_fileName.c_str());
97  if (!m_file) {
98  B2ERROR(m_fileName << ": can't open file");
99  return;
100  }
101  m_tree = (TTree*) m_file->Get(m_treeName.c_str());
102  if (!m_tree) {
103  B2ERROR(m_treeName << ": no such TTree in the SAD file");
104  return;
105  }
106 
107  m_tree->SetBranchAddress("s", &m_sad.s);
108  m_tree->SetBranchAddress("x", &m_sad.x);
109  m_tree->SetBranchAddress("px", &m_sad.px);
110  m_tree->SetBranchAddress("y", &m_sad.y);
111  m_tree->SetBranchAddress("py", &m_sad.py);
112  m_tree->SetBranchAddress("E", &m_sad.E);
113  m_tree->SetBranchAddress("rate", &m_sad.rate);
114 
115  // for SADMetaHit
116  m_tree->SetBranchAddress("ss", &m_sad.ss);
117  m_tree->SetBranchAddress("sraw", &m_sad.sraw);
118  m_tree->SetBranchAddress("ssraw", &m_sad.ssraw);
119  m_tree->SetBranchAddress("nturn", &m_sad.nturn);
120  m_tree->SetBranchAddress("xraw", &m_sad.xraw);
121  m_tree->SetBranchAddress("yraw", &m_sad.yraw);
122  m_tree->SetBranchAddress("r", &m_sad.r);
123  m_tree->SetBranchAddress("rr", &m_sad.rr);
124  m_tree->SetBranchAddress("dp_over_p0", &m_sad.dp_over_p0);
125  m_tree->SetBranchAddress("watt", &m_sad.watt);
126 
127  int numEntries = m_tree->GetEntries();
128  if (numEntries <= 0) {
129  B2ERROR("SAD tree is empty");
130  return;
131  }
132 
133  // construct vector of cumulative rates and determine number of events to generate
134 
135  m_rates.push_back(0);
136  for (int i = 0; i < numEntries; i++) {
137  m_tree->GetEntry(i);
138  m_rates.push_back(m_sad.rate + m_rates.back());
139  m_counters.push_back(0);
140  }
141 
142  m_numEvents = gRandom->Poisson(m_rates.back() * m_realTime / Unit::s);
143 
144  // set rotation from SAD to Belle II
145 
146  if (m_ringName == "LER") {
147  GearDir ring("/Detector/SuperKEKB/LER/");
148  m_rotation.RotateY(ring.getAngle("angle"));
149  m_ring = 2;
150  } else {
151  GearDir ring("/Detector/SuperKEKB/HER/");
152  m_rotation.RotateY(ring.getAngle("angle"));
153  m_ring = 1;
154  }
155 
156  m_sectionOrdering.insert(m_sectionOrdering.end(), {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2});
157 
158  B2RESULT("BG rate: " << m_rates.back() / 1e6 << " MHz, events to generate: "
159  << m_numEvents);
160 
161  }
162 
163  void BeamBkgGeneratorModule::beginRun()
164  {
165  }
166 
167  void BeamBkgGeneratorModule::event()
168  {
169 
170  // check event counter
171 
172  if (m_eventCounter >= m_numEvents) {
173  StoreObjPtr<EventMetaData> eventMetaData;
174  eventMetaData->setEndOfData(); // stop event processing
175  return;
176  }
177  m_eventCounter++;
178 
179  // get SAD entry
180 
181  int i = generateEntry();
182  m_tree->GetEntry(i);
183  m_counters[i]++;
184 
185  // make SADMetaHit
186  int ring_section = -1;
187  double ssraw = m_sad.ss;
188  if (m_sad.ss < 0) ssraw += 3016.;
189  int section = (int)(ssraw * 12 / 3016.);
190  if ((unsigned)section < m_sectionOrdering.size()) ring_section = m_sectionOrdering[section];
191 
192  StoreArray<SADMetaHit> SADMetaHits;
193  SADMetaHits.appendNew(SADMetaHit(m_sad.ssraw, m_sad.sraw, m_sad.ss, m_sad.s,
194  0., m_sad.nturn,
195  m_sad.x, m_sad.y, m_sad.px, m_sad.py, m_sad.xraw, m_sad.yraw,
196  m_sad.r, m_sad.rr, m_sad.dp_over_p0, m_sad.E, m_sad.rate,
197  m_sad.watt, m_ring, ring_section));
198 
199 
200  // transform to Belle II (flip sign of y and s, rotate)
201 
202  TVector3 position(m_sad.x * Unit::m, -m_sad.y * Unit::m, -m_sad.s * Unit::m);
203  position.Transform(m_rotation);
204 
205  double pz = sqrt(m_sad.E * m_sad.E - m_sad.px * m_sad.px - m_sad.py * m_sad.py);
206  if (m_ringName == "LER") pz = -pz;
207  TVector3 momentum(m_sad.px, -m_sad.py, pz);
208  momentum.Transform(m_rotation);
209 
210  // PDG code and mass
211 
212  int pdgCode = Const::electron.getPDGCode();
213  double mass = Const::electron.getMass();
214  if (m_ringName == "LER") pdgCode = -pdgCode;
215 
216  // append SAD particle to MCParticles
217 
218  StoreArray<MCParticle> MCParticles;
219  MCParticle* part = MCParticles.appendNew();
220  part->setPDG(pdgCode);
221  part->setMass(mass);
222  part->setProductionVertex(position);
223  part->setProductionTime(0);
224  part->setMomentum(momentum);
225  part->setEnergy(sqrt(momentum.Mag2() + mass * mass));
226  part->setValidVertex(true);
227  part->setStatus(MCParticle::c_PrimaryParticle);
228  part->addStatus(MCParticle::c_StableInGenerator);
229 
230  }
231 
232 
233  void BeamBkgGeneratorModule::endRun()
234  {
235  }
236 
237  void BeamBkgGeneratorModule::terminate()
238  {
239  // close SAD file
240  m_file->Close();
241 
242  B2RESULT("BG rate: " << m_rates.back() / 1e6 << " MHz, equivalent time: "
243  << m_realTime / Unit::us << " us, "
244  << m_eventCounter << " events generated");
245  if (m_eventCounter != m_numEvents)
246  B2ERROR("Number of generated events does not match the equivalent running time: "
247  << m_numEvents << " events needed, but "
248  << m_eventCounter << " generated");
249 
250  int imax = 0;
251  double average = 0;
252  for (unsigned i = 0; i < m_counters.size(); i++) {
253  if (m_counters[i] > m_counters[imax]) imax = i;
254  average += m_counters[i];
255  }
256  average /= m_counters.size();
257  B2RESULT("SAD particle usage: on average " << average << " times, max "
258  << m_counters[imax] << " times (entry " << imax << ")");
259 
260  }
261 
262 
263  int BeamBkgGeneratorModule::generateEntry() const
264  {
265  double rate = gRandom->Uniform(m_rates.back());
266  int i1 = 0;
267  int i2 = m_rates.size() - 1;
268  while (i2 - i1 > 1) {
269  int i = (i1 + i2) / 2;
270  if (rate > m_rates[i]) {
271  i1 = i;
272  } else {
273  i2 = i;
274  }
275  }
276  return i1;
277  }
278 
279 
281 } // end Belle2 namespace
282 
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::gearbox::Interface::getAngle
double getAngle(const std::string &path="") const noexcept(false)
Get the parameter path as a double converted to the standard angle unit.
Definition: Interface.h:301
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::SADMetaHit
ClassSADMetaHit - digitization simulated metahit for the SAD detector.
Definition: SADMetaHit.h:35
Belle2::BeamBkgGeneratorModule
Beam background generator based on SAD files.
Definition: BeamBkgGeneratorModule.h:37
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::GearDir
GearDir is the basic class used for accessing the parameter store.
Definition: GearDir.h:41
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::StoreArray< MCParticle >