Belle II Software light-2406-ragdoll
SoftGaussMassConstraint.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/SoftGaussMassConstraint.h"
18#include "analysis/OrcaKinFit/ParticleFitObject.h"
19#include <framework/logging/Logger.h>
20
21#include<iostream>
22#include<cmath>
23#undef NDEBUG
24#include<cassert>
25
26
27namespace Belle2 {
32 namespace OrcaKinFit {
33
34// 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 totE[index] += fitobjects[i]->getE();
53 totpx[index] += fitobjects[i]->getPx();
54 totpy[index] += fitobjects[i]->getPy();
55 totpz[index] += fitobjects[i]->getPz();
56 }
57 double result = -mass;
58 result += std::sqrt(std::abs(totE[0] * totE[0] - totpx[0] * totpx[0] - totpy[0] * totpy[0] - totpz[0] * totpz[0]));
59 result -= std::sqrt(std::abs(totE[1] * totE[1] - totpx[1] * totpx[1] - totpy[1] * totpy[1] - totpz[1] * totpz[1]));
60 return result;
61 }
62
63// calculate vector/array of derivatives of this constraint
64// w.r.t. to ALL parameters of all fitobjects
65// here: d M /d par(j)
66// = d M /d p(i) * d p(i) /d par(j)
67// = +-1/M * p(i) * d p(i) /d par(j)
68 void SoftGaussMassConstraint::getDerivatives(int idim, double der[]) const
69 {
70 double totE[2] = {0, 0};
71 double totpx[2] = {0, 0};
72 double totpy[2] = {0, 0};
73 double totpz[2] = {0, 0};
74 bool valid[2] = {false, false};
75 for (unsigned int i = 0; i < fitobjects.size(); i++) {
76 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
77 valid[index] = true;
78 totE[index] += fitobjects[i]->getE();
79 totpx[index] += fitobjects[i]->getPx();
80 totpy[index] += fitobjects[i]->getPy();
81 totpz[index] += fitobjects[i]->getPz();
82 }
83 double m2[2];
84 double m_inv[2] = {0, 0};
85 for (int index = 0; index < 2; ++index) {
86 m2[index] = totE[index] * totE[index] - totpx[index] * totpx[index]
87 - totpy[index] * totpy[index] - totpz[index] * totpz[index];
88 if (m2[index] < 0 && m2[index] > -1E-9) m2[index] = 0;
89 if (m2[index] < 0 && valid[index]) {
90 B2ERROR("SoftGaussMassConstraint::getDerivatives: m2<0!");
91 for (unsigned int j = 0; j < fitobjects.size(); j++) {
92 int jndex = (flags[j] == 1) ? 0 : 1;
93 if (jndex == index) {
94 B2ERROR(fitobjects[j]->getName() << ": E=" << fitobjects[j]->getE() << ", px=" << fitobjects[j]->getPx()
95 << ", py=" << fitobjects[j]->getPy() << ", pz=" << fitobjects[j]->getPz());
96 }
97 }
98 B2ERROR("sum: E=" << totE[index] << ", px=" << totpx[index]
99 << ", py=" << totpy[index] << ", pz=" << totpz[index] << ", m2=" << m2[index]);
100 }
101 if (m2[index] != 0) m_inv[index] = 1 / std::sqrt(std::abs(m2[index]));
102 }
103
104 for (unsigned int i = 0; i < fitobjects.size(); i++) {
105 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
106 for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
107 if (!fitobjects[i]->isParamFixed(ilocal)) {
108 int iglobal = fitobjects[i]->getGlobalParNum(ilocal);
109 assert(iglobal >= 0 && iglobal < idim);
110 if (m2[index] != 0) {
111 der[iglobal] = totE[index] * fitobjects[i]->getDE(ilocal)
112 - totpx[index] * fitobjects[i]->getDPx(ilocal)
113 - totpy[index] * fitobjects[i]->getDPy(ilocal)
114 - totpz[index] * fitobjects[i]->getDPz(ilocal);
115 der[iglobal] *= m_inv[index];
116 } else der[iglobal] = 1;
117 if (index == 1) der[iglobal] *= -1.;
118 }
119 }
120 }
121 }
122
124 {
125 double totE = 0;
126 double totpx = 0;
127 double totpy = 0;
128 double totpz = 0;
129 for (unsigned int i = 0; i < fitobjects.size(); i++) {
130 if (flags[i] == flag) {
131 totE += fitobjects[i]->getE();
132 totpx += fitobjects[i]->getPx();
133 totpy += fitobjects[i]->getPy();
134 totpz += fitobjects[i]->getPz();
135 }
136 }
137 return std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
138 }
139
141 {
142 mass = mass_;
143 }
144
145 bool SoftGaussMassConstraint::secondDerivatives(int i, int j, double* dderivatives) const
146 {
147 B2DEBUG(14, "SoftGaussMassConstraint::secondDerivatives: i=" << i << ", j=" << j);
148 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
149 int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
150 if (index != jndex) return false;
151 double totE = 0;
152 double totpx = 0;
153 double totpy = 0;
154 double totpz = 0;
155 for (unsigned int k = 0; k < fitobjects.size(); ++k) {
156 int kndex = (flags[k] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
157 const ParticleFitObject* fok = fitobjects[k];
158 assert(fok);
159 if (index == kndex) {
160 totE += fok->getE();
161 totpx += fok->getPx();
162 totpy += fok->getPy();
163 totpz += fok->getPz();
164 }
165 }
166
167 if (totE <= 0) {
168 B2ERROR("SoftGaussMassConstraint::secondDerivatives: totE = " << totE);
169 }
170
171 double m2 = std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz);
172 double m = std::sqrt(m2);
173 if (index) m = -m;
174 double minv3 = 1 / (m * m * m);
175
176 assert(dderivatives);
177 for (int k = 0; k < 16; ++k) dderivatives[k] = 0;
178 dderivatives[4 * 0 + 0] = (m2 - totE * totE) * minv3;
179 dderivatives[4 * 0 + 1] = dderivatives[4 * 1 + 0] = totE * totpx * minv3;
180 dderivatives[4 * 0 + 2] = dderivatives[4 * 2 + 0] = totE * totpy * minv3;
181 dderivatives[4 * 0 + 3] = dderivatives[4 * 3 + 0] = totE * totpz * minv3;
182 dderivatives[4 * 1 + 1] = -(m2 + totpx * totpx) * minv3;
183 dderivatives[4 * 1 + 2] = dderivatives[4 * 2 + 1] = -totpx * totpy * minv3;
184 dderivatives[4 * 1 + 3] = dderivatives[4 * 3 + 1] = -totpx * totpz * minv3;
185 dderivatives[4 * 2 + 2] = -(m2 + totpy * totpy) * minv3;
186 dderivatives[4 * 2 + 3] = dderivatives[4 * 3 + 2] = -totpy * totpz * minv3;
187 dderivatives[4 * 3 + 3] = -(m2 + totpz * totpz) * minv3;
188 return true;
189 }
190
191 bool SoftGaussMassConstraint::firstDerivatives(int i, double* dderivatives) const
192 {
193 double totE = 0;
194 double totpx = 0;
195 double totpy = 0;
196 double totpz = 0;
197 int index = (flags[i] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
198 for (unsigned int j = 0; j < fitobjects.size(); ++j) {
199 int jndex = (flags[j] == 1) ? 0 : 1; // default is 1, but 2 may indicate fitobjects for a second W -> equal mass constraint!
200 const ParticleFitObject* foj = fitobjects[j];
201 assert(foj);
202 if (index == jndex) {
203 totE += foj->getE();
204 totpx += foj->getPx();
205 totpy += foj->getPy();
206 totpz += foj->getPz();
207 }
208 }
209
210 if (totE <= 0) {
211 B2ERROR("SoftGaussMassConstraint::firstDerivatives: totE = " << totE);
212 }
213
214 double m = std::sqrt(std::abs(totE * totE - totpx * totpx - totpy * totpy - totpz * totpz));
215 if (index) m = -m;
216
217 dderivatives[0] = totE / m;
218 dderivatives[1] = -totpx / m;
219 dderivatives[2] = -totpy / m;
220 dderivatives[3] = -totpz / m;
221 return true;
222 }
223 }// end OrcaKinFit namespace
225} // end Belle2 namespace
virtual const char * getName() const
Returns the name of the constraint.
virtual double getPx() const
Return px.
virtual double getPz() const
Return pz.
virtual double getPy() const
Return py.
virtual double getE() const
Return E.
virtual void getDerivatives(int idim, double der[]) const override
Get first order derivatives.
double mass
The mass difference between object sets 1 and 2.
virtual ~SoftGaussMassConstraint()
Virtual destructor.
virtual double getValue() const override
Returns the value of the constraint function.
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.
SoftGaussMassConstraint(double sigma_, double mass_=0.)
Constructor.
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 ...
Abstract base class for constraints of kinematic fits.
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...
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24