| File: | genfit2/code2/trackReps/src/RKTrackRep.cc |
| Warning: | line 2002, column 7 Value stored to 'Way' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | /* Copyright 2008-2013, Technische Universitaet Muenchen, Ludwig-Maximilians-Universität München |
| 2 | Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch & Tobias Schlüter |
| 3 | |
| 4 | This file is part of GENFIT. |
| 5 | |
| 6 | GENFIT is free software: you can redistribute it and/or modify |
| 7 | it under the terms of the GNU Lesser General Public License as published |
| 8 | by the Free Software Foundation, either version 3 of the License, or |
| 9 | (at your option) any later version. |
| 10 | |
| 11 | GENFIT is distributed in the hope that it will be useful, |
| 12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 14 | GNU Lesser General Public License for more details. |
| 15 | |
| 16 | You should have received a copy of the GNU Lesser General Public License |
| 17 | along with GENFIT. If not, see <http://www.gnu.org/licenses/>. |
| 18 | */ |
| 19 | |
| 20 | #include "RKTrackRep.h" |
| 21 | #include "IO.h" |
| 22 | |
| 23 | #include <Exception.h> |
| 24 | #include <FieldManager.h> |
| 25 | #include <MaterialEffects.h> |
| 26 | #include <MeasuredStateOnPlane.h> |
| 27 | #include <MeasurementOnPlane.h> |
| 28 | |
| 29 | #include <TMatrixDSymEigen.h> |
| 30 | #include <TBuffer.h> |
| 31 | #include <TDecompLU.h> |
| 32 | #include <TMath.h> |
| 33 | |
| 34 | #include <algorithm> |
| 35 | |
| 36 | #define MINSTEP0.001 0.001 // minimum step [cm] for Runge Kutta and iteration to POCA |
| 37 | |
| 38 | namespace { |
| 39 | // Use fast inversion instead of LU decomposition? |
| 40 | const bool useInvertFast = false; |
| 41 | } |
| 42 | |
| 43 | namespace genfit { |
| 44 | |
| 45 | |
| 46 | RKTrackRep::RKTrackRep() : |
| 47 | AbsTrackRep(), |
| 48 | lastStartState_(this), |
| 49 | lastEndState_(this), |
| 50 | RKStepsFXStart_(0), |
| 51 | RKStepsFXStop_(0), |
| 52 | fJacobian_(5,5), |
| 53 | fNoise_(5), |
| 54 | useCache_(false), |
| 55 | cachePos_(0) |
| 56 | { |
| 57 | initArrays(); |
| 58 | } |
| 59 | |
| 60 | |
| 61 | RKTrackRep::RKTrackRep(int pdgCode, char propDir) : |
| 62 | AbsTrackRep(pdgCode, propDir), |
| 63 | lastStartState_(this), |
| 64 | lastEndState_(this), |
| 65 | RKStepsFXStart_(0), |
| 66 | RKStepsFXStop_(0), |
| 67 | fJacobian_(5,5), |
| 68 | fNoise_(5), |
| 69 | useCache_(false), |
| 70 | cachePos_(0) |
| 71 | { |
| 72 | initArrays(); |
| 73 | } |
| 74 | |
| 75 | |
| 76 | RKTrackRep::~RKTrackRep() { |
| 77 | ; |
| 78 | } |
| 79 | |
| 80 | |
| 81 | double RKTrackRep::extrapolateToPlane(StateOnPlane& state, |
| 82 | const SharedPlanePtr& plane, |
| 83 | bool stopAtBoundary, |
| 84 | bool calcJacobianNoise) const { |
| 85 | |
| 86 | if (debugLvl_ > 0) { |
| 87 | debugOut << "RKTrackRep::extrapolateToPlane()\n"; |
| 88 | } |
| 89 | |
| 90 | |
| 91 | if (state.getPlane() == plane) { |
| 92 | if (debugLvl_ > 0) { |
| 93 | debugOut << "state is already defined at plane. Do nothing! \n"; |
| 94 | } |
| 95 | return 0; |
| 96 | } |
| 97 | |
| 98 | checkCache(state, &plane); |
| 99 | |
| 100 | // to 7D |
| 101 | M1x7 state7 = {{0, 0, 0, 0, 0, 0, 0}}; |
| 102 | getState7(state, state7); |
| 103 | |
| 104 | TMatrixDSym* covPtr(nullptr); |
| 105 | bool fillExtrapSteps(false); |
| 106 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 107 | covPtr = &(static_cast<MeasuredStateOnPlane*>(&state)->getCov()); |
| 108 | fillExtrapSteps = true; |
| 109 | } |
| 110 | else if (calcJacobianNoise) |
| 111 | fillExtrapSteps = true; |
| 112 | |
| 113 | // actual extrapolation |
| 114 | bool isAtBoundary(false); |
| 115 | double flightTime( 0. ); |
| 116 | double coveredDistance( Extrap(*(state.getPlane()), *plane, getCharge(state), getMass(state), isAtBoundary, state7, flightTime, fillExtrapSteps, covPtr, false, stopAtBoundary) ); |
| 117 | |
| 118 | if (stopAtBoundary && isAtBoundary) { |
| 119 | state.setPlane(SharedPlanePtr(new DetPlane(TVector3(state7[0], state7[1], state7[2]), |
| 120 | TVector3(state7[3], state7[4], state7[5])))); |
| 121 | } |
| 122 | else { |
| 123 | state.setPlane(plane); |
| 124 | } |
| 125 | |
| 126 | // back to 5D |
| 127 | getState5(state, state7); |
| 128 | setTime(state, getTime(state) + flightTime); |
| 129 | |
| 130 | lastEndState_ = state; |
| 131 | |
| 132 | return coveredDistance; |
| 133 | } |
| 134 | |
| 135 | |
| 136 | double RKTrackRep::extrapolateToLine(StateOnPlane& state, |
| 137 | const TVector3& linePoint, |
| 138 | const TVector3& lineDirection, |
| 139 | bool stopAtBoundary, |
| 140 | bool calcJacobianNoise) const { |
| 141 | |
| 142 | if (debugLvl_ > 0) { |
| 143 | debugOut << "RKTrackRep::extrapolateToLine()\n"; |
| 144 | } |
| 145 | |
| 146 | checkCache(state, nullptr); |
| 147 | |
| 148 | static const unsigned int maxIt(1000); |
| 149 | |
| 150 | // to 7D |
| 151 | M1x7 state7; |
| 152 | getState7(state, state7); |
| 153 | |
| 154 | bool fillExtrapSteps(false); |
| 155 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 156 | fillExtrapSteps = true; |
| 157 | } |
| 158 | else if (calcJacobianNoise) |
| 159 | fillExtrapSteps = true; |
| 160 | |
| 161 | // cppcheck-suppress unreadVariable |
| 162 | double step(0.), lastStep(0.), maxStep(1.E99), angle(0), distToPoca(0), tracklength(0); |
| 163 | double charge = getCharge(state); |
| 164 | double mass = getMass(state); |
| 165 | double flightTime = 0; |
| 166 | TVector3 dir(state7[3], state7[4], state7[5]); |
| 167 | TVector3 lastDir(0,0,0); |
| 168 | TVector3 poca, poca_onwire; |
| 169 | bool isAtBoundary(false); |
| 170 | |
| 171 | DetPlane startPlane(*(state.getPlane())); |
| 172 | SharedPlanePtr plane(new DetPlane(linePoint, dir.Cross(lineDirection), lineDirection)); |
| 173 | unsigned int iterations(0); |
| 174 | |
| 175 | while(true){ |
| 176 | if(++iterations == maxIt) { |
| 177 | Exception exc("RKTrackRep::extrapolateToLine ==> extrapolation to line failed, maximum number of iterations reached",__LINE__177,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 178 | exc.setFatal(); |
| 179 | throw exc; |
| 180 | } |
| 181 | |
| 182 | lastStep = step; |
| 183 | lastDir = dir; |
| 184 | |
| 185 | step = this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, maxStep); |
| 186 | tracklength += step; |
| 187 | |
| 188 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 189 | poca.SetXYZ(state7[0], state7[1], state7[2]); |
| 190 | poca_onwire = pocaOnLine(linePoint, lineDirection, poca); |
| 191 | |
| 192 | // check break conditions |
| 193 | if (stopAtBoundary && isAtBoundary) { |
| 194 | plane->setON(dir, poca); |
| 195 | break; |
| 196 | } |
| 197 | |
| 198 | angle = fabs(dir.Angle((poca_onwire-poca))-TMath::PiOver2()); // angle between direction and connection to point - 90 deg |
| 199 | distToPoca = (poca_onwire-poca).Mag(); |
| 200 | if (angle*distToPoca < 0.1*MINSTEP0.001) break; |
| 201 | |
| 202 | // if lastStep and step have opposite sign, the real normal vector lies somewhere between the last two normal vectors (i.e. the directions) |
| 203 | // -> try mean value of the two (normalization not needed) |
| 204 | if (lastStep*step < 0){ |
| 205 | dir += lastDir; |
| 206 | maxStep = 0.5*fabs(lastStep); // make it converge! |
| 207 | } |
| 208 | |
| 209 | startPlane = *plane; |
| 210 | plane->setU(dir.Cross(lineDirection)); |
| 211 | } |
| 212 | |
| 213 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 214 | // make use of the cache |
| 215 | lastEndState_.setPlane(plane); |
| 216 | getState5(lastEndState_, state7); |
| 217 | |
| 218 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 219 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 220 | } |
| 221 | else { |
| 222 | state.setPlane(plane); |
| 223 | getState5(state, state7); |
| 224 | state.getAuxInfo()(1) += flightTime; |
| 225 | } |
| 226 | |
| 227 | if (debugLvl_ > 0) { |
| 228 | debugOut << "RKTrackRep::extrapolateToLine(): Reached POCA after " << iterations+1 << " iterations. Distance: " << (poca_onwire-poca).Mag() << " cm. Angle deviation: " << dir.Angle((poca_onwire-poca))-TMath::PiOver2() << " rad \n"; |
| 229 | } |
| 230 | |
| 231 | lastEndState_ = state; |
| 232 | |
| 233 | return tracklength; |
| 234 | } |
| 235 | |
| 236 | |
| 237 | double RKTrackRep::extrapToPoint(StateOnPlane& state, |
| 238 | const TVector3& point, |
| 239 | const TMatrixDSym* G, |
| 240 | bool stopAtBoundary, |
| 241 | bool calcJacobianNoise) const { |
| 242 | |
| 243 | if (debugLvl_ > 0) { |
| 244 | debugOut << "RKTrackRep::extrapolateToPoint()\n"; |
| 245 | } |
| 246 | |
| 247 | checkCache(state, nullptr); |
| 248 | |
| 249 | static const unsigned int maxIt(1000); |
| 250 | |
| 251 | // to 7D |
| 252 | M1x7 state7; |
| 253 | getState7(state, state7); |
| 254 | |
| 255 | bool fillExtrapSteps(false); |
| 256 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 257 | fillExtrapSteps = true; |
| 258 | } |
| 259 | else if (calcJacobianNoise) |
| 260 | fillExtrapSteps = true; |
| 261 | |
| 262 | // cppcheck-suppress unreadVariable |
| 263 | double step(0.), lastStep(0.), maxStep(1.E99), angle(0), distToPoca(0), tracklength(0); |
| 264 | TVector3 dir(state7[3], state7[4], state7[5]); |
| 265 | if (G != nullptr) { |
| 266 | if(G->GetNrows() != 3) { |
| 267 | Exception exc("RKTrackRep::extrapolateToLine ==> G is not 3x3",__LINE__267,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 268 | exc.setFatal(); |
| 269 | throw exc; |
| 270 | } |
| 271 | dir = TMatrix(*G) * dir; |
| 272 | } |
| 273 | TVector3 lastDir(0,0,0); |
| 274 | |
| 275 | TVector3 poca; |
| 276 | bool isAtBoundary(false); |
| 277 | |
| 278 | DetPlane startPlane(*(state.getPlane())); |
| 279 | SharedPlanePtr plane(new DetPlane(point, dir)); |
| 280 | unsigned int iterations(0); |
| 281 | double charge = getCharge(state); |
| 282 | double mass = getMass(state); |
| 283 | double flightTime = 0; |
| 284 | |
| 285 | while(true){ |
| 286 | if(++iterations == maxIt) { |
| 287 | Exception exc("RKTrackRep::extrapolateToPoint ==> extrapolation to point failed, maximum number of iterations reached",__LINE__287,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 288 | exc.setFatal(); |
| 289 | throw exc; |
| 290 | } |
| 291 | |
| 292 | lastStep = step; |
| 293 | lastDir = dir; |
| 294 | |
| 295 | step = this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, maxStep); |
| 296 | tracklength += step; |
| 297 | |
| 298 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 299 | if (G != nullptr) { |
| 300 | dir = TMatrix(*G) * dir; |
| 301 | } |
| 302 | poca.SetXYZ(state7[0], state7[1], state7[2]); |
| 303 | |
| 304 | // check break conditions |
| 305 | if (stopAtBoundary && isAtBoundary) { |
| 306 | plane->setON(dir, poca); |
| 307 | break; |
| 308 | } |
| 309 | |
| 310 | angle = fabs(dir.Angle((point-poca))-TMath::PiOver2()); // angle between direction and connection to point - 90 deg |
| 311 | distToPoca = (point-poca).Mag(); |
| 312 | if (angle*distToPoca < 0.1*MINSTEP0.001) break; |
| 313 | |
| 314 | // if lastStep and step have opposite sign, the real normal vector lies somewhere between the last two normal vectors (i.e. the directions) |
| 315 | // -> try mean value of the two |
| 316 | if (lastStep*step < 0){ |
| 317 | if (G != nullptr) { // after multiplication with G, dir has not length 1 anymore in general |
| 318 | dir.SetMag(1.); |
| 319 | lastDir.SetMag(1.); |
| 320 | } |
| 321 | dir += lastDir; |
| 322 | maxStep = 0.5*fabs(lastStep); // make it converge! |
| 323 | } |
| 324 | |
| 325 | startPlane = *plane; |
| 326 | plane->setNormal(dir); |
| 327 | } |
| 328 | |
| 329 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 330 | // make use of the cache |
| 331 | lastEndState_.setPlane(plane); |
| 332 | getState5(lastEndState_, state7); |
| 333 | |
| 334 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 335 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 336 | } |
| 337 | else { |
| 338 | state.setPlane(plane); |
| 339 | getState5(state, state7); |
| 340 | state.getAuxInfo()(1) += flightTime; |
| 341 | } |
| 342 | |
| 343 | |
| 344 | if (debugLvl_ > 0) { |
| 345 | debugOut << "RKTrackRep::extrapolateToPoint(): Reached POCA after " << iterations+1 << " iterations. Distance: " << (point-poca).Mag() << " cm. Angle deviation: " << dir.Angle((point-poca))-TMath::PiOver2() << " rad \n"; |
| 346 | } |
| 347 | |
| 348 | lastEndState_ = state; |
| 349 | |
| 350 | return tracklength; |
| 351 | } |
| 352 | |
| 353 | |
| 354 | double RKTrackRep::extrapolateToCylinder(StateOnPlane& state, |
| 355 | double radius, |
| 356 | const TVector3& linePoint, |
| 357 | const TVector3& lineDirection, |
| 358 | bool stopAtBoundary, |
| 359 | bool calcJacobianNoise) const { |
| 360 | |
| 361 | if (debugLvl_ > 0) { |
| 362 | debugOut << "RKTrackRep::extrapolateToCylinder()\n"; |
| 363 | } |
| 364 | |
| 365 | checkCache(state, nullptr); |
| 366 | |
| 367 | static const unsigned int maxIt(1000); |
| 368 | |
| 369 | // to 7D |
| 370 | M1x7 state7; |
| 371 | getState7(state, state7); |
| 372 | |
| 373 | bool fillExtrapSteps(false); |
| 374 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 375 | fillExtrapSteps = true; |
| 376 | } |
| 377 | else if (calcJacobianNoise) |
| 378 | fillExtrapSteps = true; |
| 379 | |
| 380 | double tracklength(0.), maxStep(1.E99); |
| 381 | |
| 382 | TVector3 dest, pos, dir; |
| 383 | |
| 384 | bool isAtBoundary(false); |
| 385 | |
| 386 | DetPlane startPlane(*(state.getPlane())); |
| 387 | SharedPlanePtr plane(new DetPlane()); |
| 388 | unsigned int iterations(0); |
| 389 | double charge = getCharge(state); |
| 390 | double mass = getMass(state); |
| 391 | double flightTime = 0; |
| 392 | |
| 393 | while(true){ |
| 394 | if(++iterations == maxIt) { |
| 395 | Exception exc("RKTrackRep::extrapolateToCylinder ==> maximum number of iterations reached",__LINE__395,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 396 | exc.setFatal(); |
| 397 | throw exc; |
| 398 | } |
| 399 | |
| 400 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 401 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 402 | |
| 403 | // solve quadratic equation |
| 404 | TVector3 AO = (pos - linePoint); |
| 405 | TVector3 AOxAB = (AO.Cross(lineDirection)); |
| 406 | TVector3 VxAB = (dir.Cross(lineDirection)); |
| 407 | float ab2 = (lineDirection * lineDirection); |
| 408 | float a = (VxAB * VxAB); |
| 409 | float b = 2 * (VxAB * AOxAB); |
| 410 | float c = (AOxAB * AOxAB) - (radius*radius * ab2); |
| 411 | double arg = b*b - 4.*a*c; |
| 412 | if(arg < 0) { |
| 413 | Exception exc("RKTrackRep::extrapolateToCylinder ==> cannot solve",__LINE__413,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 414 | exc.setFatal(); |
| 415 | throw exc; |
| 416 | } |
| 417 | double term = sqrt(arg); |
| 418 | double k1, k2; |
| 419 | if (b<0) { |
| 420 | k1 = (-b + term)/(2.*a); |
| 421 | k2 = 2.*c/(-b + term); |
| 422 | } |
| 423 | else { |
| 424 | k1 = 2.*c/(-b - term); |
| 425 | k2 = (-b - term)/(2.*a); |
| 426 | } |
| 427 | |
| 428 | // select smallest absolute solution -> closest cylinder surface |
| 429 | double k = k1; |
| 430 | if (fabs(k2)<fabs(k)) |
| 431 | k = k2; |
| 432 | |
| 433 | if (debugLvl_ > 0) { |
| 434 | debugOut << "RKTrackRep::extrapolateToCylinder(); k = " << k << "\n"; |
| 435 | } |
| 436 | |
| 437 | dest = pos + k * dir; |
| 438 | |
| 439 | plane->setO(dest); |
| 440 | plane->setUV((dest-linePoint).Cross(lineDirection), lineDirection); |
| 441 | |
| 442 | tracklength += this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, maxStep); |
| 443 | |
| 444 | // check break conditions |
| 445 | if (stopAtBoundary && isAtBoundary) { |
| 446 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 447 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 448 | plane->setO(pos); |
| 449 | plane->setUV((pos-linePoint).Cross(lineDirection), lineDirection); |
| 450 | break; |
| 451 | } |
| 452 | |
| 453 | if(fabs(k)<MINSTEP0.001) break; |
| 454 | |
| 455 | startPlane = *plane; |
| 456 | |
| 457 | } |
| 458 | |
| 459 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 460 | // make use of the cache |
| 461 | lastEndState_.setPlane(plane); |
| 462 | getState5(lastEndState_, state7); |
| 463 | |
| 464 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 465 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 466 | } |
| 467 | else { |
| 468 | state.setPlane(plane); |
| 469 | getState5(state, state7); |
| 470 | state.getAuxInfo()(1) += flightTime; |
| 471 | } |
| 472 | |
| 473 | lastEndState_ = state; |
| 474 | |
| 475 | return tracklength; |
| 476 | } |
| 477 | |
| 478 | |
| 479 | double RKTrackRep::extrapolateToCone(StateOnPlane& state, |
| 480 | double openingAngle, |
| 481 | const TVector3& conePoint, |
| 482 | const TVector3& coneDirection, |
| 483 | bool stopAtBoundary, |
| 484 | bool calcJacobianNoise) const { |
| 485 | |
| 486 | if (debugLvl_ > 0) { |
| 487 | debugOut << "RKTrackRep::extrapolateToCone()\n"; |
| 488 | } |
| 489 | |
| 490 | checkCache(state, nullptr); |
| 491 | |
| 492 | static const unsigned int maxIt(1000); |
| 493 | |
| 494 | // to 7D |
| 495 | M1x7 state7; |
| 496 | getState7(state, state7); |
| 497 | |
| 498 | bool fillExtrapSteps(false); |
| 499 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 500 | fillExtrapSteps = true; |
| 501 | } |
| 502 | else if (calcJacobianNoise) |
| 503 | fillExtrapSteps = true; |
| 504 | |
| 505 | double tracklength(0.), maxStep(1.E99); |
| 506 | |
| 507 | TVector3 dest, pos, dir; |
| 508 | |
| 509 | bool isAtBoundary(false); |
| 510 | |
| 511 | DetPlane startPlane(*(state.getPlane())); |
| 512 | SharedPlanePtr plane(new DetPlane()); |
| 513 | unsigned int iterations(0); |
| 514 | double charge = getCharge(state); |
| 515 | double mass = getMass(state); |
| 516 | double flightTime = 0; |
| 517 | |
| 518 | while(true){ |
| 519 | if(++iterations == maxIt) { |
| 520 | Exception exc("RKTrackRep::extrapolateToCone ==> maximum number of iterations reached",__LINE__520,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 521 | exc.setFatal(); |
| 522 | throw exc; |
| 523 | } |
| 524 | |
| 525 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 526 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 527 | |
| 528 | // solve quadratic equation a k^2 + 2 b k + c = 0 |
| 529 | // a = (U . D)^2 - cos^2 alpha * U^2 |
| 530 | // b = (Delta . D) * (U . D) - cos^2 alpha * (U . Delta) |
| 531 | // c = (Delta . D)^2 - cos^2 alpha * Delta^2 |
| 532 | // Delta = P - V, P track point, U track direction, V cone point, D cone direction, alpha opening angle of cone |
| 533 | TVector3 cDirection = coneDirection.Unit(); |
| 534 | TVector3 Delta = (pos - conePoint); |
| 535 | double DirDelta = cDirection * Delta; |
| 536 | double Delta2 = Delta*Delta; |
| 537 | double UDir = dir * cDirection; |
| 538 | double UDelta = dir * Delta; |
| 539 | double U2 = dir * dir; |
| 540 | double cosAngle2 = cos(openingAngle)*cos(openingAngle); |
| 541 | double a = UDir*UDir - cosAngle2*U2; |
| 542 | double b = UDir*DirDelta - cosAngle2*UDelta; |
| 543 | double c = DirDelta*DirDelta - cosAngle2*Delta2; |
| 544 | |
| 545 | double arg = b*b - a*c; |
| 546 | if(arg < -1e-9) { |
| 547 | Exception exc("RKTrackRep::extrapolateToCone ==> cannot solve",__LINE__547,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 548 | exc.setFatal(); |
| 549 | throw exc; |
| 550 | } else if(arg < 0) { |
| 551 | arg = 0; |
| 552 | } |
| 553 | |
| 554 | double term = sqrt(arg); |
| 555 | double k1, k2; |
| 556 | k1 = (-b + term) / a; |
| 557 | k2 = (-b - term) / a; |
| 558 | |
| 559 | // select smallest absolute solution -> closest cone surface |
| 560 | double k = k1; |
| 561 | if(fabs(k2) < fabs(k)) { |
| 562 | k = k2; |
| 563 | } |
| 564 | |
| 565 | if (debugLvl_ > 0) { |
| 566 | debugOut << "RKTrackRep::extrapolateToCone(); k = " << k << "\n"; |
| 567 | } |
| 568 | |
| 569 | dest = pos + k * dir; |
| 570 | // debugOut << "In cone extrapolation "; |
| 571 | // dest.Print(); |
| 572 | |
| 573 | plane->setO(dest); |
| 574 | plane->setUV((dest-conePoint).Cross(coneDirection), dest-conePoint); |
| 575 | |
| 576 | tracklength += this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, maxStep); |
| 577 | |
| 578 | // check break conditions |
| 579 | if (stopAtBoundary && isAtBoundary) { |
| 580 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 581 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 582 | plane->setO(pos); |
| 583 | plane->setUV((pos-conePoint).Cross(coneDirection), pos-conePoint); |
| 584 | break; |
| 585 | } |
| 586 | |
| 587 | if(fabs(k)<MINSTEP0.001) break; |
| 588 | |
| 589 | startPlane = *plane; |
| 590 | |
| 591 | } |
| 592 | |
| 593 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 594 | // make use of the cache |
| 595 | lastEndState_.setPlane(plane); |
| 596 | getState5(lastEndState_, state7); |
| 597 | |
| 598 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 599 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 600 | } |
| 601 | else { |
| 602 | state.setPlane(plane); |
| 603 | getState5(state, state7); |
| 604 | state.getAuxInfo()(1) += flightTime; |
| 605 | } |
| 606 | |
| 607 | lastEndState_ = state; |
| 608 | |
| 609 | return tracklength; |
| 610 | } |
| 611 | |
| 612 | |
| 613 | double RKTrackRep::extrapolateToSphere(StateOnPlane& state, |
| 614 | double radius, |
| 615 | const TVector3& point, // center |
| 616 | bool stopAtBoundary, |
| 617 | bool calcJacobianNoise) const { |
| 618 | |
| 619 | if (debugLvl_ > 0) { |
| 620 | debugOut << "RKTrackRep::extrapolateToSphere()\n"; |
| 621 | } |
| 622 | |
| 623 | checkCache(state, nullptr); |
| 624 | |
| 625 | static const unsigned int maxIt(1000); |
| 626 | |
| 627 | // to 7D |
| 628 | M1x7 state7; |
| 629 | getState7(state, state7); |
| 630 | |
| 631 | bool fillExtrapSteps(false); |
| 632 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 633 | fillExtrapSteps = true; |
| 634 | } |
| 635 | else if (calcJacobianNoise) |
| 636 | fillExtrapSteps = true; |
| 637 | |
| 638 | double tracklength(0.), maxStep(1.E99); |
| 639 | |
| 640 | TVector3 dest, pos, dir; |
| 641 | |
| 642 | bool isAtBoundary(false); |
| 643 | |
| 644 | DetPlane startPlane(*(state.getPlane())); |
| 645 | SharedPlanePtr plane(new DetPlane()); |
| 646 | unsigned int iterations(0); |
| 647 | double charge = getCharge(state); |
| 648 | double mass = getMass(state); |
| 649 | double flightTime = 0; |
| 650 | |
| 651 | while(true){ |
| 652 | if(++iterations == maxIt) { |
| 653 | Exception exc("RKTrackRep::extrapolateToSphere ==> maximum number of iterations reached",__LINE__653,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 654 | exc.setFatal(); |
| 655 | throw exc; |
| 656 | } |
| 657 | |
| 658 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 659 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 660 | |
| 661 | // solve quadratic equation |
| 662 | TVector3 AO = (pos - point); |
| 663 | double dirAO = dir * AO; |
| 664 | double arg = dirAO*dirAO - AO*AO + radius*radius; |
| 665 | if(arg < 0) { |
| 666 | Exception exc("RKTrackRep::extrapolateToSphere ==> cannot solve",__LINE__666,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 667 | exc.setFatal(); |
| 668 | throw exc; |
| 669 | } |
| 670 | double term = sqrt(arg); |
| 671 | double k1, k2; |
| 672 | k1 = -dirAO + term; |
| 673 | k2 = -dirAO - term; |
| 674 | |
| 675 | // select smallest absolute solution -> closest cylinder surface |
| 676 | double k = k1; |
| 677 | if (fabs(k2)<fabs(k)) |
| 678 | k = k2; |
| 679 | |
| 680 | if (debugLvl_ > 0) { |
| 681 | debugOut << "RKTrackRep::extrapolateToSphere(); k = " << k << "\n"; |
| 682 | } |
| 683 | |
| 684 | dest = pos + k * dir; |
| 685 | |
| 686 | plane->setON(dest, dest-point); |
| 687 | |
| 688 | tracklength += this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, maxStep); |
| 689 | |
| 690 | // check break conditions |
| 691 | if (stopAtBoundary && isAtBoundary) { |
| 692 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 693 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 694 | plane->setON(pos, pos-point); |
| 695 | break; |
| 696 | } |
| 697 | |
| 698 | if(fabs(k)<MINSTEP0.001) break; |
| 699 | |
| 700 | startPlane = *plane; |
| 701 | |
| 702 | } |
| 703 | |
| 704 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 705 | // make use of the cache |
| 706 | lastEndState_.setPlane(plane); |
| 707 | getState5(lastEndState_, state7); |
| 708 | |
| 709 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 710 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 711 | } |
| 712 | else { |
| 713 | state.setPlane(plane); |
| 714 | getState5(state, state7); |
| 715 | state.getAuxInfo()(1) += flightTime; |
| 716 | } |
| 717 | |
| 718 | lastEndState_ = state; |
| 719 | |
| 720 | return tracklength; |
| 721 | } |
| 722 | |
| 723 | |
| 724 | double RKTrackRep::extrapolateBy(StateOnPlane& state, |
| 725 | double step, |
| 726 | bool stopAtBoundary, |
| 727 | bool calcJacobianNoise) const { |
| 728 | |
| 729 | if (debugLvl_ > 0) { |
| 730 | debugOut << "RKTrackRep::extrapolateBy()\n"; |
| 731 | } |
| 732 | |
| 733 | checkCache(state, nullptr); |
| 734 | |
| 735 | static const unsigned int maxIt(1000); |
| 736 | |
| 737 | // to 7D |
| 738 | M1x7 state7; |
| 739 | getState7(state, state7); |
| 740 | |
| 741 | bool fillExtrapSteps(false); |
| 742 | if (dynamic_cast<MeasuredStateOnPlane*>(&state) != nullptr) { |
| 743 | fillExtrapSteps = true; |
| 744 | } |
| 745 | else if (calcJacobianNoise) |
| 746 | fillExtrapSteps = true; |
| 747 | |
| 748 | double tracklength(0.); |
| 749 | |
| 750 | TVector3 dest, pos, dir; |
| 751 | |
| 752 | bool isAtBoundary(false); |
| 753 | |
| 754 | DetPlane startPlane(*(state.getPlane())); |
| 755 | SharedPlanePtr plane(new DetPlane()); |
| 756 | unsigned int iterations(0); |
| 757 | double charge = getCharge(state); |
| 758 | double mass = getMass(state); |
| 759 | double flightTime = 0; |
| 760 | |
| 761 | while(true){ |
| 762 | if(++iterations == maxIt) { |
| 763 | Exception exc("RKTrackRep::extrapolateBy ==> maximum number of iterations reached",__LINE__763,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 764 | exc.setFatal(); |
| 765 | throw exc; |
| 766 | } |
| 767 | |
| 768 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 769 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 770 | |
| 771 | dest = pos + 1.5*(step-tracklength) * dir; |
| 772 | |
| 773 | plane->setON(dest, dir); |
| 774 | |
| 775 | tracklength += this->Extrap(startPlane, *plane, charge, mass, isAtBoundary, state7, flightTime, false, nullptr, true, stopAtBoundary, (step-tracklength)); |
| 776 | |
| 777 | // check break conditions |
| 778 | if (stopAtBoundary && isAtBoundary) { |
| 779 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 780 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 781 | plane->setON(pos, dir); |
| 782 | break; |
| 783 | } |
| 784 | |
| 785 | if (fabs(tracklength-step) < MINSTEP0.001) { |
| 786 | if (debugLvl_ > 0) { |
| 787 | debugOut << "RKTrackRep::extrapolateBy(): reached after " << iterations << " iterations. \n"; |
| 788 | } |
| 789 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 790 | dir.SetXYZ(state7[3], state7[4], state7[5]); |
| 791 | plane->setON(pos, dir); |
| 792 | break; |
| 793 | } |
| 794 | |
| 795 | startPlane = *plane; |
| 796 | |
| 797 | } |
| 798 | |
| 799 | if (fillExtrapSteps) { // now do the full extrapolation with covariance matrix |
| 800 | // make use of the cache |
| 801 | lastEndState_.setPlane(plane); |
| 802 | getState5(lastEndState_, state7); |
| 803 | |
| 804 | tracklength = extrapolateToPlane(state, plane, false, true); |
| 805 | lastEndState_.getAuxInfo()(1) = state.getAuxInfo()(1); // Flight time |
| 806 | } |
| 807 | else { |
| 808 | state.setPlane(plane); |
| 809 | getState5(state, state7); |
| 810 | state.getAuxInfo()(1) += flightTime; |
| 811 | } |
| 812 | |
| 813 | lastEndState_ = state; |
| 814 | |
| 815 | return tracklength; |
| 816 | } |
| 817 | |
| 818 | |
| 819 | TVector3 RKTrackRep::getPos(const StateOnPlane& state) const { |
| 820 | M1x7 state7; |
| 821 | getState7(state, state7); |
| 822 | |
| 823 | return TVector3(state7[0], state7[1], state7[2]); |
| 824 | } |
| 825 | |
| 826 | |
| 827 | TVector3 RKTrackRep::getMom(const StateOnPlane& state) const { |
| 828 | M1x7 state7 = {{0, 0, 0, 0, 0, 0, 0}}; |
| 829 | getState7(state, state7); |
| 830 | |
| 831 | TVector3 mom(state7[3], state7[4], state7[5]); |
| 832 | mom.SetMag(getCharge(state)/state7[6]); |
| 833 | return mom; |
| 834 | } |
| 835 | |
| 836 | |
| 837 | void RKTrackRep::getPosMom(const StateOnPlane& state, TVector3& pos, TVector3& mom) const { |
| 838 | M1x7 state7 = {{0, 0, 0, 0, 0, 0, 0}}; |
| 839 | getState7(state, state7); |
| 840 | |
| 841 | pos.SetXYZ(state7[0], state7[1], state7[2]); |
| 842 | mom.SetXYZ(state7[3], state7[4], state7[5]); |
| 843 | mom.SetMag(getCharge(state)/state7[6]); |
| 844 | } |
| 845 | |
| 846 | |
| 847 | void RKTrackRep::getPosMomCov(const MeasuredStateOnPlane& state, TVector3& pos, TVector3& mom, TMatrixDSym& cov) const { |
| 848 | getPosMom(state, pos, mom); |
| 849 | cov.ResizeTo(6,6); |
| 850 | transformPM6(state, *((M6x6*) cov.GetMatrixArray())); |
| 851 | } |
| 852 | |
| 853 | |
| 854 | TMatrixDSym RKTrackRep::get6DCov(const MeasuredStateOnPlane& state) const { |
| 855 | TMatrixDSym cov(6); |
| 856 | transformPM6(state, *((M6x6*) cov.GetMatrixArray())); |
| 857 | |
| 858 | return cov; |
| 859 | } |
| 860 | |
| 861 | |
| 862 | double RKTrackRep::getCharge(const StateOnPlane& state) const { |
| 863 | |
| 864 | if (dynamic_cast<const MeasurementOnPlane*>(&state) != nullptr) { |
| 865 | Exception exc("RKTrackRep::getCharge - cannot get charge from MeasurementOnPlane",__LINE__865,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 866 | exc.setFatal(); |
| 867 | throw exc; |
| 868 | } |
| 869 | |
| 870 | double pdgCharge( this->getPDGCharge() ); |
| 871 | |
| 872 | // return pdgCharge with sign of q/p |
| 873 | if (state.getState()(0) * pdgCharge < 0) |
| 874 | return -pdgCharge; |
| 875 | else |
| 876 | return pdgCharge; |
| 877 | } |
| 878 | |
| 879 | |
| 880 | double RKTrackRep::getMomMag(const StateOnPlane& state) const { |
| 881 | // p = q / qop |
| 882 | double p = getCharge(state)/state.getState()(0); |
| 883 | assert (p>=0)(static_cast <bool> (p>=0) ? void (0) : __assert_fail ("p>=0", "genfit2/code2/trackReps/src/RKTrackRep.cc", 883 , __extension__ __PRETTY_FUNCTION__)); |
| 884 | return p; |
| 885 | } |
| 886 | |
| 887 | |
| 888 | double RKTrackRep::getMomVar(const MeasuredStateOnPlane& state) const { |
| 889 | |
| 890 | if (dynamic_cast<const MeasurementOnPlane*>(&state) != nullptr) { |
| 891 | Exception exc("RKTrackRep::getMomVar - cannot get momVar from MeasurementOnPlane",__LINE__891,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 892 | exc.setFatal(); |
| 893 | throw exc; |
| 894 | } |
| 895 | |
| 896 | // p(qop) = q/qop |
| 897 | // dp/d(qop) = - q / (qop^2) |
| 898 | // (delta p) = (delta qop) * |dp/d(qop)| = delta qop * |q / (qop^2)| |
| 899 | // (var p) = (var qop) * q^2 / (qop^4) |
| 900 | |
| 901 | // delta means sigma |
| 902 | // cov(0,0) is sigma^2 |
| 903 | |
| 904 | return state.getCov()(0,0) * pow(getCharge(state), 2) / pow(state.getState()(0), 4); |
| 905 | } |
| 906 | |
| 907 | |
| 908 | double RKTrackRep::getSpu(const StateOnPlane& state) const { |
| 909 | |
| 910 | if (dynamic_cast<const MeasurementOnPlane*>(&state) != nullptr) { |
| 911 | Exception exc("RKTrackRep::getSpu - cannot get spu from MeasurementOnPlane",__LINE__911,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 912 | exc.setFatal(); |
| 913 | throw exc; |
| 914 | } |
| 915 | |
| 916 | const TVectorD& auxInfo = state.getAuxInfo(); |
| 917 | if (auxInfo.GetNrows() == 2 |
| 918 | || auxInfo.GetNrows() == 1) // backwards compatibility with old RKTrackRep |
| 919 | return state.getAuxInfo()(0); |
| 920 | else |
| 921 | return 1.; |
| 922 | } |
| 923 | |
| 924 | double RKTrackRep::getTime(const StateOnPlane& state) const { |
| 925 | |
| 926 | const TVectorD& auxInfo = state.getAuxInfo(); |
| 927 | if (auxInfo.GetNrows() == 2) |
| 928 | return state.getAuxInfo()(1); |
| 929 | else |
| 930 | return 0.; |
| 931 | } |
| 932 | |
| 933 | |
| 934 | void RKTrackRep::calcForwardJacobianAndNoise(const M1x7& startState7, const DetPlane& startPlane, |
| 935 | const M1x7& destState7, const DetPlane& destPlane) const { |
| 936 | |
| 937 | if (debugLvl_ > 0) { |
| 938 | debugOut << "RKTrackRep::calcForwardJacobianAndNoise " << std::endl; |
| 939 | } |
| 940 | |
| 941 | if (ExtrapSteps_.size() == 0) { |
| 942 | Exception exc("RKTrackRep::calcForwardJacobianAndNoise ==> cache is empty. Extrapolation must run with a MeasuredStateOnPlane.",__LINE__942,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 943 | throw exc; |
| 944 | } |
| 945 | |
| 946 | // The Jacobians returned from RKutta are transposed. |
| 947 | TMatrixD jac(TMatrixD::kTransposed, TMatrixD(7, 7, ExtrapSteps_.back().jac7_.begin())); |
| 948 | TMatrixDSym noise(7, ExtrapSteps_.back().noise7_.begin()); |
| 949 | for (int i = ExtrapSteps_.size() - 2; i >= 0; --i) { |
| 950 | noise += TMatrixDSym(7, ExtrapSteps_[i].noise7_.begin()).Similarity(jac); |
| 951 | jac *= TMatrixD(TMatrixD::kTransposed, TMatrixD(7, 7, ExtrapSteps_[i].jac7_.begin())); |
| 952 | } |
| 953 | |
| 954 | // Project into 5x5 space. |
| 955 | M1x3 pTilde = {{startState7[3], startState7[4], startState7[5]}}; |
| 956 | const TVector3& normal = startPlane.getNormal(); |
| 957 | double pTildeW = pTilde[0] * normal.X() + pTilde[1] * normal.Y() + pTilde[2] * normal.Z(); |
| 958 | double spu = pTildeW > 0 ? 1 : -1; |
| 959 | for (unsigned int i=0; i<3; ++i) { |
| 960 | pTilde[i] *= spu/pTildeW; // | pTilde * W | has to be 1 (definition of pTilde) |
| 961 | } |
| 962 | M5x7 J_pM; |
| 963 | calcJ_pM_5x7(J_pM, startPlane.getU(), startPlane.getV(), pTilde, spu); |
| 964 | M7x5 J_Mp; |
| 965 | calcJ_Mp_7x5(J_Mp, destPlane.getU(), destPlane.getV(), destPlane.getNormal(), *((M1x3*) &destState7[3])); |
| 966 | jac.Transpose(jac); // Because the helper function wants transposed input. |
| 967 | RKTools::J_pMTTxJ_MMTTxJ_MpTT(J_Mp, *(M7x7 *)jac.GetMatrixArray(), |
| 968 | J_pM, *(M5x5 *)fJacobian_.GetMatrixArray()); |
| 969 | RKTools::J_MpTxcov7xJ_Mp(J_Mp, *(M7x7 *)noise.GetMatrixArray(), |
| 970 | *(M5x5 *)fNoise_.GetMatrixArray()); |
| 971 | |
| 972 | if (debugLvl_ > 0) { |
| 973 | debugOut << "total jacobian : "; fJacobian_.Print(); |
| 974 | debugOut << "total noise : "; fNoise_.Print(); |
| 975 | } |
| 976 | |
| 977 | } |
| 978 | |
| 979 | void RKTrackRep::getForwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise) const { |
| 980 | |
| 981 | jacobian.ResizeTo(5,5); |
| 982 | jacobian = fJacobian_; |
| 983 | |
| 984 | noise.ResizeTo(5,5); |
| 985 | noise = fNoise_; |
| 986 | |
| 987 | } |
| 988 | |
| 989 | void RKTrackRep::getForwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const { |
| 990 | |
| 991 | jacobian.ResizeTo(5,5); |
| 992 | jacobian = fJacobian_; |
| 993 | |
| 994 | noise.ResizeTo(5,5); |
| 995 | noise = fNoise_; |
| 996 | |
| 997 | // lastEndState_ = jacobian * lastStartState_ + deltaState |
| 998 | deltaState.ResizeTo(5); |
| 999 | // Calculate this without temporaries: |
| 1000 | //deltaState = lastEndState_.getState() - jacobian * lastStartState_.getState() |
| 1001 | deltaState = lastStartState_.getState(); |
| 1002 | deltaState *= jacobian; |
| 1003 | deltaState -= lastEndState_.getState(); |
| 1004 | deltaState *= -1; |
| 1005 | |
| 1006 | |
| 1007 | if (debugLvl_ > 0) { |
| 1008 | debugOut << "delta state : "; deltaState.Print(); |
| 1009 | } |
| 1010 | } |
| 1011 | |
| 1012 | |
| 1013 | void RKTrackRep::getBackwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const { |
| 1014 | |
| 1015 | if (debugLvl_ > 0) { |
| 1016 | debugOut << "RKTrackRep::getBackwardJacobianAndNoise " << std::endl; |
| 1017 | } |
| 1018 | |
| 1019 | if (ExtrapSteps_.size() == 0) { |
| 1020 | Exception exc("RKTrackRep::getBackwardJacobianAndNoise ==> cache is empty. Extrapolation must run with a MeasuredStateOnPlane.",__LINE__1020,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1021 | throw exc; |
| 1022 | } |
| 1023 | |
| 1024 | jacobian.ResizeTo(5,5); |
| 1025 | jacobian = fJacobian_; |
| 1026 | if (!useInvertFast) { |
| 1027 | bool status = TDecompLU::InvertLU(jacobian, 0.0); |
| 1028 | if(status == 0){ |
| 1029 | Exception e("cannot invert matrix, status = 0", __LINE__1029,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1030 | e.setFatal(); |
| 1031 | throw e; |
| 1032 | } |
| 1033 | } else { |
| 1034 | double det; |
| 1035 | jacobian.InvertFast(&det); |
| 1036 | if(det < 1e-80){ |
| 1037 | Exception e("cannot invert matrix, status = 0", __LINE__1037,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1038 | e.setFatal(); |
| 1039 | throw e; |
| 1040 | } |
| 1041 | } |
| 1042 | |
| 1043 | noise.ResizeTo(5,5); |
| 1044 | noise = fNoise_; |
| 1045 | noise.Similarity(jacobian); |
| 1046 | |
| 1047 | // lastStartState_ = jacobian * lastEndState_ + deltaState |
| 1048 | deltaState.ResizeTo(5); |
| 1049 | deltaState = lastStartState_.getState() - jacobian * lastEndState_.getState(); |
| 1050 | } |
| 1051 | |
| 1052 | |
| 1053 | std::vector<genfit::MatStep> RKTrackRep::getSteps() const { |
| 1054 | |
| 1055 | // Todo: test |
| 1056 | |
| 1057 | if (RKSteps_.size() == 0) { |
| 1058 | Exception exc("RKTrackRep::getSteps ==> cache is empty.",__LINE__1058,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1059 | throw exc; |
| 1060 | } |
| 1061 | |
| 1062 | std::vector<MatStep> retVal; |
| 1063 | retVal.reserve(RKSteps_.size()); |
| 1064 | |
| 1065 | for (unsigned int i = 0; i<RKSteps_.size(); ++i) { |
| 1066 | retVal.push_back(RKSteps_[i].matStep_); |
| 1067 | } |
| 1068 | |
| 1069 | return retVal; |
| 1070 | } |
| 1071 | |
| 1072 | |
| 1073 | double RKTrackRep::getRadiationLenght() const { |
| 1074 | |
| 1075 | // Todo: test |
| 1076 | |
| 1077 | if (RKSteps_.size() == 0) { |
| 1078 | Exception exc("RKTrackRep::getRadiationLenght ==> cache is empty.",__LINE__1078,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1079 | throw exc; |
| 1080 | } |
| 1081 | |
| 1082 | double radLen(0); |
| 1083 | |
| 1084 | for (unsigned int i = 0; i<RKSteps_.size(); ++i) { |
| 1085 | radLen += RKSteps_.at(i).matStep_.stepSize_ / RKSteps_.at(i).matStep_.material_.radiationLength; |
| 1086 | } |
| 1087 | |
| 1088 | return radLen; |
| 1089 | } |
| 1090 | |
| 1091 | |
| 1092 | |
| 1093 | void RKTrackRep::setPosMom(StateOnPlane& state, const TVector3& pos, const TVector3& mom) const { |
| 1094 | |
| 1095 | if (state.getRep() != this){ |
| 1096 | Exception exc("RKTrackRep::setPosMom ==> state is defined wrt. another TrackRep",__LINE__1096,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1097 | throw exc; |
| 1098 | } |
| 1099 | |
| 1100 | if (dynamic_cast<MeasurementOnPlane*>(&state) != nullptr) { |
| 1101 | Exception exc("RKTrackRep::setPosMom - cannot set pos/mom of a MeasurementOnPlane",__LINE__1101,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1102 | exc.setFatal(); |
| 1103 | throw exc; |
| 1104 | } |
| 1105 | |
| 1106 | if (mom.Mag2() == 0) { |
| 1107 | Exception exc("RKTrackRep::setPosMom - momentum is 0",__LINE__1107,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1108 | exc.setFatal(); |
| 1109 | throw exc; |
| 1110 | } |
| 1111 | |
| 1112 | // init auxInfo if that has not yet happened |
| 1113 | TVectorD& auxInfo = state.getAuxInfo(); |
| 1114 | if (auxInfo.GetNrows() != 2) { |
| 1115 | bool alreadySet = auxInfo.GetNrows() == 1; // backwards compatibility: don't overwrite old setting |
| 1116 | auxInfo.ResizeTo(2); |
| 1117 | if (!alreadySet) |
| 1118 | setSpu(state, 1.); |
| 1119 | } |
| 1120 | |
| 1121 | if (state.getPlane() != nullptr && state.getPlane()->distance(pos) < MINSTEP0.001) { // pos is on plane -> do not change plane! |
| 1122 | |
| 1123 | M1x7 state7; |
| 1124 | |
| 1125 | state7[0] = pos.X(); |
| 1126 | state7[1] = pos.Y(); |
| 1127 | state7[2] = pos.Z(); |
| 1128 | |
| 1129 | state7[3] = mom.X(); |
| 1130 | state7[4] = mom.Y(); |
| 1131 | state7[5] = mom.Z(); |
| 1132 | |
| 1133 | // normalize dir |
| 1134 | double norm = 1. / sqrt(state7[3]*state7[3] + state7[4]*state7[4] + state7[5]*state7[5]); |
| 1135 | for (unsigned int i=3; i<6; ++i) |
| 1136 | state7[i] *= norm; |
| 1137 | |
| 1138 | state7[6] = getCharge(state) * norm; |
| 1139 | |
| 1140 | getState5(state, state7); |
| 1141 | |
| 1142 | } |
| 1143 | else { // pos is not on plane -> create new plane! |
| 1144 | |
| 1145 | // TODO: Raise Warning that a new plane has been created! |
| 1146 | SharedPlanePtr plane(new DetPlane(pos, mom)); |
| 1147 | state.setPlane(plane); |
| 1148 | |
| 1149 | TVectorD& state5(state.getState()); |
| 1150 | |
| 1151 | state5(0) = getCharge(state)/mom.Mag(); // q/p |
| 1152 | state5(1) = 0.; // u' |
| 1153 | state5(2) = 0.; // v' |
| 1154 | state5(3) = 0.; // u |
| 1155 | state5(4) = 0.; // v |
| 1156 | |
| 1157 | setSpu(state, 1.); |
| 1158 | } |
| 1159 | |
| 1160 | } |
| 1161 | |
| 1162 | |
| 1163 | void RKTrackRep::setPosMom(StateOnPlane& state, const TVectorD& state6) const { |
| 1164 | if (state6.GetNrows()!=6){ |
| 1165 | Exception exc("RKTrackRep::setPosMom ==> state has to be 6d (x, y, z, px, py, pz)",__LINE__1165,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1166 | throw exc; |
| 1167 | } |
| 1168 | setPosMom(state, TVector3(state6(0), state6(1), state6(2)), TVector3(state6(3), state6(4), state6(5))); |
| 1169 | } |
| 1170 | |
| 1171 | |
| 1172 | void RKTrackRep::setPosMomErr(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TVector3& posErr, const TVector3& momErr) const { |
| 1173 | |
| 1174 | // TODO: test! |
| 1175 | |
| 1176 | setPosMom(state, pos, mom); |
| 1177 | |
| 1178 | const TVector3& U(state.getPlane()->getU()); |
| 1179 | const TVector3& V(state.getPlane()->getV()); |
| 1180 | TVector3 W(state.getPlane()->getNormal()); |
| 1181 | |
| 1182 | double pw = mom * W; |
| 1183 | double pu = mom * U; |
| 1184 | double pv = mom * V; |
| 1185 | |
| 1186 | TMatrixDSym& cov(state.getCov()); |
| 1187 | |
| 1188 | cov(0,0) = pow(getCharge(state), 2) / pow(mom.Mag(), 6) * |
| 1189 | (mom.X()*mom.X() * momErr.X()*momErr.X()+ |
| 1190 | mom.Y()*mom.Y() * momErr.Y()*momErr.Y()+ |
| 1191 | mom.Z()*mom.Z() * momErr.Z()*momErr.Z()); |
| 1192 | |
| 1193 | cov(1,1) = pow((U.X()/pw - W.X()*pu/(pw*pw)),2.) * momErr.X()*momErr.X() + |
| 1194 | pow((U.Y()/pw - W.Y()*pu/(pw*pw)),2.) * momErr.Y()*momErr.Y() + |
| 1195 | pow((U.Z()/pw - W.Z()*pu/(pw*pw)),2.) * momErr.Z()*momErr.Z(); |
| 1196 | |
| 1197 | cov(2,2) = pow((V.X()/pw - W.X()*pv/(pw*pw)),2.) * momErr.X()*momErr.X() + |
| 1198 | pow((V.Y()/pw - W.Y()*pv/(pw*pw)),2.) * momErr.Y()*momErr.Y() + |
| 1199 | pow((V.Z()/pw - W.Z()*pv/(pw*pw)),2.) * momErr.Z()*momErr.Z(); |
| 1200 | |
| 1201 | cov(3,3) = posErr.X()*posErr.X() * U.X()*U.X() + |
| 1202 | posErr.Y()*posErr.Y() * U.Y()*U.Y() + |
| 1203 | posErr.Z()*posErr.Z() * U.Z()*U.Z(); |
| 1204 | |
| 1205 | cov(4,4) = posErr.X()*posErr.X() * V.X()*V.X() + |
| 1206 | posErr.Y()*posErr.Y() * V.Y()*V.Y() + |
| 1207 | posErr.Z()*posErr.Z() * V.Z()*V.Z(); |
| 1208 | |
| 1209 | } |
| 1210 | |
| 1211 | |
| 1212 | |
| 1213 | |
| 1214 | void RKTrackRep::setPosMomCov(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TMatrixDSym& cov6x6) const { |
| 1215 | |
| 1216 | if (cov6x6.GetNcols()!=6 || cov6x6.GetNrows()!=6){ |
| 1217 | Exception exc("RKTrackRep::setPosMomCov ==> cov has to be 6x6 (x, y, z, px, py, pz)",__LINE__1217,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1218 | throw exc; |
| 1219 | } |
| 1220 | |
| 1221 | setPosMom(state, pos, mom); // charge does not change! |
| 1222 | |
| 1223 | M1x7 state7; |
| 1224 | getState7(state, state7); |
| 1225 | |
| 1226 | const M6x6& cov6x6_( *((M6x6*) cov6x6.GetMatrixArray()) ); |
| 1227 | |
| 1228 | transformM6P(cov6x6_, state7, state); |
| 1229 | |
| 1230 | } |
| 1231 | |
| 1232 | void RKTrackRep::setPosMomCov(MeasuredStateOnPlane& state, const TVectorD& state6, const TMatrixDSym& cov6x6) const { |
| 1233 | |
| 1234 | if (state6.GetNrows()!=6){ |
| 1235 | Exception exc("RKTrackRep::setPosMomCov ==> state has to be 6d (x, y, z, px, py, pz)",__LINE__1235,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1236 | throw exc; |
| 1237 | } |
| 1238 | |
| 1239 | if (cov6x6.GetNcols()!=6 || cov6x6.GetNrows()!=6){ |
| 1240 | Exception exc("RKTrackRep::setPosMomCov ==> cov has to be 6x6 (x, y, z, px, py, pz)",__LINE__1240,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1241 | throw exc; |
| 1242 | } |
| 1243 | |
| 1244 | TVector3 pos(state6(0), state6(1), state6(2)); |
| 1245 | TVector3 mom(state6(3), state6(4), state6(5)); |
| 1246 | setPosMom(state, pos, mom); // charge does not change! |
| 1247 | |
| 1248 | M1x7 state7; |
| 1249 | getState7(state, state7); |
| 1250 | |
| 1251 | const M6x6& cov6x6_( *((M6x6*) cov6x6.GetMatrixArray()) ); |
| 1252 | |
| 1253 | transformM6P(cov6x6_, state7, state); |
| 1254 | |
| 1255 | } |
| 1256 | |
| 1257 | |
| 1258 | void RKTrackRep::setChargeSign(StateOnPlane& state, double charge) const { |
| 1259 | |
| 1260 | if (dynamic_cast<MeasurementOnPlane*>(&state) != nullptr) { |
| 1261 | Exception exc("RKTrackRep::setChargeSign - cannot set charge of a MeasurementOnPlane",__LINE__1261,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1262 | exc.setFatal(); |
| 1263 | throw exc; |
| 1264 | } |
| 1265 | |
| 1266 | if (state.getState()(0) * charge < 0) { |
| 1267 | state.getState()(0) *= -1.; |
| 1268 | } |
| 1269 | } |
| 1270 | |
| 1271 | |
| 1272 | void RKTrackRep::setSpu(StateOnPlane& state, double spu) const { |
| 1273 | state.getAuxInfo().ResizeTo(2); |
| 1274 | (state.getAuxInfo())(0) = spu; |
| 1275 | } |
| 1276 | |
| 1277 | void RKTrackRep::setTime(StateOnPlane& state, double time) const { |
| 1278 | state.getAuxInfo().ResizeTo(2); |
| 1279 | (state.getAuxInfo())(1) = time; |
| 1280 | } |
| 1281 | |
| 1282 | |
| 1283 | |
| 1284 | double RKTrackRep::RKPropagate(M1x7& state7, |
| 1285 | M7x7* jacobianT, |
| 1286 | M1x3& SA, |
| 1287 | double S, |
| 1288 | bool varField, |
| 1289 | bool calcOnlyLastRowOfJ) const { |
| 1290 | // The algorithm is |
| 1291 | // E Lund et al 2009 JINST 4 P04001 doi:10.1088/1748-0221/4/04/P04001 |
| 1292 | // "Track parameter propagation through the application of a new adaptive Runge-Kutta-Nyström method in the ATLAS experiment" |
| 1293 | // http://inspirehep.net/search?ln=en&ln=en&p=10.1088/1748-0221/4/04/P04001&of=hb&action_search=Search&sf=earliestdate&so=d&rm=&rg=25&sc=0 |
| 1294 | // where the transport of the Jacobian is described in |
| 1295 | // L. Bugge, J. Myrheim Nucl.Instrum.Meth. 160 (1979) 43-48 |
| 1296 | // "A Fast Runge-kutta Method For Fitting Tracks In A Magnetic Field" |
| 1297 | // http://inspirehep.net/record/145692 |
| 1298 | // and |
| 1299 | // L. Bugge, J. Myrheim Nucl.Instrum.Meth. 179 (1981) 365-381 |
| 1300 | // "Tracking And Track Fitting" |
| 1301 | // http://inspirehep.net/record/160548 |
| 1302 | |
| 1303 | // important fixed numbers |
| 1304 | static const double EC ( 0.000149896229 ); // c/(2*10^12) resp. c/2Tera |
| 1305 | static const double P3 ( 1./3. ); // 1/3 |
| 1306 | static const double DLT ( .0002 ); // max. deviation for approximation-quality test |
| 1307 | // Aux parameters |
| 1308 | M1x3& R = *((M1x3*) &state7[0]); // Start coordinates [cm] (x, y, z) |
| 1309 | M1x3& A = *((M1x3*) &state7[3]); // Start directions (ax, ay, az); ax^2+ay^2+az^2=1 |
| 1310 | double S3(0), S4(0), PS2(0); |
| 1311 | M1x3 H0 = {{0.,0.,0.}}, H1 = {{0.,0.,0.}}, H2 = {{0.,0.,0.}}; |
| 1312 | M1x3 r = {{0.,0.,0.}}; |
| 1313 | // Variables for Runge Kutta solver |
| 1314 | double A0(0), A1(0), A2(0), A3(0), A4(0), A5(0), A6(0); |
| 1315 | double B0(0), B1(0), B2(0), B3(0), B4(0), B5(0), B6(0); |
| 1316 | double C0(0), C1(0), C2(0), C3(0), C4(0), C5(0), C6(0); |
| 1317 | |
| 1318 | // |
| 1319 | // Runge Kutta Extrapolation |
| 1320 | // |
| 1321 | S3 = P3*S; |
| 1322 | S4 = 0.25*S; |
| 1323 | PS2 = state7[6]*EC * S; |
| 1324 | |
| 1325 | // First point |
| 1326 | r[0] = R[0]; r[1] = R[1]; r[2]=R[2]; |
| 1327 | FieldManager::getInstance()->getFieldVal(r[0], r[1], r[2], H0[0], H0[1], H0[2]); // magnetic field in 10^-1 T = kGauss |
| 1328 | H0[0] *= PS2; H0[1] *= PS2; H0[2] *= PS2; // H0 is PS2*(Hx, Hy, Hz) @ R0 |
| 1329 | A0 = A[1]*H0[2]-A[2]*H0[1]; B0 = A[2]*H0[0]-A[0]*H0[2]; C0 = A[0]*H0[1]-A[1]*H0[0]; // (ax, ay, az) x H0 |
| 1330 | A2 = A[0]+A0 ; B2 = A[1]+B0 ; C2 = A[2]+C0 ; // (A0, B0, C0) + (ax, ay, az) |
| 1331 | A1 = A2+A[0] ; B1 = B2+A[1] ; C1 = C2+A[2] ; // (A0, B0, C0) + 2*(ax, ay, az) |
| 1332 | |
| 1333 | // Second point |
| 1334 | if (varField) { |
| 1335 | r[0] += A1*S4; r[1] += B1*S4; r[2] += C1*S4; |
| 1336 | FieldManager::getInstance()->getFieldVal(r[0], r[1], r[2], H1[0], H1[1], H1[2]); |
| 1337 | H1[0] *= PS2; H1[1] *= PS2; H1[2] *= PS2; // H1 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * [(A0, B0, C0) + 2*(ax, ay, az)] |
| 1338 | } |
| 1339 | else { H1 = H0; }; |
| 1340 | A3 = B2*H1[2]-C2*H1[1]+A[0]; B3 = C2*H1[0]-A2*H1[2]+A[1]; C3 = A2*H1[1]-B2*H1[0]+A[2]; // (A2, B2, C2) x H1 + (ax, ay, az) |
| 1341 | A4 = B3*H1[2]-C3*H1[1]+A[0]; B4 = C3*H1[0]-A3*H1[2]+A[1]; C4 = A3*H1[1]-B3*H1[0]+A[2]; // (A3, B3, C3) x H1 + (ax, ay, az) |
| 1342 | A5 = A4-A[0]+A4 ; B5 = B4-A[1]+B4 ; C5 = C4-A[2]+C4 ; // 2*(A4, B4, C4) - (ax, ay, az) |
| 1343 | |
| 1344 | // Last point |
| 1345 | if (varField) { |
| 1346 | r[0]=R[0]+S*A4; r[1]=R[1]+S*B4; r[2]=R[2]+S*C4; //setup.Field(r,H2); |
| 1347 | FieldManager::getInstance()->getFieldVal(r[0], r[1], r[2], H2[0], H2[1], H2[2]); |
| 1348 | H2[0] *= PS2; H2[1] *= PS2; H2[2] *= PS2; // H2 is PS2*(Hx, Hy, Hz) @ (x, y, z) + 0.25*S * (A4, B4, C4) |
| 1349 | } |
| 1350 | else { H2 = H0; }; |
| 1351 | A6 = B5*H2[2]-C5*H2[1]; B6 = C5*H2[0]-A5*H2[2]; C6 = A5*H2[1]-B5*H2[0]; // (A5, B5, C5) x H2 |
| 1352 | |
| 1353 | |
| 1354 | // |
| 1355 | // Derivatives of track parameters |
| 1356 | // |
| 1357 | if(jacobianT != nullptr){ |
| 1358 | |
| 1359 | // jacobianT |
| 1360 | // 1 0 0 0 0 0 0 x |
| 1361 | // 0 1 0 0 0 0 0 y |
| 1362 | // 0 0 1 0 0 0 0 z |
| 1363 | // x x x x x x 0 a_x |
| 1364 | // x x x x x x 0 a_y |
| 1365 | // x x x x x x 0 a_z |
| 1366 | // x x x x x x 1 q/p |
| 1367 | M7x7& J = *jacobianT; |
| 1368 | |
| 1369 | double dA0(0), dA2(0), dA3(0), dA4(0), dA5(0), dA6(0); |
| 1370 | double dB0(0), dB2(0), dB3(0), dB4(0), dB5(0), dB6(0); |
| 1371 | double dC0(0), dC2(0), dC3(0), dC4(0), dC5(0), dC6(0); |
| 1372 | |
| 1373 | int start(0); |
| 1374 | |
| 1375 | if (!calcOnlyLastRowOfJ) { |
| 1376 | |
| 1377 | if (!varField) { |
| 1378 | // d(x, y, z)/d(x, y, z) submatrix is unit matrix |
| 1379 | J(0, 0) = 1; J(1, 1) = 1; J(2, 2) = 1; |
| 1380 | // d(ax, ay, az)/d(ax, ay, az) submatrix is 0 |
| 1381 | // start with d(x, y, z)/d(ax, ay, az) |
| 1382 | start = 3; |
| 1383 | } |
| 1384 | |
| 1385 | for(int i=start; i<6; ++i) { |
| 1386 | |
| 1387 | //first point |
| 1388 | dA0 = H0[2]*J(i, 4)-H0[1]*J(i, 5); // dA0/dp } |
| 1389 | dB0 = H0[0]*J(i, 5)-H0[2]*J(i, 3); // dB0/dp } = dA x H0 |
| 1390 | dC0 = H0[1]*J(i, 3)-H0[0]*J(i, 4); // dC0/dp } |
| 1391 | |
| 1392 | dA2 = dA0+J(i, 3); // } |
| 1393 | dB2 = dB0+J(i, 4); // } = (dA0, dB0, dC0) + dA |
| 1394 | dC2 = dC0+J(i, 5); // } |
| 1395 | |
| 1396 | //second point |
| 1397 | dA3 = J(i, 3)+dB2*H1[2]-dC2*H1[1]; // dA3/dp } |
| 1398 | dB3 = J(i, 4)+dC2*H1[0]-dA2*H1[2]; // dB3/dp } = dA + (dA2, dB2, dC2) x H1 |
| 1399 | dC3 = J(i, 5)+dA2*H1[1]-dB2*H1[0]; // dC3/dp } |
| 1400 | |
| 1401 | dA4 = J(i, 3)+dB3*H1[2]-dC3*H1[1]; // dA4/dp } |
| 1402 | dB4 = J(i, 4)+dC3*H1[0]-dA3*H1[2]; // dB4/dp } = dA + (dA3, dB3, dC3) x H1 |
| 1403 | dC4 = J(i, 5)+dA3*H1[1]-dB3*H1[0]; // dC4/dp } |
| 1404 | |
| 1405 | //last point |
| 1406 | dA5 = dA4+dA4-J(i, 3); // } |
| 1407 | dB5 = dB4+dB4-J(i, 4); // } = 2*(dA4, dB4, dC4) - dA |
| 1408 | dC5 = dC4+dC4-J(i, 5); // } |
| 1409 | |
| 1410 | dA6 = dB5*H2[2]-dC5*H2[1]; // dA6/dp } |
| 1411 | dB6 = dC5*H2[0]-dA5*H2[2]; // dB6/dp } = (dA5, dB5, dC5) x H2 |
| 1412 | dC6 = dA5*H2[1]-dB5*H2[0]; // dC6/dp } |
| 1413 | |
| 1414 | // this gives the same results as multiplying the old with the new Jacobian |
| 1415 | J(i, 0) += (dA2+dA3+dA4)*S3; J(i, 3) = ((dA0+2.*dA3)+(dA5+dA6))*P3; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)] |
| 1416 | J(i, 1) += (dB2+dB3+dB4)*S3; J(i, 4) = ((dB0+2.*dB3)+(dB5+dB6))*P3; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)] |
| 1417 | J(i, 2) += (dC2+dC3+dC4)*S3; J(i, 5) = ((dC0+2.*dC3)+(dC5+dC6))*P3; |
| 1418 | } |
| 1419 | |
| 1420 | } // end if (!calcOnlyLastRowOfJ) |
| 1421 | |
| 1422 | J(6, 3) *= state7[6]; J(6, 4) *= state7[6]; J(6, 5) *= state7[6]; |
| 1423 | |
| 1424 | //first point |
| 1425 | dA0 = H0[2]*J(6, 4)-H0[1]*J(6, 5) + A0; // dA0/dp } |
| 1426 | dB0 = H0[0]*J(6, 5)-H0[2]*J(6, 3) + B0; // dB0/dp } = dA x H0 + (A0, B0, C0) |
| 1427 | dC0 = H0[1]*J(6, 3)-H0[0]*J(6, 4) + C0; // dC0/dp } |
| 1428 | |
| 1429 | dA2 = dA0+J(6, 3); // } |
| 1430 | dB2 = dB0+J(6, 4); // } = (dA0, dB0, dC0) + dA |
| 1431 | dC2 = dC0+J(6, 5); // } |
| 1432 | |
| 1433 | //second point |
| 1434 | dA3 = J(6, 3)+dB2*H1[2]-dC2*H1[1] + (A3-A[0]); // dA3/dp } |
| 1435 | dB3 = J(6, 4)+dC2*H1[0]-dA2*H1[2] + (B3-A[1]); // dB3/dp } = dA + (dA2, dB2, dC2) x H1 |
| 1436 | dC3 = J(6, 5)+dA2*H1[1]-dB2*H1[0] + (C3-A[2]); // dC3/dp } |
| 1437 | |
| 1438 | dA4 = J(6, 3)+dB3*H1[2]-dC3*H1[1] + (A4-A[0]); // dA4/dp } |
| 1439 | dB4 = J(6, 4)+dC3*H1[0]-dA3*H1[2] + (B4-A[1]); // dB4/dp } = dA + (dA3, dB3, dC3) x H1 |
| 1440 | dC4 = J(6, 5)+dA3*H1[1]-dB3*H1[0] + (C4-A[2]); // dC4/dp } |
| 1441 | |
| 1442 | //last point |
| 1443 | dA5 = dA4+dA4-J(6, 3); // } |
| 1444 | dB5 = dB4+dB4-J(6, 4); // } = 2*(dA4, dB4, dC4) - dA |
| 1445 | dC5 = dC4+dC4-J(6, 5); // } |
| 1446 | |
| 1447 | dA6 = dB5*H2[2]-dC5*H2[1] + A6; // dA6/dp } |
| 1448 | dB6 = dC5*H2[0]-dA5*H2[2] + B6; // dB6/dp } = (dA5, dB5, dC5) x H2 + (A6, B6, C6) |
| 1449 | dC6 = dA5*H2[1]-dB5*H2[0] + C6; // dC6/dp } |
| 1450 | |
| 1451 | // this gives the same results as multiplying the old with the new Jacobian |
| 1452 | J(6, 0) += (dA2+dA3+dA4)*S3/state7[6]; J(6, 3) = ((dA0+2.*dA3)+(dA5+dA6))*P3/state7[6]; // dR := dR + S3*[(dA2, dB2, dC2) + (dA3, dB3, dC3) + (dA4, dB4, dC4)] |
| 1453 | J(6, 1) += (dB2+dB3+dB4)*S3/state7[6]; J(6, 4) = ((dB0+2.*dB3)+(dB5+dB6))*P3/state7[6]; // dA := 1/3*[(dA0, dB0, dC0) + 2*(dA3, dB3, dC3) + (dA5, dB5, dC5) + (dA6, dB6, dC6)] |
| 1454 | J(6, 2) += (dC2+dC3+dC4)*S3/state7[6]; J(6, 5) = ((dC0+2.*dC3)+(dC5+dC6))*P3/state7[6]; |
| 1455 | |
| 1456 | } |
| 1457 | |
| 1458 | // |
| 1459 | // Track parameters in last point |
| 1460 | // |
| 1461 | R[0] += (A2+A3+A4)*S3; A[0] += (SA[0]=((A0+2.*A3)+(A5+A6))*P3-A[0]); // R = R0 + S3*[(A2, B2, C2) + (A3, B3, C3) + (A4, B4, C4)] |
| 1462 | R[1] += (B2+B3+B4)*S3; A[1] += (SA[1]=((B0+2.*B3)+(B5+B6))*P3-A[1]); // A = 1/3*[(A0, B0, C0) + 2*(A3, B3, C3) + (A5, B5, C5) + (A6, B6, C6)] |
| 1463 | R[2] += (C2+C3+C4)*S3; A[2] += (SA[2]=((C0+2.*C3)+(C5+C6))*P3-A[2]); // SA = A_new - A_old |
| 1464 | |
| 1465 | // normalize A |
| 1466 | double CBA ( 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]) ); // 1/|A| |
| 1467 | A[0] *= CBA; A[1] *= CBA; A[2] *= CBA; |
| 1468 | |
| 1469 | |
| 1470 | // Test approximation quality on given step |
| 1471 | double EST ( fabs((A1+A6)-(A3+A4)) + |
| 1472 | fabs((B1+B6)-(B3+B4)) + |
| 1473 | fabs((C1+C6)-(C3+C4)) ); // EST = ||(ABC1+ABC6)-(ABC3+ABC4)||_1 = ||(axzy x H0 + ABC5 x H2) - (ABC2 x H1 + ABC3 x H1)||_1 |
| 1474 | if (debugLvl_ > 0) { |
| 1475 | debugOut << " RKTrackRep::RKPropagate. Step = "<< S << "; quality EST = " << EST << " \n"; |
| 1476 | } |
| 1477 | |
| 1478 | // Prevent the step length increase from getting too large, this is |
| 1479 | // just the point where it becomes 10. |
| 1480 | if (EST < DLT*1e-5) |
| 1481 | return 10; |
| 1482 | |
| 1483 | // Step length increase for a fifth order Runge-Kutta, see e.g. 17.2 |
| 1484 | // in Numerical Recipes. FIXME: move to caller. |
| 1485 | return pow(DLT/EST, 1./5.); |
| 1486 | } |
| 1487 | |
| 1488 | |
| 1489 | |
| 1490 | void RKTrackRep::initArrays() const { |
| 1491 | std::fill(noiseArray_.begin(), noiseArray_.end(), 0); |
| 1492 | std::fill(noiseProjection_.begin(), noiseProjection_.end(), 0); |
| 1493 | for (unsigned int i=0; i<7; ++i) // initialize as diagonal matrix |
| 1494 | noiseProjection_[i*8] = 1; |
| 1495 | std::fill(J_MMT_.begin(), J_MMT_.end(), 0); |
| 1496 | |
| 1497 | fJacobian_.UnitMatrix(); |
| 1498 | fNoise_.Zero(); |
| 1499 | limits_.reset(); |
| 1500 | |
| 1501 | RKSteps_.reserve(100); |
| 1502 | ExtrapSteps_.reserve(100); |
| 1503 | |
| 1504 | lastStartState_.getAuxInfo().ResizeTo(2); |
| 1505 | lastEndState_.getAuxInfo().ResizeTo(2); |
| 1506 | } |
| 1507 | |
| 1508 | |
| 1509 | void RKTrackRep::getState7(const StateOnPlane& state, M1x7& state7) const { |
| 1510 | |
| 1511 | if (dynamic_cast<const MeasurementOnPlane*>(&state) != nullptr) { |
| 1512 | Exception exc("RKTrackRep::getState7 - cannot get pos or mom from a MeasurementOnPlane",__LINE__1512,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1513 | exc.setFatal(); |
| 1514 | throw exc; |
| 1515 | } |
| 1516 | |
| 1517 | const TVector3& U(state.getPlane()->getU()); |
| 1518 | const TVector3& V(state.getPlane()->getV()); |
| 1519 | const TVector3& O(state.getPlane()->getO()); |
| 1520 | const TVector3& W(state.getPlane()->getNormal()); |
| 1521 | |
| 1522 | assert(state.getState().GetNrows() == 5)(static_cast <bool> (state.getState().GetNrows() == 5) ? void (0) : __assert_fail ("state.getState().GetNrows() == 5" , "genfit2/code2/trackReps/src/RKTrackRep.cc", 1522, __extension__ __PRETTY_FUNCTION__)); |
| 1523 | const double* state5 = state.getState().GetMatrixArray(); |
| 1524 | |
| 1525 | double spu = getSpu(state); |
| 1526 | |
| 1527 | state7[0] = O.X() + state5[3]*U.X() + state5[4]*V.X(); // x |
| 1528 | state7[1] = O.Y() + state5[3]*U.Y() + state5[4]*V.Y(); // y |
| 1529 | state7[2] = O.Z() + state5[3]*U.Z() + state5[4]*V.Z(); // z |
| 1530 | |
| 1531 | state7[3] = spu * (W.X() + state5[1]*U.X() + state5[2]*V.X()); // a_x |
| 1532 | state7[4] = spu * (W.Y() + state5[1]*U.Y() + state5[2]*V.Y()); // a_y |
| 1533 | state7[5] = spu * (W.Z() + state5[1]*U.Z() + state5[2]*V.Z()); // a_z |
| 1534 | |
| 1535 | // normalize dir |
| 1536 | double norm = 1. / sqrt(state7[3]*state7[3] + state7[4]*state7[4] + state7[5]*state7[5]); |
| 1537 | for (unsigned int i=3; i<6; ++i) state7[i] *= norm; |
| 1538 | |
| 1539 | state7[6] = state5[0]; // q/p |
| 1540 | } |
| 1541 | |
| 1542 | |
| 1543 | void RKTrackRep::getState5(StateOnPlane& state, const M1x7& state7) const { |
| 1544 | |
| 1545 | // state5: (q/p, u', v'. u, v) |
| 1546 | |
| 1547 | double spu(1.); |
| 1548 | |
| 1549 | const TVector3& O(state.getPlane()->getO()); |
| 1550 | const TVector3& U(state.getPlane()->getU()); |
| 1551 | const TVector3& V(state.getPlane()->getV()); |
| 1552 | const TVector3& W(state.getPlane()->getNormal()); |
| 1553 | |
| 1554 | // force A to be in normal direction and set spu accordingly |
| 1555 | double AtW( state7[3]*W.X() + state7[4]*W.Y() + state7[5]*W.Z() ); |
| 1556 | if (AtW < 0.) { |
| 1557 | //fDir *= -1.; |
| 1558 | //AtW *= -1.; |
| 1559 | spu = -1.; |
| 1560 | } |
| 1561 | |
| 1562 | double* state5 = state.getState().GetMatrixArray(); |
| 1563 | |
| 1564 | state5[0] = state7[6]; // q/p |
| 1565 | state5[1] = (state7[3]*U.X() + state7[4]*U.Y() + state7[5]*U.Z()) / AtW; // u' = (A * U) / (A * W) |
| 1566 | state5[2] = (state7[3]*V.X() + state7[4]*V.Y() + state7[5]*V.Z()) / AtW; // v' = (A * V) / (A * W) |
| 1567 | state5[3] = ((state7[0]-O.X())*U.X() + |
| 1568 | (state7[1]-O.Y())*U.Y() + |
| 1569 | (state7[2]-O.Z())*U.Z()); // u = (pos - O) * U |
| 1570 | state5[4] = ((state7[0]-O.X())*V.X() + |
| 1571 | (state7[1]-O.Y())*V.Y() + |
| 1572 | (state7[2]-O.Z())*V.Z()); // v = (pos - O) * V |
| 1573 | |
| 1574 | setSpu(state, spu); |
| 1575 | |
| 1576 | } |
| 1577 | |
| 1578 | void RKTrackRep::calcJ_pM_5x7(M5x7& J_pM, const TVector3& U, const TVector3& V, const M1x3& pTilde, double spu) const { |
| 1579 | /*if (debugLvl_ > 1) { |
| 1580 | debugOut << "RKTrackRep::calcJ_pM_5x7 \n"; |
| 1581 | debugOut << " U = "; U.Print(); |
| 1582 | debugOut << " V = "; V.Print(); |
| 1583 | debugOut << " pTilde = "; RKTools::printDim(pTilde, 3,1); |
| 1584 | debugOut << " spu = " << spu << "\n"; |
| 1585 | }*/ |
| 1586 | |
| 1587 | std::fill(J_pM.begin(), J_pM.end(), 0); |
| 1588 | |
| 1589 | const double pTildeMag = sqrt(pTilde[0]*pTilde[0] + pTilde[1]*pTilde[1] + pTilde[2]*pTilde[2]); |
| 1590 | const double pTildeMag2 = pTildeMag*pTildeMag; |
| 1591 | |
| 1592 | const double utpTildeOverpTildeMag2 = (U.X()*pTilde[0] + U.Y()*pTilde[1] + U.Z()*pTilde[2]) / pTildeMag2; |
| 1593 | const double vtpTildeOverpTildeMag2 = (V.X()*pTilde[0] + V.Y()*pTilde[1] + V.Z()*pTilde[2]) / pTildeMag2; |
| 1594 | |
| 1595 | //J_pM matrix is d(x,y,z,ax,ay,az,q/p) / d(q/p,u',v',u,v) (out is 7x7) |
| 1596 | |
| 1597 | // d(x,y,z)/d(u) |
| 1598 | J_pM[21] = U.X(); // [3][0] |
| 1599 | J_pM[22] = U.Y(); // [3][1] |
| 1600 | J_pM[23] = U.Z(); // [3][2] |
| 1601 | // d(x,y,z)/d(v) |
| 1602 | J_pM[28] = V.X(); // [4][0] |
| 1603 | J_pM[29] = V.Y(); // [4][1] |
| 1604 | J_pM[30] = V.Z(); // [4][2] |
| 1605 | // d(q/p)/d(q/p) |
| 1606 | J_pM[6] = 1.; // not needed for array matrix multiplication |
| 1607 | // d(ax,ay,az)/d(u') |
| 1608 | double fact = spu / pTildeMag; |
| 1609 | J_pM[10] = fact * ( U.X() - pTilde[0]*utpTildeOverpTildeMag2 ); // [1][3] |
| 1610 | J_pM[11] = fact * ( U.Y() - pTilde[1]*utpTildeOverpTildeMag2 ); // [1][4] |
| 1611 | J_pM[12] = fact * ( U.Z() - pTilde[2]*utpTildeOverpTildeMag2 ); // [1][5] |
| 1612 | // d(ax,ay,az)/d(v') |
| 1613 | J_pM[17] = fact * ( V.X() - pTilde[0]*vtpTildeOverpTildeMag2 ); // [2][3] |
| 1614 | J_pM[18] = fact * ( V.Y() - pTilde[1]*vtpTildeOverpTildeMag2 ); // [2][4] |
| 1615 | J_pM[19] = fact * ( V.Z() - pTilde[2]*vtpTildeOverpTildeMag2 ); // [2][5] |
| 1616 | |
| 1617 | /*if (debugLvl_ > 1) { |
| 1618 | debugOut << " J_pM_5x7_ = "; RKTools::printDim(J_pM_5x7_, 5,7); |
| 1619 | }*/ |
| 1620 | } |
| 1621 | |
| 1622 | |
| 1623 | void RKTrackRep::transformPM6(const MeasuredStateOnPlane& state, |
| 1624 | M6x6& out6x6) const { |
| 1625 | |
| 1626 | // get vectors and aux variables |
| 1627 | const TVector3& U(state.getPlane()->getU()); |
| 1628 | const TVector3& V(state.getPlane()->getV()); |
| 1629 | const TVector3& W(state.getPlane()->getNormal()); |
| 1630 | |
| 1631 | const TVectorD& state5(state.getState()); |
| 1632 | double spu(getSpu(state)); |
| 1633 | |
| 1634 | M1x3 pTilde; |
| 1635 | pTilde[0] = spu * (W.X() + state5(1)*U.X() + state5(2)*V.X()); // a_x |
| 1636 | pTilde[1] = spu * (W.Y() + state5(1)*U.Y() + state5(2)*V.Y()); // a_y |
| 1637 | pTilde[2] = spu * (W.Z() + state5(1)*U.Z() + state5(2)*V.Z()); // a_z |
| 1638 | |
| 1639 | const double pTildeMag = sqrt(pTilde[0]*pTilde[0] + pTilde[1]*pTilde[1] + pTilde[2]*pTilde[2]); |
| 1640 | const double pTildeMag2 = pTildeMag*pTildeMag; |
| 1641 | |
| 1642 | const double utpTildeOverpTildeMag2 = (U.X()*pTilde[0] + U.Y()*pTilde[1] + U.Z()*pTilde[2]) / pTildeMag2; |
| 1643 | const double vtpTildeOverpTildeMag2 = (V.X()*pTilde[0] + V.Y()*pTilde[1] + V.Z()*pTilde[2]) / pTildeMag2; |
| 1644 | |
| 1645 | //J_pM matrix is d(x,y,z,px,py,pz) / d(q/p,u',v',u,v) (out is 6x6) |
| 1646 | |
| 1647 | const double qop = state5(0); |
| 1648 | const double p = getCharge(state)/qop; // momentum |
| 1649 | |
| 1650 | M5x6 J_pM_5x6; |
| 1651 | std::fill(J_pM_5x6.begin(), J_pM_5x6.end(), 0); |
| 1652 | |
| 1653 | // d(px,py,pz)/d(q/p) |
| 1654 | double fact = -1. * p / (pTildeMag * qop); |
| 1655 | J_pM_5x6[3] = fact * pTilde[0]; // [0][3] |
| 1656 | J_pM_5x6[4] = fact * pTilde[1]; // [0][4] |
| 1657 | J_pM_5x6[5] = fact * pTilde[2]; // [0][5] |
| 1658 | // d(px,py,pz)/d(u') |
| 1659 | fact = p * spu / pTildeMag; |
| 1660 | J_pM_5x6[9] = fact * ( U.X() - pTilde[0]*utpTildeOverpTildeMag2 ); // [1][3] |
| 1661 | J_pM_5x6[10] = fact * ( U.Y() - pTilde[1]*utpTildeOverpTildeMag2 ); // [1][4] |
| 1662 | J_pM_5x6[11] = fact * ( U.Z() - pTilde[2]*utpTildeOverpTildeMag2 ); // [1][5] |
| 1663 | // d(px,py,pz)/d(v') |
| 1664 | J_pM_5x6[15] = fact * ( V.X() - pTilde[0]*vtpTildeOverpTildeMag2 ); // [2][3] |
| 1665 | J_pM_5x6[16] = fact * ( V.Y() - pTilde[1]*vtpTildeOverpTildeMag2 ); // [2][4] |
| 1666 | J_pM_5x6[17] = fact * ( V.Z() - pTilde[2]*vtpTildeOverpTildeMag2 ); // [2][5] |
| 1667 | // d(x,y,z)/d(u) |
| 1668 | J_pM_5x6[18] = U.X(); // [3][0] |
| 1669 | J_pM_5x6[19] = U.Y(); // [3][1] |
| 1670 | J_pM_5x6[20] = U.Z(); // [3][2] |
| 1671 | // d(x,y,z)/d(v) |
| 1672 | J_pM_5x6[24] = V.X(); // [4][0] |
| 1673 | J_pM_5x6[25] = V.Y(); // [4][1] |
| 1674 | J_pM_5x6[26] = V.Z(); // [4][2] |
| 1675 | |
| 1676 | |
| 1677 | // do the transformation |
| 1678 | // out = J_pM^T * in5x5 * J_pM |
| 1679 | const M5x5& in5x5_ = *((M5x5*) state.getCov().GetMatrixArray()); |
| 1680 | RKTools::J_pMTxcov5xJ_pM(J_pM_5x6, in5x5_, out6x6); |
| 1681 | |
| 1682 | } |
| 1683 | |
| 1684 | void RKTrackRep::calcJ_Mp_7x5(M7x5& J_Mp, const TVector3& U, const TVector3& V, const TVector3& W, const M1x3& A) const { |
| 1685 | |
| 1686 | /*if (debugLvl_ > 1) { |
| 1687 | debugOut << "RKTrackRep::calcJ_Mp_7x5 \n"; |
| 1688 | debugOut << " U = "; U.Print(); |
| 1689 | debugOut << " V = "; V.Print(); |
| 1690 | debugOut << " W = "; W.Print(); |
| 1691 | debugOut << " A = "; RKTools::printDim(A, 3,1); |
| 1692 | }*/ |
| 1693 | |
| 1694 | std::fill(J_Mp.begin(), J_Mp.end(), 0); |
| 1695 | |
| 1696 | const double AtU = A[0]*U.X() + A[1]*U.Y() + A[2]*U.Z(); |
| 1697 | const double AtV = A[0]*V.X() + A[1]*V.Y() + A[2]*V.Z(); |
| 1698 | const double AtW = A[0]*W.X() + A[1]*W.Y() + A[2]*W.Z(); |
| 1699 | |
| 1700 | // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,ax,ay,az,q/p) (in is 7x7) |
| 1701 | |
| 1702 | // d(u')/d(ax,ay,az) |
| 1703 | double fact = 1./(AtW*AtW); |
| 1704 | J_Mp[16] = fact * (U.X()*AtW - W.X()*AtU); // [3][1] |
| 1705 | J_Mp[21] = fact * (U.Y()*AtW - W.Y()*AtU); // [4][1] |
| 1706 | J_Mp[26] = fact * (U.Z()*AtW - W.Z()*AtU); // [5][1] |
| 1707 | // d(v')/d(ax,ay,az) |
| 1708 | J_Mp[17] = fact * (V.X()*AtW - W.X()*AtV); // [3][2] |
| 1709 | J_Mp[22] = fact * (V.Y()*AtW - W.Y()*AtV); // [4][2] |
| 1710 | J_Mp[27] = fact * (V.Z()*AtW - W.Z()*AtV); // [5][2] |
| 1711 | // d(q/p)/d(q/p) |
| 1712 | J_Mp[30] = 1.; // [6][0] - not needed for array matrix multiplication |
| 1713 | //d(u)/d(x,y,z) |
| 1714 | J_Mp[3] = U.X(); // [0][3] |
| 1715 | J_Mp[8] = U.Y(); // [1][3] |
| 1716 | J_Mp[13] = U.Z(); // [2][3] |
| 1717 | //d(v)/d(x,y,z) |
| 1718 | J_Mp[4] = V.X(); // [0][4] |
| 1719 | J_Mp[9] = V.Y(); // [1][4] |
| 1720 | J_Mp[14] = V.Z(); // [2][4] |
| 1721 | |
| 1722 | /*if (debugLvl_ > 1) { |
| 1723 | debugOut << " J_Mp_7x5_ = "; RKTools::printDim(J_Mp, 7,5); |
| 1724 | }*/ |
| 1725 | |
| 1726 | } |
| 1727 | |
| 1728 | |
| 1729 | void RKTrackRep::transformM6P(const M6x6& in6x6, |
| 1730 | const M1x7& state7, |
| 1731 | MeasuredStateOnPlane& state) const { // plane and charge must already be set! |
| 1732 | |
| 1733 | // get vectors and aux variables |
| 1734 | const TVector3& U(state.getPlane()->getU()); |
| 1735 | const TVector3& V(state.getPlane()->getV()); |
| 1736 | const TVector3& W(state.getPlane()->getNormal()); |
| 1737 | |
| 1738 | const double AtU = state7[3]*U.X() + state7[4]*U.Y() + state7[5]*U.Z(); |
| 1739 | const double AtV = state7[3]*V.X() + state7[4]*V.Y() + state7[5]*V.Z(); |
| 1740 | const double AtW = state7[3]*W.X() + state7[4]*W.Y() + state7[5]*W.Z(); |
| 1741 | |
| 1742 | // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,px,py,pz) (in is 6x6) |
| 1743 | |
| 1744 | const double qop = state7[6]; |
| 1745 | const double p = getCharge(state)/qop; // momentum |
| 1746 | |
| 1747 | M6x5 J_Mp_6x5; |
| 1748 | std::fill(J_Mp_6x5.begin(), J_Mp_6x5.end(), 0); |
| 1749 | |
| 1750 | //d(u)/d(x,y,z) |
| 1751 | J_Mp_6x5[3] = U.X(); // [0][3] |
| 1752 | J_Mp_6x5[8] = U.Y(); // [1][3] |
| 1753 | J_Mp_6x5[13] = U.Z(); // [2][3] |
| 1754 | //d(v)/d(x,y,z) |
| 1755 | J_Mp_6x5[4] = V.X(); // [0][4] |
| 1756 | J_Mp_6x5[9] = V.Y(); // [1][4] |
| 1757 | J_Mp_6x5[14] = V.Z(); // [2][4] |
| 1758 | // d(q/p)/d(px,py,pz) |
| 1759 | double fact = (-1.) * qop / p; |
| 1760 | J_Mp_6x5[15] = fact * state7[3]; // [3][0] |
| 1761 | J_Mp_6x5[20] = fact * state7[4]; // [4][0] |
| 1762 | J_Mp_6x5[25] = fact * state7[5]; // [5][0] |
| 1763 | // d(u')/d(px,py,pz) |
| 1764 | fact = 1./(p*AtW*AtW); |
| 1765 | J_Mp_6x5[16] = fact * (U.X()*AtW - W.X()*AtU); // [3][1] |
| 1766 | J_Mp_6x5[21] = fact * (U.Y()*AtW - W.Y()*AtU); // [4][1] |
| 1767 | J_Mp_6x5[26] = fact * (U.Z()*AtW - W.Z()*AtU); // [5][1] |
| 1768 | // d(v')/d(px,py,pz) |
| 1769 | J_Mp_6x5[17] = fact * (V.X()*AtW - W.X()*AtV); // [3][2] |
| 1770 | J_Mp_6x5[22] = fact * (V.Y()*AtW - W.Y()*AtV); // [4][2] |
| 1771 | J_Mp_6x5[27] = fact * (V.Z()*AtW - W.Z()*AtV); // [5][2] |
| 1772 | |
| 1773 | // do the transformation |
| 1774 | // out5x5 = J_Mp^T * in * J_Mp |
| 1775 | M5x5& out5x5_ = *((M5x5*) state.getCov().GetMatrixArray()); |
| 1776 | RKTools::J_MpTxcov6xJ_Mp(J_Mp_6x5, in6x6, out5x5_); |
| 1777 | |
| 1778 | } |
| 1779 | |
| 1780 | |
| 1781 | // |
| 1782 | // Runge-Kutta method for tracking a particles through a magnetic field. |
| 1783 | // Uses Nystroem algorithm (See Handbook Nat. Bur. of Standards, procedure 25.5.20) |
| 1784 | // in the way described in |
| 1785 | // E Lund et al 2009 JINST 4 P04001 doi:10.1088/1748-0221/4/04/P04001 |
| 1786 | // "Track parameter propagation through the application of a new adaptive Runge-Kutta-Nyström method in the ATLAS experiment" |
| 1787 | // http://inspirehep.net/search?ln=en&ln=en&p=10.1088/1748-0221/4/04/P04001&of=hb&action_search=Search&sf=earliestdate&so=d&rm=&rg=25&sc=0 |
| 1788 | // |
| 1789 | // Input parameters: |
| 1790 | // SU - plane parameters |
| 1791 | // SU[0] - direction cosines normal to surface Ex |
| 1792 | // SU[1] - ------- Ey |
| 1793 | // SU[2] - ------- Ez; Ex*Ex+Ey*Ey+Ez*Ez=1 |
| 1794 | // SU[3] - distance to surface from (0,0,0) > 0 cm |
| 1795 | // |
| 1796 | // state7 - initial parameters (coordinates(cm), direction, |
| 1797 | // charge/momentum (Gev-1) |
| 1798 | // cov and derivatives this parameters (7x7) |
| 1799 | // |
| 1800 | // X Y Z Ax Ay Az q/P |
| 1801 | // state7[0] state7[1] state7[2] state7[3] state7[4] state7[5] state7[6] |
| 1802 | // |
| 1803 | // dX/dp dY/dp dZ/dp dAx/dp dAy/dp dAz/dp d(q/P)/dp |
| 1804 | // cov[ 0] cov[ 1] cov[ 2] cov[ 3] cov[ 4] cov[ 5] cov[ 6] d()/dp1 |
| 1805 | // |
| 1806 | // cov[ 7] cov[ 8] cov[ 9] cov[10] cov[11] cov[12] cov[13] d()/dp2 |
| 1807 | // ............................................................................ d()/dpND |
| 1808 | // |
| 1809 | // Authors: R.Brun, M.Hansroul, V.Perevoztchikov (Geant3) |
| 1810 | // |
| 1811 | bool RKTrackRep::RKutta(const M1x4& SU, |
| 1812 | const DetPlane& plane, |
| 1813 | double charge, |
| 1814 | double mass, |
| 1815 | M1x7& state7, |
| 1816 | M7x7* jacobianT, |
| 1817 | M1x7* J_MMT_unprojected_lastRow, |
| 1818 | double& coveredDistance, |
| 1819 | double& flightTime, |
| 1820 | bool& checkJacProj, |
| 1821 | M7x7& noiseProjection, |
| 1822 | StepLimits& limits, |
| 1823 | bool onlyOneStep, |
| 1824 | bool calcOnlyLastRowOfJ) const { |
| 1825 | |
| 1826 | // limits, check-values, etc. Can be tuned! |
| 1827 | static const double Wmax ( 3000. ); // max. way allowed [cm] |
| 1828 | static const double AngleMax ( 6.3 ); // max. total angle change of momentum. Prevents extrapolating a curler round and round if no active plane is found. |
| 1829 | static const double Pmin ( 4.E-3 ); // minimum momentum for propagation [GeV] |
| 1830 | static const unsigned int maxNumIt ( 1000 ); // maximum number of iterations in main loop |
| 1831 | // Aux parameters |
| 1832 | M1x3& R ( *((M1x3*) &state7[0]) ); // Start coordinates [cm] (x, y, z) |
| 1833 | M1x3& A ( *((M1x3*) &state7[3]) ); // Start directions (ax, ay, az); ax^2+ay^2+az^2=1 |
| 1834 | M1x3 SA = {{0.,0.,0.}}; // Start directions derivatives dA/S |
| 1835 | double Way ( 0. ); // Sum of absolute values of all extrapolation steps [cm] |
| 1836 | double momentum ( fabs(charge/state7[6]) ); // momentum [GeV] |
| 1837 | double relMomLoss ( 0 ); // relative momentum loss in RKutta |
| 1838 | double deltaAngle ( 0. ); // total angle by which the momentum has changed during extrapolation |
| 1839 | // cppcheck-suppress unreadVariable |
| 1840 | double An(0), S(0), Sl(0), CBA(0); |
| 1841 | |
| 1842 | |
| 1843 | if (debugLvl_ > 0) { |
| 1844 | debugOut << "RKTrackRep::RKutta \n"; |
| 1845 | debugOut << "position: "; TVector3(R[0], R[1], R[2]).Print(); |
| 1846 | debugOut << "direction: "; TVector3(A[0], A[1], A[2]).Print(); |
| 1847 | debugOut << "momentum: " << momentum << " GeV\n"; |
| 1848 | debugOut << "destination: "; plane.Print(); |
| 1849 | } |
| 1850 | |
| 1851 | checkJacProj = false; |
| 1852 | |
| 1853 | // check momentum |
| 1854 | if(momentum < Pmin){ |
| 1855 | std::ostringstream sstream; |
| 1856 | sstream << "RKTrackRep::RKutta ==> momentum too low: " << momentum*1000. << " MeV"; |
| 1857 | Exception exc(sstream.str(),__LINE__1857,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1858 | exc.setFatal(); |
| 1859 | throw exc; |
| 1860 | } |
| 1861 | |
| 1862 | unsigned int counter(0); |
| 1863 | |
| 1864 | // Step estimation (signed) |
| 1865 | S = estimateStep(state7, SU, plane, charge, relMomLoss, limits); |
| 1866 | |
| 1867 | // |
| 1868 | // Main loop of Runge-Kutta method |
| 1869 | // |
| 1870 | while (fabs(S) >= MINSTEP0.001 || counter == 0) { |
| 1871 | |
| 1872 | if(++counter > maxNumIt){ |
| 1873 | Exception exc("RKTrackRep::RKutta ==> maximum number of iterations exceeded",__LINE__1873,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1874 | exc.setFatal(); |
| 1875 | throw exc; |
| 1876 | } |
| 1877 | |
| 1878 | if (debugLvl_ > 0) { |
| 1879 | debugOut << "------ RKutta main loop nr. " << counter-1 << " ------\n"; |
| 1880 | } |
| 1881 | |
| 1882 | M1x3 ABefore = {{ A[0], A[1], A[2] }}; |
| 1883 | RKPropagate(state7, jacobianT, SA, S, true, calcOnlyLastRowOfJ); // the actual Runge Kutta propagation |
| 1884 | |
| 1885 | // update paths |
| 1886 | coveredDistance += S; // add stepsize to way (signed) |
| 1887 | Way += fabs(S); |
| 1888 | |
| 1889 | double beta = 1/hypot(1, mass*state7[6]/charge); |
| 1890 | flightTime += S / beta / 29.9792458; // in ns |
| 1891 | |
| 1892 | // check way limit |
| 1893 | if(Way > Wmax){ |
| 1894 | std::ostringstream sstream; |
| 1895 | sstream << "RKTrackRep::RKutta ==> Total extrapolation length is longer than length limit : " << Way << " cm !"; |
| 1896 | Exception exc(sstream.str(),__LINE__1896,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1897 | exc.setFatal(); |
| 1898 | throw exc; |
| 1899 | } |
| 1900 | |
| 1901 | if (onlyOneStep) return(true); |
| 1902 | |
| 1903 | // if stepsize has been limited by material, break the loop and return. No linear extrapolation! |
| 1904 | if (limits.getLowestLimit().first == EStepLimitType::stp_momLoss) { |
| 1905 | if (debugLvl_ > 0) { |
| 1906 | debugOut<<" momLossExceeded -> return(true); \n"; |
| 1907 | } |
| 1908 | return(true); |
| 1909 | } |
| 1910 | |
| 1911 | // if stepsize has been limited by material boundary, break the loop and return. No linear extrapolation! |
| 1912 | if (limits.getLowestLimit().first == EStepLimitType::stp_boundary) { |
| 1913 | if (debugLvl_ > 0) { |
| 1914 | debugOut<<" at boundary -> return(true); \n"; |
| 1915 | } |
| 1916 | return(true); |
| 1917 | } |
| 1918 | |
| 1919 | |
| 1920 | // estimate Step for next loop or linear extrapolation |
| 1921 | Sl = S; // last S used |
| 1922 | limits.removeLimit(EStepLimitType::stp_fieldCurv); |
| 1923 | limits.removeLimit(EStepLimitType::stp_momLoss); |
| 1924 | limits.removeLimit(EStepLimitType::stp_boundary); |
| 1925 | limits.removeLimit(EStepLimitType::stp_plane); |
| 1926 | S = estimateStep(state7, SU, plane, charge, relMomLoss, limits); |
| 1927 | |
| 1928 | if (limits.getLowestLimit().first == EStepLimitType::stp_plane && |
| 1929 | fabs(S) < MINSTEP0.001) { |
| 1930 | if (debugLvl_ > 0) { |
| 1931 | debugOut<<" (at Plane && fabs(S) < MINSTEP) -> break and do linear extrapolation \n"; |
| 1932 | } |
| 1933 | break; |
| 1934 | } |
| 1935 | if (limits.getLowestLimit().first == EStepLimitType::stp_momLoss && |
| 1936 | fabs(S) < MINSTEP0.001) { |
| 1937 | if (debugLvl_ > 0) { |
| 1938 | debugOut<<" (momLossExceeded && fabs(S) < MINSTEP) -> return(true), no linear extrapolation; \n"; |
| 1939 | } |
| 1940 | RKSteps_.erase(RKSteps_.end()-1); |
| 1941 | --RKStepsFXStop_; |
| 1942 | return(true); // no linear extrapolation! |
| 1943 | } |
| 1944 | |
| 1945 | // check if total angle is bigger than AngleMax. Can happen if a curler should be fitted and it does not hit the active area of the next plane. |
| 1946 | double arg = ABefore[0]*A[0] + ABefore[1]*A[1] + ABefore[2]*A[2]; |
| 1947 | arg = arg > 1 ? 1 : arg; |
| 1948 | arg = arg < -1 ? -1 : arg; |
| 1949 | deltaAngle += acos(arg); |
| 1950 | if (fabs(deltaAngle) > AngleMax){ |
| 1951 | std::ostringstream sstream; |
| 1952 | sstream << "RKTrackRep::RKutta ==> Do not get to an active plane! Already extrapolated " << deltaAngle * 180 / TMath::Pi() << "°."; |
| 1953 | Exception exc(sstream.str(),__LINE__1953,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1954 | exc.setFatal(); |
| 1955 | throw exc; |
| 1956 | } |
| 1957 | |
| 1958 | // check if we went back and forth multiple times -> we don't come closer to the plane! |
| 1959 | if (counter > 3){ |
| 1960 | if (S *RKSteps_.at(counter-1).matStep_.stepSize_ < 0 && |
| 1961 | RKSteps_.at(counter-1).matStep_.stepSize_*RKSteps_.at(counter-2).matStep_.stepSize_ < 0 && |
| 1962 | RKSteps_.at(counter-2).matStep_.stepSize_*RKSteps_.at(counter-3).matStep_.stepSize_ < 0){ |
| 1963 | Exception exc("RKTrackRep::RKutta ==> Do not get closer to plane!",__LINE__1963,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 1964 | exc.setFatal(); |
| 1965 | throw exc; |
| 1966 | } |
| 1967 | } |
| 1968 | |
| 1969 | } //end of main loop |
| 1970 | |
| 1971 | |
| 1972 | // |
| 1973 | // linear extrapolation to plane |
| 1974 | // |
| 1975 | if (limits.getLowestLimit().first == EStepLimitType::stp_plane) { |
| 1976 | |
| 1977 | if (fabs(Sl) > 0.001*MINSTEP0.001){ |
| 1978 | if (debugLvl_ > 0) { |
| 1979 | debugOut << " RKutta - linear extrapolation to surface\n"; |
| 1980 | } |
| 1981 | Sl = 1./Sl; // Sl = inverted last Stepsize Sl |
| 1982 | |
| 1983 | // normalize SA |
| 1984 | SA[0]*=Sl; SA[1]*=Sl; SA[2]*=Sl; // SA/Sl = delta A / delta way; local derivative of A with respect to the length of the way |
| 1985 | |
| 1986 | // calculate A |
| 1987 | A[0] += SA[0]*S; // S = distance to surface |
| 1988 | A[1] += SA[1]*S; // A = A + S * SA*Sl |
| 1989 | A[2] += SA[2]*S; |
| 1990 | |
| 1991 | // normalize A |
| 1992 | CBA = 1./sqrt(A[0]*A[0]+A[1]*A[1]+A[2]*A[2]); // 1/|A| |
| 1993 | A[0] *= CBA; A[1] *= CBA; A[2] *= CBA; |
| 1994 | |
| 1995 | R[0] += S*(A[0]-0.5*S*SA[0]); // R = R + S*(A - 0.5*S*SA); approximation for final point on surface |
| 1996 | R[1] += S*(A[1]-0.5*S*SA[1]); |
| 1997 | R[2] += S*(A[2]-0.5*S*SA[2]); |
| 1998 | |
| 1999 | |
| 2000 | coveredDistance += S; |
| 2001 | // cppcheck-suppress unreadVariable |
| 2002 | Way += fabs(S); |
Value stored to 'Way' is never read | |
| 2003 | |
| 2004 | double beta = 1/hypot(1, mass*state7[6]/charge); |
| 2005 | flightTime += S / beta / 29.9792458; // in ns; |
| 2006 | } |
| 2007 | else if (debugLvl_ > 0) { |
| 2008 | debugOut << " RKutta - last stepsize too small -> can't do linear extrapolation! \n"; |
| 2009 | } |
| 2010 | |
| 2011 | // |
| 2012 | // Project Jacobian of extrapolation onto destination plane |
| 2013 | // |
| 2014 | if (jacobianT != nullptr) { |
| 2015 | |
| 2016 | // projected jacobianT |
| 2017 | // x x x x x x 0 |
| 2018 | // x x x x x x 0 |
| 2019 | // x x x x x x 0 |
| 2020 | // x x x x x x 0 |
| 2021 | // x x x x x x 0 |
| 2022 | // x x x x x x 0 |
| 2023 | // x x x x x x 1 |
| 2024 | |
| 2025 | if (checkJacProj && RKSteps_.size()>0){ |
| 2026 | Exception exc("RKTrackRep::Extrap ==> covariance is projected onto destination plane again",__LINE__2026,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 2027 | throw exc; |
| 2028 | } |
| 2029 | |
| 2030 | if (debugLvl_ > 0) { |
| 2031 | //debugOut << " Jacobian^T of extrapolation before Projection:\n"; |
| 2032 | //RKTools::printDim(*jacobianT, 7,7); |
| 2033 | debugOut << " Project Jacobian of extrapolation onto destination plane\n"; |
| 2034 | } |
| 2035 | An = A[0]*SU[0] + A[1]*SU[1] + A[2]*SU[2]; |
| 2036 | An = (fabs(An) > 1.E-7 ? 1./An : 0); // 1/A_normal |
| 2037 | double norm; |
| 2038 | int i=0; |
| 2039 | if (calcOnlyLastRowOfJ) |
| 2040 | i = 42; |
| 2041 | |
| 2042 | double* jacPtr = jacobianT->begin(); |
| 2043 | |
| 2044 | for(unsigned int j=42; j<49; j+=7) { |
| 2045 | (*J_MMT_unprojected_lastRow)[j-42] = jacPtr[j]; |
| 2046 | } |
| 2047 | |
| 2048 | for(; i<49; i+=7) { |
| 2049 | norm = (jacPtr[i]*SU[0] + jacPtr[i+1]*SU[1] + jacPtr[i+2]*SU[2]) * An; // dR_normal / A_normal |
| 2050 | jacPtr[i] -= norm*A [0]; jacPtr[i+1] -= norm*A [1]; jacPtr[i+2] -= norm*A [2]; |
| 2051 | jacPtr[i+3] -= norm*SA[0]; jacPtr[i+4] -= norm*SA[1]; jacPtr[i+5] -= norm*SA[2]; |
| 2052 | } |
| 2053 | checkJacProj = true; |
| 2054 | |
| 2055 | |
| 2056 | if (debugLvl_ > 0) { |
| 2057 | //debugOut << " Jacobian^T of extrapolation after Projection:\n"; |
| 2058 | //RKTools::printDim(*jacobianT, 7,7); |
| 2059 | } |
| 2060 | |
| 2061 | if (!calcOnlyLastRowOfJ) { |
| 2062 | for (int iRow = 0; iRow < 3; ++iRow) { |
| 2063 | for (int iCol = 0; iCol < 3; ++iCol) { |
| 2064 | noiseProjection[iRow*7 + iCol] = (iRow == iCol) - An * SU[iCol] * A[iRow]; |
| 2065 | noiseProjection[(iRow + 3)*7 + iCol] = - An * SU[iCol] * SA[iRow]; |
| 2066 | } |
| 2067 | } |
| 2068 | |
| 2069 | // noiseProjection will look like this: |
| 2070 | // x x x 0 0 0 0 |
| 2071 | // x x x 0 0 0 0 |
| 2072 | // x x x 0 0 0 0 |
| 2073 | // x x x 1 0 0 0 |
| 2074 | // x x x 0 1 0 0 |
| 2075 | // x x x 0 0 1 0 |
| 2076 | // 0 0 0 0 0 0 1 |
| 2077 | } |
| 2078 | |
| 2079 | } |
| 2080 | } // end of linear extrapolation to surface |
| 2081 | |
| 2082 | return(true); |
| 2083 | |
| 2084 | } |
| 2085 | |
| 2086 | |
| 2087 | double RKTrackRep::estimateStep(const M1x7& state7, |
| 2088 | const M1x4& SU, |
| 2089 | const DetPlane& plane, |
| 2090 | const double& charge, |
| 2091 | double& relMomLoss, |
| 2092 | StepLimits& limits) const { |
| 2093 | |
| 2094 | if (useCache_) { |
| 2095 | if (cachePos_ >= RKSteps_.size()) { |
| 2096 | useCache_ = false; |
| 2097 | } |
| 2098 | else { |
| 2099 | if (RKSteps_.at(cachePos_).limits_.getLowestLimit().first == EStepLimitType::stp_plane) { |
| 2100 | // we need to step exactly to the plane, so don't use the cache! |
| 2101 | useCache_ = false; |
| 2102 | RKSteps_.erase(RKSteps_.begin() + cachePos_, RKSteps_.end()); |
| 2103 | } |
| 2104 | else { |
| 2105 | if (debugLvl_ > 0) { |
| 2106 | debugOut << " RKTrackRep::estimateStep: use stepSize " << cachePos_ << " from cache: " << RKSteps_.at(cachePos_).matStep_.stepSize_ << "\n"; |
| 2107 | } |
| 2108 | //for(int n = 0; n < 1*7; ++n) RKSteps_[cachePos_].state7_[n] = state7[n]; |
| 2109 | ++RKStepsFXStop_; |
| 2110 | limits = RKSteps_.at(cachePos_).limits_; |
| 2111 | return RKSteps_.at(cachePos_++).matStep_.stepSize_; |
| 2112 | } |
| 2113 | } |
| 2114 | } |
| 2115 | |
| 2116 | limits.setLimit(EStepLimitType::stp_sMax, 25.); // max. step allowed [cm] |
| 2117 | |
| 2118 | if (debugLvl_ > 0) { |
| 2119 | debugOut << " RKTrackRep::estimateStep \n"; |
| 2120 | debugOut << " position: "; TVector3(state7[0], state7[1], state7[2]).Print(); |
| 2121 | debugOut << " direction: "; TVector3(state7[3], state7[4], state7[5]).Print(); |
| 2122 | } |
| 2123 | |
| 2124 | // calculate SL distance to surface |
| 2125 | double Dist ( SU[3] - (state7[0]*SU[0] + |
| 2126 | state7[1]*SU[1] + |
| 2127 | state7[2]*SU[2]) ); // Distance between start coordinates and surface |
| 2128 | double An ( state7[3]*SU[0] + |
| 2129 | state7[4]*SU[1] + |
| 2130 | state7[5]*SU[2] ); // An = dir * N; component of dir normal to surface |
| 2131 | |
| 2132 | double SLDist; // signed |
| 2133 | if (fabs(An) > 1.E-10) |
| 2134 | SLDist = Dist/An; |
| 2135 | else { |
| 2136 | SLDist = Dist*1.E10; |
| 2137 | if (An<0) SLDist *= -1.; |
| 2138 | } |
| 2139 | |
| 2140 | limits.setLimit(EStepLimitType::stp_plane, SLDist); |
| 2141 | limits.setStepSign(SLDist); |
| 2142 | |
| 2143 | if (debugLvl_ > 0) { |
| 2144 | debugOut << " Distance to plane: " << Dist << "\n"; |
| 2145 | debugOut << " SL distance to plane: " << SLDist << "\n"; |
| 2146 | if (limits.getStepSign()>0) |
| 2147 | debugOut << " Direction is pointing towards surface.\n"; |
| 2148 | else |
| 2149 | debugOut << " Direction is pointing away from surface.\n"; |
| 2150 | } |
| 2151 | // DONE calculate SL distance to surface |
| 2152 | |
| 2153 | // |
| 2154 | // Limit according to curvature and magnetic field inhomogenities |
| 2155 | // and improve stepsize estimation to reach plane |
| 2156 | // |
| 2157 | double fieldCurvLimit( limits.getLowestLimitSignedVal() ); // signed |
| 2158 | std::pair<double, double> distVsStep (9.E99, 9.E99); // first: smallest straight line distances to plane; second: RK steps |
| 2159 | |
| 2160 | static const unsigned int maxNumIt = 10; |
| 2161 | unsigned int counter(0); |
| 2162 | |
| 2163 | while (fabs(fieldCurvLimit) > MINSTEP0.001) { |
| 2164 | |
| 2165 | if(++counter > maxNumIt){ |
| 2166 | // if max iterations are reached, take a safe value |
| 2167 | // (in previous iteration, fieldCurvLimit has been not more than doubled) |
| 2168 | // and break. |
| 2169 | fieldCurvLimit *= 0.5; |
| 2170 | break; |
| 2171 | } |
| 2172 | |
| 2173 | M1x7 state7_temp = state7; |
| 2174 | M1x3 SA = {{0, 0, 0}}; |
| 2175 | |
| 2176 | double q ( RKPropagate(state7_temp, nullptr, SA, fieldCurvLimit, true) ); |
| 2177 | if (debugLvl_ > 0) { |
| 2178 | debugOut << " maxStepArg = " << fieldCurvLimit << "; q = " << q << " \n"; |
| 2179 | } |
| 2180 | |
| 2181 | // remember steps and resulting SL distances to plane for stepsize improvement |
| 2182 | // calculate distance to surface |
| 2183 | Dist = SU[3] - (state7_temp[0] * SU[0] + |
| 2184 | state7_temp[1] * SU[1] + |
| 2185 | state7_temp[2] * SU[2]); // Distance between position and surface |
| 2186 | |
| 2187 | An = state7_temp[3] * SU[0] + |
| 2188 | state7_temp[4] * SU[1] + |
| 2189 | state7_temp[5] * SU[2]; // An = dir * N; component of dir normal to surface |
| 2190 | |
| 2191 | if (fabs(Dist/An) < fabs(distVsStep.first)) { |
| 2192 | distVsStep.first = Dist/An; |
| 2193 | distVsStep.second = fieldCurvLimit; |
| 2194 | } |
| 2195 | |
| 2196 | // resize limit according to q never grow step size more than |
| 2197 | // two-fold to avoid infinite grow-shrink loops with strongly |
| 2198 | // inhomogeneous fields. |
| 2199 | if (q>2) { |
| 2200 | fieldCurvLimit *= 2; |
| 2201 | break; |
| 2202 | } |
| 2203 | |
| 2204 | fieldCurvLimit *= q * 0.95; |
| 2205 | |
| 2206 | if (fabs(q-1) < 0.25 || // good enough! |
| 2207 | fabs(fieldCurvLimit) > limits.getLowestLimitVal()) // other limits are lower! |
| 2208 | break; |
| 2209 | } |
| 2210 | if (fabs(fieldCurvLimit) < MINSTEP0.001) |
| 2211 | limits.setLimit(EStepLimitType::stp_fieldCurv, MINSTEP0.001); |
| 2212 | else |
| 2213 | limits.setLimit(EStepLimitType::stp_fieldCurv, fieldCurvLimit); |
| 2214 | |
| 2215 | double stepToPlane(limits.getLimitSigned(EStepLimitType::stp_plane)); |
| 2216 | if (fabs(distVsStep.first) < 8.E99) { |
| 2217 | stepToPlane = distVsStep.first + distVsStep.second; |
| 2218 | } |
| 2219 | limits.setLimit(EStepLimitType::stp_plane, stepToPlane); |
| 2220 | |
| 2221 | |
| 2222 | // |
| 2223 | // Select direction |
| 2224 | // |
| 2225 | // auto select |
| 2226 | if (propDir_ == 0 || !plane.isFinite()){ |
| 2227 | if (debugLvl_ > 0) { |
| 2228 | debugOut << " auto select direction"; |
| 2229 | if (!plane.isFinite()) debugOut << ", plane is not finite"; |
| 2230 | debugOut << ".\n"; |
| 2231 | } |
| 2232 | } |
| 2233 | // see if straight line approximation is ok |
| 2234 | else if ( limits.getLimit(EStepLimitType::stp_plane) < 0.2*limits.getLimit(EStepLimitType::stp_fieldCurv) ){ |
| 2235 | if (debugLvl_ > 0) { |
| 2236 | debugOut << " straight line approximation is fine.\n"; |
| 2237 | } |
| 2238 | |
| 2239 | // if direction is pointing to active part of surface |
| 2240 | if( plane.isInActive(state7[0], state7[1], state7[2], state7[3], state7[4], state7[5]) ) { |
| 2241 | if (debugLvl_ > 0) { |
| 2242 | debugOut << " direction is pointing to active part of surface. \n"; |
| 2243 | } |
| 2244 | } |
| 2245 | // if we are near the plane, but not pointing to the active area, make a big step! |
| 2246 | else { |
| 2247 | limits.removeLimit(EStepLimitType::stp_plane); |
| 2248 | limits.setStepSign(propDir_); |
| 2249 | if (debugLvl_ > 0) { |
| 2250 | debugOut << " we are near the plane, but not pointing to the active area. make a big step! \n"; |
| 2251 | } |
| 2252 | } |
| 2253 | } |
| 2254 | // propDir_ is set and we are not pointing to an active part of a plane -> propDir_ decides! |
| 2255 | else { |
| 2256 | if (limits.getStepSign() * propDir_ < 0){ |
| 2257 | limits.removeLimit(EStepLimitType::stp_plane); |
| 2258 | limits.setStepSign(propDir_); |
| 2259 | if (debugLvl_ > 0) { |
| 2260 | debugOut << " invert Step according to propDir_ and make a big step. \n"; |
| 2261 | } |
| 2262 | } |
| 2263 | } |
| 2264 | |
| 2265 | |
| 2266 | // call stepper and reduce stepsize if step not too small |
| 2267 | static const RKStep defaultRKStep; |
| 2268 | RKSteps_.push_back( defaultRKStep ); |
| 2269 | std::vector<RKStep>::iterator lastStep = RKSteps_.end() - 1; |
| 2270 | lastStep->state7_ = state7; |
| 2271 | ++RKStepsFXStop_; |
| 2272 | |
| 2273 | if(limits.getLowestLimitVal() > MINSTEP0.001){ // only call stepper if step estimation big enough |
| 2274 | M1x7 state7_temp = {{ state7[0], state7[1], state7[2], state7[3], state7[4], state7[5], state7[6] }}; |
| 2275 | |
| 2276 | MaterialEffects::getInstance()->stepper(this, |
| 2277 | state7_temp, |
| 2278 | charge/state7[6], // |p| |
| 2279 | relMomLoss, |
| 2280 | pdgCode_, |
| 2281 | lastStep->matStep_.material_, |
| 2282 | limits, |
| 2283 | true); |
| 2284 | } else { //assume material has not changed |
| 2285 | if (RKSteps_.size()>1) { |
| 2286 | lastStep->matStep_.material_ = (lastStep - 1)->matStep_.material_; |
| 2287 | } |
| 2288 | } |
| 2289 | |
| 2290 | if (debugLvl_ > 0) { |
| 2291 | debugOut << " final limits:\n"; |
| 2292 | limits.Print(); |
| 2293 | } |
| 2294 | |
| 2295 | double finalStep = limits.getLowestLimitSignedVal(); |
| 2296 | |
| 2297 | lastStep->matStep_.stepSize_ = finalStep; |
| 2298 | lastStep->limits_ = limits; |
| 2299 | |
| 2300 | if (debugLvl_ > 0) { |
| 2301 | debugOut << " --> Step to be used: " << finalStep << "\n"; |
| 2302 | } |
| 2303 | |
| 2304 | return finalStep; |
| 2305 | |
| 2306 | } |
| 2307 | |
| 2308 | |
| 2309 | TVector3 RKTrackRep::pocaOnLine(const TVector3& linePoint, const TVector3& lineDirection, const TVector3& point) const { |
| 2310 | |
| 2311 | TVector3 retVal(lineDirection); |
| 2312 | |
| 2313 | double t = 1./(retVal.Mag2()) * ((point*retVal) - (linePoint*retVal)); |
| 2314 | retVal *= t; |
| 2315 | retVal += linePoint; |
| 2316 | return retVal; // = linePoint + t*lineDirection |
| 2317 | |
| 2318 | } |
| 2319 | |
| 2320 | |
| 2321 | double RKTrackRep::Extrap(const DetPlane& startPlane, |
| 2322 | const DetPlane& destPlane, |
| 2323 | double charge, |
| 2324 | double mass, |
| 2325 | bool& isAtBoundary, |
| 2326 | M1x7& state7, |
| 2327 | double& flightTime, |
| 2328 | bool fillExtrapSteps, |
| 2329 | TMatrixDSym* cov, // 5D |
| 2330 | bool onlyOneStep, |
| 2331 | bool stopAtBoundary, |
| 2332 | double maxStep) const |
| 2333 | { |
| 2334 | |
| 2335 | static const unsigned int maxNumIt(500); |
| 2336 | unsigned int numIt(0); |
| 2337 | |
| 2338 | double coveredDistance(0.); |
| 2339 | // cppcheck-suppress unreadVariable |
| 2340 | double dqop(0.); |
| 2341 | |
| 2342 | const TVector3 W(destPlane.getNormal()); |
| 2343 | M1x4 SU = {{W.X(), W.Y(), W.Z(), destPlane.distance(0., 0., 0.)}}; |
| 2344 | |
| 2345 | // make SU vector point away from origin |
| 2346 | if (W*destPlane.getO() < 0) { |
| 2347 | SU[0] *= -1; |
| 2348 | SU[1] *= -1; |
| 2349 | SU[2] *= -1; |
| 2350 | } |
| 2351 | |
| 2352 | |
| 2353 | M1x7 startState7 = state7; |
| 2354 | |
| 2355 | while(true){ |
| 2356 | |
| 2357 | if (debugLvl_ > 0) { |
| 2358 | debugOut << "\n============ RKTrackRep::Extrap loop nr. " << numIt << " ============\n"; |
| 2359 | debugOut << "Start plane: "; startPlane.Print(); |
| 2360 | debugOut << "fillExtrapSteps " << fillExtrapSteps << "\n"; |
| 2361 | } |
| 2362 | |
| 2363 | if(++numIt > maxNumIt){ |
| 2364 | Exception exc("RKTrackRep::Extrap ==> maximum number of iterations exceeded",__LINE__2364,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 2365 | exc.setFatal(); |
| 2366 | throw exc; |
| 2367 | } |
| 2368 | |
| 2369 | // initialize jacobianT with unit matrix |
| 2370 | for(int i = 0; i < 7*7; ++i) J_MMT_[i] = 0; |
| 2371 | for(int i=0; i<7; ++i) J_MMT_[8*i] = 1.; |
| 2372 | |
| 2373 | M7x7* noise = nullptr; |
| 2374 | isAtBoundary = false; |
| 2375 | |
| 2376 | // propagation |
| 2377 | bool checkJacProj = false; |
| 2378 | limits_.reset(); |
| 2379 | limits_.setLimit(EStepLimitType::stp_sMaxArg, maxStep-fabs(coveredDistance)); |
| 2380 | |
| 2381 | M1x7 J_MMT_unprojected_lastRow = {{0, 0, 0, 0, 0, 0, 1}}; |
| 2382 | |
| 2383 | if( ! RKutta(SU, destPlane, charge, mass, state7, &J_MMT_, &J_MMT_unprojected_lastRow, |
| 2384 | coveredDistance, flightTime, checkJacProj, noiseProjection_, |
| 2385 | limits_, onlyOneStep, !fillExtrapSteps) ) { |
| 2386 | Exception exc("RKTrackRep::Extrap ==> Runge Kutta propagation failed",__LINE__2386,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 2387 | exc.setFatal(); |
| 2388 | throw exc; |
| 2389 | } |
| 2390 | |
| 2391 | bool atPlane(limits_.getLowestLimit().first == EStepLimitType::stp_plane); |
| 2392 | if (limits_.getLowestLimit().first == EStepLimitType::stp_boundary) |
| 2393 | isAtBoundary = true; |
| 2394 | |
| 2395 | |
| 2396 | if (debugLvl_ > 0) { |
| 2397 | debugOut<<"RKSteps \n"; |
| 2398 | for (std::vector<RKStep>::iterator it = RKSteps_.begin(); it != RKSteps_.end(); ++it){ |
| 2399 | debugOut << "stepSize = " << it->matStep_.stepSize_ << "\t"; |
| 2400 | it->matStep_.material_.Print(); |
| 2401 | } |
| 2402 | debugOut<<"\n"; |
| 2403 | } |
| 2404 | |
| 2405 | |
| 2406 | |
| 2407 | // call MatFX |
| 2408 | if(fillExtrapSteps) { |
| 2409 | noise = &noiseArray_; |
| 2410 | for(int i = 0; i < 7*7; ++i) noiseArray_[i] = 0; // set noiseArray_ to 0 |
| 2411 | } |
| 2412 | |
| 2413 | unsigned int nPoints(RKStepsFXStop_ - RKStepsFXStart_); |
| 2414 | if (/*!fNoMaterial &&*/ nPoints>0){ |
| 2415 | // momLoss has a sign - negative loss means momentum gain |
| 2416 | double momLoss = MaterialEffects::getInstance()->effects(RKSteps_, |
| 2417 | RKStepsFXStart_, |
| 2418 | RKStepsFXStop_, |
| 2419 | fabs(charge/state7[6]), // momentum |
| 2420 | pdgCode_, |
| 2421 | noise); |
| 2422 | |
| 2423 | RKStepsFXStart_ = RKStepsFXStop_; |
| 2424 | |
| 2425 | if (debugLvl_ > 0) { |
| 2426 | debugOut << "momLoss: " << momLoss << " GeV; relative: " << momLoss/fabs(charge/state7[6]) |
| 2427 | << "; coveredDistance = " << coveredDistance << "\n"; |
| 2428 | if (debugLvl_ > 1 && noise != nullptr) { |
| 2429 | debugOut << "7D noise: \n"; |
| 2430 | RKTools::printDim(noise->begin(), 7, 7); |
| 2431 | } |
| 2432 | } |
| 2433 | |
| 2434 | // do momLoss only for defined 1/momentum .ne.0 |
| 2435 | if(fabs(state7[6])>1.E-10) { |
| 2436 | |
| 2437 | if (debugLvl_ > 0) { |
| 2438 | debugOut << "correct state7 with dx/dqop, dy/dqop ...\n"; |
| 2439 | } |
| 2440 | |
| 2441 | dqop = charge/(fabs(charge/state7[6])-momLoss) - state7[6]; |
| 2442 | |
| 2443 | // Correct coveredDistance and flightTime and momLoss if checkJacProj == true |
| 2444 | // The idea is to calculate the state correction (based on the mometum loss) twice: |
| 2445 | // Once with the unprojected Jacobian (which preserves coveredDistance), |
| 2446 | // and once with the projected Jacobian (which is constrained to the plane and does NOT preserve coveredDistance). |
| 2447 | // The difference of these two corrections can then be used to calculate a correction factor. |
| 2448 | if (checkJacProj && fabs(coveredDistance) > MINSTEP0.001) { |
| 2449 | M1x3 state7_correction_unprojected = {{0, 0, 0}}; |
| 2450 | for (unsigned int i=0; i<3; ++i) { |
| 2451 | state7_correction_unprojected[i] = 0.5 * dqop * J_MMT_unprojected_lastRow[i]; |
| 2452 | //debugOut << "J_MMT_unprojected_lastRow[i] " << J_MMT_unprojected_lastRow[i] << "\n"; |
| 2453 | //debugOut << "state7_correction_unprojected[i] " << state7_correction_unprojected[i] << "\n"; |
| 2454 | } |
| 2455 | |
| 2456 | M1x3 state7_correction_projected = {{0, 0, 0}}; |
| 2457 | for (unsigned int i=0; i<3; ++i) { |
| 2458 | state7_correction_projected[i] = 0.5 * dqop * J_MMT_[6*7 + i]; |
| 2459 | //debugOut << "J_MMT_[6*7 + i] " << J_MMT_[6*7 + i] << "\n"; |
| 2460 | //debugOut << "state7_correction_projected[i] " << state7_correction_projected[i] << "\n"; |
| 2461 | } |
| 2462 | |
| 2463 | // delta distance |
| 2464 | M1x3 delta_state = {{0, 0, 0}}; |
| 2465 | for (unsigned int i=0; i<3; ++i) { |
| 2466 | delta_state[i] = state7_correction_unprojected[i] - state7_correction_projected[i]; |
| 2467 | } |
| 2468 | |
| 2469 | double Dist( sqrt(delta_state[0]*delta_state[0] |
| 2470 | + delta_state[1]*delta_state[1] |
| 2471 | + delta_state[2]*delta_state[2] ) ); |
| 2472 | |
| 2473 | // sign: delta * a |
| 2474 | if (delta_state[0]*state7[3] + delta_state[1]*state7[4] + delta_state[2]*state7[5] > 0) |
| 2475 | Dist *= -1.; |
| 2476 | |
| 2477 | double correctionFactor( 1. + Dist / coveredDistance ); |
| 2478 | flightTime *= correctionFactor; |
| 2479 | momLoss *= correctionFactor; |
| 2480 | coveredDistance = coveredDistance + Dist; |
| 2481 | |
| 2482 | if (debugLvl_ > 0) { |
| 2483 | debugOut << "correctionFactor-1 = " << correctionFactor-1. << "; Dist = " << Dist << "\n"; |
| 2484 | debugOut << "corrected momLoss: " << momLoss << " GeV; relative: " << momLoss/fabs(charge/state7[6]) |
| 2485 | << "; corrected coveredDistance = " << coveredDistance << "\n"; |
| 2486 | } |
| 2487 | } |
| 2488 | |
| 2489 | // correct state7 with dx/dqop, dy/dqop ... Greatly improves extrapolation accuracy! |
| 2490 | double qop( charge/(fabs(charge/state7[6])-momLoss) ); |
| 2491 | dqop = qop - state7[6]; |
| 2492 | state7[6] = qop; |
| 2493 | |
| 2494 | for (unsigned int i=0; i<6; ++i) { |
| 2495 | state7[i] += 0.5 * dqop * J_MMT_[6*7 + i]; |
| 2496 | } |
| 2497 | // normalize direction, just to make sure |
| 2498 | double norm( 1. / sqrt(state7[3]*state7[3] + state7[4]*state7[4] + state7[5]*state7[5]) ); |
| 2499 | for (unsigned int i=3; i<6; ++i) |
| 2500 | state7[i] *= norm; |
| 2501 | |
| 2502 | } |
| 2503 | } // finished MatFX |
| 2504 | |
| 2505 | |
| 2506 | // fill ExtrapSteps_ |
| 2507 | if (fillExtrapSteps) { |
| 2508 | static const ExtrapStep defaultExtrapStep; |
| 2509 | ExtrapSteps_.push_back(defaultExtrapStep); |
| 2510 | std::vector<ExtrapStep>::iterator lastStep = ExtrapSteps_.end() - 1; |
| 2511 | |
| 2512 | // Store Jacobian of this step for final calculation. |
| 2513 | lastStep->jac7_ = J_MMT_; |
| 2514 | |
| 2515 | if( checkJacProj == true ){ |
| 2516 | //project the noise onto the destPlane |
| 2517 | RKTools::Np_N_NpT(noiseProjection_, noiseArray_); |
| 2518 | |
| 2519 | if (debugLvl_ > 1) { |
| 2520 | debugOut << "7D noise projected onto plane: \n"; |
| 2521 | RKTools::printDim(noiseArray_.begin(), 7, 7); |
| 2522 | } |
| 2523 | } |
| 2524 | |
| 2525 | // Store this step's noise for final calculation. |
| 2526 | lastStep->noise7_ = noiseArray_; |
| 2527 | |
| 2528 | if (debugLvl_ > 2) { |
| 2529 | debugOut<<"ExtrapSteps \n"; |
| 2530 | for (std::vector<ExtrapStep>::iterator it = ExtrapSteps_.begin(); it != ExtrapSteps_.end(); ++it){ |
| 2531 | debugOut << "7D Jacobian: "; RKTools::printDim((it->jac7_.begin()), 5,5); |
| 2532 | debugOut << "7D noise: "; RKTools::printDim((it->noise7_.begin()), 5,5); |
| 2533 | } |
| 2534 | debugOut<<"\n"; |
| 2535 | } |
| 2536 | } |
| 2537 | |
| 2538 | |
| 2539 | |
| 2540 | // check if at boundary |
| 2541 | if (stopAtBoundary and isAtBoundary) { |
| 2542 | if (debugLvl_ > 0) { |
| 2543 | debugOut << "stopAtBoundary -> break; \n "; |
| 2544 | } |
| 2545 | break; |
| 2546 | } |
| 2547 | |
| 2548 | if (onlyOneStep) { |
| 2549 | if (debugLvl_ > 0) { |
| 2550 | debugOut << "onlyOneStep -> break; \n "; |
| 2551 | } |
| 2552 | break; |
| 2553 | } |
| 2554 | |
| 2555 | //break if we arrived at destPlane |
| 2556 | if(atPlane) { |
| 2557 | if (debugLvl_ > 0) { |
| 2558 | debugOut << "arrived at destPlane with a distance of " << destPlane.distance(state7[0], state7[1], state7[2]) << " cm left. "; |
| 2559 | if (destPlane.isInActive(state7[0], state7[1], state7[2], state7[3], state7[4], state7[5])) |
| 2560 | debugOut << "In active area of destPlane. \n"; |
| 2561 | else |
| 2562 | debugOut << "NOT in active area of plane. \n"; |
| 2563 | |
| 2564 | debugOut << " position: "; TVector3(state7[0], state7[1], state7[2]).Print(); |
| 2565 | debugOut << " direction: "; TVector3(state7[3], state7[4], state7[5]).Print(); |
| 2566 | } |
| 2567 | break; |
| 2568 | } |
| 2569 | |
| 2570 | } |
| 2571 | |
| 2572 | if (fillExtrapSteps) { |
| 2573 | // propagate cov and add noise |
| 2574 | calcForwardJacobianAndNoise(startState7, startPlane, state7, destPlane); |
| 2575 | |
| 2576 | if (cov != nullptr) { |
| 2577 | cov->Similarity(fJacobian_); |
| 2578 | *cov += fNoise_; |
| 2579 | } |
| 2580 | |
| 2581 | if (debugLvl_ > 0) { |
| 2582 | if (cov != nullptr) { |
| 2583 | debugOut << "final covariance matrix after Extrap: "; cov->Print(); |
| 2584 | } |
| 2585 | } |
| 2586 | } |
| 2587 | |
| 2588 | return coveredDistance; |
| 2589 | } |
| 2590 | |
| 2591 | |
| 2592 | void RKTrackRep::checkCache(const StateOnPlane& state, const SharedPlanePtr* plane) const { |
| 2593 | |
| 2594 | if (state.getRep() != this){ |
| 2595 | Exception exc("RKTrackRep::checkCache ==> state is defined wrt. another TrackRep",__LINE__2595,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 2596 | exc.setFatal(); |
| 2597 | throw exc; |
| 2598 | } |
| 2599 | |
| 2600 | if (dynamic_cast<const MeasurementOnPlane*>(&state) != nullptr) { |
| 2601 | Exception exc("RKTrackRep::checkCache - cannot extrapolate MeasurementOnPlane",__LINE__2601,__FILE__"genfit2/code2/trackReps/src/RKTrackRep.cc"); |
| 2602 | exc.setFatal(); |
| 2603 | throw exc; |
| 2604 | } |
| 2605 | |
| 2606 | cachePos_ = 0; |
| 2607 | RKStepsFXStart_ = 0; |
| 2608 | RKStepsFXStop_ = 0; |
| 2609 | ExtrapSteps_.clear(); |
| 2610 | initArrays(); |
| 2611 | |
| 2612 | |
| 2613 | if (plane && |
| 2614 | lastStartState_.getPlane() && |
| 2615 | lastEndState_.getPlane() && |
| 2616 | state.getPlane() == lastStartState_.getPlane() && |
| 2617 | state.getState() == lastStartState_.getState() && |
| 2618 | (*plane)->distance(getPos(lastEndState_)) <= MINSTEP0.001) { |
| 2619 | useCache_ = true; |
| 2620 | |
| 2621 | // clean up cache. Only use steps with same sign. |
| 2622 | double firstStep(0); |
| 2623 | for (unsigned int i=0; i<RKSteps_.size(); ++i) { |
| 2624 | if (i == 0) { |
| 2625 | firstStep = RKSteps_.at(0).matStep_.stepSize_; |
| 2626 | continue; |
| 2627 | } |
| 2628 | if (RKSteps_.at(i).matStep_.stepSize_ * firstStep < 0) { |
| 2629 | if (RKSteps_.at(i-1).matStep_.material_ == RKSteps_.at(i).matStep_.material_) { |
| 2630 | RKSteps_.at(i-1).matStep_.stepSize_ += RKSteps_.at(i).matStep_.stepSize_; |
| 2631 | } |
| 2632 | RKSteps_.erase(RKSteps_.begin()+i, RKSteps_.end()); |
| 2633 | } |
| 2634 | } |
| 2635 | |
| 2636 | if (debugLvl_ > 0) { |
| 2637 | debugOut << "RKTrackRep::checkCache: use cached material and step values.\n"; |
| 2638 | } |
| 2639 | } |
| 2640 | else { |
| 2641 | |
| 2642 | if (debugLvl_ > 0) { |
| 2643 | debugOut << "RKTrackRep::checkCache: can NOT use cached material and step values.\n"; |
| 2644 | |
| 2645 | if (plane != nullptr) { |
| 2646 | if (state.getPlane() != lastStartState_.getPlane()) { |
| 2647 | debugOut << "state.getPlane() != lastStartState_.getPlane()\n"; |
| 2648 | } |
| 2649 | else { |
| 2650 | if (! (state.getState() == lastStartState_.getState())) { |
| 2651 | debugOut << "state.getState() != lastStartState_.getState()\n"; |
| 2652 | } |
| 2653 | else if (lastEndState_.getPlane().get() != nullptr) { |
| 2654 | debugOut << "distance " << (*plane)->distance(getPos(lastEndState_)) << "\n"; |
| 2655 | } |
| 2656 | } |
| 2657 | } |
| 2658 | } |
| 2659 | |
| 2660 | useCache_ = false; |
| 2661 | RKSteps_.clear(); |
| 2662 | |
| 2663 | lastStartState_.setStatePlane(state.getState(), state.getPlane()); |
| 2664 | } |
| 2665 | } |
| 2666 | |
| 2667 | |
| 2668 | double RKTrackRep::momMag(const M1x7& state7) const { |
| 2669 | // FIXME given this interface this function cannot work for charge =/= +-1 |
| 2670 | return fabs(1/state7[6]); |
| 2671 | } |
| 2672 | |
| 2673 | |
| 2674 | bool RKTrackRep::isSameType(const AbsTrackRep* other) { |
| 2675 | if (dynamic_cast<const RKTrackRep*>(other) == nullptr) |
| 2676 | return false; |
| 2677 | |
| 2678 | return true; |
| 2679 | } |
| 2680 | |
| 2681 | |
| 2682 | bool RKTrackRep::isSame(const AbsTrackRep* other) { |
| 2683 | if (getPDG() != other->getPDG()) |
| 2684 | return false; |
| 2685 | |
| 2686 | return isSameType(other); |
| 2687 | } |
| 2688 | |
| 2689 | |
| 2690 | void RKTrackRep::Streamer(TBuffer &R__b) |
| 2691 | { |
| 2692 | // Stream an object of class genfit::RKTrackRep. |
| 2693 | |
| 2694 | //This works around a msvc bug and should be harmless on other platforms |
| 2695 | typedef ::genfit::RKTrackRep thisClass; |
| 2696 | UInt_t R__s, R__c; |
| 2697 | if (R__b.IsReading()) { |
| 2698 | ::genfit::AbsTrackRep::Streamer(R__b); |
| 2699 | Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { } |
| 2700 | R__b.CheckByteCount(R__s, R__c, thisClass::IsA()); |
| 2701 | lastStartState_.setRep(this); |
| 2702 | lastEndState_.setRep(this); |
| 2703 | } else { |
| 2704 | ::genfit::AbsTrackRep::Streamer(R__b); |
| 2705 | R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE); |
| 2706 | R__b.SetByteCount(R__c, kTRUE); |
| 2707 | } |
| 2708 | } |
| 2709 | |
| 2710 | } /* End of namespace genfit */ |