SpectrumLoaderBase.cxx

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoaderBase.h"

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

#include "CAFAna/Core/CutApplier.h"
#include "CAFAna/Core/SystApplier.h"
#include "CAFAna/Core/WeightApplier.h"

#include "StandardRecord/Proxy/SRProxy.h"

#include "ifdh.h"
#include "WebAPI.h"

#include "TFile.h"
#include "TH1.h"
#include "TTree.h"

#include <algorithm>
#include <cassert>
#include <iostream>
#include "boost/algorithm/string.hpp"

namespace ana
{
  //----------------------------------------------------------------------
  SpectrumLoaderBase::SpectrumLoaderBase()
    : fGone(false),
      fPOT(0), fRunPOT(0), fSinglePOT(0), fBeamNue(0), fSingleNue(0), fWeightedPOT(0)
  {
  }

  //----------------------------------------------------------------------
  SpectrumLoaderBase::SpectrumLoaderBase(const std::string& wildcard)
    : SpectrumLoaderBase()
  {
    fWildcard = wildcard;
    fFileSource = std::unique_ptr<IFileSource>(WildcardOrSAMQuery(wildcard));
  }

  //----------------------------------------------------------------------
  SpectrumLoaderBase::SpectrumLoaderBase(const std::vector<std::string>& fnames)
    : SpectrumLoaderBase()
  {
    fWildcard = "file list";
    fFileSource = std::unique_ptr<IFileSource>(new FileListSource(fnames));

    assert(!fnames.empty());
    std::cout << "Loading from " << fnames.size() << " files" << std::endl;
  }

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

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::SetSpillCut(const SpillCut& cut)
  {
    assert(!fGone);
    fSpillCut = std::make_unique<SpillCut>(cut);
  }

  //----------------------------------------------------------------------
  IFileSource* SpectrumLoaderBase::
  WildcardOrSAMQuery(const std::string& str) const
  {
    // stat() blows up on strings with spaces
    if(str.find(' ') == std::string::npos){
      WildcardSource* ret = new WildcardSource(str, Stride(), Offset(), Limit());
      if(ret->NFiles() > 0) return ret;
      delete ret;
    }

    // Maybe this the name of a SAM project?
    {
      IFDHSilent silent; // the usual case is for this to fail
      ifdh_ns::ifdh i;

      // findProject always gives back an address just by gluing bits together.
      // (a tad annoying, because it _does_ go and look for the project,
      //  and even would print its 'didn't-find-this-project' error out to stderr
      //  if not for the IFDHSilent, but without scraping its stderr
      //  there's no way to know whether the project is there or not--
      //  you still get the URL.)
      // however, the WebAPI call looking for the /status will return a 404 if
      // the project doesn't exist.  (suggested by Robert I. in INC000000925362.)
      try
      {
        ifdh_util_ns::WebAPI webapi(i.findProject(str, getenv("SAM_STATION")) + "/status");
        return new SAMProjectSource(str);
      }
      catch (ifdh_util_ns::WebAPIException &e)
      {
        ;
      }
    }

    // Maybe this is a SAM dataset or query?
    return new SAMQuerySource(str, Stride(), Offset(), Limit());
  }

