SpectrumLoader.cxx

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"

#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/SAMProjectSource.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/Utilities.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include <cassert>
#include <iostream>
#include <cmath>

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

namespace ana
{
  //----------------------------------------------------------------------
  SpectrumLoader::SpectrumLoader(const std::string& wildcard, DataSource src)
    : SpectrumLoaderBase(wildcard, src)
  {
  }

  //----------------------------------------------------------------------
  SpectrumLoader::SpectrumLoader(const std::vector<std::string>& fnames,
                                 DataSource src)
    : SpectrumLoaderBase(fnames, src)
  {
  }

  //----------------------------------------------------------------------
  SpectrumLoader::SpectrumLoader(DataSource src)
    : SpectrumLoaderBase(src)
  {
  }

  //----------------------------------------------------------------------
  SpectrumLoader SpectrumLoader::FromSAMProject(const std::string& proj,
                                                DataSource src,
                                                int fileLimit)
  {
    SpectrumLoader ret;
    ret.fSource = src;
    ret.fWildcard = "project "+proj;
    ret.fFileSource = std::unique_ptr<IFileSource>(new SAMProjectSource(proj, fileLimit));
    return ret;
  }

  //----------------------------------------------------------------------
  SpectrumLoader::~SpectrumLoader()
  {
  }

  struct CompareByID
  {
    bool operator()(const Cut& a, const Cut& b)
    {
      return a.ID() < b.ID();
    }
  };

  //----------------------------------------------------------------------
  void SpectrumLoader::Go()
  {
    if(fGone){
      std::cerr << "Error: can only call Go() once on a SpectrumLoader" << std::endl;
      abort();
    }
    fGone = true;

    // Find all the unique cuts
    std::set<Cut, CompareByID> cuts;
    for(auto& shiftdef: fHistDefs)
      for(auto& cutdef: shiftdef.second)
        cuts.insert(cutdef.first);
    for(const Cut& cut: cuts) fAllCuts.push_back(cut);

    fLivetimeByCut.resize(fAllCuts.size());
    fPOTByCut.resize(fAllCuts.size());


    const int Nfiles = NFiles();

    Progress* prog = 0;

    caf::SRBranchRegistry::clear();

    int fileIdx = -1;
    while(TFile* f = GetNextFile()){
      ++fileIdx;

      if(Nfiles >= 0 && !prog){
        std::string sum = TString::Format("Filling %lu spectra", fHistDefs.TotalSize()).Data();
        if(Nfiles > 1){
          if(!fSpillHistDefs.empty()){
            sum += TString::Format(", %lu spillTree spectra", fSpillHistDefs.size()).Data();
          }
          if(!fNuHistDefs.empty()){
            sum += TString::Format(", %lu nuTree spectra", fNuHistDefs.size()).Data();
          }
        }

        sum += TString::Format(" from %d files matching '%s'", Nfiles, fWildcard.c_str()).Data();
        prog = new Progress(sum);
      }

      HandleFile(f, Nfiles == 1 ? prog : 0);

      if(Nfiles > 1 && prog) prog->SetProgress((fileIdx+1.)/Nfiles);
    } // end for fileIdx

    StoreExposures();

    if(prog){
      prog->Done();
      delete prog;
    }

    ReportExposures();

    fHistDefs.RemoveLoader(this);
    fHistDefs.Clear();
  }

  //----------------------------------------------------------------------
  void SpectrumLoader::HandleFile(TFile* f, Progress* prog)
  {
    assert(!f->IsZombie());

    // Test for flat (trees are inside a directory) or nested (tree is at top
    // level) cases.
    TDirectory* dir = 0;
    TTree* tr = 0;

    TObject* obj = f->Get("recTree");
    assert(obj); // Must exist in one form or the other

    // It might seem like you could use GetObject() to do the type-checking
    // directly, but that method seems to have a memory leak in the case that
    // the types don't match.
    if(obj->ClassName() == std::string("TTree")){
      // nested case
      tr = (TTree*)obj;
    }
    else{
      // FlatCAF case
      dir = (TDirectory*)obj;
      tr = (TTree*)dir->Get("rec");
      assert(tr);
    }

    long n;
    caf::SRProxy sr(dir, tr, "rec", n, 0);

    //    FloatingExceptionOnNaN fpnan;

    const long Nentries = tr->GetEntries();
    for(n = 0; n < Nentries; ++n){
      if(!dir) tr->LoadTree(n); // nested mode

      HandleRecord(&sr);

      if(prog && n%100 == 0) prog->SetProgress(double(n)/Nentries);
    } // end for n
  }

