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