FindMatches.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    FindMatches.cxx
/// \brief   Core of the LEM match-finding algorithm
/// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "LEM/func/FindMatches.h"

#include "LEM/func/DistanceMap.h"
#include "LEM/func/Heads.h"
#include "LEM/func/MatchableEvent.h"

#include <algorithm>
#include <cfloat>
#include <iostream>
#include <set>

namespace lem
{
  //......................................................................
  float CalcEnergy(const Potential& Va, const MatchableEvent& a, const MatchableEvent& b)
  {
    float energy = a.selfEnergy + b.selfEnergy;

    const LiteHit* hb = b.hits;
    const unsigned int Nb = b.nhits;

    for(unsigned int i = 0; i < Nb; ++i){
      energy -= hb->pecorr*Va.V[hb->cellIdx];
      ++hb;
    }

    return energy;
  }

  //......................................................................
  double FracChargeMatched(const EventSummary& a, const EventSummary& b,
                           bool flipEven, bool flipOdd)
  {
    // TODO: Duplicate hits definitely screw this up

    float qmatched = 0;

    // Quicker to spend this (not much) memory than the more obvious
    // nested-loops implementation.
    float acharge[256*256] = {0,};

    const unsigned int Na = a.nhits;
    const unsigned int Nb = b.nhits;
    for(unsigned int ia = 0; ia < Na; ++ia){
      const LiteHit& ha = a.hits[ia];
      acharge[ha.cellIdx] += ha.pecorr;
    } // end for ia

    for(unsigned int ib = 0; ib < Nb; ++ib){
      const LiteHit& hb = b.hits[ib];

      int hbc = hb.Cell();

      const int v = abs(hb.Plane())%2; // Matches definition in FillPotential
      if((v == 0 && flipEven) ||
         (v == 1 && flipOdd)) hbc = 2*lem::kVertexCell-hbc;

      const int bidx = 256*hb.Plane()+hbc;

      // This is the MINOS definition
      //      if(acharge[bidx]) qmatched += acharge[bidx]+hb.pecorr;

      qmatched += 2*std::min(acharge[bidx], hb.pecorr);
    } // end for ib

    return qmatched/(a.totalPE+b.totalPE);
  }

