Belle II Software light-2406-ragdoll
SoftBWMassConstraint.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#ifdef MARLIN_USE_ROOT
18
19#include "analysis/OrcaKinFit/SoftBWMassConstraint.h"
20#include "analysis/OrcaKinFit/ParticleFitObject.h"
21#include <framework/logging/Logger.h>
22
23#include<iostream>
24#include<cmath>
25#include<limits>
26#undef NDEBUG
27#include<cassert>
28
29
30namespace Belle2 {
35 namespace OrcaKinFit {
36
37// constructor
38 SoftBWMassConstraint::SoftBWMassConstraint(double gamma_, double mass_, double massmin_, double massmax_)
39 : SoftBWParticleConstraint(gamma_, massmin_ - mass_, massmax_ - mass_),
40 mass(mass_)
41 {}
42
43// destructor
44 SoftBWMassConstraint::~SoftBWMassConstraint()
45 {}
46
47// calculate current value of constraint function
48 double SoftBWMassConstraint::getValue() const
49 {
50 double totE[2] = {0, 0};
51 double totpx[2] = {0, 0};
52 double totpy[2] = {0, 0};
53 double totpz[2] = {0, 0};
54 for (unsigned int i = 0; i < fitobjects.size(); i++) {
55 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
56 totE[index] += fitobjects[i]->getE();
57 totpx[index] += fitobjects[i]->getPx();
58 totpy[index] += fitobjects[i]->getPy();
59 totpz[index] += fitobjects[i]->getPz();
60 }
61 double m1 = std::sqrt(std::abs(totE[0] * totE[0] - totpx[0] * totpx[0] - totpy[0] * totpy[0] - totpz[0] * totpz[0]));
62 if (!std::isfinite(m1))
63 B2INFO("SoftBWMassConstraint::getValue(): m1 is nan: "
64 << "p[0]=(" << totE[0] << ", " << totpx[0] << ", " << totpy[0] << ", " << totpz[0]
65 << "), p^2=" << totE[0]*totE[0] - totpx[0]*totpx[0] - totpy[0]*totpy[0] - totpz[0]*totpz[0]);
66
67 assert(std::isfinite(m1));
68 double m2 = std::sqrt(std::abs(totE[1] * totE[1] - totpx[1] * totpx[1] - totpy[1] * totpy[1] - totpz[1] * totpz[1]));
69 if (!std::isfinite(m2))
70 B2INFO("SoftBWMassConstraint::getValue(): m2 is nan: "
71 << "p[1]=(" << totE[1] << ", " << totpx[1] << ", " << totpy[1] << ", " << totpz[1]
72 << "), p^2=" << totE[1]*totE[1] - totpx[1]*totpx[1] - totpy[1]*totpy[1] - totpz[1]*totpz[1]);
73 assert(std::isfinite(m2));
74 double result = m1 - m2 - mass;
75 assert(std::isfinite(result));
76 return result;
77 }
78
79// calculate vector/array of derivatives of this constraint
80// w.r.t. to ALL parameters of all fitobjects
81// here: d M /d par(j)
82// = d M /d p(i) * d p(i) /d par(j)
83// = +-1/M * p(i) * d p(i) /d par(j)
84 void SoftBWMassConstraint::getDerivatives(int idim, double der[]) const
85 {
86 double totE[2] = {0, 0};
87 double totpx[2] = {0, 0};
88 double totpy[2] = {0, 0};
89 double totpz[2] = {0, 0};
90 bool valid[2] = {false, false};
91 for (unsigned int i = 0; i < fitobjects.size(); i++) {
92 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
93 valid[index] = true;
94 totE[index] += fitobjects[i]->getE();
95 totpx[index] += fitobjects[i]->getPx();
96 totpy[index] += fitobjects[i]->getPy();
97 totpz[index] += fitobjects[i]->getPz();
98 }
99 double m2[2];
100 double m_inv[2] = {0, 0};
101 for (int index = 0; index < 2; ++index) {
102 m2[index] = totE[index] * totE[index] - totpx[index] * totpx[index]
103 - totpy[index] * totpy[index] - totpz[index] * totpz[index];
104 if (m2[index] < 0 && m2[index] > -1E-9) m2[index] = 0;
105 if (m2[index] < 0 && valid[index]) {
106 B2ERROR("SoftBWMassConstraint::getDerivatives: m2<0!");
107 for (unsigned int j = 0; j < fitobjects.size(); j++) {
108 int jndex = (flags[j] == 1) ? 0 : 1;
109 if (jndex == index) {
110 B2ERROR(fitobjects[j]->getName() << ": E=" << fitobjects[j]->getE() << ", px=" << fitobjects[j]->getPx()
111 << ", py=" << fitobjects[j]->getPy() << ", pz=" << fitobjects[j]->getPz());
112 }
113 }
114 B2ERROR("sum: E=" << totE[index] << ", px=" << totpx[index]
115 << ", py=" << totpy[index] << ", pz=" << totpz[index] << ", m2=" << m2[index]);
116 }
117 if (m2[index] != 0) m_inv[index] = 1 / std::sqrt(std::abs(m2[index]));
118 }
119
120 for (unsigned int i = 0; i < fitobjects.size(); i++) {
121 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
122 for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
123 if (!fitobjects[i]->isParamFixed(ilocal)) {
124 int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
125 assert(iglobal >= 0 && iglobal < idim);
126 if (m2[index] != 0) {
127 der[iglobal] = totE[index] * fitobjects[i]->getDE(ilocal)
128 - totpx[index] * fitobjects[i]->getDPx(ilocal)
129 - totpy[index] * fitobjects[i]->getDPy(ilocal)
130 - totpz[index] * fitobjects[i]->getDPz(ilocal);
131 der[iglobal] *= m_inv[index];
132 } else der[iglobal] = 1;
133 if (index == 1) der[iglobal] *= -1.;
134 }
135 }
136 }
137 }
138
139 double SoftBWMassConstraint::getMass(int flag)
140 {
141 double totE = 0;
142 double totpx = 0;
143 double totpy = 0;
144 double totpz = 0;
145 for (unsigned int i = 0; i < fitobjects.size(); i++) {
146 if (flags[i] == flag) {
147 totE += fitobjects[i]->getE();
148 totpx += fitobjects[i]->getPx();
149 totpy += fitobjects[i]->getPy();
150 totpz += fitobjects[i]->getPz();
151 }
152 }
153 return std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
154 }
155
156 void SoftBWMassConstraint::setMass(double mass_)
157 {
158 mass = mass_;
159 }
160
161 bool SoftBWMassConstraint::secondDerivatives(int i, int j, double* dderivatives) const
162 {
163 B2DEBUG(14, "SoftBWMassConstraint::secondDerivatives: i=" << i << ", j=" << j);
164 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
165 int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
166 if (index != jndex) return false;
167 double totE = 0;
168 double totpx = 0;
169 double totpy = 0;
170 double totpz = 0;
171 for (unsigned int k = 0; k < fitobjects.size(); ++k) {
172 int kndex = (flags[k] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
173 B2DEBUG(14, "SoftBWMassConstraint::secondDerivatives: i=" << i
174 << ", j=" << j << ", k=" << k << ", index=" << index << ", kndex=" << kndex);
175 const ParticleFitObject* fok = fitobjects[k];
176 assert(fok);
177 B2DEBUG(14, "fok: E=" << fok->getE());
178 if (index == kndex) {
179 totE += fok->getE();
180 totpx += fok->getPx();
181 totpy += fok->getPy();
182 totpz += fok->getPz();
183 }
184 }
185
186 if (totE <= 0) {
187 B2ERROR("SoftBWMassConstraint::secondDerivatives: totE = " << totE);
188 }
189
190 double m2 = std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz);
191 double m = std::sqrt(m2);
192 if (index) m = -m;
193 double minv3 = 1 / (m * m * m);
194
195 assert(dderivatives);
196 for (int k = 0; k < 16; ++k) dderivatives[k] = 0;
197 dderivatives[4 * 0 + 0] = (m2 - totE * totE) * minv3;
198 dderivatives[4 * 0 + 1] = dderivatives[4 * 1 + 0] = totE * totpx * minv3;
199 dderivatives[4 * 0 + 2] = dderivatives[4 * 2 + 0] = totE * totpy * minv3;
200 dderivatives[4 * 0 + 3] = dderivatives[4 * 3 + 0] = totE * totpz * minv3;
201 dderivatives[4 * 1 + 1] = -(m2 + totpx * totpx) * minv3;
202 dderivatives[4 * 1 + 2] = dderivatives[4 * 2 + 1] = -totpx * totpy * minv3;
203 dderivatives[4 * 1 + 3] = dderivatives[4 * 3 + 1] = -totpx * totpz * minv3;
204 dderivatives[4 * 2 + 2] = -(m2 + totpy * totpy) * minv3;
205 dderivatives[4 * 2 + 3] = dderivatives[4 * 3 + 2] = -totpy * totpz * minv3;
206 dderivatives[4 * 3 + 3] = -(m2 + totpz * totpz) * minv3;
207 return true;
208 }
209
210 bool SoftBWMassConstraint::firstDerivatives(int i, double* dderivatives) const
211 {
212 double totE = 0;
213 double totpx = 0;
214 double totpy = 0;
215 double totpz = 0;
216 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
217 for (unsigned int j = 0; j < fitobjects.size(); ++j) {
218 int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
219 B2DEBUG(14, "SoftBWMassConstraint::firstDerivatives: i=" << i
220 << ", j=" << j << ", index=" << index << ", jndex=" << jndex);
221 const ParticleFitObject* foj = fitobjects[j];
222 assert(foj);
223 B2DEBUG(14, "foj: E=" << foj->getE());
224 if (index == jndex) {
225 totE += foj->getE();
226 totpx += foj->getPx();
227 totpy += foj->getPy();
228 totpz += foj->getPz();
229 }
230 }
231
232 if (totE <= 0) {
233 B2INFO("SoftBWMassConstraint::firstDerivatives: totE = " << totE);
234 }
235
236 double m = std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
237 if (index) m = -m;
238
239 dderivatives[0] = totE / m;
240 dderivatives[1] = -totpx / m;
241 dderivatives[2] = -totpy / m;
242 dderivatives[3] = -totpz / m;
243 return true;
244 }
245 }// end OrcaKinFit namespace
247} // end Belle2 namespace
248
249#endif // MARLIN_USE_ROOT
250
Belle2
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24