MCShowerRecoPart.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
//  MCShowerRecoPart source
//
////////////////////////////////////////////////////////////////////////

#include "fhiclcpp/ParameterSet.h"

#include "MCShowerRecoPart.h"
#include "larsim/MCSTReco/MCRecoPart.h"

namespace sim {

  const unsigned int MCShowerRecoPart::kINVALID_UINT = std::numeric_limits<unsigned int>::max();
  const int MCShowerRecoPart::kINVALID_INT = std::numeric_limits<int>::max();

  //##################################################################
  MCShowerRecoPart::MCShowerRecoPart(fhicl::ParameterSet const& pset)
  //##################################################################
  {
    _debug_mode = pset.get<bool>("DebugMode");
  }

  //-----------------------------------------------------------------------------------------------
  void MCShowerRecoPart::ConstructShower(const MCRecoPart& part_v)
  //-----------------------------------------------------------------------------------------------
  {
    _shower_id.clear();
    _shower_id.resize(part_v.size(),-1);
    _shower_index.clear();
    _shower_daughters.clear();

    if(!part_v.size()) return;

    // Construct MCShower
    std::vector<std::multimap<double,unsigned int> > daughter_map;
    for(size_t i=0; i<part_v.size(); ++i) {

      auto const& mcp = part_v[i];

      int candidate_mom_index=-1;
      if( mcp._pdgcode == 22 ||
          mcp._pdgcode == 11 ||
          mcp._pdgcode == -11 )
        candidate_mom_index = i;

      unsigned int mom_track = mcp._mother;
      auto mom_iter = part_v._track_index.find(mom_track);
      while(mom_iter != part_v._track_index.end()) {

        unsigned int mom_index = (*mom_iter).second;

        if( part_v.at(mom_index)._pdgcode == 22 || part_v.at(mom_index)._pdgcode == 11 || part_v.at(mom_index)._pdgcode == -11 )

          candidate_mom_index = mom_index;

        mom_iter = part_v._track_index.find(part_v.at(mom_index)._mother);

      }

      if(candidate_mom_index >= 0) {

        auto candidate_mom_iter = _shower_index.find(candidate_mom_index);
        if(candidate_mom_iter == _shower_index.end()) {
          _shower_index.insert(std::make_pair((unsigned int)candidate_mom_index, (unsigned int)_shower_index.size()));
          daughter_map.push_back(std::multimap<double,unsigned int>());
        }
        unsigned int shower_index = (*_shower_index.find(candidate_mom_index)).second;
        daughter_map.at(shower_index).insert(std::make_pair((double)(mcp._start_vtx[3]),(unsigned int)i));
        _shower_id.at(i) = shower_index;

      } else if(_debug_mode) {

        std::cout
          << "Found a particle that does not belong to a shower!" << std::endl
          << Form(" PDGID: %d ... Track %d @ (%g,%g,%g,%g) with (%g,%g,%g,%g)",
                  mcp._pdgcode,
                  mcp._track_id,
                  mcp._start_vtx[0],
                  mcp._start_vtx[1],
                  mcp._start_vtx[2],
                  mcp._start_vtx[3],
                  mcp._start_mom[0],
                  mcp._start_mom[1],
                  mcp._start_mom[2],
                  mcp._start_mom[3])
          << std::endl << std::endl;

      }
    }


    if(_debug_mode)
      std::cout
        << Form("Found %zu MCShowers....",_shower_index.size()) << std::endl;

    _shower_daughters.resize(_shower_index.size(),std::vector<unsigned int>());
    for(auto const &mom : _shower_index) {

      _shower_daughters.at(mom.second).reserve(daughter_map.at(mom.second).size());
      for(auto const &part_index : daughter_map.at(mom.second))

        _shower_daughters.at(mom.second).push_back(part_index.second);

      auto const& mcp = part_v.at(mom.first);
      if(_debug_mode)
        std::cout
          << Form("PDGID: %d ... Track %d @ (%g,%g,%g,%g) with (%g,%g,%g,%g) ... %zu daughters!",
                  mcp._pdgcode,
                  mcp._track_id,
                  mcp._start_vtx[0],
                  mcp._start_vtx[1],
                  mcp._start_vtx[2],
                  mcp._start_vtx[3],
                  mcp._start_mom[0],
                  mcp._start_mom[1],
                  mcp._start_mom[2],
                  mcp._start_mom[3],
                  _shower_daughters.at(mom.second).size())
          << std::endl;
    }

    if(_debug_mode)
      std::cout<<std::endl;
  }

}