ana::SpectrumLoader Class Reference

Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files. More...

#include <SpectrumLoader.h>

Inheritance diagram for ana::SpectrumLoader:

Public Member Functions
	SpectrumLoader (const std::string &wildcard, DataSource src=kBeam, int max=0)
	SpectrumLoader (const std::vector< std::string > &fnames, DataSource src=kBeam, int max=0)
virtual	~SpectrumLoader ()
virtual void	Go () override
	Load all the registered spectra. More...
Public Member Functions inherited from ana::SpectrumLoaderBase
virtual	~SpectrumLoaderBase ()
virtual void	AddSpectrum (Spectrum &spect, const Var &var, const SpillCut &spillcut, const Cut &cut, const SystShifts &shift, const Var &wei=kUnweighted)
	For use by the Spectrum constructor. More...
virtual void	AddSpectrum (Spectrum &spect, const MultiVar &var, const SpillCut &spillcut, const Cut &cut, const SystShifts &shift, const Var &wei=kUnweighted)
	For use by the Spectrum constructor. More...
virtual void	AddSpectrum (Spectrum &spect, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
	For use by the Spectrum constructor. More...
virtual void	AddSpectrum (Spectrum &spect, const SpillMultiVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
	For use by the Spectrum constructor. More...
virtual void	AddReweightableSpectrum (ReweightableSpectrum &spect, const Var &var, const Cut &cut, const SystShifts &shift, const Var &wei)
	For use by the constructors of ReweightableSpectrum subclasses. More...
virtual void	AddReweightableSpectrum (ReweightableSpectrum &spect, const Var &var, const SpillCut &spillcut, const SliceCut &slicecut, const SystShifts &shift, const Var &wei)
	For use by the constructors of ReweightableSpectrum subclasses. More...
virtual bool	Gone () const
	Indicate whether or not Go has been called. More...

Static Public Member Functions
static SpectrumLoader	FromSAMProject (const std::string &proj, DataSource src=kBeam, int fileLimit=-1)
	Named constructor for SAM projects. More...

Protected Member Functions
	SpectrumLoader (DataSource src=kBeam)
	SpectrumLoader (SpectrumLoader &&)=default
SpectrumLoader &	operator= (SpectrumLoader &&)=default
	SpectrumLoader (const SpectrumLoader &)=delete
SpectrumLoader &	operator= (const SpectrumLoader &)=delete
virtual void	HandleFile (TFile *f, Progress *prog=0)
virtual void	HandleRecord (caf::SRSpillProxy *sr)
virtual void	StoreExposures ()
	Save accumulated exposures into the individual spectra. More...
Protected Member Functions inherited from ana::SpectrumLoaderBase
	SpectrumLoaderBase (DataSource src=kBeam)
	Component of other constructors. More...
	SpectrumLoaderBase (const std::string &wildcard, DataSource src=kBeam)
	Construct from a filename, wildcard, SAM definition, or SAM query. More...
	SpectrumLoaderBase (const std::vector< std::string > &fnames, DataSource src=kBeam)
	Construct from an explicit list of files. More...
	SpectrumLoaderBase (SpectrumLoaderBase &&)=default
SpectrumLoaderBase &	operator= (SpectrumLoaderBase &&)=default
	SpectrumLoaderBase (const SpectrumLoaderBase &)=delete
SpectrumLoaderBase &	operator= (const SpectrumLoaderBase &)=delete
IFileSource *	WildcardOrSAMQuery (const std::string &str) const
	Figure out if str is a wildcard or SAM query and return a source. More...
virtual void	RemoveSpectrum (Spectrum *)
virtual void	RemoveReweightableSpectrum (ReweightableSpectrum *)
int	NFiles () const
	Forwards to fFileSource. More...
TFile *	GetNextFile ()

Protected Attributes
int	max_entries
	All unique cuts contained in fHistDefs. More...
Protected Attributes inherited from ana::SpectrumLoaderBase
std::string	fWildcard
std::unique_ptr< IFileSource >	fFileSource
DataSource	fSource
bool	fGone
	Has Go() been called? Can't add more histograms after that. More...
double	fPOT
double	fPOTFromHist
	Accumulated by calls to GetNextFile. More...
int	fNReadouts
IDMap< SpillCut, IDMap < SystShifts, IDMap< Cut, IDMap< Var, IDMap < VarOrMultiVar, SpectList > > > > >	fHistDefs
	All the spectra that need to be filled. More...
IDMap< SpillCut, IDMap < SpillVar, IDMap < SpillVarOrMultiVar, SpectList > > >	fSpillHistDefs
	[spillcut][spillwei][spillvar] More...

Additional Inherited Members
Protected Types inherited from ana::SpectrumLoaderBase
typedef _VarOrMultiVar < caf::SRSliceProxy >	VarOrMultiVar
typedef _VarOrMultiVar < caf::SRSpillProxy >	SpillVarOrMultiVar

Detailed Description

Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.

The outside user should just pass a filename, wildcard, SAM dataset name, or query to the constructor. Then construct all the spectra you need. They will register with this loader. Finally, calling Go will cause all the spectra to be filled at once. After this the loader may not be used again.

Definition at line 19 of file SpectrumLoader.h.

Constructor & Destructor Documentation

ana::SpectrumLoader::SpectrumLoader	(	const std::string &	wildcard,
		DataSource	src = kBeam,
		int	max = 0
	)

Definition at line 24 of file SpectrumLoader.cxx.

    : SpectrumLoaderBase(wildcard, src), max_entries(max)
  {
  }

ana::SpectrumLoader::SpectrumLoader	(	const std::vector< std::string > &	fnames,
		DataSource	src = kBeam,
		int	max = 0
	)

Definition at line 30 of file SpectrumLoader.cxx.

    : SpectrumLoaderBase(fnames, src), max_entries(max)
  {
  }

ana::SpectrumLoader::~SpectrumLoader ( )

virtual

Definition at line 55 of file SpectrumLoader.cxx.

  {
  }

ana::SpectrumLoader::SpectrumLoader ( DataSource  src = kBeam )

protected

Definition at line 37 of file SpectrumLoader.cxx.

    : SpectrumLoaderBase(src)
  {
  }

ana::SpectrumLoader::SpectrumLoader ( SpectrumLoader &&  )

protecteddefault

ana::SpectrumLoader::SpectrumLoader ( const SpectrumLoader &  )

protecteddelete

Member Function Documentation

SpectrumLoader ana::SpectrumLoader::FromSAMProject ( const std::string &  proj, DataSource  src = kBeam, int  fileLimit = -1  )

static

Named constructor for SAM projects.

Definition at line 43 of file SpectrumLoader.cxx.

  {
    SpectrumLoader ret;
    ret.fSource = src;
    ret.fWildcard = "project "+proj;
    ret.fFileSource = std::unique_ptr<IFileSource>(new SAMProjectSource(proj, fileLimit));
    return ret;
  }

void ana::SpectrumLoader::Go ( )

overridevirtual

Load all the registered spectra.

Implements ana::SpectrumLoaderBase.

Definition at line 68 of file SpectrumLoader.cxx.

  {
    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) prog = new Progress(TString::Format("Filling %lu spectra from %d files matching '%s'", fHistDefs.TotalSize(), Nfiles, fWildcard.c_str()).Data());

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

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

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

    StoreExposures(); // also triggers the POT printout

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

void ana::SpectrumLoader::HandleFile ( TFile *  f, Progress *  prog = 0  )

protectedvirtual

Definition at line 116 of file SpectrumLoader.cxx.

  {
    assert(!f->IsZombie());

    TTree* tr = (TTree*)f->Get("recTree");
    assert(tr);

    // We try to access this field for every record. It was only added to the
    // files in late 2021, and we don't want to render all earlier files
    // unusable at a stroke. This logic can safely be removed once all extant
    // files have such a field (estimate mid-2022?)
    const bool has_husk = tr->GetLeaf("rec.hdr.husk");

    caf::SRSpillProxy sr(tr, "rec");

    //    FloatingExceptionOnNaN fpnan;

    long Nentries = tr->GetEntries();
    if (max_entries != 0 && max_entries < Nentries) Nentries = max_entries;

    for(long n = 0; n < Nentries; ++n){
      tr->LoadTree(n);

      // If there is no husk field there is no concept of husk events
      if(!has_husk) sr.hdr.husk = false;

      HandleRecord(&sr);

      if(prog) prog->SetProgress(double(n)/Nentries);
    } // end for n
  }

void ana::SpectrumLoader::HandleRecord ( caf::SRSpillProxy *  sr )

protectedvirtual

Definition at line 177 of file SpectrumLoader.cxx.

  {
    if(sr->hdr.first_in_subrun){
      fPOT += sr->hdr.pot;
      // TODO think about if this should be gated behind first_in_file. At the
      // moment I think these will be synonymous. And despite the comment on
      // hdr.pot, I think it may be file-based in practice too.
      if(sr->hdr.ismc) fNReadouts += sr->hdr.ngenevt;
    }

    if(!sr->hdr.ismc){
      const int nbnb  = sr->hdr.bnbinfo.size();
      const int nnumi = sr->hdr.numiinfo.size();
      if(nbnb > 0 && nnumi > 0){
        std::cout << "SpectrumLoader: nonzero number of both BNB (" << nbnb
                  << ") and NuMI (" << nnumi << ") triggers. I'm confused"
                  << std::endl;
        abort();
      }

      fNReadouts += nbnb + nnumi;
    }

    // This record was only kept as a receptacle for exposure information. It
    // shouldn't be included in any selected spectra.
    if(sr->hdr.husk) return;

    // Do the spill-level spectra first. Keep this very simple because we
    // intend to change it.
    for(auto& spillcutdef: fSpillHistDefs){
      if(!spillcutdef.first(sr)) continue;
      for(auto& spillweidef: spillcutdef.second){
        const double wei = spillweidef.first(sr);
        if(wei == 0) continue;
        for(auto spillvardef: spillweidef.second){
          if(spillvardef.first.IsMulti()){
            for(double val: spillvardef.first.GetMultiVar()(sr)){
              for(Spectrum* s: spillvardef.second.spects) s->Fill(val, wei);
            }
          }
          else{
            const double val = spillvardef.first.GetVar()(sr);
            for(Spectrum* s: spillvardef.second.spects) s->Fill(val, wei);
          }
        }
      }
    }


    for(auto& spillcutdef: fHistDefs){
      const SpillCut& spillcut = spillcutdef.first;

      const bool spillpass = spillcut(sr); // nomSpillCutCache.Get(spillcut, sr);
      // Cut failed, skip all the histograms that depended on it
      if(!spillpass) continue;

      for(caf::SRSliceProxy& slc: sr->slc){

        // 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, caf::SRSliceProxy> nomCutCache;
        CutVarCache<double, Var, caf::SRSliceProxy> nomWeiCache;
        CutVarCache<double, Var, caf::SRSliceProxy> nomVarCache;

        for(auto& shiftdef: spillcutdef.second){
          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.

          caf::SRProxySystController::BeginTransaction();

          bool shifted = false;

          double systWeight = 1;
          // Can special-case nominal to not pay cost of Shift()
          if(!shift.IsNominal()){
            shift.Shift(&slc, 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(&slc) : nomCutCache.Get(cut, &slc);
            // 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(&slc) : nomWeiCache.Get(weivar, &slc);

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

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

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

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

                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()(&slc);

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

                  rw->fHist->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
      } // end for slc
    } // end for spillcutdef
  }

SpectrumLoader& ana::SpectrumLoader::operator= ( SpectrumLoader &&  )

protecteddefault

SpectrumLoader& ana::SpectrumLoader::operator= ( const SpectrumLoader &  )

protecteddelete

void ana::SpectrumLoader::StoreExposures ( )

protectedvirtual

Save accumulated exposures into the individual spectra.

Definition at line 326 of file SpectrumLoader.cxx.

  {
    if(fabs(fPOT - fPOTFromHist)/std::min(fPOT, fPOTFromHist) > 0.001){
      std::cout << fPOT << " POT from hdr differs from " << fPOTFromHist << " POT from the TotalPOT histogram!" << std::endl;
      abort();
    }

    std::cout << fPOT << " POT over " << fNReadouts << " readouts" << std::endl;

    for(auto& shiftdef: fHistDefs){
      for(auto& spillcutdef: shiftdef.second){
        for(auto& cutdef: spillcutdef.second){
          for(auto& weidef: cutdef.second){
            for(auto& vardef: weidef.second){
              for(Spectrum* s: vardef.second.spects){
                s->fPOT += fPOT;
                s->fLivetime += fNReadouts;
              }
              for(ReweightableSpectrum* rw: vardef.second.rwSpects){
                rw->fPOT += fPOT;
                rw->fLivetime += fNReadouts;
              }
            }
          }
        }
      }
    }


    for(auto& spillcutdef: fSpillHistDefs){
      for(auto& spillweidef: spillcutdef.second){
        for(auto spillvardef: spillweidef.second){
          for(Spectrum* s: spillvardef.second.spects){
            s->fPOT += fPOT;
            s->fLivetime += fNReadouts;
          }
        }
      }
    }

  }

Member Data Documentation

int ana::SpectrumLoader::max_entries

protected

All unique cuts contained in fHistDefs.

Definition at line 55 of file SpectrumLoader.h.

---

The documentation for this class was generated from the following files:

• SpectrumLoader.h
• SpectrumLoader.cxx