NDDataMCSystBandLoad.C

Go to the documentation of this file.

// Macro to plot ND Data/MC spectra with systematic error band
// To properly include NusFlatSysts (MC Stat, Norm, ND Cont),
// comment out the line in NusSysts.cxx, NusFlatSyst::Shift
// that only shifts the FD

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/Cuts.h"
#include "NuXAna/Cuts/NusCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "NuXAna/Systs/NusSysts.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/HistAxes.h"
#include "NuXAna/Vars/HistAxes.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"

#include "TFile.h"

#include <iostream>
#include <string>
#include <vector>

using namespace ana;

void NDDataMCSystBandLoad()
{
  TH1::AddDirectory(0);

  SpectrumLoader floadnear(fnamenear_concat);
  SpectrumLoader floadneardata(fnameneardata_concat);

  floadnear    .SetSpillCut(kStandardSpillCuts);
  floadneardata.SetSpillCut(kStandardSpillCuts);

  std::string labelHits    = "Number of Hits";
  Binning kHitBins  = Binning::Simple(220, 0., 220.);
  const HistAxis kAxisNHit(labelHits,    kHitBins,      kNHit);

  // Set up nominal spectra
  Spectrum sMCCalE(floadnear,     kNCAxis, kNusND, kNoShift, kTuftsWeightCC);
  Spectrum sDtCalE(floadneardata, kNCAxis, kNusND);
  CheatDecomp dcCalE(floadnear, kNCAxis, kNusND, kNoShift, kTuftsWeightCC);

  Spectrum sMCNHit(floadnear,     kAxisNHit, kNusND, kNoShift, kTuftsWeightCC);
  Spectrum sDtNHit(floadneardata, kAxisNHit, kNusND);
  CheatDecomp dcNHit(floadnear, kAxisNHit, kNusND, kNoShift, kTuftsWeightCC);

  // Create vector of systematics
  std::vector<const ISyst*> systs;
  systs.push_back(&kNusOscParamSysts);
  systs.push_back(&kNusNDBeamSysts);
  systs.push_back(&kNusNDBirksSyst);
  systs.push_back(&kNusNDCalibFlatSyst);
  systs.push_back(&kNusNDCalibSlopeXSyst);
  systs.push_back(&kNusNDCalibSlopeYSyst);
  systs.push_back(&kNusNDGENIESmallSysts);
  systs.push_back(&kNusNDNDRockSyst);
  systs.push_back(&kNusNDNueCCSyst);
  systs.push_back(&kNusNDNumuCCSyst);
  systs.push_back(&kNusMCStatsSyst);
  systs.push_back(&kNusNDContSyst);
  systs.push_back(&kNusNormSyst);

  // Declare vectors for shifted spectra
  std::vector<Spectrum> upsCalE, dnsCalE;
  std::vector<Spectrum> upsNHit, dnsNHit;

  // These are important... without this the entire vector gets copied when 
  // resized and triggers an assertion in the Spectrum copy constructor
  upsCalE.reserve(systs.size());
  dnsCalE.reserve(systs.size());
  upsNHit.reserve(systs.size());
  dnsNHit.reserve(systs.size());

  // Set up the shifted spectra
  for(const ISyst* syst:systs)
  {
    upsCalE.emplace_back(floadnear, kNCAxis,   kNusND, SystShifts(syst, +1), kTuftsWeightCC);
    dnsCalE.emplace_back(floadnear, kNCAxis,   kNusND, SystShifts(syst, -1), kTuftsWeightCC);
    upsNHit.emplace_back(floadnear, kAxisNHit, kNusND, SystShifts(syst, +1), kTuftsWeightCC);
    dnsNHit.emplace_back(floadnear, kAxisNHit, kNusND, SystShifts(syst, -1), kTuftsWeightCC);
  }

  floadnear.Go();
  floadneardata.Go();


  TFile* rootF = new TFile("NDDataMCSystBand.root", "RECREATE");

  TDirectory* tmp = gDirectory;
  TDirectory* saveDir = gDirectory;

  // Save nominal spectra
  saveDir = rootF->mkdir("sMCCalE");
  sMCCalE.SaveTo(saveDir);
  saveDir = rootF->mkdir("sDtCalE");
  sDtCalE.SaveTo(saveDir);
  saveDir = rootF->mkdir("dcCalE");
  dcCalE.SaveTo(saveDir);
  saveDir = rootF->mkdir("sMCNHit");
  sMCNHit.SaveTo(saveDir);
  saveDir = rootF->mkdir("sDtNHit");
  sDtNHit.SaveTo(saveDir);
  saveDir = rootF->mkdir("dcNHit");
  dcNHit.SaveTo(saveDir);

  // Save shifted spectra
  std::string specName = "";
  for(unsigned int i = 0, n = systs.size(); i < n; ++i) {
    specName = "sUpCalE" + systs[i]->ShortName();
    saveDir = rootF->mkdir(specName.c_str());
    upsCalE[i].SaveTo(saveDir);
    specName = "sDnCalE" + systs[i]->ShortName();
    saveDir = rootF->mkdir(specName.c_str());
    dnsCalE[i].SaveTo(saveDir);

    specName = "sUpNHit" + systs[i]->ShortName();
    saveDir = rootF->mkdir(specName.c_str());
    upsNHit[i].SaveTo(saveDir);
    specName = "sDnNHit" + systs[i]->ShortName();
    saveDir = rootF->mkdir(specName.c_str());
    dnsNHit[i].SaveTo(saveDir);
  }

  rootF->Close();
}