| File: | analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc |
| Warning: | line 177, column 39 The left operand of '*' is a garbage value |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 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 | #include "analysis/OrcaKinFit/SoftGaussParticleConstraint.h" | |||
| 18 | #include "analysis/OrcaKinFit/ParticleFitObject.h" | |||
| 19 | #include <framework/logging/Logger.h> | |||
| 20 | #include <iostream> | |||
| 21 | #include <cmath> | |||
| 22 | using namespace std; | |||
| 23 | ||||
| 24 | namespace Belle2 { | |||
| 25 | namespace OrcaKinFit { | |||
| 26 | ||||
| 27 | SoftGaussParticleConstraint::SoftGaussParticleConstraint(double sigma_) | |||
| 28 | : sigma(sigma_) | |||
| 29 | { | |||
| 30 | invalidateCache(); | |||
| 31 | } | |||
| 32 | ||||
| 33 | double SoftGaussParticleConstraint::getSigma() const | |||
| 34 | { | |||
| 35 | return sigma; | |||
| 36 | } | |||
| 37 | ||||
| 38 | double SoftGaussParticleConstraint::setSigma(double sigma_) | |||
| 39 | { | |||
| 40 | if (sigma_ != 0) sigma = std::abs(sigma_); | |||
| 41 | return sigma; | |||
| 42 | } | |||
| 43 | ||||
| 44 | double SoftGaussParticleConstraint::getChi2() const | |||
| 45 | { | |||
| 46 | double r = getValue() / getSigma(); | |||
| 47 | return r * r; | |||
| 48 | } | |||
| 49 | ||||
| 50 | double SoftGaussParticleConstraint::getError() const | |||
| 51 | { | |||
| 52 | double dgdpi[4]; | |||
| 53 | double error2 = 0; | |||
| 54 | for (unsigned int i = 0; i < fitobjects.size(); ++i) { | |||
| 55 | const ParticleFitObject* foi = fitobjects[i]; | |||
| 56 | assert(foi)(static_cast <bool> (foi) ? void (0) : __assert_fail ("foi" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 56 , __extension__ __PRETTY_FUNCTION__)); | |||
| 57 | if (firstDerivatives(i, dgdpi)) { | |||
| 58 | error2 += foi->getError2(dgdpi, getVarBasis()); | |||
| 59 | } | |||
| 60 | } | |||
| 61 | return std::sqrt(std::abs(error2)); | |||
| 62 | } | |||
| 63 | ||||
| 64 | /** | |||
| 65 | * Calculates the second derivative of the constraint g w.r.t. the various parameters | |||
| 66 | * and adds it to the global covariance matrix | |||
| 67 | * | |||
| 68 | * We denote with P_i the 4-vector of the i-th ParticleFitObject, | |||
| 69 | * then | |||
| 70 | * \f[ | |||
| 71 | * \frac{\partial ^2 g}{\partial a_k \partial a_l} | |||
| 72 | * = \sum_i \sum_j \frac{\partial ^2 g}{\partial P_i \partial P_j} \cdot | |||
| 73 | * \frac{\partial P_i}{\partial a_k} \cdot \frac{\partial P_j}{\partial a_l} | |||
| 74 | * + \sum_i \frac{\partial g}{\partial P_i} \cdot | |||
| 75 | * \frac{\partial^2 P_i}{\partial a_k \partial a_l} | |||
| 76 | * \f] | |||
| 77 | * Here, \f$\frac{\partial P_i}{\partial a_k}\f$ is a \f$4 \times n_i\f$ Matrix, where | |||
| 78 | * \f$n_i\f$ is the number of parameters of FitObject i; | |||
| 79 | * Correspondingly, \f$\frac{\partial^2 P_i}{\partial a_k \partial a_l}\f$ is a | |||
| 80 | * \f$4 \times n_i \times n_i\f$ matrix. | |||
| 81 | * Also, \f$\frac{\partial ^2 g}{\partial P_i \partial P_j}\f$ is a \f$4\times 4\f$ matrix | |||
| 82 | * for a given i and j, and \f$\frac{\partial g}{\partial P_i}\f$ is a 4-vector | |||
| 83 | * (though not a Lorentz-vector!). | |||
| 84 | * | |||
| 85 | */ | |||
| 86 | ||||
| 87 | ||||
| 88 | void SoftGaussParticleConstraint::add2ndDerivativesToMatrix(double* M, int idim) const | |||
| 89 | { | |||
| 90 | ||||
| 91 | /** First, treat the part | |||
| 92 | * \f[ | |||
| 93 | * \frac{\partial ^2 g}{\partial P_i \partial P_j} \cdot | |||
| 94 | * \frac{\partial P_i}{\partial a_k} \cdot \frac{\partial P_j}{\partial a_l} | |||
| 95 | * \f] | |||
| 96 | */ | |||
| 97 | double s = getSigma(); | |||
| 98 | double fact = 2 * getValue() / (s * s); | |||
| 99 | ||||
| 100 | // Derivatives \f$\frac{\partial ^2 g}{\partial P_i \partial P_j}\f$ at fixed i, j | |||
| 101 | // d2GdPidPj[4*ii+jj] is derivative w.r.t. P_i,ii and P_j,jj, where ii=0,1,2,3 for E,px,py,pz | |||
| 102 | double d2GdPidPj[16]; | |||
| 103 | // Derivatives \f$\frac {\partial P_i}{\partial a_k}\f$ for all i; | |||
| 104 | // k is local parameter number | |||
| 105 | // dPidAk[KMAX*4*i + 4*k + ii] is \f$\frac {\partial P_{i,ii}}{\partial a_k}\f$, | |||
| 106 | // with ii=0, 1, 2, 3 for E, px, py, pz | |||
| 107 | const int KMAX = 4; | |||
| 108 | const int n = fitobjects.size(); | |||
| 109 | auto* dPidAk = new double[n * KMAX * 4]; | |||
| 110 | bool* dPidAkval = new bool[n]; | |||
| 111 | ||||
| 112 | for (int i = 0; i
| |||
| ||||
| 113 | ||||
| 114 | // Derivatives \f$\frac{\partial ^2 g}{\partial P_i \partial a_l}\f$ at fixed i | |||
| 115 | // d2GdPdAl[4*l + ii] is \f$\frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l}\f$ | |||
| 116 | double d2GdPdAl[4 * KMAX]; | |||
| 117 | // Derivatives \f$\frac{\partial ^2 g}{\partial a_k \partial a_l}\f$ | |||
| 118 | double d2GdAkdAl[KMAX * KMAX] = {0}; | |||
| 119 | ||||
| 120 | // Global parameter numbers: parglobal[KMAX*i+klocal] | |||
| 121 | // is global parameter number of local parameter klocal of i-th Fit object | |||
| 122 | int* parglobal = new int[KMAX * n]; | |||
| 123 | ||||
| 124 | for (int i = 0; i < n; ++i) { | |||
| 125 | const ParticleFitObject* foi = fitobjects[i]; | |||
| 126 | assert(foi)(static_cast <bool> (foi) ? void (0) : __assert_fail ("foi" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 126 , __extension__ __PRETTY_FUNCTION__)); | |||
| 127 | for (int klocal = 0; klocal < foi->getNPar(); ++klocal) { | |||
| 128 | parglobal [KMAX * i + klocal] = foi->getGlobalParNum(klocal); | |||
| 129 | } | |||
| 130 | } | |||
| 131 | ||||
| 132 | ||||
| 133 | for (int i = 0; i < n; ++i) { | |||
| 134 | const ParticleFitObject* foi = fitobjects[i]; | |||
| 135 | assert(foi)(static_cast <bool> (foi) ? void (0) : __assert_fail ("foi" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 135 , __extension__ __PRETTY_FUNCTION__)); | |||
| 136 | for (int j = 0; j < n; ++j) { | |||
| 137 | const ParticleFitObject* foj = fitobjects[j]; | |||
| 138 | assert(foj)(static_cast <bool> (foj) ? void (0) : __assert_fail ("foj" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 138 , __extension__ __PRETTY_FUNCTION__)); | |||
| 139 | if (secondDerivatives(i, j, d2GdPidPj)) { | |||
| 140 | if (!dPidAkval[i]) { | |||
| 141 | foi->getDerivatives(dPidAk + i * (KMAX * 4), KMAX * 4); | |||
| 142 | dPidAkval[i] = true; | |||
| 143 | } | |||
| 144 | if (!dPidAkval[j]) { | |||
| 145 | foj->getDerivatives(dPidAk + j * (KMAX * 4), KMAX * 4); | |||
| 146 | dPidAkval[j] = true; | |||
| 147 | } | |||
| 148 | // Now sum over E/px/Py/Pz for object j: | |||
| 149 | // \f[ | |||
| 150 | // \frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l} | |||
| 151 | // = (sum_{j}) sum_{jj} frac{\partial ^2 g}{\partial P_{i,ii} \partial P_{j,jj}} | |||
| 152 | // \cdot \frac{\partial P_{j,jj}}{\partial a_l} | |||
| 153 | // \f] | |||
| 154 | // We're summing over jj here | |||
| 155 | for (int llocal = 0; llocal < foj->getNPar(); ++llocal) { | |||
| 156 | for (int ii = 0; ii < 4; ++ii) { | |||
| 157 | int ind1 = 4 * ii; | |||
| 158 | int ind2 = (KMAX * 4) * j + 4 * llocal; | |||
| 159 | double& r = d2GdPdAl[4 * llocal + ii]; | |||
| 160 | r = d2GdPidPj[ ind1] * dPidAk[ ind2]; // E | |||
| 161 | r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // px | |||
| 162 | r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // py | |||
| 163 | r += d2GdPidPj[++ind1] * dPidAk[++ind2]; // pz | |||
| 164 | } | |||
| 165 | } | |||
| 166 | // Now sum over E/px/Py/Pz for object i, i.e. sum over ii: | |||
| 167 | // \f[ | |||
| 168 | // \frac{\partial ^2 g}{\partial a_k \partial a_l} | |||
| 169 | // = \sum_{ii} \frac{\partial ^2 g}{\partial P_{i,ii} \partial a_l} \cdot | |||
| 170 | // \frac{\partial P_{i,ii}}{\partial a_k} | |||
| 171 | // \f] | |||
| 172 | for (int klocal = 0; klocal < foi->getNPar(); ++klocal) { | |||
| 173 | for (int llocal = 0; llocal < foj->getNPar(); ++llocal) { | |||
| 174 | int ind1 = 4 * llocal; | |||
| 175 | int ind2 = (KMAX * 4) * i + 4 * klocal; | |||
| 176 | double& r = d2GdAkdAl[KMAX * klocal + llocal]; | |||
| 177 | r = d2GdPdAl[ ind1] * dPidAk[ ind2]; //E | |||
| ||||
| 178 | r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // px | |||
| 179 | r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // py | |||
| 180 | r += d2GdPdAl[++ind1] * dPidAk[++ind2]; // pz | |||
| 181 | } | |||
| 182 | } | |||
| 183 | // Now expand the local parameter numbers to global ones | |||
| 184 | for (int klocal = 0; klocal < foi->getNPar(); ++klocal) { | |||
| 185 | int kglobal = parglobal [KMAX * i + klocal]; | |||
| 186 | for (int llocal = 0; llocal < foj->getNPar(); ++llocal) { | |||
| 187 | int lglobal = parglobal [KMAX * j + llocal]; | |||
| 188 | M [idim * kglobal + lglobal] += fact * d2GdAkdAl[KMAX * klocal + llocal]; | |||
| 189 | } | |||
| 190 | } | |||
| 191 | } | |||
| 192 | } | |||
| 193 | } | |||
| 194 | /** Second, treat the parts | |||
| 195 | * \f[ | |||
| 196 | * \sum_i \frac{\partial g}{\partial P_i} \cdot | |||
| 197 | * \frac{\partial^2 P_i}{\partial a_k \partial a_l} | |||
| 198 | * \f] | |||
| 199 | * and | |||
| 200 | * \f[ | |||
| 201 | * \frac{\partial^2 h}{\partial g^2} | |||
| 202 | * \sum_i \frac{\partial g}{\partial P_i} \cdot | |||
| 203 | * \frac{\partial P_i}{\partial a_k} | |||
| 204 | * \sum_j \frac{\partial g}{\partial P_j} \cdot | |||
| 205 | * \frac{\partial P_j}{\partial a_l} | |||
| 206 | * \f] | |||
| 207 | * | |||
| 208 | * Here, \f$\frac{\partial g}{\partial P_i}\f$ is a 4-vector, which we pass on to | |||
| 209 | * the FitObject | |||
| 210 | */ | |||
| 211 | ||||
| 212 | auto* v = new double[idim]; | |||
| 213 | for (int i = 0; i < idim; ++i) v[i] = 0; | |||
| 214 | double sqrtfact2 = sqrt(2.0) / s; | |||
| 215 | ||||
| 216 | double dgdpi[4]; | |||
| 217 | for (int i = 0; i < n; ++i) { | |||
| 218 | const ParticleFitObject* foi = fitobjects[i]; | |||
| 219 | assert(foi)(static_cast <bool> (foi) ? void (0) : __assert_fail ("foi" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 219 , __extension__ __PRETTY_FUNCTION__)); | |||
| 220 | if (firstDerivatives(i, dgdpi)) { | |||
| 221 | foi->addTo2ndDerivatives(M, idim, fact, dgdpi, getVarBasis()); | |||
| 222 | foi->addToGlobalChi2DerVector(v, idim, sqrtfact2, dgdpi, getVarBasis()); | |||
| 223 | } | |||
| 224 | } | |||
| 225 | ||||
| 226 | for (int i = 0; i < idim; ++i) { | |||
| 227 | if (double vi = v[i]) { | |||
| 228 | int ioffs = i * idim; | |||
| 229 | for (double* pvj = v; pvj < v + idim; ++pvj) { | |||
| 230 | M[ioffs++] += vi * (*pvj); | |||
| 231 | } | |||
| 232 | } | |||
| 233 | } | |||
| 234 | ||||
| 235 | ||||
| 236 | delete[] dPidAk; | |||
| 237 | delete[] dPidAkval; | |||
| 238 | delete[] parglobal; | |||
| 239 | delete[] v; | |||
| 240 | } | |||
| 241 | ||||
| 242 | void SoftGaussParticleConstraint::addToGlobalChi2DerVector(double* y, int idim) const | |||
| 243 | { | |||
| 244 | double dgdpi[4]; | |||
| 245 | double s = getSigma(); | |||
| 246 | double r = 2 * getValue() / (s * s); | |||
| 247 | for (unsigned int i = 0; i < fitobjects.size(); ++i) { | |||
| 248 | const ParticleFitObject* foi = fitobjects[i]; | |||
| 249 | assert(foi)(static_cast <bool> (foi) ? void (0) : __assert_fail ("foi" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 249 , __extension__ __PRETTY_FUNCTION__)); | |||
| 250 | if (firstDerivatives(i, dgdpi)) { | |||
| 251 | foi->addToGlobalChi2DerVector(y, idim, r, dgdpi, getVarBasis()); | |||
| 252 | } | |||
| 253 | } | |||
| 254 | } | |||
| 255 | ||||
| 256 | void SoftGaussParticleConstraint::test1stDerivatives() | |||
| 257 | { | |||
| 258 | B2INFO("SoftGaussParticleConstraint::test1stDerivatives for " << getName())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "SoftGaussParticleConstraint::test1stDerivatives for " << getName(); Belle2::LogSystem::Instance().sendMessage (Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream ), "analysis", __PRETTY_FUNCTION__, "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc" , 258, 0)); }; } } while(false); | |||
| 259 | double y[100]; | |||
| 260 | for (double& i : y) i = 0; | |||
| 261 | addToGlobalChi2DerVector(y, 100); | |||
| 262 | double eps = 0.00001; | |||
| 263 | for (unsigned int ifo = 0; ifo < fitobjects.size(); ++ifo) { | |||
| 264 | ParticleFitObject* fo = fitobjects[ifo]; | |||
| 265 | assert(fo)(static_cast <bool> (fo) ? void (0) : __assert_fail ("fo" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 265 , __extension__ __PRETTY_FUNCTION__)); | |||
| 266 | for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) { | |||
| 267 | int iglobal = fo->getGlobalParNum(ilocal); | |||
| 268 | double calc = y[iglobal]; | |||
| 269 | double num = num1stDerivative(ifo, ilocal, eps); | |||
| 270 | B2INFO("fo: " << fo->getName() << " par " << ilocal << "/"do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo: " << fo->getName() << " par " << ilocal << "/" << iglobal << " (" << fo->getParamName(ilocal) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2::LogSystem::Instance().sendMessage (Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream ), "analysis", __PRETTY_FUNCTION__, "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc" , 272, 0)); }; } } while(false) | |||
| 271 | << iglobal << " (" << fo->getParamName(ilocal)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo: " << fo->getName() << " par " << ilocal << "/" << iglobal << " (" << fo->getParamName(ilocal) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2::LogSystem::Instance().sendMessage (Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream ), "analysis", __PRETTY_FUNCTION__, "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc" , 272, 0)); }; } } while(false) | |||
| 272 | << ") calc: " << calc << " - num: " << num << " = " << calc - num)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo: " << fo->getName() << " par " << ilocal << "/" << iglobal << " (" << fo->getParamName(ilocal) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2::LogSystem::Instance().sendMessage (Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream ), "analysis", __PRETTY_FUNCTION__, "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc" , 272, 0)); }; } } while(false); | |||
| 273 | } | |||
| 274 | } | |||
| 275 | } | |||
| 276 | void SoftGaussParticleConstraint::test2ndDerivatives() | |||
| 277 | { | |||
| 278 | B2INFO("SoftGaussParticleConstraint::test2ndDerivatives for " << getName())do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "SoftGaussParticleConstraint::test2ndDerivatives for " << getName(); Belle2::LogSystem::Instance().sendMessage (Belle2::LogMessage(Belle2::LogConfig::c_Info, std::move(varStream ), "analysis", __PRETTY_FUNCTION__, "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc" , 278, 0)); }; } } while(false); | |||
| 279 | const int idim = 100; | |||
| 280 | auto* M = new double[idim * idim]; | |||
| 281 | for (int i = 0; i < idim * idim; ++i) M[i] = 0; | |||
| 282 | add2ndDerivativesToMatrix(M, idim); | |||
| 283 | double eps = 0.0001; | |||
| 284 | B2INFO("eps=" << eps)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "eps=" << eps; Belle2::LogSystem:: Instance().sendMessage(Belle2::LogMessage(Belle2::LogConfig:: c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 284 , 0)); }; } } while(false); | |||
| 285 | ||||
| 286 | for (unsigned int ifo1 = 0; ifo1 < fitobjects.size(); ++ifo1) { | |||
| 287 | ParticleFitObject* fo1 = fitobjects[ifo1]; | |||
| 288 | assert(fo1)(static_cast <bool> (fo1) ? void (0) : __assert_fail ("fo1" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 288 , __extension__ __PRETTY_FUNCTION__)); | |||
| 289 | for (unsigned int ifo2 = ifo1; ifo2 < fitobjects.size(); ++ifo2) { | |||
| 290 | ParticleFitObject* fo2 = fitobjects[ifo2]; | |||
| 291 | assert(fo2)(static_cast <bool> (fo2) ? void (0) : __assert_fail ("fo2" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 291 , __extension__ __PRETTY_FUNCTION__)); | |||
| 292 | for (int ilocal1 = 0; ilocal1 < fo1->getNPar(); ++ilocal1) { | |||
| 293 | int iglobal1 = fo1->getGlobalParNum(ilocal1); | |||
| 294 | for (int ilocal2 = (ifo1 == ifo2 ? ilocal1 : 0); ilocal2 < fo2->getNPar(); ++ilocal2) { | |||
| 295 | int iglobal2 = fo2->getGlobalParNum(ilocal2); | |||
| 296 | double calc = M[idim * iglobal1 + iglobal2]; | |||
| 297 | double num = num2ndDerivative(ifo1, ilocal1, eps, ifo2, ilocal2, eps); | |||
| 298 | B2INFO("fo1: " << fo1->getName() << " par " << ilocal1 << "/"do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo1: " << fo1->getName() << " par " << ilocal1 << "/" << iglobal1 << " (" << fo1->getParamName(ilocal1) << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/" << iglobal2 << " (" << fo2-> getParamName(ilocal2) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2 ::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2 ::LogConfig::c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 302 , 0)); }; } } while(false) | |||
| 299 | << iglobal1 << " (" << fo1->getParamName(ilocal1)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo1: " << fo1->getName() << " par " << ilocal1 << "/" << iglobal1 << " (" << fo1->getParamName(ilocal1) << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/" << iglobal2 << " (" << fo2-> getParamName(ilocal2) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2 ::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2 ::LogConfig::c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 302 , 0)); }; } } while(false) | |||
| 300 | << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/"do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo1: " << fo1->getName() << " par " << ilocal1 << "/" << iglobal1 << " (" << fo1->getParamName(ilocal1) << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/" << iglobal2 << " (" << fo2-> getParamName(ilocal2) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2 ::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2 ::LogConfig::c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 302 , 0)); }; } } while(false) | |||
| 301 | << iglobal2 << " (" << fo2->getParamName(ilocal2)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo1: " << fo1->getName() << " par " << ilocal1 << "/" << iglobal1 << " (" << fo1->getParamName(ilocal1) << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/" << iglobal2 << " (" << fo2-> getParamName(ilocal2) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2 ::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2 ::LogConfig::c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 302 , 0)); }; } } while(false) | |||
| 302 | << ") calc: " << calc << " - num: " << num << " = " << calc - num)do { if (Belle2::LogSystem::Instance().isLevelEnabled(Belle2:: LogConfig::c_Info, 0, "analysis")) { { LogVariableStream varStream ; varStream << "fo1: " << fo1->getName() << " par " << ilocal1 << "/" << iglobal1 << " (" << fo1->getParamName(ilocal1) << "), fo2: " << fo2->getName() << " par " << ilocal2 << "/" << iglobal2 << " (" << fo2-> getParamName(ilocal2) << ") calc: " << calc << " - num: " << num << " = " << calc - num; Belle2 ::LogSystem::Instance().sendMessage(Belle2::LogMessage(Belle2 ::LogConfig::c_Info, std::move(varStream), "analysis", __PRETTY_FUNCTION__ , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 302 , 0)); }; } } while(false); | |||
| 303 | } | |||
| 304 | } | |||
| 305 | } | |||
| 306 | } | |||
| 307 | delete[] M; | |||
| 308 | } | |||
| 309 | ||||
| 310 | ||||
| 311 | double SoftGaussParticleConstraint::num1stDerivative(int ifo, int ilocal, double eps) | |||
| 312 | { | |||
| 313 | ParticleFitObject* fo = fitobjects[ifo]; | |||
| 314 | assert(fo)(static_cast <bool> (fo) ? void (0) : __assert_fail ("fo" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 314 , __extension__ __PRETTY_FUNCTION__)); | |||
| 315 | double save = fo->getParam(ilocal); | |||
| 316 | fo->setParam(ilocal, save + eps); | |||
| 317 | double v1 = getChi2(); | |||
| 318 | fo->setParam(ilocal, save - eps); | |||
| 319 | double v2 = getChi2(); | |||
| 320 | double result = (v1 - v2) / (2 * eps); | |||
| 321 | fo->setParam(ilocal, save); | |||
| 322 | return result; | |||
| 323 | } | |||
| 324 | ||||
| 325 | double SoftGaussParticleConstraint::num2ndDerivative(int ifo1, int ilocal1, double eps1, | |||
| 326 | int ifo2, int ilocal2, double eps2) | |||
| 327 | { | |||
| 328 | double result; | |||
| 329 | ||||
| 330 | if (ifo1 == ifo2 && ilocal1 == ilocal2) { | |||
| 331 | ParticleFitObject* fo = fitobjects[ifo1]; | |||
| 332 | assert(fo)(static_cast <bool> (fo) ? void (0) : __assert_fail ("fo" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 332 , __extension__ __PRETTY_FUNCTION__)); | |||
| 333 | double save = fo->getParam(ilocal1); | |||
| 334 | double v0 = getChi2(); | |||
| 335 | fo->setParam(ilocal1, save + eps1); | |||
| 336 | double v1 = getChi2(); | |||
| 337 | fo->setParam(ilocal1, save - eps1); | |||
| 338 | double v2 = getChi2(); | |||
| 339 | result = (v1 + v2 - 2 * v0) / (eps1 * eps1); | |||
| 340 | fo->setParam(ilocal1, save); | |||
| 341 | } else { | |||
| 342 | ParticleFitObject* fo1 = fitobjects[ifo1]; | |||
| 343 | assert(fo1)(static_cast <bool> (fo1) ? void (0) : __assert_fail ("fo1" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 343 , __extension__ __PRETTY_FUNCTION__)); | |||
| 344 | ParticleFitObject* fo2 = fitobjects[ifo2]; | |||
| 345 | assert(fo2)(static_cast <bool> (fo2) ? void (0) : __assert_fail ("fo2" , "analysis/OrcaKinFit/src/SoftGaussParticleConstraint.cc", 345 , __extension__ __PRETTY_FUNCTION__)); | |||
| 346 | double save1 = fo1->getParam(ilocal1); | |||
| 347 | double save2 = fo2->getParam(ilocal2); | |||
| 348 | fo1->setParam(ilocal1, save1 + eps1); | |||
| 349 | fo2->setParam(ilocal2, save2 + eps2); | |||
| 350 | double v11 = getChi2(); | |||
| 351 | fo2->setParam(ilocal2, save2 - eps2); | |||
| 352 | double v12 = getChi2(); | |||
| 353 | fo1->setParam(ilocal1, save1 - eps1); | |||
| 354 | double v22 = getChi2(); | |||
| 355 | fo2->setParam(ilocal2, save2 + eps2); | |||
| 356 | double v21 = getChi2(); | |||
| 357 | result = (v11 + v22 - v12 - v21) / (4 * eps1 * eps2); | |||
| 358 | fo1->setParam(ilocal1, save1); | |||
| 359 | fo2->setParam(ilocal2, save2); | |||
| 360 | } | |||
| 361 | return result; | |||
| 362 | } | |||
| 363 | ||||
| 364 | int SoftGaussParticleConstraint::getVarBasis() | |||
| 365 | { | |||
| 366 | return VAR_BASIS; | |||
| 367 | } | |||
| 368 | }// end OrcaKinFit namespace | |||
| 369 | } // end Belle2 namespace | |||
| 370 |