  //----------------------------------------------------------------------
  WeightApplier* SpectrumLoaderBase::GetSinkSCW(const SystShifts& shift,
                                                const Cut& cut,
                                                const Var& wei)
  {
    const int sid = shift.ID();
    const int cid = cut.ID();
    const int wid = wei.ID();

    const std::pair<int, int> sc(sid, cid);
    const std::tuple<int, int, int> scw(sid, cid, wid);

    if(fSinksS.find(sid) == fSinksS.end()){
      SystApplier* sa = new SystApplier(shift);
      fSinksS[sid] = sa;
      fSink.AddSink(sa);
    }
    if(fSinksSC.find(sc) == fSinksSC.end()){
      CutApplier* ca = new CutApplier(cut);
      fSinksS[sid]->AddSink(ca);
      fSinksSC[sc] = ca;
    }
    if(fSinksSCW.find(scw) == fSinksSCW.end()){
      WeightApplier* wa = new WeightApplier(wei);
      fSinksSC[sc]->AddSink(wa);
      fSinksSCW[scw] = wa;
    }

    return fSinksSCW[scw];
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddSpectrum(Spectrum& spect,
                                       const Var& var,
                                       const Cut& cut,
                                       const SystShifts& shift,
                                       const Var& wei)
  {
    if(fGone){
      std::cerr << "Error: can't add Spectra after the call to Go()" << std::endl;
      abort();
    }

    SpectrumSink* sink = new SpectrumSink(var, &spect);
    fSpectrumSinks.insert(sink);
    GetSinkSCW(shift, cut, wei)->AddSink(sink);
    fExposureSource.AddSink(cut, sink);
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddSpectrum(Spectrum& spect,
                                       const MultiVar& var,
                                       const Cut& cut,
                                       const SystShifts& shift,
                                       const Var& wei)
  {
    if(fGone){
      std::cerr << "Error: can't add Spectra after the call to Go()" << std::endl;
      abort();
    }

    MultiVarSpectrumSink* sink = new MultiVarSpectrumSink(var, &spect);
    fSpectrumSinks.insert(sink);
    GetSinkSCW(shift, cut, wei)->AddSink(sink);
    fExposureSource.AddSink(cut, sink);
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddSpectrum(Spectrum& spect,
                                       const NuTruthVar& var,
                                       const NuTruthCut& cut,
                                       const SystShifts& shift,
                                       const NuTruthVar& wei)
  {
    // TODO - implement with the same source/sink deal as regular spectra
    fNuHistDefs.emplace_back(spect, var, cut, shift, wei);
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddSpectrum(Spectrum& spect,
                                       const NuTruthMultiVar& multivar,
                                       const NuTruthCut& cut,
                                       const SystShifts& shift,
                                       const NuTruthVar& wei)
  {
    // TODO - implement with the same source/sink deal as regular spectra
    fNuHistDefs.emplace_back(spect, multivar, cut, shift, wei);
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddReweightableSpectrum(ReweightableSpectrum& spect,
                                                   const Var& xvar,
                                                   const Var& yvar,
                                                   const Cut& cut,
                                                   const SystShifts& shift,
                                                   const Var& wei)
  {
    if(fGone){
      std::cerr << "Error: can't add Spectra after the call to Go()" << std::endl;
      abort();
    }

    ReweightableSpectrumSink* sink = new ReweightableSpectrumSink(xvar, yvar, &spect);
    fSpectrumSinks.insert(sink);
    GetSinkSCW(shift, cut, wei)->AddSink(sink);
    fExposureSource.AddSink(cut, sink);
  }

  //----------------------------------------------------------------------
  void SpectrumLoaderBase::AddSpillHistogram(TH1* h,
                                             const SpillVar& var,
                                             const SpillCut& cut,
                                             const SpillVar& wei)
  {
    // TODO - implement with the same source/sink deal as regular spectra
    fSpillHistDefs.emplace_back(h, var, cut, wei);
  }

  //----------------------------------------------------------------------
  int SpectrumLoaderBase::NFiles() const
  {
    return fFileSource->NFiles();
  }

  //----------------------------------------------------------------------
  TFile* SpectrumLoaderBase::GetNextFile()
  {
    TFile* f = fFileSource->GetNextFile();
    if(!f) return 0; // out of files

    if(f->IsZombie()){
      std::cout << "Bad file (zombie): " << f->GetName() << std::endl;
      abort();
    }

    // Crude way to count POT, at a file level without referring to spill cuts
    TH1* hPOT = (TH1*)f->Get("TotalPOT");

    if(!hPOT){
      std::cout << "Bad file (no POT tree): " << f->GetName() << std::endl;
      abort();
    }

    fRunPOT += hPOT->Integral();

    TKey *key = f->FindKey("TotalSinglePOT");

    // Initialise and only fill these histos if they exist
    // We know they exist for nue overlay files
    TH1* hSinglePOT           = 0;
    TH1* hTotalTrueSingleNue  = 0;
    TH1* hTotalTrueNonswapNue = 0;
    if(key != 0){
      hSinglePOT = (TH1*)f->Get("TotalSinglePOT");
      assert(hSinglePOT);
      fSinglePOT = hSinglePOT->Integral();
      
      hTotalTrueNonswapNue = (TH1*)f->Get("TotalTrueNonswapNue");
      assert(hTotalTrueNonswapNue);
      fBeamNue   = hTotalTrueNonswapNue->Integral();

      hTotalTrueSingleNue = (TH1*)f->Get("TotalTrueSingleNue");
      assert(hTotalTrueSingleNue);
      fSingleNue = hTotalTrueSingleNue->Integral();
    }
    else{
      delete hSinglePOT;
      delete hTotalTrueSingleNue;
      delete hTotalTrueNonswapNue;
    }

    TDirectory* spillDir = 0;
    TTree* spillTree = 0;

    TObject* obj = f->Get("spillTree");
    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")){
      spillTree = (TTree*)obj;
    }
    else{
      spillDir = (TDirectory*)obj;
      spillTree = (TTree*)spillDir->Get("spill");
      assert(spillTree);
    }

    const caf::CAFType type = caf::GetCAFType(spillDir, spillTree);

    long n = 0;
    caf::SRSpillProxy spill(spillDir, spillTree, "spill", n, 0);

    //    FloatingExceptionOnNaN fpnan;

    // With multiple files exposure counting alternates with histogram filling,
    // but with one file they deserve to be reported as two separate steps.
    Progress* prog = 0;
    if(fFileSource->NFiles() == 1)
      prog = new Progress(TString::Format("Counting exposure from 1 file matching '%s'", fWildcard.c_str()).Data());

    const int Nentries = spillTree->GetEntries();

    double totFilePOT = 0;

    for(n = 0; n < Nentries; ++n){
      if(type != caf::kFlatMultiTree) spillTree->LoadTree(n); // for all single-tree modes

      if(!fSpillCut || (*fSpillCut)(&spill)){
        AccumulateExposures(&spill);

        totFilePOT += spill.spillpot;

        for(const SpillHistDef& hd: fSpillHistDefs){
          // No cut, or cut passes
          if(hd.cut(&spill)){
            const double wei = hd.wei(&spill);
            if(wei){
              const double val = hd.var(&spill);
              hd.h->Fill(val, wei);
            }
          } // end if(cut passes)
        } // end for hd
      }

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

    // The POT we want to use is either the singles POT in case of overlay files
    // i.e the singlepot histo exists, or the standard after spill cuts
    // Note: must be sure to run spill cuts for the overlays otherwise this gets
    // reset to the fRunPOT
    if(key != 0) fPOT += fSinglePOT;
    else         fPOT += totFilePOT;
    
    delete prog;
    prog = 0;

    if(!fNuHistDefs.empty()){
      TDirectory* nuDir = 0;
      TTree* nuTree = 0;

      TObject* obj = f->Get("nuTree");
      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")){
        nuTree = (TTree*)obj;
      }
      else{
        nuDir = (TDirectory*)obj;
        nuTree = (TTree*)spillDir->Get("nu");
        assert(nuTree);
      }

      const caf::CAFType type = caf::GetCAFType(nuDir, nuTree);

      const int NNuEntries = nuTree->GetEntries();

      long n = 0;
      caf::SRNeutrinoProxy nu(nuDir, nuTree, (type == caf::kNested) ? "mc.nu" : "nu", n, 0);

      if(fFileSource->NFiles() == 1)
        prog = new Progress(TString::Format("Filling from nuTree from 1 file matching '%s'", fWildcard.c_str()).Data());


      for(n = 0; n < NNuEntries; ++n){
        if(type != caf::kFlatMultiTree) nuTree->LoadTree(n); // for all single-tree modes

        for(const NuHistDef& hd: fNuHistDefs){
          caf::SRProxySystController::BeginTransaction();

          double systWeight = 1;
          // Can special-case nominal to not pay cost of Shift()
          if(!hd.shift.IsNominal()){
            hd.shift.Shift(&nu, systWeight);
          }

          if(hd.cut(&nu)){
            const double wei = hd.wei(&nu)*systWeight;
            if(wei){
              if(hd.var){
                const double val = (*hd.var)(&nu);
                hd.spect.Fill(val, wei);
              }
              else{
                for(double val: (*hd.multivar)(&nu)) hd.spect.Fill(val, wei);
              }
            }
          }

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

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

      for(const NuHistDef& hd: fNuHistDefs) hd.spect.fPOT += totFilePOT;

      delete prog;
    } // end if NuHistDefs

    return f;
  }

  //---------------------------------------------------------------------
  double SpectrumLoaderBase::GetWeightedPOT()
  {
    fWeightedPOT = ((fBeamNue/(fBeamNue+fSingleNue))*fPOT) + ((fSingleNue/(fBeamNue+fSingleNue))*fSinglePOT);
    std::cout << " beam nue:  " << fBeamNue << "  single nue:  " << fSingleNue << "  fPOT:  " << fPOT << "  fSinglePOT:  " << fSinglePOT << std::endl;
    return fWeightedPOT;
  }

  //----------------------------------------------------------------------
  void NullLoader::Go()
  {
  }

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

} // namespace