sPhenix code displayed by LXR master/​coresoftware/​offline/​packages/​decayfinder/​DecayFinder.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-12-16 09:19:54

 /*
  * Find decay topologies in HepMC
  * Cameron Dean
  * 04/06/2021
  */
 
 // This is an array of PID's which decay faster than we can see in the detector
 // This means if we find them in the HepMC record we need to skip over them
 // If your decay descriptor contains one of these IDs then it will be removed from the list before starting the search
 
 #include "DecayFinder.h"
 
 #include "DecayFinderContainerBase.h"  // for DecayFinderContainerBase::Iter
 #include "DecayFinderContainer_v1.h"   // for DecayFinderContainer_v1
 
 #include <g4main/PHG4TruthInfoContainer.h>
 
 #include <phhepmc/PHHepMCGenEvent.h>
 #include <phhepmc/PHHepMCGenEventMap.h>
 
 #include <fun4all/Fun4AllReturnCodes.h>
 
 #include <phool/PHCompositeNode.h>
 #include <phool/PHIODataNode.h>    // for PHIODataNode
 #include <phool/PHNodeIterator.h>  // for PHNodeIterator
 #include <phool/getClass.h>
 
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
 #include <HepMC/GenEvent.h>
 #include <HepMC/GenVertex.h>  // for GenVertex::particle_iterator
 #pragma GCC diagnostic pop
 
 #include <HepMC/GenParticle.h>
 #include <HepMC/IteratorRange.h>
 #include <HepMC/SimpleVector.h>
 
 #include <TDatabasePDG.h>
 
 #include <algorithm>  // for find, for_each, sort, transform
 #include <cctype>     // for toupper
 #include <cstdlib>    // for abs, size_t
 #include <iostream>   // for operator<<, endl, basic_ostream
 #include <iterator>   // for end, begin
 #include <map>        // for map, map<>::mapped_type, _Rb...
 #include <memory>     // for allocator_traits<>::value_type
 
 int listOfResonantPIDs[] = {111, 113, 213, 333, 310, 311, 313, 323, 413, 423, 513, 523, 441, 443, 100443, 9000111, 9000211, 100111, 100211, 10111,
                             10211, 9010111, 9010211, 10113, 10213, 20113, 20213, 9000113, 9000213, 100113, 100213, 9010113, 9010213, 9020113, 9020213,
                             30113, 30213, 9030113, 9030213, 9040113, 9040213, 115, 215, 10115, 10215, 9000115, 9000215, 9010115, 9010215, 117, 217,
                             9000117, 9000217, 9010117, 9010217, 119, 219, 221, 331, 9000221, 9010221, 100221, 10221, 9020221, 100331, 9030221, 10331,
                             9040221, 9050221, 9060221, 9070221, 9080221, 223, 10223, 20223, 10333, 20333, 1000223, 9000223, 9010223, 30223, 100333, 225,
                             9000225, 335, 9010225, 9020225, 10225, 9030225, 10335, 9040225, 9050225, 9060225, 9070225, 9080225, 9090225, 227, 337, 229,
                             9000229, 9010229, 9000311, 9000321, 10311, 10321, 100311, 100321, 9010311, 9010321, 9020311, 9020321, 10313, 10323, 20313,
                             20323, 100313, 100323, 9000313, 9000323, 30313, 30323, 315, 325, 9000315, 9000325, 10315, 10325, 20315, 20325, 9010315,
                             9010325, 9020315, 9020325, 317, 327, 9010317, 9010327, 319, 329, 9000319, 9000329, 10411, 10421, 10413, 10423, 20413, 20423,
                             415, 425, 431, 10431, 433, 10433, 20433, 435, 10511, 10521, 10513, 10523, 20513, 20523, 515, 525, 10531, 533, 10533, 20533, 535,
                             10541, 543, 10543, 20543, 545, 10441, 100441, 10443, 20443, 30443, 9000443, 9010443, 9020443, 445, 100445, 551, 10551, 100551,
                             110551, 200551, 210551, 10553, 20553, 30553, 110553, 120553, 130553, 210553, 220553, 9000553, 9010553, 555, 10555, 20555, 100555,
                             110555, 120555, 200555, 557, 100557, 2224, 2214, 2114, 1114, 3212, 3224, 3214, 3114, 3324, 3314, 4222, 4212, 4112, 4224, 4214,
                             4114, 4232, 4132, 4322, 4312, 4324, 4314, 4332, 4334, 4412, 4422, 4414, 4424, 4432, 4434, 4444};
 
 DecayFinder::DecayFinder()
   : SubsysReco("DECAYFINDER")
   , m_save_dst(false)
 {
 }
 
 DecayFinder::DecayFinder(const std::string& name)
   : SubsysReco(name)
   , m_save_dst(false)
 {
 }
 
 int DecayFinder::Init(PHCompositeNode* topNode)
 {
   if (Verbosity() >= VERBOSITY_SOME)
   {
     std::cout << "DecayFinder name: " << Name() << std::endl;
     std::cout << "Decay descriptor: " << m_decayDescriptor << std::endl;
   }
 
   int canSearchDecay = parseDecayDescriptor();
 
   if (m_save_dst)
   {
     createDecayNode(topNode);
   }
 
   return canSearchDecay;
 }
 
 int DecayFinder::process_event(PHCompositeNode* topNode)
 {
   bool decayFound = findDecay(topNode);
 
   if (decayFound && m_save_dst && Verbosity() >= VERBOSITY_MORE)
   {
     printNode(topNode);
   }
 
   if (m_triggerOnDecay && !decayFound)
   {
     if (Verbosity() >= VERBOSITY_MORE)
     {
       std::cout << "The decay, " << m_decayDescriptor << " was not found or reconstructable in this event, skipping" << std::endl;
     }
     return Fun4AllReturnCodes::ABORTEVENT;
   }
   else
   {
     return Fun4AllReturnCodes::EVENT_OK;
   }
 }
 
 int DecayFinder::End(PHCompositeNode* /*topNode*/)
 {
   if (Verbosity() >= VERBOSITY_SOME)
   {
     printInfo();
   }
 
   return 0;
 }
 
 /*
  *
  *This function tries to intepret the string you wrote
  *Will be semi-independent of the rest of this class
  *(semi- because we need this to set the charges and strings)
  *
  */
 int DecayFinder::parseDecayDescriptor()
 {
   bool ddCanBeParsed = true;
 
   size_t daughterLocator;
 
   std::string mother;
   std::string intermediate;
   std::string daughter;
 
   int mother_charge = 0;
   std::vector<int> intermediates_charge;
   std::vector<int> daughters_charge;
 
   std::string decayArrow = "->";
   std::string chargeIndicator = "^";
   std::string startIntermediate = "{";
   std::string endIntermediate = "}";
 
   // These tracks require a + or - after their name for TDatabasePDG
   std::string specialTracks[] = {"e", "mu", "pi", "K"};
 
   std::string manipulateDecayDescriptor = m_decayDescriptor;
 
   // Remove all white space before we begin
   size_t pos;
   while ((pos = manipulateDecayDescriptor.find(' ')) != std::string::npos)
   {
     manipulateDecayDescriptor.replace(pos, 1, "");
   }
 
   // Check for charge conjugate requirement
   std::string checkForCC = manipulateDecayDescriptor.substr(0, 1) + manipulateDecayDescriptor.substr(manipulateDecayDescriptor.size() - 3, 3);
   std::for_each(checkForCC.begin(), checkForCC.end(), [](char& c)
                 { c = ::toupper(c); });
 
   // Remove the CC check if needed
   if (checkForCC == "[]CC")
   {
     manipulateDecayDescriptor = manipulateDecayDescriptor.substr(1, manipulateDecayDescriptor.size() - 4);
     m_getChargeConjugate = true;
   }
 
   // Try and find the initial particle
   size_t findMotherEndPoint = manipulateDecayDescriptor.find(decayArrow);
   mother = manipulateDecayDescriptor.substr(0, findMotherEndPoint);
   if (findParticle(mother))
   {
     m_mother_ID = abs(get_pdgcode(mother));
   }
   else
   {
     ddCanBeParsed = false;
   }
   manipulateDecayDescriptor.erase(0, findMotherEndPoint + decayArrow.length());
 
   // Try and find the intermediates
   while ((pos = manipulateDecayDescriptor.find(startIntermediate)) != std::string::npos)
   {
     size_t findIntermediateStartPoint = manipulateDecayDescriptor.find(startIntermediate, pos);
     size_t findIntermediateEndPoint = manipulateDecayDescriptor.find(endIntermediate, pos);
     std::string intermediateDecay = manipulateDecayDescriptor.substr(pos + 1, findIntermediateEndPoint - (pos + 1));
 
     intermediate = intermediateDecay.substr(0, intermediateDecay.find(decayArrow));
     if (findParticle(intermediate))
     {
       m_intermediates_ID.push_back(abs(get_pdgcode(intermediate)));
     }
     else
     {
       ddCanBeParsed = false;
     }
 
     // Now find the daughters associated to this intermediate
     int nDaughters = 0;
     intermediateDecay.erase(0, intermediateDecay.find(decayArrow) + decayArrow.length());
     while ((daughterLocator = intermediateDecay.find(chargeIndicator)) != std::string::npos)
     {
       daughter = intermediateDecay.substr(0, daughterLocator);
       std::string daughterChargeString = intermediateDecay.substr(daughterLocator + 1, 1);
       if (std::find(std::begin(specialTracks), std::end(specialTracks), daughter) != std::end(specialTracks))
       {
         daughter += daughterChargeString;
       }
       if (findParticle(daughter))
       {
         int daughter_ID = abs(get_pdgcode(daughter));
         m_daughters_ID.push_back(daughter_ID);
         daughters_charge.push_back(get_charge(daughter));
 
         if (daughter_ID == 22)
         {
           m_hasPhotonDaughter = m_allowPhotons = true; 
         }
       }
       else
       {
         ddCanBeParsed = false;
       }
       intermediateDecay.erase(0, daughterLocator + 2);
       ++nDaughters;
     }
     manipulateDecayDescriptor.erase(findIntermediateStartPoint, findIntermediateEndPoint + 1 - findIntermediateStartPoint);
     m_nTracksFromIntermediates.push_back(nDaughters);
     m_nTracksFromMother += 1;
   }
 
   // Now find any remaining reconstructable tracks from the mother
   while ((daughterLocator = manipulateDecayDescriptor.find(chargeIndicator)) != std::string::npos)
   {
     daughter = manipulateDecayDescriptor.substr(0, daughterLocator);
     std::string daughterChargeString = manipulateDecayDescriptor.substr(daughterLocator + 1, 1);
     if (std::find(std::begin(specialTracks), std::end(specialTracks), daughter) != std::end(specialTracks))
     {
       daughter += daughterChargeString;
     }
     if (findParticle(daughter))
     {
       m_daughters_ID.push_back(abs(get_pdgcode(daughter)));
       daughters_charge.push_back(get_charge(daughter));
     }
     else
     {
       ddCanBeParsed = false;
     }
     manipulateDecayDescriptor.erase(0, daughterLocator + 2);
     m_nTracksFromMother += 1;
   }
 
   unsigned int trackStart = 0;
   unsigned int trackEnd = 0;
   for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i)
   {
     trackStart = trackEnd;
     trackEnd = m_nTracksFromIntermediates[i] + trackStart;
 
     int vtxCharge = 0;
 
     for (unsigned int j = trackStart; j < trackEnd; ++j)
     {
       vtxCharge += daughters_charge[j];
     }
 
     intermediates_charge.push_back(vtxCharge);
   }
 
   for (unsigned int i = 0; i < m_daughters_ID.size(); ++i)
   {
     mother_charge += daughters_charge[i];
   }
 
   m_mother_ID = mother_charge == 0 ? m_mother_ID : mother_charge * m_mother_ID;
   for (unsigned int i = 0; i < m_intermediates_ID.size(); ++i)
   {
     m_intermediates_ID[i] = intermediates_charge[i] == 0 ? m_intermediates_ID[i] : intermediates_charge[i] * m_intermediates_ID[i];
     m_motherDecayProducts.push_back(m_intermediates_ID[i]);
   }
   for (unsigned int i = 0; i < m_daughters_ID.size(); ++i)
   {
     m_daughters_ID[i] = daughters_charge[i] == 0 ? m_daughters_ID[i] : daughters_charge[i] * m_daughters_ID[i];
     if (i >= trackEnd)
     {
       m_motherDecayProducts.push_back(m_daughters_ID[i]);
     }
   }
 
   if (ddCanBeParsed)
   {
     if (Verbosity() >= VERBOSITY_MORE)
     {
       std::cout << "Your decay descriptor can be parsed" << std::endl;
     }
     return 0;
   }
   else
   {
     if (Verbosity() >= VERBOSITY_SOME)
     {
       std::cout << "Your decay descriptor cannot be parsed, " << Name() << " will not be registered" << std::endl;
     }
     return Fun4AllReturnCodes::DONOTREGISTERSUBSYSTEM;
   }
 }
 
 /*
  * Main body for searching your decay
  * Slight issue if you limit the generator decay
  * as decays wont enter the HepMC record
  * need a switch to go to Geant4 record
  */
 bool DecayFinder::findDecay(PHCompositeNode* topNode)
 {
   bool decayWasFound = false;
   bool reconstructableDecayWasFound = false;
   bool aTrackFailedPT = false;
   bool aTrackFailedETA = false;
   bool aMotherHasPhoton = false;
   bool aMotherHasPi0 = false;
   std::vector<int> correctMotherProducts;
   std::vector<int> positive_motherDecayProducts;
 
   int n = sizeof(listOfResonantPIDs) / sizeof(listOfResonantPIDs[0]);
   for (int i : m_intermediates_ID)
   {
     n = deleteElement(listOfResonantPIDs, n, i);
   }
 
   n = deleteElement(listOfResonantPIDs, n, m_mother_ID);
 
   for (int m_motherDecayProduct : m_motherDecayProducts)
   {
     positive_motherDecayProducts.push_back(std::abs(m_motherDecayProduct));
   }
 
   m_truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
   if (!m_truthinfo)
   {
     std::cout << "DecayFinder: Missing node G4TruthInfo" << std::endl;
   }
 
   m_geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
   if (!m_geneventmap)
   {
     std::cout << "DecayFinder: Missing node PHHepMCGenEventMap" << std::endl;
   }
 
   if (!m_truthinfo && !m_geneventmap)
   {
     std::cout << "You have neither the PHHepMCGenEventMap or G4TruthInfo nodes" << std::endl;
     std::cout << "DecayFinder will crash, exiting now!" << std::endl;
     exit(1);
   }
 
   if (m_truthinfo && !m_geneventmap)  // This should use the truth info container if we have no HepMC record
   {
     if (Verbosity() >= VERBOSITY_SOME)
     {
       std::cout << "You have a valid G4TruthInfo but no PHHepMCGenEventMap container" << std::endl;
       std::cout << "You are probably using a particle gun" << std::endl;
     }
 
     PHG4TruthInfoContainer::ConstRange range = m_truthinfo->GetParticleRange();
     for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
     {
       PHG4Particle* g4particle = iter->second;
       int this_pid = m_getChargeConjugate ? abs(g4particle->get_pid()) : g4particle->get_pid();
       if (this_pid == m_mother_ID)
       {
         if (Verbosity() >= VERBOSITY_MAX)
         {
           std::cout << "parent->pdg_id(): " << g4particle->get_pid() << std::endl;
         }
 
         bool breakOut = false;
         correctMotherProducts.clear();
         decayChain.clear();
         decayChain.push_back(std::make_pair(std::make_pair(g4particle->get_primary_id(), g4particle->get_barcode()), g4particle->get_pid()));
 
         searchGeant4Record(g4particle->get_barcode(), g4particle->get_pid(), positive_motherDecayProducts, breakOut, aMotherHasPhoton, aMotherHasPi0, aTrackFailedPT, aTrackFailedETA, correctMotherProducts);
 
         if (breakOut)
         {
           continue;
         }
 
         decayWasFound = compareDecays(m_motherDecayProducts, correctMotherProducts);
 
         if (decayWasFound)
         {
           m_counter += 1;
           if (aTrackFailedPT && !aTrackFailedETA)
           {
             m_nCandFail_pT += 1;
           }
           else if (!aTrackFailedPT && aTrackFailedETA)
           {
             m_nCandFail_eta += 1;
           }
           else if (aTrackFailedPT && aTrackFailedETA)
           {
             m_nCandFail_pT_and_eta += 1;
           }
           else
           {
             m_nCandReconstructable += 1;
             reconstructableDecayWasFound = true;
             if (m_save_dst)
             {
               fillDecayNode(topNode, decayChain);
             }
             if (aMotherHasPhoton && !aMotherHasPi0)
             {
               m_nCandHas_Photon += 1;
             }
             else if (!aMotherHasPhoton && aMotherHasPi0)
             {
               m_nCandHas_Pi0 += 1;
             }
             else if (aMotherHasPhoton && aMotherHasPi0)
             {
               m_nCandHas_Photon_and_Pi0 += 1;
             }
             else
             {
               m_nCandHas_noPhoton_and_noPi0 += 1;
             }
           }
         }
       }
     }
     return reconstructableDecayWasFound;
   }
 
   for (PHHepMCGenEventMap::ConstReverseIter riter = m_geneventmap->rbegin();
        riter != m_geneventmap->rend(); ++riter)
   {
     m_genevt = riter->second;
     if (!m_genevt)
     {
       std::cout << "DecayFinder: Missing node PHHepMCGenEvent" << std::endl;
       return false;
     }
 
     HepMC::GenEvent* theEvent = m_genevt->getEvent();
 
     for (HepMC::GenEvent::particle_const_iterator p = theEvent->particles_begin(); p != theEvent->particles_end(); ++p)
     {
       int this_pid = m_getChargeConjugate ? abs((*p)->pdg_id()) : (*p)->pdg_id();
       if (this_pid == m_mother_ID)
       {
         if (Verbosity() >= VERBOSITY_MAX)
         {
           std::cout << "parent->pdg_id(): " << (*p)->pdg_id() << std::endl;
         }
 
         bool breakOut = false;
         correctMotherProducts.clear();
         decayChain.clear();
         decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*p)->barcode()), (*p)->pdg_id()));
 
         // Make sure that the mother has a decay in our record
         if (!(*p)->end_vertex())  // Mother has no end vertex, decay volume was limited
         {
           searchGeant4Record((*p)->barcode(), (*p)->pdg_id(), positive_motherDecayProducts, breakOut, aMotherHasPhoton, aMotherHasPi0, aTrackFailedPT, aTrackFailedETA, correctMotherProducts);
         }
         else
         {
           searchHepMCRecord((*p), positive_motherDecayProducts, breakOut, aMotherHasPhoton, aMotherHasPi0, aTrackFailedPT, aTrackFailedETA, correctMotherProducts);
         }
 
         if (breakOut)
         {
           continue;
         }
 
         decayWasFound = compareDecays(m_motherDecayProducts, correctMotherProducts);
 
         if (decayWasFound)
         {
           m_counter += 1;
           if (aTrackFailedPT && !aTrackFailedETA)
           {
             m_nCandFail_pT += 1;
           }
           else if (!aTrackFailedPT && aTrackFailedETA)
           {
             m_nCandFail_eta += 1;
           }
           else if (aTrackFailedPT && aTrackFailedETA)
           {
             m_nCandFail_pT_and_eta += 1;
           }
           else
           {
             m_nCandReconstructable += 1;
             reconstructableDecayWasFound = true;
             if (m_save_dst)
             {
               fillDecayNode(topNode, decayChain);
             }
             if (aMotherHasPhoton && !aMotherHasPi0)
             {
               m_nCandHas_Photon += 1;
             }
             else if (!aMotherHasPhoton && aMotherHasPi0)
             {
               m_nCandHas_Pi0 += 1;
             }
             else if (aMotherHasPhoton && aMotherHasPi0)
             {
               m_nCandHas_Photon_and_Pi0 += 1;
             }
             else
             {
               m_nCandHas_noPhoton_and_noPi0 += 1;
             }
           }
         }
       }
     }
   }
 
   return reconstructableDecayWasFound;
 }
 
 /*
  * Function to search HepMC record
  * Can switch to Geant4 search if needed
  */
 void DecayFinder::searchHepMCRecord(HepMC::GenParticle* particle, std::vector<int> decayProducts,
                                     bool& breakLoop, bool& hasPhoton, bool& hasPi0, bool& failedPT, bool& failedETA,
                                     std::vector<int>& actualDecayProducts)
 {
   for (HepMC::GenVertex::particle_iterator children = particle->end_vertex()->particles_begin(HepMC::children);
        children != particle->end_vertex()->particles_end(HepMC::children); ++children)
   {
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "--children->pdg_id(): " << (*children)->pdg_id() << std::endl;
     }
 
     if ((*children)->pdg_id() == 22)
     {
       hasPhoton = true;
     }
     if ((*children)->pdg_id() == 111)
     {
       hasPi0 = true;
     }
     if ((!m_allowPhotons && (*children)->pdg_id() == 22) || (!m_allowPi0 && (*children)->pdg_id() == 111))
     {
       breakLoop = true;
       break;
     }
 
     // This is one of the children we are looking for
     if (std::find(decayProducts.begin(), decayProducts.end(),
                   std::abs((*children)->pdg_id())) != decayProducts.end())
     {
       if (Verbosity() >= VERBOSITY_MAX)
       {
         std::cout << "This is a child you were looking for" << std::endl;
       }
       // Check if this is an internediate decay that didnt decay in the generator
       std::vector<int> positive_intermediates_ID;
       for (int i : m_intermediates_ID)
       {
         positive_intermediates_ID.push_back(abs(i));
       }
 
       if (!(*children)->end_vertex() && std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
                                                   abs((*children)->pdg_id())) != positive_intermediates_ID.end())
       {
         std::vector<int> requiredIntermediateDecayProducts;
         std::vector<int> positive_requiredIntermediateDecayProducts;
         std::vector<int> actualIntermediateDecayProducts;
         // Which intermediate decay list to we need
         auto it = std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(), abs((*children)->pdg_id()));
         int index = std::distance(positive_intermediates_ID.begin(), it);
 
         unsigned int trackStart = 0, trackStop = 0;
         if (index == 0)
         {
           trackStop = m_nTracksFromIntermediates[0];
         }
         else
         {
           for (int i = 0; i < index; ++i)
           {
             trackStart += m_nTracksFromIntermediates[i];
           }
           trackStop = trackStart + m_nTracksFromIntermediates[index];
         }
         for (unsigned int i = trackStart; i < trackStop; ++i)
         {
           requiredIntermediateDecayProducts.push_back(m_daughters_ID[i]);
           positive_requiredIntermediateDecayProducts.push_back(abs(m_daughters_ID[i]));
         }
 
         searchGeant4Record((*children)->barcode(), (*children)->pdg_id(), positive_requiredIntermediateDecayProducts,
                            breakLoop, hasPhoton, hasPi0, failedPT, failedETA, actualIntermediateDecayProducts);
 
         bool needThisParticle = compareDecays(requiredIntermediateDecayProducts, actualIntermediateDecayProducts);
         if (needThisParticle)
         {
           actualDecayProducts.push_back((*children)->pdg_id());
           decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*children)->barcode()), (*children)->pdg_id()));
         }
         else  // Remove what I added to the decayChain
         {
           decayChain.erase(decayChain.end() - (m_intermediate_product_counter + 1), decayChain.end());
           decayChain.pop_back();
         }
 
         m_intermediate_product_counter = 0;
       }
       else
       {
         bool needThisParticle = checkIfCorrectHepMCParticle((*children), failedPT, failedETA);
         if (needThisParticle)
         {
           actualDecayProducts.push_back((*children)->pdg_id());
           decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*children)->barcode()), (*children)->pdg_id()));
         }
       }
     }  // Now check if it's part of the resonance list
     else if (std::find(std::begin(listOfResonantPIDs), std::end(listOfResonantPIDs),
                        std::abs((*children)->pdg_id())) != std::end(listOfResonantPIDs))
     {
       if (Verbosity() >= VERBOSITY_MAX)
       {
         std::cout << "This is a resonance to investigate further" << std::endl;
       }
       if (!(*children)->end_vertex())
       {
         searchGeant4Record((*children)->barcode(), (*children)->pdg_id(), decayProducts,
                            breakLoop, hasPhoton, hasPi0, failedPT, failedETA, actualDecayProducts);
       }
       else
       {
         for (HepMC::GenVertex::particle_iterator grandchildren = (*children)->end_vertex()->particles_begin(HepMC::children);
              grandchildren != (*children)->end_vertex()->particles_end(HepMC::children); ++grandchildren)
         {
           bool needThisParticle = checkIfCorrectHepMCParticle((*grandchildren), failedPT, failedETA);
           if (needThisParticle)
           {
             actualDecayProducts.push_back((*grandchildren)->pdg_id());
             decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*grandchildren)->barcode()), (*grandchildren)->pdg_id()));
           }
         }
       }
     }
     else
     {
       breakLoop = true;  // This particle is not in the decay descriptor, stop
       break;
     }
   }
 }
 
 /*
  *Function to search Geant4 record
  * Cannot switch to HepMC search (doesn't make sense to)
  */
 void DecayFinder::searchGeant4Record(int barcode, int pid, std::vector<int> decayProducts, bool& breakLoop, bool& hasPhoton, bool& hasPi0, bool& failedPT, bool& failedETA, std::vector<int>& actualDecayProducts)
 {
   PHG4TruthInfoContainer::ConstRange range = m_truthinfo->GetParticleRange();
   for (PHG4TruthInfoContainer::ConstIterator iter = range.first; iter != range.second; ++iter)
   {
     if (decayChain.size() == 100)
     {
       breakLoop = true;  // Stuck in loop. Sympton not cause!
       break;
     }
 
     PHG4Particle* g4particle = iter->second;
     PHG4Particle* mother = nullptr;
     if (g4particle->get_parent_id() != 0)
     {
       mother = m_truthinfo->GetParticle(g4particle->get_parent_id());
     }
     else
     {
       continue;
     }
     if (mother->get_barcode() == barcode && abs(mother->get_pid()) == abs(pid))
     {
       int particleID = g4particle->get_pid();
       if (Verbosity() >= VERBOSITY_MAX)
       {
         std::cout << "--children->pdg_id(): " << particleID << std::endl;
       }
       if (particleID == 22)
       {
         hasPhoton = true;
       }
       if (particleID == 111)
       {
         hasPi0 = true;
       }
 
       if ((!m_allowPhotons && particleID == 22) || (!m_hasPhotonDaughter && particleID == 22) || (!m_allowPi0 && particleID == 111))
       {
         breakLoop = true;
         break;
       }
       // This is one of the children we are looking for
       if (std::find(decayProducts.begin(), decayProducts.end(),
                     std::abs(particleID)) != decayProducts.end())
       {
         bool needThisParticle = checkIfCorrectGeant4Particle(g4particle, hasPhoton, hasPi0, failedPT, failedETA);
         if (needThisParticle)
         {
           if (Verbosity() >= VERBOSITY_MAX)
           {
             std::cout << "This is a child you were looking for" << std::endl;
           }
           actualDecayProducts.push_back(particleID);
           int embedding_id = m_geneventmap ? m_genevt->get_embedding_id() : g4particle->get_primary_id();
           decayChain.push_back(std::make_pair(std::make_pair(embedding_id, g4particle->get_barcode()), particleID));
         }
       }  // Now check if it's part of the other resonance list
       else if ((m_allowPhotons && !m_hasPhotonDaughter && particleID == 22) || (m_allowPi0 && particleID == 111))
       {
         continue;
       }
       else if (std::find(std::begin(listOfResonantPIDs), std::end(listOfResonantPIDs),
                          std::abs(particleID)) != std::end(listOfResonantPIDs))
       {
         if (Verbosity() >= VERBOSITY_MAX)
         {
           std::cout << "This is a resonance to investigate further" << std::endl;
         }
         searchGeant4Record(g4particle->get_barcode(), g4particle->get_pid(), decayProducts,
                            breakLoop, hasPhoton, hasPi0, failedPT, failedETA, actualDecayProducts);
       }
       else
       {
         breakLoop = true;  // This particle is not in the decay descriptor, stop
         break;
       }
     }
   }
 }
 
 /*
  * Checks if this prticle really matches your request
  */
 bool DecayFinder::checkIfCorrectHepMCParticle(HepMC::GenParticle* particle, bool& trackFailedPT, bool& trackFailedETA)
 {
   bool acceptParticle = false;
 
   std::vector<int> positive_intermediates_ID;
   for (int i : m_intermediates_ID)
   {
     positive_intermediates_ID.push_back(abs(i));
   }
 
   // Check if it is an intermediate or a final track
   if (std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
                 abs(particle->pdg_id())) != positive_intermediates_ID.end())
   {
     std::vector<int> requiredIntermediateDecayProducts;
     std::vector<int> actualIntermediateDecayProducts;
     // Which intermediate decay list to we need
     auto it = std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(), abs(particle->pdg_id()));
     int index = std::distance(positive_intermediates_ID.begin(), it);
 
     unsigned int trackStart = 0, trackStop = 0;
     if (index == 0)
     {
       trackStop = m_nTracksFromIntermediates[0];
     }
     else
     {
       for (int i = 0; i < index; ++i)
       {
         trackStart += m_nTracksFromIntermediates[i];
       }
       trackStop = trackStart + m_nTracksFromIntermediates[index];
     }
 
     for (unsigned int i = trackStart; i < trackStop; ++i)
     {
       requiredIntermediateDecayProducts.push_back(m_daughters_ID[i]);
     }
 
     for (HepMC::GenVertex::particle_iterator grandchildren = particle->end_vertex()->particles_begin(HepMC::children);
          grandchildren != particle->end_vertex()->particles_end(HepMC::children); ++grandchildren)
     {
       if (Verbosity() >= VERBOSITY_MAX)
       {
         std::cout << "----grandchildren->pdg_id(): " << (*grandchildren)->pdg_id() << std::endl;
       }
 
       if (std::find(std::begin(listOfResonantPIDs), std::end(listOfResonantPIDs),
                     std::abs((*grandchildren)->pdg_id())) != std::end(listOfResonantPIDs))
       {
         for (HepMC::GenVertex::particle_iterator greatgrandchildren = (*grandchildren)->end_vertex()->particles_begin(HepMC::children);
              greatgrandchildren != (*grandchildren)->end_vertex()->particles_end(HepMC::children); ++greatgrandchildren)
         {
           if (Verbosity() >= VERBOSITY_MAX)
           {
             std::cout << "--------greatgrandchildren->pdg_id(): " << (*greatgrandchildren)->pdg_id() << std::endl;
           }
 
           if ((m_allowPhotons && !m_hasPhotonDaughter && (*greatgrandchildren)->pdg_id() == 22) || (m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111))
           {
             continue;
           }
           else if ((!m_hasPhotonDaughter && (*greatgrandchildren)->pdg_id() == 22) || (!m_allowPhotons && (*greatgrandchildren)->pdg_id() == 22) || (!m_allowPi0 && (*greatgrandchildren)->pdg_id() == 111))
           {
             break;
           }
           else
           {
             actualIntermediateDecayProducts.push_back((*greatgrandchildren)->pdg_id());
             decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*greatgrandchildren)->barcode()), (*greatgrandchildren)->pdg_id()));
             ++m_intermediate_product_counter;
 
             HepMC::FourVector myFourVector = (*greatgrandchildren)->momentum();
 
             HepMC::GenVertex* thisVtx = (*greatgrandchildren)->production_vertex();
             double vtxPos[3] = {thisVtx->point3d().x(), thisVtx->point3d().y(), thisVtx->point3d().z()};
             if (m_recalcualteEtaRange)
             {
               recalculateEta(myFourVector.py(), vtxPos);
             }
 
             if (myFourVector.perp() < m_pt_req)
             {
               trackFailedPT = true;
             }
             if (!isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req))
             {
               trackFailedETA = true;
             }
 
             if (Verbosity() >= VERBOSITY_MAX)
             {
               std::cout << "pT = " << myFourVector.perp() << ", eta = " << myFourVector.eta() << std::endl;
               std::cout << "The track " << passOrFail(myFourVector.perp() >= m_pt_req) << " the pT requirement, ";
               std::cout << "the track " << passOrFail(isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req)) << " the eta requirement." << std::endl;
             }
           }
         }
       }
       else if ((m_allowPhotons && !m_hasPhotonDaughter && (*grandchildren)->pdg_id() == 22) || (m_allowPi0 && (*grandchildren)->pdg_id() == 111))
       {
         continue;
       }
       else if ((!m_hasPhotonDaughter && (*grandchildren)->pdg_id() == 22) || (!m_allowPhotons && (*grandchildren)->pdg_id() == 22) || (!m_allowPi0 && (*grandchildren)->pdg_id() == 111))
       {
         break;
       }
       else
       {
         actualIntermediateDecayProducts.push_back((*grandchildren)->pdg_id());
         decayChain.push_back(std::make_pair(std::make_pair(m_genevt->get_embedding_id(), (*grandchildren)->barcode()), (*grandchildren)->pdg_id()));
         ++m_intermediate_product_counter;
 
         HepMC::FourVector myFourVector = (*grandchildren)->momentum();
 
         HepMC::GenVertex* thisVtx = (*grandchildren)->production_vertex();
         double vtxPos[3] = {thisVtx->point3d().x(), thisVtx->point3d().y(), thisVtx->point3d().z()};
         if (m_recalcualteEtaRange)
         {
           recalculateEta(myFourVector.py(), vtxPos);
         }
 
         if (myFourVector.perp() < m_pt_req)
         {
           trackFailedPT = true;
         }
         if (!isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req))
         {
           trackFailedETA = true;
         }
 
         if (Verbosity() >= VERBOSITY_MAX)
         {
           std::cout << "pT = " << myFourVector.perp() << ", eta = " << myFourVector.eta() << std::endl;
           std::cout << "The track " << passOrFail(myFourVector.perp() >= m_pt_req) << " the pT requirement, ";
           std::cout << "the track " << passOrFail(isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req)) << " the eta requirement." << std::endl;
         }
       }
     }
 
     acceptParticle = compareDecays(requiredIntermediateDecayProducts, actualIntermediateDecayProducts);
   }
   //else if ((particle->pdg_id() == 22) || (particle->pdg_id() == 111))
   //{
   //  return false;
   //}
   else
   {
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "This is a final state track" << std::endl;
     }
     const HepMC::FourVector& myFourVector = particle->momentum();
 
     HepMC::GenVertex* thisVtx = particle->production_vertex();
     double vtxPos[3] = {thisVtx->point3d().x(), thisVtx->point3d().y(), thisVtx->point3d().z()};
     if (m_recalcualteEtaRange)
     {
       recalculateEta(myFourVector.py(), vtxPos);
     }
 
     if (myFourVector.perp() < m_pt_req)
     {
       trackFailedPT = true;
     }
     if (!isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req))
     {
       trackFailedETA = true;
     }
 
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "pT = " << myFourVector.perp() << ", eta = " << myFourVector.eta() << std::endl;
       std::cout << "The track " << passOrFail(myFourVector.perp() >= m_pt_req) << " the pT requirement, ";
       std::cout << "the track " << passOrFail(isInRange(m_eta_low_req, myFourVector.eta(), m_eta_high_req)) << " the eta requirement." << std::endl;
     }
     acceptParticle = true;
   }
 
   return acceptParticle;
 }
 
 /*
  * Checks if this prticle really matches your request
  */
 bool DecayFinder::checkIfCorrectGeant4Particle(PHG4Particle* particle, bool& hasPhoton, bool& hasPi0, bool& trackFailedPT, bool& trackFailedETA)
 {
   bool acceptParticle = false;
 
   std::vector<int> positive_intermediates_ID;
   for (int i : m_intermediates_ID)
   {
     positive_intermediates_ID.push_back(abs(i));
   }
   // Check if it is an intermediate or a final track
   if (std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(),
                 abs(particle->get_pid())) != positive_intermediates_ID.end())
   {
     std::vector<int> requiredIntermediateDecayProducts;
     std::vector<int> actualIntermediateDecayProducts;
     // Which intermediate decay list to we need
     auto it = std::find(positive_intermediates_ID.begin(), positive_intermediates_ID.end(), abs(particle->get_pid()));
     int index = std::distance(positive_intermediates_ID.begin(), it);
 
     unsigned int trackStart = 0, trackStop = 0;
     if (index == 0)
     {
       trackStop = m_nTracksFromIntermediates[0];
     }
     else
     {
       for (int i = 0; i < index; ++i)
       {
         trackStart += m_nTracksFromIntermediates[i];
       }
       trackStop = trackStart + m_nTracksFromIntermediates[index];
     }
 
     std::vector<int> positive_intermediateDecayProducts;
     for (unsigned int i = trackStart; i < trackStop; ++i)
     {
       positive_intermediateDecayProducts.push_back(abs(m_daughters_ID[i]));
       requiredIntermediateDecayProducts.push_back(m_daughters_ID[i]);
     }
 
     bool fakeBreak = false;
     searchGeant4Record(particle->get_barcode(), particle->get_pid(), positive_intermediateDecayProducts, fakeBreak,
                        hasPhoton, hasPi0, trackFailedPT, trackFailedETA, actualIntermediateDecayProducts);
 
     acceptParticle = compareDecays(requiredIntermediateDecayProducts, actualIntermediateDecayProducts);
   }
   //else if ((particle->get_pid() == 22) || (particle->get_pid() == 111))
   //{
   //  return false;
   //}
   else
   {
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "This is a final state track" << std::endl;
     }
     double px = particle->get_px();
     double py = particle->get_py();
     double pz = particle->get_pz();
 
     double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
     double p = std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2));
     double eta = 0.5 * std::log((p + pz) / (p - pz));
 
     PHG4VtxPoint* thisVtx = m_truthinfo->GetVtx(particle->get_vtx_id());
     double vtxPos[3] = {thisVtx->get_x(), thisVtx->get_y(), thisVtx->get_z()};
     if (m_recalcualteEtaRange)
     {
       recalculateEta(py, vtxPos);
     }
 
     if (pt < m_pt_req)
     {
       trackFailedPT = true;
     }
     if (!isInRange(m_eta_low_req, eta, m_eta_high_req))
     {
       trackFailedETA = true;
     }
 
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "pT = " << pt << ", eta = " << eta << std::endl;
       std::cout << "The track " << passOrFail(pt >= m_pt_req) << " the pT requirement, ";
       std::cout << "the track " << passOrFail(isInRange(m_eta_low_req, eta, m_eta_high_req)) << " the eta requirement." << std::endl;
     }
 
     acceptParticle = true;
   }
 
   return acceptParticle;
 }
 
 /*
  * Check if the decay matches what we need
  */
 bool DecayFinder::compareDecays(std::vector<int> required, std::vector<int> actual)
 {
   bool accept = false;
   multiplyVectorByScalarAndSort(required, +1);
   multiplyVectorByScalarAndSort(actual, +1);
 
   if (Verbosity() >= VERBOSITY_MAX)
   {
     std::cout << "Printing required decay products: ";
     for (int i : required)
     {
       std::cout << i << ", ";
     }
     std::cout << std::endl;
     std::cout << "Printing actual decay products: ";
     for (int i : actual)
     {
       std::cout << i << ", ";
     }
     std::cout << std::endl;
     if (required == actual)
     {
       std::cout << "*\n* These vectors match\n*\n"
                 << std::endl;
     }
     else
     {
       std::cout << "*\n* These vectors DONT match\n*\n"
                 << std::endl;
     }
   }
 
   if (required == actual)
   {
     accept = true;
   }
 
   if (m_getChargeConjugate && !accept)
   {
     multiplyVectorByScalarAndSort(required, -1);
     if (required == actual)
     {
       accept = true;
     }
 
     if (Verbosity() >= VERBOSITY_MAX)
     {
       std::cout << "Checking CC state" << std::endl;
       if (required == actual)
       {
         std::cout << "*\n* These vectors match\n*\n"
                   << std::endl;
       }
       else
       {
         std::cout << "*\n* These vectors DONT match\n*\n"
                   << std::endl;
       }
     }
   }
   return accept;
 }
 
 /*
  * Everything below this is helper functions
  */
 int DecayFinder::deleteElement(int arr[], int n, int x)
 {
   // https://www.geeksforgeeks.org/delete-an-element-from-array-using-two-traversals-and-one-traversal/
   // Search x in array
   int i;
   for (i = 0; i < n; i++)
   {
     if (arr[i] == x)
     {
       break;
     }
   }
 
   // If x found in array
   if (i < n)
   {
     // reduce size of array and move all
     // elements on space ahead
     n = n - 1;
     for (int j = i; j < n; j++)
     {
       arr[j] = arr[j + 1];
     }
   }
 
   return n;
 }
 
 bool DecayFinder::findParticle(const std::string& particle)
 {
   bool particleFound = true;
   if (!TDatabasePDG::Instance()->GetParticle(particle.c_str()))
   {
     if (Verbosity() >= VERBOSITY_SOME)
     {
       std::cout << "The particle, " << particle << " is not in the TDatabasePDG particle list" << std::endl;
     }
     particleFound = false;
   }
 
   return particleFound;
 }
 
 void DecayFinder::multiplyVectorByScalarAndSort(std::vector<int>& v, int k)
 {
   // https://slaystudy.com/c-multiply-vector-by-scalar/
   std::transform(v.begin(), v.end(), v.begin(), [k](const int& c)
   {
     int particlesWithNoCC[] = {111, 113, 115, 130, 220, 221, 223, 225, 310, 330, 331, 333, 335, 440, 441, 443, 
                                445, 551, 553, 555, 10111, 10113, 10221, 10223, 10331, 10333, 10441, 10443, 
                                10551, 10553, 20113, 20223, 20333, 20443, 20553, 100443, 100553};
     if (std::find(std::begin(particlesWithNoCC), std::end(particlesWithNoCC), c) != std::end(particlesWithNoCC))
     {
       return c;
     }
     else
     {
       return c * k;
     }
   });
   std::sort(v.begin(), v.end());
 }
 
 int DecayFinder::get_pdgcode(const std::string& name)
 {
   if (findParticle(name))
   {
     return TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
   }
   else
   {
     return 0;
   }
 }
 
 int DecayFinder::get_charge(const std::string& name)
 {
   if (findParticle(name))
   {
     return TDatabasePDG::Instance()->GetParticle(name.c_str())->Charge() / 3;
   }
   else
   {
     return -99;
   }
 }
 
 bool DecayFinder::isInRange(float min, float value, float max)
 {
   return min <= value && value <= max;
 }
 
 void DecayFinder::calculateEffectiveTPCradius(double vertex[3], double& effective_top_r, double& effective_bottom_r)
 {
   double r_sq = std::pow(m_tpc_r, 2);
   double x_sq = std::pow(vertex[0], 2);
 
   effective_top_r = std::sqrt(r_sq - x_sq) - vertex[1];
   effective_bottom_r = std::sqrt(r_sq - x_sq) + vertex[1];
 }
 
 void DecayFinder::recalculateEta(double py, double vertex[3])
 {
   calculateEffectiveTPCradius(vertex, m_effective_top_tpc_r, m_effective_bottom_tpc_r);
 
   double z_north = m_tpc_z - vertex[2];  // Distance between north end of TPC and vertex z position
   double z_south = m_tpc_z + vertex[2];  // Same but for south of TPC
 
   double top_north_theta = atan(m_effective_top_tpc_r / z_north);
   double bottom_north_theta = atan(m_effective_bottom_tpc_r / z_north);
   double top_south_theta = atan(m_effective_top_tpc_r / z_south);
   double bottom_south_theta = atan(m_effective_bottom_tpc_r / z_south);
 
   double top_north_eta = -1 * std::log(tan(top_north_theta / 2));
   double bottom_north_eta = -1 * std::log(tan(bottom_north_theta / 2));
   double top_south_eta = std::log(tan(top_south_theta / 2));
   double bottom_south_eta = std::log(tan(bottom_south_theta / 2));
 
   m_eta_high_req = py < 0 ? bottom_north_eta : top_north_eta;
   m_eta_low_req = py < 0 ? bottom_south_eta : top_south_eta;
 }
 
 int DecayFinder::createDecayNode(PHCompositeNode* topNode)
 {
   PHNodeIterator iter(topNode);
 
   PHCompositeNode* lowerNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
   if (!lowerNode)
   {
     lowerNode = new PHCompositeNode("DST");
     topNode->addNode(lowerNode);
     std::cout << "DST node added" << std::endl;
   }
 
   std::string baseName;
 
   if (m_container_name.empty())
   {
     // baseName = "decay";
     baseName = Name();
   }
   else
   {
     baseName = m_container_name;
   }
 
   // Cant have forward slashes in DST or else you make a subdirectory on save!!!
   size_t pos;
   std::string undrscr = "_";
   std::string nothing = "";
   std::map<std::string, std::string> forbiddenStrings;
   forbiddenStrings["/"] = undrscr;
   forbiddenStrings["("] = undrscr;
   forbiddenStrings[")"] = nothing;
   forbiddenStrings["+"] = "plus";
   forbiddenStrings["-"] = "minus";
   forbiddenStrings["*"] = "star";
   forbiddenStrings[" "] = nothing;
   for (auto const& [badString, goodString] : forbiddenStrings)
   {
     while ((pos = baseName.find(badString)) != std::string::npos)
     {
       baseName.replace(pos, 1, goodString);
     }
   }
 
   m_nodeName = baseName + "_DecayMap";
 
   m_decayMap = new DecayFinderContainer_v1();
   PHIODataNode<PHObject>* decayNode = new PHIODataNode<PHObject>(m_decayMap, m_nodeName.c_str(), "PHObject");
   lowerNode->addNode(decayNode);
   std::cout << m_nodeName << " node added" << std::endl;
 
   return Fun4AllReturnCodes::EVENT_OK;
 }
 
 void DecayFinder::fillDecayNode(PHCompositeNode* topNode, Decay& decay)
 {
   m_decayMap = findNode::getClass<DecayFinderContainer_v1>(topNode, m_nodeName.c_str());
   m_decayMap->insert(decay);
 }
 
 void DecayFinder::printInfo()
 {
   std::cout << "\n---------------DecayFinder information---------------" << std::endl;
   std::cout << "Module name: " << Name() << std::endl;
   std::cout << "Decay descriptor: " << m_decayDescriptor << std::endl;
   std::cout << "Number of generated decays: " << m_counter << std::endl;
   std::cout << "  Number of decays that failed pT requirement: " << m_nCandFail_pT << std::endl;
   std::cout << "  Number of decays that failed eta requirement: " << m_nCandFail_eta << std::endl;
   std::cout << "  Number of decays that failed pT and eta requirements: " << m_nCandFail_pT_and_eta << std::endl;
   std::cout << "Number of decays that could be reconstructed: " << m_nCandReconstructable << std::endl;
   std::cout << "  Number of reconstructable mothers with associated photon: " << m_nCandHas_Photon << std::endl;
   std::cout << "  Number of reconstructable mothers with associated Pi0: " << m_nCandHas_Pi0 << std::endl;
   std::cout << "  Number of reconstructable mothers with associated photon and Pi0: " << m_nCandHas_Photon_and_Pi0 << std::endl;
   std::cout << "  Number of reconstructable mothers with no associated photons or Pi0s: " << m_nCandHas_noPhoton_and_noPi0 << std::endl;
   std::cout << "-----------------------------------------------------\n";
   std::cout << std::endl;
 }
 
 void DecayFinder::printNode(PHCompositeNode* topNode)
 {
   std::cout << "----------------";
   std::cout << " DecayFinderNode: " << m_nodeName << " information ";
   std::cout << "----------------" << std::endl;
   DecayFinderContainer_v1* map = findNode::getClass<DecayFinderContainer_v1>(topNode, m_nodeName.c_str());
   for (auto& iter : *map)
   {
     Decay decay = iter.second;
     for (auto& i : decay)
     {
       std::cout << "Particle embedding ID: " << i.first.first << ", particle barcode: " << i.first.second << ", particle PDG ID: " << i.second << std::endl;
     }
   }
   std::cout << "--------------------------------------------------------------------------------------------------" << std::endl;
 }
 
 std::string DecayFinder::passOrFail(bool condition)
 {
   return condition ? "passed" : "failed";
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team