FindLEMMatches_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    FindMatches_module.cc
// \brief   Module to find best LEM matches
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "LEM/func/DistanceMap.h"
#include "LEM/func/FindMatchesAlg.h"
#include "LEM/func/Library.h"
#include "LEM/func/MatchIndices.h"
#include "LEM/func/MatchList.h"

#include "LEM/LEMInput.h"
#include "LEM/Util.h"

// ROOT includes
#include "TH2.h"

// NOvASoft includes
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/Vertex.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/func/MathUtil.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Principal/Handle.h"
#include "fhiclcpp/ParameterSet.h"

/// \f$ \nu_e \f$ %PID
namespace lem
{
  /// Module to find best %LEM matches
  class FindLEMMatches: public art::EDProducer // TODO rename once FindMatches is gone
  {
  public:
    explicit FindLEMMatches(const fhicl::ParameterSet& pset);
    ~FindLEMMatches();

    virtual void reconfigure(const fhicl::ParameterSet& pset);
    virtual void beginRun(art::Run& run);
    virtual void produce(art::Event& evt);
  protected:
    void AddHists(std::vector<TH2F>* th2col, const Match& bestMatch,
                  const EventSummary& trial,
                  const EventSummary* head,
                  bool headFlipEven, bool headFlipOdd) const;

    rb::Vertex MatchToVertex(const Match& match,
                             int vtxPlane, int vtxCell, int vtxCellOther) const;

    FindMatchesAlg* fAlg;

    std::string fInputLabel;

    bool fSaveEventImages, fSaveEventPotentials;
  };

  //......................................................................
  FindLEMMatches::FindLEMMatches(const fhicl::ParameterSet& pset)
  {
    reconfigure(pset);

    produces<std::vector<lem::MatchList> >();
    // MatchLists are associated with the slice they are against
    produces<art::Assns<lem::MatchList, rb::Cluster> >();

    produces<std::vector<lem::MatchIndices> >();
    // MatchLists are associated with the slice they are against
    produces<art::Assns<lem::MatchIndices, rb::Cluster> >();

    // True vertex of best match. Associate to slice
    produces<std::vector<rb::Vertex> >();
    produces<art::Assns<rb::Vertex, rb::Cluster> >();

    // Really should be file level, not run level, but don't know how to do
    // that
    produces<lem::LibrarySummary, art::InRun>();

    // Just to mark where the vertex is so you can check it's right
    produces<std::vector<geo::OfflineChan> >();

    // Images of best match
    if(fSaveEventImages || fSaveEventPotentials)
      produces<std::vector<TH2F> >();
  }

  //......................................................................
  FindLEMMatches::~FindLEMMatches()
  {
    delete fAlg;
  }

  //......................................................................
  void FindLEMMatches::reconfigure(const fhicl::ParameterSet& pset)
  {
    fAlg = new FindMatchesAlg(util::EnvExpansion(pset.get<std::string>("LibraryPath")),
                              pset.get<bool>("PreloadLib"),
                              pset.get<bool>("UseHeads"));

    fInputLabel = pset.get<std::string>("InputLabel");

    // This needs to be before anyone asks the DistanceMap anything, should be
    // safe enough here.
    DistanceMap::SetPlaneScale(pset.get<double>("DistPlaneScale"));
    DistanceMap::SetCellScale(pset.get<double>("DistCellScale"));
    DistanceMap::SetDecayPower(pset.get<double>("DistDecayPower"));
    DistanceMap::SetExpMode(pset.get<bool>("DistExpMode"));

    EventSummary::fgDecayRate = pset.get<double>("HitDecayRate");
    EventSummary::fgChargePower = pset.get<double>("HitChargePower");

    fSaveEventImages = pset.get<bool>("SaveEventImages", true);
    fSaveEventPotentials = pset.get<bool>("SaveEventPotentials", false);
  }

  //......................................................................
  void FindLEMMatches::beginRun(art::Run& run)
  {
    std::unique_ptr<lem::LibrarySummary> ls(new lem::LibrarySummary);
    *ls = fAlg->GetLibrary()->Summary();
    run.put(std::move(ls));
  }

