g4simulation/g4eval/FillTruthRecoMatchMap.cc

0001 #include "FillTruthRecoMatchMap.h"
0002
0003 #include <g4tracking/EmbRecoMatch.h>
0004 #include <g4tracking/EmbRecoMatchContainer.h>
0005
0006 #include <trackbase_historic/PHG4ParticleSvtxMap.h>
0007 #include <trackbase_historic/PHG4ParticleSvtxMap_v1.h>
0008 #include <trackbase_historic/SvtxPHG4ParticleMap_v1.h>
0009
0010 #include <fun4all/Fun4AllReturnCodes.h>
0011
0012 #include <phool/PHCompositeNode.h>
0013 #include <phool/PHDataNode.h>
0014 #include <phool/PHIODataNode.h>
0015 #include <phool/PHNode.h>
0016 #include <phool/PHNodeIterator.h>
0017 #include <phool/PHObject.h>  // for PHObject
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>  // for PHWHERE
0020
0021 #include <boost/format.hpp>
0022
0023 #include <iostream>
0024
0025 //____________________________________________________________________________..
0026 FillTruthRecoMatchMap::FillTruthRecoMatchMap(const std::string &name)
0027   : SubsysReco(name)
0028 {
0029 }
0030
0031 //____________________________________________________________________________..
0032 int FillTruthRecoMatchMap::Init(PHCompositeNode *topNode)
0033 {
0034   if (Verbosity() > 1)
0035   {
0036     topNode->print();
0037   }
0038
0039   return Fun4AllReturnCodes::EVENT_OK;
0040 }
0041
0042 //____________________________________________________________________________..
0043 int FillTruthRecoMatchMap::InitRun(PHCompositeNode *topNode)
0044 {
0045   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0046   {
0047     return Fun4AllReturnCodes::ABORTEVENT;
0048   }
0049
0050   return Fun4AllReturnCodes::EVENT_OK;
0051 }
0052
0053 int FillTruthRecoMatchMap::createNodes(PHCompositeNode *topNode)
0054 {
0055   PHNodeIterator iter(topNode);
0056
0057   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0058
0059   if (!dstNode)
0060   {
0061     std::cout << PHWHERE << " DST node is missing, quitting" << std::endl;
0062     throw std::runtime_error("Failed to find DST node in FillTruthRecoMatchMap::createNodes");
0063   }
0064
0065   m_EmbRecoMatchContainer = findNode::getClass<EmbRecoMatchContainer>(topNode, "TRKR_EMBRECOMATCHCONTAINER");
0066   if (!m_EmbRecoMatchContainer)
0067   {
0068     std::cout << PHWHERE << " Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting " << std::endl;
0069     throw std::runtime_error(" Cannot find node TRKR_EMBRECOMATCHCONTAINER on node tree; quitting");
0070   }
0071
0072   PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "SVTX"));
0073   if (!svtxNode)
0074   {
0075     svtxNode = new PHCompositeNode("SVTX");
0076     dstNode->addNode(svtxNode);
0077   }
0078
0079   m_PHG4ParticleSvtxMap = findNode::getClass<PHG4ParticleSvtxMap>(topNode, "PHG4ParticleToRecoMap");
0080   if (!m_PHG4ParticleSvtxMap)
0081   {
0082     m_PHG4ParticleSvtxMap = new PHG4ParticleSvtxMap_v1;
0083     PHIODataNode<PHObject> *truthNode =
0084         new PHIODataNode<PHObject>(m_PHG4ParticleSvtxMap, "PHG4ParticleToRecoMap", "PHObject");
0085     svtxNode->addNode(truthNode);
0086   }
0087
0088   m_SvtxPHG4ParticleMap = findNode::getClass<SvtxPHG4ParticleMap>(topNode, "RecoToPHG4ParticleMap");
0089   if (!m_SvtxPHG4ParticleMap)
0090   {
0091     m_SvtxPHG4ParticleMap = new SvtxPHG4ParticleMap_v1;
0092     PHIODataNode<PHObject> *recoNode =
0093         new PHIODataNode<PHObject>(m_SvtxPHG4ParticleMap, "RecoToPHG4ParticleMap", "PHObject");
0094     svtxNode->addNode(recoNode);
0095   }
0096
0097   return Fun4AllReturnCodes::EVENT_OK;
0098 }
0099
0100 int FillTruthRecoMatchMap::process_event(PHCompositeNode * /*topNode*/)
0101 {
0102   if (Verbosity() > 5)
0103   {
0104     std::cout << " FillTruthRecoMatchMap::process_event() " << std::endl;
0105   }
0106   // make maps to all the matches for both truth and reco and fill the
0107   /* auto matches = m_EmbRecoMatchContainer->getMatches(); */
0108
0109   /* map<unsigned short, set<int>> truth_entries; */
0110   /* map<unsigned short, set<int>> m_reco_entries; */
0111   /* int cnt {0}; */
0112   /* for (auto& match : m_EmbRecoMatchContainer->getMatches()) { */
0113
0114   /* } */
0115
0116   SvtxPHG4ParticleMap::Map map_RtoT{};  // RtoT = "Reco to Truth"
0117   PHG4ParticleSvtxMap::Map map_TtoR{};  // TtoR = "Truth to Reco"
0118
0119   for (auto &match : m_EmbRecoMatchContainer->getMatches())
0120   {
0121     // every match is a unique match from truth to embedded particles
0122     const int gtrackID = match->idTruthTrack();
0123     const unsigned int id_reco = match->idRecoTrack();
0124     const unsigned short n_match = match->nClustersMatched();
0125     const unsigned short n_truth = match->nClustersTruth();
0126     const unsigned short n_reco = match->nClustersReco();
0127
0128     // fill the map_TtoR
0129     /* PHG4ParticleSvtxMap::WeightedRecoTrackMap*  entry_TtoR; */
0130     if (map_TtoR.find(gtrackID) == map_TtoR.end())
0131     {
0132       map_TtoR[gtrackID] = PHG4ParticleSvtxMap::WeightedRecoTrackMap{};
0133     }
0134     auto &entry_TtoR = map_TtoR[gtrackID];
0135     float weight_TtoR = (float) n_match + (float) n_truth / 100.;
0136     if (entry_TtoR.find(weight_TtoR) == entry_TtoR.end())
0137     {
0138       entry_TtoR[weight_TtoR] = {id_reco};  // i.e. std::set<unsigned int> { id_reco };
0139     }
0140     else
0141     {
0142       entry_TtoR[weight_TtoR].insert(id_reco);
0143     }
0144
0145     // fill the map_RtoT
0146     /* SvtxPHG4ParticleMap::WeightedTruthTrackMap*  entry_RtoT; */
0147     if (map_RtoT.find(id_reco) == map_RtoT.end())
0148     {
0149       map_RtoT[id_reco] = SvtxPHG4ParticleMap::WeightedTruthTrackMap{};
0150     }
0151     auto &entry_RtoT = map_RtoT[id_reco];
0152     float weight_RtoT = (float) n_match + (float) n_reco / 100.;
0153     if (entry_RtoT.find(weight_RtoT) == entry_RtoT.end())
0154     {
0155       entry_RtoT[weight_RtoT] = {gtrackID};  // i.e. std::set<int> { gtrackID }
0156     }
0157     else
0158     {
0159       entry_RtoT[weight_RtoT].insert(gtrackID);
0160     }
0161
0162     if (Verbosity() > 20)
0163     {
0164       std::cout << (boost::format("EmbRecoMatch: gtrackID(%2i) id_reco(%2i) nclusters:match(%i),gtrack(%2i),reco(%2i)") % gtrackID % ((int) id_reco) % ((int) n_match) % ((int) n_truth) % ((int) n_reco)).str()
0165                 << std::endl;
0166       std::cout << (boost::format("   -> in SvtxPHG4ParticleMap {id_reco->{weight->id_true}} = {%2i->{%5.2f->%2i}}") % ((int) id_reco) % weight_RtoT % ((int) gtrackID)).str() << std::endl;
0167       std::cout << (boost::format("   -> in PHG4ParticleSvtxMap {id_true->{weight->id_reco}} = {%2i->{%5.2f->%2i}}") % gtrackID % weight_TtoR % ((int) id_reco)).str() << std::endl;
0168     }
0169   }
0170
0171   // fill the output maps
0172   for (auto &entry : map_TtoR)
0173   {
0174     m_PHG4ParticleSvtxMap->insert(entry.first, entry.second);
0175   }
0176   for (auto &entry : map_RtoT)
0177   {
0178     m_SvtxPHG4ParticleMap->insert(entry.first, entry.second);
0179   }
0180
0181   if (Verbosity() > 15)
0182   {
0183     std::cout << PHWHERE << " Print out: " << std::endl
0184               << " Contents of SvtxPHG4ParticleMap (node \"RecoToPHG4ParticleMap\")" << std::endl
0185               << " and PHG4ParticleSvtxMap (node \"PHG4ParticleToRecoMap\")" << std::endl;
0186
0187     std::cout << " --BEGIN-- Contents of SvtxPHG4ParticleMap: " << std::endl;
0188     for (auto &iter : *m_SvtxPHG4ParticleMap)
0189     {
0190       std::cout << (boost::format("    { %2i ") % ((int) iter.first)).str();  // id_reco
0191       auto n_matches = iter.second.size();
0192       long unsigned int cnt_matches = 0;
0193       for (const auto &matches : iter.second)
0194       {
0195         if (cnt_matches == 0)
0196         {
0197           std::cout << (boost::format("-> { %5.2f -> ") % matches.first).str();
0198         }
0199         else
0200         {
0201           std::cout << (boost::format("         -> { %5.2f -> ") % matches.first).str();
0202         }
0203         auto size = matches.second.size();
0204         long unsigned int i = 0;
0205         for (auto id_true : matches.second)
0206         {
0207           if (i < size - 1)
0208           {
0209             std::cout << (boost::format("%2i, ") % id_true).str();
0210           }
0211           else
0212           {
0213             std::cout << (boost::format("%2i }") % id_true).str();
0214           }
0215           ++i;
0216         }  // end matches cnt
0217         if (cnt_matches < (n_matches - 1))
0218         {
0219           std::cout << "," << std::endl;
0220         }
0221         else
0222         {
0223           std::cout << "}" << std::endl;
0224         }
0225         ++cnt_matches;
0226       }
0227     }
0228     std::cout << " --END-- Contents of SvtxPHG4ParticleMap: " << std::endl
0229               << std::endl;
0230
0231     std::cout << " --BEGIN-- Contents of PHG4ParticleToRecoMap: " << std::endl;
0232     for (auto &iter : *m_PHG4ParticleSvtxMap)
0233     {
0234       std::cout << (boost::format("    { %2i ") % iter.first).str();  // id_true
0235       auto n_matches = iter.second.size();
0236       long unsigned int cnt_matches = 0;
0237       for (const auto &matches : iter.second)
0238       {
0239         if (cnt_matches == 0)
0240         {
0241           std::cout << (boost::format("-> { %5.2f -> ") % matches.first).str();
0242         }
0243         else
0244         {
0245           std::cout << (boost::format("         -> { %5.2f -> ") % matches.first).str();
0246         }
0247         auto size = matches.second.size();
0248         long unsigned int i = 0;
0249         for (auto id_reco : matches.second)
0250         {
0251           if (i < size - 1)
0252           {
0253             std::cout << (boost::format("%2i, ") % id_reco).str();
0254           }
0255           else
0256           {
0257             std::cout << (boost::format("%2i }") % id_reco).str();
0258           }
0259           ++i;
0260         }  // end matches cnt
0261         if (cnt_matches < (n_matches - 1))
0262         {
0263           std::cout << "," << std::endl;
0264         }
0265         else
0266         {
0267           std::cout << "}" << std::endl;
0268         }
0269         ++cnt_matches;
0270       }
0271     }
0272     std::cout << " --END-- Contents of SvtxPHG4ParticleMap: " << std::endl
0273               << std::endl;
0274   }
0275
0276   return Fun4AllReturnCodes::EVENT_OK;
0277 }
0278
0279 int FillTruthRecoMatchMap::End(PHCompositeNode * /*unused*/)
0280 {
0281   return Fun4AllReturnCodes::EVENT_OK;
0282 }