Belle II Software  release-05-01-25
MassConstraint.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * See https://github.com/tferber/OrcaKinfit, forked from *
4  * https://github.com/iLCSoft/MarlinKinfit *
5  * *
6  * Further information about the fit engine and the user interface *
7  * provided in MarlinKinfit can be found at *
8  * https://www.desy.de/~blist/kinfit/doc/html/ *
9  * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available *
10  * from http://www-flc.desy.de/lcnotes/ *
11  * *
12  * Adopted by: Torben Ferber (torben.ferber@desy.de) (TF) *
13  * *
14  * This software is provided "as is" without any warranty. *
15  **************************************************************************/
16 
17 #include "analysis/OrcaKinFit/MassConstraint.h"
18 #include "analysis/OrcaKinFit/ParticleFitObject.h"
19 #include <framework/logging/Logger.h>
20 
21 #include<iostream>
22 #include<cmath>
23 
24 #undef NDEBUG
25 #include<cassert>
26 
27 namespace Belle2 {
32  namespace OrcaKinFit {
33 
34 
35 // constructor
36  MassConstraint::MassConstraint(double mass_)
37  : mass(mass_)
38  {}
39 
40 // destructor
41  MassConstraint::~MassConstraint() = default;
42 
43 // calulate current value of constraint function
44  double MassConstraint::getValue() const
45  {
46  double totE[2] = {0, 0};
47  double totpx[2] = {0, 0};
48  double totpy[2] = {0, 0};
49  double totpz[2] = {0, 0};
50  for (unsigned int i = 0; i < fitobjects.size(); i++) {
51  int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
52  auto* pfo = dynamic_cast < ParticleFitObject* >(fitobjects[i]);
53  assert(pfo);
54  totE[index] += pfo->getE();
55  totpx[index] += pfo->getPx();
56  totpy[index] += pfo->getPy();
57  totpz[index] += pfo->getPz();
58  }
59  double result = -mass;
60  result += std::sqrt(std::abs(totE[0] * totE[0] - totpx[0] * totpx[0] - totpy[0] * totpy[0] - totpz[0] * totpz[0]));
61  result -= std::sqrt(std::abs(totE[1] * totE[1] - totpx[1] * totpx[1] - totpy[1] * totpy[1] - totpz[1] * totpz[1]));
62  return result;
63  }
64 
65 // calculate vector/array of derivatives of this contraint
66 // w.r.t. to ALL parameters of all fitobjects
67 // here: d M /d par(j)
68 // = d M /d p(i) * d p(i) /d par(j)
69 // = +-1/M * p(i) * d p(i) /d par(j)
70  void MassConstraint::getDerivatives(int idim, double der[]) const
71  {
72  double totE[2] = {0, 0};
73  double totpx[2] = {0, 0};
74  double totpy[2] = {0, 0};
75  double totpz[2] = {0, 0};
76  bool valid[2] = {false, false};
77  for (unsigned int i = 0; i < fitobjects.size(); i++) {
78  int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
79  valid[index] = true;
80  auto* pfo = dynamic_cast < ParticleFitObject* >(fitobjects[i]);
81  assert(pfo);
82  totE[index] += pfo->getE();
83  totpx[index] += pfo->getPx();
84  totpy[index] += pfo->getPy();
85  totpz[index] += pfo->getPz();
86  }
87  double m2[2];
88  double m_inv[2] = {0, 0};
89  for (int index = 0; index < 2; ++index) {
90  m2[index] = totE[index] * totE[index] - totpx[index] * totpx[index]
91  - totpy[index] * totpy[index] - totpz[index] * totpz[index];
92  if (m2[index] < 0 && m2[index] > -1E-9) m2[index] = 0;
93  if (m2[index] < 0 && valid[index]) {
94  B2ERROR("MassConstraint::getDerivatives: m2<0!");
95  for (unsigned int j = 0; j < fitobjects.size(); j++) {
96  int jndex = (flags[j] == 1) ? 0 : 1;
97  if (jndex == index) {
98 
99  auto* pfo = dynamic_cast < ParticleFitObject* >(fitobjects[j]);
100  assert(pfo);
101  B2ERROR(pfo->getName() <<
102  ": E =" << pfo->getE() <<
103  ", px=" << pfo->getPx() <<
104  ", py=" << pfo->getPy() <<
105  ", pz=" << pfo->getPz());
106  }
107  }
108  B2ERROR("sum: E=" << totE[index] << ", px=" << totpx[index]
109  << ", py=" << totpy[index] << ", pz=" << totpz[index] << ", m2=" << m2[index]);
110  }
111  if (m2[index] != 0) m_inv[index] = 1 / std::sqrt(std::abs(m2[index]));
112  }
113 
114  for (unsigned int i = 0; i < fitobjects.size(); i++) {
115  int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
116  for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
117  if (!fitobjects[i]->isParamFixed(ilocal)) {
118  int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
119  assert(iglobal >= 0 && iglobal < idim);
120  if (m2[index] != 0) {
121 
122  auto* pfo = dynamic_cast < ParticleFitObject* >(fitobjects[i]);
123  assert(pfo);
124 
125  der[iglobal] = totE[index] * pfo->getDE(ilocal)
126  - totpx[index] * pfo->getDPx(ilocal)
127  - totpy[index] * pfo->getDPy(ilocal)
128  - totpz[index] * pfo->getDPz(ilocal);
129  der[iglobal] *= m_inv[index];
130  } else der[iglobal] = 1;
131  if (index == 1) der[iglobal] *= -1.;
132  }
133  }
134  }
135  }
136 
137  double MassConstraint::getMass(int flag)
138  {
139  double totE = 0;
140  double totpx = 0;
141  double totpy = 0;
142  double totpz = 0;
143  for (unsigned int i = 0; i < fitobjects.size(); i++) {
144  if (flags[i] == flag) {
145 
146  const ParticleFitObject* fok = dynamic_cast < ParticleFitObject* >(fitobjects[i]);
147  assert(fok);
148  totE += fok->getE();
149  totpx += fok->getPx();
150  totpy += fok->getPy();
151  totpz += fok->getPz();
152  }
153  }
154  return std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
155  }
156 
157  void MassConstraint::setMass(double mass_)
158  {
159  mass = mass_;
160  }
161 
162  bool MassConstraint::secondDerivatives(int i, int j, double* dderivatives) const
163  {
164  B2DEBUG(14, "MassConstraint::secondDerivatives: i=" << i << ", j=" << j);
165  int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
166  int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
167  if (index != jndex) return false;
168  double totE = 0;
169  double totpx = 0;
170  double totpy = 0;
171  double totpz = 0;
172  for (unsigned int k = 0; k < fitobjects.size(); ++k) {
173  int kndex = (flags[k] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
174  const ParticleFitObject* fok = dynamic_cast < ParticleFitObject* >(fitobjects[k]);
175  assert(fok);
176  if (index == kndex) {
177  totE += fok->getE();
178  totpx += fok->getPx();
179  totpy += fok->getPy();
180  totpz += fok->getPz();
181  }
182  }
183 
184  if (totE <= 0) {
185  B2ERROR("MassConstraint::secondDerivatives: totE = " << totE);
186  }
187 
188  double m2 = std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz);
189  double m = std::sqrt(m2);
190  if (index) m = -m;
191  double minv3 = 1 / (m * m * m);
192 
193  assert(dderivatives);
194  for (int k = 0; k < 16; ++k) dderivatives[k] = 0;
195  dderivatives[4 * 0 + 0] = (m2 - totE * totE) * minv3;
196  dderivatives[4 * 0 + 1] = dderivatives[4 * 1 + 0] = totE * totpx * minv3;
197  dderivatives[4 * 0 + 2] = dderivatives[4 * 2 + 0] = totE * totpy * minv3;
198  dderivatives[4 * 0 + 3] = dderivatives[4 * 3 + 0] = totE * totpz * minv3;
199  dderivatives[4 * 1 + 1] = -(m2 + totpx * totpx) * minv3;
200  dderivatives[4 * 1 + 2] = dderivatives[4 * 2 + 1] = -totpx * totpy * minv3;
201  dderivatives[4 * 1 + 3] = dderivatives[4 * 3 + 1] = -totpx * totpz * minv3;
202  dderivatives[4 * 2 + 2] = -(m2 + totpy * totpy) * minv3;
203  dderivatives[4 * 2 + 3] = dderivatives[4 * 3 + 2] = -totpy * totpz * minv3;
204  dderivatives[4 * 3 + 3] = -(m2 + totpz * totpz) * minv3;
205  return true;
206  }
207 
208  bool MassConstraint::firstDerivatives(int i, double* dderivatives) const
209  {
210  double totE = 0;
211  double totpx = 0;
212  double totpy = 0;
213  double totpz = 0;
214  int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
215  for (unsigned int j = 0; j < fitobjects.size(); ++j) {
216  int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
217  const ParticleFitObject* foj = dynamic_cast < ParticleFitObject* >(fitobjects[j]);
218  assert(foj);
219  if (index == jndex) {
220  totE += foj->getE();
221  totpx += foj->getPx();
222  totpy += foj->getPy();
223  totpz += foj->getPz();
224  }
225  }
226 
227  if (totE <= 0) {
228  B2ERROR("MassConstraint::firstDerivatives: totE = " << totE);
229  }
230 
231  double m = std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
232  if (index) m = -m;
233 
234  dderivatives[0] = totE / m;
235  dderivatives[1] = -totpx / m;
236  dderivatives[2] = -totpy / m;
237  dderivatives[3] = -totpz / m;
238  return true;
239  }
240 
241  int MassConstraint::getVarBasis() const
242  {
243  return VAR_BASIS;
244  }
245  }// end OrcaKinFit namespace
247 } // end Belle2 namespace
Belle2::OrcaKinFit::ParticleFitObject::getE
virtual double getE() const
Return E.
Definition: ParticleFitObject.cc:106
Belle2::OrcaKinFit::MassConstraint::~MassConstraint
virtual ~MassConstraint()
Virtual destructor.
Belle2::OrcaKinFit::MassConstraint::getMass
virtual double getMass(int flag=1)
Get the actual invariant mass of the fit objects with a given flag.
Definition: MassConstraint.cc:151
Belle2::OrcaKinFit::ParticleFitObject::getPx
virtual double getPx() const
Return px.
Definition: ParticleFitObject.cc:110
Belle2::OrcaKinFit::MassConstraint::getValue
virtual double getValue() const override
Returns the value of the constraint.
Definition: MassConstraint.cc:58
Belle2::OrcaKinFit::MassConstraint::setMass
virtual void setMass(double mass_)
Sets the target mass of the constraint.
Definition: MassConstraint.cc:171
Belle2::OrcaKinFit::BaseHardConstraint::fitobjects
FitObjectContainer fitobjects
The FitObjectContainer.
Definition: BaseHardConstraint.h:191
Belle2::OrcaKinFit::MassConstraint::firstDerivatives
virtual bool firstDerivatives(int i, double *derivatives) const override
First derivatives with respect to the 4-vector of Fit objects i; result false if all derivatives are ...
Definition: MassConstraint.cc:222
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::OrcaKinFit::MassConstraint::MassConstraint
MassConstraint(double mass_=0.)
Constructor.
Definition: MassConstraint.cc:50
Belle2::OrcaKinFit::BaseHardConstraint::flags
std::vector< int > flags
The flags can be used to divide the FitObjectContainer into several subsets used for example to imple...
Definition: BaseHardConstraint.h:196
Belle2::OrcaKinFit::ParticleFitObject::getPy
virtual double getPy() const
Return py.
Definition: ParticleFitObject.cc:114
Belle2::OrcaKinFit::MassConstraint::mass
double mass
The mass difference between object sets 1 and 2.
Definition: MassConstraint.h:90
Belle2::OrcaKinFit::MassConstraint::getDerivatives
virtual void getDerivatives(int idim, double der[]) const override
Get first order derivatives.
Definition: MassConstraint.cc:84
Belle2::OrcaKinFit::MassConstraint::secondDerivatives
virtual bool secondDerivatives(int i, int j, double *derivatives) const override
Second derivatives with respect to the 4-vectors of Fit objects i and j; result false if all derivati...
Definition: MassConstraint.cc:176
Belle2::OrcaKinFit::ParticleFitObject
Definition: ParticleFitObject.h:93
Belle2::OrcaKinFit::ParticleFitObject::getPz
virtual double getPz() const
Return pz.
Definition: ParticleFitObject.cc:118