12 #include <background/modules/BeamBkgGenerator/BeamBkgGeneratorModule.h>
15 #include <framework/datastore/StoreArray.h>
16 #include <framework/datastore/StoreObjPtr.h>
19 #include <framework/gearbox/Unit.h>
20 #include <framework/gearbox/Const.h>
21 #include <framework/logging/Logger.h>
22 #include <framework/gearbox/GearDir.h>
25 #include <framework/dataobjects/EventMetaData.h>
26 #include <mdst/dataobjects/MCParticle.h>
27 #include <generators/SAD/dataobjects/SADMetaHit.h>
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.");
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].");
70 BeamBkgGeneratorModule::~BeamBkgGeneratorModule()
74 void BeamBkgGeneratorModule::initialize()
79 mcParticles.registerInDataStore();
82 sadHits.registerInDataStore();
87 if (m_ringName !=
"LER" and m_ringName !=
"HER") {
88 B2ERROR(
"ringName can only be LER or HER");
90 if (m_realTime <= 0) {
91 B2ERROR(
"realTime must be positive");
96 m_file = TFile::Open(m_fileName.c_str());
98 B2ERROR(m_fileName <<
": can't open file");
101 m_tree = (TTree*) m_file->Get(m_treeName.c_str());
103 B2ERROR(m_treeName <<
": no such TTree in the SAD file");
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);
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);
127 int numEntries = m_tree->GetEntries();
128 if (numEntries <= 0) {
129 B2ERROR(
"SAD tree is empty");
135 m_rates.push_back(0);
136 for (
int i = 0; i < numEntries; i++) {
138 m_rates.push_back(m_sad.rate + m_rates.back());
139 m_counters.push_back(0);
142 m_numEvents = gRandom->Poisson(m_rates.back() * m_realTime / Unit::s);
146 if (m_ringName ==
"LER") {
147 GearDir ring(
"/Detector/SuperKEKB/LER/");
148 m_rotation.RotateY(ring.
getAngle(
"angle"));
151 GearDir ring(
"/Detector/SuperKEKB/HER/");
152 m_rotation.RotateY(ring.
getAngle(
"angle"));
156 m_sectionOrdering.insert(m_sectionOrdering.end(), {1, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2});
158 B2RESULT(
"BG rate: " << m_rates.back() / 1e6 <<
" MHz, events to generate: "
163 void BeamBkgGeneratorModule::beginRun()
167 void BeamBkgGeneratorModule::event()
172 if (m_eventCounter >= m_numEvents) {
174 eventMetaData->setEndOfData();
181 int i = generateEntry();
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];
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));
202 TVector3 position(m_sad.x * Unit::m, -m_sad.y * Unit::m, -m_sad.s * Unit::m);
203 position.Transform(m_rotation);
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);
212 int pdgCode = Const::electron.getPDGCode();
213 double mass = Const::electron.getMass();
214 if (m_ringName ==
"LER") pdgCode = -pdgCode;
220 part->setPDG(pdgCode);
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);
233 void BeamBkgGeneratorModule::endRun()
237 void BeamBkgGeneratorModule::terminate()
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");
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];
256 average /= m_counters.size();
257 B2RESULT(
"SAD particle usage: on average " << average <<
" times, max "
258 << m_counters[imax] <<
" times (entry " << imax <<
")");
263 int BeamBkgGeneratorModule::generateEntry()
const
265 double rate = gRandom->Uniform(m_rates.back());
267 int i2 = m_rates.size() - 1;
268 while (i2 - i1 > 1) {
269 int i = (i1 + i2) / 2;
270 if (rate > m_rates[i]) {