  //......................................................................
  void FindLEMMatches::produce(art::Event& evt)
  {
    std::unique_ptr<std::vector<lem::MatchList> > matchcol(new std::vector<lem::MatchList>);
    std::unique_ptr<std::vector<lem::MatchIndices>> matchidxs(new std::vector<lem::MatchIndices>);

    // Associate matches with the slice they came from
    std::unique_ptr<art::Assns<lem::MatchList, rb::Cluster> > assns(new art::Assns<lem::MatchList, rb::Cluster>);
    std::unique_ptr<art::Assns<lem::MatchIndices, rb::Cluster> > idxassns(new art::Assns<lem::MatchIndices, rb::Cluster>);

    std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
    std::unique_ptr<art::Assns<rb::Vertex, rb::Cluster> > vtxassns(new art::Assns<rb::Vertex, rb::Cluster>);

    std::unique_ptr<std::vector<geo::OfflineChan> > chancol(new std::vector<geo::OfflineChan>);

    std::unique_ptr<std::vector<TH2F>> th2col(new std::vector<TH2F>);


    art::Handle<std::vector<LEMInput>> inputs;
    evt.getByLabel(fInputLabel, inputs);
    art::FindManyP<rb::Cluster> fmt(inputs, evt, fInputLabel);

    const int inputMax = inputs->size();
    for(int inputIdx = 0; inputIdx < inputMax; ++inputIdx){
      const LEMInput& input = (*inputs)[inputIdx];

      EventSummary sum(input.hits, -1, std::vector<int>());

      sum.run = evt.run();
      sum.subrun = evt.subRun();
      sum.event = evt.event();

      sum.origPdg = 0;

      sum.pdg = 0;
      sum.ccnc = 0;
      sum.y = 0;
      sum.trueEVis = 0;


      //      chancol->push_back(geo::OfflineChan(vtxPlane, vtxCell));
      // Should really be nextPlane, but that's lost inside of Util.cxx. This
      // is almost always right, and only for display purposes anyway.
      //      chancol->push_back(geo::OfflineChan(vtxPlane+1, vtxCellOther));

      MatchableEvent trial(sum);

      const std::vector<Match> matches = fAlg->FindMatches(trial, kMaxNumMatches, false);
      const std::vector<Match> matchesEnrich = fAlg->FindMatches(trial, kMaxNumMatches, true);

/*
      if(!matches.empty()){
        vtxcol->push_back(MatchToVertex(matches[0], vtxPlane, vtxCell, vtxCellOther));
        vtxcol->back().SetT(slice.MinTNS());
        util::CreateAssn(*this, evt, *vtxcol,
                         art::Ptr<rb::Cluster>(slices, sliceIdx), *vtxassns);
      }
*/

      // TODO TODO TODO
/*
      if(!matches.empty()){
        AddHists(th2col.get(), matches[0], sum,
                 (bestHead >= 0) ? &fLib->Event(fHeads->HeadIdx(bestHead)) : 0,
                 headFlipEven, headFlipOdd);
      }
*/

      matchcol->push_back(MatchList());
      MatchList& ml = matchcol->back();

      matchidxs->push_back(MatchIndices());
      MatchIndices& mi = matchidxs->back();
      // TODO TODO TODO
      //      mi.headIdx = (bestHead >= 0) ? fHeads->HeadIdx(bestHead) : -1;

      for(unsigned int i = 0; i < matches.size(); ++i){
        ml.matches.push_back(MatchSummary(matches[i]));
        mi.matches.push_back(matches[i].evt - &fAlg->GetLibrary()->Event(0));
      } // end for i

      for(unsigned int i = 0; i < matchesEnrich.size(); ++i){
        ml.matches.push_back(MatchSummary(matchesEnrich[i]));
        mi.matchesEnrich.push_back(matchesEnrich[i].evt - &fAlg->GetLibrary()->Event(0));
      } // end for i

      ml.trial = MatchSummary(Match(0, &trial, false, false));

      // TODO TODO TODO
      //      ml.headIdx = (bestHead >= 0) ? fHeads->HeadIdx(bestHead) : -1;

      // If the original input hits were associated to a cluster, then so
      // should our matches be.
      try{
        std::vector<art::Ptr<rb::Cluster>> slices = fmt.at(inputIdx);
        if(slices.size() == 1){
          util::CreateAssn(*this, evt, *matchcol, slices[0], *assns);
          util::CreateAssn(*this, evt, *matchidxs, slices[0], *idxassns);
        }
      }
      catch(...){
        // ugh
      }
    } // end for sliceIdx

    evt.put(std::move(matchcol));
    evt.put(std::move(assns));

    evt.put(std::move(vtxcol));
    evt.put(std::move(vtxassns));

    evt.put(std::move(chancol));
    if(fSaveEventImages || fSaveEventPotentials) evt.put(std::move(th2col));
  }

