FindMatchesAlg.cxx

Go to the documentation of this file.

#include "LEM/func/FindMatchesAlg.h"

#include "LEM/func/FindMatches.h" // For FillPotentials

#include <cassert>
#include <cfloat>

namespace lem
{
  // TODO investigate other data structures
  class SortedMatchList: public std::multiset<Match>
  {
  public:
    SortedMatchList(unsigned int numMatches)
      : fNumMatches(numMatches), fTarget(FLT_MAX) {}

    template<class... T> void emplace(double E, T... args)
    {
      if(E >= fTarget) return;

      std::multiset<Match>::emplace(E, args...);

      if(size() > fNumMatches){
        // Always drop the last one, to keep list correct size
        auto itLargest = end();
        --itLargest;
        erase(itLargest);
        fTarget = rbegin()->energy;
      }
    }

  protected:
    void insert(Match m); // not implemented

    unsigned int fNumMatches;
    double fTarget;
  };

  //......................................................................
  FindMatchesAlg::FindMatchesAlg(const std::string& libDir,
                                 bool preload,
                                 bool useHeads)
  {
    fLib = Library::FromMMap(libDir, preload);
    assert(fLib);

    if(useHeads){
      fHeads = Heads::FromMMap(libDir, preload);
      assert(fHeads);
    }
    else{
      fHeads = 0;
    }
  }

  //......................................................................
  FindMatchesAlg::~FindMatchesAlg()
  {
    delete fLib;
    delete fHeads;
  }

  //......................................................................
  std::vector<Match> FindMatchesAlg::
  PackageMatches(const std::multiset<Match>& matches,
                 const MatchableEvent& trial) const
  {
    // Convert back to a vector, and calculate the qfrac variable for surviving
    // matches.
    std::vector<Match> ret;
    ret.reserve(matches.size());
    for(const Match& m: matches){
      ret.push_back(m);
      ret.back().qfrac = FracChargeMatched(trial, *m.evt, m.flipEven, m.flipOdd);
      ret.back().energy /= trial.selfEnergy; // Normalize to approx 0-1
    }

    return ret;
  }

  //......................................................................
  std::vector<Match> FindMatchesAlg::
  FindMatches(const MatchableEvent& trial,
              unsigned int numMatches,
              int enrich) const
  {
    if(fHeads)
      return FindMatchesHeads  (trial, numMatches, enrich);
    else
      return FindMatchesNoHeads(trial, numMatches, enrich);
  }

  //......................................................................
  std::vector<Match> FindMatchesAlg::
  FindMatchesNoHeads(const MatchableEvent& trial,
                     unsigned int numMatches,
                     int enrich) const
  {
    FlippedPotentials trialVs;
    FillPotential(trial, trialVs);

    // Candidates for return, remains sorted
    SortedMatchList Elist(numMatches);

    // This is the quickest ordering of the loops (locality into lib?)
    const unsigned int I = fLib->NEvents();
    for(unsigned int i = 0; i < I; ++i){
      const MatchableEvent& libEvt = fLib->Event(i);

      if(enrich != 2){
        if(libEvt.enrich != bool(enrich)) continue;
      }

      // There's beam nue in the library, but we don't want to match to it
      if(libEvt.weight  == 0) continue;

      for(int flipEven: {0, 1}){
        for(int flipOdd: {0, 1}){

          const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial, libEvt);
          Elist.emplace(E, &libEvt, flipEven, flipOdd, i);
        } // end for flipOdd
      } // end for flipEven
    } // end for i

    return PackageMatches(Elist, trial);
  }

  //......................................................................
  int FindMatchesAlg::BestHeadFor(const FlippedPotentials& trialVs,
                                  const MatchableEvent& trial,
                                  bool& bestFlipEven, bool& bestFlipOdd,
                                  double& bestEnergy) const
  {
    int bestHead = -1;
    bestEnergy = FLT_MAX;

    const int nHeads = fHeads->NHeads();
    for(int iHead = 0; iHead < nHeads; ++iHead){
      const MatchableEvent& headEvt = fLib->Event(fHeads->HeadIdx(iHead));

      for(int flipEven: {0, 1}){
        for(int flipOdd: {0, 1}){
          const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial,
                                     headEvt);
          if(E < bestEnergy){
            bestEnergy = E;
            bestHead = iHead;
            bestFlipEven = flipEven;
            bestFlipOdd = flipOdd;
          }
        }
      }
    }
    assert(bestHead >= 0);

    return bestHead;
  }

  //......................................................................
  std::vector<Match> FindMatchesAlg::
  FindMatchesHeads(const MatchableEvent& trial,
                   unsigned int numMatches,
                   int enrich) const
  {
    FlippedPotentials trialVs;
    FillPotential(trial, trialVs);

    bool headFlipEven, headFlipOdd;
    double headE;
    const int headIdx = BestHeadFor(trialVs, trial,
                                    headFlipEven, headFlipOdd, headE);
    //    const MatchableEvent& head = fLib->Event(fHeads->HeadIdx(headIdx));

    // Candidates for return, remains sorted
    SortedMatchList Elist(numMatches);

    // This is the quickest ordering of the loops (locality into lib?)
    const unsigned int J = fHeads->HeadSeqsLen();
    for(unsigned int j = 0; j < J; ++j){
      const int i = fHeads->ChildIdxFor(headIdx, j);

      const MatchableEvent& libEvt = fLib->Event(i);

      if(enrich != 2){
        if(libEvt.enrich != bool(enrich)) continue;
      }

      // There's beam nue in the library, but we don't want to match to it
      if(libEvt.weight  == 0) continue;

      const bool flipEven = (fHeads->FlipEvenFor(headIdx, j) ^ headFlipEven);
      const bool flipOdd  = (fHeads->FlipOddFor (headIdx, j) ^ headFlipOdd );

      const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial, libEvt);
      Elist.emplace(E, &libEvt, flipEven, flipOdd, j);
    } // end for i/j

    return PackageMatches(Elist, trial);
  }
} // namespace