Belle II Software light-2406-ragdoll
B2BIIConvertBeamParamsModule.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#include <b2bii/modules/B2BIIMdstInput/B2BIIConvertBeamParamsModule.h>
10#include <belle_legacy/benergy/BeamEnergy.h>
11#include <belle_legacy/ip/IpProfile.h>
12#include <belle_legacy/panther/panther.h>
13#include <belle_legacy/tables/belletdf.h>
14
15#include <framework/database/Database.h>
16#include <framework/database/EventDependency.h>
17#include <framework/database/IntervalOfValidity.h>
18#include <framework/dbobjects/BeamParameters.h>
19#include <framework/geometry/B2Vector3.h>
20#include <framework/utilities/FileSystem.h>
21#include <mdst/dbobjects/BeamSpot.h>
22#include <mdst/dbobjects/CollisionAxisCMS.h>
23#include <mdst/dbobjects/CollisionBoostVector.h>
24#include <mdst/dbobjects/CollisionInvariantMass.h>
25
26#include <Math/VectorUtil.h>
27
28#include <cmath>
29#include <list>
30
31namespace Belle {
35 int check_beginrun(bool) { return 0; }
36}
37
38namespace {
40 TVector3 CLHEPtoROOT(const HepPoint3D& point)
41 {
42 return TVector3(point.x(), point.y(), point.z());
43 }
45 TMatrixDSym CLHEPtoROOT(const HepSymMatrix& matrix)
46 {
47 TMatrixDSym result(matrix.num_col());
48 for (int i = 0; i < matrix.num_row(); ++i) {
49 for (int j = 0; j < matrix.num_col(); ++j) {
50 result(i, j) = matrix(i + 1, j + 1);
51 }
52 }
53 return result;
54 }
55}
56
57namespace Belle2 {
63 REG_MODULE(B2BIIConvertBeamParams);
64
66 {
67 addParam("mcFlag", m_mcFlag, "which mc flag to use", m_mcFlag);
68 addParam("missingBenergy", m_missingBenergy, "where to store information about runs with missing database info", m_missingBenergy);
69 addParam("missingIp", m_missingIp, "where to store information about runs with missing IP profile info", m_missingIp);
70 addParam("smearEnergy", m_SmearEnergy, "Allow beam energy smearing for MC generation", true);
71 addParam("smearDirection", m_SmearDirection, "Allow beam direction smearing for MC generation", true);
72 addParam("smearVertex", m_SmearVertex, "Allow IP position smearing for MC generation", true);
73 addParam("generateCMS", m_GenerateCMS, "Generate events in CMS, not lab system.", false);
74 addParam("storeBeamParameters", m_storeBeamParameters, "Store the BeamParameters payloads in the localDB", true);
75 addParam("storeCollisionAxisCMS", m_storeCollisionAxisCMS,
76 "Store the CollisionAxisCMS payloads in the localDB",
77 true);
78 addParam("storeCollisionInvariantMass", m_storeCollisionInvariantMass, "Store the CollisionInvariantMass payloads in the localDB",
79 true);
80 addParam("storeCollisionBoostVector", m_storeCollisionBoostVector, "Store the CollisionBoostVector payloads in the localDB", true);
81 addParam("storeBeamSpot", m_storeBeamSpot, "Store the BeamSpot payloads in the localDB", true);
82 }
83
85 {
86 BsInit(0);
87 m_event.isRequired();
88 }
89
91 {
92 // Fudge a runhead record with all the info we need to be able to get the
93 // Beamenergy/IpProfile from the database
94 BsClrTab(BBS_CLEAR);
95 Belle::Belle_runhead_Manager& runMgr = Belle::Belle_runhead_Manager::get_manager();
96 runMgr.add();
97 Belle::Belle_runhead& runhead = runMgr[0];
98 runhead.ExpMC(m_mcFlag);
99 runhead.ExpNo(m_event->getExperiment());
100 runhead.RunNo(m_event->getRun());
101 B2INFO("Obtaining values for exp " << m_event->getExperiment() << ", run " << m_event->getRun());
102 // and then get the values from the database
103 Belle::BeamEnergy::flag_database(0);
104 Belle::IpProfile::begin_run();
105 Belle::Ip_profile_general_Manager& ipMgr = Belle::Ip_profile_general_Manager::get_manager();
106
107 if (!Belle::BeamEnergy::usable()) {
109 if (!lock.lock()) {
110 B2ERROR("No BeamEnergy for exp " << m_event->getExperiment() << ", run " << m_event->getRun());
111 } else {
112 std::ofstream file(m_missingBenergy.c_str(), std::ios::app);
113 file << m_event->getExperiment() << "," << m_event->getRun() << std::endl;
114 }
115 return;
116 }
117
118 const double Eher = Belle::BeamEnergy::E_HER();
119 const double Eler = Belle::BeamEnergy::E_LER();
120 const double crossingAngle = Belle::BeamEnergy::Cross_angle();
121 const double angleLer = M_PI; //parallel to negative z axis (different from Belle II!)
122 const double angleHer = crossingAngle; //in positive z and x direction, verified to be consistent with Upsilon(4S) momentum
123
124 // Beam energy spread taken from the hard-coded values in the Belle evtgen:
125 // /sw/belle/belle/b20090127_0910/src/util/evtgenutil/basf_if/evtgen.cc
126 // The beam energy spread seems not have been logged in the belle DB at all.
127 // If you are reading this and you know how to access the run-dependent beam energy spread,
128 // please fix this (April 2021)!
129 std::vector<double> covarianceHer = {0.00513 * 0.00513}; // Energy spread only. No idea about the direction spread
130 std::vector<double> covarianceLer = {0.002375 * 0.002375}; // Energy spread only. No idea about the direction spread
131
132 IntervalOfValidity iov(m_event->getExperiment(), m_event->getRun(), m_event->getExperiment(), m_event->getRun());
133
134 BeamParameters beamParams;
135 beamParams.setLER(Eler, angleLer, 0, covarianceLer);
136 beamParams.setHER(Eher, angleHer, 0, covarianceHer);
137
138 // set the flags according to the module settings.
139 int flags = 0;
140 if (m_GenerateCMS)
142 if (m_SmearEnergy)
146 if (m_SmearVertex)
148 beamParams.setGenerationFlags(flags);
149
150 CollisionAxisCMS collisionAxisCMS;
151 CollisionBoostVector collisionBoostVector;
152 CollisionInvariantMass collisionInvM;
153 ROOT::Math::PxPyPzEVector momentumHER = beamParams.getHER();
154 ROOT::Math::PxPyPzEVector cms = momentumHER + beamParams.getLER();
155 ROOT::Math::XYZVector boost = -cms.BoostToCM();
156 ROOT::Math::VectorUtil::boost(momentumHER, boost);
157 double angleXZ = std::atan(momentumHER.X() / momentumHER.Z());
158 double angleYZ = std::atan(momentumHER.Y() / momentumHER.Z());
159 collisionAxisCMS.setAngles(angleXZ, angleYZ, TMatrixTSym<double>(2));
160 collisionAxisCMS.setSpread(TMatrixTSym<double>(2), 0, 0, 0);
161 collisionBoostVector.setBoost(B2Vector3D(boost), TMatrixTSym<double>(3));
162 //note: maybe we could use Belle::BeamEnergy::E_beam_corr(), Belle::BeamEnergy::E_beam_err()
163 collisionInvM.setMass(cms.M(), 0.0, 0.0);
164
165 // Boost vector and invariant mass are not intra-run dependent, store now
168 "CollisionAxisCMS", &collisionAxisCMS, iov);
169 }
172 "CollisionBoostVector", &collisionBoostVector, iov);
173 }
176 "CollisionInvariantMass", &collisionInvM, iov);
177 }
178
179 // and now we continue with the vertex
180 if (!Belle::IpProfile::usable()) {
182 if (!lock.lock()) {
183 B2ERROR("No IpProfile for exp " << m_event->getExperiment() << ", run " << m_event->getRun());
184 } else {
185 std::ofstream file(m_missingIp.c_str(), std::ios::app);
186 file << m_event->getExperiment() << "," << m_event->getRun() << std::endl;
187 }
188 std::vector<double> covariance; // 0 = no error
189 beamParams.setVertex(ROOT::Math::XYZVector(
190 std::numeric_limits<double>::quiet_NaN(),
191 std::numeric_limits<double>::quiet_NaN(),
192 std::numeric_limits<double>::quiet_NaN()), covariance);
194 Database::Instance().storeData("BeamParameters", &beamParams, iov);
195
196 BeamSpot beamSpot;
197 beamSpot.setIP(
198 TVector3(std::numeric_limits<double>::quiet_NaN(),
199 std::numeric_limits<double>::quiet_NaN(),
200 std::numeric_limits<double>::quiet_NaN()
201 ), TMatrixTSym<double>(3)
202 );
203 if (m_storeBeamSpot) {
204 Database::Instance().storeData("BeamSpot", &beamSpot, iov);
205 B2INFO("No IpProfile, created BeamSpot Payload");
206 }
207 return;
208 }
209
210 int nbins = std::max(ipMgr.count() - 1, 1);
211 B2INFO("exp " << m_event->getExperiment() << ", run " << m_event->getRun() << ": " << nbins << " IpProfile bins");
212
213 // need to keep the objects alive until we created the payloads
214 std::list<BeamParameters> beamparamsList;
215 std::list<BeamSpot> beamspotList;
216 // and we want them to be intra-run dependent
217 std::unique_ptr<EventDependency> beamparams;
218 std::unique_ptr<EventDependency> beamspots;
219
220 Belle::Belle_event_Manager& evtMgr = Belle::Belle_event_Manager::get_manager();
221 if (!evtMgr.count()) evtMgr.add();
222 for (int i = 0; i < nbins; ++i) {
223 int evtNr = i * BIN_EVENTS;
224 // fudge belle event record to get the correct run dependent ip profile
225 Belle::Belle_event& evt = evtMgr[0];
226 evt.ExpMC(m_mcFlag);
227 evt.ExpNo(m_event->getExperiment());
228 evt.RunNo(m_event->getRun());
229 evt.EvtNo(evtNr);
230 Belle::IpProfile::set_evtbin_number();
231 B2ASSERT("something wrong: " << Belle::IpProfile::EvtBinNo() << "!=" << i, Belle::IpProfile::EvtBinNo() == i);
232
233 // and convert
234 HepPoint3D ip = Belle::IpProfile::e_position();
235 HepSymMatrix ipErr = Belle::IpProfile::e_position_err();
236
237 beamParams.setVertex(ROOT::Math::XYZVector(ip.x(), ip.y(), ip.z()));
238 beamParams.setCovVertex(CLHEPtoROOT(ipErr));
239 beamparamsList.emplace_back(beamParams);
240
241 BeamSpot beamSpot;
242 beamSpot.setIP(CLHEPtoROOT(ip), CLHEPtoROOT(ipErr));
243 beamspotList.emplace_back(beamSpot);
244
245 if (!beamparams) {
246 beamparams = std::make_unique<EventDependency>(&beamparamsList.back(), false);
247 beamspots = std::make_unique<EventDependency>(&beamspotList.back(), false);
248 } else {
249 beamparams->add(evtNr, &beamparamsList.back());
250 beamspots->add(evtNr, &beamspotList.back());
251 }
252 }
253 if (beamparamsList.size() < 1) {
254 B2ERROR("Something is wrong with exp " << m_event->getExperiment() << ", run " << m_event->getRun() << ": no bins found");
255 return;
256 }
257 if (beamparamsList.size() == 1) {
258 // just one bin? no need to store event dependency
260 Database::Instance().storeData("BeamParameters", &beamparamsList.back(), iov);
261 if (m_storeBeamSpot)
262 Database::Instance().storeData("BeamSpot", &beamspotList.back(), iov);
263 B2INFO("Created event independent payload");
264 } else {
265 // otherwise store full information
267 Database::Instance().storeData("BeamParameters", beamparams.get(), iov);
268 if (m_storeBeamSpot)
269 Database::Instance().storeData("BeamSpot", beamspots.get(), iov);
270 B2INFO("Created event dependent payload with " << beamparamsList.size() << " entries");
271 }
272 }
274} //Belle2 namespace
bool m_SmearDirection
Smear beam direction when generating initial events.
std::string m_missingBenergy
Where to store information about runs without beam energy information.
bool m_storeCollisionAxisCMS
Store the CollisionAxisCMS payloads in the local database.
std::string m_missingIp
Where to store information about runs without IP profile information.
bool m_storeCollisionInvariantMass
Store the CollisionInvariantMass payloads in the local database.
StoreObjPtr< EventMetaData > m_event
Event metadata.
bool m_SmearVertex
Smear vertex position when generating initial events.
bool m_storeBeamParameters
Store the BeamParameters payloads in the local database.
bool m_GenerateCMS
Generate events in CMS, not lab system.
bool m_storeBeamSpot
Store the BeamSpot payloads in the local database.
bool m_SmearEnergy
Smear energy when generating initial events.
bool m_storeCollisionBoostVector
Store the CollisionBoostVector payloads in the local database.
int m_mcFlag
mcFlag to use when getting belle database content
This class contains the nominal beam parameters and the parameters used for smearing of the primary v...
void setHER(double energy, double angleX, double angleY, const std::vector< double > &cov)
Set the HER FourVector and error matrix from beam energy, angle and covariance entries.
void setLER(double energy, double angleX, double angleY, const std::vector< double > &cov)
Set the LER FourVector and error matrix from beam energy, angle and covariance entries.
void setVertex(const ROOT::Math::XYZVector &vertex, const std::vector< double > &cov)
Set the vertex position and error matrix.
void setGenerationFlags(int flags) override
Set the generation flags to be used for event generation (ORed combination of EGenerationFlags).
void setCovVertex(const TMatrixDSym &cov)
Set the covariance matrix of the vertex position.
This class contains the beam spot position and size modeled as a gaussian distribution in space.
Definition: BeamSpot.h:22
void setIP(const TVector3 &ipPosition, const TMatrixDSym &covariance)
Set the IP position and its error matrix.
Definition: BeamSpot.h:59
This class contains the measured values of the orientation of the collision axis in the CM system obt...
void setSpread(const TMatrixDSym &spreadCovariance, double spreadXZunc, double spreadYZunc, double spreadPhiUnc)
Set spread covariance and uncertainties of the eigenvalues of this matrix.
void setAngles(double angleXZ, double angleYZ, const TMatrixDSym &centerCovariance)
Set the central values and uncertainty of them.
This class contains the measured average boost vector vec(beta) = (beta_x, beta_y,...
void setBoost(const TVector3 &boost, const TMatrixDSym &covariance)
Set the boost vector and its error matrix.
This class contains the measured average center-of-mass energy, which is equal to the invariant mass ...
void setMass(double mass, double error, double spread)
Set the CMS energy and its uncertainty.
Helper class for locking a file.
Definition: FileSystem.h:97
bool lock(int timeout=300, bool ignoreErrors=false)
Try to lock the file.
Definition: FileSystem.cc:189
A class that describes the interval of experiments/runs for which an object in the database is valid.
const ROOT::Math::PxPyPzEVector & getHER() const
Get 4vector of the high energy beam.
const ROOT::Math::PxPyPzEVector & getLER() const
Get 4vector of the low energy beam.
@ c_generateCMS
generate initial event in CMS instead of lab
@ c_smearBeamEnergy
smear energy of HER and LER (but not direction)
@ c_smearBeamDirection
smear direction of HER and LER (but not energy)
void initialize() override
Initialize phanter banks.
void beginRun() override
Set run info in panther and load IPProfile/Benergy and convert the values to payloads.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
static Database & Instance()
Instance of a singleton Database.
Definition: Database.cc:42
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition: Database.cc:141
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24