Belle II Software development
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | Related Functions | List of all members
MassConstraint Class Reference

Implements constraint 0 = mass1 - mass2 - m. More...

#include <MassConstraint.h>

Inheritance diagram for MassConstraint:
ParticleConstraint BaseHardConstraint BaseConstraint

Public Member Functions

 MassConstraint (double mass_=0.)
 Constructor.
 
virtual ~MassConstraint ()
 Virtual destructor.
 
virtual double getValue () const override
 Returns the value of the constraint.
 
virtual void getDerivatives (int idim, double der[]) const override
 Get first order derivatives.
 
virtual double getMass (int flag=1)
 Get the actual invariant mass of the fit objects with a given flag.
 
virtual void setMass (double mass_)
 Sets the target mass of the constraint.
 
virtual int getVarBasis () const override
 
virtual void setFOList (std::vector< ParticleFitObject * > *fitobjects_)
 Adds several ParticleFitObject objects to the list.
 
virtual void addToFOList (ParticleFitObject &fitobject, int flag=1)
 Adds one ParticleFitObject objects to the list.
 
virtual void resetFOList ()
 Resests ParticleFitObject list.
 
virtual void invalidateCache () const
 Invalidates any cached values for the next event.
 
virtual void add1stDerivativesToMatrix (double *M, int idim) const
 Adds first order derivatives to global covariance matrix M.
 
virtual void add2ndDerivativesToMatrix (double *M, int idim, double lambda) const
 Adds second order derivatives to global covariance matrix M.
 
virtual void addToGlobalChi2DerVector (double *y, int idim, double lambda) const
 Add lambda times derivatives of chi squared to global derivative vector.
 
virtual double dirDer (double *p, double *w, int idim, double mu=1)
 Calculate directional derivative.
 
virtual double dirDerAbs (double *p, double *w, int idim, double mu=1)
 Calculate directional derivative for abs(c)
 
virtual double getError () const override
 Returns the error on the value of the constraint.
 
virtual int getGlobalNum () const
 Accesses position of constraint in global constraint list.
 
virtual void setGlobalNum (int iglobal)
 Sets position of constraint in global constraint list.
 
virtual void printFirstDerivatives () const
 
virtual void printSecondDerivatives () const
 
virtual void test1stDerivatives ()
 
virtual void test2ndDerivatives ()
 
virtual double num1stDerivative (int ifo, int ilocal, double eps)
 Evaluates numerically the 1st derivative w.r.t. a parameter.
 
virtual double num2ndDerivative (int ifo1, int ilocal1, double eps1, int ifo2, int ilocal2, double eps2)
 Evaluates numerically the 2nd derivative w.r.t. 2 parameters.
 
virtual const char * getName () const
 Returns the name of the constraint.
 
void setName (const char *name_)
 Set object's name.
 
virtual std::ostream & print (std::ostream &os) const
 print object to ostream
 

Protected Types

enum  { VAR_BASIS = 0 }
 
typedef std::vector< BaseFitObject * > FitObjectContainer
 Vector of pointers to ParticleFitObjects.
 
typedef FitObjectContainer::iterator FitObjectIterator
 Iterator through vector of pointers to ParticleFitObjects.
 
typedef FitObjectContainer::const_iterator ConstFitObjectIterator
 Constant iterator through vector of pointers to ParticleFitObjects.
 

Protected Member Functions

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 derivatives are zero.
 
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 zero.
 

Protected Attributes

double mass
 The mass difference between object sets 1 and 2.
 
FitObjectContainer fitobjects
 The FitObjectContainer.
 
std::vector< double > derivatives
 The derivatives.
 
std::vector< int > flags
 The flags can be used to divide the FitObjectContainer into several subsets used for example to implement an equal mass constraint (see MassConstraint).
 
int globalNum
 Position of constraint in global constraint list.
 
char * name
 

Related Functions

(Note that these are not member functions.)

std::ostream & operator<< (std::ostream &os, const BaseConstraint &bc)
 Prints out a BaseConstraint, using its print method.
 

Detailed Description

Implements constraint 0 = mass1 - mass2 - m.

This class implements different mass constraints:

  • the invariant mass of several objects should be m
  • the difference of the invariant masses between two sets of objects should be m (normally m=0 in this case).

Author: Jenny List, Benno List Last update:

