MCShowerRecoAlg.cxx

Go to the documentation of this file.

//
//  MCShowerRecoAlg source
//
//

#include "MCShowerRecoAlg.h"

#include "RtypesCore.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TString.h"
#include "TVector3.h"

#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Utilities/Exception.h"
#include "fhiclcpp/ParameterSet.h"

#include "larcore/Geometry/Geometry.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/MCBase/MCLimits.h"
#include "lardataobj/MCBase/MCShower.h"
#include "lardataobj/MCBase/MCStep.h"
#include "larsim/MCSTReco/MCRecoEdep.h"
#include "larsim/MCSTReco/MCRecoPart.h"
#include "larsim/MCSTReco/MCShowerRecoPart.h"

#include <memory>
#include <vector>

namespace sim {

  //##################################################################
  MCShowerRecoAlg::MCShowerRecoAlg(fhicl::ParameterSet const& pset)
    : fPartAlg(pset.get<fhicl::ParameterSet>("MCShowerRecoPart"))
    , fDebugMode(pset.get<bool>("DebugMode"))
    , fMinShowerEnergy(pset.get<double>("MinShowerEnergy"))
    , fMinNumDaughters(pset.get<unsigned int>("MinNumDaughters"))
  //##################################################################
  {}

  std::unique_ptr<std::vector<sim::MCShower>> MCShowerRecoAlg::Reconstruct(MCRecoPart& part_v,
                                                                           MCRecoEdep& edep_v)
  {

    art::ServiceHandle<geo::Geometry const> geo;

    auto pindex = details::createPlaneIndexMap();

    fPartAlg.ConstructShower(part_v);
    auto result = std::make_unique<std::vector<sim::MCShower>>();
    auto& mcshower = *result;
    //std::vector<sim::MCShower> mcshower;
    // Get shower info from grouped particles
    const std::vector<unsigned int> shower_index_v = fPartAlg.ShowerMothers();
    mcshower.reserve(shower_index_v.size());
    std::vector<size_t> mcs_to_spart_v;
    mcs_to_spart_v.reserve(shower_index_v.size());

    bool daughter_stored = false;
    for (size_t shower_index = 0; shower_index < shower_index_v.size(); ++shower_index) {

      unsigned int shower_candidate = shower_index_v.at(shower_index);
      auto const& shower_part = part_v.at(shower_candidate);

      unsigned int mother_track = part_v.MotherTrackID(shower_candidate);
      unsigned int ancestor_track = part_v.AncestorTrackID(shower_candidate);

      if (mother_track == kINVALID_UINT || ancestor_track == kINVALID_UINT)

        throw cet::exception(__FUNCTION__) << "LOGIC ERROR: mother/ancestor track ID is invalid!";

      MCMiniPart mother_part;
      MCMiniPart ancestor_part;

      unsigned int mother_index = part_v.TrackToParticleIndex(mother_track);
      unsigned int ancestor_index = part_v.TrackToParticleIndex(ancestor_track);

      if (mother_index != kINVALID_UINT)
        mother_part = part_v[mother_index];
      else
        mother_part._track_id = mother_track;

      if (ancestor_index != kINVALID_UINT)
        ancestor_part = part_v[ancestor_index];
      else
        ancestor_part._track_id = ancestor_track;

      double shower_g4_energy = shower_part._start_mom[3];

      if (fDebugMode)

        std::cout << "Found MCShower with mother energy: " << shower_g4_energy << " MeV";

      // Skip if mother energy is less than the enery threshold
      if (shower_g4_energy < fMinShowerEnergy) {
        if (fDebugMode) std::cout << " ... below energy threshold: skipping!" << std::endl;

        continue;
      }
      else if (shower_part._daughters.size() < fMinNumDaughters) {
        if (fDebugMode)
          std::cout << " ... below # daughter particle count threshold: skipping!" << std::endl;

        continue;
      }
      else if (fDebugMode) {
        std::cout << " ... condition matched. Storing this MCShower..." << std::endl;
      }

      // Record this MCShower
      mcs_to_spart_v.push_back(shower_index);

      if (fDebugMode)

        std::cout << " Storage index " << mcshower.size() << " => Shower index " << shower_index
                  << std::endl;

      ::sim::MCShower shower_prof;

      shower_prof.Origin(shower_part._origin);
      shower_prof.PdgCode(shower_part._pdgcode);
      shower_prof.TrackID(shower_part._track_id);
      shower_prof.Process(shower_part._process);

      shower_prof.MotherPdgCode(mother_part._pdgcode);
      shower_prof.MotherTrackID(mother_part._track_id);
      shower_prof.MotherProcess(mother_part._process);

      shower_prof.AncestorPdgCode(ancestor_part._pdgcode);
      shower_prof.AncestorTrackID(ancestor_part._track_id);
      shower_prof.AncestorProcess(ancestor_part._process);

      shower_prof.Start(MCStep(shower_part._start_vtx, shower_part._start_mom));
      shower_prof.End(MCStep(shower_part._end_vtx, shower_part._end_mom));
      shower_prof.MotherStart(MCStep(mother_part._start_vtx, mother_part._start_mom));
      shower_prof.MotherEnd(MCStep(mother_part._end_vtx, mother_part._end_mom));
      shower_prof.AncestorStart(MCStep(ancestor_part._start_vtx, ancestor_part._start_mom));
      shower_prof.AncestorEnd(MCStep(ancestor_part._end_vtx, ancestor_part._end_mom));

      // Daughter list
      std::vector<unsigned int> daughter_track_id;
      daughter_track_id.reserve(fPartAlg.ShowerDaughters(shower_index).size());

      for (auto const& index : fPartAlg.ShowerDaughters(shower_index))

        daughter_track_id.push_back(part_v.at(index)._track_id);

      shower_prof.DaughterTrackID(daughter_track_id);

      if (!daughter_stored && daughter_track_id.size() > 1) daughter_stored = true;

      mcshower.push_back(shower_prof);
    }

    if (fDebugMode)
      std::cout << " Found " << mcshower.size()
                << " MCShowers. Now computing DetProfile position..." << std::endl;

    //
    // Daughter vtx
    //
    std::vector<TLorentzVector> mcs_daughter_vtx_v(
      mcshower.size(),
      TLorentzVector(
        sim::kINVALID_DOUBLE, sim::kINVALID_DOUBLE, sim::kINVALID_DOUBLE, sim::kINVALID_DOUBLE));
    std::vector<TLorentzVector> mcs_daughter_mom_v(mcshower.size(), TLorentzVector());

    std::vector<std::vector<double>> plane_charge_v(mcshower.size(), std::vector<double>(3, 0));
    std::vector<std::vector<double>> plane_dqdx_v(mcshower.size(), std::vector<double>(3, 0));

    //For dEdx Calculation
    std::vector<double> mcs_daughter_dedx_v(mcshower.size(), 0);
    std::vector<double> mcs_daughter_dedxRAD_v(mcshower.size(), 0);
    std::vector<TVector3> mcs_daughter_dir_v(mcshower.size(), TVector3());

    for (size_t mcs_index = 0; mcs_index < mcshower.size(); ++mcs_index) {

      auto& mcs_daughter_vtx = mcs_daughter_vtx_v[mcs_index];
      auto& mcs_daughter_mom = mcs_daughter_mom_v[mcs_index];
      auto& mcs_daughter_dedx = mcs_daughter_dedx_v[mcs_index];
      auto& mcs_daughter_dedxRAD = mcs_daughter_dedxRAD_v[mcs_index];
      auto& mcs_daughter_dir = mcs_daughter_dir_v[mcs_index];
      auto& plane_charge = plane_charge_v[mcs_index];
      auto& plane_dqdx = plane_dqdx_v[mcs_index];

      for (auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {

        auto const daughter_part_index = part_v.TrackToParticleIndex(daughter_trk_id);

        auto const& daughter_part = part_v[daughter_part_index];

        auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);

        if (daughter_edep_index < 0) continue;

        auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);

        if (!(daughter_edep.size())) continue;

        // Record first daughter's vtx point
        double min_dist = sim::kINVALID_DOUBLE;
        for (auto const& edep : daughter_edep) {

          double dist = sqrt(pow(edep.pos._x - daughter_part._start_vtx[0], 2) +
                             pow(edep.pos._y - daughter_part._start_vtx[1], 2) +
                             pow(edep.pos._z - daughter_part._start_vtx[2], 2));

          if (dist < min_dist) {
            min_dist = dist;
            mcs_daughter_vtx[0] = edep.pos._x;
            mcs_daughter_vtx[1] = edep.pos._y;
            mcs_daughter_vtx[2] = edep.pos._z;
            mcs_daughter_vtx[3] = (dist / 100. / 2.998e8) * 1.e9 + daughter_part._start_vtx[3];
          }
        }
        if (!daughter_stored) {
          // If daughter is not stored, and shower id energetic enough, attempt to include angle info
          std::vector<double> shower_dir(3, 0);
          shower_dir[0] = mcshower[mcs_index].Start().Px();
          shower_dir[1] = mcshower[mcs_index].Start().Py();
          shower_dir[2] = mcshower[mcs_index].Start().Pz();
          double magnitude = 0;
          for (size_t i = 0; i < 3; ++i)
            magnitude += pow(shower_dir[i], 2);

          magnitude = sqrt(magnitude);

          if (magnitude > 1.e-10) {
            // If enough momentum, include angle info
            min_dist = sim::kINVALID_DOUBLE;

            for (auto& v : shower_dir)
              v /= magnitude;

            for (auto const& edep : daughter_edep) {
              std::vector<double> shower_dep_dir(3, 0);
              shower_dep_dir[0] = edep.pos._x - mcshower[mcs_index].Start().X();
              shower_dep_dir[1] = edep.pos._y - mcshower[mcs_index].Start().Y();
              shower_dep_dir[2] = edep.pos._z - mcshower[mcs_index].Start().Z();

              double dist = sqrt(pow(shower_dep_dir[0], 2) + pow(shower_dep_dir[1], 2) +
                                 pow(shower_dep_dir[2], 2));
              for (auto& v : shower_dep_dir)
                v /= dist;

              double angle =
                acos(shower_dep_dir[0] * shower_dir[0] + shower_dep_dir[1] * shower_dir[1] +
                     shower_dep_dir[2] * shower_dir[2]) /
                TMath::Pi() * 180.;

              if (dist < min_dist && angle < 10) {

                min_dist = dist;
                mcs_daughter_vtx[0] = edep.pos._x;
                mcs_daughter_vtx[1] = edep.pos._y;
                mcs_daughter_vtx[2] = edep.pos._z;
                mcs_daughter_vtx[3] =
                  (dist / 100. / 2.998e8) * 1.e9 + mcshower[mcs_index].Start().T();
              }
            }
          }
        }
        break;
      }
      // Now take care of momentum & plane charge

      std::vector<double> mom(3, 0);
      for (auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {

        //auto const daughter_part_index = part_v.TrackToParticleIndex(daughter_trk_id);

        // for c2: daughter_part is unused
        //auto const& daughter_part = part_v[daughter_part_index];

        auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);

        if (daughter_edep_index < 0) continue;

        auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);

        if (!(daughter_edep.size())) continue;

        //bool first=true;  // unused
        for (auto const& edep : daughter_edep) {

          // Compute unit vector to this energy deposition
          mom[0] = edep.pos._x - mcs_daughter_vtx[0];
          mom[1] = edep.pos._y - mcs_daughter_vtx[1];
          mom[2] = edep.pos._z - mcs_daughter_vtx[2];

          // Weight by energy (momentum)
          double magnitude = sqrt(pow(mom.at(0), 2) + pow(mom.at(1), 2) + pow(mom.at(2), 2));

          double energy = 0;
          double npid = 0;
          for (auto const& pid_energy : edep.deps) {
            npid++;
            energy += pid_energy.energy;
          }
          energy /= npid;
          if (magnitude > 1.e-10) {
            mom.at(0) = mom.at(0) * energy / magnitude;
            mom.at(1) = mom.at(1) * energy / magnitude;
            mom.at(2) = mom.at(2) * energy / magnitude;
            mcs_daughter_mom[0] += mom.at(0);
            mcs_daughter_mom[1] += mom.at(1);
            mcs_daughter_mom[2] += mom.at(2);
          }
          //Determine the direction of the shower right at the start point
          double E = 0;
          double N = 0;
          if (sqrt(pow(edep.pos._x - mcs_daughter_vtx[0], 2) +
                   pow(edep.pos._y - mcs_daughter_vtx[1], 2) +
                   pow(edep.pos._z - mcs_daughter_vtx[2], 2)) < 2.4 &&
              magnitude > 1.e-10) {

            mcs_daughter_dir[0] += mom.at(0);
            mcs_daughter_dir[1] += mom.at(1);
            mcs_daughter_dir[2] += mom.at(2);
            E += energy;
            N += 1;
          }

          if (E > 0) E /= N;
          mcs_daughter_dedxRAD += E;

          mcs_daughter_mom[3] += energy;

          // Charge
          auto const pid = edep.pid;
          auto q_i = pindex.find(pid);
          if (q_i != pindex.end())
            plane_charge[pid.Plane] += (double)(edep.deps[pindex[pid]].charge);

        } 

      } 
      mcs_daughter_dedxRAD /= 2.4;

