11 #include <alignment/modules/MillepedeCollector/MillepedeCollectorModule.h>
13 #include <alignment/dataobjects/MilleData.h>
14 #include <alignment/GblMultipleScatteringController.h>
15 #include <alignment/GlobalDerivatives.h>
16 #include <alignment/GlobalLabel.h>
17 #include <alignment/GlobalParam.h>
18 #include <alignment/GlobalTimeLine.h>
19 #include <alignment/Manager.h>
20 #include <alignment/reconstruction/AlignableCDCRecoHit.h>
21 #include <alignment/reconstruction/AlignablePXDRecoHit.h>
22 #include <alignment/reconstruction/AlignableSVDRecoHit.h>
23 #include <alignment/reconstruction/AlignableSVDRecoHit2D.h>
24 #include <alignment/reconstruction/AlignableBKLMRecoHit.h>
25 #include <alignment/reconstruction/AlignableEKLMRecoHit.h>
26 #include <analysis/dataobjects/ParticleList.h>
27 #include <analysis/utility/ReferenceFrame.h>
28 #include <framework/core/FileCatalog.h>
29 #include <framework/database/DBObjPtr.h>
30 #include <framework/dataobjects/EventT0.h>
31 #include <framework/dataobjects/FileMetaData.h>
32 #include <framework/datastore/StoreArray.h>
33 #include <framework/datastore/StoreObjPtr.h>
34 #include <framework/dbobjects/BeamParameters.h>
35 #include <framework/particledb/EvtGenDatabasePDG.h>
36 #include <framework/pcore/ProcHandler.h>
37 #include <mdst/dbobjects/BeamSpot.h>
38 #include <mdst/dataobjects/Track.h>
39 #include <tracking/trackFitting/fitter/base/TrackFitter.h>
40 #include <tracking/trackFitting/measurementCreator/adder/MeasurementAdder.h>
42 #include <genfit/FullMeasurement.h>
43 #include <genfit/GblFitter.h>
44 #include <genfit/KalmanFitterInfo.h>
45 #include <genfit/PlanarMeasurement.h>
46 #include <genfit/Track.h>
51 #include <TDecompSVD.h>
55 using namespace alignment;
68 setPropertyFlags(c_ParallelProcessingCertified);
69 setDescription(
"Calibration data collector for Millepede Algorithm");
72 addParam(
"tracks", m_tracks,
"Names of collections of RecoTracks (already fitted with DAF) for calibration", vector<string>({
""}));
73 addParam(
"particles", m_particles,
"Names of particle list of single particles", vector<string>());
74 addParam(
"vertices", m_vertices,
75 "Name of particle list of (mother) particles with daughters for calibration using vertex constraint", vector<string>());
76 addParam(
"primaryVertices", m_primaryVertices,
77 "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile constraint",
79 addParam(
"twoBodyDecays", m_twoBodyDecays,
80 "Name of particle list of (mother) particles with daughters for calibration using vertex + mass constraint",
82 addParam(
"primaryTwoBodyDecays", m_primaryTwoBodyDecays,
83 "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile + kinematics constraint",
85 addParam(
"primaryMassTwoBodyDecays", m_primaryMassTwoBodyDecays,
86 "Name of particle list of (mother) particles with daughters for calibration using vertex + mass constraint",
88 addParam(
"primaryMassVertexTwoBodyDecays", m_primaryMassVertexTwoBodyDecays,
89 "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile + mass constraint",
92 addParam(
"stableParticleWidth", m_stableParticleWidth,
93 "Width (in GeV/c/c) to use for invariant mass constraint for 'stable' particles (like K short). Temporary until proper solution is found.",
96 addParam(
"doublePrecision", m_doublePrecision,
"Use double (=true) or single/float (=false) precision for writing binary files",
98 addParam(
"useGblTree", m_useGblTree,
"Store GBL trajectories in a tree instead of output to binary files",
100 addParam(
"absFilePaths", m_absFilePaths,
"Use absolute paths to remember binary files. Only applies if useGblTree=False",
104 addParam(
"components", m_components,
105 "Specify which DB objects are calibrated, like ['BeamSpot', 'CDCTimeWalks'] or leave empty to use all components available.",
107 addParam(
"calibrateVertex", m_calibrateVertex,
108 "For primary vertices / two body decays, beam spot vertex calibration derivatives are added",
110 addParam(
"calibrateKinematics", m_calibrateKinematics,
111 "For primary two body decays, beam spot kinematics calibration derivatives are added",
115 addParam(
"externalIterations", m_externalIterations,
"Number of external iterations of GBL fitter",
117 addParam(
"internalIterations", m_internalIterations,
"String defining internal GBL iterations for outlier down-weighting",
119 addParam(
"recalcJacobians", m_recalcJacobians,
"Up to which external iteration propagation Jacobians should be re-calculated",
122 addParam(
"minPValue", m_minPValue,
"Minimum p-value to write out a (combined) trajectory. Set <0 to write out all.",
126 addParam(
"fitTrackT0", m_fitTrackT0,
"Add local parameter for track T0 fit in GBL",
128 addParam(
"updateCDCWeights", m_updateCDCWeights,
"Update L/R weights from previous DAF fit result",
130 addParam(
"minCDCHitWeight", m_minCDCHitWeight,
"Minimum (DAF) CDC hit weight for usage by GBL",
132 addParam(
"minUsedCDCHitFraction", m_minUsedCDCHitFraction,
"Minimum used CDC hit fraction to write out a trajectory",
135 addParam(
"hierarchyType", m_hierarchyType,
"Type of (VXD only now) hierarchy: 0 = None, 1 = Flat, 2 = Half-Shells, 3 = Full",
137 addParam(
"enablePXDHierarchy", m_enablePXDHierarchy,
"Enable PXD in hierarchy (flat or full)",
139 addParam(
"enableSVDHierarchy", m_enableSVDHierarchy,
"Enable SVD in hierarchy (flat or full)",
142 addParam(
"enableWireByWireAlignment", m_enableWireByWireAlignment,
"Enable global derivatives for wire-by-wire alignment",
144 addParam(
"enableWireSagging", m_enableWireSagging,
"Enable global derivatives for wire sagging",
148 addParam(
"events", m_eventNumbers,
149 "List of (event, run, exp) with event numbers at which payloads can change for timedep calibration.",
152 addParam(
"timedepConfig", m_timedepConfig,
153 "list{ {list{param1, param2, ...}, list{(ev1, run1, exp1), ...}}, ... }.",
157 addParam(
"customMassConfig", m_customMassConfig,
158 "dict{ list_name: (mass, width), ... } with custom mass and width to use as external measurement.",
162 void MillepedeCollectorModule::prepare()
170 if (m_tracks.empty() &&
171 m_particles.empty() &&
172 m_vertices.empty() &&
173 m_primaryVertices.empty() &&
174 m_twoBodyDecays.empty() &&
175 m_primaryTwoBodyDecays.empty() &&
176 m_primaryMassTwoBodyDecays.empty() &&
177 m_primaryMassVertexTwoBodyDecays.empty())
178 B2ERROR(
"You have to specify either arrays of single tracks or particle lists of single single particles or mothers with vertex constrained daughters.");
180 if (!m_tracks.empty()) {
181 for (
auto arrayName : m_tracks)
186 if (!m_particles.empty() || !m_vertices.empty() || !m_primaryVertices.empty()) {
196 for (
auto listName : m_particles) {
201 for (
auto listName : m_vertices) {
206 for (
auto listName : m_primaryVertices) {
212 registerObject<MilleData>(
"mille",
new MilleData(m_doublePrecision, m_absFilePaths));
214 auto gblDataTree =
new TTree(
"GblDataTree",
"GblDataTree");
215 gblDataTree->Branch<std::vector<gbl::GblData>>(
"GblData", &m_currentGblData, 32000, 99);
216 registerObject<TTree>(
"GblDataTree", gblDataTree);
218 registerObject<TH1I>(
"ndf",
new TH1I(
"ndf",
"ndf", 200, 0, 200));
219 registerObject<TH1F>(
"chi2_per_ndf",
new TH1F(
"chi2_per_ndf",
"chi2 divided by ndf", 200, 0., 50.));
220 registerObject<TH1F>(
"pval",
new TH1F(
"pval",
"pval", 100, 0., 1.));
222 registerObject<TH1F>(
"cdc_hit_fraction",
new TH1F(
"cdc_hit_fraction",
"cdc_hit_fraction", 100, 0., 1.));
223 registerObject<TH1F>(
"evt0",
new TH1F(
"evt0",
"evt0", 400, -100., 100.));
226 if (m_hierarchyType == 0)
228 else if (m_hierarchyType == 1)
230 else if (m_hierarchyType == 2)
232 else if (m_hierarchyType == 3)
238 std::vector<EventMetaData> events;
239 for (
auto& ev_run_exp : m_eventNumbers) {
240 events.push_back(
EventMetaData(std::get<0>(ev_run_exp), std::get<1>(ev_run_exp), std::get<2>(ev_run_exp)));
244 if (!m_timedepConfig.empty() && m_eventNumbers.empty()) {
245 auto autoEvents = Belle2::alignment::timeline::setupTimedepGlobalLabels(m_timedepConfig);
247 }
else if (m_timedepConfig.empty() && !m_eventNumbers.empty()) {
249 }
else if (m_timedepConfig.empty() && m_eventNumbers.empty()) {
252 B2ERROR(
"Cannot set both, event list and timedep config.");
257 AlignableCDCRecoHit::s_enableTrackT0LocalDerivative = m_fitTrackT0;
258 AlignableCDCRecoHit::s_enableWireSaggingGlobalDerivative = m_enableWireSagging;
259 AlignableCDCRecoHit::s_enableWireByWireAlignmentGlobalDerivatives = m_enableWireByWireAlignment;
262 void MillepedeCollectorModule::collect()
265 alignment::GlobalCalibrationManager::getInstance().preCollect(*emd);
270 auto mille = getObjectPtr<MilleData>(
"mille");
271 if (!mille->isOpen())
272 mille->open(getUniqueMilleName());
277 double lostWeight = -1.;
281 for (
auto arrayName : m_tracks) {
286 for (
auto& recoTrack : recoTracks) {
288 if (!fitRecoTrack(recoTrack))
291 auto& track = RecoTrackGenfitAccess::getGenfitTrack(recoTrack);
292 if (!track.hasFitStatus())
298 if (!fs->isFittedWithReferenceTrack())
302 GblTrajectory trajectory(
gbl->collectGblPoints(&track, track.getCardinalRep()), fs->hasCurvature());
304 trajectory.
fit(chi2, ndf, lostWeight);
305 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
306 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
307 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
308 if (eventT0.isValid() && eventT0->hasEventT0()) {
309 evt0 = eventT0->getEventT0();
310 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
313 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(trajectory);
319 for (
auto listName : m_particles) {
324 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
325 for (
auto& track : getParticlesTracks({list->getParticle(iParticle)},
false)) {
328 gbl::GblTrajectory trajectory(
gbl->collectGblPoints(track, track->getCardinalRep()), gblfs->hasCurvature());
330 trajectory.
fit(chi2, ndf, lostWeight);
331 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
332 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
333 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
334 if (eventT0.isValid() && eventT0->hasEventT0()) {
335 evt0 = eventT0->getEventT0();
336 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
339 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(trajectory);
345 for (
auto listName : m_vertices) {
350 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
351 auto mother = list->getParticle(iParticle);
352 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
354 for (
auto& track : getParticlesTracks(mother->getDaughters()))
355 daughters.push_back({
356 gbl->collectGblPoints(track, track->getCardinalRep()),
357 getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2)
360 if (daughters.size() > 1) {
363 combined.
fit(chi2, ndf, lostWeight);
364 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
365 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
366 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
367 if (eventT0.isValid() && eventT0->hasEventT0()) {
368 evt0 = eventT0->getEventT0();
369 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
373 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
375 B2RESULT(
"Vertex-constrained fit NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
381 for (
auto listName : m_primaryVertices) {
386 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
387 auto mother = list->getParticle(iParticle);
388 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
390 TMatrixD extProjection(5, 3);
391 TMatrixD locProjection(3, 5);
394 for (
auto& track : getParticlesTracks(mother->getDaughters())) {
397 extProjection = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2);
398 locProjection = getLocalToGlobalTransform(track->getFittedState()).GetSub(0, 2, 0, 4);
401 daughters.push_back({
402 gbl->collectGblPoints(track, track->getCardinalRep()),
403 getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2)
407 if (daughters.size() > 1) {
408 auto beam = getPrimaryVertexAndCov();
410 TMatrixDSym vertexCov(get<TMatrixDSym>(beam));
411 TMatrixDSym vertexPrec(get<TMatrixDSym>(beam).Invert());
412 B2Vector3D vertexResidual = - (mother->getVertex() - get<B2Vector3D>(beam));
414 TVectorD extMeasurements(3);
415 extMeasurements[0] = vertexResidual[0];
416 extMeasurements[1] = vertexResidual[1];
417 extMeasurements[2] = vertexResidual[2];
419 TMatrixD extDeriv(3, 3);
426 if (m_calibrateVertex) {
427 TMatrixD derivatives(3, 3);
429 derivatives(0, 0) = 1.;
430 derivatives(1, 1) = 1.;
431 derivatives(2, 2) = 1.;
433 std::vector<int> labels;
434 GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
435 labels.push_back(label.setParameterId(1));
436 labels.push_back(label.setParameterId(2));
437 labels.push_back(label.setParameterId(3));
444 std::vector<int> lab(globals); TMatrixD der(globals);
470 TMatrixD dLocal_dExt = extProjection;
471 TMatrixD dExt_dLocal = locProjection;
473 TVectorD locRes = dLocal_dExt * extMeasurements;
475 TMatrixD locCov = dLocal_dExt * vertexCov * dExt_dLocal;
477 TMatrixD locPrec = locCov.GetSub(3, 4, 3, 4).Invert();
478 TMatrixDSym locPrec2D(2); locPrec2D.Zero();
479 for (
int i = 0; i < 2; ++i)
480 for (
int j = 0; j < 2; ++j)
481 locPrec2D(i, j) = locPrec(i, j);
486 TVectorD locRes2D = locRes.GetSub(3, 4);
487 TMatrixD locDerivs2D = (extProjection * der).GetSub(3, 4, 0, 2);
491 daughters[0].first[0].addMeasurement(locRes2D, locPrec2D);
493 daughters[0].first[0].addGlobals(lab, locDerivs2D);
500 combined.
fit(chi2, ndf, lostWeight);
501 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
502 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
503 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
504 if (eventT0.isValid() && eventT0->hasEventT0()) {
505 evt0 = eventT0->getEventT0();
506 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
509 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
510 B2RESULT(
"Beam vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
516 combined.
fit(chi2, ndf, lostWeight);
517 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
518 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
519 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
520 if (eventT0.isValid() && eventT0->hasEventT0()) {
521 evt0 = eventT0->getEventT0();
522 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
525 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
527 B2RESULT(
"Beam vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
534 for (
auto listName : m_twoBodyDecays) {
539 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
541 auto mother = list->getParticle(iParticle);
542 auto track12 = getParticlesTracks(mother->getDaughters());
543 if (track12.size() != 2) {
544 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
548 auto pdgdb = EvtGenDatabasePDG::Instance();
549 double motherMass = mother->getPDGMass();
550 double motherWidth = pdgdb->GetParticle(mother->getPDGCode())->Width();
552 updateMassWidthIfSet(listName, motherMass, motherWidth);
555 if (motherWidth == 0.) {
556 motherWidth = m_stableParticleWidth * Unit::GeV;
557 B2WARNING(
"Using artificial width for " << pdgdb->GetParticle(mother->getPDGCode())->GetName() <<
" : " << motherWidth <<
" GeV");
560 auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
561 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
563 daughters.push_back({
gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
564 daughters.push_back({
gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
566 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
567 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
569 TVectorD extMeasurements(1);
570 extMeasurements[0] = massResidual[0];
572 TMatrixD extDeriv(1, 9);
578 combined.
fit(chi2, ndf, lostWeight);
581 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
582 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
583 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
584 if (eventT0.isValid() && eventT0->hasEventT0()) {
585 evt0 = eventT0->getEventT0();
586 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
590 B2RESULT(
"Mass(PDG) + vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
592 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
597 for (
auto listName : m_primaryMassTwoBodyDecays) {
604 double motherMass = beam->getMass();
605 double motherWidth = sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
607 updateMassWidthIfSet(listName, motherMass, motherWidth);
609 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
611 auto mother = list->getParticle(iParticle);
612 auto track12 = getParticlesTracks(mother->getDaughters());
613 if (track12.size() != 2) {
614 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
618 auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
619 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
621 daughters.push_back({
gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
622 daughters.push_back({
gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
624 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
625 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
627 TVectorD extMeasurements(1);
628 extMeasurements[0] = massResidual[0];
630 TMatrixD extDeriv(1, 9);
636 combined.
fit(chi2, ndf, lostWeight);
637 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
638 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
639 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
640 if (eventT0.isValid() && eventT0->hasEventT0()) {
641 evt0 = eventT0->getEventT0();
642 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
646 B2RESULT(
"Mass constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
648 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
653 for (
auto listName : m_primaryMassVertexTwoBodyDecays) {
660 double motherMass = beam->getMass();
661 double motherWidth = sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
663 updateMassWidthIfSet(listName, motherMass, motherWidth);
665 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
667 auto mother = list->getParticle(iParticle);
668 auto track12 = getParticlesTracks(mother->getDaughters());
669 if (track12.size() != 2) {
670 B2ERROR(
"Did not get 2 fitted tracks. Skipping this mother.");
674 auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
675 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
677 daughters.push_back({
gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
678 daughters.push_back({
gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
680 TMatrixDSym vertexPrec(get<TMatrixDSym>(getPrimaryVertexAndCov()).Invert());
681 B2Vector3D vertexResidual = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()));
683 TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
684 TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
686 TMatrixDSym extPrec(4); extPrec.Zero();
687 extPrec.SetSub(0, 0, vertexPrec);
688 extPrec(3, 3) = massPrec(0, 0);
690 TVectorD extMeasurements(4);
691 extMeasurements[0] = vertexResidual[0];
692 extMeasurements[1] = vertexResidual[1];
693 extMeasurements[2] = vertexResidual[2];
694 extMeasurements[3] = massResidual[0];
696 TMatrixD extDeriv(4, 9);
705 combined.
fit(chi2, ndf, lostWeight);
706 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
707 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
708 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
709 if (eventT0.isValid() && eventT0->hasEventT0()) {
710 evt0 = eventT0->getEventT0();
711 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
715 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
717 B2RESULT(
"Mass + vertex constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
722 for (
auto listName : m_primaryTwoBodyDecays) {
723 B2WARNING(
"This should NOT be used for production of calibration constants for the real detector (yet)!");
732 double M = beam->getMass();
733 double E_HER = beam->getHER().E();
734 double E_LER = beam->getLER().E();
736 double pz = (beam->getHER().Vect() + beam->getLER().Vect())[2];
737 double E = (beam->getHER() + beam->getLER()).E();
739 double motherMass = beam->getMass();
740 double motherWidth = sqrt((E_HER / M) * (E_HER / M) * beam->getCovLER()(0, 0) + (E_LER / M) * (E_LER / M) * beam->getCovHER()(0,
743 updateMassWidthIfSet(listName, motherMass, motherWidth);
745 for (
unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
747 B2WARNING(
"Two body decays with full kinematic constraint not yet correct - need to resolve strange covariance provided by BeamParameters!");
749 auto mother = list->getParticle(iParticle);
751 auto track12 = getParticlesTracks(mother->getDaughters());
752 if (track12.size() != 2) {
753 B2ERROR(
"Did not get exactly 2 fitted tracks. Skipping this mother in list " << listName);
757 auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
758 std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
760 daughters.push_back({
gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
761 daughters.push_back({
gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
763 TMatrixDSym extCov(7); extCov.Zero();
766 extCov.SetSub(0, 0, get<TMatrixDSym>(getPrimaryVertexAndCov()));
772 TMatrixD dBoost_dVect(3, 3);
773 dBoost_dVect(0, 0) = 0.; dBoost_dVect(0, 1) = 1. / pz; dBoost_dVect(0, 2) = 0.;
774 dBoost_dVect(1, 0) = 0.; dBoost_dVect(1, 1) = 0.; dBoost_dVect(1, 2) = 1. / pz;
775 dBoost_dVect(2, 0) = pz / E; dBoost_dVect(2, 1) = 0.; dBoost_dVect(2, 2) = 0.;
777 TMatrixD dVect_dBoost(3, 3);
778 dVect_dBoost(0, 0) = 0.; dVect_dBoost(0, 1) = 0.; dVect_dBoost(0, 2) = E / pz;
779 dVect_dBoost(1, 0) = pz; dVect_dBoost(1, 1) = 0.; dVect_dBoost(1, 2) = 0.;
780 dVect_dBoost(2, 0) = 0.; dVect_dBoost(2, 1) = pz; dVect_dBoost(2, 2) = 0.;
782 TMatrixD covBoost(3, 3);
783 for (
int i = 0; i < 3; ++i) {
784 for (
int j = i; j < 3; ++j) {
785 covBoost(j, i) = covBoost(i, j) = (beam->getCovHER() + beam->getCovLER())(i, j);
791 if (covBoost(1, 1) == 0.) covBoost(1, 1) = 1.e-4;
792 if (covBoost(2, 2) == 0.) covBoost(2, 2) = 1.e-4;
794 TMatrixD covVect = dBoost_dVect * covBoost * dVect_dBoost;
796 extCov.SetSub(3, 3, covVect);
798 extCov(6, 6) = motherWidth * motherWidth;
799 auto extPrec = extCov; extPrec.Invert();
801 TVectorD extMeasurements(7);
802 extMeasurements[0] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[0];
803 extMeasurements[1] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[1];
804 extMeasurements[2] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[2];
805 extMeasurements[3] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[0];
806 extMeasurements[4] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[1];
807 extMeasurements[5] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[2];
808 extMeasurements[6] = - (mother->getMass() - motherMass);
810 B2INFO(
"mother mass = " << mother->getMass() <<
" and beam mass = " << beam->getMass());
812 TMatrixD extDeriv(7, 9);
825 if (m_calibrateVertex || m_calibrateKinematics) {
826 B2WARNING(
"Primary vertex+kinematics calibration not (yet?) fully implemented!");
827 B2WARNING(
"This code is highly experimental and has (un)known issues!");
830 TMatrixD derivatives(9, 6);
831 std::vector<int> labels;
834 if (m_calibrateVertex) {
835 derivatives(0, 0) = 1.;
836 derivatives(1, 1) = 1.;
837 derivatives(2, 2) = 1.;
838 GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
839 labels.push_back(label.setParameterId(1));
840 labels.push_back(label.setParameterId(2));
841 labels.push_back(label.setParameterId(3));
848 if (m_calibrateKinematics) {
849 derivatives(3, 3) = mother->getMomentumMagnitude();
850 derivatives(4, 4) = mother->getMomentumMagnitude();
851 derivatives(8, 5) = (beam->getLER().E() + beam->getHER().E()) / beam->getMass();
853 GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
854 labels.push_back(label.setParameterId(4));
855 labels.push_back(label.setParameterId(5));
856 labels.push_back(label.setParameterId(6));
870 std::vector<int> lab(globals); TMatrixD der(globals);
873 TMatrixD dTwoBody_dExt(9, 7);
874 dTwoBody_dExt.Zero();
876 dTwoBody_dExt(0, 0) = 1.;
877 dTwoBody_dExt(1, 1) = 1.;
878 dTwoBody_dExt(2, 2) = 1.;
880 dTwoBody_dExt(3, 3) = 1.;
881 dTwoBody_dExt(4, 4) = 1.;
882 dTwoBody_dExt(5, 5) = 1.;
884 dTwoBody_dExt(8, 6) = 1.;
886 const TMatrixD dLocal_dExt = dfdextPlusMinus.first * dTwoBody_dExt;
887 TMatrixD dLocal_dExt_T = dLocal_dExt; dLocal_dExt_T.T();
896 TDecompSVD svd(dLocal_dExt_T);
897 TMatrixD dExt_dLocal = svd.Invert().T();
914 for (
int i = 0; i < 7; ++i) {
915 for (
int j = 0; j < 5; ++j) {
916 if (fabs(dExt_dLocal(i, j)) < 1.e-6)
917 dExt_dLocal(i, j) = 0.;
920 const TVectorD locRes = dLocal_dExt * extMeasurements;
921 const TMatrixD locPrec = dLocal_dExt * extPrec * dExt_dLocal;
923 TMatrixDSym locPrecSym(5); locPrecSym.Zero();
924 for (
int i = 0; i < 5; ++i) {
925 for (
int j = i; j < 5; ++j) {
927 locPrecSym(j, i) = locPrecSym(i, j) = (fabs(locPrec(i, j)) > 1.e-6) ? locPrec(i, j) : 0.;
931 daughters[0].first[0].addMeasurement(locRes, locPrecSym);
933 daughters[0].first[0].addGlobals(lab, dfdextPlusMinus.first * der);
945 combined.
fit(chi2, ndf, lostWeight);
946 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
947 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
948 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
949 if (eventT0.isValid() && eventT0->hasEventT0()) {
950 evt0 = eventT0->getEventT0();
951 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
955 B2RESULT(
"Full kinematic-constrained fit (calibration version) results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
957 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
965 combined.
fit(chi2, ndf, lostWeight);
966 getObjectPtr<TH1I>(
"ndf")->Fill(ndf);
967 getObjectPtr<TH1F>(
"chi2_per_ndf")->Fill(chi2 /
double(ndf));
968 getObjectPtr<TH1F>(
"pval")->Fill(TMath::Prob(chi2, ndf));
969 if (eventT0.isValid() && eventT0->hasEventT0()) {
970 evt0 = eventT0->getEventT0();
971 getObjectPtr<TH1F>(
"evt0")->Fill(evt0);
975 B2RESULT(
"Full kinematic-constrained fit results NDF = " << ndf <<
" Chi2/NDF = " << chi2 /
double(ndf));
977 if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
983 void MillepedeCollectorModule::closeRun()
988 auto mille = getObjectPtr<MilleData>(
"mille");
993 void MillepedeCollectorModule::finish()
999 B2ERROR(
"Cannot register binaries in FileCatalog.");
1004 const std::vector<string> parents = {fileMetaData->getLfn()};
1005 for (
auto binary : getObjectPtr<MilleData>(
"mille")->getFiles()) {
1008 milleMetaData.
setLfn(
"");
1010 FileCatalog::Instance().registerFile(binary, milleMetaData);
1019 m_currentGblData = trajectory.getData();
1021 m_currentGblData.clear();
1023 if (!m_currentGblData.empty())
1024 getObjectPtr<TTree>(
"GblDataTree")->Fill();
1026 getObjectPtr<MilleData>(
"mille")->fill(trajectory);
1030 std::string MillepedeCollectorModule::getUniqueMilleName()
1033 string name = getName();
1035 name +=
"-e" + to_string(emd->getExperiment());
1036 name +=
"-r" + to_string(emd->getRun());
1037 name +=
"-ev" + to_string(emd->getEvent());
1039 if (ProcHandler::parallelProcessingUsed())
1040 name +=
"-pid" + to_string(ProcHandler::EvtProcID());
1055 auto relatedRecoHitInformation =
1060 if (recoHitInformation.getFlag() == RecoHitInformation::c_pruned) {
1061 B2FATAL(
"Found pruned point in RecoTrack. Pruned tracks cannot be used in MillepedeCollector.");
1064 if (recoHitInformation.getTrackingDetector() != RecoHitInformation::c_CDC)
continue;
1071 if (not kalmanFitterInfo) {
1074 std::vector<double> weights = kalmanFitterInfo->getWeights();
1075 if (weights.size() == 2) {
1076 if (weights.at(0) > weights.at(1))
1077 recoHitInformation.setRightLeftInformation(RecoHitInformation::c_left);
1078 else if (weights.at(0) < weights.at(1))
1079 recoHitInformation.setRightLeftInformation(RecoHitInformation::c_right);
1081 double weightLR = weights.at(0) + weights.at(1);
1082 if (weightLR < m_minCDCHitWeight) recoHitInformation.setUseInFit(
false);
1083 sumCDCWeights += weightLR - 1.;
1090 getObjectPtr<TH1F>(
"cdc_hit_fraction")->Fill(usedCDCHitFraction);
1091 if (usedCDCHitFraction < m_minUsedCDCHitFraction)
1095 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 gbl->setOptions(m_internalIterations,
true,
true, m_externalIterations, m_recalcJacobians);
1116 if (pxdHits.isOptional()) {
1119 genfitMeasurementFactory.
addProducer(Const::PXD, PXDProducer);
1122 if (svdHits.isOptional()) {
1125 genfitMeasurementFactory.
addProducer(Const::SVD, SVDProducer);
1128 if (cdcHits.isOptional()) {
1131 genfitMeasurementFactory.
addProducer(Const::CDC, CDCProducer);
1134 if (bklmHits.isOptional()) {
1137 genfitMeasurementFactory.
addProducer(Const::BKLM, BKLMProducer);
1140 if (eklmHits.isOptional()) {
1143 genfitMeasurementFactory.
addProducer(Const::EKLM, EKLMProducer);
1148 std::vector<std::shared_ptr<PXDBaseMeasurementCreator>> pxdMeasurementCreators = { std::shared_ptr<PXDBaseMeasurementCreator>(
new PXDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1149 std::vector<std::shared_ptr<SVDBaseMeasurementCreator>> svdMeasurementCreators = { std::shared_ptr<SVDBaseMeasurementCreator>(
new SVDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1152 std::vector<std::shared_ptr<CDCBaseMeasurementCreator>> cdcMeasurementCreators = { std::shared_ptr<CDCBaseMeasurementCreator>(
new CDCCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1153 std::vector<std::shared_ptr<BKLMBaseMeasurementCreator>> bklmMeasurementCreators = { std::shared_ptr<BKLMBaseMeasurementCreator>(
new BKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1154 std::vector<std::shared_ptr<EKLMBaseMeasurementCreator>> eklmMeasurementCreators = { std::shared_ptr<EKLMBaseMeasurementCreator>(
new EKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1157 std::vector<std::shared_ptr<BaseMeasurementCreator>> additionalMeasurementCreators = {};
1158 factory.
resetMeasurementCreators(pxdMeasurementCreators, svdMeasurementCreators, cdcMeasurementCreators, bklmMeasurementCreators,
1159 eklmMeasurementCreators, additionalMeasurementCreators);
1162 auto& gfTrack = RecoTrackGenfitAccess::getGenfitTrack(recoTrack);
1164 int currentPdgCode = TrackFitter::createCorrectPDGCodeForChargedStable(Const::muon, recoTrack);
1166 currentPdgCode = particle->getPDGCode();
1168 genfit::AbsTrackRep* trackRep = RecoTrackGenfitAccess::createOrReturnRKTrackRep(recoTrack, currentPdgCode);
1169 gfTrack.setCardinalRep(gfTrack.getIdForRep(trackRep));
1172 B2Vector3D vertexPos = particle->getVertex();
1173 B2Vector3D vertexMom = particle->getMomentum();
1174 gfTrack.setStateSeed(vertexPos, vertexMom);
1177 B2Vector3D vertexRPhiDir(vertexPos[0], vertexPos[1], 0);
1184 vertexSOP.setPlane(vertexPlane);
1185 vertexSOP.setPosMom(vertexPos, vertexMom);
1186 TMatrixDSym vertexCov(5);
1187 vertexCov.UnitMatrix();
1193 gfTrack.insertMeasurement(vertex, 0);
1197 for (
unsigned int i = 0; i < gfTrack.getNumPoints() - 1; ++i) {
1201 i)->getRawMeasurement(0));
1203 i + 1)->getRawMeasurement(0));
1205 if (planarMeas1 != NULL && planarMeas2 != NULL &&
1206 planarMeas1->getDetId() == planarMeas2->getDetId() &&
1207 planarMeas1->getPlaneId() != -1 &&
1208 planarMeas1->getPlaneId() == planarMeas2->getPlaneId()) {
1215 if (hit1->
isU() && !hit2->
isU()) {
1218 }
else if (!hit1->
isU() && hit2->
isU()) {
1226 gfTrack.insertMeasurement(hit, i);
1228 gfTrack.deletePoint(i + 1);
1229 gfTrack.deletePoint(i + 1);
1233 }
catch (std::exception& e) {
1235 B2ERROR(
"SVD Cluster combination failed. This is symptomatic of pruned tracks. MillepedeCollector cannot process pruned tracks.");
1240 gbl->processTrackWithRep(&gfTrack, gfTrack.getCardinalRep(),
true);
1245 B2ERROR(
"GBL fit failed.");
1252 std::vector< genfit::Track* > MillepedeCollectorModule::getParticlesTracks(std::vector<Particle*> particles,
bool addVertexPoint)
1254 std::vector< genfit::Track* > tracks;
1255 for (
auto particle : particles) {
1256 auto belle2Track = particle->getTrack();
1258 B2WARNING(
"No Belle2::Track for particle (particle->X");
1267 auto recoTrack = belle2Track->getRelatedTo<
RecoTrack>();
1270 B2WARNING(
"No related RecoTrack for Belle2::Track (particle->Track->X)");
1275 if (!fitRecoTrack(*recoTrack, (addVertexPoint) ? particle :
nullptr))
1278 auto& track = RecoTrackGenfitAccess::getGenfitTrack(*recoTrack);
1280 if (!track.hasFitStatus()) {
1281 B2WARNING(
"Track has no fit status");
1286 B2WARNING(
"Track FitStatus is not GblFitStatus.");
1289 if (!fs->isFittedWithReferenceTrack()) {
1290 B2WARNING(
"Track is not fitted with reference track.");
1294 tracks.push_back(&track);
1300 std::pair<TMatrixD, TMatrixD> MillepedeCollectorModule::getTwoBodyToLocalTransform(
Particle& mother,
1303 std::vector<TMatrixD> result;
1305 double px = mother.getMomentum()[0];
1306 double py = mother.getMomentum()[1];
1307 double pz = mother.getMomentum()[2];
1308 double pt = sqrt(px * px + py * py);
1309 double p = mother.getMomentumMagnitude();
1310 double M = motherMass;
1311 double m = mother.getDaughter(0)->getPDGMass();
1313 if (mother.getNDaughters() != 2
1314 || m != mother.getDaughter(1)->getPDGMass()) B2FATAL(
"Only two same-mass daughters (V0->f+f- decays) allowed.");
1317 TMatrixD mother2lab(3, 3);
1318 mother2lab(0, 0) = px * pz / pt / p; mother2lab(0, 1) = - py / pt; mother2lab(0, 2) = px / p;
1319 mother2lab(1, 0) = py * pz / pt / p; mother2lab(1, 1) = px / pt; mother2lab(1, 2) = py / p;
1320 mother2lab(2, 0) = - pt / p; mother2lab(2, 1) = 0; mother2lab(2, 2) = pz / p;
1321 auto lab2mother = mother2lab; lab2mother.Invert();
1326 TLorentzVector fourVector1(lab2mother * mother.getDaughter(0)->get4Vector().Vect(), mother.getDaughter(0)->get4Vector().T());
1327 TLorentzVector fourVector2(lab2mother * mother.getDaughter(1)->get4Vector().Vect(), mother.getDaughter(1)->get4Vector().T());
1329 auto mom1 = boostedFrame.
getMomentum(fourVector1).Vect();
1330 auto mom2 = boostedFrame.
getMomentum(fourVector2).Vect();
1333 auto avgMom = 0.5 * (mom1 - mom2);
1334 if (avgMom[2] < 0.) {
1340 double theta = atan2(avgMom.Perp(), avgMom[2]);
1341 double phi = atan2(avgMom[1], avgMom[0]);
1342 if (phi < 0.) phi += 2. * TMath::Pi();
1344 double alpha = M / 2. / m;
1345 double c1 = m * sqrt(alpha * alpha - 1.);
1346 double c2 = 0.5 * sqrt((alpha * alpha - 1.) / alpha / alpha * (p * p + M * M));
1348 double p3 = p * p * p;
1349 double pt3 = pt * pt * pt;
1352 for (
auto& track : getParticlesTracks(mother.getDaughters())) {
1355 TMatrixD R = mother2lab;
1356 B2Vector3D P(sign * c1 * sin(theta) * cos(phi),
1357 sign * c1 * sin(theta) * sin(phi),
1358 p / 2. + sign * c2 * cos(theta));
1360 TMatrixD dRdpx(3, 3);
1361 dRdpx(0, 0) = - pz * (pow(px, 4.) - pow(py, 4.) - py * py * pz * pz) / pt3 / p3;
1362 dRdpx(0, 1) = px * py / pt3;
1363 dRdpx(0, 2) = (py * py + pz * pz) / p3;
1365 dRdpx(1, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1366 dRdpx(1, 1) = - py * py / pt3;
1367 dRdpx(1, 2) = px * py / p3;
1369 dRdpx(2, 0) = - px * pz * pz / pt / p3;
1371 dRdpx(2, 2) = - px * pz / p3;
1373 TMatrixD dRdpy(3, 3);
1374 dRdpy(0, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1375 dRdpy(0, 1) = - px * px / pt3;
1376 dRdpy(0, 2) = px * pz / p3;
1378 dRdpy(1, 0) = - pz * (- pow(px, 4.) - px * px * pz * pz + pow(py, 4.)) / pt3 / p3;
1379 dRdpy(1, 1) = px * py / pt3;
1380 dRdpy(1, 2) = (px * px + pz * pz) / p3;
1382 dRdpy(2, 0) = - py * pz * pz / pt / p3;
1384 dRdpy(2, 2) = - py * pz / p3;
1386 TMatrixD dRdpz(3, 3);
1387 dRdpz(0, 0) = px * pt / p3;
1389 dRdpz(0, 2) = - px * pz / p3;
1391 dRdpz(1, 0) = py * pt / p3;
1393 dRdpz(1, 2) = py * pz / p3;
1395 dRdpz(2, 0) = pz * pt / p3;
1397 dRdpz(2, 2) = (px * px + py * py) / p3;
1404 sign * c1 * cos(theta) * sin(phi),
1405 p / 2. + sign * c2 * (- sin(phi)));
1409 sign * c1 * sin(theta) * cos(phi),
1412 double dc1dM = m * M / (2. * sqrt(M * M - 4. * m * m));
1413 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)));
1416 sign * sin(theta) * sin(phi) * dc1dM,
1417 sign * cos(theta) * dc2dM);
1419 TMatrixD dpdz(3, 6);
1420 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);
1421 dpdz(0, 5) = dpdM(0);
1422 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);
1423 dpdz(1, 5) = dpdM(1);
1424 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);
1425 dpdz(2, 5) = dpdM(2);
1427 TMatrixD dqdv = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2);
1428 TMatrixD dqdp = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 3, 5);
1429 TMatrixD dfdvz(5, 9);
1430 dfdvz.SetSub(0, 0, dqdv);
1431 dfdvz.SetSub(0, 3, dqdp * dpdz);
1433 result.push_back(dfdvz);
1439 return {result[0], result[1]};
1445 const B2Vector3D& U(state.getPlane()->getU());
1446 const B2Vector3D& V(state.getPlane()->getV());
1447 const B2Vector3D& O(state.getPlane()->getO());
1448 const B2Vector3D& W(state.getPlane()->getNormal());
1450 const double* state5 = state.getState().GetMatrixArray();
1454 const TVectorD& auxInfo = state.getAuxInfo();
1455 if (auxInfo.GetNrows() == 2
1456 || auxInfo.GetNrows() == 1)
1457 spu = state.getAuxInfo()(0);
1461 state7[0] = O.
X() + state5[3] * U.
X() + state5[4] * V.
X();
1462 state7[1] = O.
Y() + state5[3] * U.
Y() + state5[4] * V.
Y();
1463 state7[2] = O.
Z() + state5[3] * U.
Z() + state5[4] * V.
Z();
1465 state7[3] = spu * (W.
X() + state5[1] * U.
X() + state5[2] * V.
X());
1466 state7[4] = spu * (W.
Y() + state5[1] * U.
Y() + state5[2] * V.
Y());
1467 state7[5] = spu * (W.
Z() + state5[1] * U.
Z() + state5[2] * V.
Z());
1470 double norm = 1. / sqrt(state7[3] * state7[3] + state7[4] * state7[4] + state7[5] * state7[5]);
1471 for (
unsigned int i = 3; i < 6; ++i) state7[i] *= norm;
1473 state7[6] = state5[0];
1475 const double AtU = state7[3] * U.
X() + state7[4] * U.
Y() + state7[5] * U.
Z();
1476 const double AtV = state7[3] * V.
X() + state7[4] * V.
Y() + state7[5] * V.
Z();
1477 const double AtW = state7[3] * W.
X() + state7[4] * W.
Y() + state7[5] * W.
Z();
1481 const double qop = state7[6];
1482 const double p = state.getCharge() / qop;
1484 TMatrixD J_Mp_6x5(6, 5);
1488 J_Mp_6x5(0, 3) = U.
X();
1489 J_Mp_6x5(1, 3) = U.
Y();
1490 J_Mp_6x5(2, 3) = U.
Z();
1492 J_Mp_6x5(0, 4) = V.
X();
1493 J_Mp_6x5(1, 4) = V.
Y();
1494 J_Mp_6x5(2, 4) = V.
Z();
1497 double fact = (-1.) * qop / p;
1498 J_Mp_6x5(3, 0) = fact * state7[3];
1499 J_Mp_6x5(4, 0) = fact * state7[4];
1500 J_Mp_6x5(5, 0) = fact * state7[5];
1502 fact = 1. / (p * AtW * AtW);
1503 J_Mp_6x5(3, 1) = fact * (U.
X() * AtW - W.
X() * AtU);
1504 J_Mp_6x5(4, 1) = fact * (U.
Y() * AtW - W.
Y() * AtU);
1505 J_Mp_6x5(5, 1) = fact * (U.
Z() * AtW - W.
Z() * AtU);
1507 J_Mp_6x5(3, 2) = fact * (V.
X() * AtW - W.
X() * AtV);
1508 J_Mp_6x5(4, 2) = fact * (V.
Y() * AtW - W.
Y() * AtV);
1509 J_Mp_6x5(5, 2) = fact * (V.
Z() * AtW - W.
Z() * AtV);
1511 return J_Mp_6x5.T();
1518 const B2Vector3D& U(state.getPlane()->getU());
1519 const B2Vector3D& V(state.getPlane()->getV());
1520 const B2Vector3D& W(state.getPlane()->getNormal());
1522 const TVectorD& state5(state.getState());
1525 const TVectorD& auxInfo = state.getAuxInfo();
1526 if (auxInfo.GetNrows() == 2
1527 || auxInfo.GetNrows() == 1)
1528 spu = state.getAuxInfo()(0);
1531 pTilde[0] = spu * (W.
X() + state5(1) * U.
X() + state5(2) * V.
X());
1532 pTilde[1] = spu * (W.
Y() + state5(1) * U.
Y() + state5(2) * V.
Y());
1533 pTilde[2] = spu * (W.
Z() + state5(1) * U.
Z() + state5(2) * V.
Z());
1535 const double pTildeMag = sqrt(pTilde[0] * pTilde[0] + pTilde[1] * pTilde[1] + pTilde[2] * pTilde[2]);
1536 const double pTildeMag2 = pTildeMag * pTildeMag;
1538 const double utpTildeOverpTildeMag2 = (U.
X() * pTilde[0] + U.
Y() * pTilde[1] + U.
Z() * pTilde[2]) / pTildeMag2;
1539 const double vtpTildeOverpTildeMag2 = (V.
X() * pTilde[0] + V.
Y() * pTilde[1] + V.
Z() * pTilde[2]) / pTildeMag2;
1543 const double qop = state5(0);
1544 const double p = state.getCharge() / qop;
1546 TMatrixD J_pM_5x6(5, 6);
1550 double fact = -1. * p / (pTildeMag * qop);
1551 J_pM_5x6(0, 3) = fact * pTilde[0];
1552 J_pM_5x6(0, 4) = fact * pTilde[1];
1553 J_pM_5x6(0, 5) = fact * pTilde[2];
1555 fact = p * spu / pTildeMag;
1556 J_pM_5x6(1, 3) = fact * (U.
X() - pTilde[0] * utpTildeOverpTildeMag2);
1557 J_pM_5x6(1, 4) = fact * (U.
Y() - pTilde[1] * utpTildeOverpTildeMag2);
1558 J_pM_5x6(1, 5) = fact * (U.
Z() - pTilde[2] * utpTildeOverpTildeMag2);
1560 J_pM_5x6(2, 3) = fact * (V.
X() - pTilde[0] * vtpTildeOverpTildeMag2);
1561 J_pM_5x6(2, 4) = fact * (V.
Y() - pTilde[1] * vtpTildeOverpTildeMag2);
1562 J_pM_5x6(2, 5) = fact * (V.
Z() - pTilde[2] * vtpTildeOverpTildeMag2);
1564 J_pM_5x6(3, 0) = U.
X();
1565 J_pM_5x6(3, 1) = U.
Y();
1566 J_pM_5x6(3, 2) = U.
Z();
1568 J_pM_5x6(4, 0) = V.
X();
1569 J_pM_5x6(4, 1) = V.
Y();
1570 J_pM_5x6(4, 2) = V.
Z();
1572 return J_pM_5x6.T();
1576 tuple<B2Vector3D, TMatrixDSym> MillepedeCollectorModule::getPrimaryVertexAndCov()
const
1579 return {beam->getIPPosition(), beam->getSizeCovMatrix()};
1582 void MillepedeCollectorModule::updateMassWidthIfSet(
string listName,
double& mass,
double& width)
1584 if (m_customMassConfig.find(listName) != m_customMassConfig.end()) {
1585 auto massWidth = m_customMassConfig.at(listName);
1586 mass = std::get<0>(massWidth);
1587 width = std::get<1>(massWidth);