Date
2008/02/12 11:03:32

by:

Author
blist

Definition at line 46 of file MassConstraint.h.

Member Typedef Documentation

◆ ConstFitObjectIterator

typedef FitObjectContainer::const_iterator ConstFitObjectIterator
protectedinherited

Constant iterator through vector of pointers to ParticleFitObjects.

Definition at line 175 of file BaseHardConstraint.h.

◆ FitObjectContainer

typedef std::vector<BaseFitObject*> FitObjectContainer
protectedinherited

Vector of pointers to ParticleFitObjects.

Definition at line 171 of file BaseHardConstraint.h.

◆ FitObjectIterator

typedef FitObjectContainer::iterator FitObjectIterator
protectedinherited

Iterator through vector of pointers to ParticleFitObjects.

Definition at line 173 of file BaseHardConstraint.h.

Member Enumeration Documentation

◆ anonymous enum

anonymous enum
protected

Definition at line 89 of file MassConstraint.h.

89{ VAR_BASIS = 0 }; // this means that the constraint knows about E,px,py,pz

Constructor & Destructor Documentation

◆ MassConstraint()

MassConstraint ( double  mass_ = 0.)
explicit

Constructor.

Parameters
mass_The mass difference between object sets 1 and 2

Definition at line 36 of file MassConstraint.cc.

37 : mass(mass_)
38 {}
Belle2::OrcaKinFit::MassConstraint::mass
double mass
The mass difference between object sets 1 and 2.
Definition: MassConstraint.h:76

Member Function Documentation

◆ add1stDerivativesToMatrix()

void add1stDerivativesToMatrix ( double *  M,
int  idim 
) const
virtualinherited

Adds first order derivatives to global covariance matrix M.

Parameters
MGlobal covariance matrix, dimension at least idim x idim
idimFirst dimension of array der

Definition at line 37 of file BaseHardConstraint.cc.

38 {
39 double dgdpi[BaseDefs::MAXINTERVARS];
40 for (unsigned int i = 0; i < fitobjects.size(); ++i) {
41 const BaseFitObject* foi = fitobjects[i];
42 assert(foi);
43 if (firstDerivatives(i, dgdpi)) {
44 foi->addTo1stDerivatives(M, idim, dgdpi, getGlobalNum(), getVarBasis());
45 }
46 }
47 }
Belle2::OrcaKinFit::BaseHardConstraint::getGlobalNum
virtual int getGlobalNum() const
Accesses position of constraint in global constraint list.
Definition: BaseHardConstraint.h:137
Belle2::OrcaKinFit::BaseHardConstraint::fitobjects
FitObjectContainer fitobjects
The FitObjectContainer.
Definition: BaseHardConstraint.h:177
Belle2::OrcaKinFit::BaseHardConstraint::firstDerivatives
virtual bool firstDerivatives(int i, double *derivatives) const =0
First derivatives with respect to the meta-variables of Fit objects i; result false if all derivative...

◆ add2ndDerivativesToMatrix()

void add2ndDerivativesToMatrix ( double *  M,
int  idim,
double  lambda 
) const
virtualinherited

Adds second order derivatives to global covariance matrix M.

Calculates the second derivative of the constraint g w.r.t.

the various parameters, multiplies it by lambda and adds it to the global covariance matrix

in case of particlefitobject: We denote with P_i the 4-vector of the i-th ParticleFitObject, then

\[ \frac{\partial ^2 g}{\partial a_k \partial a_l} = \sum_i \sum_j \frac{\partial ^2 g}{\partial P_i \partial P_j} \cdot \frac{\partial P_i}{\partial a_k} \cdot \frac{\partial P_j}{\partial a_l} + \sum_i \frac{\partial g}{\partial P_i} \cdot \frac{\partial^2 P_i}{\partial a_k \partial a_l} \]

Here, $\frac{\partial P_i}{\partial a_k}$ is a $4 \times n_i$ Matrix, where $n_i$ is the number of parameters of FitObject i; Correspondingly, $\frac{\partial^2 P_i}{\partial a_k \partial a_l}$ is a $4 \times n_i \times n_i$ matrix. Also, $\frac{\partial ^2 g}{\partial P_i \partial P_j}$ is a $4\times 4$ matrix for a given i and j, and $\frac{\partial g}{\partial P_i}$ is a 4-vector (though not a Lorentz-vector!).