  //......................................................................
  void FindLEMMatches::AddHists(std::vector<TH2F>* th2col,
                                const Match& bestMatch,
                                const EventSummary& trial,
                                const EventSummary* head,
                                bool headFlipEven, bool headFlipOdd) const
  {
    // Best
    const EventSummary* bestEvt = bestMatch.evt;

    const std::string bestTitle = bestEvt->Description()+";Plane;Cell";

    if(fSaveEventImages){
      TH2F evenBest("bestEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192);
      TH2F oddBest ("bestOdd",  bestTitle.c_str(), 64, 0, 128, 128, 64, 192);

      bestEvt->FillHists(&evenBest, &oddBest, bestMatch.flipEven, bestMatch.flipOdd);

      th2col->push_back(evenBest);
      th2col->push_back(oddBest);
    }

    if(fSaveEventPotentials){
      TH2F hVsBest[2] = {TH2F("bestVEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192),
                         TH2F("bestVOdd",  bestTitle.c_str(), 64, 0, 128, 128, 64, 192)};
      Potential VBest;
      FillPotential(*bestEvt, VBest, bestMatch.flipEven, bestMatch.flipOdd);
      for(int plane = 0; plane < 256; ++plane){
        for(int cell = 0; cell < 256; ++cell){
          const double p = VBest.V[256*plane+cell];
          hVsBest[plane%2].Fill(plane, cell, p);
        }
      }
      th2col->push_back(hVsBest[0]);
      th2col->push_back(hVsBest[1]);
    }


    // Head
    if(head){
      const char* headTitle = (head->Description()+";Plane;Cell").c_str();

      if(fSaveEventImages){
        TH2F evenHead("headEven", headTitle, 64, 0, 128, 128, 64, 192);
        TH2F oddHead ("headOdd",  headTitle, 64, 0, 128, 128, 64, 192);

        head->FillHists(&evenHead, &oddHead, headFlipEven, headFlipOdd);

        th2col->push_back(evenHead);
        th2col->push_back(oddHead);
      }

      if(fSaveEventPotentials){
        TH2F hVsHead[2] = {TH2F("headVEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192),
                           TH2F("headVOdd",  bestTitle.c_str(), 64, 0, 128, 128, 64, 192)};
        Potential VHead;
        FillPotential(*head, VHead, headFlipEven, headFlipOdd);
        for(int plane = 0; plane < 256; ++plane){
          for(int cell = 0; cell < 256; ++cell){
            const double p = VHead.V[256*plane+cell];
            hVsHead[plane%2].Fill(plane, cell, p);
          }
        }
        th2col->push_back(hVsHead[0]);
        th2col->push_back(hVsHead[1]);
      }
    }


    // Trial
    const char* trialTitle = "Trial event;Plane;Cell";

    if(fSaveEventImages){
      TH2F evenTrial("trialEven", trialTitle, 64, 0, 128, 128, 64, 192);
      TH2F oddTrial ("trialOdd",  trialTitle, 64, 0, 128, 128, 64, 192);

      trial.FillHists(&evenTrial, &oddTrial);

      th2col->push_back(evenTrial);
      th2col->push_back(oddTrial);
    }

    if(fSaveEventPotentials){
      TH2F hVsTrial[2] = {TH2F("trialVEven", trialTitle, 64, 0, 128, 128, 64, 192),
                          TH2F("trialVOdd",  trialTitle, 64, 0, 128, 128, 64, 192)};
      Potential VTrial;
      FillPotential(trial, VTrial, false, false);
      for(int plane = 0; plane < 256; ++plane){
        for(int cell = 0; cell < 256; ++cell){
          const double p = VTrial.V[256*plane+cell];
          hVsTrial[plane%2].Fill(plane, cell, p);
        }
      }
      th2col->push_back(hVsTrial[0]);
      th2col->push_back(hVsTrial[1]);
    }
  }

  //......................................................................
  rb::Vertex FindLEMMatches::MatchToVertex(const Match& match,
                                           int vtxPlane, int vtxCell,
                                           int vtxCellOther) const
  {
    art::ServiceHandle<geo::Geometry> geom;

    int plane = vtxPlane+match.evt->trueVtxPlane-kVertexPlane;

    if(plane < 0) plane = 0;
    if(plane > int(geom->NPlanes())) plane = geom->NPlanes()-1;

    int cell, cellOther;
    if(abs(match.evt->trueVtxPlane-kVertexPlane)%2 == 0){
      cell = vtxCell+match.evt->trueVtxCell-kVertexCell;
      cellOther = vtxCellOther+match.evt->trueVtxCellOther-kVertexCell;
    }
    else{
      cell = vtxCellOther+match.evt->trueVtxCell-kVertexCell;
      cellOther = vtxCell+match.evt->trueVtxCellOther-kVertexCell;
    }

    int planeOther = 0;
    if(geom->Plane(plane)){
      geom->NextPlaneOtherView(plane, +1);
      if(planeOther == geo::kPLANE_NOT_FOUND) planeOther = geom->NextPlaneOtherView(plane, -1);
    }

    if(cell < 0) cell = 0;
    if(cellOther < 0) cellOther = 0;
    if(cell >= int(geom->Plane(plane)->Ncells())) cell = geom->Plane(plane)->Ncells()-1;
    if(cellOther >= int(geom->Plane(planeOther)->Ncells())) cellOther = geom->Plane(planeOther)->Ncells()-1;

    double xyz[3];
    geom->Plane(plane)->Cell(cell)->GetCenter(xyz);

    double xyz2[3];
    geom->Plane(planeOther)->Cell(cellOther)->GetCenter(xyz2);

    const geo::View_t v = geom->Plane(plane)->View();

    return rb::Vertex((v == geo::kX) ? xyz[0] : xyz2[0],
                      (v == geo::kY) ? xyz[1] : xyz2[1],
                      xyz[2],
                      -1); // Don't bother to fill the time
  }

  DEFINE_ART_MODULE(FindLEMMatches)

} //namespace lem
////////////////////////////////////////////////////////////////////////