Bug Summary

File:include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h
Warning:line 185, column 39
The result of left shift is undefined because the right operand '32' is not smaller than 32, the capacity of 'int'

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -O3 -analyze -disable-free -clear-ast-before-backend -disable-llvm-verifier -discard-value-names -main-file-name BeamSpotAlgorithm.cc -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -ffp-contract=on -fno-rounding-math -mconstructor-aliases -funwind-tables=2 -target-cpu x86-64 -tune-cpu generic -debugger-tuning=gdb -fdebug-compilation-dir=/data/b2soft/buildbot/development/build -fcoverage-compilation-dir=/data/b2soft/buildbot/development/build -resource-dir /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/lib/clang/21 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/c++ -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/c++/x86_64-redhat-linux -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/c++/backward -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/python3.12 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/include/CLHEP -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/Geant4 -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/include/root -isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/include/belle_legacy -I include/ -D GBL_EIGEN_SUPPORT_ROOT -D _PACKAGE_="reconstruction" -D G4UI_USE_TCSH -D RaveDllExport= -D HAS_SQLITE -D HAS_CALLGRIND -I include -I /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/include/libxml2 -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++ -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/x86_64-redhat-linux -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../include/c++/backward -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/lib/clang/21/include -internal-isystem /usr/local/include -internal-isystem /cvmfs/belle.cern.ch/el9/externals/v02-04-01/Linux_x86_64/common/bin/../lib64/gcc/x86_64-redhat-linux/15.2.0/../../../../x86_64-redhat-linux/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -Wno-missing-braces -Wno-unused-command-line-argument -std=c++20 -fdeprecated-macro -ferror-limit 19 -fgnuc-version=4.2.1 -fno-implicit-modules -fskip-odr-check-in-gmf -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -D__GCC_HAVE_DWARF2_CFI_ASM=1 -o /scan_build/2026-07-11-002547-1017033-1 -x c++ reconstruction/calibration/BeamSpotBoostInvMass/src/BeamSpotAlgorithm.cc

reconstruction/calibration/BeamSpotBoostInvMass/src/BeamSpotAlgorithm.cc

→
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <mdst/dbobjects/BeamSpot.h>
10#include <reconstruction/calibration/BeamSpotBoostInvMass/BeamSpotAlgorithm.h>
11#include <reconstruction/calibration/BeamSpotBoostInvMass/BeamSpotStandAlone.h>
12#include <reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h>
13
14#include <Eigen/Dense>
15
16using Eigen::Vector3d;
17using Eigen::Matrix3d;
18
19using namespace Belle2;
20
21using Belle2::BeamSpotCalib::getEvents;
22using Belle2::BeamSpotCalib::runBeamSpotAnalysis;
23
24BeamSpotAlgorithm::BeamSpotAlgorithm() : CalibrationAlgorithm("BeamSpotCollector")
25{
26 setDescription("BeamSpot calibration algorithm");
27}
28
29
30/** Create BS object */
31static TObject* getBeamSpotObj(Vector3d ipVtx, Matrix3d ipVtxUnc, Matrix3d sizeMat)
32{
33 auto payload = new BeamSpot();
34 payload->setIP(toB2Vector3(ipVtx), toTMatrixDSym(ipVtxUnc));
35 payload->setSizeCovMatrix(toTMatrixDSym(sizeMat));
36 TObject* obj = static_cast<TObject*>(payload);
37 return obj;
38}
39
40
41
42
43/* Main calibration method calling dedicated functions */
44CalibrationAlgorithm::EResult BeamSpotAlgorithm::calibrate()
45{
46 TTree* tracks = getObjectPtr<TTree>("events").get();
47 return runCalibration(tracks, "BeamSpot", getEvents,
1
Calling 'runCalibration<std::vector<Belle2::BeamSpotCalib::Event> (*)(TTree *), std::tuple<std::vector<Eigen::Matrix<double, -1, 1>>, std::vector<Eigen::Matrix<double, -1, -1>>, Eigen::Matrix<double, -1, -1>> (*)(std::vector<Belle2::BeamSpotCalib::Event>, const std::vector<double> &)>'
→
48 runBeamSpotAnalysis, getBeamSpotObj,
49 m_lossFunctionOuter, m_lossFunctionInner);
50}
←

include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h