but here it's been generalised

First, treat the part

\[ \frac{\partial ^2 g}{\partial P_i \partial P_j} \cdot \frac{\partial P_i}{\partial a_k} \cdot \frac{\partial P_j}{\partial a_l} \]

Second, treat the part

\[ \sum_i \frac{\partial g}{\partial P_i} \cdot \frac{\partial^2 P_i}{\partial a_k \partial a_l} \]

Here, $\frac{\partial g}{\partial P_i}$ is a 4-vector, which we pass on to the FitObject

Parameters
MGlobal covariance matrix, dimension at least idim x idim
idimFirst dimension of array der
lambdaLagrange multiplier for this constraint

Definition at line 76 of file BaseHardConstraint.cc.

77 {
78
85 // Derivatives \f$\frac{\partial ^2 g}{\partial P_i \partial P_j}\f$ at fixed i, j
86 // d2GdPidPj[4*ii+jj] is derivative w.r.t. P_i,ii and P_j,jj, where ii=0,1,2,3 for E,px,py,pz
87 double d2GdPidPj[BaseDefs::MAXINTERVARS * BaseDefs::MAXINTERVARS];
88
89 // Derivatives \f$\frac {\partial P_i}{\partial a_k}\f$ for all i;
90 // k is local parameter number
91 // dPidAk[KMAX*4*i + 4*k + ii] is \f$\frac {\partial P_{i,ii}}{\partial a_k}\f$,
92 // with ii=0, 1, 2, 3 for E, px, py, pz
93 const int n = fitobjects.size();
94 auto* dPidAk = new double[n * BaseDefs::MAXPAR * BaseDefs::MAXINTERVARS];
95 bool* dPidAkval = new bool[n];
96
97 for (int i = 0; i < n; ++i) dPidAkval[i] = false;
98
99 // Derivatives \f$\frac{\partial ^2 g}{\partial P_i \partial a_l}\f$ at fixed i
100 // d2GdPdAl[4*l + ii] is \f$\frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l}\f$
101 double d2GdPdAl[static_cast<int>(BaseDefs::MAXINTERVARS) * BaseDefs::MAXPAR];
102 // Derivatives \f$\frac{\partial ^2 g}{\partial a_k \partial a_l}\f$
103 double d2GdAkdAl[BaseDefs::MAXPAR * BaseDefs::MAXPAR] = {0};
104
105 // Global parameter numbers: parglobal[BaseDefs::MAXPAR*i+klocal]
106 // is global parameter number of local parameter klocal of i-th Fit object
107 int* parglobal = new int[BaseDefs::MAXPAR * n];
108
109 for (int i = 0; i < n; ++i) {
110 const BaseFitObject* foi = fitobjects[i];
111 assert(foi);
112 for (int klocal = 0; klocal < foi->getNPar(); ++klocal) {
113 parglobal [BaseDefs::MAXPAR * i + klocal] = foi->getGlobalParNum(klocal);
114 }
115 }
116
117
118 for (int i = 0; i < n; ++i) {
119 const BaseFitObject* foi = fitobjects[i];
120 assert(foi);
121 for (int j = 0; j < n; ++j) {
122 const BaseFitObject* foj = fitobjects[j];
123 assert(foj);
124 if (secondDerivatives(i, j, d2GdPidPj)) {
125 if (!dPidAkval[i]) {
126 foi->getDerivatives(dPidAk + i * (static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS),
127 static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS);
128 dPidAkval[i] = true;
129 }
130 if (!dPidAkval[j]) {
131 foj->getDerivatives(dPidAk + j * (static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS),
132 static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS);
133 dPidAkval[j] = true;
134 }
135 // Now sum over E/px/Py/Pz for object j:
136 // \f[
137 // \frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l}
138 // = (sum_{j}) sum_{jj} frac{\partial ^2 g}{\partial P_{i,ii} \partial P_{j,jj}}
139 // \cdot \frac{\partial P_{j,jj}}{\partial a_l}
140 // \f]
141 // We're summing over jj here
142 for (int llocal = 0; llocal < foj->getNPar(); ++llocal) {
143 for (int ii = 0; ii < BaseDefs::MAXINTERVARS; ++ii) {
144 int ind1 = BaseDefs::MAXINTERVARS * ii;
145 int ind2 = (static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS) * j + BaseDefs::MAXINTERVARS * llocal;
146 double& r = d2GdPdAl[BaseDefs::MAXINTERVARS * llocal + ii];
147 r = d2GdPidPj[ ind1] * dPidAk[ ind2]; // E
148 r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // px
149 r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // py
150 r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // pz
151 }
152 }
153 // Now sum over E/px/Py/Pz for object i, i.e. sum over ii:
154 // \f[
155 // \frac{\partial ^2 g}{\partial a_k \partial a_l}
156 // = \sum_{ii} \frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l} \cdot
157 // \frac{\partial P_{i,ii}}{\partial a_k}
158 // \f]
159 for (int klocal = 0; klocal < foi->getNPar(); ++klocal) {
160 for (int llocal = 0; llocal < foj->getNPar(); ++llocal) {
161 int ind1 = BaseDefs::MAXINTERVARS * llocal;
162 int ind2 = (static_cast<int>(BaseDefs::MAXPAR) * BaseDefs::MAXINTERVARS) * i + BaseDefs::MAXINTERVARS * klocal;
163 double& r = d2GdAkdAl[BaseDefs::MAXPAR * klocal + llocal];
164 r = d2GdPdAl[ ind1] * dPidAk[ ind2]; //E
165 r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // px
166 r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // py
167 r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // pz
168 }
169 }
170 // Now expand the local parameter numbers to global ones
171 for (int klocal = 0; klocal < foi->getNPar(); ++klocal) {
172 int kglobal = parglobal [BaseDefs::MAXPAR * i + klocal];
173 for (int llocal = 0; llocal < foj->getNPar(); ++llocal) {
174 int lglobal = parglobal [BaseDefs::MAXPAR * j + llocal];
175 M [idim * kglobal + lglobal] += lambda * d2GdAkdAl[BaseDefs::MAXPAR * klocal + llocal];
176 }
177 }
178 }
179 }
180 }
190 double dgdpi[BaseDefs::MAXINTERVARS];
191 for (int i = 0; i < n; ++i) {
192 const BaseFitObject* foi = fitobjects[i];
193 assert(foi);
194 if (firstDerivatives(i, dgdpi)) {
195 foi->addTo2ndDerivatives(M, idim, lambda, dgdpi, getVarBasis());
196 }
197 }
198
199 delete[] dPidAk;
200 delete[] dPidAkval;
201 delete[] parglobal;
202 }
Belle2::OrcaKinFit::BaseHardConstraint::secondDerivatives
virtual bool secondDerivatives(int i, int j, double *derivatives) const =0
Second derivatives with respect to the meta-variables of Fit objects i and j; result false if all der...