  //......................................................................
  void FillPotential(const EventSummary& trial, Potential& V,
                     bool flipEven, bool flipOdd)
  {
    const DistanceMap& dm = DistanceMap::Instance();

    const unsigned int N = trial.nhits;

    for(unsigned int i = 0; i < N; ++i){
      const LiteHit& h = trial.hits[i];

      // Doesn't matter? Has to match in qfrac definition
      const int view = abs(h.Plane())%2;

      int hc = h.Cell();
      if((flipEven && view == 0) || (flipOdd && view == 1))
        hc = 2*lem::kVertexCell-hc;

      for(int plane = 0; plane < 256; ++plane){
        // Only in the same view
        if(abs(plane-h.Plane())%2) continue;
        for(int cell = 0; cell < 256; ++cell){

          const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
                                                     cell, hc);

          V.V[256*plane+cell] += contrib;
        }
      }
    }
  }

  //......................................................................
  void FillPotential(const EventSummary& trial, FlippedPotentials& Vs)
  {
    // We could just call FillPotential 4 times with each flip variant like so
    /*
    for(bool flipEven: {false, true}){
      for(bool flipOdd: {false, true}){
        FillPotential(trial, Vs.V[flipEven][flipOdd], flipEven, flipOdd);
      }
    }
    */

    // But this function actually takes noticeable runtime, so instead just
    // fill it once, and then copy into the over three flips.
    const DistanceMap& dm = DistanceMap::Instance();

    const unsigned int N = trial.nhits;

    // There's a complication because 256 is odd there isn't a perfect
    // centreline to flip round. The way things are defined, we end up needing
    // to know what cell 256 in the original map would have been, which we
    // don't have. Make a place to store it here.
    float overflow[256] = {0,};

    for(unsigned int i = 0; i < N; ++i){
      const LiteHit& h = trial.hits[i];

      for(int plane = 0; plane < 256; ++plane){
        // Only in the same view
        if(abs(plane-h.Plane())%2) continue;
        for(int cell = 0; cell < 256; ++cell){

          const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
                                                     cell,  h.Cell());

          Vs.V[false][false].V[256*plane+cell] += contrib;
        }
        // And one for cell 256, stored seperately
        overflow[plane] += h.pecorr*dm.InvDist(plane, h.Plane(),
                                               256,   h.Cell());
      }
    }

    // Now copy over into the flipped variants

    for(bool flipEven: {false, true}){
      for(bool flipOdd: {false, true}){
        if(!flipEven && !flipOdd) continue;

        for(int plane = 0; plane < 256; ++plane){
          const int view = plane%2;

          for(int cell = 0; cell < 256; ++cell){

            int fc = cell; // flipped cell
            if((flipEven && view == 0) || (flipOdd && view == 1))
              fc = 2*lem::kVertexCell-fc;

            if(fc >= 0 && fc < 256)
              Vs.V[flipEven][flipOdd].V[plane*256+cell] = Vs.V[false][false].V[plane*256+fc];
            if(fc == 256) Vs.V[flipEven][flipOdd].V[plane*256+cell] = overflow[plane];
          }
        }
      }
    }
  }

  //......................................................................
  std::vector<Match> FindMatches(const MatchableEvent& trial,
                                 unsigned int libSize,
                                 const MatchableEvent* lib,
                                 const MatchableEvent* libDownsample,
                                 int factor,
                                 const std::multiset<float>& alreadyInput,
                                 unsigned int numMatches,
                                 int enrich,
                                 const Heads* heads,
                                 int headIdx,
                                 bool flipHeadEven,
                                 bool flipHeadOdd,
                                 bool useDownsample)
  {
    //    assert(libsize >= kMaxNumMatches);

    assert(!heads || headIdx >= 0);

    float target = FLT_MAX; // infinity

    // "De-normalize" values back into the range most of this function deals
    // with. Need to copy to a new container, since set iterators are const
    std::multiset<float> already;
    for(auto a: alreadyInput) already.insert(a*trial.selfEnergy);

    // Truncate to the right length. Update target if possible.
    while(already.size() > numMatches){
      auto itLargest = already.end();
      --itLargest;
      already.erase(itLargest);
      target = *already.rbegin();
    }

    // Precalculate potential
    // Try it with either/both of x&y flipped. Effectively gain factor of 4 library stats.
    // TODO: This isn't entirely right because of attenuation. It /would/
    // actually be possible to introduce this properly during library building
    // when truth information is available.

    FlippedPotentials Vs;
    FillPotential(trial, Vs);

    // Make downsampled versions of the four potential cases
    FlippedPotentials Vs_down;
    if(useDownsample){
      for(int flipEven = 0; flipEven <= 1; ++flipEven){
        for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
          Vs_down.V[flipEven][flipOdd] = Vs.V[flipEven][flipOdd].Downsampled(factor);
        }
      }
    }

    // Candidates for return, remains sorted
    std::multiset<Match> Elist;

    // This is the quickest ordering of the loops (locality into lib?)
    for(unsigned int j = 0; j < libSize; ++j){
      int i = j;
      if(heads){
        if(int(j) >= heads->HeadSeqsLen()) break;
        i = heads->ChildIdxFor(headIdx, j);
      }

      if(enrich != 2){
        if(lib[i].enrich != bool(enrich)) continue;
      }

      const MatchableEvent& downEvt = libDownsample[i];

      // There's beam nue in the library, but we don't want to match to it
      if( useDownsample && downEvt.weight == 0) continue;
      if(!useDownsample && lib[i].weight  == 0) continue;

      for(int flipEven = 0; flipEven <= 1; ++flipEven){
        if(heads && (heads->FlipEvenFor(headIdx, j) ^ flipHeadEven) != flipEven) continue;
        for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
          if(heads && (heads->FlipOddFor(headIdx, j) ^ flipHeadOdd) != flipOdd) continue;

          // Compare the downsampled event to the downsampled potential. The
          // resulting energy is a lower-bound on what the real match would
          // have done.
          float E = useDownsample ? CalcEnergy(Vs_down.V[flipEven][flipOdd], trial, downEvt) : 0;

          if(E < target){
            // Calculate the energy fully
            E = CalcEnergy(Vs.V[flipEven][flipOdd], trial, lib[i]);

            // If it's still good, use it
            if(E < target){
              Elist.insert(Match(E, &lib[i], flipEven, flipOdd, j));
              already.insert(E);

              if(Elist.size() > numMatches){
                // Always drop the last one, to keep list correct size
                auto itLargest = Elist.end();
                --itLargest;
                Elist.erase(itLargest);
                // In case "already" isn't long enough
                target = Elist.rbegin()->energy;
              }
              if(already.size() > numMatches){
                // Always drop the last one, to keep list correct size
                auto itLargest = already.end();
                --itLargest;
                already.erase(itLargest);
                target = *already.rbegin();
              }
            }
          }
        } // end for flipOdd
      } // end for flipEven
    } // end for i

    // Convert back to a vector, and calculate the qfrac variable for surviving
    // matches.
    std::vector<Match> ret;
    ret.reserve(numMatches);
    for(const auto& m: Elist){
      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;
  }

} // namespace