1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9
10
11#pragma once
12
13#include <reconstruction/calibration/BeamSpotBoostInvMass/Splitter.h>
14#include <TMatrixDSym.h>
15#include <functional>
16#include <map>
17
18#include <framework/database/EventDependency.h>
19#include <framework/datastore/StoreArray.h>
20#include <framework/geometry/B2Vector3.h>
21#include <calibration/CalibrationAlgorithm.h>
22
23#include <Eigen/Dense>
24
25namespace Belle2 {
26
27 // General functions to perform the calibration
28 // Notice that the goal of the calibration is to estimate the parameters
29 // of the Gaussian distribution: center + covariance matrix describing the spread.
30 // In general it requires more data to determine the spread, so there can be
31 // several calib. subintervals with different values of
32 // the center position of the Gaussian (mean) but with identical spread parameters.
33 // The longer intervals of constant spread are called "intervals"
34 // The shorter intervals of constant mean value are called "subintervals"
35 // In general there are several subintervals within single interval
36 // By definition a subinterval always belongs only to single interval.
37
38
39 /** Function that converts Eigen symmetric matrix to ROOT matrix */
40 inline TMatrixDSym toTMatrixDSym(Eigen::MatrixXd mIn)
41 {
42 TMatrixDSym mOut(mIn.rows());
43 for (int i = 0; i < mIn.rows(); ++i)
44 for (int j = 0; j < mIn.cols(); ++j)
45 mOut(i, j) = (mIn(i, j) + mIn(j, i)) / 2.;
46 return mOut;
47 }
48
49 /** Function that converts Eigen vector to ROOT vector */
50 inline B2Vector3D toB2Vector3(Eigen::VectorXd vIn)
51 {
52 return B2Vector3D(vIn(0), vIn(1), vIn(2));
53 }
54
55 /** get id of the time point t */
56 inline int getID(const std::vector<double>& breaks, double t)
57 {
58 for (int i = 0; i < int(breaks.size()) + 1; ++i) {
59 double s = (i == 0) ? 0 : breaks[i - 1];
60 double e = (i == int(breaks.size())) ? 1e20 : breaks[i];
61 if (s <= t && t < e)
62 return i;
63 }
64 return -1;
65 }
66
67 /** The parameters related to single calibration interval */
68 struct CalibPars {
69 std::vector<Eigen::VectorXd> cnt; ///< vector of means for each calib. subinterval
70 std::vector<Eigen::MatrixXd> cntUnc; ///< vector of uncertainties of means for each calib. subinterval
71 Eigen::MatrixXd spreadMat; ///< spread CovMatrix
72
73 double spreadUnc = std::numeric_limits<double>::quiet_NaN(); ///< stat uncertainty of the spread (for eCMS)
74 double shift =
75 std::numeric_limits<double>::quiet_NaN(); ///< difference between eCMS for hadronic B decay method and mumu method, i.e. hadB - mumu
76 double shiftUnc = std::numeric_limits<double>::quiet_NaN(); ///< stat uncertainty of the shift
77 std::vector<double> pulls; ///< vector of pulls between mumu and hadB methods (for eCMS)
78 int size() const {return cnt.size();} ///< number of the subintervals
79 };
80
81
82 /** Parameters and data relevant for single calibration interval */
83 struct CalibrationData {
84 /** vector of the start and end times of the calibration subintervals */
85 std::vector<std::map<ExpRun, std::pair<double, double>>> subIntervals;
86
87 std::vector<ExpRunEvt> breakPoints; ///< vector with break points positions
88
89 CalibPars pars; ///< The parameters of the calibration itself
90
91 bool isCalibrated = false; ///< true if calibration run was successful
92
93 };
94
95
96
97 /** Extrapolate calibration to intervals where it failed */
98 inline void extrapolateCalibration(std::vector<CalibrationData>& calVec)
99 {
100 //put closest neighbor, where the statistic was low or algo failed
101 for (unsigned i = 0; i < calVec.size(); ++i) {
102 if (calVec[i].pars.cnt.size() != 0) continue;
103 const auto& r = calVec[i].subIntervals;
104 double Start, End;
105 std::tie(Start, End) = Splitter::getStartEnd(r);
106
107 Eigen::Vector3d ipNow;
108 Eigen::MatrixXd ipeNow;
109 Eigen::MatrixXd sizeMatNow;
110
111 double distMin = 1e20;
112 //Find the closest calibrated interval
113 for (unsigned j = 0; j < calVec.size(); ++j) {
114 if (calVec[j].isCalibrated == false) continue; //skip not-calibrated intervals
115 const auto& rJ = calVec[j].subIntervals;
116 for (unsigned jj = 0; jj < rJ.size(); ++jj) { //loop over subintervals
117 const auto& rNow = rJ[jj];
118 double s = rNow.begin()->second.first;
119 double e = rNow.rbegin()->second.second;
120
121 double dist1 = (s - End >= 0) ? (s - End) : 1e20;
122 double dist2 = (Start - e >= 0) ? (Start - e) : 1e20;
123 double dist = std::min(dist1, dist2);
124
125 if (dist < distMin) {
126 ipNow = calVec[j].pars.cnt.at(jj);
127 ipeNow = calVec[j].pars.cntUnc.at(jj);
128 sizeMatNow = calVec[j].pars.spreadMat;
129 distMin = dist;
130 }
131 }
132 }
133
134 //Store it to vectors
135 calVec[i].pars.cnt.resize(r.size());
136 calVec[i].pars.cntUnc.resize(r.size());
137 for (unsigned ii = 0; ii < r.size(); ++ii) {
138 calVec[i].pars.cnt.at(ii) = ipNow;
139 calVec[i].pars.cntUnc.at(ii) = ipeNow;
140 }
141 calVec[i].pars.spreadMat = sizeMatNow;
142 }
143
144 }
145
146 /** Extrapolate calibration to the very short runs which were filtered before */
147 inline void addShortRun(std::vector<CalibrationData>& calVec, std::pair<ExpRun, std::pair<double, double>> shortRun)
148 {
149 double shortStart = shortRun.second.first;
150 double shortEnd = shortRun.second.second;
151
152 double distMin = 1e20;
153 int iMin = -1, jMin = -1;
154
155 for (unsigned i = 0; i < calVec.size(); ++i) {
156 if (calVec[i].isCalibrated == false)
157 continue;
158 for (unsigned j = 0; j < calVec[i].subIntervals.size(); ++j) {
159 for (auto I : calVec[i].subIntervals[j]) {
160 double s = I.second.first;
161 double e = I.second.second;
162
163 double dist1 = (s - shortEnd >= 0) ? (s - shortEnd) : 1e20;
164 double dist2 = (shortStart - e >= 0) ? (shortStart - e) : 1e20;
165 double dist = std::min(dist1, dist2);
166
167 if (dist < distMin) {
168 distMin = dist;
169 iMin = i;
170 jMin = j;
171 }
172 }
173 }
174 }
175
176 B2ASSERT("Must be found", iMin != -1 && jMin != -1)do { if (!(iMin != -1 && jMin != -1)) { do { { LogVariableStream
varStream; varStream << "Must be found"; Belle2::LogSystem
::Instance().sendMessage(Belle2::LogMessage(Belle2::LogConfig
::c_Fatal, std::move(varStream), "reconstruction", __PRETTY_FUNCTION__
, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 176, 0)); }; exit(1); } while(false); } } while(false);
177 calVec[iMin].subIntervals[jMin].insert(shortRun);
178 }
179
180 /** Encode integer num into double val such that val is nearly not changed (maximally by a relative shift 1e-6).
181 * It is use to store time information to the payloads */
182 inline double encodeNumber(double val, unsigned num)
183 {
184 double factor = pow(FLT_RADIX2, DBL_MANT_DIG53);
185 static const long long fEnc = 1 << 32; // pow(2, 32), 32 binary digits for encoded number
17
←
The result of left shift is undefined because the right operand '32' is not smaller than 32, the capacity of 'int'
186
187 int e; //exponent of the number
188 double mantisa = std::frexp(val, &e);
189 long long mantisaI = mantisa * factor; //mantissa as integer
190
191 if (val != 0)
192 mantisaI = (mantisaI / fEnc) * fEnc + num; //adding encoded number to last digits of mantissa
193 else {
194 mantisaI = factor / 2 + num;
195 e = -100; //if the val is zero, ensure very small number by the exponent
196 }
197
198 double newVal = ldexp(mantisaI / factor, e);
199
200 return newVal;
201 }
202
203 /** Decode the integer number encoded in val */
204 inline unsigned decodeNumber(double val)
205 {
206 double factor = pow(FLT_RADIX2, DBL_MANT_DIG53);
207 static const long long fEnc = 1 << 32; // pow(2, 32), 32 binary digits for encoded number
208
209 int e;
210 double mantisa = std::frexp(val, &e);
211 long long mantisaI = mantisa * factor;
212
213 return (mantisaI % fEnc);
214 }
215
216
217
218
219 /** Store payloads to files */
220 template<typename Evt>
221 inline void storePayloads(const std::vector<Evt>& evts, const std::vector<CalibrationData>& calVecConst, std::string objName,
222 std::function<TObject*(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj)
223 {
224 auto calVec = calVecConst;
225
226 // Loop to store payloads
227 ExpRun exprunLast(-1, -1); //last exprun
228 EventDependency* intraRun = nullptr;
229
230 // Loop over calibration intervals
231 for (unsigned i = 0; i < calVec.size(); ++i) {
10
←
Assuming the condition is true
→
11
←
Loop condition is true. Entering loop body
→
232 const auto& r = calVec[i].subIntervals; // splits[i];
233 // Loop over calibration subintervals
234 for (int k = 0; k < int(r.size()); ++k) {
12
←
Assuming the condition is true
→
13
←
Loop condition is true. Entering loop body
→
235
236 for (auto I : r[k]) { //interval required to be within single run
237 ExpRun exprun = I.first;
238
239 //Encode Start+End time in seconds of the payload
240 if (calVec[i].pars.cntUnc.at(k).rows() == 3) {
14
←
Assuming the condition is false
→
15
←
Taking false branch
→
241 calVec[i].pars.cntUnc.at(k)(0, 1) = calVec[i].pars.cntUnc.at(k)(1, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 1),
242 round(I.second.first * 3600));
243 calVec[i].pars.cntUnc.at(k)(0, 2) = calVec[i].pars.cntUnc.at(k)(2, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 2),
244 round(I.second.second * 3600));
245 } else {
246 calVec[i].pars.cntUnc.at(k)(0, 0) = encodeNumber(calVec[i].pars.cntUnc.at(k)(0, 0), round(I.second.first * 3600));
16
←
Calling 'encodeNumber'
→
247 calVec[i].pars.spreadMat(0, 0) = encodeNumber(calVec[i].pars.spreadMat(0, 0), round(I.second.second * 3600));
248 }
249
250 TObject* obj = getCalibObj(calVec[i].pars.cnt.at(k), calVec[i].pars.cntUnc.at(k), calVec[i].pars.spreadMat);
251 if (exprun != exprunLast) { //if new run
252 if (intraRun) { //if not first -> store
253 auto m_iov = IntervalOfValidity(exprunLast.exp, exprunLast.run, exprunLast.exp, exprunLast.run);
254 Database::Instance().storeData(objName, intraRun, m_iov);
255 }
256
257 intraRun = new EventDependency(obj);
258 } else {
259 int breakPoint;
260 if (k - 1 >= 0) {
261 breakPoint = calVec[i].breakPoints.at(k - 1).evt;
262 B2ASSERT("Payload saving consistency", calVec[i].breakPoints.at(k - 1).run == exprun.run)do { if (!(calVec[i].breakPoints.at(k - 1).run == exprun.run)
) { do { { LogVariableStream varStream; varStream << "Payload saving consistency"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 262, 0)); }; exit(1); } while(false); } } while(false);
263 } else {
264 B2ASSERT("Payload saving consistency", i != 0)do { if (!(i != 0)) { do { { LogVariableStream varStream; varStream
<< "Payload saving consistency"; Belle2::LogSystem::Instance
().sendMessage(Belle2::LogMessage(Belle2::LogConfig::c_Fatal,
std::move(varStream), "reconstruction", __PRETTY_FUNCTION__,
"include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 264, 0)); }; exit(1); } while(false); } } while(false);
265 double rStart, rEnd;
266 std::tie(rStart, rEnd) = Splitter::getStartEnd(r);
267 auto pos = getPosition(evts, rStart);
268 breakPoint = pos.evt;
269 B2ASSERT("Payload saving consistency", pos.run == exprun.run)do { if (!(pos.run == exprun.run)) { do { { LogVariableStream
varStream; varStream << "Payload saving consistency"; Belle2
::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2
::LogConfig::c_Fatal, std::move(varStream), "reconstruction",
__PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 269, 0)); }; exit(1); } while(false); } } while(false);
270 }
271 intraRun->add(breakPoint, obj);
272 }
273 exprunLast = exprun;
274 }
275 } //end loop over calibration subintervals
276
277 } //end loop over calibration intervals
278
279 //Store the last entry
280 auto m_iov = IntervalOfValidity(exprunLast.exp, exprunLast.run, exprunLast.exp, exprunLast.run);
281 Database::Instance().storeData(objName, intraRun, m_iov);
282 }
283
284
285 /** Store payloads to files, where calib data have no intra-run dependence */
286 inline void storePayloadsNoIntraRun(const std::vector<CalibrationData>& calVecConst, std::string objName,
287 std::function<TObject*(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd) > getCalibObj)
288 {
289 auto calVec = calVecConst;
290
291 // Check that there is no intra-run dependence
292 std::set<ExpRun> existingRuns;
293 for (unsigned i = 0; i < calVec.size(); ++i) {
294 const auto& r = calVec[i].subIntervals;
295 // Loop over calibration subintervals
296 for (int k = 0; k < int(r.size()); ++k) {
297
298 for (auto I : r[k]) {
299 ExpRun exprun = I.first;
300 // make sure that the run isn't already in the list, to avoid duplicity
301 if (existingRuns.count(exprun) != 0)
302 B2FATAL("Intra-run dependence exists")do { { LogVariableStream varStream; varStream << "Intra-run dependence exists"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 302, 0)); }; exit(1); } while(false);
303 existingRuns.insert(exprun);
304 }
305 }
306 }
307
308
309 // Loop over calibration intervals
310 for (unsigned i = 0; i < calVec.size(); ++i) {
311 const auto& r = calVec[i].subIntervals; // splits[i];
312 // Loop over calibration subintervals
313 for (unsigned k = 0; k < r.size(); ++k) {
314
315 TObject* obj = getCalibObj(calVec[i].pars.cnt.at(k), calVec[i].pars.cntUnc.at(k), calVec[i].pars.spreadMat);
316
317 ExpRun start = (r[k].cbegin()->first);
318 ExpRun last = (r[k].crbegin()->first);
319
320 auto iov = IntervalOfValidity(start.exp, start.run, last.exp, last.run);
321 Database::Instance().storeData(objName, obj, iov);
322
323
324 } //end loop over calibration subintervals
325 } //end loop over calibration intervals
326
327 }
328
329 /** run calibration algorithm for single calibration interval */
330 template<typename Evt, typename Fun>
331 inline CalibrationData runAlgorithm(const std::vector<Evt>& evts, std::vector<std::map<ExpRun, std::pair<double, double>>> range,
332 Fun runCalibAnalysis
333 )
334 {
335 CalibrationData calD;
336 auto& r = range;
337 double rStart, rEnd;
338 std::tie(rStart, rEnd) = Splitter::getStartEnd(r);
339 B2INFO("Start of loop startTime endTime : " << rStart << " " << rEnd)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Start of loop startTime endTime : "
<< rStart << " " << rEnd; Belle2::LogSystem
::Instance().sendMessage(Belle2::LogMessage(Belle2::LogConfig
::c_Info, std::move(varStream), "reconstruction", __PRETTY_FUNCTION__
, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 339, 0)); }; } } while(false);
340
341 auto breaks = Splitter::getBreaks(r);
342
343 std::vector<Evt> evtsNow;
344
345 std::vector<int> Counts(breaks.size() + 1, 0);
346 // Select events belonging to the interval
347 for (const auto& ev : evts) {
348 if (rStart <= ev.t && ev.t < rEnd) {
349 evtsNow.push_back(ev);
350 ++Counts.at(getID(breaks, ev.t));
351 }
352 }
353
354 B2ASSERT("Number of intervals vs number of breakPoints", r.size() == breaks.size() + 1)do { if (!(r.size() == breaks.size() + 1)) { do { { LogVariableStream
varStream; varStream << "Number of intervals vs number of breakPoints"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 354, 0)); }; exit(1); } while(false); } } while(false);
355
356 //Merge smallest interval if with low stat (try it 10times)
357 for (int k = 0; k < 10; ++k) {
358 int iMin = min_element(Counts.begin(), Counts.end()) - Counts.begin();
359 if (Counts.size() >= 2 && Counts[iMin] < 50) { //merge with neighbor if possible
360 auto iM = -1;
361 if (iMin == 0)
362 iM = iMin + 1;
363 else if (iMin == int(Counts.size()) - 1)
364 iM = iMin - 1;
365 else {
366 if (Counts[iMin + 1] < Counts[iMin - 1])
367 iM = iMin + 1;
368 else
369 iM = iMin - 1;
370 }
371 B2ASSERT("Number of intervals equal to size of counters", r.size() == Counts.size())do { if (!(r.size() == Counts.size())) { do { { LogVariableStream
varStream; varStream << "Number of intervals equal to size of counters"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 371, 0)); }; exit(1); } while(false); } } while(false);
372
373 r.at(iM) = Splitter::mergeIntervals(r[iM], r[iMin]);
374 r.erase(r.begin() + iMin);
375 breaks = Splitter::getBreaks(r);
376 Counts[iM] += Counts[iMin];
377 Counts.erase(Counts.begin() + iMin);
378 }
379 }
380
381 B2INFO("#events " << " : " << evtsNow.size())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "#events " << " : " <<
evtsNow.size(); Belle2::LogSystem::Instance().sendMessage(Belle2
::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream),
"reconstruction", __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 381, 0)); }; } } while(false);
382 B2INFO("Breaks size " << " : " << breaks.size())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Breaks size " << " : " <<
breaks.size(); Belle2::LogSystem::Instance().sendMessage(Belle2
::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream),
"reconstruction", __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 382, 0)); }; } } while(false);
383
384 calD.breakPoints = convertSplitPoints(evtsNow, breaks);
385
386 calD.subIntervals = r;
387
388 if (breaks.size() > 0)
389 B2INFO("StartOfCalibInterval (run,evtNo,vtxIntervalsSize) " << calD.breakPoints.at(0).run << " " <<do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "StartOfCalibInterval (run,evtNo,vtxIntervalsSize) "
<< calD.breakPoints.at(0).run << " " << calD
.breakPoints.at(0).evt << " " << calD.breakPoints
.size(); Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Info, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 390, 0)); }; } } while(false)
390 calD.breakPoints.at(0).evt << " " << calD.breakPoints.size())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "StartOfCalibInterval (run,evtNo,vtxIntervalsSize) "
<< calD.breakPoints.at(0).run << " " << calD
.breakPoints.at(0).evt << " " << calD.breakPoints
.size(); Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Info, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 390, 0)); }; } } while(false);
391
392
393 //If too few events, let have the output empty
394 //Will be filled with the closest neighbor at the next stage
395 if (evtsNow.size() < 50) {
396 return calD;
397 }
398
399 // Run the calibration
400 B2INFO("Start of running calibration over calibration interval")do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Start of running calibration over calibration interval"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Info, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 400, 0)); }; } } while(false);
401 tie(calD.pars.cnt, calD.pars.cntUnc, calD.pars.spreadMat) = runCalibAnalysis(evtsNow, breaks);
402 calD.pars.pulls.resize(calD.pars.cnt.size());
403 B2INFO("End of running analysis - SpreadMatX : " << sqrt(abs(calD.pars.spreadMat(0, 0))))do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "End of running analysis - SpreadMatX : "
<< sqrt(abs(calD.pars.spreadMat(0, 0))); Belle2::LogSystem
::Instance().sendMessage(Belle2::LogMessage(Belle2::LogConfig
::c_Info, std::move(varStream), "reconstruction", __PRETTY_FUNCTION__
, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 403, 0)); }; } } while(false);
404 B2ASSERT("All subintervals have calibration of the mean value", calD.pars.cnt.size() == r.size())do { if (!(calD.pars.cnt.size() == r.size())) { do { { LogVariableStream
varStream; varStream << "All subintervals have calibration of the mean value"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 404, 0)); }; exit(1); } while(false); } } while(false);
405 B2ASSERT("All subintervals have calibration of the unc. of mean", calD.pars.cntUnc.size() == r.size())do { if (!(calD.pars.cntUnc.size() == r.size())) { do { { LogVariableStream
varStream; varStream << "All subintervals have calibration of the unc. of mean"
; Belle2::LogSystem::Instance().sendMessage(Belle2::LogMessage
(Belle2::LogConfig::c_Fatal, std::move(varStream), "reconstruction"
, __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 405, 0)); }; exit(1); } while(false); } } while(false);
406
407 calD.isCalibrated = true;
408
409 return calD;
410 }
411
412
413 /** Run the the calibration over the whole event sample
414 @param tracks: TTree object with mu-mu events
415 @param calibName: name of the calibration payload
416 @param GetEvents: function that transforms TTree to std::vector
417 @param calibAnalysis: function that performs the calibration on a single calibration interval
418 @param calibObjCreator: function that stores results to the payload class which inherits from TObject
419 @param m_lossFunctionOuter: Lost function for the calibration intervals of the spread parameters
420 @param m_lossFunctionInner: Lost function for the calibration subintervals (for the mean value parameters)
421 @return State of the calibration run, i.e. EResult::c_OK if everything OK
422 */
423 template<typename Fun1, typename Fun2>
424 CalibrationAlgorithm::EResult runCalibration(TTree* tracks, const std::string& calibName, Fun1 GetEvents, Fun2 calibAnalysis,
425 std::function<TObject*(Eigen::VectorXd, Eigen::MatrixXd, Eigen::MatrixXd)> calibObjCreator,
426 TString m_lossFunctionOuter, TString m_lossFunctionInner)
427 {
428 // Check that there are at least some data
429 if (!tracks || tracks->GetEntries() < 15) {
2
←
Assuming 'tracks' is non-null
→
3
←
Assuming the condition is false
→
4
←
Taking false branch
→
430 if (tracks)
431 B2WARNING("Too few data : " << tracks->GetEntries())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Warning, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Too few data : " << tracks
->GetEntries(); Belle2::LogSystem::Instance().sendMessage(
Belle2::LogMessage(Belle2::LogConfig::c_Warning, std::move(varStream
), "reconstruction", __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 431, 0)); }; } } while(false);
432 return CalibrationAlgorithm::EResult::c_NotEnoughData;
433 }
434 B2INFO("Number of tracks: " << tracks->GetEntries())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Info, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Number of tracks: " << tracks
->GetEntries(); Belle2::LogSystem::Instance().sendMessage(
Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream
), "reconstruction", __PRETTY_FUNCTION__, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 434, 0)); }; } } while(false);
5
←
Taking false branch
→
6
←
Loop condition is false. Exiting loop
→
435
436 // Tree to vector of Events
437 auto evts = GetEvents(tracks);
438
439 //Time range for each ExpRun
440 std::map<ExpRun, std::pair<double, double>> runsInfoOrg = getRunInfo(evts);
441 std::map<ExpRun, std::pair<double, double>> runsRemoved; //map with time intervals of very short runs
442 auto runsInfo = filter(runsInfoOrg, 2. / 60, runsRemoved); //include only runs longer than 2mins
443
444 // If nothing remains
445 if (runsInfo.size() == 0) {
7
←
Assuming the condition is false
→
8
←
Taking false branch
→
446 B2WARNING("Too short run")do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2::
LogConfig::c_Warning, 0, "reconstruction")) { { LogVariableStream
varStream; varStream << "Too short run"; Belle2::LogSystem
::Instance().sendMessage(Belle2::LogMessage(Belle2::LogConfig
::c_Warning, std::move(varStream), "reconstruction", __PRETTY_FUNCTION__
, "include/reconstruction/calibration/BeamSpotBoostInvMass/calibTools.h"
, 446, 0)); }; } } while(false);
447 return CalibrationAlgorithm::EResult::c_NotEnoughData;
448 }
449
450 // Get intervals based on the input loss functions
451 Splitter splt;
452 auto splits = splt.getIntervals(runsInfo, evts, m_lossFunctionOuter, m_lossFunctionInner);
453
454 //Loop over all calibration intervals
455 std::vector<CalibrationData> calVec;
456 for (auto s : splits) {
457 CalibrationData calD = runAlgorithm(evts, s, calibAnalysis); // run the calibration over the interval s
458 calVec.push_back(calD);
459 }
460
461 // extrapolate results to the low-stat intervals
462 extrapolateCalibration(calVec);
463
464 // Include removed short runs
465 for (auto shortRun : runsRemoved) {
466 addShortRun(calVec, shortRun);
467 }
468
469 // Store Payloads to files
470 storePayloads(evts, calVec, calibName, calibObjCreator);
9
←
Calling 'storePayloads<Belle2::BeamSpotCalib::Event>'
→
471
472 return CalibrationAlgorithm::EResult::c_OK;
473 }
474
475
476}