◆ addToFOList()

virtual void addToFOList ( ParticleFitObject &  fitobject,
int  flag = 1 
)
inlinevirtualinherited

Adds one ParticleFitObject objects to the list.

Definition at line 90 of file ParticleConstraint.h.

92 {
93 fitobjects.push_back(reinterpret_cast < BaseFitObject* >(&fitobject));
94 flags.push_back(flag);
95 };
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:182

◆ addToGlobalChi2DerVector()

void addToGlobalChi2DerVector ( double *  y,
int  idim,
double  lambda 
) const
virtualinherited

Add lambda times derivatives of chi squared to global derivative vector.

Parameters
yVector of chi2 derivatives
idimVector size

Definition at line 204 of file BaseHardConstraint.cc.

205 {
206 double dgdpi[BaseDefs::MAXINTERVARS];
207 for (unsigned int i = 0; i < fitobjects.size(); ++i) {
208 const BaseFitObject* foi = fitobjects[i];
209 assert(foi);
210 if (firstDerivatives(i, dgdpi)) {
211 foi->addToGlobalChi2DerVector(y, idim, lambda, dgdpi, getVarBasis());
212 }
213 }
214 }

◆ dirDer()

double dirDer ( double *  p,
double *  w,
int  idim,
double  mu = 1 
)
virtualinherited

Calculate directional derivative.

Parameters
pVector of direction
wWork vector
idimVector size
muoptional multiplier

Definition at line 232 of file BaseHardConstraint.cc.

