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/SelectionUtility.h"
#include "CovarianceMatrixFit/utilities/ParameterUtility.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)
  {

    LOG_VERBATIM("EventListManipulator") 
      << "\n-------------------------------------------------------"
      << "\nConfiguring production 5 event lists"
      << "\n-------------------------------------------------------";

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

    auto detectors = pset.get<std::vector<std::string>>("DetectorsToUse", std::vector<std::string>({"Near", "Far"}));
    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);

    // Comment out using event caps for now
    //this->SetEventCaps(pset);

  }

  // limits for the fraction of events to use from each MD type
  // the configuration uses a set fraction for each MD type initially,
  // but this number can be adjusted when the caps are applied if a given
  // MD type has too few events that way.  Too few is hardcoded to be 2.e4
  //...........................................................................
  void EventListManipulator::SetEventCaps(const fhicl::ParameterSet &pset)
  {
    auto caps = pset.get<fhicl::ParameterSet>("EventCap", fhicl::ParameterSet());
    std::string keyStr;

    // the event caps are parameter sets named by the detector and file type, the
    // periods are stored in the parameters for each detector/file type combo
    for(auto const& name : caps.get_names()){

      LOG_DEBUG("ELMEventCap")
          << "cap parameter: "
          << name;

      // no flux or tau files for the ND
      if(name.find("Near") != std::string::npos &&
         name.find("Beam") == std::string::npos &&
         name.find("Data") == std::string::npos ) continue;

      for(auto const& per : fPeriodsToUse){
        keyStr = name + "P" + std::to_string(per);

        fEventCaps[keyStr] = caps.get<fhicl::ParameterSet>(name).get<float>("P" + std::to_string(per), 1.);

        LOG_DEBUG("ELMEventCap")
            << keyStr
            << " "
            << caps.get<fhicl::ParameterSet>(name).get<float>("P" + std::to_string(per), 1.);

      } // end loop over periods
    } // end loop over caps

  }

  // 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>() );
    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  = (psetDet == "NearDet") ? cmf::kNEARDET : cmf::kFARDET;
      cmf::BeamType_t beam = PeriodToBeamType(std::to_string(psetPeriod));

      // 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)
    //  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
  //.....................................................................
  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)
    {

      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){
      LOG_VERBATIM("EventListManipulator") 
        << "\n-------------------------------------------------------"
        << "\nLIVETIME SUMMARY"
        << "\n-------------------------------------------------------";
    }
    else{
      LOG_VERBATIM("EventListManipulator") 
        << "\n-------------------------------------------------------"
        << "\nPOT SUMMARY"
        << "\n-------------------------------------------------------";
    }

    std::sort(ps.begin(), ps.end());
    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){
        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){
      LOG_DEBUG("ELMUEFMMC")
          << md.ToString()
          << " failed cmf::kMC and !md.isMC";
      return false;
    }
    else if(dataMC == cmf::kData && md.isMC){
      LOG_DEBUG("ELMUEFMRealData")
      << md.ToString()
      << " failed cmf::kData and md.isMC";
      return false;
    }
    else if(dataMC == cmf::kFakeData && !md.isMC){
      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){
      LOG_DEBUG("ELMUEFMPeriod")
      << md.ToString()
      << " failed count of periods to use "
      << md.Period();
      return false;
    }

    // Skip if uncategorized interactions
    if(md.interactionType == cmf::kUncategorised){
      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)){
      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
  //......................................................................
  float EventListManipulator::PeriodEventCap(cmf::MetaData const& md,
                                             long                 treeEvents)
  {
    //auto key = md.DetectorFilePeriodString();

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

    // If we are going to have more than fMaxEventsPerTree in a tree,
    // we want to cap it at that value to keep the memory footprint down

    // Capping events from trees with fewer events than fMinEventsPerTree
    // doesn't make a lot of sense, so return 1 in that case

    // if(treeEvents < fMinEventsPerTree + 1){
    //   cap = 1.;
    // }
    // else if(treeEvents > fMaxEventsPerTree - 1){
    //   cap = fMaxEventsPerTree / (1. * treeEvents);
    // }
    // else if(treeEvents > fMinEventsPerTree &&
    //         treeEvents < fMaxEventsPerTree){

    //   // now we are in the in between situation - fewer than fMaxEventsPerTree
    //   // but more than fMinEventsPerTree.  

    //   if(fEventCaps.count(key) > 0){
    //     float testNum = fEventCaps[key] * treeEvents;

    //     // If the capped number would be smaller than the
    //     // minimum number of events we want, set the cap to
    //     // be equivalent to the minimum events to keep.
    //     // Otherwise just return the found cap
    //     if(testNum < fMinEventsPerTree){
    //       cap = (1. * fMinEventsPerTree) / (1. * treeEvents);
    //     }
    //     else{
    //       cap = fEventCaps[key];
    //     }
    //   } // end if there is a defined cap for this key
    // } // end if there are more than the minimum number of events to apply a cap

    LOG_DEBUG("EVentListManipulator")
        << "Check event cap for "
        << md.ToString()
        << ": "
        << treeEvents
        << " / "
        << fMinEventsPerTree
        << " / "
        << fMaxEventsPerTree
      //<< " cap fraction: "
      //<< fEventCaps[key]
        << " 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)
  {
    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()));
      }

      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){
    //   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())
        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);

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

      // map the detector/file/beam + general selection type to the minimum POT for that
      // combination to do the POT accounting for concatenation
      std::map<long, double> keyToMinExposure;
      long   key;
      double cappedExposure;

      for(auto & itr : extractVec){
        key            = itr.Key();
        cappedExposure = (itr.metadata.IsCosmicMuon()) ? itr.CappedLiveTime() : itr.CappedPOT();

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

        LOG_DEBUG("EventListManipulator")
            << "input key: "
            << cmf::KeyToString(itr.metadata.Key())
            << " concat key: "
            << cmf::KeyToString(key)
            << " POT: "
            << itr.spillSummary.goodPOT
            << " / "
            << keyToMinExposure[key]
            << " LiveTime: "
            << itr.spillSummary.liveTime
            << " / "
            << keyToMinExposure[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
        ss = itr.spillSummary;

        // Rescale POT based on the event cap, which is a fraction of the total
        // events available
        ss.goodPOT  *= itr.eventCap;
        ss.totalPOT *= itr.eventCap;
        ss.liveTime *= itr.eventCap;

        // 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()) ? ss.liveTime : ss.goodPOT;

        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 = keyToMinExposure[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 = keyToMinExposure[itr.Key()] / exposure;
        }
        //if (itr.metadata.fileType == kCosmicBackgroundFile)
        //  normalization = fPOTSum[md.Period()].liveTime / ss.liveTime;

        //if(md.IsNuMuSelected())
        LOG_DEBUG("EventListManipulator")
            << "Filling "
            << itr.metadata.ToString()
            << " / "
            << md.ToString()
            << " with cap "
            << itr.eventCap
            << " Exposure "
            << exposure / itr.eventCap
            //<< " : Events "
            //<< itr.eventTree->GetEntries()
            << " cap POT/livetime "
            << ss.goodPOT
                  << " / "
                  << ss.liveTime
            //<< " : Events "
            //<< itr.eventTree->GetEntries() * itr.eventCap
            << " with normalization of "
            << 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->size();
      } // end loop over extracted trees

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

    LOG_VERBATIM("EventListManipulator")
        << "There are "
        << eventLists.size()
        << " event lists and "
        << totalEvents
        << " events";

    for(auto const& itr : eventLists){

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

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

  //--------------------------------------------------------------------------------------------
  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 += ".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.fHad_RecoE + dv.fLep_RecoE
        << " "
        << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
        << " ";
      if (md.isMC){
        ofile
          << tv.fTrueE
          << " "
          << tv.fTruePDG
          << " "
          << tv.fTrueCCNC
          << " "
          << wv.fXSecCVPPFX_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
      << " "
      << cmf::RecoEnergy(dv.fLep_RecoE, dv.fHad_RecoE, md)
      << " ";
    if (md.isMC){
      ofile
        << tv.fTrueE
        << " "
        << tv.fTruePDG
        << " "
        << tv.fTrueCCNC
        << " "
        << wv.fXSecCVPPFX_Weight;
    }

    ofile.close();
  }

  //......................................................................
  void EventListManipulator::FillEventList(TFile               * eventFile,
                                           std::string    const& treeName,
                                           double         const& eventCap,
                                           double         const& normalization,
                                           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)
    //   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;

    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;

    // 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);

      // we have to fill the data vals first to check the
      // the energy of this event
      dataVars.fFake_Weight = normalization;
      cmf::DataVarVals dataVals(dataVars);

      // don't bother filling events with greater than the maximum
      // energy for the analysis in the trees - it just wastes time
      // in the fitting
      if(dataVals.val_at(cmf::kNu_RecoE, md) > fMaxNuEnergy.find(selection)->second) continue;

      // 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)
              //        LOG_VERBATIM("EventListManipulator")
              //          << md.ToString();
              
              for(size_t v = 0; v < systVarKeys->size(); ++v){
                loc = (*systVarKeys)[v];

                // 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(std::find(sysKeys.begin(), sysKeys.end(), 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]);
                
                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 FD
              if(md.detector == cmf::kFARDET){
      
                std::string selectionStr("numu");
                if     (md.IsNuESelected()) selectionStr = "nue";
                else if(md.IsNCSelected() ) selectionStr = "nc";
                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){
    //   LOG_VERBATIM("EventListManipulator")
    //     << "event has systematic uncertainties with values ";
    //   for(auto const& sitr : evt->MCVals().SystematicShifts())
    //     LOG_VERBATIM("EventListManipulator")
    //       << "\t" 
    //       << sitr; 
    // }

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

} // end namespace