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 "CAFAna/Core/CutApplier.h"
#include "CAFAna/Core/SystApplier.h"
#include "CAFAna/Core/WeightApplier.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)
    : SpectrumLoaderBase(wildcard)
  {
  }

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

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

  //----------------------------------------------------------------------
  SpectrumLoader SpectrumLoader::FromSAMProject(const std::string& proj,
                                                int fileLimit)
  {
    SpectrumLoader ret;
    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;

    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 %d spectra", fSink.NSinks()).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

    fExposureSource.Finalize();

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

    std::cout << "Total POT before spill cuts: " << fRunPOT << std::endl;
    fExposureSource.ReportExposures(fRunPOT);

    // Delete the tree of syst/cut/weight appliers
    for(auto it: fSinksS) delete it.second;
    for(auto it: fSinksSC) delete it.second;
    for(auto it: fSinksSCW) delete it.second;
    fSinksS.clear();
    fSinksSC.clear();
    fSinksSCW.clear();

    // Delete the actual spectrum sinks themselves
    for(IRecordSink* s: fSpectrumSinks) delete s;
    fSpectrumSinks.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
  }

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

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

    fSink.HandleRecord(sr, 1);
  }

  //----------------------------------------------------------------------
  void SpectrumLoader::AccumulateExposures(const caf::SRSpillProxy* spill)
  {
    fExposureSource.HandleSpill(spill);
  }

} // namespace