233 {
234 double* pw, *pp;
235 for (pw = w; pw < w + idim; * (pw++) = 0);
236 addToGlobalChi2DerVector(w, idim, mu);
237 double result = 0;
238 for (pw = w, pp = p; pw < w + idim; result += *(pp++) * *(pw++));
239 return mu * result;
240 }
Belle2::OrcaKinFit::BaseHardConstraint::addToGlobalChi2DerVector
virtual void addToGlobalChi2DerVector(double *y, int idim, double lambda) const
Add lambda times derivatives of chi squared to global derivative vector.
Definition: BaseHardConstraint.cc:204

◆ dirDerAbs()

double dirDerAbs ( double *  p,
double *  w,
int  idim,
double  mu = 1 
)
virtualinherited

Calculate directional derivative for abs(c)

Parameters
pVector of direction
wWork vector
idimVector size
muoptional multiplier

Definition at line 242 of file BaseHardConstraint.cc.

243 {
244 double val = getValue();
245 if (val == 0) return mu * std::fabs(dirDer(p, w, idim, 1));
246 else if (val > 0) return mu * dirDer(p, w, idim, 1);
247 else return -mu * dirDer(p, w, idim, 1);
248 }
Belle2::OrcaKinFit::BaseHardConstraint::getValue
virtual double getValue() const override=0
Returns the value of the constraint.
Belle2::OrcaKinFit::BaseHardConstraint::dirDer
virtual double dirDer(double *p, double *w, int idim, double mu=1)
Calculate directional derivative.
Definition: BaseHardConstraint.cc:232

◆ firstDerivatives()

bool firstDerivatives ( int  i,
double *  derivatives 
) const
overrideprotectedvirtual

First derivatives with respect to the 4-vector of Fit objects i; result false if all derivatives are zero.

Parameters
inumber of 1st FitObject
derivativesThe result 4-vector

Implements BaseHardConstraint.

Definition at line 208 of file MassConstraint.cc.

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 }

◆ getDerivatives()

void getDerivatives ( int  idim,
double  der[] 
) const
overridevirtual

Get first order derivatives.

Call this with a predefined array "der" with the necessary number of entries!

Parameters
idimFirst dimension of the array
derArray of derivatives, at least idim x idim

Implements BaseHardConstraint.

Definition at line 70 of file MassConstraint.cc.

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 }
E
R E
internal precision of FFTW codelets
Definition: DiscreteCosineTransform_31points.cc:36

◆ getError()

double getError ( ) const
overridevirtualinherited

Returns the error on the value of the constraint.

Reimplemented from BaseConstraint.

Definition at line 217 of file BaseHardConstraint.cc.

218 {
219 double dgdpi[BaseDefs::MAXINTERVARS];
220 double error2 = 0;
221 for (unsigned int i = 0; i < fitobjects.size(); ++i) {
222 const BaseFitObject* foi = fitobjects[i];
223 assert(foi);
224 if (firstDerivatives(i, dgdpi)) {
225 error2 += foi->getError2(dgdpi, getVarBasis());
226 }
227 }
228 return std::sqrt(std::abs(error2));
229 }

◆ getGlobalNum()

virtual int getGlobalNum ( ) const
inlinevirtualinherited

Accesses position of constraint in global constraint list.

Definition at line 137 of file BaseHardConstraint.h.

138 {return globalNum;}
Belle2::OrcaKinFit::BaseHardConstraint::globalNum
int globalNum
Position of constraint in global constraint list.
Definition: BaseHardConstraint.h:185

◆ getMass()

double getMass ( int  flag = 1)
virtual

Get the actual invariant mass of the fit objects with a given flag.

Parameters
flagThe flag

Definition at line 137 of file MassConstraint.cc.

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 }

◆ getName()

const char * getName ( ) const
virtualinherited

Returns the name of the constraint.

Definition at line 56 of file BaseConstraint.cc.

57 {
58 return name;
59 }

◆ getValue()

double getValue ( ) const
overridevirtual

Returns the value of the constraint.

Implements BaseHardConstraint.

Definition at line 44 of file MassConstraint.cc.

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 }

◆ getVarBasis()

int getVarBasis ( ) const
overridevirtual

Implements BaseHardConstraint.

Definition at line 241 of file MassConstraint.cc.

242 {
243 return VAR_BASIS;
244 }

◆ invalidateCache()

virtual void invalidateCache ( ) const
inlinevirtualinherited

