SBNSoftware · cjbacchus · Mar 30, 2021 · Mar 30, 2021 · Apr 13, 2021 · Apr 13, 2021
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -36,7 +36,7 @@ cet_set_compiler_flags(DIAGS CAUTIOUS
   WERROR
   NO_UNDEFINED
   ALLOW_DEPRECATIONS
-  EXTRA_FLAGS -pedantic -Wno-unused-local-typedefs -Wno-undefined-var-template -I $ENV{CLANG_FQ_DIR}/include/c++/v1
+  EXTRA_FLAGS -pedantic -Wno-unused-local-typedefs -Wno-undefined-var-template -I $ENV{CLANG_FQ_DIR}/include/c++/v1 -DTBB_INTERFACE_NEW
 )
 
 cet_report_compiler_flags()
@@ -49,6 +49,8 @@ find_ups_root(v6_00_00)
 find_ups_product(ifdhc)
 find_ups_product( srproxy )
 find_ups_product( osclib )
+find_ups_product( cafanacore )
+find_ups_product( stan_math )
 find_ups_product( pygccxml )
 
 #ifdhc

diff --git a/sbnana/CAFAna/Analysis/CMakeLists.txt b/sbnana/CAFAna/Analysis/CMakeLists.txt
@@ -1,6 +1,3 @@
-# Allow "StandardRecord/" to find the right directory
-include_directories(..)
-
 cet_make(LIBRARY_NAME CAFAnaAnalysis
          LIBRARIES CAFAnaCore CAFAnaExperiment
                    ${ROOT_MINUIT2}

diff --git a/sbnana/CAFAna/Analysis/Calcs.cxx b/sbnana/CAFAna/Analysis/Calcs.cxx
@@ -25,7 +25,7 @@ namespace ana
     calc->SetTh23(M_PI/4);
 
     // NH from PDG2014 + 2015 update
-    //http://pdg.lbl.gov/2015/tables/rpp2015-sum-leptons.pdf 
+    //http://pdg.lbl.gov/2015/tables/rpp2015-sum-leptons.pdf
     calc->SetDmsq32(2.44e-3); // NB: this is normal hierarchy
 
     // Reactor average from PDG2014 + 2015 update
@@ -46,11 +46,11 @@ namespace ana
   //----------------------------------------------------------------------
   void ResetOscCalcToDefaultIH(osc::IOscCalcAdjustable* calc)
   {
-    //Share most defaults 
+    //Share most defaults
     ResetOscCalcToDefault(calc);
-    // IH from PDG2014 + 2015 update 
+    // IH from PDG2014 + 2015 update
     // http://pdg.lbl.gov/2015/tables/rpp2015-sum-leptons.pdf
-    calc->SetDmsq32(-2.49e-3); 
+    calc->SetDmsq32(-2.49e-3);
 
   }
 

diff --git a/sbnana/CAFAna/Analysis/CalcsNuFit_cdr.cxx b/sbnana/CAFAna/Analysis/CalcsNuFit_cdr.cxx
@@ -49,7 +49,7 @@ namespace ana
 
     // Uniform throws within +/-3 sigma
     if(hie > 0){
-      ret->SetDmsq32(gRandom->Uniform(kNuFitDmsq32CVNH-3*kNuFitDmsq32ErrNH, 
+      ret->SetDmsq32(gRandom->Uniform(kNuFitDmsq32CVNH-3*kNuFitDmsq32ErrNH,
 				      kNuFitDmsq32CVNH+3*kNuFitDmsq32ErrNH));
       ret->SetTh23(gRandom->Uniform(kNuFitTh23CVNH-3*kNuFitTh23ErrNH,
 				    kNuFitTh23CVNH+3*kNuFitTh23ErrNH));
@@ -64,7 +64,7 @@ namespace ana
                                     kNuFitTh13CVIH+3*kNuFitTh13ErrIH));
     }
     ret->SetdCP(gRandom->Uniform(-1*TMath::Pi(), TMath::Pi()));
-		  
+
     return ret;
   }
 
@@ -145,9 +145,9 @@ namespace ana
     fTh23Err = (hietrue > 0) ? 0.058*fTh23 : 0.048*fTh23;
 
     fRhoErr = 0.02*fRho;
-    
+
   }
