Add muon decay background

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

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

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,13 +774,16 @@ 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);
       tdper[istage].CorrLi9EvtsSpec[idet]->SetDirectory(0);
       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,

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;
 

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;

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:

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

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

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

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' )]
 
 

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' )]
 

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

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");

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