Belle II Software  release-08-01-10
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 
39 using namespace std;
40 
41 namespace 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