Invalidates any cached values for the next event.

Reimplemented in MomentumConstraint.

Definition at line 104 of file ParticleConstraint.h.

105 {}

◆ num1stDerivative()

double num1stDerivative ( int  ifo,
int  ilocal,
double  eps 
)
virtualinherited

Evaluates numerically the 1st derivative w.r.t. a parameter.

Parameters
ifoNumber of FitObject
ilocalLocal parameter number
epsvariation of local parameter

Definition at line 307 of file BaseHardConstraint.cc.

308 {
309 BaseFitObject* fo = fitobjects[ifo];
310 assert(fo);
311 double save = fo->getParam(ilocal);
312 fo->setParam(ilocal, save + eps);
313 double v1 = getValue();
314 fo->setParam(ilocal, save - eps);
315 double v2 = getValue();
316 double result = (v1 - v2) / (2 * eps);
317 fo->setParam(ilocal, save);
318 return result;
319 }
VXDTFFilterTest::v2
const std::vector< double > v2
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:42
VXDTFFilterTest::v1
const std::vector< double > v1
MATLAB generated random vector.
Definition: decorrelationMatrix.cc:28

◆ num2ndDerivative()

double num2ndDerivative ( int  ifo1,
int  ilocal1,
double  eps1,
int  ifo2,
int  ilocal2,
double  eps2 
)
virtualinherited

Evaluates numerically the 2nd derivative w.r.t. 2 parameters.

Parameters
ifo1Number of 1st FitObject
ilocal11st local parameter number
eps1variation of 1st local parameter
ifo2Number of 1st FitObject
ilocal21st local parameter number
eps2variation of 2nd local parameter

Definition at line 321 of file BaseHardConstraint.cc.

323 {
324 double result;
325
326 if (ifo1 == ifo2 && ilocal1 == ilocal2) {
327 BaseFitObject* fo = fitobjects[ifo1];
328 assert(fo);
329 double save = fo->getParam(ilocal1);
330 double v0 = getValue();
331 fo->setParam(ilocal1, save + eps1);
332 double v1 = getValue();
333 fo->setParam(ilocal1, save - eps1);
334 double v2 = getValue();
335 result = (v1 + v2 - 2 * v0) / (eps1 * eps1);
336 fo->setParam(ilocal1, save);
337 } else {
338 BaseFitObject* fo1 = fitobjects[ifo1];
339 assert(fo1);
340 BaseFitObject* fo2 = fitobjects[ifo2];
341 assert(fo2);
342 double save1 = fo1->getParam(ilocal1);
343 double save2 = fo2->getParam(ilocal2);
344 fo1->setParam(ilocal1, save1 + eps1);
345 fo2->setParam(ilocal2, save2 + eps2);
346 double v11 = getValue();
347 fo2->setParam(ilocal2, save2 - eps2);
348 double v12 = getValue();
349 fo1->setParam(ilocal1, save1 - eps1);
350 double v22 = getValue();
351 fo2->setParam(ilocal2, save2 + eps2);
352 double v21 = getValue();
353 result = (v11 + v22 - v12 - v21) / (4 * eps1 * eps2);
354 fo1->setParam(ilocal1, save1);
355 fo2->setParam(ilocal2, save2);
356 }
357 return result;
358 }

◆ print()

std::ostream & print ( std::ostream &  os) const
virtualinherited

print object to ostream

Parameters
osThe output stream

Definition at line 76 of file BaseConstraint.cc.

77 {
78 os << getName() << "=" << getValue();
79 return os;
80 }
Belle2::OrcaKinFit::BaseConstraint::getValue
virtual double getValue() const =0
Returns the value of the constraint function.
Belle2::OrcaKinFit::BaseConstraint::getName
virtual const char * getName() const
Returns the name of the constraint.
Definition: BaseConstraint.cc:56

◆ printFirstDerivatives()

void printFirstDerivatives ( ) const
virtualinherited

Definition at line 360 of file BaseHardConstraint.cc.

361 {
362
363 B2INFO("BaseHardConstraint::printFirstDerivatives " << fitobjects.size());
364
365 double dgdpi[BaseDefs::MAXINTERVARS];
366 for (unsigned int i = 0; i < fitobjects.size(); ++i) {
367 const BaseFitObject* foi = fitobjects[i];
368 assert(foi);
369 if (firstDerivatives(i, dgdpi)) {
370 B2INFO("first derivs for obj " << i << " : ");
371 for (double j : dgdpi)
372 B2INFO(j << " ");
373 }
374 }
375
376 return;
377 }

