SideBandLoad.C

Go to the documentation of this file.

// This macro creates spectra for OBSERVED
// and predicted event rates in sideband regions
// This uses UNBLIND DATA. Edit with extreme caution.

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "3FlavorAna/Vars/NumuVars.h"

#include <string>
#include <utility>

using namespace ana;

struct Sample {
  HistAxis* axnc;
  HistAxis* axnumu;
  Cut* nd;
  Cut* fd;
};

Sample MakeSample(const HistAxis& axnc, const HistAxis& axnm,
                  const Cut& nd, const Cut& fd);

std::vector<std::string> MakeUnblindList();

void SideBandLoad(std::string outfile)
{
  TH1::AddDirectory(0);

  // ND DATA + MC
  SpectrumLoader loaderNDMC(fnamenear_concat);
  SpectrumLoader loaderNDdata(fnameneardata_concat);
  loaderNDMC.SetSpillCut(kStandardSpillCuts);
  loaderNDdata.SetSpillCut(kStandardSpillCuts);
  
    // PERIOD 1, 2 + EPOCH 3B
  SpectrumLoader loaderFDMC_non(fFDMC_non);
  SpectrumLoader loaderFDMC_swp(fFDMC_swp);
  SpectrumLoader loaderFDMC_tau(fFDMC_tau);
  loaderFDMC_non.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_swp.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_tau.SetSpillCut(kStandardSpillCuts);

  // PERIOD 3C (and implied 3D)
  SpectrumLoader loaderFDMC_non_3c(fFDMC_non_3c);
  SpectrumLoader loaderFDMC_swp_3c(fFDMC_swp_3c);
  SpectrumLoader loaderFDMC_tau_3c(fFDMC_tau_3c);
  loaderFDMC_non_3c.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_swp_3c.SetSpillCut(kStandardSpillCuts);
  loaderFDMC_tau_3c.SetSpillCut(kStandardSpillCuts);

  // FD DATA
  const std::vector<std::string> fnameunblind(MakeUnblindList());
  SpectrumLoader loaderFDdata(fnameunblind);
  loaderFDdata.SetSpillCut(kStandardSpillCuts);

  const Binning kNCHighEBins = Binning::Simple(8, 4, 6);
  const HistAxis kNCHighEAxis( "Calorimetric Energy (GeV)",
                               kNCHighEBins, kCaloE );
  const HistAxis kHighENumuCCAxis( "Reconstructed Neutrino Energy (GeV)",
                                   kNCHighEBins, kCCE );

  std::string labelRecoE = "Calorimetric Energy (GeV)";
  const Cut kLowCVN = kNHit < 200 && kCVNnc < 0.2;
  const Cut kMidBDT({
      "sel.cosrej.numucontpid", "sel.nuecosrej.partptp",
      "vtx.elastic.fuzzyk.nshwlid", "vtx.elastic.fuzzyk.png.shwlid.lid.ismuon",
      "slc.calE", "slc.nhit","sel.nuecosrej.starttop",
      "sel.nuecosrej.stoptop"
    },
    [](const caf::SRProxy* sr)
    {
      double numucontpid = sr->sel.cosrej.numucontpid;
      double partptp = sr->sel.nuecosrej.partptp;
      if(sr->vtx.elastic.fuzzyk.nshwlid == 0) return false;
      if(numucontpid > 0.5) return false;
      if(numucontpid <= 0.42) return false;
      if(partptp >= 0.8) return false;
      if(sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon == 1) return false;
      double nhit = (double)sr->slc.nhit;
      if(nhit <= 0.) return false;
      if(sr->slc.calE/nhit <= 0.018) return false;
      if(std::min(sr->sel.nuecosrej.starttop, sr->sel.nuecosrej.stoptop) < 480) return false;
      return true;
    });

  std::map<std::string, Sample> sidebands;
  sidebands["HighE"] = MakeSample(
    kNCHighEAxis, kHighENumuCCAxis,
    kNusNDAllE && kCaloE > 4.,
    kNusFDAllE && kCaloE > 4.
  );
  sidebands["LowCVN"] = MakeSample(
    kNCAxis, kNCBinsNumuCCAxis,
    kNusNDPresel && kNusCosRejMod && kLowCVN && kNCCalERange,
    kNusFDPresel && kNusCosRej    && kLowCVN && kNCCalERange
  );
  sidebands["MidBDT"] = MakeSample(
    kNCAxis, kNCBinsNumuCCAxis,
    kNusND,
    kNusFDPresel && kNusNCSel && kMidBDT && kNCCalERange
  );

  std::map<std::string, IDecomp*> ncdecomps;
  std::map<std::string, IDecomp*> numudecomps;
  std::map<std::string, ModularExtrapSterile*> extrap123bs;
  std::map<std::string, ModularExtrapSterile*> extrap3c3ds;
  std::map<std::string, PredictionSterile*> predst123bs;
  std::map<std::string, PredictionSterile*> predst3c3ds;
  std::map<std::string, PredictionCombinePeriods*> preds;
  std::map<std::string, Spectrum*> cosmics;
  std::map<std::string, Spectrum*> dataspecs;

  for(const auto& sideband : sidebands) {
    std::string samplelabel = sideband.first;
    Sample sample = sideband.second;

    ncdecomps[samplelabel] = new ProportionalDecomp(
      loaderNDMC, loaderNDdata,
      *sample.axnc, *sample.nd,
      kNoShift, kNoShift, kTuftsWeightCC
    );

    numudecomps[samplelabel] = new ProportionalDecomp(
      loaderNDMC, loaderNDdata,
      *sample.axnumu, kNumuND,
      kNoShift, kNoShift, kTuftsWeightCC
    );

    extrap123bs[samplelabel] = new ModularExtrapSterile(
      ModularExtrapSterile::NCDisappearance(
        loaderNDMC, loaderFDMC_swp, loaderFDMC_non, loaderFDMC_tau,
        *ncdecomps[samplelabel], *numudecomps[samplelabel],
        *sample.axnc, *sample.axnumu,
        *sample.fd, *sample.nd, kNumuND,
        kNoShift, kTuftsWeightCC
      )
    );

    extrap3c3ds[samplelabel] = new ModularExtrapSterile(
      ModularExtrapSterile::NCDisappearance(
        loaderNDMC, loaderFDMC_swp_3c, loaderFDMC_non_3c, loaderFDMC_tau_3c,
        *ncdecomps[samplelabel], *numudecomps[samplelabel],
        *sample.axnc, *sample.axnumu,
        *sample.fd, *sample.nd, kNumuND,
        kNoShift, kTuftsWeightCC
      )
    );

    predst123bs[samplelabel] = new PredictionSterile(extrap123bs[samplelabel]);
    predst3c3ds[samplelabel] = new PredictionSterile(extrap3c3ds[samplelabel]);

    cosmics[samplelabel] = new Spectrum(loaderFDdata, *sample.axnc,
                                        kInTimingSideband && *sample.fd,
                                        kNoShift, kTimingSidebandWeight);

    dataspecs[samplelabel] = new Spectrum(loaderFDdata, *sample.axnc, kInBeamSpill && *sample.fd);
  }

  // Fill the spectra!
  loaderNDMC.Go();
  loaderNDdata.Go();
  loaderFDMC_non.Go();
  loaderFDMC_swp.Go();
  loaderFDMC_tau.Go();
  loaderFDMC_non_3c.Go();
  loaderFDMC_swp_3c.Go();
  loaderFDMC_tau_3c.Go();
  loaderFDdata.Go();

  // Two scale factors for scaling spectra, 1 data yr and 3 data yr
  const double pot_p1p2e3b = kSecondAnaPeriod1POT+kSecondAnaPeriod2POT+kSecondAnaEpoch3bPOT;
  const double pot_e3ce3d  = kSecondAnaEpoch3cPOT+kSecondAnaEpoch3dPOT;

  for(const auto& sideband : sidebands) {
    std::string samplelabel = sideband.first;
    preds[samplelabel] = new PredictionCombinePeriods({
      std::make_pair(predst123bs[samplelabel], pot_p1p2e3b),
      std::make_pair(predst3c3ds[samplelabel], pot_e3ce3d)
    });
  }

  // Set up output filename
  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");

  // Save all of the objects
  TDirectory* tmp = gDirectory;
  TDirectory* saveDir = gDirectory;
  std::string dir;
  std::string sep = "__";

  for(const auto& cafanaobj : ncdecomps) {
    dir = "ncdecomp" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : numudecomps) {
    dir = "numudecomp" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : extrap123bs) {
    dir = "extrap123b" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : extrap3c3ds) {
    dir = "extrap3c3d" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : predst123bs) {
    dir = "predst123b" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : predst3c3ds) {
    dir = "predst3c3d" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : preds) {
    dir = "pred" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : cosmics) {
    dir = "cosmic" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  for(const auto& cafanaobj : dataspecs) {
    dir = "dataspec" + sep + cafanaobj.first; // sample name
    saveDir = rootF->mkdir(dir.c_str());
    cafanaobj.second->SaveTo(saveDir);
  }

  tmp->cd();
  rootF->Close(); // Close the file
}

Sample MakeSample(const HistAxis& axnc, const HistAxis& axnm,
                  const Cut& nd, const Cut& fd)
{
  Sample ret;
  ret.axnc   = new HistAxis(axnc.label, axnc.bins, axnc.var);
  ret.axnumu = new HistAxis(axnm.label, axnm.bins, axnm.var);
  ret.nd = new Cut(nd);
  ret.fd = new Cut(fd);

  return ret;
}

std::vector<std::string> MakeUnblindList()
{
  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
  std::cout << "            WARNING: YOU ARE USING UNBLIND DATA." << std::endl;
  std::cout << "  Be very cautious about what distributions you are making," << std::endl;
  std::cout << "     where the results are going, and who will see them." << std::endl;
  std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;

  std::vector<std::string> fname1 = Wildcard(
    diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot/"
    + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_period1_nus_contain_v1_goodruns_prod2-snapshot*.root"
  );
  std::vector<std::string> fname2 = Wildcard(
    diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot/"
    + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_period2_nus_contain_v1_goodruns_prod2-snapshot*.root"
  );
  std::vector<std::string> fname3b = Wildcard(
    diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3*_nus_contain_v1_goodruns_prod2-snapshot.root"
  );
  std::vector<std::string> fname3c = Wildcard(
    diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot/"
    + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3c_nus_contain_v1_goodruns_prod2-snapshot*.root"
  );
  std::vector<std::string> fname3d = Wildcard(
    diskdirFD + "by_period/prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot/"
    + "prod_restricteddecaf_S16-05-20_fd_numi_fhc_epoch3d_nus_contain_v1_goodruns_prod2-snapshot*.root"
  );
  std::vector<std::string> ret;
  for(const auto& file : fname1)  { ret.push_back(file); }
  for(const auto& file : fname2)  { ret.push_back(file); }
  for(const auto& file : fname3b) { ret.push_back(file); }
  for(const auto& file : fname3c) { ret.push_back(file); }
  for(const auto& file : fname3d) { ret.push_back(file); }
  return ret;
}