diff --git a/src/samples/MCTools/NuisToCAFTruth.cxx b/src/samples/MCTools/NuisToCAFTruth.cxx
index 5ba6935..2377685 100644
--- a/src/samples/MCTools/NuisToCAFTruth.cxx
+++ b/src/samples/MCTools/NuisToCAFTruth.cxx
@@ -1,253 +1,253 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 //********************************************************************
 
 #include "samples/SimpleMCStudy.hxx"
 
 #include "utility/EventTopologyUtility.hxx"
 #include "utility/FullEventUtility.hxx"
 #include "utility/KinematicUtility.hxx"
 #include "utility/PDGCodeUtility.hxx"
 #include "utility/ROOTUtility.hxx"
 
 #include "persistency/ROOTOutput.hxx"
 
 using namespace nuis::event;
 using namespace nuis::utility;
 
 class NuisToCAFTruth : public SimpleMCStudy {
 
   TreeFile outputTree;
 
   // int mode;
 
   int isFD;
   int isFHC;
   int isCC;
 
   int nuPDG;
   int nuPDGunosc;
   int LepPDG;
 
   double NuMomX;
   double NuMomY;
   double NuMomZ;
   double Ev;
   double LepMomX;
   double LepMomY;
   double LepMomZ;
   double LepE;
   double LepNuAngle;
 
   double Q2;
   double W;
   double X;
   double Y;
 
   int nP;
   int nN;
   int nipip;
   int nipim;
   int nipi0;
   int nikp;
   int nikm;
   int nik0;
   int niem;
   int niother;
   int nNucleus;
   int nUNKNOWN;
 
   double eP;
   double eN;
   double ePip;
   double ePim;
   double ePi0;
   double eOther;
 
   double xsecWeight;
 
   bool isOscSwap;
 
 public:
   NuisToCAFTruth() { ReadGlobalConfigDefaults(); }
 
   void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
 
-    outputTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(),
-                                   "cafTruthTree");
-
     isFD = instance_sample_configuration.get<bool>("is_fd", false);
     isFHC = instance_sample_configuration.get<bool>("is_numode", true);
     isOscSwap = instance_sample_configuration.get<bool>("is_nueswap", false);
 
     // Perform any per-sample configuration in the base class
     SimpleMCStudy::Initialize(instance_sample_configuration);
 