      for (auto const& daughter_trk_id : mcshower[mcs_index].DaughterTrackID()) {

        //auto const daughter_part_index = part_v.TrackToParticleIndex(daughter_trk_id);

        // for c2: daughter_part is unused
        //auto const& daughter_part = part_v[daughter_part_index];

        auto const daughter_edep_index = edep_v.TrackToEdepIndex(daughter_trk_id);

        if (daughter_edep_index < 0) continue;

        auto const& daughter_edep = edep_v.GetEdepArrayAt(daughter_edep_index);

        if (!(daughter_edep.size())) continue;

        for (auto const& edep : daughter_edep) {

          //Defining dEdx
          //Need to define a plane through the shower start point (x_0, y_0, z_0) with a normal along the momentum vector of the shower
          //The plane will be defined in the typical way:
          // a*x + b*y + c*z + d = 0
          // where, a = dir_x, b = dir_y, c = dir_z, d = - (a*x_0+b*y_0+c*z_0)
          // then the *signed* distance of any point (x_1, y_1, z_1) from this plane is:
          // D = (a*x_1 + b*y_1 + c*z_1 + d )/sqrt( pow(a,2) + pow(b,2) + pow(c,2))

          double p_mag = sqrt(pow(mcs_daughter_dir[0], 2) + pow(mcs_daughter_dir[1], 2) +
                              pow(mcs_daughter_dir[2], 2));
          double a = 0, b = 0, c = 0, d = 0;
          if (p_mag > 1.e-10) {
            a = mcs_daughter_dir[0] / p_mag;
            b = mcs_daughter_dir[1] / p_mag;
            c = mcs_daughter_dir[2] / p_mag;
            d = -1 * (a * mcs_daughter_vtx[0] + b * mcs_daughter_vtx[1] + c * mcs_daughter_vtx[2]);
          }
          else {
            mcs_daughter_dedx += 0;
            continue;
          }
          //Radial Distance
          if ((a * edep.pos._x + b * edep.pos._y + c * edep.pos._z + d) /
                  sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) <
                2.4 &&
              (a * edep.pos._x + b * edep.pos._y + c * edep.pos._z + d) /
                  sqrt(pow(a, 2) + pow(b, 2) + pow(c, 2)) >
                0) {

            double E = 0;
            double N = 0;

            for (auto const& pid_energy : edep.deps) {
              N += 1;
              E += pid_energy.energy;
            }

            if (N > 0) { E /= N; }
            else {
              E = 0;
            }

            mcs_daughter_dedx += E;

            // Charge
            auto const pid = edep.pid;
            auto q_i = pindex.find(pid);
            if (q_i != pindex.end())
              plane_dqdx[pid.Plane] += (double)(edep.deps[pindex[pid]].charge);
          }
        }
      }
      mcs_daughter_dedx /= 2.4;
      plane_dqdx.at(0) /= 2.4;
      plane_dqdx.at(1) /= 2.4;
      plane_dqdx.at(2) /= 2.4;

    } 

    if (fDebugMode)
      std::cout << " Found " << mcshower.size() << " MCShowers. Now storing..." << std::endl;

    // Store plane charge & daughter momentum
    for (size_t mcs_index = 0; mcs_index < mcshower.size(); ++mcs_index) {

      auto& daughter_vtx = mcs_daughter_vtx_v[mcs_index];
      auto& daughter_mom = mcs_daughter_mom_v[mcs_index];
      auto& daughter_dedx = mcs_daughter_dedx_v[mcs_index];
      auto& daughter_dedxRAD = mcs_daughter_dedxRAD_v[mcs_index];
      auto& daughter_dir = mcs_daughter_dir_v[mcs_index];
      auto& plane_charge = plane_charge_v[mcs_index];
      auto& plane_dqdx = plane_dqdx_v[mcs_index];

      double magnitude =
        sqrt(pow(daughter_mom[0], 2) + pow(daughter_mom[1], 2) + pow(daughter_mom[2], 2));

      if (daughter_mom[3] > 1.e-10) {
        daughter_mom[0] *= daughter_mom[3] / magnitude;
        daughter_mom[1] *= daughter_mom[3] / magnitude;
        daughter_mom[2] *= daughter_mom[3] / magnitude;
      }
      else
        for (size_t i = 0; i < 4; ++i)
          daughter_mom[i] = 0;

      mcshower.at(mcs_index).DetProfile(MCStep(daughter_vtx, daughter_mom));
      mcshower.at(mcs_index).Charge(plane_charge);
      mcshower.at(mcs_index).dQdx(plane_dqdx);
      mcshower.at(mcs_index).dEdx(daughter_dedx);
      mcshower.at(mcs_index).dEdxRAD(daughter_dedxRAD);
      mcshower.at(mcs_index).StartDir(daughter_dir);
    }

    if (fDebugMode) {

      for (auto const& prof : mcshower) {

        std::cout

          << Form("  Shower particle:     PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
                  "(%g,%g,%g,%g)",
                  prof.PdgCode(),
                  prof.TrackID(),
                  prof.Start().X(),
                  prof.Start().Y(),
                  prof.Start().Z(),
                  prof.Start().T(),
                  prof.Start().Px(),
                  prof.Start().Py(),
                  prof.Start().Pz(),
                  prof.Start().E())
          << std::endl
          << Form("    Mother particle:   PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
                  "(%g,%g,%g,%g)",
                  prof.MotherPdgCode(),
                  prof.MotherTrackID(),
                  prof.MotherStart().X(),
                  prof.MotherStart().Y(),
                  prof.MotherStart().Z(),
                  prof.MotherStart().T(),
                  prof.MotherStart().Px(),
                  prof.MotherStart().Py(),
                  prof.MotherStart().Pz(),
                  prof.MotherStart().E())
          << std::endl
          << Form("    Ancestor particle: PDG=%d : Track ID=%d Start @ (%g,%g,%g,%g) with Momentum "
                  "(%g,%g,%g,%g)",
                  prof.AncestorPdgCode(),
                  prof.AncestorTrackID(),
                  prof.AncestorStart().X(),
                  prof.AncestorStart().Y(),
                  prof.AncestorStart().Z(),
                  prof.AncestorStart().T(),
                  prof.AncestorStart().Px(),
                  prof.AncestorStart().Py(),
                  prof.AncestorStart().Pz(),
                  prof.AncestorStart().E())
          << std::endl
          << Form("    ... with %zu daughters: Start @ (%g,%g,%g,%g) with Momentum (%g,%g,%g,%g)",
                  prof.DaughterTrackID().size(),
                  prof.DetProfile().X(),
                  prof.DetProfile().Y(),
                  prof.DetProfile().Z(),
                  prof.DetProfile().T(),
                  prof.DetProfile().Px(),
                  prof.DetProfile().Py(),
                  prof.DetProfile().Pz(),
                  prof.DetProfile().E())
          << std::endl
          << "    Charge per plane: ";
        size_t const nPlanes = prof.Charge().size();
        for (size_t i = 0; i < nPlanes; ++i) {

          std::cout << " | Plane " << i << std::flush;
          std::cout << " ... Q = " << prof.Charge(i) << std::flush;
        }
        std::cout << std::endl << std::endl;
      }
    }
    return result;
  }
}