Belle II Software light-2406-ragdoll
RecoilMassConstraint.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include "analysis/OrcaKinFit/RecoilMassConstraint.h"
10#include "analysis/OrcaKinFit/ParticleFitObject.h"
11
12#include<iostream>
13#include<cmath>
14
15#undef NDEBUG
16#include<cassert>
17
18namespace Belle2 {
23 namespace OrcaKinFit {
24
25// constructor
26 RecoilMassConstraint::RecoilMassConstraint(double recoilmass, double beampx, double beampy, double beampz, double beampe)
27 : m_recoilMass(recoilmass), m_beamPx(beampx), m_beamPy(beampy), m_beamPz(beampz), m_beamE(beampe)
28 {}
29
30// destructor
32 {
33 ;
34 }
35
36// calculate current value of constraint function
38 {
39 double totE = 0.;
40 double totpx = 0.;
41 double totpy = 0.;
42 double totpz = 0.;
43
44 for (auto fitobject : fitobjects) {
45 auto* pfo = dynamic_cast < ParticleFitObject* >(fitobject);
46 assert(pfo);
47 totE += pfo->getE();
48 totpx += pfo->getPx();
49 totpy += pfo->getPy();
50 totpz += pfo->getPz();
51 }
52
53 const double recoilE = (m_beamE - totE);
54 const double recoilE2 = recoilE * recoilE;
55 const double recoilpx = (m_beamPx - totpx);
56 const double recoilpx2 = recoilpx * recoilpx;
57 const double recoilpy = (m_beamPy - totpy);
58 const double recoilpy2 = recoilpy * recoilpy;
59 const double recoilpz = (m_beamPz - totpz);
60 const double recoilpz2 = recoilpz * recoilpz;
61 const double recoil2 = recoilE2 - recoilpx2 - recoilpy2 - recoilpz2;
62 const double result = recoil2 - m_recoilMass * m_recoilMass;
63 return result;
64 }
65
66// calculate vector/array of derivatives of this constraint
67// w.r.t. to ALL parameters of all fitobjects
68// here: d RM2 /d par(j)
69// = d RM2 /d p(i) * d p(i) /d par(j)
70// = -+2 * recoil_p(i) * d p(i) /d par(j)
71
72 void RecoilMassConstraint::getDerivatives(int idim, double der[]) const
73 {
74
75 double totE = 0.;
76 double totpx = 0.;
77 double totpy = 0.;
78 double totpz = 0.;
79
80 for (auto fitobject : fitobjects) {
81 auto* pfo = dynamic_cast < ParticleFitObject* >(fitobject);
82 assert(pfo);
83 totE += pfo->getE();
84 totpx += pfo->getPx();
85 totpy += pfo->getPy();
86 totpz += pfo->getPz();
87 }
88
89 const double recoilE = (m_beamE - totE);
90 const double recoilpx = (m_beamPx - totpx);
91 const double recoilpy = (m_beamPy - totpy);
92 const double recoilpz = (m_beamPz - totpz);
93
94 for (auto fitobject : fitobjects) {
95 for (int ilocal = 0; ilocal < fitobject->getNPar(); ilocal++) {
96 if (!fitobject->isParamFixed(ilocal)) {
97 int iglobal = fitobject->getGlobalParNum(ilocal);
98 assert(iglobal >= 0 && iglobal < idim);
99
100 auto* pfo = dynamic_cast < ParticleFitObject* >(fitobject);
101 assert(pfo);
102 der[iglobal] = - recoilE * pfo->getDE(ilocal)
103 + recoilpx * pfo->getDPx(ilocal)
104 + recoilpy * pfo->getDPy(ilocal)
105 + recoilpz * pfo->getDPz(ilocal);
106 der[iglobal] *= 2;
107 }
108 }
109 }
110 }
111
112
114 {
115 double totE = 0.;
116 double totpx = 0.;
117 double totpy = 0.;
118 double totpz = 0.;
119
120 for (auto& fitobject : fitobjects) {
121 auto* pfo = dynamic_cast < ParticleFitObject* >(fitobject);
122 assert(pfo);
123 totE += pfo->getE();
124 totpx += pfo->getPx();
125 totpy += pfo->getPy();
126 totpz += pfo->getPz();
127 }
128
129 const double recoilE = (m_beamE - totE);
130 const double recoilE2 = recoilE * recoilE;
131 const double recoilpx = (m_beamPx - totpx);
132 const double recoilpx2 = recoilpx * recoilpx;
133 const double recoilpy = (m_beamPy - totpy);
134 const double recoilpy2 = recoilpy * recoilpy;
135 const double recoilpz = (m_beamPz - totpz);
136 const double recoilpz2 = recoilpz * recoilpz;
137 const double recoil2 = recoilE2 - recoilpx2 - recoilpy2 - recoilpz2;
138 const double recoil = std::sqrt(std::fabs(recoil2));
139
140 return recoil;
141 }
142
143 void RecoilMassConstraint::setRecoilMass(double recoilmass_)
144 {
145 m_recoilMass = recoilmass_;
146 }
147
148 bool RecoilMassConstraint::secondDerivatives(int i, int j, double* dderivatives) const
149 {
150 (void) i;
151 (void) j;
152
153 assert(dderivatives);
154 for (int k = 0; k < 16; ++k) dderivatives[k] = 0;
155
156
157 dderivatives[4 * 0 + 0] = 2; //dE^2
158 dderivatives[4 * 0 + 1] = 0; //dEdpx
159 dderivatives[4 * 0 + 2] = 0; //dEdpy
160 dderivatives[4 * 0 + 3] = 0; //dEdpz
161 dderivatives[4 * 1 + 1] = -2;//dpx^2
162 dderivatives[4 * 1 + 2] = 0; //dpxdpy
163 dderivatives[4 * 1 + 3] = 0; //dpxdpz
164 dderivatives[4 * 2 + 2] = -2;//dpy^2
165 dderivatives[4 * 2 + 3] = 0; //dpydpz
166 dderivatives[4 * 3 + 3] = -2;//dpz^2
167
168 return true;
169
170 }
171
172 bool RecoilMassConstraint::firstDerivatives(int i, double* dderivatives) const
173 {
174 (void) i;
175
176 double totE = 0.;
177 double totpx = 0.;
178 double totpy = 0.;
179 double totpz = 0.;
180
181 for (auto fitobject : fitobjects) {
182 auto* pfo = dynamic_cast < ParticleFitObject* >(fitobject);
183 assert(pfo);
184 totE += pfo->getE();
185 totpx += pfo->getPx();
186 totpy += pfo->getPy();
187 totpz += pfo->getPz();
188 }
189
190 const double recoilE = (m_beamE - totE);
191 const double recoilpx = (m_beamPx - totpx);
192 const double recoilpy = (m_beamPy - totpy);
193 const double recoilpz = (m_beamPz - totpz);
194
195 dderivatives[0] = -2 * recoilE;
196 dderivatives[1] = 2 * recoilpx;
197 dderivatives[2] = 2 * recoilpy;
198 dderivatives[3] = 2 * recoilpz;
199
200 return true;
201 }
202
203 int RecoilMassConstraint::getVarBasis() const
204 {
205 return VAR_BASIS;
206 }
207 }// end OrcaKinFit namespace
209} // end Belle2 namespace
210
FitObjectContainer fitobjects
The FitObjectContainer.
virtual void setRecoilMass(double recoilmass)
Sets the target recoil mass of the constraint.
virtual void getDerivatives(int idim, double der[]) const override
Get first order derivatives.
virtual double getRecoilMass()
Get the actual recoil mass of the fit objects.
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...
RecoilMassConstraint(double recoilmass=0., double beampx=0., double beampy=0., double beampz=0, double beampe=0.)
Constructor.
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 ~RecoilMassConstraint()
Virtual destructor.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24