Belle II Software light-2406-ragdoll
IterationScanner.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/IterationScanner.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
30#include <gsl/gsl_vector.h>
31
32#undef NDEBUG
33#include <cassert>
34#include <iostream>
35
36using namespace std;
37
38namespace Belle2 {
43 namespace OrcaKinFit {
44
45 IterationScanner::IterationScanner(BaseFitter& fitter_)
46 : fitter(fitter_)
47 {}
48
49 void IterationScanner::doScan(int xglobal,
50 int nx,
51 double xstart,
52 double xstop,
53 int yglobal,
54 int ny,
55 double ystart,
56 double ystop,
57 const char* idprefix,
58 const char* titleprefix)
59 {
60
61 NewFitterGSL* newfitter = dynamic_cast<NewFitterGSL*>(&fitter);
62
63 FitObjectContainer* fitobjects = fitter.getFitObjects();
64 if (fitobjects == 0) return;
65 ConstraintContainer* constraints = fitter.getConstraints();
66 if (constraints == 0) return;
67
68 FitObjectContainer fitobjects_backup(fitobjects->size());
69 for (unsigned int i = 0; i < fitobjects->size(); ++i) {
70 BaseFitObject* fo = (*fitobjects)[i];
71 assert(fo);
72 fitobjects_backup[i] = fo->copy();
73 }
74 // Get largest global parameter number,
75 // find parameter names
76 TString xname("");
77 TString yname("");
78
79 int idim = 1;
80 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
81 BaseFitObject* fo = *i;
82 assert(fo);
83 for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
84 int iglobal = fo->getGlobalParNum(ilocal);
85 if (iglobal >= idim) idim = iglobal + 1;
86 if (iglobal == xglobal) xname = fo->getParamName(ilocal);
87 if (iglobal == yglobal) yname = fo->getParamName(ilocal);
88 }
89 }
90 for (ConstraintIterator i = constraints->begin(); i != constraints->end(); ++i) {
91 BaseHardConstraint* c = *i;
92 assert(c);
93 int iglobal = c->getGlobalNum();
94 if (iglobal >= idim) idim = iglobal + 1;
95 }
96
97 if (xglobal >= idim) return;
98 if (yglobal >= idim) return;
99
100 double* parsave = new double [idim];
101 double* par = new double [idim];
102
103 // Get starting values
104 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
105 BaseFitObject* fo = *i;
106 assert(fo);
107 for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
108 int iglobal = fo->getGlobalParNum(ilocal);
109 assert(iglobal >= 0 && iglobal < idim);
110 parsave[iglobal] = fo->getParam(ilocal);
111 }
112 }
113
114 // Book Histograms
115
116 TString idpostfix("");
117 idpostfix += xglobal;
118 idpostfix += "_";
119 idpostfix += yglobal;
120
121 TString titlepostfix(" vs ");
122 titlepostfix += xname;
123 titlepostfix += " and ";
124 titlepostfix += yname;
125 titlepostfix += ";";
126 titlepostfix += xname;
127 titlepostfix += ";";
128 titlepostfix += yname;
129
130 // Chi2 Histogram
131
132 TString id(idprefix);
133 id += "nit_";
134 id += idpostfix;
135
136 TString title(titleprefix);
137 title += "Iterations ";
138 title += titlepostfix;
139
140 TH2F* hnit = new TH2F(id, title, nx, xstart, xstop, ny, ystart, ystop);
141 B2INFO("Booking Histo '" << id << "': '" << title << "'");
142
143
144 // Do the scan
145
146 for (int ix = 1; ix <= nx; ++ix) {
147 double x = (ix - 0.5) * (xstop - xstart) / nx + xstart;
148 for (int iy = 1; iy <= ny; ++iy) {
149 double y = (iy - 0.5) * (ystop - ystart) / ny + ystart;
150
151 // Set parameters
152 for (int i = 0; i < idim; ++i) par[i] = parsave[i];
153 par[xglobal] = x;
154 par[yglobal] = y;
155 for (FitObjectIterator i = fitobjects->begin(); i != fitobjects->end(); ++i) {
156 BaseFitObject* fo = *i;
157 assert(fo);
158 fo->updateParams(par, idim);
159 }
160
161 for (unsigned int i = 0; i < fitobjects->size(); ++i) {
162 BaseFitObject* fo = (*fitobjects)[i];
163 BaseFitObject* fobu = fitobjects_backup[i];
164 assert(fo);
165 assert(fobu);
166 for (int j = 0; j < fobu->getNPar(); ++j) {
167 for (int k = 0; k < fobu->getNPar(); ++k) {
168 fo->setCov(j, k, fobu->getCov(j, k));
169 }
170 }
171 }
172
173 // Calculate nit
174
175 double fprob = fitter.fit();
176 double nit = fitter.getIterations();
177
178 B2INFO("IterationScanner::doScan: x=" << x);
179 if (newfitter) B2INFO(" -> " << gsl_vector_get(newfitter->x, xglobal));
180 B2INFO(", y=" << y);
181 if (newfitter) B2INFO(" -> " << gsl_vector_get(newfitter->x, yglobal));
182 B2INFO(", fitprob=" << fprob << ", nit=" << nit);
183
184 // Fill Histos
185 hnit->SetBinContent(ix, iy, nit);
186 }
187 }
188
189 // Write histos;
190 hnit->Write();
191
192 for (unsigned int i = 0; i < fitobjects_backup.size(); ++i) {
193 delete fitobjects_backup[i];
194 }
195
196 delete[] par;
197 delete[] parsave;
198 }
199
200 }// end OrcaKinFit namespace
202} // end Belle2 namespace
203
204#endif // MARLIN_USE_ROOT
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.