+    outputTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(),
+                                   "cafTruthTree", write_directory);
+
     // outputTree->Branch("mode", &mode, "mode/I");
     outputTree->Branch("isFD", &isFD, "isFD/I");
     outputTree->Branch("isFHC", &isFHC, "isFHC/I");
     outputTree->Branch("isCC", &isCC, "isCC/I");
     outputTree->Branch("nuPDG", &nuPDG, "nuPDG/I");
     outputTree->Branch("nuPDGunosc", &nuPDGunosc, "nuPDGunosc/I");
     outputTree->Branch("LepPDG", &LepPDG, "LepPDG/I");
     outputTree->Branch("NuMomX", &NuMomX, "NuMomX/D");
     outputTree->Branch("NuMomY", &NuMomY, "NuMomY/D");
     outputTree->Branch("NuMomZ", &NuMomZ, "NuMomZ/D");
     outputTree->Branch("Ev", &Ev, "Ev/D");
     outputTree->Branch("LepMomX", &LepMomX, "LepMomX/D");
     outputTree->Branch("LepMomY", &LepMomY, "LepMomY/D");
     outputTree->Branch("LepMomZ", &LepMomZ, "LepMomZ/D");
     outputTree->Branch("LepE", &LepE, "LepE/D");
     outputTree->Branch("LepNuAngle", &LepNuAngle, "LepNuAngle/D");
     outputTree->Branch("Q2", &Q2, "Q2/D");
     outputTree->Branch("W", &W, "W/D");
     outputTree->Branch("X", &X, "X/D");
     outputTree->Branch("Y", &Y, "Y/D");
     outputTree->Branch("nP", &nP, "nP/I");
     outputTree->Branch("nN", &nN, "nN/I");
     outputTree->Branch("nipip", &nipip, "nipip/I");
     outputTree->Branch("nipim", &nipim, "nipim/I");
     outputTree->Branch("nipi0", &nipi0, "nipi0/I");
     outputTree->Branch("nikp", &nikp, "nikp/I");
     outputTree->Branch("nikm", &nikm, "nikm/I");
     outputTree->Branch("nik0", &nik0, "nik0/I");
     outputTree->Branch("niem", &niem, "niem/I");
     outputTree->Branch("niother", &niother, "niother/I");
     outputTree->Branch("nNucleus", &nNucleus, "nNucleus/I");
     outputTree->Branch("nUNKNOWN", &nUNKNOWN, "nUNKNOWN/I");
     outputTree->Branch("eP", &eP, "eP/D");
     outputTree->Branch("eN", &eN, "eN/D");
     outputTree->Branch("ePip", &ePip, "ePip/D");
     outputTree->Branch("ePim", &ePim, "ePim/D");
     outputTree->Branch("ePi0", &ePi0, "ePi0/D");
     outputTree->Branch("eOther", &eOther, "eOther/D");
     outputTree->Branch("xsecWeight", &xsecWeight, "xsecWeight/D");
 
     //! Here you do whatever you want to do on a per-event basis.
     //! This method will be run for every event given by the input handler,
     //! the weight allows differential xsecs to be plotted easily and also
     //! contains any requested reweightable parameter variations.
     ProcessEventFunction = [&](nuis::event::FullEvent const &ev,
                                double weight) -> void {
       xsecWeight = weight;
 
       nP = 0;
       nN = 0;
       nipip = 0;
       nipim = 0;
       nipi0 = 0;
       nikp = 0;
       nikm = 0;
       nik0 = 0;
       niem = 0;
       niother = 0;
       nNucleus = 0;
       nUNKNOWN = 0;
       eP = 0;
       eN = 0;
       ePip = 0;
       ePim = 0;
       ePi0 = 0;
       eOther = 0;
 
       Particle ISLep = GetHMISNeutralLepton(ev);
       Particle FSLep = GetHMFSLepton(ev);
 
       isCC = IsChargedLepton(FSLep.pdg);
 
       nuPDG = ISLep.pdg;
       nuPDGunosc = isOscSwap
                        ? (ISLep.pdg < 0 ? pdgcodes::kNuMuBar : pdgcodes::kNuMu)
                        : ISLep.pdg;
       LepPDG = FSLep.pdg;
 
       NuMomX = ISLep.P4.X();
       NuMomY = ISLep.P4.Y();
       NuMomZ = ISLep.P4.Z();
       Ev = ISLep.P4.E();
 
       LepMomX = FSLep.P4.X();
       LepMomY = FSLep.P4.Y();
       LepMomZ = FSLep.P4.Z();
       LepE = FSLep.P4.E();
 
       LepNuAngle = FSLep.Dir().Dot(ISLep.Dir());
 
       TLorentzVector FourMomTransfer = GetEnergyMomentumTransfer(ev);
 
       Q2 = -FourMomTransfer.Mag2();
       W = sqrt(-Q2 +
                2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E() +
                pdgmasses::kNucleonAverageMass_MeV *
                    pdgmasses::kNucleonAverageMass_MeV);
       X = Q2 / (2 * pdgmasses::kNucleonAverageMass_MeV * FourMomTransfer.E());
       Y = 1 - LepE / Ev;
 
       for (Particle const &p : GetNuclearLeavingParticles(ev)) {
         if (p == FSLep) {
           continue;
         }
 
         if (IsProton(p.pdg)) {
           nP++;
           eP += p.KE();
         } else if (IsNeutron(p.pdg)) {
           nN++;
           eN += p.KE();
         } else if (IsPositivePion(p.pdg)) {
           nipip++;
           ePip += p.KE();
         } else if (IsNegativePion(p.pdg)) {
           nipim++;
           ePim += p.KE();
         } else if (IsNeutralPion(p.pdg)) {
           nipi0++;
           ePi0 += p.KE();
         } else if (IsPositiveKaon(p.pdg)) {
           nikp++;
           eOther += p.KE();
         } else if (IsNegativeKaon(p.pdg)) {
           nikm++;
           eOther += p.KE();
         } else if (IsNeutralKaon(p.pdg)) {
           nik0++;
           eOther += p.KE();
         } else if (IsGamma(p.pdg)) {
           niem++;
           eOther += p.KE();
         } else if (IsNuclearPDG(p.pdg)) {
           nNucleus++;
         } else {
           niother++;
           eOther += p.KE();
         }
       }
 
       outputTree->Fill();
     };
   }
 
   std::string Name() { return "NuisToCAFTruth"; }
 
   void Write() {}
 };
 
 //! These declarations allow your class to be loaded dynamically by NUISANCE
 DECLARE_PLUGIN(IEventProcessor, NuisToCAFTruth);
diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx
index 60551e8..75b0d47 100644
--- a/src/utility/ROOTUtility.hxx
+++ b/src/utility/ROOTUtility.hxx
@@ -1,155 +1,164 @@
 // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
 
 /*******************************************************************************
  *    This file is part of NUISANCE.
  *
  *    NUISANCE is free software: you can redistribute it and/or modify
  *    it under the terms of the GNU General Public License as published by
  *    the Free Software Foundation, either version 3 of the License, or
  *    (at your option) any later version.
  *
  *    NUISANCE is distributed in the hope that it will be useful,
  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *    GNU General Public License for more details.
  *
  *    You should have received a copy of the GNU General Public License
  *    along with NUISANCE.  If not, see <http://www.gnu.org/licenses/>.
  *******************************************************************************/
 
 #ifndef UTILITY_ROOTUTILITY_HXX_SEEN
 #define UTILITY_ROOTUTILITY_HXX_SEEN
 
 #include "exception/exception.hxx"
 
 #include "TFile.h"
 #include "TTree.h"
 
 #include <iomanip>
 #include <iostream>
 #include <memory>
 #include <string>
 
 namespace nuis {
 namespace utility {
 
 NEW_NUIS_EXCEPT(failed_to_open_TFile);
 NEW_NUIS_EXCEPT(failed_to_get_TTree);
 
 inline void CloseTFile(TFile *f = nullptr) {
   if (f) {
     std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f
               << std::endl;
     f->Close();
   }
   delete f;
 }
 
 using TFile_ptr = std::unique_ptr<TFile, decltype(&CloseTFile)>;
 inline TFile_ptr make_unique_TFile(TFile *f) {
   return TFile_ptr(f, &CloseTFile);
 }
 
 inline TFile_ptr CheckOpenTFile(std::string const &fname,
                                 char const *opts = "") {
 
   TDirectory *ogDir = gDirectory;
 
   TFile *inpF = new TFile(fname.c_str(), opts);
   if (!inpF || !inpF->IsOpen()) {
     throw failed_to_open_TFile()
         << "[ERROR]: Couldn't open input file: " << std::quoted(fname);
   }
 
   if (ogDir) {
     ogDir->cd();
   }
 
   return make_unique_TFile(inpF);
 }
 
 struct TreeFile {
   TFile_ptr file;
   TTree *tree;
   bool file_owned;
   TreeFile()
       : file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {}
   TreeFile(TreeFile const &other) = delete;
   TreeFile &operator=(TreeFile &&other) {
     file = std::move(other.file);
     tree = other.tree;
     file_owned = other.file_owned;
 
     other.file = nullptr;
     other.tree = nullptr;
     other.file_owned = false;
 
     return *this;
   }
 
   TreeFile(TreeFile &&other)
       : file(std::move(other.file)), tree(other.tree),
         file_owned(other.file_owned) {
     other.file = nullptr;
     other.tree = nullptr;
     other.file_owned = false;
   }
 
   TTree *operator->() { return tree; }
 
   ~TreeFile() {
     if (!file_owned) {
       file.release();
     }
   }
 };
 
 inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) {
   TreeFile tf;
   tf.file = make_unique_TFile(file);
   tf.tree = dynamic_cast<TTree *>(tf.file->Get(tname.c_str()));
   tf.file_owned = false;
 
   if (!tf.tree) {
     throw failed_to_get_TTree()
         << "[ERROR]: Failed to get TTree named: " << std::quoted(tname)
         << " from TFile: " << std::quoted(file->GetName());
   }
   std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file
             << std::endl;
   return tf;
 }
 
 inline TreeFile CheckGetTTree(std::string const &fname,
                               std::string const &tname, char const *opts = "") {
   TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname);
   tf.file_owned = true;
   return tf;
 }
 
 inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname,
                              char const *opts = "") {
   TreeFile tf;
   tf.file = CheckOpenTFile(fname, opts);
   tf.tree = new TTree(tname.c_str(), "");
   tf.tree->SetDirectory(tf.file.get());
   tf.file_owned = true;
   return tf;
 }
 
-inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname) {
+inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname,
+                                  std::string const &dirname = "") {
   TreeFile tf;
   // This is pretty dumb...
   tf.file = make_unique_TFile(f);
   tf.tree = new TTree(tname.c_str(), "");
-  tf.tree->SetDirectory(tf.file.get());
+  if (dirname.size()) {
+    TDirectory *dir = tf.file->GetDirectory(dirname.c_str());
+    if (!dir) {
+      dir = tf.file->mkdir(dirname.c_str());
+    }
+    tf.tree->SetDirectory(dir);
+  } else {
+    tf.tree->SetDirectory(tf.file.get());
+  }
   tf.file_owned = false;
   return tf;
 }
 
 inline std::string SanitizeROOTObjectName(std::string name) { return name; }
 
 } // namespace utility
 } // namespace nuis
 
 #endif