Belle II Software light-2406-ragdoll
MassConstraint.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * Forked from https://github.com/iLCSoft/MarlinKinfit *
6 * *
7 * Further information about the fit engine and the user interface *
8 * provided in MarlinKinfit can be found at *
9 * https://www.desy.de/~blist/kinfit/doc/html/ *
10 * and in the LCNotes LC-TOOL-2009-001 and LC-TOOL-2009-004 available *
11 * from http://www-flc.desy.de/lcnotes/ *
12 * *
13 * See git log for contributors and copyright holders. *
14 * This file is licensed under LGPL-3.0, see LICENSE.md. *
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
27namespace Belle2 {
32 namespace OrcaKinFit {
33
34
35// constructor
37 : mass(mass_)
38 {}
39
40// destructor
42
43// calculate current value of constraint function
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 constraint
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
FitObjectContainer fitobjects
The FitObjectContainer.
std::vector< int > flags
The flags can be used to divide the FitObjectContainer into several subsets used for example to imple...
virtual ~MassConstraint()
Virtual destructor.
virtual void getDerivatives(int idim, double der[]) const override
Get first order derivatives.
double mass
The mass difference between object sets 1 and 2.
MassConstraint(double mass_=0.)
Constructor.
virtual double getValue() const override
Returns the value of the constraint.
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...
virtual void setMass(double mass_)
Sets the target mass of the constraint.
virtual double getMass(int flag=1)
Get the actual invariant mass of the fit objects with a given flag.
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 ...
virtual double getPx() const
Return px.
virtual double getPz() const
Return pz.
virtual double getPy() const
Return py.
virtual double getE() const
Return E.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24