261 if (!mille->isOpen())
265 std::shared_ptr<genfit::GblFitter> gbl(
new genfit::GblFitter());
267 double lostWeight = -1.;
276 for (
auto& recoTrack : recoTracks) {
282 if (!track.hasFitStatus())
284 genfit::GblFitStatus* fs =
dynamic_cast<genfit::GblFitStatus*
>(track.getFitStatus());
288 if (!fs->isFittedWithReferenceTrack())
292 GblTrajectory trajectory(gbl->collectGblPoints(&track, track.getCardinalRep()), fs->hasCurvature());
294 trajectory.fit(chi2, ndf, lostWeight);
314 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
316 auto gblfs =
dynamic_cast<genfit::GblFitStatus*
>(track->getFitStatus());
318 gbl::GblTrajectory trajectory(gbl->collectGblPoints(track, track->getCardinalRep()), gblfs->hasCurvature());
320 trajectory.fit(chi2, ndf, lostWeight);
340 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
341 auto mother = list->getParticle(iParticle);
342 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
345 daughters.push_back({
346 gbl->collectGblPoints(track, track->getCardinalRep()),
350 if (daughters.size() > 1) {
351 gbl::GblTrajectory combined(daughters);
353 combined.fit(chi2, ndf, lostWeight);
365 B2RESULT(
"Vertex-constrained fit NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
376 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
377 auto mother = list->getParticle(iParticle);
378 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
380 TMatrixD extProjection(5, 3);
381 TMatrixD locProjection(3, 5);
384 TMatrixD innerTrafo(5, 3 + 3 * mother->getDaughters().size());
385 unsigned int iCol(3);
398 innerTrafo[0][iCol++] = 1.;
399 innerTrafo[1][iCol++] = 1.;
400 innerTrafo[2][iCol++] = 1.;
403 daughters.push_back({
404 gbl->collectGblPoints(track, track->getCardinalRep()),
409 if (daughters.size() > 1) {
412 TMatrixDSym vertexCov(get<TMatrixDSym>(beam));
413 TMatrixDSym vertexPrec(get<TMatrixDSym>(beam).Invert());
416 TVectorD extMeasurements(3);
417 extMeasurements[0] = vertexResidual[0];
418 extMeasurements[1] = vertexResidual[1];
419 extMeasurements[2] = vertexResidual[2];
421 TMatrixD extDeriv(3, 9);
429 TMatrixD derivatives(3, 3);
431 derivatives(0, 0) = 1.;
432 derivatives(1, 1) = 1.;
433 derivatives(2, 2) = 1.;
435 std::vector<int> labels;
437 labels.push_back(label.setParameterId(1));
438 labels.push_back(label.setParameterId(2));
439 labels.push_back(label.setParameterId(3));
446 std::vector<int> lab(globals); TMatrixD der(globals);
472 TMatrixD dLocal_dExt = extProjection;
473 TMatrixD dExt_dLocal = locProjection;
475 TVectorD locRes = dLocal_dExt * extMeasurements;
477 TMatrixD locCov = dLocal_dExt * vertexCov * dExt_dLocal;
479 TMatrixD locPrec = locCov.GetSub(3, 4, 3, 4).Invert();
480 TMatrixDSym locPrec2D(2); locPrec2D.Zero();
481 for (
int i = 0; i < 2; ++i)
482 for (
int j = 0; j < 2; ++j)
483 locPrec2D(i, j) = locPrec(i, j);
488 TVectorD locRes2D = locRes.GetSub(3, 4);
489 TMatrixD locDerivs2D = (extProjection * der).GetSub(3, 4, 0, 2);
493 daughters[0].first[0].addMeasurement(locRes2D, locPrec2D);
495 daughters[0].first[0].addGlobals(lab, locDerivs2D);
498 gbl::GblTrajectory combined(daughters);
502 combined.fit(chi2, ndf, lostWeight);
512 B2RESULT(
"Beam vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
516 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, vertexPrec);
518 combined.fit(chi2, ndf, lostWeight);
529 B2RESULT(
"Beam vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
541 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
543 auto mother = list->getParticle(iParticle);
545 if (track12.size() != 2) {
546 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
551 double motherMass = mother->getPDGMass();
552 double motherWidth = pdgdb->GetParticle(mother->getPDGCode())->Width();
557 if (motherWidth == 0.) {
559 B2WARNING(
"Using artificial width for " << pdgdb->GetParticle(mother->getPDGCode())->GetName() <<
" : " << motherWidth <<
" GeV");
563 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
565 daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
566 daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
568 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
569 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
571 TVectorD extMeasurements(1);
572 extMeasurements[0] = massResidual[0];
574 TMatrixD extDeriv(1, 9);
578 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, massPrec);
580 combined.fit(chi2, ndf, lostWeight);
592 B2RESULT(
"Mass(PDG) + vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
606 double motherMass = beam->getMass();
607 double motherWidth =
sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
611 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
613 auto mother = list->getParticle(iParticle);
615 if (track12.size() != 2) {
616 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
621 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
623 daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
624 daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
626 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
627 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
629 TVectorD extMeasurements(1);
630 extMeasurements[0] = massResidual[0];
632 TMatrixD extDeriv(1, 9);
636 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, massPrec);
638 combined.fit(chi2, ndf, lostWeight);
648 B2RESULT(
"Mass constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
662 double motherMass = beam->getMass();
663 double motherWidth =
sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
667 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
669 auto mother = list->getParticle(iParticle);
671 if (track12.size() != 2) {
672 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
677 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
679 daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
680 daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
685 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
686 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
688 TMatrixDSym extPrec(4); extPrec.Zero();
689 extPrec.SetSub(0, 0, vertexPrec);
690 extPrec(3, 3) = massPrec(0, 0);
692 TVectorD extMeasurements(4);
693 extMeasurements[0] = vertexResidual[0];
694 extMeasurements[1] = vertexResidual[1];
695 extMeasurements[2] = vertexResidual[2];
696 extMeasurements[3] = massResidual[0];
698 TMatrixD extDeriv(4, 9);
705 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
707 combined.fit(chi2, ndf, lostWeight);
719 B2RESULT(
"Mass + vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
725 B2WARNING(
"This should NOT be used for production of calibration constants for the real detector (yet)!");
734 double M = beam->getMass();
735 double E_HER = beam->getHER().E();
736 double E_LER = beam->getLER().E();
738 double pz = beam->getHER().Pz() + beam->getLER().Pz();
739 double E = (beam->getHER() + beam->getLER()).E();
741 double motherMass = beam->getMass();
742 double motherWidth =
sqrt((E_HER / M) * (E_HER / M) * beam->getCovLER()(0, 0) + (E_LER / M) * (E_LER / M) * beam->getCovHER()(0,
747 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
749 B2WARNING(
"Two body decays with full kinematic constraint not yet correct - need to resolve strange covariance provided by BeamParameters!");
751 auto mother = list->getParticle(iParticle);
754 if (track12.size() != 2) {
755 B2ERROR(
"Did not get exactly 2 fitted tracks. Skipping this mother in list " << listName);
760 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
762 daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
763 daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
765 TMatrixDSym extCov(7); extCov.Zero();
774 TMatrixD dBoost_dVect(3, 3);
775 dBoost_dVect(0, 0) = 0.; dBoost_dVect(0, 1) = 1. / pz; dBoost_dVect(0, 2) = 0.;
776 dBoost_dVect(1, 0) = 0.; dBoost_dVect(1, 1) = 0.; dBoost_dVect(1, 2) = 1. / pz;
777 dBoost_dVect(2, 0) = pz /
E; dBoost_dVect(2, 1) = 0.; dBoost_dVect(2, 2) = 0.;
779 TMatrixD dVect_dBoost(3, 3);
780 dVect_dBoost(0, 0) = 0.; dVect_dBoost(0, 1) = 0.; dVect_dBoost(0, 2) =
E / pz;
781 dVect_dBoost(1, 0) = pz; dVect_dBoost(1, 1) = 0.; dVect_dBoost(1, 2) = 0.;
782 dVect_dBoost(2, 0) = 0.; dVect_dBoost(2, 1) = pz; dVect_dBoost(2, 2) = 0.;
784 TMatrixD covBoost(3, 3);
785 for (
int i = 0; i < 3; ++i) {
786 for (
int j = i; j < 3; ++j) {
787 covBoost(j, i) = covBoost(i, j) = (beam->getCovHER() + beam->getCovLER())(i, j);
793 if (covBoost(1, 1) == 0.) covBoost(1, 1) = 1.e-4;
794 if (covBoost(2, 2) == 0.) covBoost(2, 2) = 1.e-4;
796 TMatrixD covVect = dBoost_dVect * covBoost * dVect_dBoost;
798 extCov.SetSub(3, 3, covVect);
800 extCov(6, 6) = motherWidth * motherWidth;
801 auto extPrec = extCov; extPrec.Invert();
803 TVectorD extMeasurements(7);
807 extMeasurements[3] = - (
B2Vector3D(mother->getMomentum()) - (beam->getHER().Vect() + beam->getLER().Vect()))[0];
808 extMeasurements[4] = - (
B2Vector3D(mother->getMomentum()) - (beam->getHER().Vect() + beam->getLER().Vect()))[1];
809 extMeasurements[5] = - (
B2Vector3D(mother->getMomentum()) - (beam->getHER().Vect() + beam->getLER().Vect()))[2];
810 extMeasurements[6] = - (mother->getMass() - motherMass);
812 B2INFO(
"mother mass = " << mother->getMass() <<
" and beam mass = " << beam->getMass());
814 TMatrixD extDeriv(7, 9);
828 B2WARNING(
"Primary vertex+kinematics calibration not (yet?) fully implemented!");
829 B2WARNING(
"This code is highly experimental and has (un)known issues!");
832 TMatrixD derivatives(9, 6);
833 std::vector<int> labels;
837 derivatives(0, 0) = 1.;
838 derivatives(1, 1) = 1.;
839 derivatives(2, 2) = 1.;
841 labels.push_back(label.setParameterId(1));
842 labels.push_back(label.setParameterId(2));
843 labels.push_back(label.setParameterId(3));
851 derivatives(3, 3) = mother->getMomentumMagnitude();
852 derivatives(4, 4) = mother->getMomentumMagnitude();
853 derivatives(8, 5) = (beam->getLER().
E() + beam->getHER().E()) / beam->getMass();
856 labels.push_back(label.setParameterId(4));
857 labels.push_back(label.setParameterId(5));
858 labels.push_back(label.setParameterId(6));
872 std::vector<int> lab(globals); TMatrixD der(globals);
875 TMatrixD dTwoBody_dExt(9, 7);
876 dTwoBody_dExt.Zero();
878 dTwoBody_dExt(0, 0) = 1.;
879 dTwoBody_dExt(1, 1) = 1.;
880 dTwoBody_dExt(2, 2) = 1.;
882 dTwoBody_dExt(3, 3) = 1.;
883 dTwoBody_dExt(4, 4) = 1.;
884 dTwoBody_dExt(5, 5) = 1.;
886 dTwoBody_dExt(8, 6) = 1.;
888 const TMatrixD dLocal_dExt = dfdextPlusMinus.first * dTwoBody_dExt;
889 TMatrixD dLocal_dExt_T = dLocal_dExt; dLocal_dExt_T.T();
898 TDecompSVD svd(dLocal_dExt_T);
899 TMatrixD dExt_dLocal = svd.Invert().T();
916 for (
int i = 0; i < 7; ++i) {
917 for (
int j = 0; j < 5; ++j) {
918 if (fabs(dExt_dLocal(i, j)) < 1.e-6)
919 dExt_dLocal(i, j) = 0.;
922 const TVectorD locRes = dLocal_dExt * extMeasurements;
923 const TMatrixD locPrec = dLocal_dExt * extPrec * dExt_dLocal;
925 TMatrixDSym locPrecSym(5); locPrecSym.Zero();
926 for (
int i = 0; i < 5; ++i) {
927 for (
int j = i; j < 5; ++j) {
929 locPrecSym(j, i) = locPrecSym(i, j) = (fabs(locPrec(i, j)) > 1.e-6) ? locPrec(i, j) : 0.;
933 daughters[0].first[0].addMeasurement(locRes, locPrecSym);
935 daughters[0].first[0].addGlobals(lab, dfdextPlusMinus.first * der);
943 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
947 combined.fit(chi2, ndf, lostWeight);
957 B2RESULT(
"Full kinematic-constrained fit (calibration version) results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
963 gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
967 combined.fit(chi2, ndf, lostWeight);
977 B2RESULT(
"Full kinematic-constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
1056 auto relatedRecoHitInformation =
1061 if (recoHitInformation.getFlag() == RecoHitInformation::c_pruned) {
1062 B2FATAL(
"Found pruned point in RecoTrack. Pruned tracks cannot be used in MillepedeCollector.");
1065 if (recoHitInformation.getTrackingDetector() != RecoHitInformation::c_CDC)
continue;
1071 auto kalmanFitterInfo =
dynamic_cast<genfit::KalmanFitterInfo*
>(trackPoint->getFitterInfo());
1072 if (not kalmanFitterInfo) {
1075 std::vector<double> weights = kalmanFitterInfo->getWeights();
1076 if (weights.size() == 2) {
1077 if (weights.at(0) > weights.at(1))
1078 recoHitInformation.setRightLeftInformation(RecoHitInformation::c_left);
1079 else if (weights.at(0) < weights.at(1))
1080 recoHitInformation.setRightLeftInformation(RecoHitInformation::c_right);
1082 double weightLR = weights.at(0) + weights.at(1);
1084 sumCDCWeights += weightLR - 1.;
1096 B2ERROR(
"Error in checking DAF weights from previous fit to resolve hit ambiguity. Why? Failed fit points in DAF? Skip track to be sure.");
1100 std::shared_ptr<genfit::GblFitter> gbl(
new genfit::GblFitter());
1114 genfit::MeasurementFactory<genfit::AbsMeasurement> genfitMeasurementFactory;
1118 genfit::MeasurementProducer <RecoHitInformation::UsedPXDHit, AlignablePXDRecoHit>* PXDProducer =
new genfit::MeasurementProducer
1120 genfitMeasurementFactory.addProducer(Const::PXD, PXDProducer);
1124 genfit::MeasurementProducer <RecoHitInformation::UsedSVDHit, AlignableSVDRecoHit>* SVDProducer =
new genfit::MeasurementProducer
1126 genfitMeasurementFactory.addProducer(Const::SVD, SVDProducer);
1130 genfit::MeasurementProducer <RecoHitInformation::UsedCDCHit, AlignableCDCRecoHit>* CDCProducer =
new genfit::MeasurementProducer
1132 genfitMeasurementFactory.addProducer(Const::CDC, CDCProducer);
1136 genfit::MeasurementProducer <RecoHitInformation::UsedBKLMHit, AlignableBKLMRecoHit>* BKLMProducer =
new genfit::MeasurementProducer
1138 genfitMeasurementFactory.addProducer(Const::BKLM, BKLMProducer);
1142 genfit::MeasurementProducer <RecoHitInformation::UsedEKLMHit, AlignableEKLMRecoHit>* EKLMProducer =
new genfit::MeasurementProducer
1144 genfitMeasurementFactory.addProducer(Const::EKLM, EKLMProducer);
1149 std::vector<std::shared_ptr<PXDBaseMeasurementCreator>> pxdMeasurementCreators = { std::shared_ptr<PXDBaseMeasurementCreator>(
new PXDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1150 std::vector<std::shared_ptr<SVDBaseMeasurementCreator>> svdMeasurementCreators = { std::shared_ptr<SVDBaseMeasurementCreator>(
new SVDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1153 std::vector<std::shared_ptr<CDCBaseMeasurementCreator>> cdcMeasurementCreators = { std::shared_ptr<CDCBaseMeasurementCreator>(
new CDCCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1154 std::vector<std::shared_ptr<BKLMBaseMeasurementCreator>> bklmMeasurementCreators = { std::shared_ptr<BKLMBaseMeasurementCreator>(
new BKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1155 std::vector<std::shared_ptr<EKLMBaseMeasurementCreator>> eklmMeasurementCreators = { std::shared_ptr<EKLMBaseMeasurementCreator>(
new EKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1158 std::vector<std::shared_ptr<BaseMeasurementCreator>> additionalMeasurementCreators = {};
1159 factory.
resetMeasurementCreators(pxdMeasurementCreators, svdMeasurementCreators, cdcMeasurementCreators, bklmMeasurementCreators,
1160 eklmMeasurementCreators, additionalMeasurementCreators);
1167 currentPdgCode = particle->getPDGCode();
1170 gfTrack.setCardinalRep(gfTrack.getIdForRep(trackRep));
1173 B2Vector3D vertexPos = particle->getVertex();
1174 B2Vector3D vertexMom = particle->getMomentum();
1175 gfTrack.setStateSeed(vertexPos, vertexMom);
1177 genfit::StateOnPlane vertexSOP(gfTrack.getCardinalRep());
1178 B2Vector3D vertexRPhiDir(vertexPos[0], vertexPos[1], 0);
1183 genfit::SharedPlanePtr vertexPlane(
new genfit::DetPlane(vertexPos, vertexMom));
1185 vertexSOP.setPlane(vertexPlane);
1186 vertexSOP.setPosMom(vertexPos, vertexMom);
1187 TMatrixDSym vertexCov(5);
1188 vertexCov.UnitMatrix();
1192 genfit::MeasuredStateOnPlane mop(vertexSOP, vertexCov);
1193 genfit::FullMeasurement* vertex =
new genfit::FullMeasurement(mop, Const::IR);
1194 gfTrack.insertMeasurement(vertex, 0);
1198 for (
unsigned int i = 0; i < gfTrack.getNumPoints() - 1; ++i) {
1201 genfit::PlanarMeasurement* planarMeas1 =
dynamic_cast<genfit::PlanarMeasurement*
>(gfTrack.getPointWithMeasurement(
1202 i)->getRawMeasurement(0));
1203 genfit::PlanarMeasurement* planarMeas2 =
dynamic_cast<genfit::PlanarMeasurement*
>(gfTrack.getPointWithMeasurement(
1204 i + 1)->getRawMeasurement(0));
1206 if (planarMeas1 != NULL && planarMeas2 != NULL &&
1207 planarMeas1->getDetId() == planarMeas2->getDetId() &&
1208 planarMeas1->getPlaneId() != -1 &&
1209 planarMeas1->getPlaneId() == planarMeas2->getPlaneId()) {
1216 if (hit1->
isU() && !hit2->
isU()) {
1219 }
else if (!hit1->
isU() && hit2->
isU()) {
1227 gfTrack.insertMeasurement(hit, i);
1229 gfTrack.deletePoint(i + 1);
1230 gfTrack.deletePoint(i + 1);
1234 }
catch (std::exception& e) {
1236 B2ERROR(
"SVD Cluster combination failed. This is symptomatic of pruned tracks. MillepedeCollector cannot process pruned tracks.");
1241 gbl->processTrackWithRep(&gfTrack, gfTrack.getCardinalRep(),
true);
1242 }
catch (genfit::Exception& e) {
1246 B2ERROR(
"GBL fit failed.");
1304 std::vector<TMatrixD> result;
1306 double px = mother.
getPx();
1307 double py = mother.
getPy();
1308 double pz = mother.
getPz();
1309 double pt =
sqrt(px * px + py * py);
1311 double M = motherMass;
1315 || m != mother.
getDaughter(1)->
getPDGMass()) B2FATAL(
"Only two same-mass daughters (V0->f+f- decays) allowed.");
1318 TMatrixD mother2lab(3, 3);
1319 mother2lab(0, 0) = px * pz / pt / p; mother2lab(0, 1) = - py / pt; mother2lab(0, 2) = px / p;
1320 mother2lab(1, 0) = py * pz / pt / p; mother2lab(1, 1) = px / pt; mother2lab(1, 2) = py / p;
1321 mother2lab(2, 0) = - pt / p; mother2lab(2, 1) = 0; mother2lab(2, 2) = pz / p;
1322 ROOT::Math::Rotation3D lab2mother;
1323 lab2mother.SetRotationMatrix(mother2lab); lab2mother.Invert();
1331 auto mom1 = lab2mother * boostedFrame.
getMomentum(fourVector1).Vect();
1332 auto mom2 = lab2mother * boostedFrame.
getMomentum(fourVector2).Vect();
1335 auto avgMom = 0.5 * (mom1 - mom2);
1336 if (avgMom.Z() < 0.) {
1342 double theta = atan2(avgMom.rho(), avgMom.Z());
1343 double phi = atan2(avgMom.Y(), avgMom.X());
1344 if (phi < 0.) phi += 2. * TMath::Pi();
1346 double alpha = M / 2. / m;
1347 double c1 = m *
sqrt(alpha * alpha - 1.);
1348 double c2 = 0.5 *
sqrt((alpha * alpha - 1.) / alpha / alpha * (p * p + M * M));
1350 double p3 = p * p * p;
1351 double pt3 = pt * pt * pt;
1357 TMatrixD
R = mother2lab;
1358 B2Vector3D P(sign * c1 * sin(theta) * cos(phi),
1359 sign * c1 * sin(theta) * sin(phi),
1360 p / 2. + sign * c2 * cos(theta));
1362 TMatrixD dRdpx(3, 3);
1363 dRdpx(0, 0) = - pz * (pow(px, 4.) - pow(py, 4.) - py * py * pz * pz) / pt3 / p3;
1364 dRdpx(0, 1) = px * py / pt3;
1365 dRdpx(0, 2) = (py * py + pz * pz) / p3;
1367 dRdpx(1, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1368 dRdpx(1, 1) = - py * py / pt3;
1369 dRdpx(1, 2) = px * py / p3;
1371 dRdpx(2, 0) = - px * pz * pz / pt / p3;
1373 dRdpx(2, 2) = - px * pz / p3;
1375 TMatrixD dRdpy(3, 3);
1376 dRdpy(0, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1377 dRdpy(0, 1) = - px * px / pt3;
1378 dRdpy(0, 2) = px * pz / p3;
1380 dRdpy(1, 0) = - pz * (- pow(px, 4.) - px * px * pz * pz + pow(py, 4.)) / pt3 / p3;
1381 dRdpy(1, 1) = px * py / pt3;
1382 dRdpy(1, 2) = (px * px + pz * pz) / p3;
1384 dRdpy(2, 0) = - py * pz * pz / pt / p3;
1386 dRdpy(2, 2) = - py * pz / p3;
1388 TMatrixD dRdpz(3, 3);
1389 dRdpz(0, 0) = px * pt / p3;
1391 dRdpz(0, 2) = - px * pz / p3;
1393 dRdpz(1, 0) = py * pt / p3;
1395 dRdpz(1, 2) = py * pz / p3;
1397 dRdpz(2, 0) = pz * pt / p3;
1399 dRdpz(2, 2) = (px * px + py * py) / p3;
1401 auto K = 1. / 2. / p + sign * cos(theta) * m * m * (M * M / 4. / m / m - 1.) / M / M /
sqrt(m * m * (M * M / 4. / m / m - 1.) *
1402 (M * M + p * p) / M / M);
1409 sign * c1 * cos(theta) * sin(phi),
1410 sign * c2 * (- sin(theta)));
1414 sign * c1 * sin(theta) * cos(phi),
1417 double dc1dM = m * M / (2. *
sqrt(M * M - 4. * m * m));
1418 double dc2dM = M * (4. * m * m * p * p + pow(M, 4)) / (2 * M * M * M *
sqrt((M * M - 4. * m * m) * (p * p + M * M)));
1421 sign * sin(theta) * sin(phi) * dc1dM,
1422 sign * cos(theta) * dc2dM);
1424 TMatrixD dpdz(3, 6);
1425 dpdz(0, 0) = dpdpx(0); dpdz(0, 1) = dpdpy(0); dpdz(0, 2) = dpdpz(0); dpdz(0, 3) = dpdtheta(0); dpdz(0, 4) = dpdphi(0);
1426 dpdz(0, 5) = dpdM(0);
1427 dpdz(1, 0) = dpdpx(1); dpdz(1, 1) = dpdpy(1); dpdz(1, 2) = dpdpz(1); dpdz(1, 3) = dpdtheta(1); dpdz(1, 4) = dpdphi(1);
1428 dpdz(1, 5) = dpdM(1);
1429 dpdz(2, 0) = dpdpx(2); dpdz(2, 1) = dpdpy(2); dpdz(2, 2) = dpdpz(2); dpdz(2, 3) = dpdtheta(2); dpdz(2, 4) = dpdphi(2);
1430 dpdz(2, 5) = dpdM(2);
1434 TMatrixD dfdvz(5, 9);
1435 dfdvz.SetSub(0, 0, dqdv);
1436 dfdvz.SetSub(0, 3, dqdp * dpdz);
1438 result.push_back(dfdvz);
1444 return {result[0], result[1]};