  //----------------------------------------------------------------------
  /// Helper for \ref SpectrumLoader::HandleRecord
  template<class T, class U> class CutVarCache
  {
  public:
    CutVarCache() : fVals(U::MaxID()+1), fValsSet(U::MaxID()+1, false) {}

    inline T Get(const U& var, const caf::SRProxy* sr)
    {
      const unsigned int id = var.ID();

      if(fValsSet[id]){
        return fVals[id];
      }
      else{
        const T val = var(sr);
        fVals[id] = val;
        fValsSet[id] = true;
        return val;
      }
    }

  protected:
    // Seems to be faster to do this than [unordered_]map
    std::vector<T> fVals;
    std::vector<bool> fValsSet;
  };

  //----------------------------------------------------------------------
  void SpectrumLoader::HandleRecord(caf::SRProxy* sr)
  {
    // First up, check this spill passes the global quality cuts
    if(fSpillCut && !(*fSpillCut)(&sr->spill)) return;

    // Some shifts only adjust the weight, so they're effectively nominal, but
    // aren't grouped with the other nominal histograms. Keep track of the
    // results for nominals in these caches to speed those systs up.
    CutVarCache<bool, Cut> nomCutCache;
    CutVarCache<double, Var> nomWeiCache;
    CutVarCache<double, Var> nomVarCache;

    // These aren't set in SRNeutrino, but we need them for flux systs, so copy
    // them over from SRSpill.
    //
    // TODO TODO - remove this hack (and the equivalent in
    // SpectrumLoaderBase.cxx) once they're available in the CAFs
    if(sr->mc.nnu > 0){
      sr->mc.nu[0].isFHC = bool(sr->spill.isFHC);
      sr->mc.nu[0].is0HC = bool(sr->spill.is0HC);
      sr->mc.nu[0].isRHC = bool(sr->spill.isRHC);
      sr->mc.nu[0].det   = caf::Det_t(sr->hdr.det);
    }

    for(auto& shiftdef: fHistDefs){
      caf::SRProxySystController::BeginTransaction();

      const SystShifts& shift = shiftdef.first;

      // Need to provide a clean slate for each new set of systematic shifts to
      // work from. Copying the whole StandardRecord is pretty expensive, so
      // modify it in place and revert it afterwards.

      bool shifted = false;

      double systWeight = 1;
      // Can special-case nominal to not pay cost of Shift()
      if(!shift.IsNominal()){
        shift.Shift(sr, systWeight);
        // If there were only weighting systs applied then the cached nominal
        // values are still valid.
        shifted = caf::SRProxySystController::AnyShifted();
      }

      for(auto& cutdef: shiftdef.second){
        const Cut& cut = cutdef.first;

        const bool pass = shifted ? cut(sr) : nomCutCache.Get(cut, sr);
        // Cut failed, skip all the histograms that depended on it
        if(!pass) continue;

        for(auto& weidef: cutdef.second){
          const Var& weivar = weidef.first;

          double wei = shifted ? weivar(sr) : nomWeiCache.Get(weivar, sr);

          if(wei < 0){
            std::cerr << "Negative weight " << wei
                      << " returned from Var";
            std::cerr << std::endl;
            abort();
          }

          wei *= systWeight;
          if(wei == 0) continue;

          for(auto& vardef: weidef.second){
            if(vardef.first.IsMulti()){
              for(double val: vardef.first.GetMultiVar()(sr)){
                for(Spectrum* s: vardef.second.spects)
                  s->Fill(val, wei);
              }
              continue;
            }

            const Var& var = vardef.first.GetVar();

            const double val = shifted ? var(sr) : nomVarCache.Get(var, sr);

            if(std::isnan(val) || std::isinf(val)){
              std::cerr << "Warning: Bad value: " << val
                        << " returned from a Var. The input variable(s) could "
                        << "be NaN in the CAF,  or perhaps your "
                        << "Var code computed 0/0?";
              std::cout << ". Not filling into this histogram for this slice." << std::endl;
              continue;
            }

            for(Spectrum* s: vardef.second.spects) s->Fill(val, wei);

            for(ReweightableSpectrum* rw: vardef.second.rwSpects){
              const double yval = rw->ReweightVar()(sr);

              if(std::isnan(yval) || std::isinf(yval)){
                std::cerr << "Warning: Bad value: " << yval
                          << " for reweighting Var";
                std::cout << ". Not filling into histogram." << std::endl;
                continue;
              }

              if(rw->fHistD) rw->fHistD->Fill(val, yval, wei);
              if(rw->fHistF) rw->fHistF->Fill(val, yval, wei);
            } // end for rw
          } // end for vardef
        } // end for weidef
      } // end for cutdef

      // Return StandardRecord to its unshifted form ready for the next
      // histogram.
      caf::SRProxySystController::Rollback();
    } // end for shiftdef
  }

