MakePID_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    MakePID_module.cc
// \brief   Calculate final PID variables from match details
// \author  Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "LEM/PIDDetails.h"

#include "LEM/LEMInput.h"
#include "LEM/func/FindMatches.h" // For num matches constant
#include "LEM/func/LibrarySummary.h"
#include "LEM/func/MatchList.h"

#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "nusimdata/SimulationBase/MCTruth.h" // mode enum
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/EnvExpand.h"
#include "Utilities/func/MathUtil.h"

//#include "doublefann.h" // from the FANN external

#include "TMVA/Reader.h"

#include "LEM/func/DecisionTree.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/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

float gTMVAVar[6];
float gJunk;

#include "TFile.h"
#include "TH2.h"

namespace lem
{
  /// Calculate final %PID variables from match details
  class MakePID: public art::EDProducer
  {
  public:
    explicit MakePID(const fhicl::ParameterSet& pset);
    ~MakePID();

    virtual void beginJob();
    virtual void beginRun(art::Run& run);
    virtual void produce(art::Event& evt);
    virtual void reconfigure(const fhicl::ParameterSet& pset);

  protected:
    bool fScaleCalE;

    enum WeightScheme{
      kUnweighted,
      kSimple,
      kExp,
      kSigOnly,
      kSigOnlyExp,
      kBkgOnly,
      kBkgOnlyExp,
      kSig50Only,
      kSig90Only,
      kPi0Only
    };

    enum MatchVar{
      kSigFrac,
      kBkgFrac,
      kNumuFrac,
      kY,
      kQFrac,
      kSigFrac50,
      kSigFrac90,
      kL0,
      kL1,
      kTrueEVis,
      kQESigFrac,
      kResSigFrac,
      kDISSigFrac,
      kQENCFrac,
      kResNCFrac,
      kDISNCFrac,
      kQENumuFrac,
      kResNumuFrac,
      kDISNumuFrac
    };

    void AvgAndFit(const std::vector<MatchSummary>& matches,
                   WeightScheme weightScheme, MatchVar matchVar,
                   double& mean, double& fit) const;

    std::vector<double> fExpTable;

    double fNCScale;
    double fSigScale;
    double fEnrichScale;

    double fExpLength, fExpPow;

    //    fann *fANN, *fANNEnrich;

    bool fTrimMatchLists;
    bool fUseFitInput;
    bool fPositionalWeight;
    std::string fMatchesLabel;
    std::string fSummarizeLabel;

    TMVA::Reader* fReader;
    dec::Forest* fForest;

    TH1* fFlatHist;
  };

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

    produces<std::vector<lem::PIDDetails>>();
    produces<art::Assns<lem::PIDDetails, rb::Cluster>>();

    const std::string path     = pset.get<std::string>("LibraryPath");
    const std::string libPath  = util::EnvExpansion(path) + "/";

    fReader = 0;
    if(pset.get<bool>("FillTMVAs")){
      fReader = new TMVA::Reader;

      fReader->AddVariable("pidExtrap", &gTMVAVar[0]); // vars.fPIDFit
      fReader->AddVariable("yexp",      &gTMVAVar[1]); // vars.fMeanYExp
      fReader->AddVariable("qfracexp",  &gTMVAVar[2]); // vars.fMeanQFracExp
      fReader->AddVariable("ediff",     &gTMVAVar[3]); // vars.fEnergyDiffExp
      fReader->AddVariable("calE",      &gTMVAVar[4]); // 2*totalRecoGeV
      fReader->AddVariable("enrichExp", &gTMVAVar[5]); // vars.fEnrichFracExp

      fReader->AddSpectator("isSig", &gJunk);
      fReader->AddSpectator("myW",   &gJunk);

      const std::string tmvaPath = libPath+"weights/";

      fReader->BookMVA("BDTG", tmvaPath+"TMVAClassification_BDTG.weights.xml");
      fReader->BookMVA("BDT",  tmvaPath+"TMVAClassification_BDT.weights.xml");
      fReader->BookMVA("BDTD", tmvaPath+"TMVAClassification_BDTD.weights.xml");
      fReader->BookMVA("BDTMitFisher", tmvaPath+"TMVAClassification_BDTMitFisher.weights.xml");
      fReader->BookMVA("MLP",  tmvaPath+"TMVAClassification_MLP.weights.xml");

      fReader->BookMVA("BDTG_CC", tmvaPath+"TMVAClassification_BDTG_CC.weights.xml");
      fReader->BookMVA("MLP_CC",  tmvaPath+"TMVAClassification_MLP_CC.weights.xml");
    }