◆ printSecondDerivatives()

void printSecondDerivatives ( ) const
virtualinherited

Definition at line 379 of file BaseHardConstraint.cc.

380 {
381
382 double d2GdPidPj[BaseDefs::MAXINTERVARS * BaseDefs::MAXINTERVARS];
383
384 int n = fitobjects.size();
385
386 for (int i = 0; i < n; ++i) {
387 const BaseFitObject* foi = fitobjects[i];
388 assert(foi);
389 for (int j = 0; j < n; ++j) {
390 const BaseFitObject* foj = fitobjects[j];
391 assert(foj);
392 if (secondDerivatives(i, j, d2GdPidPj)) {
393
394 B2INFO("second derivs for objs " << i << " " << j);
395
396 int k(0);
397 for (int k1 = 0; k1 < BaseDefs::MAXINTERVARS; k1++) {
398 for (int k2 = 0; k2 < BaseDefs::MAXINTERVARS; k2++) {
399 B2INFO(d2GdPidPj[k++] << " ");
400 }
401 }
402 }
403 }
404 }
405
406 return;
407 }

◆ resetFOList()

virtual void resetFOList ( )
inlinevirtualinherited

Resests ParticleFitObject list.

Definition at line 97 of file ParticleConstraint.h.

98 {
99 fitobjects.resize(0);
100 flags.resize(0);
101 };

◆ secondDerivatives()

bool secondDerivatives ( int  i,
int  j,
double *  derivatives 
) const
overrideprotectedvirtual

Second derivatives with respect to the 4-vectors of Fit objects i and j; result false if all derivatives are zero.

Parameters
inumber of 1st FitObject
jnumber of 2nd FitObject
derivativesThe result 4x4 matrix

Implements BaseHardConstraint.

Definition at line 162 of file MassConstraint.cc.

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 }

◆ setFOList()

virtual void setFOList ( std::vector< ParticleFitObject * > *  fitobjects_)
inlinevirtualinherited

Adds several ParticleFitObject objects to the list.

Parameters
fitobjects_A list of BaseFitObject objects

Definition at line 81 of file ParticleConstraint.h.

83 {
84 for (int i = 0; i < (int) fitobjects_->size(); i++) {
85 fitobjects.push_back(reinterpret_cast < BaseFitObject* >((*fitobjects_)[i]));
86 flags.push_back(1);
87 }
88 };

◆ setGlobalNum()

virtual void setGlobalNum ( int  iglobal)
inlinevirtualinherited

Sets position of constraint in global constraint list.

Parameters
iglobalGlobal constraint number

Definition at line 140 of file BaseHardConstraint.h.

142 {globalNum = iglobal;}

◆ setMass()

void setMass ( double  mass_)
virtual

Sets the target mass of the constraint.

Parameters
mass_The new mass

Definition at line 157 of file MassConstraint.cc.

158 {
159 mass = mass_;
160 }

◆ setName()

void setName ( const char *  name_)
inherited

Set object's name.

Definition at line 61 of file BaseConstraint.cc.

62 {
63 if (name_ == nullptr) return;
64 size_t l = strlen(name_);
65 if (name) delete[] name;
66 name = new char[l + 1];
67 strcpy(name, name_);
68 }

◆ test1stDerivatives()

void test1stDerivatives ( )
virtualinherited

Definition at line 251 of file BaseHardConstraint.cc.

252 {
253 B2INFO("BaseConstraint::test1stDerivatives for " << getName());
254 double y[100];
255 for (double& i : y) i = 0;
256 addToGlobalChi2DerVector(y, 100, 1);
257 double eps = 0.00001;
258 for (unsigned int ifo = 0; ifo < fitobjects.size(); ++ifo) {
259 BaseFitObject* fo = fitobjects[ifo];
260 assert(fo);
261 for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
262 int iglobal = fo->getGlobalParNum(ilocal);
263 double calc = y[iglobal];
264 double num = num1stDerivative(ifo, ilocal, eps);
265 B2INFO("fo: " << fo->getName() << " par " << ilocal << "/"
266 << iglobal << " (" << fo->getParamName(ilocal)
267 << ") calc: " << calc << " - num: " << num << " = " << calc - num);
268 }
269 }
270 }
Belle2::OrcaKinFit::BaseHardConstraint::num1stDerivative
virtual double num1stDerivative(int ifo, int ilocal, double eps)
Evaluates numerically the 1st derivative w.r.t. a parameter.
Definition: BaseHardConstraint.cc:307

