diff --git a/ShapeFit/Paths.cc b/ShapeFit/Paths.cc
index 7645342..537891c 100644
--- a/ShapeFit/Paths.cc
+++ b/ShapeFit/Paths.cc
@@ -91,6 +91,11 @@ const char* aln()
   return "../alpha-n-spectrum/result-DocDB9667.root";
 }
 
+const char* muon_decay()
+{
+  return "../muon_decay_spectrum/MuonDecaySpec.root";
+}
+
 const char* baselines()
 {
   return "../toySpectra/unblinded_baseline.txt";
diff --git a/ShapeFit/Paths.h b/ShapeFit/Paths.h
index 763d5c5..adb04c0 100644
--- a/ShapeFit/Paths.h
+++ b/ShapeFit/Paths.h
@@ -35,6 +35,7 @@ const char* amc();
 // ihep=false. Go figure.
 const char* fn(bool ihep=true);
 const char* aln();
+const char* muon_decay();
 
 const char* baselines();
 
diff --git a/ShapeFit/Predictor.cc b/ShapeFit/Predictor.cc
index 1ae0405..ae2007a 100644
--- a/ShapeFit/Predictor.cc
+++ b/ShapeFit/Predictor.cc
@@ -1,5 +1,6 @@
 #include "Predictor.h"
 
+#include "Paths.h"
 #include "Utils.h"
 
 #include "TF1.h"
@@ -334,7 +335,23 @@ void Predictor::LoadMainData(const Char_t* mainmatrixname)
               tdper[istage - 1].DMCEff[ii] * tdper[istage - 1].MuonVetoEff[ii];
         }
       }
-
+      //-->muon decay bg
+      if (row == 21) {
+        for (int ii = 0; ii < Ndetectors; ii++) {
+          tdper[istage - 1].MuonDecayEvts[ii] = readvals[ii];
+          // in new input format have to correct for efficiencies here
+          tdper[istage - 1].MuonDecayEvts[ii] *=
+              tdper[istage - 1].DMCEff[ii] * tdper[istage - 1].MuonVetoEff[ii];
+        }
+      }
+      if (row == 22) {
+        for (int ii = 0; ii < Ndetectors; ii++) {
+          tdper[istage - 1].MuonDecayErr[ii] = readvals[ii];
+          // in new input format have to correct for efficiencies here
+          tdper[istage - 1].MuonDecayErr[ii] *=
+              tdper[istage - 1].DMCEff[ii] * tdper[istage - 1].MuonVetoEff[ii];
+        }
+      }
     } // only lines >1
     ++linenum;
   } // Reading main file
@@ -347,6 +364,7 @@ void Predictor::LoadMainData(const Char_t* mainmatrixname)
     tdper[istage - 1].BgEvts[ii] += tdper[istage - 1].Li9Evts[ii];
     tdper[istage - 1].BgEvts[ii] += tdper[istage - 1].FnEvts[ii];
     tdper[istage - 1].BgEvts[ii] += tdper[istage - 1].AlnEvts[ii];
+    tdper[istage - 1].BgEvts[ii] += tdper[istage - 1].MuonDecayEvts[ii];
   }
 
   // for(int ii=0;ii<Nstage;++ii){
@@ -470,9 +488,22 @@ void Predictor::LoadToyMCEntry(Int_t i, bool correct)
   RecalculateCovMatrix = true; // Force to recalculate
 }
 
