FilterGenInTime_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/EDFilter.h"
#include "art/Framework/Core/ModuleMacros.h"

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"

// LArSoft Includes
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/Simulation/sim.h"
#include "nug4/ParticleNavigation/ParticleList.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// C++ Includes
#include <cstring>
#include <sys/stat.h>

// Root includes
#include <TMath.h>

namespace simfilter {

  class FilterGenInTime : public art::EDFilter {
  public:
    explicit FilterGenInTime(fhicl::ParameterSet const& pset);

    bool filter(art::Event&);

    virtual void beginJob();

  private:
    bool KeepParticle(simb::MCParticle const& part) const;

    std::vector<std::array<double, 6>> fCryostatBoundaries; 
    double fMinKE;       //<only keep based on particles with greater than this energy
    bool fKeepOnlyMuons; //keep based only on muons if enabled
    double fMinT, fMaxT; //<time range in which to keep particles
    bool fSortParticles; //create new MCTruth collections with particles sorted by their timing
    bool fAlwaysPass;    // flag to have filter always pass (to be used with sorting...)
  };

} // namespace simfilter

namespace simfilter {

  FilterGenInTime::FilterGenInTime(fhicl::ParameterSet const& pset)
    : EDFilter{pset}
    , fMinKE(pset.get<double>("MinEnergy", 0.0))
    , fKeepOnlyMuons(pset.get<bool>("KeepOnlyMuons", false))
    , fMinT(pset.get<double>("MinT", 0.0))
    , fMaxT(pset.get<double>("MaxT"))
    , fSortParticles(pset.get<bool>("SortParticles", false))
    , fAlwaysPass(pset.get<bool>("AlwaysPass", false))
  {
    if (fSortParticles) {
      produces<std::vector<simb::MCTruth>>("intime");
      produces<std::vector<simb::MCTruth>>("outtime");
    }
    std::cout << "New Filter!\n";
  }

  void FilterGenInTime::beginJob()
  {
    auto const& geom = *art::ServiceHandle<geo::Geometry const>();
    for (auto const& cryo : geom.Iterate<geo::CryostatGeo>()) {
      std::array<double, 6> this_cryo_boundaries{};
      cryo.Boundaries(&this_cryo_boundaries[0]);
      fCryostatBoundaries.push_back(this_cryo_boundaries);
    }
  }

  bool FilterGenInTime::KeepParticle(simb::MCParticle const& part) const
  {
    const TLorentzVector& v4 = part.Position();
    const TLorentzVector& p4 = part.Momentum();
    // origin of particle
    double x0[3] = {v4.X(), v4.Y(), v4.Z()};
    // normalized direction of particle
    double dx[3] = {
      p4.Px() / p4.Vect().Mag(), p4.Py() / p4.Vect().Mag(), p4.Pz() / p4.Vect().Mag()};

    // tolernace for treating number as "zero"
    double eps = 1e-5;

    //check to see if particle crosses boundary of any cryostat within appropriate time window
    std::vector<bool> intersects_cryo(fCryostatBoundaries.size(), false);
    std::vector<bool> inside_cryo(fCryostatBoundaries.size(), false);
    std::vector<double> distance_to_cryo(fCryostatBoundaries.size(), 0.);

    // Check to see if particle intersects any cryostat
    //
    // Algorithmically, this is looking for ray-box intersection. This is a common problem in
    // computer graphics. The algorithm below is taken from "Graphics Gems", Academic Press, 1990
    for (size_t i_cryo = 0; i_cryo < fCryostatBoundaries.size(); i_cryo++) {
      auto const& bound = fCryostatBoundaries[i_cryo];
      std::array<int, 3> quadrant{}; // 0 == RIGHT, 1 == LEFT, 2 == MIDDLE
      std::array<double, 3> candidatePlane{};
      std::array<double, 3> coord{};

      std::array<double, 3> bound_lo = {{bound[0], bound[2], bound[4]}};
      std::array<double, 3> bound_hi = {{bound[1], bound[3], bound[5]}};

      // First check if origin is inside box
      // Also check which of the two planes in each dimmension is the
      // "candidate" for the ray to hit
      bool inside = true;
      for (int i = 0; i < 3; i++) {
        if (x0[i] < bound_lo[i]) {
          quadrant[i] = 1; // LEFT
          candidatePlane[i] = bound_lo[i];
          inside = false;
        }
        else if (x0[i] > bound_hi[i]) {
          quadrant[i] = 0; // RIGHT
          candidatePlane[i] = bound_hi[i];
          inside = false;
        }
        else {
          quadrant[i] = 2; // MIDDLE
        }
      }

      if (inside) {
        inside_cryo[i_cryo] = true;
        // if we're inside the cryostat, then we do intersect it
        intersects_cryo[i_cryo] = true;
        continue;
      }

      // ray origin is outside the box -- calculate the distance to the cryostat and see if it intersects

      // calculate distances to candidate planes
      std::array<double, 3> maxT{};
      for (int i = 0; i < 3; i++) {
        if (quadrant[i] != 2 /* MIDDLE */ && abs(dx[i]) > eps) {
          maxT[i] = (candidatePlane[i] - x0[i]) / dx[i];
        }
        // if a ray origin is between two the two planes in a dimmension, it would never hit that plane first
        else {
          maxT[i] = -1;
        }
      }

      // The plane on the box that the ray hits is the one with the largest distance
      int whichPlane = 0;
      for (int i = 1; i < 3; i++) {
        if (maxT[whichPlane] < maxT[i]) whichPlane = i;
      }

      // check if the candidate intersection point is inside the box

      // no intersection
      if (maxT[whichPlane] < 0.) {
        intersects_cryo[i_cryo] = false;
        continue;
      }

      for (int i = 0; i < 3; i++) {
        if (whichPlane != i) { coord[i] = x0[i] + maxT[whichPlane] * dx[i]; }
        else {
          coord[i] = candidatePlane[i];
        }
      }

      // check if intersection is in box
      intersects_cryo[i_cryo] = true;
      for (int i = 0; i < 3; i++) {
        if (coord[i] < bound_lo[i] || coord[i] > bound_hi[i]) { intersects_cryo[i_cryo] = false; }
      }

      if (intersects_cryo[i_cryo]) { distance_to_cryo[i_cryo] = maxT[whichPlane]; }
    }

    // check if any cryostats are intersected in-time
    for (size_t i_cryo = 0; i_cryo < fCryostatBoundaries.size(); i_cryo++) {

      // If the particle originates inside the cryostat, then
      // we can't really say when it will leave. Thus, accept
      // the particle
      if (inside_cryo[i_cryo]) { return true; }
      // otherwise check arrival time at boundary of cryostat
      if (intersects_cryo[i_cryo]) {
        double ptime = (distance_to_cryo[i_cryo] * 1e-2 /* cm -> m */) /
                       (TMath::C() * sqrt(1 - pow(part.Mass() / part.E(), 2))) /* velocity */;
        double totT = part.T() + ptime * 1e9 /* s -> ns */;
        if (totT > fMinT && totT < fMaxT) { return true; }
      }
    }

    return false;
  }