◆ test2ndDerivatives()

void test2ndDerivatives ( )
virtualinherited

Definition at line 272 of file BaseHardConstraint.cc.

273 {
274 B2INFO("BaseConstraint::test2ndDerivatives for " << getName());
275 const int idim = 100;
276 auto* M = new double[idim * idim];
277 for (int i = 0; i < idim * idim; ++i) M[i] = 0;
278 add2ndDerivativesToMatrix(M, idim, 1);
279 double eps = 0.0001;
280 B2INFO("eps=" << eps);
281
282 for (unsigned int ifo1 = 0; ifo1 < fitobjects.size(); ++ifo1) {
283 BaseFitObject* fo1 = fitobjects[ifo1];
284 assert(fo1);
285 for (unsigned int ifo2 = ifo1; ifo2 < fitobjects.size(); ++ifo2) {
286 BaseFitObject* fo2 = fitobjects[ifo2];
287 assert(fo2);
288 for (int ilocal1 = 0; ilocal1 < fo1->getNPar(); ++ilocal1) {
289 int iglobal1 = fo1->getGlobalParNum(ilocal1);
290 for (int ilocal2 = (ifo1 == ifo2 ? ilocal1 : 0); ilocal2 < fo2->getNPar(); ++ilocal2) {
291 int iglobal2 = fo2->getGlobalParNum(ilocal2);
292 double calc = M[idim * iglobal1 + iglobal2];
293 double num = num2ndDerivative(ifo1, ilocal1, eps, ifo2, ilocal2, eps);
294 B2INFO("fo1: " << fo1->getName() << " par " << ilocal1 << "/"
295 << iglobal1 << " (" << fo1->getParamName(ilocal1)
296 << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/"
297 << iglobal2 << " (" << fo2->getParamName(ilocal2)
298 << ") calc: " << calc << " - num: " << num << " = " << calc - num);
299 }
300 }
301 }
302 }
303 delete[] M;
304 }
Belle2::OrcaKinFit::BaseHardConstraint::add2ndDerivativesToMatrix
virtual void add2ndDerivativesToMatrix(double *M, int idim, double lambda) const
Adds second order derivatives to global covariance matrix M.
Definition: BaseHardConstraint.cc:76
Belle2::OrcaKinFit::BaseHardConstraint::num2ndDerivative
virtual double num2ndDerivative(int ifo1, int ilocal1, double eps1, int ifo2, int ilocal2, double eps2)
Evaluates numerically the 2nd derivative w.r.t. 2 parameters.
Definition: BaseHardConstraint.cc:321

Friends And Related Function Documentation

◆ operator<<()

std::ostream & operator<< ( std::ostream &  os,
const BaseConstraint &  bc 
)
related

Prints out a BaseConstraint, using its print method.

Parameters
osThe output stream
bcThe object to print

Definition at line 114 of file BaseConstraint.h.

117 {
118 return bc.print(os);
119 }

Member Data Documentation

◆ derivatives

std::vector<double> derivatives
protectedinherited

The derivatives.

Definition at line 179 of file BaseHardConstraint.h.

◆ fitobjects

FitObjectContainer fitobjects
protectedinherited

The FitObjectContainer.

Definition at line 177 of file BaseHardConstraint.h.

◆ flags

std::vector<int> flags
protectedinherited

The flags can be used to divide the FitObjectContainer into several subsets used for example to implement an equal mass constraint (see MassConstraint).

Definition at line 182 of file BaseHardConstraint.h.

◆ globalNum

int globalNum
protectedinherited

Position of constraint in global constraint list.

Definition at line 185 of file BaseHardConstraint.h.

◆ mass

double mass
protected

The mass difference between object sets 1 and 2.

Definition at line 76 of file MassConstraint.h.

◆ name

char* name
protectedinherited

Definition at line 108 of file BaseConstraint.h.

The documentation for this class was generated from the following files: