EventListManipulator.cxx

Go to the documentation of this file.

//
//  EventListManipulator.cxx
//
//  Created by Brian Rebel on 3/21/16.
//

#include <stdexcept>
#include <memory>
#include <iomanip>
#include <fstream>

#include "TFile.h"

#include "messagefacility/MessageLogger/MessageLogger.h"

#include "CovarianceMatrixFit/core/Event.h"
#include "CovarianceMatrixFit/core/ShifterAndWeighter.h"
#include "CovarianceMatrixFit/core/VarVals.h"
#include "CovarianceMatrixFit/dataProducts/Constants.h"
#include "CovarianceMatrixFit/modules/EventListManipulator.h"
#include "CovarianceMatrixFit/utilities/CovarianceBinUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.h"
#include "CovarianceMatrixFit/utilities/SelectionUtility.h"

#include "nusimdata/SimulationBase/MCNeutrino.h"

namespace cmf{

  //----------------------------------------------------------------------------
  EventListManipulator::EventListManipulator(fhicl::ParameterSet const& p)
  {
    fEventCaps.clear();

    this->reconfigure(p);
  }

  //......................................................................
  EventListManipulator::~EventListManipulator()
  {
    fExposure.clear();
  }

  //......................................................................
  void EventListManipulator::reconfigure(const fhicl::ParameterSet& pset)
  {

    fRandom.SetSeed(pset.get<unsigned int>("RandomSeed", 0));

    auto detectors = pset.get<std::vector<std::string>>("DetectorsToUse", std::vector<std::string>({"NearDet", "FarDet"}));
    for(auto const& itr : detectors) fDetectors.insert(cmf::StringToDetectorType(itr));

    fCMFEventLabels     = pset.get< std::vector<std::string>         >("CMFEventLabels"           , std::vector<std::string>()         );
    fTreeDirectories    = pset.get< std::vector<std::string>         >("TreeDirectories"          , std::vector<std::string>()         );
    fUseEventId         = pset.get< bool                             >("UseEventId"               , false                              );
    fMaxNuEnergy[cmf::kNuESelection ] = pset.get< double             >("MaxNuEEnergyForAnalysis"  , 5.0                                );
    fMaxNuEnergy[cmf::kNuMuSelection] = pset.get< double             >("MaxNuMuEnergyForAnalysis" , 10.0                               );
    fMaxNuEnergy[cmf::kNCSelection  ] = pset.get< double             >("MaxNCEnergyForAnalysis"   , 20.0                               );
    fMaxEventsPerTree   = pset.get< float                            >("MaxEventsPerTree"         , 1.e6                               );
    fMinEventsPerTree   = pset.get< float                            >("MinEventsPerTree"         , 2.e4                               );
    fFillTextFile       = pset.get< bool                             >("FillTextFile"             , false                              );
    fDeserializeCosmics = pset.get< bool                             >("DeserializeCosmics"       , true                               );
    
    // if we're printing the event lists we want the event ID information!
    if (fFillTextFile){
      fUseEventId = true;
      fMaxEventsPerTree = 1.e6;
    }

    // for(auto & dir : fTreeDirectories) dir.append("/full");

    // determine which general experiments we have, ie nue or numu or both
    // allow users to force the use of all event lists if they choose
    fLoadAllEventLists = pset.get<bool>("LoadAllEventLists", false);

    this->SetExposures(pset);

  }

  // Set the exposures which are configured from the fhicl files
  //...........................................................................
  void EventListManipulator::SetExposures(const fhicl::ParameterSet &pset)
  {

    auto psetExposure   = pset.get< std::vector<fhicl::ParameterSet> >("ExposureToUse", std::vector<fhicl::ParameterSet>() );
    MF_LOG_DEBUG("EventListManipulator")
      << "There are "
      << psetExposure.size()
      << " exposures from fhicl";

    auto fhcExposureMultiplier = pset.get<double>("FHCExposureMultiplier", 1.);
    auto rhcExposureMultiplier = pset.get<double>("RHCExposureMultiplier", 1.);
    fExposure.clear();
    for(auto const& exposure : psetExposure){

      auto psetDet      = exposure.get<std::string>("Detector" );
      auto psetPeriod   = exposure.get<int        >("Period"   );
      auto psetExposure = exposure.get<float      >("Exposure" );
      auto psetLivetime = exposure.get<float      >("Livetime" );

      fPeriodsToUse.insert(psetPeriod);

      cmf::SpillSummary spillSummary(psetExposure,
                                     psetExposure,
                                     psetLivetime,
                                     0,
                                     0);

      cmf::DetType_t  det  = cmf::StringToDetectorType(psetDet);
      cmf::BeamType_t beam = PeriodToBeamType(std::to_string(psetPeriod));

      MF_LOG_DEBUG("ELM")
            << cmf::cDetType_Strings[det]
            << " "
            << cmf::cBeamType_Strings[beam]
            << " "
            << exposure.to_indented_string();

      // if we want high stats fakes, multiply the exposure by exposure 
      // multiplier - only ever applied to FD
      if(det == cmf::kFARDET)
        psetExposure *= (beam == cmf::kFHC) ? fhcExposureMultiplier : rhcExposureMultiplier;

      auto const& dbsSet = cmf::SelectionUtility::Instance()->SelectionsToUse();
      for(auto const& sel :  (*dbsSet.find(cmf::SelectionUtility::DetBeamSels(det, beam))).sels){
        cmf::MetaData md(false,
                         det,
                         cmf::kDataFile,
                         sel,
                         cmf::kUnknownInteraction,
                         psetPeriod);

        if(fExposure.count(md) < 1){
          fExposure.emplace(md, spillSummary);
        }
        else {
          fExposure[md].totalPOT += psetExposure;
          fExposure[md].goodPOT  += psetExposure;
          fExposure[md].liveTime += psetLivetime;
        }

      } // end loop over selections

    } // end if we have detectors and selections to make

    this->PrintPOTSummary     (fExposure);
    this->PrintLivetimeSummary(fExposure);

    //for(auto const& itr : fExposure)
    //  MF_LOG_VERBATIM("EventListManipulator")
    //      << "configuring exposure from fhicl"
    //      << std::setw(30)
    //      << itr.first.ToString()
    //      << " to have "
    //      << itr.second
    //      << " POT and this event list will not be deserialized";

  }

  // just print out Exposure summary
  // TODO - do we want this to work for MINOS data as well?
  //.....................................................................
  void EventListManipulator::PrintExposureSummary(cmf::ExposureMap &exposureMap, bool isLiveTime)
  {

    std::map<int, double> NDFHC;
    std::map<int, double> FDFHC;
    std::map<int, double> NDRHC;
    std::map<int, double> FDRHC;
    std::vector<int> ps;

    for (auto const& itr : exposureMap)
    {
      MF_LOG_DEBUG("ELM")
            << itr.first.ToString()
            << " "
            << itr.second;
      
      double value;
      if (isLiveTime) value = itr.second.liveTime;
      else            value = itr.second.goodPOT;

      // ND FHC
      if (itr.first.detector == cmf::kNEARDET && itr.first.BeamType() == cmf::BeamType_t::kFHC){
        if (NDFHC.find(itr.first.period) == NDFHC.end())
          NDFHC[itr.first.period] = value;
        if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
      }
      // FD FHC
      if (itr.first.detector == cmf::kFARDET && itr.first.BeamType() == cmf::BeamType_t::kFHC){
        if (FDFHC.find(itr.first.period) == FDFHC.end())
          FDFHC[itr.first.period] = value;
        if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
      }

      // ND RHC
      if (itr.first.detector == cmf::kNEARDET && itr.first.BeamType() == cmf::BeamType_t::kRHC){
        if (NDRHC.find(itr.first.period) == NDRHC.end())
          NDRHC[itr.first.period] = value;
        if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
      }

      // FD RHC
      if (itr.first.detector == cmf::kFARDET && itr.first.BeamType() == cmf::BeamType_t::kRHC){
        if (FDRHC.find(itr.first.period) == FDRHC.end())
          FDRHC[itr.first.period] = value;
        if (std::find(ps.begin(), ps.end(), itr.first.period) == ps.end()) ps.push_back(itr.first.period);
      }
    }

    if (isLiveTime){
      MF_LOG_VERBATIM("EventListManipulator")
        << "\n-------------------------------------------------------"
        << "\nLIVETIME SUMMARY"
        << "\n-------------------------------------------------------";
    }
    else{
      MF_LOG_VERBATIM("EventListManipulator")
        << "\n-------------------------------------------------------"
        << "\nPOT SUMMARY"
        << "\n-------------------------------------------------------";
    }

    std::sort(ps.begin(), ps.end());
    MF_LOG_VERBATIM("EventListManipulator")
      << std::setw(20)
      << ""
      << std::setw(20)
      << "ND FHC"
      << std::setw(20)
      << "FD FHC"
      << std::setw(20)
      << "ND RHC"
      << std::setw(20)
      << "FD RHC";
      for (size_t i = 0;  i < ps.size(); ++i){
        MF_LOG_VERBATIM("EventListManipulator")
          << std::setw(20)
          << "Period "+std::to_string(ps[i])
          << std::setw(20)
          << NDFHC[ps[i]]
          << std::setw(20)
          << FDFHC[ps[i]]
          << std::setw(20)
          << NDRHC[ps[i]]
          << std::setw(20)
          << FDRHC[ps[i]];
      }

  }

  //----------------------------------------------------------------------------
  void EventListManipulator::PrintPOTSummary(cmf::ExposureMap &exposureMap){
    this->PrintExposureSummary(exposureMap, false);
  }

  //----------------------------------------------------------------------------
  void EventListManipulator::PrintLivetimeSummary(cmf::ExposureMap &exposureMap){
    this->PrintExposureSummary(exposureMap, true);
  }

  // Skip any event list where the selection and detector combination were
  // not requested in the configuration
  //......................................................................
  bool EventListManipulator::UseEventsFromMetaData(cmf::MetaData            const& md,
                                                   cmf::DataMC_t            const  dataMC,
                                                   std::set<cmf::DetType_t> const& detectors)
  {
    if(fLoadAllEventLists) return true;

    // if we don't want cosmic muons return false right away, 
    // if we do, keep doing checks to be sure this is a selection 
    // we want
    if(md.fileType == cmf::kCosmicBackgroundFile && 
       !fDeserializeCosmics) return false;
      
    // check if this is data and we only want MC events
    // or vice versa
    if(dataMC == cmf::kMC && !md.isMC){
      MF_LOG_DEBUG("ELMUEFMMC")
          << md.ToString()
          << " failed cmf::kMC and !md.isMC";
      return false;
    }
    else if(dataMC == cmf::kData && md.isMC){
      MF_LOG_DEBUG("ELMUEFMRealData")
      << md.ToString()
      << " failed cmf::kData and md.isMC";
      return false;
    }
    else if(dataMC == cmf::kFakeData && !md.isMC){
      MF_LOG_DEBUG("ELMUEFMFakeData")
      << md.ToString()
      << " failed cmf::kFakeData and !md.isMC";
      return false;
    }

    //Skip unless the epoch is in the fPeriodsToUse set
    if(fPeriodsToUse.count(md.Period()) < 1){
      MF_LOG_DEBUG("ELMUEFMPeriod")
      << md.ToString()
      << " failed count of periods to use "
      << md.Period();
      return false;
    }

    // Skip if uncategorized interactions
    if(md.interactionType == cmf::kUncategorised){
      MF_LOG_DEBUG("ELMUEFMInteraction")
          << md.ToString()
          << " failed interaction type check to use "
          << md.interactionType;
      return false;
    }

    // check if the detector for this metadata is desired, return false if it
    // isn't in the fDetectors set.  We have to check detector separately from
    // the SelectionUtility because the selection utility is set for the entire
    // job, not just the EventListManipulator which may only want to use
    // event from a subset of detectors for the job.
    // Moreover, even though ultimately we want to use both detectors, a job
    // may want to load one at a time for memory management, so we use the
    // detectors set to check on that
    if(fDetectors.count(md.detector) < 1 || detectors.count(md.detector) < 1) return false;

    // now check that the selection is desired for the detector, beam and
    // selection of this metadata
    if(!cmf::SelectionUtility::Instance()->UsesDetBeamAndSel(md)){
      MF_LOG_DEBUG("ELMUEFMDetector")
          << md.ToString()
          << " failed to find "
          << md.DetectorString()
          << " "
          << md.PeriodString()
          << " "
          << cmf::cSelectionType_Strings[md.selectionType]
          << " in selections to use";
      return false;
    }

    // Metadata types we are going to skip
    // - ND flux swap for now
    //if(md->fileType == cmf::kSwap            &&
    //   md->detector == cmf::kNEARDET ) return false;
    
    return true;
  }

  // We want to make sure that the pre-configured cap does not result in
  // a poor sampling of events, so pass in the number of events in the
  // tree as well to check how the cap will affect the final number of events
  // used.  If too few events will remain after that cap, up it to ensure
  //......................................................................
  double EventListManipulator::PeriodEventCap(cmf::MetaData const& md,
                                              long                 treeEvents)
  {

    double cap = (treeEvents > fMaxEventsPerTree - 1) ? fMaxEventsPerTree / (1. * treeEvents) : 1.;

    MF_LOG_DEBUG("EVentListManipulator")
        << "Check event cap for "
        << md.ToString()
        << ": "
        << treeEvents
        << " / "
        << fMinEventsPerTree
        << " / "
        << fMaxEventsPerTree
        << " final cap: "
        << cap;

    // if no cap is specified, use all the events
    return cap;
  }

  //......................................................................
  void EventListManipulator::ExtractFromFile(TFile*                                metadataFile,
                                             std::string                    const& dirName,
                                             cmf::DataMC_t                  const& dataMC,
                                             std::set<cmf::DetType_t>       const& detectors,
                                             std::vector<cmf::FileExtracts>      & extractVec)
  {
    MF_LOG_DEBUG("EventListManipulator")
        << "extracting metadata from file";

    std::string mdString;
    std::string treeName;

    TTree *metadataTree = dynamic_cast<TTree *>(metadataFile->Get((dirName + "/metadata").c_str()));

    if(!metadataTree){
      throw cet::exception("EventListManipulator")
          << "cannot find metadata tree: "
          << dirName + "/metadata";
    }

    // The metadata tree may have the same metadata information in it over and over again
    // because we hadd'ed a bunch of files together.  We don't want to read in the event
    // trees multiple times, so we will make a map of the metadata and spill summary
    // information so that we only have one entry for each metadata
    cmf::ExposureMap exposureMap;

    // loop over the entries in the metadataTree
    cmf::MetaData     *md = nullptr;
    cmf::SpillSummary *ss = nullptr;
    metadataTree->SetBranchAddress("metadata",     &md);
    metadataTree->SetBranchAddress("spillsummary", &ss);

    // We have potentially hadd'ed several files together to get the
    // eventlisttree, so first loop over the entries in the metadata
    // tree and sum all the POT values for each metadata type
    for(int ctr = 0; ctr < metadataTree->GetEntries(); ++ctr){
      metadataTree->GetEntry(ctr);

      // now check if we want this metadata for our analysis
      // and find the minimum POT for all metadata having the same
      // beam/detector/file/interaction types for concatenated numu or nue
      // selections
      if(this->UseEventsFromMetaData(*md, dataMC, detectors)){
        if(exposureMap.count(*md) < 1) exposureMap[*md] = *ss;
        else exposureMap[*md] += *ss;
      }
    } // end loop over the metadata tree

    for(auto const& itr : exposureMap){

      // get the pointer to the event TTree
      mdString = itr.first.ToString();

      treeName = dirName + "/" + mdString;

      if(itr.second.goodPOT  > 0 ||
         (itr.first.fileType == kCosmicBackgroundFile && itr.second.liveTime > 0)){
        extractVec.emplace_back(itr.first,
                                itr.second,
                                treeName,
                                this->PeriodEventCap(itr.first, dynamic_cast<TTree*>(metadataFile->Get(treeName.c_str()))->GetEntriesFast()));
      }

      MF_LOG_DEBUG("EventListManipulator")
            << "deserializing: "
            << itr.first.ToString()
            << " "
            << itr.second.goodPOT
            << " "
            << itr.second.liveTime
            << " "
            << exposureMap.size();
    } // end loop over the metadata to spill summary map


    // loop over the exposureMap map one more time to check which
    // selections are being used
    // for(auto const& itr : exposureMap){
    //   MF_LOG_VERBATIM("EventListManipulator")
    //  << "Using "
    //  << itr.first.ToString()
    //  << " in the analysis";
    // }

  }

  // Use this method to load the TTrees containing the cmf::EventLists
  // until we are able to store them directly in an art file
  //
  // If onlyMC is true, only load the MC lists from the file
  //......................................................................
  void EventListManipulator::Deserialize(cmf::EventListColl            & eventLists,
                                         cmf::DataMC_t                   dataMC,
                                         std::set<cmf::DetType_t> const& detectors)
  {
    cmf::DataMC_t dataType = (dataMC == cmf::kData) ? dataMC : cmf::kMC;

    std::vector<cmf::FileExtracts> extractVec;
    size_t totalEvents = 0;

    cmf::MetaData     md;
    cmf::SpillSummary ss;
    double            normalization = 1.;

    // each input file is assumed to have a metadata tree and a cmfEvent tree
    for(auto & fileName : fCMFEventLabels){

      TFile *f = TFile::Open(fileName.c_str());

      if(f->IsOpen())
        MF_LOG_VERBATIM("EventListManipulator")
            << "successfully opened "
            << fileName;
      else
        throw cet::exception("EventListManipulator")
            << "Cannot open file "
            << fileName;

      for(auto const& dirName : fTreeDirectories){

        this->ExtractFromFile(f, dirName, dataType, detectors, extractVec);

        MF_LOG_DEBUG("EventListManipulator")
            << "There are "
            << extractVec.size()
            << " metadata categories";
      } // end loop over directories

      // map the detector/file/beam/interaction type + general selection type to the minimum exposure for that
      // combination. Use the ratio between the capped exposure and the found minimum to adjust the cap for
      // the other lists making up the general selection.
      // That is, new cap = (capped exposure for sub list) / (capped exposure for minimum exposure of sub lists)
      std::map<long, double> keyToMaxExposure;
      long   key;
      double cappedExposure;

      for(auto & itr : extractVec){
        key = itr.Key();

        // don't waste precision on the exponent for the POT
        cappedExposure = (itr.metadata.IsCosmicMuon()) ? itr.CappedLiveTime() : itr.CappedPOT() * 1.e-20;

        if(keyToMaxExposure.count(key) < 1)
          keyToMaxExposure[key] = cappedExposure;
        else
          keyToMaxExposure[key] = std::max(keyToMaxExposure.find(key)->second, cappedExposure);

        MF_LOG_DEBUG("EventListManipulator")
            << "input key: "
            << cmf::KeyToString(itr.metadata.Key())
            << " concat key: "
            << cmf::KeyToString(key)
            << " uncapped POT: "
            << itr.spillSummary.goodPOT
            << " uncapped LiveTime: "
            << itr.spillSummary.liveTime
            << " /  capped exposure (x 1e20) "
            << keyToMaxExposure[key];
      }

      // now loop over the extracted information again to fill the event lists
      // cosmic ray background livetime normalization is done when the spectra are
      // filled
      double exposure;
      for(auto & itr : extractVec){

        normalization = 1.;

        // need a non-const spill summary for filling the event list in case we have to
        // apply a cap. Any adjustment of the exposure due to capping will be done in the
        // FillEventList method
        ss = itr.spillSummary;

        // check if we are concatenating either the numu or nue selections
        // and make the appropriate metadata for that case.  We need to also
        // adjust the event cap for this particular tree and change the
        // total POT for this selection to account for normalization differences
        // between the selections to be concatenated
        md = itr.metadata;

        exposure = (md.IsCosmicMuon()) ? itr.CappedLiveTime() : itr.CappedPOT() * 1.e-20;

        if(cmf::SelectionUtility::Instance()->UsesSelection(cmf::kNuMuSelection) &&
           cmf::IsNuMuQuantiles(md.selectionType)                                ){
          md = cmf::MetaData(md.isMC,
                             md.detector,
                             md.fileType,
                             cmf::kNuMuSelection,
                             md.interactionType,
                             md.Period());

          normalization = keyToMaxExposure[itr.Key()] / exposure;
        }
        else if(cmf::SelectionUtility::Instance()->UsesSelection(cmf::kNuESelection) &&
                cmf::IsNuESelected(md.selectionType)                                 ){
          md = cmf::MetaData(md.isMC,
                             md.detector,
                             md.fileType,
                             cmf::kNuESelection,
                             md.interactionType,
                             md.Period());

          normalization = keyToMaxExposure[itr.Key()] / exposure;
        }
        
        //if(md.IsNuMuSelected())
        MF_LOG_DEBUG("EventListManipulator")
            << "Filling "
            << itr.metadata.ToString()
            << " / "
            << md.ToString()
            << " with cap "
            << itr.eventCap * normalization
            << " Exposure "
            << exposure * normalization
            //<< " : Events "
            //<< itr.eventTree->GetEntries()
            << " uncapped POT/livetime "
            << ss.goodPOT
            << " / "
            << ss.liveTime;
        
        ss *= itr.eventCap * normalization;
              
        auto listItr = cmf::FindEventList(md, eventLists);
        if(listItr == eventLists.end()){
          eventLists.emplace_back(cmf::EventList(md, ss));
          listItr = cmf::FindEventList(md, eventLists);
        }

        this->FillEventList(f, itr.treeName, itr.eventCap * normalization, *listItr);

        totalEvents += listItr->TotalEvents();
      } // end loop over extracted trees

      f->Close();
    }// end loop over files

    MF_LOG_VERBATIM("EventListManipulator")
        << "There are "
        << eventLists.size()
        << " " 
        << cmf::cDataTypeStrings[dataMC]
        << " event lists and "
        << totalEvents
        << " weighted events";

    std::map<cmf::MetaData, double> totalPOTPerSel;

    for(auto const& itr : eventLists){

      md = itr.ListMetaData();
      ss = itr.ListSpillSummary();

      if (md.fileType == kCosmicBackgroundFile){
        MF_LOG_DEBUG("EventListManipulator")
            << std::setw(55)
            << md.ToString()
            << " of period "
            << md.Period()
            << " "
            << std::setw(15)
            << " represents "
            << ss.liveTime
            << " livetime and "
            << itr.TotalEvents()
            << " events "
            << itr.TotalEvents() / ss.liveTime;
      }
      else{
        MF_LOG_DEBUG("EventListManipulator")
            << std::setw(55)
            << " "
            << std::setw(15)
            << " represents "
            << ss.goodPOT
            << " POT and "
            << itr.TotalEvents()
            << " events "
            << itr.TotalEvents() / ss.goodPOT;
      }

      cmf::SelectionType_t                                   selType = cmf::kNuESelection;
      if      (md.selectionType == cmf::kNuMuSelectionQ1)    selType = cmf::kNuMuSelection;
      else if (md.selectionType == cmf::kNuESelectionLowPID) selType = cmf::kNuESelection;
      else if (md.selectionType == cmf::kNCSelection)        selType = cmf::kNCSelection;
      else continue;

      cmf::MetaData mdtmp (md.MCKey(), 
                           md.detector, 
                           md.fileType, 
                           selType, 
                           cmf::kUnknownInteraction, 
                           md.BeamType() == kFHC ? 1 : 4);

      if (md.interactionType != cmf::kNuMuCC) continue;
      
      if (md.fileType != cmf::kBeam && 
          md.fileType != cmf::kSwap && 
          md.fileType != cmf::kTauSwap) continue;

      if (totalPOTPerSel.find(mdtmp) == totalPOTPerSel.end())
        totalPOTPerSel[mdtmp] = ss.goodPOT;
      else 
        totalPOTPerSel[mdtmp] += ss.goodPOT;

    }

    for (auto const& itr : totalPOTPerSel){
      if (itr.first.period == 1 || itr.first.period == 4){
        MF_LOG_VERBATIM("EventListManipulator")
          << "-- TOTAL "
                << cmf::cBeamType_Strings[itr.first.BeamType()]
                << " POT: \n"
                << cDetType_Strings[itr.first.detector]
                << " "
                << cFileTypeStrings[itr.first.fileType]
                << " "
                << cSelectionType_Strings[itr.first.selectionType]
                << " POT: "
                << itr.second; 
      }
    }

  }

  //--------------------------------------------------------------------------------------------
  void EventListManipulator::FillTextFile(cmf::MetaData   const& md,
                                          cmf::EventId         & ev,
                                          cmf::DataVars        & dv,
                                          cmf::TruthVars       & tv,
                                          cmf::WeightVars      & wv)
  {

    std::string fileName = md.DetectorString();
    fileName += cBeamType_Strings[md.BeamType()];
    fileName += cSelectionType_Strings[md.selectionType];
    fileName += cFileTypeStrings[md.fileType];
    fileName += ".txt";

    std::ifstream ifile;
    std::ofstream ofile;
    ifile.open(fileName);
    if(!ifile) {
      ofile.open(fileName, std::ofstream::out);
      ofile 
        << "\n"
        //<< md.ToString()
        //<< " "
        << ev.run 
        << " "
        << ev.subrun
        << " "
        << ev.event
        << " "
        << ev.cycle
        << " "
        << ev.slice
        << " "
        << dv.fLep_RecoE
        << " "
        << dv.fHad_RecoE
        << " "
        << dv.fLep_RecoE + dv.fHad_RecoE
        << " "
        << dv.fNu_RecoE
        << " "
        << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
        << " ";
      if (md.isMC){
        ofile
          << tv.fTrueE
          << " "
          << tv.fTruePDG
          << " "
          << tv.fTrueCCNC
          << " "
          << wv.fXSecCV2020_Weight
          << " "
          << wv.fPPFXFluxCV_Weight;
      }


      //ofile
      //    << "run"
      //    << " subrun"
      //    << " ev"
      //    << " cyc"
      //    << " slc"
      //    << " Ereco";
      //  if (md.isMC){
      //    ofile
      //      << " Etrue"
      //      << " pdg"
      //      << " ccnc"
      //      << " cvwgt";
      //  }
    }
    ofile.open(fileName, std::ofstream::out | std::ofstream::app);

    ofile
      << "\n"
      //<< md.ToString()
      //<< " "
      << ev.run
      << " "
      << ev.subrun
      << " "
      << ev.event
      << " "
      << ev.cycle
      << " "
      << ev.slice
      << " "
      << dv.fLep_RecoE
      << " "
      << dv.fHad_RecoE
      << " "
      << dv.fHad_RecoE + dv.fLep_RecoE
      << " "
      << dv.fNu_RecoE
      << " "
      << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
      << " ";
      if (md.isMC){
        ofile
          << tv.fTrueE
          << " "
          << tv.fTruePDG
          << " "
          << tv.fTrueCCNC
          << " "
          << wv.fXSecCV2020_Weight
          << " "
          << wv.fPPFXFluxCV_Weight;
      }

    ofile.close();
  }

  //......................................................................
  void EventListManipulator::FillEventList(TFile               * eventFile,
                                           std::string    const& treeName,
                                           double         const& eventCap,
                                           cmf::EventList      & eventList)
  {
    cmf::EventId    eventId;
    cmf::DataVars   dataVars;
    cmf::TruthVars  truthVars;
    cmf::WeightVars weightVars;

    auto const& md = eventList.ListMetaData();

    std::vector<unsigned char>* systVarKeys       = nullptr; // PROD5
    std::vector<float>*         systVarMinus2Wgts = nullptr; // PROD5
    std::vector<float>*         systVarMinus1Wgts = nullptr; // PROD5
    std::vector<float>*         systVarPlus1Wgts  = nullptr; // PROD5
    std::vector<float>*         systVarPlus2Wgts  = nullptr; // PROD5

    auto const sysKeys = cmf::ParameterUtility::Instance()->SysParKeys();
    // for(auto const& itr : sysKeys)
    //   MF_LOG_VERBATIM("EventListManipulator")
    //   << itr
    //   << " "
    //   << cmf::KeyToVarName(itr);

    TTree *eventTree = dynamic_cast<TTree*>(eventFile->Get(treeName.c_str()));
    eventTree->SetBranchAddress("dataVars", &dataVars);
    if(fUseEventId)
      eventTree->SetBranchAddress("eventId",  &eventId);

    if(md.isMC){
      eventTree->SetBranchAddress("truthVars" , &truthVars );
      eventTree->SetBranchAddress("weightVars", &weightVars);

      eventTree->SetBranchAddress("systVarKeys",       &systVarKeys);
      eventTree->SetBranchAddress("systVarMinus2Wgts", &systVarMinus2Wgts);
      eventTree->SetBranchAddress("systVarMinus1Wgts", &systVarMinus1Wgts);
      eventTree->SetBranchAddress("systVarPlus1Wgts",  &systVarPlus1Wgts);
      eventTree->SetBranchAddress("systVarPlus2Wgts",  &systVarPlus2Wgts);
    } // end if MC

    // get the generic identifier for numu or nue selection so
    // we can then determine the maximum energy of the neutrino
    // events to use
    cmf::SelectionType_t selection = cmf::kNuMuSelection;
    if     ( md.IsNuESelected() ) selection = cmf::kNuESelection;
    else if( md.IsNCSelected()  ) selection = cmf::kNCSelection;

    MF_LOG_DEBUG("EventListManipulator")
        << md.ToString()
        << " selection is "
        << selection
        << " max energy: "
        << fMaxNuEnergy.find(selection)->second;

    // check that our event list corresponds to a reasonable exposure
    if(eventList.ListSpillSummary().goodPOT < 1.e17 && md.fileType != kCosmicBackgroundFile)
      throw cet::exception("EventListManipulator")
          << "exposure for "
          << eventList.ListMetaData().ToString()
          << " is "
          << eventList.ListSpillSummary().goodPOT
          << " < 10^9 POT.  That cannot be right, bail";

    auto    systs = std::make_unique<cmf::SystVarColl>();
    uint8_t loc;

    MF_LOG_DEBUG("ELM")
      << "There are "
      << eventTree->GetEntries()
      << " events in tree before capping, "
      << eventList.size()
      << " events in list already, expect "
      << eventTree->GetEntries() * eventCap + eventList.size()
      << " events in list after capping";
        
    // use maxEvents to adjust the number of entries loaded in testing
    long maxEvents = eventTree->GetEntries();
    for(long evNum = 0; evNum < maxEvents; ++evNum){

      // If the uniformly distributed random number is greater than the
      // cap forget the current event.  This has the correct boundry
      // conditions - if the cap is 1, then all events are used; if it is
      // 0 then none are.

      // At some point we may need to worry about keeping more events in the
      // tails of the spectrum, but for now it looks like we can keep at least
      // 80% of all MetaData options, so we should not be biasing the tails
      if(fRandom.Rndm() > eventCap) continue;

      eventTree->GetEntry(evNum);

      if (fFillTextFile)
        this->FillTextFile(md, eventId, dataVars, truthVars, weightVars);

      // don't bother filling events outside the energy range for 
      // this selection - it just wastes time in the fitting
      if(dataVars.fNu_RecoE > fMaxNuEnergy.find(selection)->second || 
         dataVars.fNu_RecoE < cmf::CovarianceBinUtility::Instance()->SelectionLowEdge(md)) continue;

      // set the reco energy bin for this event
      dataVars.fRecoEBin  = cmf::CovarianceBinUtility::Instance()->EnergyToBin(dataVars.fNu_RecoE, md);
      dataVars.fRecoEBin += cmf::CovarianceBinUtility::Instance()->KeyToOffset(md.DetectorBeamSelectionKey());

      // now set the dataVals
      cmf::DataVarVals dataVals(dataVars);

      // now we can fill the rest of the event info because we are
      // keeping this event
      if(md.isMC){
              systs->clear();
      
              // if(systVarKeys->size() < 1)
              //        MF_LOG_VERBATIM("EventListManipulator")
              //          << md.ToString();
              
              for(size_t v = 0; v < systVarKeys->size(); ++v){
                loc = (*systVarKeys)[v];
                
                // MF_LOG_VERBATIM("EventListManipulator")
                //   << (*systVarKeys)[v]
                //   << " "
                //   << cmf::KeyToVarName(loc);
                
          // check if this sytematic is something we want to evaluate in the job, 
          // ignore it if not
                if(sysKeys.find(loc) == sysKeys.end()) continue;
                
                // don't bother saving a bunch of 1s in the SystVarColl, instead just
                // return all 1s if that systematic is desired for this event from MCVars.
                if((*systVarMinus2Wgts)[v] == 1 &&
                   (*systVarMinus1Wgts)[v] == 1 &&
                   (*systVarPlus1Wgts )[v] == 1 &&
                   (*systVarPlus2Wgts )[v] == 1 ) continue;
                
                systs->emplace_back(loc,
                              (*systVarMinus2Wgts)[v],
                              (*systVarMinus1Wgts)[v],
                              (*systVarPlus1Wgts )[v],
                              (*systVarPlus2Wgts )[v]);
                
                MF_LOG_DEBUG("ELMVarCheck")
                  << "check systematic vector sizes for event "
                  << evNum
                  << ": " 
                  << systVarKeys->size()
                  << " "
                  << systVarMinus2Wgts->size()
                  << " "
                  << systVarMinus1Wgts->size()
                  << " "
                  << systVarPlus1Wgts->size()
                  << " "
                  << systVarPlus2Wgts->size()
                  << " "
                  << systs->size()
                  << " "
                  << systs->back();
              }
        
              if(fUseEventId)
                eventList.AddEvent(eventId,
                                         dataVals,
                                         cmf::MCVarVals(truthVars, weightVars, *systs));
              else  
                eventList.AddEvent(dataVals,
                                         cmf::MCVarVals(truthVars, weightVars, *systs));
      } // end if md.isMC
      else{

              // Print out the event id information for each data event in the NOvA FD
              if(md.detector == cmf::kFARDET){
          
                std::string selectionStr("numu");
                if     (md.IsNuESelected()) selectionStr = "nue";
                else if(md.IsNCSelected() ) selectionStr = "nc";
                MF_LOG_DEBUG("EventListManipulator")
                  << selectionStr
                  << " "
                  << eventId;
      
              } // end if FD to print the event information
        
              if(fUseEventId)
                eventList.AddEvent(eventId,
                                   dataVals);
              else
                eventList.AddEvent(dataVals);
      }
                
    } // end loop over events

    // loop over the events we just added and check the systematics 
    // were loaded correctly
    // for(auto const& evt : eventList){
    //   MF_LOG_VERBATIM("EventListManipulator")
    //     << "event has systematic uncertainties with values ";
    //   for(auto const& sitr : evt->MCVals().SystematicShifts())
    //     MF_LOG_VERBATIM("EventListManipulator")
    //       << "\t" 
    //       << sitr; 
    // }

    systs.reset();
    eventTree = nullptr;
  } // done

} // end namespace