19#include <TGeoGlobalMagField.h>
20#include <TGeoManager.h>
21#include <TGeoMaterial.h>
35bool TrackExtrap::sExtrapV2 =
false;
36double TrackExtrap::sSimpleBValue = 0.;
37bool TrackExtrap::sFieldON =
false;
38std::size_t TrackExtrap::sNCallExtrapToZCov = 0;
39std::size_t TrackExtrap::sNCallField = 0;
46 const double x[3] = {50., 50., SSimpleBPosition};
47 double b[3] = {0., 0., 0.};
48 TGeoGlobalMagField::Instance()->Field(
x,
b);
50 sFieldON = (TMath::Abs(sSimpleBValue) > 1.e-10) ?
true :
false;
51 LOG(info) <<
"Track extrapolation with magnetic field " << (sFieldON ?
"ON" :
"OFF");
61 const double correctionFactor = 1.1;
62 if (bendingMomentum == 0.) {
65 return correctionFactor * (-0.0003 * sSimpleBValue * SSimpleBLength * SSimpleBPosition / bendingMomentum);
75 const double correctionFactor = 1.1;
76 if (impactParam == 0.) {
80 return correctionFactor * (-0.0003 * sSimpleBValue * SSimpleBLength * SSimpleBPosition / impactParam);
82 return SMostProbBendingMomentum;
92 if (trackParam.
getZ() == zEnd) {
97 double dZ = zEnd - trackParam.
getZ();
100 trackParam.
setZ(zEnd);
109 if (trackParam.
getZ() == zEnd) {
115 LOG(warning) <<
"Covariance matrix does not exist";
122 double dZ = zEnd - trackParam.
getZ();
125 trackParam.
setZ(zEnd);
128 TMatrixD jacob(5, 5);
134 TMatrixD tmp(trackParam.
getCovariances(), TMatrixD::kMultTranspose, jacob);
135 TMatrixD tmp2(jacob, TMatrixD::kMult, tmp);
139 if (updatePropagator) {
152 }
else if (sExtrapV2) {
153 return extrapToZRungekuttaV2(trackParam, zEnd);
155 return extrapToZRungekutta(trackParam, zEnd);
165 ++sNCallExtrapToZCov;
167 if (trackParam.
getZ() == zEnd) {
178 LOG(warning) <<
"Covariance matrix does not exist";
186 double zBegin = trackParamSave.
getZ();
201 TMatrixD jacob(5, 5);
203 TMatrixD dParam(5, 1);
204 double direction[5] = {-1., -1., 1., 1., -1.};
205 for (
int i = 0;
i < 5;
i++) {
207 if (kParamCov(
i,
i) <= 0.) {
212 for (
int j = 0;
j < 5;
j++) {
214 dParam(
j, 0) = TMath::Sqrt(kParamCov(
i,
i));
215 dParam(
j, 0) *= TMath::Sign(1., direction[
j] * paramSave(
j, 0));
224 trackParamSave.
setZ(zBegin);
228 LOG(warning) <<
"Bad covariance matrix";
233 TMatrixD jacobji(trackParamSave.
getParameters(), TMatrixD::kMinus, extrapParam);
234 jacobji *= 1. / dParam(
i, 0);
235 jacob.SetSub(0,
i, jacobji);
239 TMatrixD tmp(kParamCov, TMatrixD::kMultTranspose, jacob);
240 TMatrixD tmp2(jacob, TMatrixD::kMult, tmp);
244 if (updatePropagator) {
257 if (trackParam.
getZ() == SMIDZ) {
262 if (trackParam.
getZ() < SMuonFilterZBeg) {
263 LOG(warning) <<
"The track already passed the beginning of the muon filter";
271 addMCSEffect(trackParam, -SMuonFilterThickness, SMuonFilterX0);
279 double errXVtx,
double errYVtx,
bool correctForMCS,
bool correctForEnergyLoss,
280 double xUpstream,
double yUpstream, std::optional<double> zUpstream)
292 if (trackParam.
getZ() == zVtx) {
297 if (zVtx < SAbsZBeg) {
298 if (zVtx < SAbsZEnd) {
299 LOG(warning) <<
"Ending Z (" << zVtx <<
") downstream the front absorber (zAbsorberEnd = " << SAbsZEnd <<
")";
302 LOG(warning) <<
"Ending Z (" << zVtx <<
") inside the front absorber (" << SAbsZBeg <<
", " << SAbsZEnd <<
")";
309 if (!correctForMCS && zUpstream && *zUpstream < SAbsZBeg - 0.01) {
310 if (*zUpstream < SAbsZEnd) {
311 LOG(warning) <<
"Upstream Z (" << *zUpstream <<
") downstream the front absorber (zAbsorberEnd = " << SAbsZEnd <<
")";
314 LOG(warning) <<
"Upstream Z (" << *zUpstream <<
") inside the front absorber (" << SAbsZBeg <<
", " << SAbsZEnd <<
")";
320 if (trackParam.
getZ() > SAbsZEnd) {
321 if (trackParam.
getZ() > zVtx) {
322 LOG(warning) <<
"Starting Z (" << trackParam.
getZ() <<
") upstream the vertex (zVtx = " << zVtx <<
")";
324 }
else if (trackParam.
getZ() > SAbsZBeg) {
325 LOG(warning) <<
"Starting Z (" << trackParam.
getZ() <<
") upstream the front absorber (zAbsorberBegin = " << SAbsZBeg <<
")";
328 LOG(warning) <<
"Starting Z (" << trackParam.
getZ() <<
") inside the front absorber (" << SAbsZBeg <<
", " << SAbsZEnd <<
")";
341 double trackXYZIn[3] = {0., 0., 0.};
343 trackXYZIn[2] = SAbsZBeg;
344 trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
345 trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
346 }
else if (zUpstream) {
347 trackXYZIn[2] = SAbsZBeg;
348 trackXYZIn[0] = trackXYZOut[0] + (xUpstream - trackXYZOut[0]) / (*zUpstream - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
349 trackXYZIn[1] = trackXYZOut[1] + (yUpstream - trackXYZOut[1]) / (*zUpstream - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
351 TrackParam trackParamIn(trackParam);
352 if (!
extrapToZ(trackParamIn, SAbsZBeg)) {
355 trackXYZIn[0] = trackParamIn.getNonBendingCoor();
356 trackXYZIn[1] = trackParamIn.getBendingCoor();
357 trackXYZIn[2] = trackParamIn.getZ();
359 double pTot = trackParam.
p();
360 double pathLength(0.), f0(0.),
f1(0.),
f2(0.), meanRho(0.), totalELoss(0.), sigmaELoss2(0.);
361 if (!getAbsorberCorrectionParam(trackXYZIn, trackXYZOut, pTot, pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2)) {
367 if (correctForEnergyLoss) {
369 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
370 if (!correctMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, trackXYZIn[2], pathLength, f0, f1, f2)) {
373 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
376 if (!correctMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, trackXYZIn[2], pathLength, f0, f1, f2)) {
381 if (correctForEnergyLoss) {
383 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
384 addMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2);
388 correctELossEffectInAbsorber(trackParam, 0.5 * totalELoss, 0.5 * sigmaELoss2);
394 addMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2);
405bool TrackExtrap::getAbsorberCorrectionParam(
double trackXYZIn[3],
double trackXYZOut[3],
double pTotal,
406 double& pathLength,
double& f0,
double& f1,
double& f2,
407 double& meanRho,
double& totalELoss,
double& sigmaELoss2)
421 LOG(warning) <<
"geometry is missing";
426 pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0]) * (trackXYZOut[0] - trackXYZIn[0]) +
427 (trackXYZOut[1] - trackXYZIn[1]) * (trackXYZOut[1] - trackXYZIn[1]) +
428 (trackXYZOut[2] - trackXYZIn[2]) * (trackXYZOut[2] - trackXYZIn[2]));
429 if (pathLength < TGeoShape::Tolerance()) {
430 LOG(warning) <<
"path length is too small";
433 double b[3] = {(trackXYZOut[0] - trackXYZIn[0]) / pathLength, (trackXYZOut[1] - trackXYZIn[1]) / pathLength, (trackXYZOut[2] - trackXYZIn[2]) / pathLength};
434 TGeoNode* currentnode = gGeoManager->InitTrack(trackXYZIn,
b);
436 LOG(warning) <<
"starting point out of geometry";
441 f0 =
f1 =
f2 = meanRho = totalELoss = 0.;
442 double sigmaELoss(0.);
443 double zB = trackXYZIn[2];
444 double remainingPathLength = pathLength;
448 TGeoMaterial* material = currentnode->GetVolume()->GetMedium()->GetMaterial();
449 double rho = material->GetDensity();
450 double x0 = material->GetRadLen();
451 double atomicZ = material->GetZ();
452 double atomicZoverA(0.);
453 if (material->IsMixture()) {
454 TGeoMixture* mixture =
static_cast<TGeoMixture*
>(material);
456 for (
int iel = 0; iel < mixture->GetNelements(); ++iel) {
457 sum += mixture->GetWmixt()[iel];
458 atomicZoverA += mixture->GetWmixt()[iel] * mixture->GetZmixt()[iel] / mixture->GetAmixt()[iel];
462 atomicZoverA = atomicZ / material->GetA();
466 gGeoManager->FindNextBoundary(remainingPathLength);
467 double localPathLength = gGeoManager->GetStep() + 1.e-6;
469 if (localPathLength >= remainingPathLength) {
470 localPathLength = remainingPathLength;
472 currentnode = gGeoManager->Step();
474 LOG(warning) <<
"navigation failed";
477 if (!gGeoManager->IsEntering()) {
479 gGeoManager->SetStep(0.001);
480 currentnode = gGeoManager->Step();
481 if (!gGeoManager->IsEntering() || !currentnode) {
482 LOG(warning) <<
"navigation failed";
485 localPathLength += 0.001;
490 double zE =
b[2] * localPathLength + zB;
491 double dzB = zB - trackXYZIn[2];
492 double dzE = zE - trackXYZIn[2];
493 f0 += localPathLength /
x0;
494 f1 += (dzE * dzE - dzB * dzB) /
b[2] /
b[2] /
x0 / 2.;
495 f2 += (dzE * dzE * dzE - dzB * dzB * dzB) /
b[2] /
b[2] /
b[2] /
x0 / 3.;
496 meanRho += localPathLength *
rho;
497 totalELoss += betheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
498 sigmaELoss += energyLossFluctuation(pTotal, localPathLength, rho, atomicZoverA);
502 remainingPathLength -= localPathLength;
503 }
while (remainingPathLength > TGeoShape::Tolerance());
505 meanRho /= pathLength;
506 sigmaELoss2 = sigmaELoss * sigmaELoss;
517 double bendingSlope =
param.getBendingSlope();
518 double nonBendingSlope =
param.getNonBendingSlope();
519 double inverseTotalMomentum2 =
param.getInverseBendingMomentum() *
param.getInverseBendingMomentum() *
520 (1.0 + bendingSlope * bendingSlope) /
521 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
523 double pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
527 double theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength /
x0));
528 return theta02 * theta02 * inverseTotalMomentum2 * pathLength /
x0;
543 double inverseTotalMomentum2 = inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
544 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
546 double signedPathLength = dZ * TMath::Sqrt(1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
547 double pathLengthOverX0 = (
x0 > 0.) ? TMath::Abs(signedPathLength) /
x0 : TMath::Abs(signedPathLength);
551 double theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
552 theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
554 double varCoor = (
x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
555 double varSlop = theta02;
556 double covCorrSlope = (
x0 > 0.) ? signedPathLength * theta02 / 2. : 0.;
561 newParamCov(0, 0) += varCoor;
562 newParamCov(0, 1) += covCorrSlope;
563 newParamCov(1, 0) += covCorrSlope;
564 newParamCov(1, 1) += varSlop;
566 newParamCov(2, 2) += varCoor;
567 newParamCov(2, 3) += covCorrSlope;
568 newParamCov(3, 2) += covCorrSlope;
569 newParamCov(3, 3) += varSlop;
574 double dqPxydSlopeX =
575 inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
576 double dqPxydSlopeY = -inverseBendingMomentum * nonBendingSlope * nonBendingSlope * bendingSlope /
577 (1. + bendingSlope * bendingSlope) /
578 (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
580 newParamCov(4, 0) += dqPxydSlopeX * covCorrSlope;
581 newParamCov(0, 4) += dqPxydSlopeX * covCorrSlope;
582 newParamCov(4, 1) += dqPxydSlopeX * varSlop;
583 newParamCov(1, 4) += dqPxydSlopeX * varSlop;
584 newParamCov(4, 2) += dqPxydSlopeY * covCorrSlope;
585 newParamCov(2, 4) += dqPxydSlopeY * covCorrSlope;
586 newParamCov(4, 3) += dqPxydSlopeY * varSlop;
587 newParamCov(3, 4) += dqPxydSlopeY * varSlop;
588 newParamCov(4, 4) += (dqPxydSlopeX * dqPxydSlopeX + dqPxydSlopeY * dqPxydSlopeY) * varSlop;
596void TrackExtrap::addMCSEffectInAbsorber(
TrackParam&
param,
double signedPathLength,
double f0,
double f1,
double f2)
602 double bendingSlope =
param.getBendingSlope();
603 double nonBendingSlope =
param.getNonBendingSlope();
604 double inverseBendingMomentum =
param.getInverseBendingMomentum();
605 double alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
606 (1.0 + bendingSlope * bendingSlope + nonBendingSlope * nonBendingSlope);
607 double pathLength = TMath::Abs(signedPathLength);
608 double varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
609 double covCorrSlope = TMath::Sign(1., signedPathLength) * alpha2 * (pathLength * f0 - f1);
610 double varSlop = alpha2 * f0;
613 TMatrixD newParamCov(
param.getCovariances());
615 newParamCov(0, 0) += varCoor;
616 newParamCov(0, 1) += covCorrSlope;
617 newParamCov(1, 0) += covCorrSlope;
618 newParamCov(1, 1) += varSlop;
620 newParamCov(2, 2) += varCoor;
621 newParamCov(2, 3) += covCorrSlope;
622 newParamCov(3, 2) += covCorrSlope;
623 newParamCov(3, 3) += varSlop;
628 double dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
629 double dqPxydSlopeY = -inverseBendingMomentum * nonBendingSlope * nonBendingSlope * bendingSlope /
630 (1. + bendingSlope * bendingSlope) / (1. + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope);
632 newParamCov(4, 0) += dqPxydSlopeX * covCorrSlope;
633 newParamCov(0, 4) += dqPxydSlopeX * covCorrSlope;
634 newParamCov(4, 1) += dqPxydSlopeX * varSlop;
635 newParamCov(1, 4) += dqPxydSlopeX * varSlop;
636 newParamCov(4, 2) += dqPxydSlopeY * covCorrSlope;
637 newParamCov(2, 4) += dqPxydSlopeY * covCorrSlope;
638 newParamCov(4, 3) += dqPxydSlopeY * varSlop;
639 newParamCov(3, 4) += dqPxydSlopeY * varSlop;
640 newParamCov(4, 4) += (dqPxydSlopeX * dqPxydSlopeX + dqPxydSlopeY * dqPxydSlopeY) * varSlop;
644 param.setCovariances(newParamCov);
648double TrackExtrap::betheBloch(
double pTotal,
double pathLength,
double rho,
double atomicZ,
double atomicZoverA)
654 static constexpr double mK = 0.307075e-3;
655 static constexpr double me = 0.511e-3;
656 static constexpr double x0 = 0.2 * 2.303;
657 static constexpr double x1 = 3. * 2.303;
659 double bg = pTotal / SMuMass;
660 double bg2 = bg * bg;
661 double maxT = 2 * me * bg2;
664 double mI = (atomicZ < 13) ? (12. * atomicZ + 7.) * 1.e-9 : (9.76 * atomicZ + 58.8 * TMath::Power(atomicZ, -0.19)) * 1.e-9;
667 double x = TMath::Log(bg);
668 double lhwI = TMath::Log(28.816 * 1e-9 * TMath::Sqrt(rho * atomicZoverA) / mI);
677 return pathLength *
rho * (
mK * atomicZoverA * (1 +
bg2) / bg2 * (0.5 * TMath::Log(2 * me * bg2 * maxT / (mI * mI)) - bg2 / (1 + bg2) -
d2));
681double TrackExtrap::energyLossFluctuation(
double pTotal,
double pathLength,
double rho,
double atomicZoverA)
686 static constexpr double mK = 0.307075e-3;
688 double p2 = pTotal * pTotal;
689 double beta2 =
p2 / (
p2 + SMuMass * SMuMass);
691 double fwhm = 2. *
mK *
rho * pathLength * atomicZoverA /
beta2;
693 return fwhm / TMath::Sqrt(8. * log(2.));
697bool TrackExtrap::correctMCSEffectInAbsorber(TrackParam&
param,
double xVtx,
double yVtx,
double zVtx,
double errXVtx,
double errYVtx,
698 double absZBeg,
double pathLength,
double f0,
double f1,
double f2)
706 double zB = (
f1 > 0.) ? absZBeg - f2 / f1 : 0.;
709 addMCSEffectInAbsorber(
param, -pathLength, f0, f1, f2);
718 TMatrixD newParam(5, 1);
719 newParam(0, 0) = xVtx;
720 newParam(1, 0) = (
param.getNonBendingCoor() - xVtx) / (zB - zVtx);
721 newParam(2, 0) = yVtx;
722 newParam(3, 0) = (
param.getBendingCoor() - yVtx) / (zB - zVtx);
723 newParam(4, 0) =
param.getCharge() /
param.p() *
724 TMath::Sqrt(1.0 + newParam(1, 0) * newParam(1, 0) + newParam(3, 0) * newParam(3, 0)) /
725 TMath::Sqrt(1.0 + newParam(3, 0) * newParam(3, 0));
728 TMatrixD paramCovP(
param.getCovariances());
729 cov2CovP(
param.getParameters(), paramCovP);
732 TMatrixD paramCovVtx(5, 5);
734 paramCovVtx(0, 0) = errXVtx * errXVtx;
735 paramCovVtx(1, 1) = paramCovP(0, 0);
736 paramCovVtx(2, 2) = errYVtx * errYVtx;
737 paramCovVtx(3, 3) = paramCovP(2, 2);
738 paramCovVtx(4, 4) = paramCovP(4, 4);
739 paramCovVtx(1, 3) = paramCovP(0, 2);
740 paramCovVtx(3, 1) = paramCovP(2, 0);
741 paramCovVtx(1, 4) = paramCovP(0, 4);
742 paramCovVtx(4, 1) = paramCovP(4, 0);
743 paramCovVtx(3, 4) = paramCovP(2, 4);
744 paramCovVtx(4, 3) = paramCovP(4, 2);
747 TMatrixD jacob(5, 5);
749 jacob(1, 0) = -1. / (zB - zVtx);
750 jacob(1, 1) = 1. / (zB - zVtx);
751 jacob(3, 2) = -1. / (zB - zVtx);
752 jacob(3, 3) = 1. / (zB - zVtx);
755 TMatrixD tmp(paramCovVtx, TMatrixD::kMultTranspose, jacob);
756 TMatrixD newParamCov(jacob, TMatrixD::kMult, tmp);
759 covP2Cov(newParam, newParamCov);
762 param.setParameters(newParam);
764 param.setCovariances(newParamCov);
770void TrackExtrap::correctELossEffectInAbsorber(TrackParam&
param,
double eLoss,
double sigmaELoss2)
775 TMatrixD newParamCov(
param.getCovariances());
776 cov2CovP(
param.getParameters(), newParamCov);
779 double p =
param.p();
780 double e = TMath::Sqrt(p * p + SMuMass * SMuMass);
781 double eCorr = e + eLoss;
782 double pCorr = TMath::Sqrt(eCorr * eCorr - SMuMass * SMuMass);
783 double nonBendingSlope =
param.getNonBendingSlope();
784 double bendingSlope =
param.getBendingSlope();
785 param.setInverseBendingMomentum(
param.getCharge() / pCorr *
786 TMath::Sqrt(1.0 + nonBendingSlope * nonBendingSlope + bendingSlope * bendingSlope) /
787 TMath::Sqrt(1.0 + bendingSlope * bendingSlope));
790 newParamCov(4, 4) += eCorr * eCorr / pCorr / pCorr * sigmaELoss2;
793 covP2Cov(
param.getParameters(), newParamCov);
796 param.setCovariances(newParamCov);
800void TrackExtrap::cov2CovP(
const TMatrixD&
param, TMatrixD& cov)
810 TMatrixD jacob(5, 5);
815 jacob(4, 4) = -qPTot /
param(4, 0);
818 TMatrixD tmp(cov, TMatrixD::kMultTranspose, jacob);
819 cov.Mult(jacob, tmp);
823void TrackExtrap::covP2Cov(
const TMatrixD&
param, TMatrixD& covP)
833 TMatrixD jacob(5, 5);
838 jacob(4, 4) = -
param(4, 0) / qPTot;
841 TMatrixD tmp(covP, TMatrixD::kMultTranspose, jacob);
842 covP.Mult(jacob, tmp);
846void TrackExtrap::convertTrackParamForExtrap(TrackParam& trackParam,
double forwardBackward,
double*
v3)
851 v3[0] = trackParam.getNonBendingCoor();
852 v3[1] = trackParam.getBendingCoor();
853 v3[2] = trackParam.getZ();
854 double pYZ = TMath::Abs(1.0 / trackParam.getInverseBendingMomentum());
855 double pZ = pYZ / TMath::Sqrt(1.0 + trackParam.getBendingSlope() * trackParam.getBendingSlope());
856 v3[6] = TMath::Sqrt(pYZ * pYZ + pZ * pZ * trackParam.getNonBendingSlope() * trackParam.getNonBendingSlope());
857 v3[5] = -forwardBackward * pZ /
v3[6];
858 v3[3] = trackParam.getNonBendingSlope() *
v3[5];
859 v3[4] = trackParam.getBendingSlope() *
v3[5];
863void TrackExtrap::recoverTrackParam(
double*
v3,
double charge, TrackParam& trackParam)
868 trackParam.setNonBendingCoor(
v3[0]);
869 trackParam.setBendingCoor(
v3[1]);
870 trackParam.setZ(
v3[2]);
871 double pYZ =
v3[6] * TMath::Sqrt((1. -
v3[3]) * (1. +
v3[3]));
872 trackParam.setInverseBendingMomentum(
charge / pYZ);
873 trackParam.setBendingSlope(
v3[4] /
v3[5]);
874 trackParam.setNonBendingSlope(
v3[3] /
v3[5]);
878bool TrackExtrap::extrapToZRungekutta(TrackParam& trackParam,
double zEnd)
883 if (trackParam.getZ() == zEnd) {
886 double forwardBackward;
887 if (zEnd < trackParam.getZ()) {
888 forwardBackward = 1.0;
890 forwardBackward = -1.0;
894 double chargeExtrap = forwardBackward * TMath::Sign(
double(1.0), trackParam.getInverseBendingMomentum());
895 double v3[7] = {0.}, v3New[7] = {0.};
896 double dZ(0.),
step(0.);
900 double residue = zEnd - trackParam.getZ();
901 while (TMath::Abs(residue) > SRungeKuttaMaxResidue) {
903 dZ = zEnd - trackParam.getZ();
905 step = dZ * TMath::Sqrt(1.0 + trackParam.getBendingSlope() * trackParam.getBendingSlope() +
906 trackParam.getNonBendingSlope() * trackParam.getNonBendingSlope());
907 convertTrackParamForExtrap(trackParam, forwardBackward,
v3);
910 if (++stepNumber > SMaxStepNumber) {
911 LOG(warning) <<
"Too many trials";
914 step = TMath::Abs(step);
915 if (!extrapOneStepRungekutta(chargeExtrap, step,
v3, v3New)) {
918 residue = zEnd - v3New[2];
919 step *= dZ / (v3New[2] - trackParam.getZ());
920 }
while (residue * dZ < 0 && TMath::Abs(residue) > SRungeKuttaMaxResidue);
922 if (v3New[5] *
v3[5] < 0) {
923 LOG(warning) <<
"The track turned around";
927 recoverTrackParam(v3New, chargeExtrap * forwardBackward, trackParam);
931 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + residue * trackParam.getNonBendingSlope());
932 trackParam.setBendingCoor(trackParam.getBendingCoor() + residue * trackParam.getBendingSlope());
933 trackParam.setZ(zEnd);
939bool TrackExtrap::extrapToZRungekuttaV2(TrackParam& trackParam,
double zEnd)
945 if (trackParam.getZ() == zEnd) {
949 double residue = zEnd - trackParam.getZ();
950 double forwardBackward = (residue < 0) ? 1. : -1.;
952 convertTrackParamForExtrap(trackParam, forwardBackward,
v3);
953 double charge = TMath::Sign(
double(1.), trackParam.getInverseBendingMomentum());
956 double v3New[7] = {0.};
960 if (++stepNumber > SMaxStepNumber) {
961 LOG(warning) <<
"Too many trials";
966 double slopeX =
v3[3] /
v3[5];
967 double slopeY =
v3[4] /
v3[5];
968 double step = TMath::Abs(residue) * TMath::Sqrt(1.0 + slopeX * slopeX + slopeY * slopeY);
970 if (!extrapOneStepRungekutta(forwardBackward *
charge, step,
v3, v3New)) {
974 if (v3New[5] *
v3[5] < 0) {
975 LOG(warning) <<
"The track turned around";
979 residue = zEnd - v3New[2];
980 if (TMath::Abs(residue) < SRungeKuttaMaxResidueV2) {
984 for (
int i = 0;
i < 7; ++
i) {
989 if (forwardBackward * residue > 0) {
990 forwardBackward = -forwardBackward;
997 recoverTrackParam(v3New,
charge, trackParam);
1000 trackParam.setNonBendingCoor(trackParam.getNonBendingCoor() + residue * trackParam.getNonBendingSlope());
1001 trackParam.setBendingCoor(trackParam.getBendingCoor() + residue * trackParam.getBendingSlope());
1002 trackParam.setZ(zEnd);
1008bool TrackExtrap::extrapOneStepRungekutta(
double charge,
double step,
const double* vect,
double* vout)
1033 double h2(0.), h4(0.),
f[4] = {0.};
1034 double xyzt[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
1035 double a(0.),
b(0.),
c(0.), ph(0.), ph2(0.);
1036 double secxs[4] = {0.}, secys[4] = {0.}, seczs[4] = {0.}, hxp[3] = {0.};
1037 double g1(0.), g2(0.), g3(0.), g4(0.), g5(0.), g6(0.), ang2(0.), dxt(0.), dyt(0.), dzt(0.);
1038 double est(0.), at(0.), bt(0.), ct(0.), cba(0.);
1039 double f1(0.),
f2(0.),
f3(0.),
f4(0.),
rho(0.),
tet(0.), hnorm(0.), hp(0.), rho1(0.),
sint(0.), cost(0.);
1049 double maxit = 1992;
1052 const double kdlt = 1e-4;
1053 const double kdlt32 = kdlt / 32.;
1054 const double kthird = 1. / 3.;
1055 const double khalf = 0.5;
1056 const double kec = 2.9979251e-4;
1058 const double kpisqua = 9.86960440109;
1073 for (
int j = 0;
j < 7;
j++) {
1084 if (TMath::Abs(
h) > TMath::Abs(rest)) {
1088 TGeoGlobalMagField::Instance()->Field(vout,
f);
1105 secxs[0] = (
b *
f[2] -
c *
f[1]) * ph2;
1106 secys[0] = (
c *
f[0] -
a *
f[2]) * ph2;
1107 seczs[0] = (
a *
f[1] -
b *
f[0]) * ph2;
1108 ang2 = (secxs[0] * secxs[0] + secys[0] * secys[0] + seczs[0] * seczs[0]);
1109 if (ang2 > kpisqua) {
1113 dxt = h2 *
a + h4 * secxs[0];
1114 dyt = h2 *
b + h4 * secys[0];
1115 dzt = h2 *
c + h4 * seczs[0];
1123 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1125 if (ncut++ > maxcut) {
1137 TGeoGlobalMagField::Instance()->Field(xyzt,
f);
1144 secxs[1] = (bt *
f[2] - ct *
f[1]) * ph2;
1145 secys[1] = (ct *
f[0] - at *
f[2]) * ph2;
1146 seczs[1] = (at *
f[1] - bt *
f[0]) * ph2;
1150 secxs[2] = (bt *
f[2] - ct *
f[1]) * ph2;
1151 secys[2] = (ct *
f[0] - at *
f[2]) * ph2;
1152 seczs[2] = (at *
f[1] - bt *
f[0]) * ph2;
1153 dxt =
h * (
a + secxs[2]);
1154 dyt =
h * (
b + secys[2]);
1155 dzt =
h * (
c + seczs[2]);
1159 at =
a + 2. * secxs[2];
1160 bt =
b + 2. * secys[2];
1161 ct =
c + 2. * seczs[2];
1163 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1164 if (est > 2. * TMath::Abs(
h)) {
1165 if (ncut++ > maxcut) {
1177 TGeoGlobalMagField::Instance()->Field(xyzt,
f);
1180 z =
z + (
c + (seczs[0] + seczs[1] + seczs[2]) * kthird) *
h;
1181 y =
y + (
b + (secys[0] + secys[1] + secys[2]) * kthird) *
h;
1182 x =
x + (
a + (secxs[0] + secxs[1] + secxs[2]) * kthird) *
h;
1184 secxs[3] = (bt *
f[2] - ct *
f[1]) * ph2;
1185 secys[3] = (ct *
f[0] - at *
f[2]) * ph2;
1186 seczs[3] = (at *
f[1] - bt *
f[0]) * ph2;
1187 a =
a + (secxs[0] + secxs[3] + 2. * (secxs[1] + secxs[2])) * kthird;
1188 b =
b + (secys[0] + secys[3] + 2. * (secys[1] + secys[2])) * kthird;
1189 c =
c + (seczs[0] + seczs[3] + 2. * (seczs[1] + seczs[2])) * kthird;
1191 est = TMath::Abs(secxs[0] + secxs[3] - (secxs[1] + secxs[2])) +
1192 TMath::Abs(secys[0] + secys[3] - (secys[1] + secys[2])) +
1193 TMath::Abs(seczs[0] + seczs[3] - (seczs[1] + seczs[2]));
1195 if (est > kdlt && TMath::Abs(
h) > 1.e-4) {
1196 if (ncut++ > maxcut) {
1205 if (iter++ > maxit) {
1213 cba = 1. / TMath::Sqrt(
a *
a +
b *
b +
c *
c);
1224 if (rest < 1.e-5 * TMath::Abs(step)) {
1231 LOG(warning) <<
"Ruge-Kutta failed: switch to helix";
1236 f4 = TMath::Sqrt(f1 * f1 + f2 * f2 + f3 * f3);
1238 LOG(warning) <<
"magnetic field at (" << xyzt[0] <<
", " << xyzt[1] <<
", " << xyzt[2] <<
") = " <<
f4 <<
": giving up";
1256 sint = TMath::Sin(tet);
1257 cost = 2. * TMath::Sin(khalf * tet) * TMath::Sin(khalf * tet);
1261 g3 = (
tet -
sint) * hp * rho1;
1266 vout[kix] =
vect[kix] + g1 *
vect[kipx] + g2 * hxp[0] + g3 *
f1;
1267 vout[kiy] =
vect[kiy] + g1 *
vect[kipy] + g2 * hxp[1] + g3 *
f2;
1268 vout[kiz] =
vect[kiz] + g1 *
vect[kipz] + g2 * hxp[2] + g3 *
f3;
1270 vout[kipx] =
vect[kipx] + g4 *
vect[kipx] + g5 * hxp[0] + g6 *
f1;
1271 vout[kipy] =
vect[kipy] + g4 *
vect[kipy] + g5 * hxp[1] + g6 *
f2;
1272 vout[kipz] =
vect[kipz] + g4 *
vect[kipz] + g5 * hxp[2] + g6 *
f3;
1281 LOG(info) <<
"number of times extrapToZCov() is called = " << sNCallExtrapToZCov;
1282 LOG(info) <<
"number of times Field() is called = " << sNCallField;
Definition of the MCH track parameters for internal use.
Class for time synchronization of RawReader instances.
track parameters for internal use
Double_t getInverseBendingMomentum() const
return inverse bending momentum (GeV/c ** -1) times the charge (assumed forward motion)
void setZ(Double_t z)
set Z coordinate (cm)
Bool_t hasCovariances() const
return kTRUE if the covariance matrix exist, kFALSE if not
Double_t getNonBendingCoor() const
return non bending coordinate (cm)
Double_t getZ() const
return Z coordinate (cm)
void addParameters(const TMatrixD ¶meters)
add track parameters
const TMatrixD & getParameters() const
return track parameters
Double_t getNonBendingSlope() const
return non bending slope (cm ** -1)
void setCovariances(const TMatrixD &covariances)
const TMatrixD & getCovariances() const
void setBendingCoor(Double_t bendingCoor)
set bending coordinate (cm)
Double_t getBendingCoor() const
return bending coordinate (cm)
void setParameters(const TMatrixD ¶meters)
set track parameters
void updatePropagator(const TMatrixD &propagator)
Double_t getBendingSlope() const
return bending slope (cm ** -1)
void setNonBendingCoor(Double_t nonBendingCoor)
set non bending coordinate (cm)
float sum(float s, o2::dcs::DataPointValue v)
GLuint GLfloat GLfloat GLfloat x1
GLboolean GLboolean GLboolean b
GLfloat GLfloat GLfloat GLfloat v3
GLboolean GLboolean GLboolean GLboolean a
GLdouble GLdouble GLdouble z
value_T gpu::gpustd::array< value_T, 7 > & vect
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"