+void Predictor::LoadBgSpec()
+{
+  // deliberate leak
+  TString* accLoc = new TString[3];
+  accLoc[0] = Paths::acc_spectra(0);
+  accLoc[1] = Paths::acc_spectra(1);
+  accLoc[2] = Paths::acc_spectra(2);
+
+  LoadBgSpec(accLoc, Paths::li9(), Paths::amc(), Paths::fn(),
+             Paths::aln(), Paths::muon_decay());
+}
+
 void Predictor::LoadBgSpec(TString* accspecname, const Char_t* li9specname,
                            const Char_t* amcspecname, const Char_t* fnspecname,
-                           const Char_t* alnspecname)
+                           const Char_t* alnspecname,
+                           const Char_t* muondecayspecname)
 {
   cout << "Loading bg spectra..." << endl;
   Char_t name[1024];
@@ -688,6 +719,45 @@ void Predictor::LoadBgSpec(TString* accspecname, const Char_t* li9specname,
   // m_alnspec->Close();
   cout << "--> loaded alpha-n spectra" << endl;
 
+  m_muondecayspec = new TFile(muondecayspecname, "READ");
+  for (int istage = 0; istage < Nstage; ++istage) {
+    for (int idet = 0; idet < Ndetectors; ++idet) {
+      sprintf(name, "MdSpec_EH%d;1", detConfigEH[idet]);
+      TH1F* htemp;
+      htemp = (TH1F*) m_muondecayspec->Get(name)->Clone();
+      sprintf(name, "CorrMuonDecayEvtsSpec_stage%d_ad%d", istage, idet);
+      tdper[istage].CorrMuonDecayEvtsSpec[idet] =
+        (TH1F*)tdper[istage].ObsEvtsSpec[idet]->Clone(name);
+      tdper[istage].CorrMuonDecayEvtsSpec[idet]->Reset();
+
+      // fill into destination histogram
+      for (Int_t ibin = 0; ibin < htemp->GetNbinsX(); ibin++) {
+        tdper[istage].CorrMuonDecayEvtsSpec[idet]->Fill(
+            htemp->GetBinCenter(ibin + 1), htemp->GetBinContent(ibin + 1));
+      }
+
+      delete htemp;
+
+      // assign zero errors
+      for (Int_t ibin = 0;
+           ibin < tdper[istage].CorrMuonDecayEvtsSpec[idet]->GetNbinsX(); ibin++) {
+        tdper[istage].CorrMuonDecayEvtsSpec[idet]->SetBinError(ibin + 1, 0);
+      }
+
+      // scale
+      if (tdper[istage].CorrMuonDecayEvtsSpec[idet]->Integral() == 0)
+        tdper[istage].CorrMuonDecayEvtsSpec[idet]->Scale(0);
+      else
+        tdper[istage].CorrMuonDecayEvtsSpec[idet]->Scale(
+            tdper[istage].MuonDecayEvts[idet] * 1. /
+            tdper[istage]
+                .CorrMuonDecayEvtsSpec[idet]
+                ->Integral()); //<---Scale this histogram to expected number of
+                               // events based on input file
+    }
+  }
+  cout << "--> loaded muon decay spectra" << endl;
+
   cout << "done loading bg spectra!" << endl;
 
   // create sum of bg's spectra
@@ -704,6 +774,8 @@ void Predictor::LoadBgSpec(TString* accspecname, const Char_t* li9specname,
           tdper[istage].CorrFnEvtsSpec[idet], 1);
       tdper[istage].CorrBgEvtsSpec[idet]->Add(
           tdper[istage].CorrAlnEvtsSpec[idet], 1);
+      tdper[istage].CorrBgEvtsSpec[idet]->Add(
+          tdper[istage].CorrMuonDecayEvtsSpec[idet], 1);
 
       tdper[istage].CorrBgEvtsSpec[idet]->SetDirectory(0);
       tdper[istage].CorrAccEvtsSpec[idet]->SetDirectory(0);
@@ -711,6 +783,7 @@ void Predictor::LoadBgSpec(TString* accspecname, const Char_t* li9specname,
       tdper[istage].CorrAmcEvtsSpec[idet]->SetDirectory(0);
       tdper[istage].CorrFnEvtsSpec[idet]->SetDirectory(0);
       tdper[istage].CorrAlnEvtsSpec[idet]->SetDirectory(0);
+      tdper[istage].CorrMuonDecayEvtsSpec[idet]->SetDirectory(0);
     }
   }
 
@@ -731,6 +804,11 @@ void Predictor::LoadBgSpec(TString* accspecname, const Char_t* li9specname,
 
 } // end of LoadBgSpec
 
+// TimePeriodData.CorrectEvts adjusts tdper.*Evts to give the daily rate
+// assuming perfect muon/DMC eff, so, assuming that GetCorr*EvtsSpec() is called
+// after TimerPeriodData.CorrectEvts, scaling by livetime/effs gives the raw
+// background expected in the data.
+
 // Extract Bkg Spectrum here
 TH1F* Predictor::GetCorrAccEvtsSpec(Int_t istage, Int_t idet)
 {
@@ -777,6 +855,16 @@ TH1F* Predictor::GetCorrAlnEvtsSpec(Int_t istage, Int_t idet)
   // return
   // tdper[istage].CorrAlnEvtsSpec[idet]->Scale(tdper[istage].Livetime[idet]);
 }
+TH1F* Predictor::GetCorrMuonDecayEvtsSpec(Int_t istage, Int_t idet)
+{
+  TH1F* dummy = (TH1F*)tdper[istage].CorrMuonDecayEvtsSpec[idet]->Clone();
+  dummy->Scale(tdper[istage].Livetime[idet] * tdper[istage].MuonVetoEff[idet] *
+               tdper[istage].DMCEff[idet]);
+  return dummy;
+  // return
+  // tdper[istage].CorrMuonDecayEvtsSpec[idet]->Scale(tdper[istage].Livetime[idet]);
+}
+
 // Done!
 
 void Predictor::EnterObsEvts(double obsad1, double obsad2, double obsad3,
diff --git a/ShapeFit/Predictor.h b/ShapeFit/Predictor.h
index 66e9ff1..380c48b 100644
--- a/ShapeFit/Predictor.h
+++ b/ShapeFit/Predictor.h
@@ -31,9 +31,10 @@ class Predictor : public TObject {
   Int_t LoadToyIBDSpec(const Char_t* toyibdspecname);
   void LoadToyMCEntry(Int_t i, bool correct = true);
   void LoadToyMCNominalSpec();
+  void LoadBgSpec();
   void LoadBgSpec(TString* accspecname, const Char_t* li9specname,
                   const Char_t* amcspecname, const Char_t* fnspecname,
-                  const Char_t* alnspecname);
+                  const Char_t* alnspecname, const Char_t* muondecayspecname);
   FluxCalculator* GetFluxCalculator() { return fluxcalc; }
 
   void MakePrediction(double sin2theta13, double dm2, Double_t sin22t14,
@@ -64,6 +65,7 @@ class Predictor : public TObject {
   TH1F* GetCorrAmcEvtsSpec(int istage, int idet);
   TH1F* GetCorrFnEvtsSpec(int istage, int idet);
   TH1F* GetCorrAlnEvtsSpec(int istage, int idet);
+  TH1F* GetCorrMuonDecayEvtsSpec(int istage, int idet);
 
   void LoadEvisToEnuMatrix(const Char_t* evis_to_enu_matrix_name);
 
@@ -271,6 +273,7 @@ class Predictor : public TObject {
   TFile* m_fnspec;
   TFile* m_amcspec;
   TFile* m_alnspec;
+  TFile* m_muondecayspec;
 
   PredSet* superpred;
 
diff --git a/ShapeFit/TimePeriodData.cc b/ShapeFit/TimePeriodData.cc
index 27760f4..e5aa0e8 100644
--- a/ShapeFit/TimePeriodData.cc
+++ b/ShapeFit/TimePeriodData.cc
@@ -11,11 +11,13 @@ TimePeriodData::TimePeriodData()
     AmcEvts[idet1] = 0;
     AlnEvts[idet1] = 0;
     FnEvts[idet1] = 0;
+    MuonDecayEvts[idet1] = 0;
     AccErr[idet1] = 0;
     Li9Err[idet1] = 0;
     AmcErr[idet1] = 0;
     AlnErr[idet1] = 0;
     FnErr[idet1] = 0;
+    MuonDecayErr[idet1] = 0;
     CorrBgEvts[idet1] = 0;
     ErrBgEvts[idet1] = 0;
     MuonVetoEff[idet1] = 0;
@@ -51,11 +53,13 @@ void TimePeriodData::CorrectEvts(bool print)
     AlnEvts[idet] *= Livetime[idet];
     Li9Evts[idet] *= Livetime[idet];
     FnEvts[idet] *= Livetime[idet];
+    MuonDecayEvts[idet] *= Livetime[idet];
     AccErr[idet] *= Livetime[idet];
     AmcErr[idet] *= Livetime[idet];
     AlnErr[idet] *= Livetime[idet];
     Li9Err[idet] *= Livetime[idet];
     FnErr[idet] *= Livetime[idet];
+    MuonDecayErr[idet] *= Livetime[idet];
   }
 
   // Scale back to 1 day of livetime and correct for muon efficiencies
@@ -81,12 +85,14 @@ void TimePeriodData::CorrectEvts(bool print)
     Li9Evts[idet] *= factor_bg;
     FnEvts[idet] *= factor_bg;
     AlnEvts[idet] *= factor_bg;
+    MuonDecayEvts[idet] *= factor_bg;
 
     AccErr[idet] *= factor_bg;
     AmcErr[idet] *= factor_bg;
     Li9Err[idet] *= factor_bg;
     FnErr[idet] *= factor_bg;
     AlnErr[idet] *= factor_bg;
+    MuonDecayErr[idet] *= factor_bg;
 
     ErrEvts[idet] = ErrEvts[idet] * factor;
     CorrBgEvts[idet] = BgEvtsLiv[idet] * factor_bg;
diff --git a/ShapeFit/TimePeriodData.h b/ShapeFit/TimePeriodData.h
index b241afd..d0a870f 100644
--- a/ShapeFit/TimePeriodData.h
+++ b/ShapeFit/TimePeriodData.h
@@ -34,11 +34,13 @@ class TimePeriodData : public TObject {
   Double_t AmcEvts[Ndetectors] = {};
   Double_t AlnEvts[Ndetectors] = {};
   Double_t FnEvts[Ndetectors] = {};
+  Double_t MuonDecayEvts[Ndetectors] = {};
   Double_t AccErr[Ndetectors] = {};
   Double_t Li9Err[Ndetectors] = {};
   Double_t AmcErr[Ndetectors] = {};
   Double_t AlnErr[Ndetectors] = {};
   Double_t FnErr[Ndetectors] = {};
+  Double_t MuonDecayErr[Ndetectors] = {};
   Double_t CorrBgEvts[Ndetectors] = {};
   Double_t ErrBgEvts[Ndetectors] =
       {}; // statistical uncertainties for the backgrounds
@@ -56,6 +58,7 @@ class TimePeriodData : public TObject {
   TH1F* CorrAmcEvtsSpec[Ndetectors] = {};
   TH1F* CorrFnEvtsSpec[Ndetectors] = {};
   TH1F* CorrAlnEvtsSpec[Ndetectors] = {};
+  TH1F* CorrMuonDecayEvtsSpec[Ndetectors] = {};
   TH1F* CorrEvtsSpec[Ndetectors] = {};
 
 private:
diff --git a/ShapeFit/build_covmatrix.C b/ShapeFit/build_covmatrix.C
index 7bf9e48..12e3e5d 100644
--- a/ShapeFit/build_covmatrix.C
+++ b/ShapeFit/build_covmatrix.C
@@ -66,11 +66,7 @@ void build_covmatrix(const Char_t* toymc_filename,
 
     myPred->LoadToyMCEntry(0, false);
 
-    TString AccidentalSpectrumLocation[3] = {
-      Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
-
-    myPred->LoadBgSpec(AccidentalSpectrumLocation,
-                       Paths::li9(), Paths::amc(), Paths::fn(), Paths::aln());
+    myPred->LoadBgSpec();
 
     myPred->SetEvisBins(n_evis_bins, evis_bins);
     myPred->SetEnuBins(n_enu_bins, enu_bins);
diff --git a/ShapeFit/fit_shape_2d_P17B.C b/ShapeFit/fit_shape_2d_P17B.C
index ff07372..32c3a0d 100644
--- a/ShapeFit/fit_shape_2d_P17B.C
+++ b/ShapeFit/fit_shape_2d_P17B.C
@@ -59,8 +59,6 @@ void fit_shape_2d_P17B(
 
   TString sig_spectra_filename[3] = {
       Paths::sig_spectra(0), Paths::sig_spectra(1), Paths::sig_spectra(2)};
-  TString AccidentalSpectrumLocation[3] = {
-      Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
   Char_t Theta13InputsLocation[3][1024];
 
   strcpy(Theta13InputsLocation[0], Paths::input(0));
@@ -106,8 +104,7 @@ void fit_shape_2d_P17B(
     pred->LoadIBDSpec(sig_spectra_filename);
 
     // load bg afterwards since here is when correct events
-    pred->LoadBgSpec(AccidentalSpectrumLocation,
-                     Paths::li9(), Paths::amc(), Paths::fn(), Paths::aln());
+    pred->LoadBgSpec();
 
 
     // pred->SetStatFactor(stat_factor);
@@ -992,6 +989,9 @@ void fit_shape_2d_P17B(
 
       sprintf(name, "CorrAlnEvtsSpec_stage%i_AD%i", istage, idet + 1);
       pred->GetCorrAlnEvtsSpec(istage, idet)->Clone(name)->Write();
+
+      sprintf(name, "CorrMuonDecayEvtsSpec_stage%i_AD%i", istage, idet + 1);
+      pred->GetCorrMuonDecayEvtsSpec(istage, idet)->Clone(name)->Write();
     }
   }
 
diff --git a/muon_decay_spectrum/MuonDecaySpec.root b/muon_decay_spectrum/MuonDecaySpec.root
new file mode 100644
index 0000000..c6dfa3e
Binary files /dev/null and b/muon_decay_spectrum/MuonDecaySpec.root differ
diff --git a/toySpectra/data_file/dyb_data_v1_nominal.txt b/toySpectra/data_file/dyb_data_v1_nominal.txt
index 43abfdf..33ea4cb 100644
--- a/toySpectra/data_file/dyb_data_v1_nominal.txt
+++ b/toySpectra/data_file/dyb_data_v1_nominal.txt
@@ -163,7 +163,8 @@ varyAccBg 0
 varyAmcBg 0
 varyFnBg 0
 varyLi9Bg 0
-varyAlnBg 0 
+varyAlnBg 0
+varyMuonDecayBg 0
 
 #These distort the shape of the background spectra by a line 
 #Note: the value is the amount of variation that you want to introduce
@@ -172,6 +173,7 @@ distortAmcBg 0
 distortFnBg 0
 distortLi9Bg null
 distortAlnBg 0
+distortMuonDecayBg 0
 
 #Finally, flag for statistical flactuation
 statisticalFluctuation 0
diff --git a/toySpectra/data_file/generate_data_file.py b/toySpectra/data_file/generate_data_file.py
index 9e99666..fabbe0d 100755
--- a/toySpectra/data_file/generate_data_file.py
+++ b/toySpectra/data_file/generate_data_file.py
@@ -20,9 +20,11 @@
               ('varyFnBg', 'vary_fn', '1'),
               ('varyLi9Bg', 'vary_li9', '1'),
               ('varyAlnBg', 'vary_aln', '1' ),
+              ('varyMuonDecayBg', 'vary_muon_decay', '1' ),
               ('distortAmcBg', 'distort_amc', '0.15'),
               ('distortFnBg', 'distort_fn', '0.2'),
               ('distortLi9Bg', 'distort_li9', '../li9_spectrum/8he9li_distort_neutron100_alpha100_frac0.055_N250.root'),
+              # ('distortMuonDecayBg', 'distort_muon_decay', '0.12345'), # TODO
               ('statisticalFluctuation', 'stat', '1' )]
 
 
diff --git a/toySpectra/data_file/generate_data_file_extra.py b/toySpectra/data_file/generate_data_file_extra.py
index 30fdc40..1170da8 100755
--- a/toySpectra/data_file/generate_data_file_extra.py
+++ b/toySpectra/data_file/generate_data_file_extra.py
@@ -20,9 +20,11 @@
               ('varyFnBg', 'vary_fn', '1'),
               ('varyLi9Bg', 'vary_li9', '1'),
               ('varyAlnBg', 'vary_aln', '1' ),
+              ('varyMuonDecayBg', 'vary_muon_decay', '1' ),
               ('distortAmcBg', 'distort_amc', '0.15'),
               ('distortFnBg', 'distort_fn', '0.2'),
               ('distortLi9Bg', 'distort_li9', '../li9_spectrum/8he9li_distort_neutron100_alpha100_frac0.1_N250.root'),
+              # ('distortMuonDecayBg', 'distort_muon_decay', '0.12345'), # TODO
               ('statisticalFluctuation', 'stat', '1' ),
               ('useBcwFluxUncertainty', 'bcwflux', '1' )]
 
diff --git a/toySpectra/genEvisToEnuMatrix.C b/toySpectra/genEvisToEnuMatrix.C
index 1c7cf34..84e57a3 100644
--- a/toySpectra/genEvisToEnuMatrix.C
+++ b/toySpectra/genEvisToEnuMatrix.C
@@ -90,12 +90,7 @@ void genEvisToEnuMatrix(Double_t s2t13 = -1, Double_t dm2ee = -1,
     spectrumNormNominal->loadDistances(Paths::baselines());
     spectrumNormNominal->initialize(mydata_nominal);
 
-    TString AccidentalSpectrumLocation[3] = {
-      Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
-
-    spectrumNormNominal->loadBgSpecForToy(AccidentalSpectrumLocation,
-                                          Paths::li9(), Paths::amc(), Paths::fn(),
-                                          Paths::aln());
+    spectrumNormNominal->loadBgSpecForToy();
 
     // Generate nominal spectrum
     spectrumNormNominal->updateAntinu();
diff --git a/toySpectra/genPredictedIBD.C b/toySpectra/genPredictedIBD.C
index af9c04e..cde7570 100644
--- a/toySpectra/genPredictedIBD.C
+++ b/toySpectra/genPredictedIBD.C
@@ -89,12 +89,7 @@ void genPredictedIBD(Double_t s2t13 = 0, Double_t dm2ee = -1,
   spectrumNormNominal->loadDistances(Paths::baselines());
   spectrumNormNominal->initialize(mydata_nominal);
 
-  TString AccidentalSpectrumLocation[3] = {
-    Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
-
-  spectrumNormNominal->loadBgSpecForToy(AccidentalSpectrumLocation,
-                                        Paths::li9(), Paths::amc(), Paths::fn(),
-                                        Paths::aln());
+  spectrumNormNominal->loadBgSpecForToy();
 
   // Prepare destination file
   TFile* outfile = new TFile(Paths::predicted_ibd(), "RECREATE");
diff --git a/toySpectra/genSuperHistograms.C b/toySpectra/genSuperHistograms.C
index 5a6e756..3cea9f7 100644
--- a/toySpectra/genSuperHistograms.C
+++ b/toySpectra/genSuperHistograms.C
@@ -46,12 +46,7 @@ void genSuperHistograms()
   // TString AccidentalSpectrumLocation[2] =
   // {"../ShapeFit/Spectra/accidental_eprompt_shapes_6ad.root","../ShapeFit/Spectra/accidental_eprompt_shapes_8ad_p14a.root"};
 
-  TString AccidentalSpectrumLocation[3] = {
-    Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
-
-  spectrumNorm->loadBgSpecForToy(AccidentalSpectrumLocation,
-                                 Paths::li9(), Paths::amc(), Paths::fn(),
-                                 Paths::aln());
+  spectrumNorm->loadBgSpecForToy();
 
   // load distances hare as well to construct traditional supermatrix
 
diff --git a/toySpectra/genToySpectraTree.C b/toySpectra/genToySpectraTree.C
index c074a77..6fd59a0 100644
--- a/toySpectra/genToySpectraTree.C
+++ b/toySpectra/genToySpectraTree.C
@@ -108,12 +108,7 @@ void genToySpectraTree(TString dataset_filename, TString output_filename,
   spectrumNormNominal->loadDistances(Paths::baselines());
   spectrumNormNominal->initialize(mydata_nominal);
 
-  TString AccidentalSpectrumLocation[3] = {
-      Paths::acc_spectra(0), Paths::acc_spectra(1), Paths::acc_spectra(2)};
-
-  spectrumNormNominal->loadBgSpecForToy(AccidentalSpectrumLocation,
-                                        Paths::li9(), Paths::amc(), Paths::fn(),
-                                        Paths::aln());
+  spectrumNormNominal->loadBgSpecForToy();
 
 
   // Create Spectrum for toy
@@ -127,9 +122,7 @@ void genToySpectraTree(TString dataset_filename, TString output_filename,
     // to avoid duplication
     spectrumNorm->loadDistances(Paths::baselines());
     spectrumNorm->initialize(mydata);
-    spectrumNorm->loadBgSpecForToy(AccidentalSpectrumLocation,
-                                   Paths::li9(), Paths::amc(), Paths::fn(),
-                                   Paths::aln());
+    spectrumNorm->loadBgSpecForToy();
     // NB: Adding 1 is important because 0 means random (UUID) seed
     spectrumNorm->setRandomSeed(1 + omp_get_thread_num());
   } // parallel
diff --git a/toySpectra/reactor/Spectrum.C b/toySpectra/reactor/Spectrum.C
index bc94f2d..d01fcb9 100644
--- a/toySpectra/reactor/Spectrum.C
+++ b/toySpectra/reactor/Spectrum.C
@@ -993,10 +993,23 @@ void Spectrum::loadDistances(const char* distancematrixname)
   }
 }
 
+void Spectrum::loadBgSpecForToy()
+{
+  // deliberate leak
+  TString* accLoc = new TString[3];
+  accLoc[0] = Paths::acc_spectra(0);
+  accLoc[1] = Paths::acc_spectra(1);
+  accLoc[2] = Paths::acc_spectra(2);
+
+  loadBgSpecForToy(accLoc, Paths::li9(), Paths::amc(), Paths::fn(),
+                   Paths::aln(), Paths::muon_decay());
+}
+
 void Spectrum::loadBgSpecForToy(TString* accspecname, const Char_t* li9specname,
                                 const Char_t* amcspecname,
                                 const Char_t* fnspecname,
-                                const Char_t* alnspecname)
+                                const Char_t* alnspecname,
+                                const Char_t* muondecayspecname)
 {
   cout << "Loading bg spectra..." << endl;
   Char_t name[1024];
@@ -1160,6 +1173,29 @@ void Spectrum::loadBgSpecForToy(TString* accspecname, const Char_t* li9specname,
   // m_alnspec->Close();
   cout << "--> loaded alpha-n spectra" << endl;
 
+  //(muon decay)
+  m_muondecayspec = new TFile(muondecayspecname, "READ");
+  dir->cd();
+  for (int istage = 0; istage < Nstage; ++istage) {
+    for (int idet = 0; idet < Ndetectors; ++idet) {
+      sprintf(name, "CorrMuonDecayEvtsSpec_ad%d", idet);
+      Char_t nameMd[1024];
+      sprintf(nameMd, "MdSpec_EH%d;1", detConfigEH[idet]);
+
+      CorrMuonDecayEvtsSpec[istage][idet] = (TH1F*)m_muondecayspec->Get(nameMd)->Clone(name);
+
+      for (Int_t ibin = 0; ibin < CorrMuonDecayEvtsSpec[istage][idet]->GetNbinsX();
+           ibin++) {
+        CorrMuonDecayEvtsSpec[istage][idet]->SetBinError(ibin + 1, 0);
+      }
+      CorrMuonDecayEvtsSpec[istage][idet]->Scale(
+          pred->tdper[istage].MuonDecayEvts[idet] /
+          CorrMuonDecayEvtsSpec[istage][idet]->Integral());
+    }
+  }
+  // m_muondecayspec->Close();
+  cout << "--> loaded muon decay spectra" << endl;
+
   cout << "done loading bg spectra!" << endl;
 
   // un-correcting the spectrum normalization....
@@ -1177,6 +1213,7 @@ void Spectrum::loadBgSpecForToy(TString* accspecname, const Char_t* li9specname,
       CorrLi9EvtsSpec[istage][idet]->Scale(factor);
       CorrFnEvtsSpec[istage][idet]->Scale(factor);
       CorrAlnEvtsSpec[istage][idet]->Scale(factor);
+      CorrMuonDecayEvtsSpec[istage][idet]->Scale(factor);
     }
   }
   /*
@@ -1204,13 +1241,14 @@ void Spectrum::loadBgSpecForToy(TString* accspecname, const Char_t* li9specname,
 } // end of LoadBgSpec
 
 void Spectrum::setBgRemoveFlag(bool acc_flag, bool li9_flag, bool fn_flag,
-                               bool amc_flag, bool aln_flag)
+                               bool amc_flag, bool aln_flag, bool muon_decay_flag)
 {
   m_removeAccBg = acc_flag;
   m_removeLi9Bg = li9_flag;
   m_removeFnBg = fn_flag;
   m_removeAmcBg = amc_flag;
   m_removeAlnBg = aln_flag;
+  m_removeMuonDecayBg = muon_decay_flag;
 
   this->updateBgDetected();
 }
@@ -1224,6 +1262,7 @@ void Spectrum::updateBgDetected()
   TH1F* CopyAlnEvtsSpec[Nstage][Ndetectors];
   TH1F* CopyLi9EvtsSpec[Nstage][Ndetectors];
   TH1F* CopyFnEvtsSpec[Nstage][Ndetectors];
+  TH1F* CopyMuonDecayEvtsSpec[Nstage][Ndetectors];
   Char_t name[1024];
   for (int istage = 0; istage < Nstage; ++istage) {
     for (int idet = 0; idet < Ndetectors; ++idet) {
@@ -1242,6 +1281,9 @@ void Spectrum::updateBgDetected()
       sprintf(name, "h_aln_copy_%i", idet);
       CopyAlnEvtsSpec[istage][idet] =
           (TH1F*)CorrAlnEvtsSpec[istage][idet]->Clone(name);
+      sprintf(name, "h_muondecay_copy_%i", idet);
+      CopyMuonDecayEvtsSpec[istage][idet] =
+          (TH1F*)CorrMuonDecayEvtsSpec[istage][idet]->Clone(name);
     }
   }
   // ****************************************************
@@ -1253,6 +1295,7 @@ void Spectrum::updateBgDetected()
   double scale_factor_li9[Nstage][Ndetectors];
   double scale_factor_fn[Nstage][Ndetectors];
   double scale_factor_aln[Nstage][Ndetectors];
+  double scale_factor_muondecay[Nstage][Ndetectors];
   double corr_random[Nhalls] = {0};
 
   for (int istage = 0; istage < Nstage; ++istage) {
@@ -1262,11 +1305,12 @@ void Spectrum::updateBgDetected()
       scale_factor_li9[istage][idet] = 1;
       scale_factor_fn[istage][idet] = 1;
       scale_factor_aln[istage][idet] = 1;
+      scale_factor_muondecay[istage][idet] = 1;
     }
   }
   // Fluctuate each background depending on how it makes the most sense
   // 1) For accidentals, fluctuate each AD independenlty
-  // 2) For li9 and fn, fluctuate ADs in same site with same factor
+  // 2) For li9 and fn and muon decay, fluctuate ADs in same site with same factor
   // 3) For amc, fluctuate all ADs in a correlated way
   // 4) For alpha-n, for now fluctuating all ADs independently
 
@@ -1346,6 +1390,26 @@ void Spectrum::updateBgDetected()
     }
   }
 
+  // muon decays are correlated between detectors in the same hall(?)
+  for (int ihall = 0; ihall < Nhalls; ++ihall) {
+    corr_random[ihall] = ran->Gaus(0, 1);
+  }
+
+  for (int istage = 0; istage < Nstage; ++istage) {
+    for (int idet = 0; idet < Ndetectors; ++idet) {
+      if (m_varyMuonDecayBg > 0) {
+        if (pred->tdper[istage].MuonDecayEvts[idet] > 0)
+          scale_factor_muondecay[istage][idet] =
+              (1 + pred->tdper[istage].MuonDecayErr[idet] * 1. /
+                       pred->tdper[istage].MuonDecayEvts[idet] *
+                       corr_random[detConfigEH[idet] - 1]);
+        else
+          scale_factor_muondecay[istage][idet] = 0;
+      }
+    }
+  }
+
+
   // force set background contribution to 0 if the flag is set
   for (int istage = 0; istage < Nstage; ++istage) {
     for (int idet = 0; idet < Ndetectors; ++idet) {
@@ -1359,6 +1423,8 @@ void Spectrum::updateBgDetected()
         scale_factor_li9[istage][idet] = 0;
       if (m_removeFnBg)
         scale_factor_fn[istage][idet] = 0;
+      if (m_removeMuonDecayBg)
+        scale_factor_muondecay[istage][idet] = 0;
     }
   }
 
@@ -1370,6 +1436,7 @@ void Spectrum::updateBgDetected()
   // - li9: fluctuate shape of all ADs together
   // - fn: fluctuate shape of ADs in same site together
   // - aln: fluctuate shape of all ADs together
+  // - muon decay: TODO
   // *******************************************************
   if (m_distortAccBg > 0) {
     TF1* func_acc = getDistortionCurve(m_distortAccBg);
@@ -1474,6 +1541,10 @@ void Spectrum::updateBgDetected()
     delete func_aln;
   }
 
+  if (m_distortMuonDecayBg > 0) {
+    // TODO
+  }
+
 
   // ********************************************************
   // Returning the bg spectra
@@ -1498,6 +1569,9 @@ void Spectrum::updateBgDetected()
         m_bgDetectedSpectrum[istage][idet][ibin] +=
             scale_factor_aln[istage][idet] *
             CopyAlnEvtsSpec[istage][idet]->GetBinContent(ibin + 1);
+        m_bgDetectedSpectrum[istage][idet][ibin] +=
+            scale_factor_muondecay[istage][idet] *
+            CopyMuonDecayEvtsSpec[istage][idet]->GetBinContent(ibin + 1);
       }
 
       delete CopyAccEvtsSpec[istage][idet];
@@ -1505,6 +1579,7 @@ void Spectrum::updateBgDetected()
       delete CopyLi9EvtsSpec[istage][idet];
       delete CopyFnEvtsSpec[istage][idet];
       delete CopyAlnEvtsSpec[istage][idet];
+      delete CopyMuonDecayEvtsSpec[istage][idet];
     }
   }
 
@@ -1779,6 +1854,7 @@ void Spectrum::initialize(DataSet* data)
   m_varyFnBg = data->getDouble("varyFnBg");
   m_varyLi9Bg = data->getDouble("varyLi9Bg");
   m_varyAlnBg = data->getDouble("varyAlnBg");
+  m_varyMuonDecayBg = data->getDouble("varyMuonDecayBg");
 
   m_distortAccBg = data->getDouble("distortAccBg");
   m_distortAmcBg = data->getDouble("distortAmcBg");
@@ -1793,6 +1869,7 @@ void Spectrum::initialize(DataSet* data)
   }
   m_distortFnBg = data->getDouble("distortFnBg");
   m_distortAlnBg = data->getDouble("distortAlnBg");
+  m_distortMuonDecayBg = data->getDouble("distortMuonDecayBg");
 
   m_statisticalFluctuation = data->getDouble("statisticalFluctuation");
 
diff --git a/toySpectra/reactor/Spectrum.h b/toySpectra/reactor/Spectrum.h
index f044440..0ca77b9 100644
--- a/toySpectra/reactor/Spectrum.h
+++ b/toySpectra/reactor/Spectrum.h
@@ -86,9 +86,10 @@ class Spectrum {
   // which are loaded by Predictor)
   void loadBgSpecForToy(TString* accspecname, const Char_t* li9specname,
                         const Char_t* amcspecname, const Char_t* fnspecname,
-                        const Char_t* alnspecname);
+                        const Char_t* alnspecname, const Char_t* muondecayspecname);
+  void loadBgSpecForToy();
   void setBgRemoveFlag(bool acc_flag, bool li9_flag, bool fn_flag,
-                       bool amc_flag, bool aln_flag);
+                       bool amc_flag, bool aln_flag, bool muon_decay_flag);
 
   // Update Oscillation
   // void setDeltaMSqee(double deltaMSqee);
@@ -137,6 +138,7 @@ class Spectrum {
   TH1F* CorrAlnEvtsSpec[Nstage][Ndetectors];
   TH1F* CorrLi9EvtsSpec[Nstage][Ndetectors];
   TH1F* CorrFnEvtsSpec[Nstage][Ndetectors];
+  TH1F* CorrMuonDecayEvtsSpec[Nstage][Ndetectors];
 
   void setRandomSolarOscPars();
   void setRandomDm2ee();
@@ -343,6 +345,7 @@ class Spectrum {
   double m_varyFnBg;
   double m_varyLi9Bg;
   double m_varyAlnBg;
+  double m_varyMuonDecayBg;
 
   double m_distortAccBg;
   double m_distortAmcBg;
@@ -353,6 +356,7 @@ class Spectrum {
   TH1F* m_func_distortLi9Bg;
   int m_entries_distortLi9Bg;
   double m_distortAlnBg;
+  double m_distortMuonDecayBg;
 
   double m_statisticalFluctuation;
 
@@ -362,7 +366,7 @@ class Spectrum {
   bool m_removeFnBg;
   bool m_removeAmcBg;
   bool m_removeAlnBg;
-
+  bool m_removeMuonDecayBg;
 
   void extractPredictorData();
   void loadIavCorrection(const char* iavcorrectionname);
@@ -411,6 +415,7 @@ class Spectrum {
   TFile* m_fnspec;
   TFile* m_amcspec;
   TFile* m_alnspec;
+  TFile* m_muondecayspec;
 
   TFile* m_iavCorrFile;
 };