Belle II Software development
ParameterScanner.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/ParameterScanner.h"
20
21#include "analysis/OrcaKinFit/BaseFitter.h"
22#include "analysis/OrcaKinFit/NewFitterGSL.h"
23#include "analysis/OrcaKinFit/BaseFitObject.h"
24#include "analysis/OrcaKinFit/BaseHardConstraint.h"
25#include <framework/logging/Logger.h>
26
27#include <TString.h>
28#include <TH2F.h>
29#include <TMultiGraph.h>
30#include <TGraph.h>
31
32#include <gsl/gsl_vector.h>
33
34#undef NDEBUG
35#include <cassert>
36#include <cmath>
37#include <iostream>
38
39using namespace std;
40
41namespace Belle2 {
46 namespace OrcaKinFit {
47
48 ParameterScanner::ParameterScanner(BaseFitter& fitter_)
49 : fitter(fitter_)
50 {}
51
52 void ParameterScanner::doScan(int xglobal,
53 int nx,
54 double xstart,
55 double xstop,
56 int yglobal,
57 int ny,
58 double ystart,
59 double ystop,
60 const char* idprefix,
61 const char* titleprefix,
62 double mumerit)
63 {
64
65 NewFitterGSL* newfitter = dynamic_cast<NewFitterGSL*>(&fitter);
66
67 FitObjectContainer* fitobjects = fitter.getFitObjects();
68 if (fitobjects == 0) return;
69 ConstraintContainer* constraints = fitter.getConstraints();
70 if (constraints == 0) return;
71
72 // Get largest global parameter number,
73 // find parameter names
74 TString xname("");
75 TString yname("");
76
77 int idim = 1;
78 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
79 BaseFitObject* fo = *i;
80 assert(fo);
81 for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
82 int iglobal = fo->getGlobalParNum(ilocal);
83 if (iglobal >= idim) idim = iglobal + 1;
84 if (iglobal == xglobal) xname = fo->getParamName(ilocal);
85 if (iglobal == yglobal) yname = fo->getParamName(ilocal);
86 }
87 }
88 for (ConstraintIterator i = constraints->begin(); i != constraints->end(); ++i) {
89 BaseHardConstraint* c = *i;
90 assert(c);
91 int iglobal = c->getGlobalNum();
92 if (iglobal >= idim) idim = iglobal + 1;
93 }
94
95 if (xglobal >= idim) return;
96 if (yglobal >= idim) return;
97
98 double* parsave = new double [idim];
99 double* par = new double [idim];
100
101 // Get starting values
102 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
103 BaseFitObject* fo = *i;
104 assert(fo);
105 for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
106 int iglobal = fo->getGlobalParNum(ilocal);
107 assert(iglobal >= 0 && iglobal < idim);
108 parsave[iglobal] = fo->getParam(ilocal);
109 }
110 }
111
112 // Book Histograms
113
114 TString idpostfix("");
115 idpostfix += xglobal;
116 idpostfix += "_";
117 idpostfix += yglobal;
118
119 TString titlepostfix(" vs ");
120 titlepostfix += xname;
121 titlepostfix += " and ";
122 titlepostfix += yname;
123 titlepostfix += ";";
124 titlepostfix += xname;
125 titlepostfix += ";";
126 titlepostfix += yname;
127
128 // Chi2 Histogram
129
130 TString id(idprefix);
131 id += "chi2_";
132 id += idpostfix;
133
134 TString title(titleprefix);
135 title += "Chi2 ";
136 title += titlepostfix;
137
138 TH2F* hchi2 = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
139 B2INFO("Booking Histo '" << id << "': '" << title << "'");
140
141 TH2F* halpha = 0;
142 TH2F* hlog2alpha = 0;
143 TH2F* hmu = 0;
144 TH2F* hphi1 = 0;
145
146 if (newfitter) {
147 id = idprefix;
148 id += "alpha_";
149 id += idpostfix;
150
151 title = titleprefix;
152 title += "#alpha ";
153 title += titlepostfix;
154
155 halpha = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
156 B2INFO("Booking Histo '" << id << "': '" << title << "'");
157
158 id = idprefix;
159 id += "log2alpha_";
160 id += idpostfix;
161
162 title = titleprefix;
163 title += "log_{2} (#alpha) ";
164 title += titlepostfix;
165
166 hlog2alpha = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
167 B2INFO("Booking Histo '" << id << "': '" << title << "'");
168
169 id = idprefix;
170 id += "mu_";
171 id += idpostfix;
172
173 title = titleprefix;
174 title += "#mu ";
175 title += titlepostfix;
176
177 hmu = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
178 B2INFO("Booking Histo '" << id << "': '" << title << "'");
179
180 id = idprefix;
181 id += "phi1_";
182 id += idpostfix;
183
184 title = titleprefix;
185 title += "Merit function #phi_{1} ";
186 title += titlepostfix;
187
188 hphi1 = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
189 B2INFO("Booking Histo '" << id << "': '" << title << "'");
190
191
192 }
193
194 TMultiGraph* mgstepsfull = 0;
195 TMultiGraph* mgsteps = 0;
196 TGraph* gstep0 = 0;
197 TGraph* gstep1 = 0;
198 TGraph* gstep2 = 0;
199 if (newfitter) {
200 id = idprefix;
201 id += "stepsfull_";
202 id += idpostfix;
203
204 title = titleprefix;
205 title += "Full Steps ";
206 title += titlepostfix;
207
208 mgstepsfull = new TMultiGraph(id, title);
209 B2INFO("Booking Multigraph '" << id << "': '" << title << "'");
210
211 id = idprefix;
212 id += "steps_";
213 id += idpostfix;
214
215 title = titleprefix;
216 title += "Steps ";
217 title += titlepostfix;
218
219 mgsteps = new TMultiGraph(id, title);
220 B2INFO("Booking Multigraph '" << id << "': '" << title << "'");
221
222 id = idprefix;
223 id += "stepstart_";
224 id += idpostfix;
225
226 title = titleprefix;
227 title += "Step Start Points ";
228 title += titlepostfix;
229
230 gstep0 = new TGraph(nx * ny);
231 gstep0 ->SetName(id);
232 gstep0 ->SetTitle(title);
233 gstep0->SetMarkerStyle(20);
234 gstep0->SetMarkerColor(kBlack);
235 gstep0->SetMarkerSize(0.5);
236
237 id = idprefix;
238 id += "stependfull_";
239 id += idpostfix;
240
241 title = titleprefix;
242 title += "Full Step End Points ";
243 title += titlepostfix;
244
245 gstep1 = new TGraph(nx * ny);
246 gstep1 ->SetName(id);
247 gstep1 ->SetTitle(title);
248 gstep1->SetMarkerStyle(20);
249 gstep1->SetMarkerColor(kRed);
250 gstep1->SetMarkerSize(0.5);
251
252 id = idprefix;
253 id += "stepend_";
254 id += idpostfix;
255
256 title = titleprefix;
257 title += "Step End Points ";
258 title += titlepostfix;
259
260 gstep2 = new TGraph(nx * ny);
261 gstep2 ->SetName(id);
262 gstep2 ->SetTitle(title);
263 gstep2->SetMarkerStyle(21);
264 gstep2->SetMarkerColor(kGreen);
265 gstep2->SetMarkerSize(0.5);
266
267
268 }
269
270 unsigned int ncon = constraints->size();
271 if (ncon > NCONMAX) ncon = NCONMAX;
272 TH2F* hcon[NCONMAX];
273 TH2F* hlambda[NCONMAX];
274 for (unsigned int icon = 0; icon < ncon; ++icon) {
275 hcon[icon] = 0;
276 hlambda[icon] = 0;
277 BaseConstraint* c = (*constraints)[icon];
278 if (c) {
279 id = idprefix;
280 id += "con";
281 id += icon;
282 id += "_";
283 id += idpostfix;
284
285 title = titleprefix;
286 title += "Constraint ";
287 title += icon;
288 title += ": ";
289 title += c->getName();
290 title += titlepostfix;
291
292 hcon[icon] = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
293 B2INFO("Booking Histo '" << id << "': '" << title << "'");
294
295 if (newfitter) {
296 id = idprefix;
297 id += "lambda";
298 id += icon;
299 id += "_";
300 id += idpostfix;
301
302 title = titleprefix;
303 title += "Lambda ";
304 title += icon;
305 title += ": ";
306 title += c->getName();
307 title += titlepostfix;
308
309 hlambda[icon] = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
310 B2INFO("Booking Histo '" << id << "': '" << title << "'");
311
312 }
313
314
315 }
316 }
317
318
319 // Do the scan
320
321 for (int ix = 1; ix <= nx; ++ix) {
322 double x = (ix - 0.5) * (xstop - xstart) / nx + xstart;
323 for (int iy = 1; iy <= ny; ++iy) {
324 double y = (iy - 0.5) * (ystop - ystart) / ny + ystart;
325
326 // Set parameters
327 for (int i = 0; i < idim; ++i) par[i] = parsave[i];
328 par[xglobal] = x;
329 par[yglobal] = y;
330 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
331 BaseFitObject* fo = *i;
332 assert(fo);
333 fo->updateParams(par, idim);
334 }
335
336
337 // Calculate chi2
338 double chi2 = 0;
339 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
340 BaseFitObject* fo = *i;
341 assert(fo);
342 chi2 += fo->getChi2();
343 }
344
345 //B2INFO( "Chi2 for x=" << x << " and y=" << y << ": " << chi2);
346
347 if (newfitter) {
348 newfitter->fillx(newfitter->x);
349 newfitter->assembleConstDer(newfitter->M);
350 newfitter->determineLambdas(newfitter->x, newfitter->M, newfitter->x, newfitter->W, newfitter->v1);
351 }
352
353 // Fill Histos
354 hchi2->SetBinContent(ix, iy, chi2);
355 for (unsigned int icon = 0; icon < ncon; ++icon) {
356 BaseHardConstraint* c = (*constraints)[icon];
357 TH2F* h = hcon[icon];
358 if (c && h) h->SetBinContent(ix, iy, c->getValue());
359 if (c) B2INFO("x=" << x << ", y=" << y << ": Constraint " << c->getName() << ": " << c->getValue())
360 h = hlambda[icon];
361 if (c && h && newfitter) {
362 int kglobal = c->getGlobalNum();
363 h->SetBinContent(ix, iy, gsl_vector_get(newfitter->x, kglobal));;
364 }
365
366 }
367
368 if (newfitter) {
369 newfitter->setDebug((ix == 1) ? 4 : 0);
370
371 double xval[2], yval[2];
372 xval[0] = x;
373 yval[0] = y;
374 gstep0->SetPoint((ix - 1)*ny + (iy - 1), xval[0], yval[0]);
375 newfitter->fillx(newfitter->x);
376 newfitter->fillperr(newfitter->perr);
377
378 newfitter->assembleConstDer(newfitter->M);
379 newfitter->determineLambdas(newfitter->x, newfitter->M, newfitter->x, newfitter->W, newfitter->v1);
380
381 double phi1 = newfitter->meritFunction(mumerit, newfitter->x, newfitter->perr);
382 if (hphi1) hphi1->SetBinContent(ix, iy, phi1);
383
384 newfitter->calcNewtonDx(newfitter->dx, newfitter->dxscal, newfitter->x,
385 newfitter->perr, newfitter->M, newfitter->Mscal,
386 newfitter->y, newfitter->yscal, newfitter->W, newfitter->W2,
387 newfitter->permW, newfitter->v1);
388 newfitter->add(newfitter->xnew, newfitter->x, 1, newfitter->dx);
389
390 xval[1] = gsl_vector_get(newfitter->xnew, xglobal);
391 yval[1] = gsl_vector_get(newfitter->xnew, yglobal);
392
393 B2INFO("ParameterScanner::doScan: full step from (" << xval[0] << ", " << yval[0]
394 << ") -> (" << xval[1] << ", " << yval[1] << ")");
395
396
397 gstep1->SetPoint((ix - 1)*ny + (iy - 1), xval[1], yval[1]);
398
399 TGraph* g = new TGraph(2, xval, yval);
400 g->SetLineColor(kRed);
401
402 if (xval[1] > xstart - 10 * (xstop - xstart) &&
403 xval[1] < xstop + 10 * (xstop - xstart) &&
404 yval[1] > ystart - 10 * (ystop - ystart) &&
405 yval[1] < ystop + 10 * (ystop - ystart)) {
406
407 mgstepsfull->Add(g, "L");
408 B2INFO(" -> added\n");
409 }
410 double alpha = 1;
411 double mu = 0;
412 int imode = 2;
413
414 newfitter->calcLimitedDx(alpha, mu, newfitter->xnew, imode,
415 newfitter->x, newfitter->v2, newfitter->dx, newfitter->dxscal,
416 newfitter->perr, newfitter->M, newfitter->Mscal,
417 newfitter->W, newfitter->v1);
418
419 xval[1] = gsl_vector_get(newfitter->xnew, xglobal);
420 yval[1] = gsl_vector_get(newfitter->xnew, yglobal);
421
422 B2INFO(" limited step from (" << xval[0] << ", " << yval[0]
423 << ") -> (" << xval[1] << ", " << yval[1] << "), alpha=" << alpha);
424
425 gstep2->SetPoint((ix - 1)*ny + (iy - 1), xval[1], yval[1]);
426
427 g = new TGraph(2, xval, yval);
428 g->SetLineColor(kGreen);
429
430 if (xval[1] > xstart - 10 * (xstop - xstart) &&
431 xval[1] < xstop + 10 * (xstop - xstart) &&
432 yval[1] > ystart - 10 * (ystop - ystart) &&
433 yval[1] < ystop + 10 * (ystop - ystart)) {
434
435 mgsteps->Add(g, "L");
436 B2INFO(" -> added\n");
437 }
438 if (halpha) halpha->SetBinContent(ix, iy, alpha);
439 if (hlog2alpha) hlog2alpha->SetBinContent(ix, iy, std::log(alpha) / std::log(2.));
440 if (hmu) hmu->SetBinContent(ix, iy, mu);
441
442
443 }
444
445 }
446 }
447
448 // Write histos;
449 hchi2->Write();
450 for (unsigned int icon = 0; icon < ncon; ++icon) {
451 if (hcon[icon]) hcon[icon]->Write();
452 if (hlambda[icon]) hcon[icon]->Write();
453 }
454
455 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
456 BaseFitObject* fo = *i;
457 assert(fo);
458 fo->updateParams(parsave, idim);
459 }
460
461 if (mgstepsfull) mgstepsfull->Write();
462 if (mgsteps) mgsteps->Write();
463 if (gstep0) gstep0->Write();
464 if (gstep1) gstep1->Write();
465 if (gstep2) gstep2->Write();
466 if (halpha) halpha->Write();
467 if (hlog2alpha) hlog2alpha->Write();
468 if (hmu) hmu->Write();
469 if (hphi1) hphi1->Write();
470
471 delete[] par;
472 delete[] parsave;
473 }
474
475 }// end OrcaKinFit namespace
477} // end Belle2 namespace
478
479#endif // MARLIN_USE_ROOT
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17
std
STL namespace.