244 B2DEBUG(23,
"Skipping MC Track Matcher as there are no MC Particles registered in the DataStore.");
248 B2DEBUG(23,
"########## MCRecoTracksMatcherModule ############");
253 B2DEBUG(23,
"Number pattern recognition tracks is " << nPRRecoTracks);
254 B2DEBUG(23,
"Number Monte-Carlo tracks is " << nMCRecoTracks);
256 if (not nMCRecoTracks or not nPRRecoTracks) {
264 std::multimap<DetHitIdPair, WeightedRecoTrackId > mcId_by_hitId;
272 std::set<std::pair<DetHitIdPair, WeightedRecoTrackId>> prId_by_hitId;
278 std::map<DetId, int> nHits_by_detId;
292 nHits_by_detId[Const::CDC] =
m_CDCHits.getEntries();
299 Eigen::MatrixXd confusionMatrix = Eigen::MatrixXd::Zero(nPRRecoTracks, nMCRecoTracks + 1);
300 Eigen::MatrixXd weightedConfusionMatrix = Eigen::MatrixXd::Zero(nPRRecoTracks, nMCRecoTracks + 1);
304 Eigen::RowVectorXd totalNDF_by_mcId = Eigen::RowVectorXd::Zero(nMCRecoTracks + 1);
305 Eigen::RowVectorXd totalWeight_by_mcId = Eigen::RowVectorXd::Zero(nMCRecoTracks + 1);
309 Eigen::VectorXd totalNDF_by_prId = Eigen::VectorXd::Zero(nPRRecoTracks);
312 const int mcBkgId = nMCRecoTracks;
316 for (
const std::pair<const DetId, NDF>& detId_nHits_pair : nHits_by_detId) {
318 DetId detId = detId_nHits_pair.first;
319 int nHits = detId_nHits_pair.second;
322 for (HitId hitId = 0; hitId < nHits; ++hitId) {
323 DetHitIdPair detId_hitId_pair(detId, hitId);
327 const CDCHit* cdcHit = cdcHits[hitId];
335 const auto mcIds_for_detId_hitId_pair =
336 as_range(mcId_by_hitId.equal_range(detId_hitId_pair));
339 const auto prIds_for_detId_hitId_pair =
340 as_range(std::equal_range(prId_by_hitId.begin(),
343 CompDetHitIdPair()));
347 if (mcIds_for_detId_hitId_pair.empty()) {
350 RecoTrackId mcId = mcBkgId;
352 totalNDF_by_mcId(mcId) += ndfForOneHit;
353 totalWeight_by_mcId(mcId) += ndfForOneHit * mcWeight;
355 for (
const auto& detId_hitId_pair_and_mcId : mcIds_for_detId_hitId_pair) {
356 int mcId = detId_hitId_pair_and_mcId.second;
357 double mcWeight = detId_hitId_pair_and_mcId.second.weight;
358 totalNDF_by_mcId(mcId) += ndfForOneHit;
359 totalWeight_by_mcId(mcId) += ndfForOneHit * mcWeight;
365 for (
const auto& detId_hitId_pair_and_prId : prIds_for_detId_hitId_pair) {
366 RecoTrackId prId = detId_hitId_pair_and_prId.second;
367 totalNDF_by_prId(prId) += ndfForOneHit;
371 for (
const auto& detId_hitId_pair_and_prId : prIds_for_detId_hitId_pair) {
372 int prId = detId_hitId_pair_and_prId.second;
373 if (mcIds_for_detId_hitId_pair.empty()) {
376 confusionMatrix(prId, mcId) += ndfForOneHit;
377 weightedConfusionMatrix(prId, mcId) += ndfForOneHit * mcWeight;
379 for (
const auto& detId_hitId_pair_and_mcId : mcIds_for_detId_hitId_pair) {
380 int mcId = detId_hitId_pair_and_mcId.second;
381 double mcWeight = detId_hitId_pair_and_mcId.second.weight;
382 confusionMatrix(prId, mcId) += ndfForOneHit;
383 weightedConfusionMatrix(prId, mcId) += ndfForOneHit * mcWeight;
390 B2DEBUG(24,
"Confusion matrix of the event : " << std::endl << confusionMatrix);
391 B2DEBUG(24,
"Weighted confusion matrix of the event : " << std::endl << weightedConfusionMatrix);
393 B2DEBUG(24,
"totalNDF_by_mcId : " << std::endl << totalNDF_by_mcId);
394 B2DEBUG(24,
"totalWeight_by_mcId : " << std::endl << totalWeight_by_mcId);
396 B2DEBUG(24,
"totalNDF_by_prId : " << std::endl << totalNDF_by_prId);
398 Eigen::MatrixXd purityMatrix = confusionMatrix.array().colwise() / totalNDF_by_prId.array();
399 Eigen::MatrixXd efficiencyMatrix = confusionMatrix.array().rowwise() / totalNDF_by_mcId.array();
400 Eigen::MatrixXd weightedEfficiencyMatrix = weightedConfusionMatrix.array().rowwise() / totalWeight_by_mcId.array();
402 B2DEBUG(23,
"Purities");
403 B2DEBUG(23, purityMatrix);
405 B2DEBUG(23,
"Efficiencies");
406 B2DEBUG(23, efficiencyMatrix);
408 B2DEBUG(23,
"Weighted efficiencies");
409 B2DEBUG(23, weightedEfficiencyMatrix);
413 using Efficiency = float;
414 using Purity = float;
416 struct MostWeightEfficientPRId {
418 Efficiency weightedEfficiency;
420 Efficiency efficiency;
422 std::vector<MostWeightEfficientPRId> mostWeightEfficientPRId_by_mcId(nMCRecoTracks);
423 for (RecoTrackId mcId = 0; mcId < nMCRecoTracks; ++mcId) {
424 Eigen::VectorXd efficiencyCol = efficiencyMatrix.col(mcId);
425 Eigen::VectorXd weightedEfficiencyCol = weightedEfficiencyMatrix.col(mcId);
427 RecoTrackId bestPrId = 0;
428 Efficiency bestWeightedEfficiency = weightedEfficiencyCol(0);
429 Efficiency bestEfficiency = efficiencyCol(0);
430 Purity bestPurity = purityMatrix.row(0)(mcId);
434 bestWeightedEfficiency = 0;
438 for (RecoTrackId prId = 1; prId < nPRRecoTracks; ++prId) {
439 Eigen::RowVectorXd purityRow = purityMatrix.row(prId);
441 Efficiency currentWeightedEfficiency = weightedEfficiencyCol(prId);
442 Efficiency currentEfficiency = efficiencyCol(prId);
443 Purity currentPurity = purityRow(mcId);
447 currentWeightedEfficiency = 0;
450 if (std::tie(currentWeightedEfficiency, currentEfficiency, currentPurity) >
451 std::tie(bestWeightedEfficiency, bestEfficiency, bestPurity)) {
453 bestEfficiency = currentEfficiency;
454 bestWeightedEfficiency = currentWeightedEfficiency;
455 bestPurity = currentPurity;
459 bestWeightedEfficiency = weightedEfficiencyCol(bestPrId);
460 bestEfficiency = efficiencyCol(bestPrId);
461 mostWeightEfficientPRId_by_mcId[mcId] = {bestPrId, bestWeightedEfficiency, bestEfficiency};
466 struct MostPureMCId {
471 std::vector<MostPureMCId> mostPureMCId_by_prId(nPRRecoTracks);
472 for (
int prId = 0; prId < nPRRecoTracks; ++prId) {
473 Eigen::RowVectorXd purityRow = purityMatrix.row(prId);
476 Purity highestPurity = purityRow.maxCoeff(&mcId);
478 mostPureMCId_by_prId[prId] = {mcId, highestPurity};
484 RecoTrackId mcId = -1;
485 B2DEBUG(24,
"MCTrack to highest weighted efficiency PRTrack relation");
486 for (
const auto& mostWeightEfficientPRId_for_mcId : mostWeightEfficientPRId_by_mcId) {
488 const Efficiency& weightedEfficiency = mostWeightEfficientPRId_for_mcId.weightedEfficiency;
489 const RecoTrackId& prId = mostWeightEfficientPRId_for_mcId.id;
491 "mcId : " << mcId <<
" -> prId : " << prId <<
" with weighted efficiency "
492 << weightedEfficiency);
499 RecoTrackId prId = -1;
500 B2DEBUG(24,
"PRTrack to highest purity MCTrack relation");
501 for (
const auto& mostPureMCId_for_prId : mostPureMCId_by_prId) {
503 const RecoTrackId& mcId = mostPureMCId_for_prId.id;
504 const Purity& purity = mostPureMCId_for_prId.purity;
505 B2DEBUG(24,
"prId : " << prId <<
" -> mcId : " << mcId <<
" with purity " << purity);
510 int nMatched{}, nWrongCharge{}, nBackground{}, nClones{}, nClonesWrongCharge{}, nGhost{};
516 for (RecoTrackId prId = 0; prId < nPRRecoTracks; ++prId) {
519 const MostPureMCId& mostPureMCId = mostPureMCId_by_prId[prId];
521 const RecoTrackId& mcId = mostPureMCId.id;
522 const Purity& purity = mostPureMCId.purity;
528 B2DEBUG(23,
"Stored PRTrack " << prId <<
" as ghost because of too low purity");
534 if (mcId == mcBkgId) {
537 B2DEBUG(23,
"Stored PRTrack " << prId <<
" as background because of too low purity.");
550 B2ASSERT(
"No relation from MCRecoTrack to MCParticle.", mcParticle);
552 const MostWeightEfficientPRId& mostWeightEfficientPRId_for_mcId =
553 mostWeightEfficientPRId_by_mcId[mcId];
555 const RecoTrackId& mostWeightEfficientPRId = mostWeightEfficientPRId_for_mcId.id;
556 const Efficiency& weightedEfficiency = mostWeightEfficientPRId_for_mcId.weightedEfficiency;
560 const short MCParticleTrackCharge = mcParticle->
getCharge() > 0 ? 1 : -1;
564 short nPositiveCharges = 0;
565 short nNegativeCharges = 0;
569 if (fittedTracks.
size() > 0) {
571 for (
const auto& fittedTrack : fittedTracks) {
577 trackFitResult->
getChargeSign() > 0 ? nPositiveCharges++ : nNegativeCharges++;
580 trackFitResult->
getChargeSign() > 0 ? nPositiveCharges++ : nNegativeCharges++;
584 if (nPositiveCharges > 0 and nNegativeCharges > 0) {
586 "There are different charges attributed to the same track, this shouldn't happen. Continue with the majority of positive or negative charges");
590 if (fittedTracks.
size() > 0 and (nPositiveCharges > 0 or nNegativeCharges > 0)) {
591 foundTrackCharge = nPositiveCharges > nNegativeCharges ? 1 : -1;
599 if (prId == mostWeightEfficientPRId) {
600 if (foundTrackCharge != MCParticleTrackCharge) {
614 B2DEBUG(23,
"Stored PRTrack " << prId <<
" as matched.");
615 B2DEBUG(23,
"MC Match prId " << prId <<
" to mcPartId " << mcParticle->
getArrayIndex());
616 B2DEBUG(23,
"Purity rel: prId " << prId <<
" -> mcId " << mcId <<
" : " << purity);
626 B2DEBUG(23,
"Stored PRTrack " << prId <<
" as ghost because of too low efficiency.");
633 if (foundTrackCharge != MCParticleTrackCharge) {
635 ++nClonesWrongCharge;
646 B2DEBUG(23,
"Stored PRTrack " << prId <<
" as clone.");
647 B2DEBUG(23,
"MC Match prId " << prId <<
" to mcPartId " << mcParticle->
getArrayIndex());
648 B2DEBUG(23,
"Purity rel: prId " << prId <<
" -> mcId " << mcId <<
" : " << -purity);
653 B2DEBUG(23,
"Number of matches " << nMatched);
654 B2DEBUG(23,
"Number of wrongCharge " << nWrongCharge);
655 B2DEBUG(23,
"Number of clones " << nClones);
656 B2DEBUG(23,
"Number of clones wrongCharge " << nClonesWrongCharge);
657 B2DEBUG(23,
"Number of bkg " << nBackground);
658 B2DEBUG(23,
"Number of ghost " << nGhost);
662 for (RecoTrackId mcId = 0; mcId < nMCRecoTracks; ++mcId) {
666 const MostWeightEfficientPRId& mostWeightEfficiencyPRId = mostWeightEfficientPRId_by_mcId[mcId];
668 const RecoTrackId& prId = mostWeightEfficiencyPRId.id;
669 const Efficiency& weightedEfficiency = mostWeightEfficiencyPRId.weightedEfficiency;
672 B2ASSERT(
"Index of pattern recognition tracks out of range.", prId < nPRRecoTracks and prId >= 0);
676 const MostPureMCId& mostPureMCId_for_prId = mostPureMCId_by_prId[prId];
677 const RecoTrackId& mostPureMCId = mostPureMCId_for_prId.id;
680 if (mcId == mostPureMCId and
682 prRecoTrack->
getMatchingStatus() == RecoTrack::MatchingStatus::c_matchedWrongCharge or
684 prRecoTrack->
getMatchingStatus() == RecoTrack::MatchingStatus::c_cloneWrongCharge)) {
688 B2DEBUG(23,
"Efficiency rel: mcId " << mcId <<
" -> prId " << prId <<
" : " << weightedEfficiency);
696 bool isMergedMCRecoTrack =
699 prRecoTrack->
getMatchingStatus() == RecoTrack::MatchingStatus::c_matchedWrongCharge or
701 prRecoTrack->
getMatchingStatus() == RecoTrack::MatchingStatus::c_cloneWrongCharge);
703 if (isMergedMCRecoTrack) {
704 mcRecoTrack->
addRelationTo(prRecoTrack, -weightedEfficiency);
706 B2DEBUG(23,
"Efficiency rel: mcId " << mcId <<
" -> prId " << prId <<
" : " << -weightedEfficiency);
713 B2DEBUG(23,
"mcId " << mcId <<
" is missing. No relation created.");
714 B2DEBUG(23,
"is Primary " <<
m_MCRecoTracks[mcId]->getRelatedTo<MCParticle>()->isPrimaryParticle());
715 B2DEBUG(23,
"best prId " << prId <<
" with purity " << mostPureMCId_for_prId.purity <<
" -> " << mostPureMCId);
716 B2DEBUG(23,
"MC Total ndf " << totalNDF_by_mcId[mcId]);
717 B2DEBUG(23,
"MC Total weight " << totalWeight_by_mcId[mcId]);
718 B2DEBUG(23,
"MC Overlap ndf\n " << confusionMatrix.col(mcId).transpose());
719 B2DEBUG(23,
"MC Overlap weight\n " << weightedConfusionMatrix.col(mcId).transpose());
720 B2DEBUG(23,
"MC Efficiencies for the track\n" << efficiencyMatrix.col(mcId).transpose());
721 B2DEBUG(23,
"MC Weighted efficiencies for the track\n" << weightedEfficiencyMatrix.col(mcId).transpose());
724 B2DEBUG(23,
"########## End MCRecoTracksMatcherModule ############");