Belle II Software development
BeamParametersFitter.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 header. */
10#include <mdst/calibration/BeamParametersFitter.h>
11
12/* Basf2 headers. */
13#include <framework/database/Database.h>
14#include <framework/database/DBImportObjPtr.h>
15#include <framework/database/DBStore.h>
16#include <framework/gearbox/Const.h>
17#include <framework/logging/Logger.h>
18
19/* ROOT headers. */
20#include <TMinuit.h>
21#include <TVectorD.h>
22#include <Math/Vector3D.h>
23#include <Math/Vector4D.h>
24#include <Math/VectorUtil.h>
25#include <Math/RotationY.h>
26
27using namespace Belle2;
28
30static double s_InvariantMass;
31
33static double s_InvariantMassError;
34
36static TVector3 s_BoostVector;
37
39static TMatrixDSym s_BoostVectorInverseCovariance(3);
40
42static TVector3 s_DirectionHER;
43
45static TVector3 s_DirectionLER;
46
48static double s_AngleError;
49
53static ROOT::Math::PxPyPzEVector getMomentum(double energy, double thetaX, double thetaY,
54 bool ler)
55{
56 const double pz = std::sqrt(energy * energy -
58 const double sx = sin(thetaX);
59 const double cx = cos(thetaX);
60 const double sy = sin(thetaY);
61 const double cy = cos(thetaY);
62 const double px = sy * cx * pz;
63 const double py = -sx * pz;
64 ROOT::Math::PxPyPzEVector result(px, py, cx * cy * pz, energy);
65 if (ler) {
66 ROOT::Math::RotationY rotationY(M_PI);
67 result = rotationY(result);
68 }
69 return result;
70}
71
72static void fcn(int& npar, double* grad, double& fval, double* par, int iflag)
73{
74 (void)npar;
75 (void)grad;
76 (void)iflag;
77 ROOT::Math::PxPyPzEVector pHER, pLER;
78 pHER = getMomentum(par[0], par[1], par[2], false);
79 pLER = getMomentum(par[3], par[4], par[5], true);
80 ROOT::Math::PxPyPzEVector pBeam = pHER + pLER;
81 ROOT::Math::XYZVector beamBoost = pBeam.BoostToCM();
82 TVectorD boostDifference(3);
83 boostDifference[0] = beamBoost.X() - s_BoostVector.X();
84 boostDifference[1] = beamBoost.Y() - s_BoostVector.Y();
85 boostDifference[2] = beamBoost.Z() - s_BoostVector.Z();
86 double boostChi2 = s_BoostVectorInverseCovariance.Similarity(boostDifference);
87 double invariantMass = pBeam.M();
88 double massChi2 = pow((invariantMass - s_InvariantMass) /
89 s_InvariantMassError, 2);
90 double angleHER = ROOT::Math::VectorUtil::Angle(pHER.Vect(), s_DirectionHER);
91 double angleLER = ROOT::Math::VectorUtil::Angle(pLER.Vect(), s_DirectionLER);
92 double angleChi2 = pow(angleHER / s_AngleError, 2) +
93 pow(angleLER / s_AngleError, 2);
94 fval = boostChi2 + massChi2 + angleChi2;
95}
96
98{
99 /* DataStore. */
101 StoreObjPtr<EventMetaData> eventMetaData;
102 eventMetaData.registerInDataStore();
104 /* Database. */
105 if (eventMetaData.isValid()) {
106 eventMetaData->setExperiment(m_IntervalOfValidity.getExperimentLow());
107 eventMetaData->setRun(m_IntervalOfValidity.getRunLow());
108 } else {
109 eventMetaData.construct(1, m_IntervalOfValidity.getRunLow(),
111 }
112 DBStore& dbStore = DBStore::Instance();
113 dbStore.update();
114 dbStore.updateEvent();
115}
116
117/*
118
119 In BeamParameters, the momenta are represented as (E, theta_x, theta_y).
120The cartesian coordinates of the momenta are given by
121
122 |p_x| |cos(theta_y) 0 sin(theta_y)| | 1 0 0 |
123 |p_y| = |0 1 0 | * | 0 cos(theta_x) -sin(theta_x) |
124 |p_x| |-sin(theta_y) 0 cos(theta_y)| | 0 sin(theta_x) cos(theta_x) |
125
126 | 0 |
127 * | 0 | ,
128 | p |
129
130or, after, matrix multiplication,
131
132 p_x = p * cos(theta_x) * sin(theta_y) ,
133 p_y = - p * sin(theta_x) ,
134 p_z = p * cos(theta_x) * cos(theta_y) .
135
136 The block form of the full covariance matrix is
137
138 | V_HER | 0 |
139V = |-------|-------| ,
140 | 0 | V_LER |
141
142where V_HER and V_LER are the covariance matrices of high-energy and
143low-energy beam momenta. The directions are assumed to be known exactly,
144thus, the covariance matrix has only two non-zero elements:
145V_11 = sigma_{E_HER}^2 and V_44 = sigma_{E_LER}^2.
146
147 It is necessary to reproduce the measured variance (note that error represents
148energy spread rather than the uncertainty of the mean energy) of collision
149invariant mass. The corresponding 1x1 error matrix is given by
150
151 | \partial \sqrt{s} | | \partial \sqrt{s} | ^ T
152V_{\sqrt{s}} = | ----------------- | * V * | ----------------- | .
153 | \partial p_i | | \partial p_i |
154
155Since there are only two non-zero elements, this formula reduces to
156
157 | \partial \sqrt{s} | ^ 2
158\sigma^2_{\sqrt{s}} = | ----------------- | * sigma_{E_HER}^2
159 | \partial E_HER |
160
161 | \partial \sqrt{s} | ^ 2
162 + | ----------------- | * sigma_{E_LER}^2 .
163 | \partial E_LER |
164
165 The derivatives are given by
166
167\partial \sqrt{s}
168----------------- =
169 \partial E_HER
170
171 1 E_HER
172= -------- [ E_beam - (p_beam)_x * cos(theta_x) * sin(theta_y) * -----
173 \sqrt{s} p_HER
174
175 E_HER
176 + (p_beam)_y * sin(theta_x) * -----
177 p_HER
178
179 E_HER
180 - (p_beam)_z * cos(theta_x) * cos(theta_y) * ----- ],
181 p_HER
182
183and by similar formula for E_LER.
184
185 Now it is necessary to make some assumption about the relation between
186sigma_{E_HER} and sigma_{E_LER}. It is assumed that sigma_{E_HER} = k E_HER
187and sigma_{E_LER} = k E_LER.
188
189*/
191{
192 int minuitResult;
194 /* Get p_HER and p_LER from a fit. */
195 double herMomentum, herThetaX, herThetaY;
196 double lerMomentum, lerThetaX, lerThetaY;
197 s_BoostVector = m_CollisionBoostVector->getBoost();
198 s_BoostVectorInverseCovariance = m_CollisionBoostVector->getBoostCovariance();
199 if (s_BoostVectorInverseCovariance.Determinant() == 0) {
200 B2WARNING("Determinant of boost covariance matrix is 0, "
201 "using generic inverse covariance matrix for fit.");
202 s_BoostVectorInverseCovariance[0][0] = 1.0 / (m_BoostError * m_BoostError);
203 s_BoostVectorInverseCovariance[0][1] = 0;
204 s_BoostVectorInverseCovariance[0][2] = 0;
205 s_BoostVectorInverseCovariance[1][0] = 0;
206 s_BoostVectorInverseCovariance[1][1] = 1.0 / (m_BoostError * m_BoostError);
207 s_BoostVectorInverseCovariance[1][2] = 0;
208 s_BoostVectorInverseCovariance[2][0] = 0;
209 s_BoostVectorInverseCovariance[2][1] = 0;
210 s_BoostVectorInverseCovariance[2][2] = 1.0 / (m_BoostError * m_BoostError);
211 } else {
212 s_BoostVectorInverseCovariance.Invert();
213 }
214 s_InvariantMass = m_CollisionInvariantMass->getMass();
215 s_InvariantMassError = m_CollisionInvariantMass->getMassError();
216 if (s_InvariantMassError == 0) {
217 B2WARNING("Invariant-mass errror is 0, using generic error for fit.");
218 s_InvariantMassError = m_InvariantMassError;
219 }
220 s_DirectionHER.SetX(0);
221 s_DirectionHER.SetY(0);
222 s_DirectionHER.SetZ(1);
223 s_DirectionHER.RotateY(m_AngleHER);
224 s_DirectionLER.SetX(0);
225 s_DirectionLER.SetY(0);
226 s_DirectionLER.SetZ(1);
227 s_DirectionLER.RotateY(m_AngleLER + M_PI);
228 s_AngleError = m_AngleError;
229 TMinuit minuit(6);
230 if (!m_Verbose)
231 minuit.SetPrintLevel(-1);
232 minuit.SetFCN(fcn);
233 minuit.mnparm(0, "PHER_E", 7, 0.01, 0, 0, minuitResult);
234 minuit.mnparm(1, "PHER_TX", 0, 0.01, 0, 0, minuitResult);
235 minuit.mnparm(2, "PHER_TY", 0, 0.01, 0, 0, minuitResult);
236 minuit.mnparm(3, "PLER_E", 4, 0.01, 0, 0, minuitResult);
237 minuit.mnparm(4, "PLER_TX", 0, 0.01, 0, 0, minuitResult);
238 minuit.mnparm(5, "PLER_TY", 0, 0.01, 0, 0, minuitResult);
239 minuit.mncomd("FIX 2 3 5 6", minuitResult);
240 minuit.mncomd("MIGRAD 10000", minuitResult);
241 minuit.mncomd("RELEASE 2 3 5 6", minuitResult);
242 minuit.mncomd("MIGRAD 10000", minuitResult);
243 double error;
244 minuit.GetParameter(0, herMomentum, error);
245 minuit.GetParameter(1, herThetaX, error);
246 minuit.GetParameter(2, herThetaY, error);
247 minuit.GetParameter(3, lerMomentum, error);
248 minuit.GetParameter(4, lerThetaX, error);
249 minuit.GetParameter(5, lerThetaY, error);
250 /* Calculate error. */
251 ROOT::Math::PxPyPzEVector pHER = getMomentum(herMomentum, herThetaX, herThetaY, false);
252 ROOT::Math::PxPyPzEVector pLER = getMomentum(lerMomentum, lerThetaX, lerThetaY, true);
253 ROOT::Math::PxPyPzEVector pBeam = pHER + pLER;
254 double fittedInvariantMass = pBeam.M();
255 B2RESULT("Initial invariant mass: " << s_InvariantMass <<
256 "; fitted invariant mass: " << fittedInvariantMass);
257 double cosThetaX = cos(herThetaX);
258 double sinThetaX = sin(herThetaX);
259 double cosThetaY = cos(herThetaY);
260 double sinThetaY = sin(herThetaY);
261 double herPartial =
262 (pBeam.E() - pHER.E() / pHER.P() *
263 (pBeam.Px() * cosThetaX * sinThetaY - pBeam.Py() * sinThetaX +
264 pBeam.Pz() * cosThetaX * cosThetaY)) / fittedInvariantMass;
265 cosThetaX = cos(lerThetaX);
266 sinThetaX = sin(lerThetaX);
267 cosThetaY = cos(lerThetaY + M_PI);
268 sinThetaY = sin(lerThetaY + M_PI);
269 double lerPartial =
270 (pBeam.E() - pLER.E() / pLER.P() *
271 (pBeam.Px() * cosThetaX * sinThetaY - pBeam.Py() * sinThetaX +
272 pBeam.Pz() * cosThetaX * cosThetaY)) / fittedInvariantMass;
273 double sigmaInvariantMass = m_CollisionInvariantMass->getMassSpread();
274 double k = sqrt(sigmaInvariantMass * sigmaInvariantMass /
275 (pow(herPartial * pHER.E(), 2) +
276 pow(lerPartial * pLER.E(), 2)));
277 double herSpread = k * pHER.E();
278 double lerSpread = k * pLER.E();
279 B2INFO("Invariant mass spread: " << sigmaInvariantMass);
280 B2RESULT("HER energy spread: " << herSpread <<
281 "; LER energy spread: " << lerSpread);
282 /* Fill beam parameters. */
285 TMatrixDSym covariance(3);
286 for (int i = 0; i < 3; ++i) {
287 for (int j = 0; j < 3; ++j)
288 covariance[i][j] = 0;
289 }
290 covariance[0][0] = herSpread * herSpread;
291 m_BeamParameters.setCovHER(covariance);
292 covariance[0][0] = lerSpread * lerSpread;
293 m_BeamParameters.setCovLER(covariance);
294}
295
297 double covarianceXX, double covarianceYY)
298{
300 m_BeamParameters.setVertex(ROOT::Math::XYZVector(m_BeamSpot->getIPPosition()));
301 TMatrixDSym beamSize = m_BeamSpot->getSizeCovMatrix();
302 double xScale, yScale;
303 if (covarianceXX < 0)
304 xScale = 1;
305 else
306 xScale = sqrt(covarianceXX / beamSize[0][0]);
307 if (covarianceYY < 0)
308 yScale = 1;
309 else
310 yScale = sqrt(covarianceYY / beamSize[1][1]);
311 for (int i = 0; i < 3; ++i) {
312 beamSize[0][i] *= xScale;
313 beamSize[i][0] *= xScale;
314 beamSize[1][i] *= yScale;
315 beamSize[i][1] *= yScale;
316 }
318}
319
321{
322 DBImportObjPtr<BeamParameters> beamParametersImport;
323 beamParametersImport.construct(m_BeamParameters);
324 beamParametersImport.import(m_IntervalOfValidity);
325}
BeamParameters m_BeamParameters
Beam parameters.
void setupDatabase()
Setup database.
double m_InvariantMassError
Invariant-mass error (use only if error is 0).
DBObjPtr< CollisionBoostVector > m_CollisionBoostVector
Collision boost vector.
bool m_Verbose
Whether to be verbose (print Minuit output).
void importBeamParameters()
Import beam parameters.
DBObjPtr< CollisionInvariantMass > m_CollisionInvariantMass
Collision invariant mass.
IntervalOfValidity m_IntervalOfValidity
Interval of validity.
DBObjPtr< BeamSpot > m_BeamSpot
Beam spot.
double m_BoostError
Boost error (use only if inverse error matrix is not available).
void fillVertexData(double covarianceXX, double covarianceYY)
Fill beam spot (vertex) data.
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 setCovLER(const TMatrixDSym &cov)
Set the covariance matrix for LER (E, theta_x, theta_y) where E is the energy, theta_x is the horizon...
void setCovHER(const TMatrixDSym &cov)
Set the covariance matrix for HER (E, theta_x, theta_y) where E is the energy, theta_x is the horizon...
void setVertex(const ROOT::Math::XYZVector &vertex, const std::vector< double > &cov)
Set the vertex position and error matrix.
void setCovVertex(const TMatrixDSym &cov)
Set the covariance matrix of the vertex position.
static const double electronMass
electron mass
Definition: Const.h:685
bool import(const IntervalOfValidity &iov)
Import the object to database.
Definition: DBImportBase.cc:36
Class for importing a single object to the database.
void construct(Args &&... params)
Construct an object of type T in this DBImportObjPtr using the provided constructor arguments.
Singleton class to cache database objects.
Definition: DBStore.h:31
static DataStore & Instance()
Instance of singleton Store.
Definition: DataStore.cc:54
void setInitializeActive(bool active)
Setter for m_initializeActive.
Definition: DataStore.cc:94
int getRunLow() const
Getter for lowest run number.
int getExperimentLow() const
Getter for lowest experiment number.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:96
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:111
bool construct(Args &&... params)
Construct an object of type T in this StoreObjPtr, using the provided constructor arguments.
Definition: StoreObjPtr.h:119
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:28
void updateEvent()
Updates all intra-run dependent objects.
Definition: DBStore.cc:142
void update()
Updates all objects that are outside their interval of validity.
Definition: DBStore.cc:79
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
Abstract base class for different kinds of events.