  //----------------------------------------------------------------------
  void SpectrumLoader::ReportExposures()
  {
    // The POT member variables we use here were filled as part of
    // SpectrumLoaderBase::GetNextFile() as we looped through the input files.

    // Let's just assume no-one is using the Cut::POT() function yet, so this
    // printout remains relevant...

    if(fSpillCut && fRunPOT > 0){
      std::cout << "Total " << fRunPOT << " POT, of which "
                << fPOT << " (" << int(100*fPOT/fRunPOT+.5)
                << "%) passes quality cuts." << std::endl;
    }
    else{
      std::cout << "Total " << fRunPOT << " POT." << std::endl;
      fPOT = fRunPOT;
    }

    std::vector<Spectrum*> spectra;
    fHistDefs.GetSpectra(spectra);
    std::vector<ReweightableSpectrum*> rwSpectra;
    fHistDefs.GetReweightableSpectra(rwSpectra);
    // Normally using a floating-point number as a key is a bad idea. But this
    // seems to work, and if it goes wrong all it breaks is this printout.
    std::map<double, int> livetimes;
    for(Spectrum* s: spectra) ++livetimes[s->Livetime()];
    for(ReweightableSpectrum* rs: rwSpectra) ++livetimes[rs->fLivetime];

    if(livetimes.size() == 1){
      const double t = livetimes.begin()->first;
      if(t > 0)
        std::cout << "Total livetime " << t << " seconds." << std::endl;
      else
        std::cout << "Unknown livetime." << std::endl;
    }
    else{
      std::cout << "Livetimes:" << std::endl;
      for(auto it: livetimes){
        std::cout << "  " << it.second << " spectra: ";
        if(it.first > 0)
          std::cout << it.first << " seconds" << std::endl;
        else
          std::cout << "unknown" << std::endl;
      }
    }
  }

  //----------------------------------------------------------------------
  void SpectrumLoader::AccumulateExposures(const caf::SRSpillProxy* spill)
  {
    const unsigned int N = fAllCuts.size();
    for(unsigned int i = 0; i < N; ++i){
      const Cut& cut = fAllCuts[i];

      // FD MC is 1 neutrino-per-event, so the livetime count will be all
      // screwed up there. Just go with 0, ie "unknown".
      double livetime = 0;
      if(spill->det != caf::kFARDET || !spill->ismc || fSource == kCosmic)
        livetime = cut.Livetime(spill);

      // Not implemented -> no good way to count
      if(livetime < 0){
        if(spill->det == caf::kFARDET && !spill->ismc){
          // People might well expect FD data to have a livetime, flag it
          // so that we can make a warning about it at the end.
          fLivetimeByCut[i] = -1;
        }
        // Otherwise the user probably doesn't care. Just don't sum it into the
        // livetime
      }
      else{
        // If a previous spill had no livetime and set the flag, and this one
        // then does have good livetime, we're in big trouble.
        assert(fLivetimeByCut[i] >= 0);

        // Livetime was fine, accumulate it
        fLivetimeByCut[i] += livetime;
      }


      double pot = cut.POT(spill);
      // Not implemented -> fall back to usual counting
      if(pot < 0) pot = spill->spillpot;

      fPOTByCut[i] += pot;
    }
  }

  //----------------------------------------------------------------------
  void SpectrumLoader::StoreExposures()
  {
    std::map<int, double> livetime;
    std::map<int, double> pot;

    for(unsigned int i = 0; i < fAllCuts.size(); ++i){
      const int id = fAllCuts[i].ID();
      if(fLivetimeByCut[i] < 0){
        fLivetimeByCut[i] = 0;
        std::cout << "WARNING: no way to compute livetime for FD data spectrum. If you want a livetime you need to be applying one of the cuts from TimingCuts.h or similar. You probably should be anyway to remove bad data near the spill ends." << std::endl;
      }
      livetime.emplace(id, fLivetimeByCut[i]);
      pot.emplace(id, fPOTByCut[i]);
    }


    for(auto& shiftdef: fHistDefs){
      for(auto& cutdef: shiftdef.second){
        const Cut& cut = cutdef.first;
        const int id = cut.ID();

        for(auto& weidef: cutdef.second){
          for(auto& vardef: weidef.second){
            for(Spectrum* s: vardef.second.spects){
              s->fPOT += pot[id];
              s->fLivetime += livetime[id];
            }

            for(ReweightableSpectrum* rw: vardef.second.rwSpects){
              rw->fPOT += pot[id];
              rw->fLivetime += livetime[id];
            }
          }
        }
      }
    }
  }

} // namespace