-    
+
   double Penalizer_GlbLikeCDR::ChiSq(osc::IOscCalcAdjustable* calc,
 				  const SystShifts& /*syst*/) const {
 
@@ -156,7 +156,7 @@ namespace ana
 
     double ret =
       util::sqr((calc->GetDmsq21() - fDmsq21)/fDmsq21Err) +
-      util::sqr((calc->GetTh12() - fTh12)/fTh12Err) + 
+      util::sqr((calc->GetTh12() - fTh12)/fTh12Err) +
       util::sqr((calc->GetRho() - fRho)/fRhoErr);
 
     // if fWeakOnly is set, only apply a constraint to the parameter we can only weakly constrain in DUNE

diff --git a/sbnana/CAFAna/Analysis/CutOptimizer.cxx b/sbnana/CAFAna/Analysis/CutOptimizer.cxx
@@ -1,13 +1,15 @@
 #include "sbnana/CAFAna/Analysis/CutOptimizer.h"
 
+#include "sbnana/CAFAna/Core/SpectrumLoader.h"
+
 #include "TH1.h"
 
 #include <iostream>
 
 namespace ana
 {
   // --------------------------------------------------------------------------
-  void MakeNMinusOneSpectra(SpectrumLoader& loader,
+  void MakeNMinusOneSpectra(ISliceSource& src,
                             const Cut& sigcut,
                             const Cut& presel,
                             const std::vector<HistAxis>& axes,
@@ -28,8 +30,8 @@ namespace ana
         nminusone = nminusone && axes[j].GetVars()[0] > cut_pos[j];
       } // end for j
 
-      sigs.emplace_back(loader, axes[i], nminusone && sigcut);
-      bkgs.emplace_back(loader, axes[i], nminusone && !sigcut);
+      sigs.emplace_back(src[nminusone][ sigcut], axes[i]);
+      bkgs.emplace_back(src[nminusone][!sigcut], axes[i]);
     } // end for i
   }
 
@@ -64,7 +66,7 @@ namespace ana
   {
     SpectrumLoader loader(wildcard);
     std::vector<Spectrum> sigs, bkgs;
-    MakeNMinusOneSpectra(loader, sigcut, presel, axes, cut_pos, sigs, bkgs);
+    MakeNMinusOneSpectra(loader.Slices(), sigcut, presel, axes, cut_pos, sigs, bkgs);
     loader.Go();
 
     double best_fom = 0;

diff --git a/sbnana/CAFAna/Analysis/CutOptimizer.h b/sbnana/CAFAna/Analysis/CutOptimizer.h
@@ -2,8 +2,9 @@
 
 #include "sbnana/CAFAna/Core/HistAxis.h"
 #include "sbnana/CAFAna/Core/Cut.h"
-#include "sbnana/CAFAna/Core/Spectrum.h"
-#include "sbnana/CAFAna/Core/SpectrumLoader.h"
+#include "cafanacore/Spectrum.h"
+
+#include "sbnana/CAFAna/Core/IRecordSource.h"
 
 #include <string>
 
@@ -14,14 +15,14 @@ namespace ana
   /// This function does *not* call Go() on the loader, allowing you to combine
   /// it with other spectra-filling if desired.
   ///
-  /// \param loader  SpectrumLoader to associate spectra with
+  /// \param src     Source to associate spectra with
   /// \param sigcut  Definition of signal (inverse defines background)
   /// \param presel  Cuts that will always be applied
   /// \param axes    Variables that will be cut on
   /// \param cut_pos Cut position for each variable
   /// \param[out] sigs Resulting signal spectra
   /// \param[out] bkgs Resulting background spectra
-  void MakeNMinusOneSpectra(SpectrumLoader& loader,
+  void MakeNMinusOneSpectra(ISliceSource& src,
                             const Cut& sigcut,
                             const Cut& presel,
                             const std::vector<HistAxis>& axes,

diff --git a/sbnana/CAFAna/Analysis/ExpInfo.h b/sbnana/CAFAna/Analysis/ExpInfo.h
@@ -7,7 +7,7 @@ namespace ana
   const double kPOTuBoone = 1.3e21;
   const double kBaselineSBND = 0.11;
   const double kBaselineMicroBoone = 0.47;
-  const double kBaselineIcarus = 0.6; 
+  const double kBaselineIcarus = 0.6;
   const int kSBND = 0;
   const int kMicroBoone = 1;
   const int kICARUS = 2;

diff --git a/sbnana/CAFAna/Analysis/Fit.cxx b/sbnana/CAFAna/Analysis/Fit.cxx
@@ -1,6 +1,6 @@
 #include "sbnana/CAFAna/Analysis/Fit.h"
 
-#include "sbnana/CAFAna/Core/Progress.h"
+#include "cafanacore/Progress.h"
 #include "sbnana/CAFAna/Core/IFitVar.h"
 #include "sbnana/CAFAna/Core/Utilities.h"
 #include "sbnana/CAFAna/Experiment/IExperiment.h"
@@ -166,7 +166,7 @@ namespace ana
 	mnMin->GetMinosError(i, errLow, errHigh);
 	std::cout << i << "/" << mnMin->NDim() << " " << fParamNames[i] << ": " << errLow << ", +" << errHigh << " (" << mnMin->Errors()[i] << ")" << std::endl;
 	fTempMinosErrors.push_back(std::make_pair(errLow,errHigh));
-      }      
+      }
     }
 
     return mnMin;
@@ -191,8 +191,8 @@ namespace ana
       // Need to deal with parameters that are not fit values!
       SystShifts shift = pt.shift;
 
-      // Need to copy over syst values into the seed for this fit 
-      // that were in the input syst shifts, but are not being fit for, 
+      // Need to copy over syst values into the seed for this fit
+      // that were in the input syst shifts, but are not being fit for,
       // and therefore not part of the seedPts list
       for(auto s: bestSysts.ActiveSysts()) {
         auto fit_systs = shift.ActiveSysts();

diff --git a/sbnana/CAFAna/Analysis/Fit.h b/sbnana/CAFAna/Analysis/Fit.h
@@ -1,11 +1,14 @@
 #pragma once
 
-#include "sbnana/CAFAna/Core/Spectrum.h"
+#include "cafanacore/Spectrum.h"
 #include "sbnana/CAFAna/Prediction/IPrediction.h"
 #include "sbnana/CAFAna/Core/SystShifts.h"
+#include "sbnana/CAFAna/Core/Utilities.h"
 
 #include "Math/Minimizer.h"
 
+#include "TMatrixDSym.h"
+
 #include <memory>
 
 class TGraph;

diff --git a/sbnana/CAFAna/Analysis/FitAxis.h b/sbnana/CAFAna/Analysis/FitAxis.h
@@ -13,7 +13,7 @@ namespace ana
     FitAxis(const IFitVar* var,
             int nbins, double min, double max,
             bool islog = false);
-    
+
     const IFitVar* var;
     int nbins;
     double min;

diff --git a/sbnana/CAFAna/Analysis/GradientDescent.cxx b/sbnana/CAFAna/Analysis/GradientDescent.cxx
@@ -119,7 +119,7 @@ namespace ana
     const unsigned int N = pt.size();
     ROOT::Minuit2::MnAlgebraicVector vec(N);
     for(unsigned int i = 0; i < N; ++i) vec(i) = pt[i];
-    ROOT::Minuit2::MinimumParameters params(vec, chi);    
+    ROOT::Minuit2::MinimumParameters params(vec, chi);
     ROOT::Minuit2::MinimumState state(params, 0, ncalls);
     ROOT::Minuit2::MnUserTransformation trans(pt, std::vector<double>(N));
     ROOT::Minuit2::MinimumSeed seed(state, trans);

diff --git a/sbnana/CAFAna/Analysis/MedianSurface.cxx b/sbnana/CAFAna/Analysis/MedianSurface.cxx
@@ -120,22 +120,32 @@ namespace ana
   }
 
   //----------------------------------------------------------------------
-  void MedianSurface::SaveTo(TDirectory* dir) const
+  void MedianSurface::SaveTo(TDirectory* dir, const std::string& name) const
   {
     TDirectory* tmp = gDirectory;
+
+    dir = dir->mkdir(name.c_str()); // switch to subdir
     dir->cd();
+
     TObjString("MedianSurface").Write("type");
 
     for(unsigned int i = 0; i < fThrows.size(); ++i){
-      fThrows[i].SaveTo(dir->mkdir(TString::Format("surf%d", i)));
+      fThrows[i].SaveTo(dir, TString::Format("surf%d", i).Data());
     }
 
+    dir->Write();
+    delete dir;
+
     tmp->cd();
   }
 
   //----------------------------------------------------------------------
-  std::unique_ptr<MedianSurface> MedianSurface::LoadFrom(TDirectory* dir)
+  std::unique_ptr<MedianSurface> MedianSurface::LoadFrom(TDirectory* dir,
+                                                         const std::string& name)
   {
+    dir = dir->GetDirectory(name.c_str()); // switch to subdir
+    assert(dir);
+
     DontAddDirectory guard;
 
     TObjString* tag = (TObjString*)dir->Get("type");
@@ -144,9 +154,11 @@ namespace ana
 
     std::vector<Surface> surfs;
     for(unsigned int i = 0; ; ++i){
-      TDirectory* surfdir = dir->GetDirectory(TString::Format("surf%d", i));
-      if(!surfdir) break; // we got all of them
-      surfs.push_back(*ana::LoadFrom<Surface>(surfdir));
+      const std::string subname = TString::Format("surf%d", i).Data();
+      TDirectory* subdir = dir->GetDirectory(subname.c_str());
+      if(!subdir) break; // we got all of them
+      delete subdir;
+      surfs.push_back(*ana::LoadFrom<Surface>(dir, subname));
     }
 
     return std::make_unique<MedianSurface>(surfs);

diff --git a/sbnana/CAFAna/Analysis/MedianSurface.h b/sbnana/CAFAna/Analysis/MedianSurface.h
@@ -14,8 +14,9 @@ namespace ana
     void DrawEnsemble(TH2* fc, Color_t color = kGray);
     void DrawBand(TH2* fc);
 
-    void SaveTo(TDirectory * dir) const;
-    static std::unique_ptr<MedianSurface> LoadFrom(TDirectory * dir);
+    void SaveTo(TDirectory* dir, const std::string& name) const;
+    static std::unique_ptr<MedianSurface> LoadFrom(TDirectory* dir,
+                                                   const std::string& name);
   protected:
     std::vector<Surface> fThrows;