  bool FilterGenInTime::filter(art::Event& evt)
  {

    std::unique_ptr<std::vector<simb::MCTruth>> truthInTimePtr(new std::vector<simb::MCTruth>(1));
    std::unique_ptr<std::vector<simb::MCTruth>> truthOutOfTimePtr(
      new std::vector<simb::MCTruth>(1));

    simb::MCTruth& truthInTime = truthInTimePtr->at(0);
    simb::MCTruth& truthOutOfTime = truthOutOfTimePtr->at(0);

    //get the list of particles from this event
    art::ServiceHandle<geo::Geometry const> geom;

    //std::vector< art::Handle< std::vector<simb::MCTruth> > > allmclists;
    //evt.getManyByType(allmclists);
    auto allmclists = evt.getMany<std::vector<simb::MCTruth>>();
    bool keepEvent = false;
    for (size_t mcl = 0; mcl < allmclists.size(); ++mcl) {
      art::Handle<std::vector<simb::MCTruth>> mclistHandle = allmclists[mcl];
      for (size_t m = 0; m < mclistHandle->size(); ++m) {
        art::Ptr<simb::MCTruth> mct(mclistHandle, m);

        for (int ipart = 0; ipart < mct->NParticles(); ipart++) {
          bool kp = KeepParticle(mct->GetParticle(ipart));

          if (kp && (!fKeepOnlyMuons || abs(mct->GetParticle(ipart).PdgCode()) == 13) &&
              mct->GetParticle(ipart).E() - mct->GetParticle(ipart).Mass() > fMinKE) {
            keepEvent = true;
            if (!fSortParticles) break;
          }

          if (fSortParticles) {
            simb::MCParticle particle = mct->GetParticle(ipart);
            if (kp) truthInTime.Add(particle);
            if (!kp) truthOutOfTime.Add(particle);
          }

        } //end loop over particles

        if (!fSortParticles && keepEvent) break;

      } //end loop over mctruth col

      if (!fSortParticles && keepEvent) break;

    } //end loop over all mctruth lists

    if (fSortParticles) {
      evt.put(std::move(truthInTimePtr), "intime");
      evt.put(std::move(truthOutOfTimePtr), "outtime");
    }

    return (keepEvent || fAlwaysPass);
  }

} // namespace simfilter

namespace simfilter {

  DEFINE_ART_MODULE(FilterGenInTime)

} // namespace simfilter