    // fANN = fANNEnrich = 0;
    // if(pset.get<bool>("FillANNs")){
    //   fANN = fann_create_from_file((libPath+"fann.net").c_str());
    //   assert(fANN);
    //   fann_print_connections(fANN);

    //   fANNEnrich = fann_create_from_file((libPath+"fann_enrich.net").c_str());
    //   assert(fANNEnrich);
    //   fann_print_connections(fANNEnrich);
    // }

    fForest = 0;
    if(pset.get<bool>("FillDecTree")){
      std::cout << "Loading LEM dectree from " << libPath << std::endl;
      fForest = dec::Forest::FromFile(libPath+"dectree.root");
    }

    fFlatHist = 0;
    if(pset.get<bool>("Flatten")){
      TFile* fFlat = new TFile((libPath+"flatten.root").c_str());
      assert(!fFlat->IsZombie());
      fFlatHist = (TH1*)fFlat->Get("conv");
      assert(fFlatHist);
    }
  }

  //......................................................................
  MakePID::~MakePID()
  {
    delete fReader;
    delete fForest;
  }

  //......................................................................
  void MakePID::reconfigure(const fhicl::ParameterSet& pset)
  {
    fScaleCalE = pset.get<bool>("ScaleCalE");

    fTrimMatchLists = pset.get<bool>("TrimMatchLists");
    fUseFitInput = pset.get<bool>("UseFitInput");
    fPositionalWeight = pset.get<bool>("PositionalWeight");
    fMatchesLabel = pset.get<std::string>("MatchesLabel");
    fSummarizeLabel = pset.get<std::string>("SummarizeLabel");
    fExpLength = pset.get<double>("ExpLength");
    fExpPow = pset.get<double>("ExpPow");
  }

  //......................................................................
  void MakePID::beginJob()
  {
    fExpTable.resize(kMaxNumMatches);
    for(unsigned int i = 0; i < kMaxNumMatches; ++i){
      fExpTable[i] = exp(-pow(i/fExpLength, fExpPow));
    }
  }

  //......................................................................
  void MakePID::beginRun(art::Run& run)
  {
    art::Handle<lem::LibrarySummary> libsum;
    run.getByLabel(fMatchesLabel, libsum);

    assert(!libsum.failedToGet());

    fNCScale = libsum->nTrueNC/double(libsum->nTrueNC+libsum->nSwapNC);
    fSigScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totSig;
    fEnrichScale = double(libsum->totBkg-libsum->nSwapNC)/libsum->totEnrich;

    std::cout << "Reweighted to sig, bkg, enrich = "
              << libsum->totSig*fSigScale << ", "
              << libsum->totBkg+(fNCScale-1)*libsum->nTrueNC-libsum->nSwapNC << ", "
              << libsum->totEnrich*fEnrichScale;
  }

  //......................................................................
  void MakePID::AvgAndFit(const std::vector<MatchSummary>& matches,
                          WeightScheme weightScheme, MatchVar matchVar,
                          double& mean, double& fit) const
  {
    std::vector<double> xs, ys, ws;
    double tot = 0;
    mean = 0;

    const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;

    for(unsigned int i = 0; i < matches.size(); ++i){
      const MatchSummary& match = matches[i];

      double w = match.weight;

      if(weightScheme == kUnweighted) w = 1;

      if(weightScheme == kExp ||
         weightScheme == kSigOnlyExp ||
         weightScheme == kBkgOnlyExp){
        if(fPositionalWeight){
          w *= fExpTable[i];
        }
        else{
          w *= exp(-pow(match.potential/maxPotential, fExpPow)/fExpLength);
        }
      }

      if(weightScheme == kSigOnly || weightScheme == kSigOnlyExp){
        if(!match.isSig) w = 0;
      }

      if(weightScheme == kBkgOnly || weightScheme == kBkgOnlyExp){
        if(match.isSig) w = 0;
      }

      if(weightScheme == kSig50Only){
        if(!match.isSig || match.y > 0.50) w = 0;
      }
      if(weightScheme == kSig90Only){
        if(!match.isSig || match.y > 0.90) w = 0;
      }
      if(weightScheme == kPi0Only) w = (match.photL0 > 0 || match.photL1 > 0 || match.photE0 > 0 || match.photE1 > 0);

      double var = 0;
      switch(matchVar){
      case kSigFrac: var = int(match.isSig); break;
      case kBkgFrac: var = 1-int(match.isSig); break;
      case kNumuFrac: var = int(match.ccnc == simb::kCC && abs(match.pdg) == 14); break;
      case kY: var = match.y; break;
      case kQFrac: var = match.qFrac; break;
      case kSigFrac50: var = int(match.isSig && match.y < .50); break;
      case kSigFrac90: var = int(match.isSig && match.y < .90); break;
      case kL0: var = match.photL0; break;
      case kL1: var = match.photL1; break;
      case kTrueEVis: var = match.trueEVis; break;
      case kQESigFrac:  var = (match.isSig && match.mode == simb::kQE);  break;
      case kResSigFrac: var = (match.isSig && match.mode == simb::kRes); break;
      case kDISSigFrac: var = (match.isSig && match.mode == simb::kDIS); break;
      case kQENCFrac:  var = (match.ccnc == simb::kNC && match.mode == simb::kQE);  break;
      case kResNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kRes); break;
      case kDISNCFrac: var = (match.ccnc == simb::kNC && match.mode == simb::kDIS); break;
      case kQENumuFrac:  var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kQE);  break;
      case kResNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kRes); break;
      case kDISNumuFrac: var = (match.ccnc == simb::kCC && abs(match.pdg) == 14 && match.mode == simb::kDIS); break;
      }

      xs.push_back(match.potential);
      ys.push_back(var);
      ws.push_back(w);

      tot += w;
      mean += w*var;
    } // end for i

    mean /= tot;

    if(ws.size() < 2){
      fit = -1;
    }
    else{
      double m, c;
      if(weightScheme == kUnweighted)
        util::LinFitUnweighted(xs, ys, m, c);
      else
        util::LinFit(xs, ys, ws, m, c);
      fit = m*xs[0]+c;
    }
  }

  //......................................................................
  void MakePID::produce(art::Event& evt)
  {
    std::unique_ptr<std::vector<lem::PIDDetails> > pidcol(new std::vector<lem::PIDDetails>);
    std::unique_ptr<art::Assns<lem::PIDDetails, rb::Cluster>> assns(new art::Assns<lem::PIDDetails, rb::Cluster>);

    art::Handle<std::vector<lem::MatchList> > matchlists;
    evt.getByLabel(fMatchesLabel, matchlists);

    art::Handle<std::vector<LEMInput> > inputs;
    evt.getByLabel(fSummarizeLabel, inputs);

    assert(inputs->size() == matchlists->size());

    art::FindManyP<rb::Cluster> fmt(matchlists, evt, fMatchesLabel);

    const unsigned int matchMax = matchlists->size();
    for(unsigned int matchIdx = 0; matchIdx < matchMax; ++matchIdx){
      MatchList matchlist = (*matchlists)[matchIdx];

      // Get all the slices associated with this match list
      std::vector<art::Ptr<rb::Cluster> > slices = fmt.at(matchIdx);
      // Associated slices can be empty when we're running over lemsum files,
      // where the slices have been stripped out. In that case, don't make
      // associations, LEMAssociator will sort it out later.
      assert(slices.empty() || slices.size() == 1);

      const double totalRecoGeV = (*inputs)[matchIdx].totalGeV;

      std::vector<MatchSummary> matches, enrichedMatches;
      for(unsigned int i = 0; i < matchlist.matches.size(); ++i){
        if(matchlist.matches[i].enrich){
          enrichedMatches.push_back(matchlist.matches[i]);
        }
        else{
          matches.push_back(matchlist.matches[i]);
        }
      }

      if(matches.empty()){
        lem::PIDDetails pid(-1);
        pidcol->push_back(pid);
        if(!slices.empty())
          util::CreateAssn(evt, *pidcol, slices[0], *assns);
        continue;
      }

      // Equalize stats
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig) matches[n].weight *= fSigScale;
        if(matches[n].ccnc == 1) matches[n].weight *= fNCScale;
      }
      for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
        enrichedMatches[n].weight *= fEnrichScale;
      }


      PIDExtraVars vars;

      std::vector<double> xs, ys, ws, ws2;

      double meanESig = 0;
      double meanEBkg = 0;
      double meanESigExp = 0;
      double meanEBkgExp = 0;
      double totSig = 0;
      double totBkg = 0;
      double totSigExp = 0;
      double totBkgExp = 0;
      double totW = 0;

      double invE = 0, invEDenom = 0;
      double expE = 0, expEDenom = 0;
      double sigE = 0, sigEDenom = 0;

      const double maxPotential = matches.empty() ? 1 : matches[matches.size()-1].potential;
      for(unsigned int i = 0; i < matches.size(); ++i){
        const double w = matches[i].weight;

        const double invEVal = w/(matches[i].potential+1e-5);
        if(matches[i].isSig) invE += invEVal; else invE -= invEVal;
        invEDenom += invEVal;

        const double expEVal = w*exp(-matches[i].potential);
        if(matches[i].isSig) expE += expEVal; else expE -= expEVal;
        expEDenom += expEVal;

        const double sigEVal = w/(1+matches[i].potential);
        if(matches[i].isSig) sigE += sigEVal; else sigE -= sigEVal;
        sigEDenom += sigEVal;

        xs.push_back(matches[i].potential);
        ys.push_back(matches[i].isSig);
        ws.push_back(w);
        double expVal = fExpTable[i];
        if(!fPositionalWeight){
          expVal = exp(-pow(matches[i].potential/maxPotential, fExpPow)/fExpLength);
        }
        ws2.push_back(w*expVal);

        if(matches[i].isSig){
          totSig += w;
          totSigExp += expVal*w;
          meanESig += w*matches[i].potential;
          meanESigExp += expVal*w*matches[i].potential;
        }
        else{
          totBkg += w;
          totBkgExp += expVal*w;
          meanEBkg += w*matches[i].potential;
          meanEBkgExp += expVal*w*matches[i].potential;
        }

        totW += w;
      } // end for i

      if(ws.size() < 2){
        vars.fChiSig = 0;
        vars.fChiBkg = 0;
        vars.fChiSigExp = 0;
        vars.fChiBkgExp = 0;
      }
      else{
        for(int c = 0; c <= 1; ++c){
          double num = 0, denom = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            num += ws[i]*(xs[i]*ys[i]-c*xs[i]);
            denom += xs[i]*xs[i]*ws[i];
          }
          const double m = num/denom;
          double chi = 0;
          double totW = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            chi += ws[i]*util::sqr(ys[i]-m*xs[i]-c);
            totW += ws[i];
          }
          chi /= totW;
          if(c == 0) vars.fChiBkg = chi;
          if(c == 1) vars.fChiSig = chi;
        } // end for c

        for(int c = 0; c <= 1; ++c){
          double num = 0, denom = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            num += ws2[i]*(xs[i]*ys[i]-c*xs[i]);
            denom += xs[i]*xs[i]*ws2[i];
          }
          const double m = num/denom;
          double chi = 0;
          double totW = 0;
          for(unsigned int i = 0; i < xs.size(); ++i){
            chi += ws2[i]*util::sqr(ys[i]-m*xs[i]-c);
            totW += ws2[i];
          }
          chi /= totW;
          if(c == 0) vars.fChiBkgExp = chi;
          if(c == 1) vars.fChiSigExp = chi;
        } // end for c
      }

      if(totSig) meanESig /= totSig;
      if(totBkg) meanEBkg /= totBkg;
      if(totSigExp) meanESigExp /= totSigExp;
      if(totBkgExp) meanEBkgExp /= totBkgExp;


      double junk;

      AvgAndFit(matches, kSimple,     kSigFrac, junk, vars.fPIDFit);
      AvgAndFit(matches, kUnweighted, kSigFrac, junk, vars.fPIDFitUnweight);
      AvgAndFit(matches, kExp,        kSigFrac, vars.fPIDExp, vars.fPIDExpFit);

      AvgAndFit(matches, kSimple, kY, vars.fMeanY, vars.fMeanYFit);
      AvgAndFit(matches, kExp,    kY, vars.fMeanYExp, junk);

      AvgAndFit(matches, kSigOnly,    kY, vars.fMeanYSig, vars.fMeanYSigFit);
      AvgAndFit(matches, kSigOnlyExp, kY, vars.fMeanYSigExp, junk);

      AvgAndFit(matches, kSimple, kQFrac, vars.fMeanQFrac, vars.fMeanQFracFit);
      AvgAndFit(matches, kExp,    kQFrac, vars.fMeanQFracExp, junk);

      AvgAndFit(matches, kSigOnly,    kQFrac, vars.fMeanQFracSig, vars.fMeanQFracSigFit);
      AvgAndFit(matches, kSigOnlyExp, kQFrac, vars.fMeanQFracSigExp, junk);

      AvgAndFit(matches, kBkgOnly,    kQFrac, vars.fMeanQFracBkg, vars.fMeanQFracBkgFit);
      AvgAndFit(matches, kBkgOnlyExp, kQFrac, vars.fMeanQFracBkgExp, junk);

      AvgAndFit(matches, kExp, kTrueEVis, vars.fMeanMatchEExp, junk);

      AvgAndFit(matches, kExp, kQESigFrac,   vars.fQESigFracExp,   junk);
      AvgAndFit(matches, kExp, kResSigFrac,  vars.fResSigFracExp,  junk);
      AvgAndFit(matches, kExp, kDISSigFrac,  vars.fDISSigFracExp,  junk);
      AvgAndFit(matches, kExp, kQENCFrac,    vars.fQENCFracExp,    junk);
      AvgAndFit(matches, kExp, kResNCFrac,   vars.fResNCFracExp,   junk);
      AvgAndFit(matches, kExp, kDISNCFrac,   vars.fDISNCFracExp,   junk);
      AvgAndFit(matches, kExp, kQENumuFrac,  vars.fQENumuFracExp,  junk);
      AvgAndFit(matches, kExp, kResNumuFrac, vars.fResNumuFracExp, junk);
      AvgAndFit(matches, kExp, kDISNumuFrac, vars.fDISNumuFracExp, junk);

      // Mix in the usual signal matches
      const int emOrig = enrichedMatches.size();
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig) enrichedMatches.push_back(matches[n]);
      }
      std::sort(enrichedMatches.begin(), enrichedMatches.end());
      enrichedMatches.resize(emOrig);

      AvgAndFit(enrichedMatches, kSimple, kBkgFrac, vars.fEnrichFrac, vars.fEnrichFracFit);
      AvgAndFit(enrichedMatches, kExp,    kBkgFrac, vars.fEnrichFracExp, junk);

      AvgAndFit(enrichedMatches, kBkgOnly,    kQFrac, vars.fEnrichQFrac, vars.fEnrichQFracFit);
      AvgAndFit(enrichedMatches, kBkgOnlyExp, kQFrac, vars.fEnrichQFracExp, junk);

      AvgAndFit(enrichedMatches, kSimple, kL0, junk, vars.fMeanL0Fit);
      AvgAndFit(enrichedMatches, kSimple, kL1, junk, vars.fMeanL1Fit);
      AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0Exp, junk);
      AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1Exp, junk);

      AvgAndFit(enrichedMatches, kPi0Only, kL0, junk, vars.fMeanL0PiOnlyFit);
      AvgAndFit(enrichedMatches, kPi0Only, kL1, junk, vars.fMeanL1PiOnlyFit);
      AvgAndFit(enrichedMatches, kExp, kL0, vars.fMeanL0PiOnlyExp, junk);
      AvgAndFit(enrichedMatches, kExp, kL1, vars.fMeanL1PiOnlyExp, junk);

      AvgAndFit(matches, kSimple, kSigFrac50, junk, vars.fPIDy50Fit);
      AvgAndFit(matches, kSimple, kSigFrac90, junk, vars.fPIDy90Fit);

      AvgAndFit(matches, kSig50Only, kQFrac, junk, vars.fMeanQFracSigY50Fit);
      AvgAndFit(matches, kSig90Only, kQFrac, junk, vars.fMeanQFracSigY90Fit);

      AvgAndFit(matches, kSig50Only, kY, junk, vars.fMeanYSigY50Fit);
      AvgAndFit(matches, kSig90Only, kY, junk, vars.fMeanYSigY90Fit);

      AvgAndFit(matches, kSimple, kNumuFrac, junk, vars.fPIDFracNumuFit);

      vars.fEnergyDiff = meanESig-meanEBkg;
      vars.fEnergyDiffExp = meanESigExp-meanEBkgExp;

      vars.fBestSigPotential = vars.fBestBkgPotential = vars.fBestEnrichPotential = 2;
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(matches[n].isSig){
          vars.fBestSigPotential = matches[n].potential;
          break;
        }
      }
      for(unsigned int n = 0; n < matches.size(); ++n){
        if(!matches[n].isSig){
          vars.fBestBkgPotential = matches[n].potential;
          break;
        }
      }
      for(unsigned int n = 0; n < enrichedMatches.size(); ++n){
        if(!enrichedMatches[n].isSig){
          vars.fBestEnrichPotential = enrichedMatches[n].potential;
          break;
        }
      }

      vars.fAvgInvE = invE/invEDenom;
      vars.fAvgExpE = expE/expEDenom;
      vars.fAvgSigE = sigE/sigEDenom;

      // if(fANN || fANNEnrich){
      //   fann_type input[6];
      //   input[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
      //   input[1] = vars.fMeanYExp;
      //   input[2] = vars.fMeanQFracExp;
      //   input[3] = vars.fEnergyDiffExp;

      //   if(fANN){
      //     fann_type* calc_out = fann_run(fANN, input);
      //     vars.fFannId = calc_out[0];
      //   }

      //   input[4] = 2*totalRecoGeV;
      //   input[5] = vars.fEnrichFracExp;

      //   if(fANNEnrich){
      //     fann_type* calc_out = fann_run(fANNEnrich, input);
      //     vars.fFannIdEnrich = calc_out[0];
      //   }
      // }

      if(fReader){
        gTMVAVar[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
        gTMVAVar[1] = vars.fMeanYExp;
        gTMVAVar[2] = vars.fMeanQFracExp;
        gTMVAVar[3] = vars.fEnergyDiffExp;
        if(fScaleCalE)
          gTMVAVar[4] = 2*totalRecoGeV;
        else
          gTMVAVar[4] = totalRecoGeV;
        gTMVAVar[5] = vars.fEnrichFracExp;

        vars.fTMVABDTG  = fReader->EvaluateMVA("BDTG");
        vars.fTMVABDT   = fReader->EvaluateMVA("BDT");
        vars.fTMVABDTD  = fReader->EvaluateMVA("BDTD");
        vars.fTMVABDTMF = fReader->EvaluateMVA("BDTMitFisher");
        vars.fTMVAMLP   = fReader->EvaluateMVA("MLP");

        vars.fTMVABDTG_CC = fReader->EvaluateMVA("BDTG_CC");
        vars.fTMVAMLP_CC  = fReader->EvaluateMVA("MLP_CC");
      }

      if(fForest){
        if(vars.fPIDExp == -1){
          vars.fDecTreePID = -1;
        }
        else{
          dec::Evt decevt;
          decevt.vars[0] = fUseFitInput ? vars.fPIDFit : vars.fPIDExp;
          decevt.vars[1] = vars.fMeanYExp;
          decevt.vars[2] = vars.fMeanQFracExp;
          decevt.vars[3] = vars.fEnergyDiffExp;
          if(fScaleCalE)
            decevt.vars[4] = 2*totalRecoGeV;
          else
            decevt.vars[4] = totalRecoGeV;
          decevt.vars[5] = vars.fEnrichFracExp;
          vars.fDecTreePID = fForest->Classify(decevt);
        }
      }

      double defaultPID = vars.fDecTreePID;

      // Flattened version of the Decision Tree variable
      if(fFlatHist) defaultPID = fFlatHist->GetBinContent(fFlatHist->FindBin(vars.fDecTreePID));


      lem::PIDDetails pid(defaultPID);

      if(fTrimMatchLists) matchlist.matches.resize(1);

      // TODO: sample size equalization doesn't survive this
      pid.matchList = matchlist;
      pid.vars = vars;

      pidcol->push_back(pid);

      if(!slices.empty()){
        util::CreateAssn(evt, *pidcol, slices[0], *assns);
      }
    } // end for matchIdx

    evt.put(std::move(pidcol));
    evt.put(std::move(assns));
  }

  DEFINE_ART_MODULE(MakePID)

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