Belle II Software light-2406-ragdoll
ISRPhotonFitObject.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#define NO_MARLIN // if defined: all output via B2INFO, Marlin inclusion not required
18#include "analysis/OrcaKinFit/ISRPhotonFitObject.h"
19#include <framework/logging/Logger.h>
20#include <cmath>
21
22#undef NDEBUG
23#include <cassert>
24
25#include <iostream>
26#ifndef NO_MARLIN
27#include "marlin/Processor.h"
28#endif
29
30using std::sqrt;
31using std::exp;
32using std::pow;
33using std::endl;
34#ifndef NO_MARLIN
35using namespace marlin;
36#endif
37
38namespace Belle2 {
43 namespace OrcaKinFit {
44
45 static const double pi_ = M_PI, // 3.14159265358979323846264338328,
46 a = 8. / 3. / pi_ * (pi_ - 3.) / (4. - pi_); // = ca. 0.140012289
47
48// constructor
49 ISRPhotonFitObject::ISRPhotonFitObject(double px, double py, double ppz,
50 double b_, double PzMaxB_, double PzMinB_)
51 : cachevalid(false),
52 pt2(0), p2(0), p(0), pz(0),
53 dpx0(0), dpy0(0), dpz0(0), dE0(0), dpx1(0), dpy1(0), dpz1(0), dE1(0),
54 dpx2(0), dpy2(0), dpz2(0), dE2(0), d2pz22(0), d2E22(0),
55 chi2(0), b(0), PzMinB(0), PzMaxB(0), dp2zFact(0)
56 {
57
58 assert(int(NPAR) <= int(BaseDefs::MAXPAR));
59
60 initCov();
61 b = b_;
62 PzMinB = PzMinB_;
63 PzMaxB = PzMaxB_;
64#ifdef DEBUG
65 B2INFO("ISRPhotonFitObject: b: " << b << " PzMinB: " << PzMinB << " PzMaxB: " << PzMaxB);
66#endif
67
68 if (b <= 0. || b >= 1.) {
69 B2INFO("ISRPhotonFitObject: b must be from ]0,1[ ");
70 }
71 assert(b > 0. && b < 1.);
72 if (PzMinB < 0. || PzMaxB <= PzMinB) {
73 B2INFO("ISRPhotonFitObject: PzMinB and PzMaxB must be chosen such that 0 <= PzMinB < PzMaxB");
74 }
75 assert(PzMinB >= 0.);
76 assert(PzMaxB > PzMinB);
77 dp2zFact = (PzMaxB - PzMinB) / b * sqrt(2. / pi_);
78 double pg = PgFromPz(ppz); // using internally Gauss-distributed parameter p_g instead of p_z
79 setParam(0, px, true, true);
80 setParam(1, py, true, true);
81 setParam(2, pg, true);
82 setMParam(0, 0.); // all measured parameters
83 setMParam(1, 0.); // are assumed to be zero
84 setMParam(2, 0.); // in this photon parametrization
85#ifdef DEBUG
86 B2INFO("ISRPhotonFitObject: Initial pg: " << pg);
87#endif
88 setError(2, 1.);
89 setMass(0.);
91 }
92
93// destructor
94 ISRPhotonFitObject::~ISRPhotonFitObject() = default;
95
96
98 : ParticleFitObject(rhs), cachevalid(false),
99 pt2(0), p2(0), p(0), pz(0),
100 dpx0(0), dpy0(0), dpz0(0), dE0(0), dpx1(0), dpy1(0), dpz1(0), dE1(0),
101 dpx2(0), dpy2(0), dpz2(0), dE2(0), d2pz22(0), d2E22(0),
102 chi2(0), b(0), PzMinB(0), PzMaxB(0), dp2zFact(0)
103 {
105 }
106
108 {
109 if (this != &rhs) {
110 assign(rhs); // calls virtual function assign of derived class
111 }
112 return *this;
113 }
114
116 {
117 return new ISRPhotonFitObject(*this);
118 }
119
121 {
122 if (const auto* psource = dynamic_cast<const ISRPhotonFitObject*>(&source)) {
123 if (psource != this) {
125 // only mutable data members, need not to be copied, if cache is invalid
126 }
127 } else {
128 assert(0);
129 }
130 return *this;
131 }
132
133 const char* ISRPhotonFitObject::getParamName(int ilocal) const
134 {
135 switch (ilocal) {
136 case 0: return "P_x";
137 case 1: return "P_y";
138 case 2: return "P_g";
139 }
140 return "undefined";
141 }
142
143 bool ISRPhotonFitObject::updateParams(double pp[], int idim)
144 {
146 int i2 = getGlobalParNum(2);
147 assert(i2 >= 0 && i2 < idim);
148 double pp2 = pp[i2];
149#ifdef DEBUG
150 B2INFO("ISRPhotonFitObject::updateParams: p2(new) = " << pp[i2] << " par[2](old) = " << par[2]);
151#endif
152 bool result = ((pp2 - par[2]) * (pp2 - par[2]) > eps2 * cov[2][2]);
153 par[2] = pp2;
154 pp[i2] = par[2];
155 return result;
156 }
157
158 double ISRPhotonFitObject::getDPx(int ilocal) const
159 {
160 assert(ilocal >= 0 && ilocal < NPAR);
161 if (!cachevalid) updateCache();
162 switch (ilocal) {
163 case 0: return dpx0;
164 case 1: return dpx1;
165 case 2: return dpx2;
166 }
167 return 0;
168 }
169
170 double ISRPhotonFitObject::getDPy(int ilocal) const
171 {
172 assert(ilocal >= 0 && ilocal < NPAR);
173 if (!cachevalid) updateCache();
174 switch (ilocal) {
175 case 0: return dpy0;
176 case 1: return dpy1;
177 case 2: return dpy2;
178 }
179 return 0;
180 }
181
182 double ISRPhotonFitObject::getDPz(int ilocal) const
183 {
184 assert(ilocal >= 0 && ilocal < NPAR);
185 if (!cachevalid) updateCache();
186 switch (ilocal) {
187 case 0: return dpz0;
188 case 1: return dpz1;
189 case 2: return dpz2;
190 }
191 return 0;
192 }
193
194 double ISRPhotonFitObject::getDE(int ilocal) const
195 {
196 assert(ilocal >= 0 && ilocal < NPAR);
197 if (!cachevalid) updateCache();
198 switch (ilocal) {
199 case 0: return dE0;
200 case 1: return dE1;
201 case 2: return dE2;
202 }
203 return 0;
204 }
205
206
207 double ISRPhotonFitObject::getFirstDerivative_Meta_Local(int iMeta, int ilocal, int metaSet) const
208 {
209 assert(metaSet == 0);
210 switch (iMeta) {
211 case 0:
212 return getDE(ilocal);
213 break;
214 case 1:
215 return getDPx(ilocal);
216 break;
217 case 2:
218 return getDPy(ilocal);
219 break;
220 case 3:
221 return getDPz(ilocal);
222 break;
223 default:
224 assert(0);
225 }
226 return -999;
227 }
228
229 double ISRPhotonFitObject::getSecondDerivative_Meta_Local(int iMeta, int ilocal, int jlocal, int metaSet) const
230 {
231 assert(metaSet == 0);
232 if (!cachevalid) updateCache();
233
234 if (jlocal < ilocal) {
235 int temp = jlocal;
236 jlocal = ilocal;
237 ilocal = temp;
238 }
239
240 // daniel hasn't checked these, copied from orig code
241 switch (iMeta) {
242
243 case 0:
244 if (ilocal == 2 && jlocal == 2) return d2E22;
245 else return 0;
246 break;
247 case 1:
248 case 2:
249 return 0;
250 break;
251 case 3:
252 if (ilocal == 2 && jlocal == 2) return d2pz22;
253 else return 0;
254 break;
255 default:
256 return 0;
257 }
258
259 }
260
261
262
263 double ISRPhotonFitObject::PgFromPz(double ppz)
264 {
265
266 int sign = (ppz > 0.) - (ppz < 0.);
267 double u = (pow(fabs(ppz), b) - PzMinB) / (PzMaxB - PzMinB);
268
269 if (u < 0.) {
270 B2INFO("ISRPhotonFitObject: Initial pz with abs(pz) < pzMin adjusted to zero.");
271 u = 0.;
272 }
273
274 if (u >= 1.) {
275 B2INFO("ISRPhotonFitObject: Initial pz with abs(pz) >= pzMax adjusted.");
276 u = 0.99999999;
277 }
278
279 double g = std::log(1. - u * u);
280 double g4pa = g + 4. / pi_ / a;
281 return sign * sqrt(-g4pa + sqrt(g4pa * g4pa - 4. / a * g)) ;
282 }
283
284
285 void ISRPhotonFitObject::updateCache() const
286 {
287 double px = par[0];
288 double py = par[1];
289 double pg = par[2];
290
291 int sign = (pg > 0.) - (pg < 0.);
292 double pg2h = pg * pg / 2.;
293 double exponent = -pg2h * (4. / pi_ + a * pg2h) / (1. + a * pg2h);
294 double u = sqrt((exponent < -1.e-14) ? 1. - exp(exponent) : -exponent); // approximation to avoid numerical problem
295 pz = sign * pow((PzMinB + (PzMaxB - PzMinB) * u), (1. / b));
296
297 pt2 = px * px + py * py;
298 p2 = pt2 + pz * pz;
299 p = std::sqrt(p2);
300
301 fourMomentum.setValues(p, px, py, pz);
302
303 dpx0 = 1.;
304 dpx1 = 0.;
305 dpx2 = 0.;
306 dpy0 = 0.;
307 dpy1 = 1.;
308 dpy2 = 0.;
309 dpz0 = 0.;
310 dpz1 = 0.;
311 dpz2 = dp2zFact * pow(fabs(pz), (1. - b)) * exp(-pg * pg / 2.);
312
313 // if p,pz==0, derivatives are zero (catch up 1/0)
314 if (pz) {
315 d2pz22 = dpz2 * ((1. - b) * dpz2 / pz - par[2]);
316 } else {
317 d2pz22 = 0.;
318 }
319 if (p) {
320 dE0 = px / p;
321 dE1 = py / p;
322 dE2 = pz / p * dpz2;
323 d2E22 = pz / p * d2pz22;
324 if (pt2) {
325 d2E22 += pt2 / p / p / p * dpz2 * dpz2; // NOT using /p2/p to avoid numerical problems
326 }
327 } else {
328 dE0 = 0.;
329 dE1 = 0.;
330 dE2 = 0.;
331 d2E22 = 0.;
332 }
333
334#ifdef DEBUG
335 B2INFO("ISRPhotonFitObject::updateCache: pg: " << pg << " pz: " << pz << " p: " << p << " p^2: " << p2 << "\n"
336 << " Dpz/Dpg: " << dpz2 << " DE/Dpg: " << dE2 << " D^2pz/Dpg^2: " << d2pz22
337 << " D^2E/Dpg^2: " << d2E22);
338#endif
339
340 cachevalid = true;
341 }
342
343 }// end OrcaKinFit namespace
345} // end Belle2 namespace
346
virtual bool setParam(int ilocal, double par_, bool measured_, bool fixed_=false)
Set value and measured flag of parameter i; return: significant change.
virtual int getGlobalParNum(int ilocal) const
Get global parameter number of parameter ilocal.
virtual bool setMParam(int ilocal, double mpar_)
Set measured value of parameter ilocal; return: success.
double par[BaseDefs::MAXPAR]
fit parameters
void invalidateCache() const
invalidate any cached quantities
double cov[BaseDefs::MAXPAR][BaseDefs::MAXPAR]
local covariance matrix
virtual bool setError(int ilocal, double err_)
Set error of parameter ilocal; return: success.
virtual double getDPx(int ilocal) const override
Return d p_x / d par_ilocal (derivative of px w.r.t. local parameter ilocal)
virtual double getDPy(int ilocal) const override
Return d p_y / d par_ilocal (derivative of py w.r.t. local parameter ilocal)
virtual double getDE(int ilocal) const override
Return d E / d par_ilocal (derivative of E w.r.t. local parameter ilocal)
ISRPhotonFitObject & operator=(const ISRPhotonFitObject &rhs)
Assignment.
virtual ISRPhotonFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
virtual const char * getParamName(int ilocal) const override
Get name of parameter ilocal.
virtual double getDPz(int ilocal) const override
Return d p_z / d par_ilocal (derivative of pz w.r.t. local parameter ilocal)
virtual bool updateParams(double p[], int idim) override
Read values from global vector, readjust vector; return: significant change.
ISRPhotonFitObject(double px, double py, double pz, double b_, double PzMaxB_, double PzMinB_=0.)
virtual ISRPhotonFitObject * copy() const override
Return a new copy of itself.
virtual bool setMass(double mass_)
Set mass of particle; return=success.
virtual ParticleFitObject & assign(const BaseFitObject &source) override
Assign from anther object, if of same type.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24