sterile_demo.C

Go to the documentation of this file.

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcSterile.h"

#include "CAFAna/Core/Binning.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "NuXAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Loaders.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/OscillatableSpectrum.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Core/SpectrumLoaderBase.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Vars/Vars.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "3FlavorAna/Vars/NumuVars.h"

using namespace ana;

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

  // Strings to file paths
  const std::string fnamenear = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*/*/*_s*0_FA*";
  const std::string fnamefar  = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_nonswap_*_r*00_s??_FA*";
  const std::string fnameswap = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_swap_*_r*00_s??_FA*";
  const std::string fnametau  = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/*/*/*_tau_*_r*00_s??_FA*";
  const std::string fnamedata = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*/*/*_s*0_FA*"; // note that this is the same as the near MC!!!

  // Set up a Loaders object
  Loaders loaders("$NOVA_PROD/mc/FA14-10-28", Loaders::kFHC);
  loaders.SetLoaderPath(fnamenear, caf::kNEARDET, Loaders::kMC,   ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnamefar,  caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kNonSwap);
  loaders.SetLoaderPath(fnameswap, caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kFluxSwap);
  loaders.SetLoaderPath(fnametau,  caf::kFARDET,  Loaders::kMC,   ana::kBeam, Loaders::kTauSwap);
  loaders.SetLoaderPath(fnamedata, caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap);

  // Create HistAxis objects for the extrapolation
  const HistAxis axisNC(
    "Calorimetric Energy (GeV)", Binning::Simple(100, 0., 10.), kCaloE);
  const HistAxis axisNumu(
    "Numu CC Energy Estimator (GeV)", Binning::Simple(100, 0., 10.), kCCE);

  // Create Cuts for NC Selection
  const Cut kND([](const caf::SRProxy* sr)
                {
                  if( !sr->vtx.elastic.IsValid) { return false; }
                  return (fabs(sr->vtx.elastic.vtx.X()) < 160. &&
                          fabs(sr->vtx.elastic.vtx.Y()) < 160. &&
                          sr->vtx.elastic.vtx.Z() > 15. && sr->vtx.elastic.vtx.Z() < 1250.);
                });

  const Cut kFD([](const caf::SRProxy* sr)
                {
                  if(!sr->vtx.elastic.IsValid ) { return false; }
                  return (fabs(sr->vtx.elastic.vtx.X()) < 700. &&
                          fabs(sr->vtx.elastic.vtx.Y()) < 700. &&
                          sr->vtx.elastic.vtx.Z() > 15. && sr->vtx.elastic.vtx.Z() < 5800.);
                });

  const Cut kGoodFracMIP([](const caf::SRProxy* sr)
                         {
                           float fracmip = ((sr->slc.ncalhit > 0)
                                            ? sr->slc.nmiphit/float(sr->slc.ncalhit)
                                            : -1);
                           return (fracmip > 0. && fracmip < 0.5);
                         });

  const Cut kIsNotMuon([](const caf::SRProxy* sr)
                       {
                         if(sr->vtx.elastic.fuzzyk.nshwlid < 1) { return true; }
                         return sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.ismuon < 1;
                       });

  // This is the important one, ultimately
  const Cut kNCPresel = kGoodFracMIP && 
                        kElecID < 0.5 && kLEM < 0.3 && kRemID < 0.7 &&
                        kIsNotMuon;

  // Package the cuts in a slightly nicer way
  const Cut cutNDNumu(kNumuQEND);
  const Cut cutND(kND && kNCPresel);
  const Cut cutFD(kFD && kNCPresel);

  // Setup an NC and a numu decomposition
  CheatDecomp decompNC(
    loaders.GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
    axisNC,   cutND
  );
  CheatDecomp decompNumu(
    loaders.GetLoader(caf::kNEARDET, Loaders::kData, ana::kBeam, Loaders::kNonSwap),
    axisNumu, cutNDNumu
  );

  // Setup the extrapolation
  ModularExtrapSterile extrap = ModularExtrapSterile::NCDisappearance(
    loaders, decompNC, decompNumu, axisNC, axisNumu, cutFD, cutND, cutNDNumu
  );
  /*ModularExtrapSterile extrap = ModularExtrapSterile::TrivialExtrapNC(
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap),
    axisNC, cutFD
  );*/

  // Make a sterile style prediction based on the extrapolation
  PredictionSterile pred(&extrap);
  // Make a standard 3 flavor prediction based on the extrapolation
  PredictionExtrap pred_3flav(&extrap);
  // Make a prediction based on FD MC directly from the loaders
  PredictionNoExtrap base(
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kNonSwap),
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kTauSwap),
    axisNC.label, axisNC.bins, axisNC.var, cutFD
  );

  // Define these spectra directly for individual component comparisons
  Spectrum sNueApp(     loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
                        axisNC, cutFD && !kIsAntiNu && kIsSig);
  Spectrum sAntiNueApp( loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
                        axisNC, cutFD &&  kIsAntiNu && kIsSig);
  Spectrum sNumuApp(    loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
                        axisNC, cutFD && !kIsAntiNu && kIsNumuApp);
  Spectrum sAntiNumuApp(loaders.GetLoader(caf::kFARDET, Loaders::kMC, ana::kBeam, Loaders::kFluxSwap),
                        axisNC, cutFD &&  kIsAntiNu && kIsNumuApp);

  // Load from files
  loaders.Go();

  // Open output file
  TFile* rootOut = new TFile("extrapsterile_demo.root", "RECREATE");
  // Save things in output file
  extrap.SaveTo(rootOut, "NCextrap");

  // A 'calculator' that applies no oscillations
  osc::NoOscillations noosc;

  // Set up oscillation calculator with default values
  osc::OscCalcPMNSOpt* calc = (osc::OscCalcPMNSOpt*)DefaultOscCalc();

  // Set up oscillation calculator that uses default 3 flavor parameters
  osc::OscCalcSterile* calcst = new osc::OscCalcSterile();
  calcst->SetNFlavors(3);
  calcst->SetL(810);
  calcst->SetRho(2.84);
  calcst->SetAngle(1, 2, asin(sqrt(.846))/2);
  calcst->SetAngle(1, 3, asin(sqrt(.084))/2);
  calcst->SetAngle(2, 3, M_PI/4);
  calcst->SetDm(2, 7.53e-5); //Dm21
  calcst->SetDm(3, 2.40e-3 + 7.53e-5); //Dm31 = Dm32 + Dm21
  calcst->SetDelta(1, 3, 0);

  double scale = 18e20;

/* These plots show the fraction of NC events
   that come from neutrinos of a given starting flavor/sign */
  new TCanvas;
  TH1* hncnue         = extrap.NCNueProportion().ToTH1();
  hncnue->SetTitle("Proportion of NC events originating from #nu_{e}");
  hncnue->Draw("hist");
  new TCanvas;
  TH1* hncantinue     = extrap.NCAntiNueProportion().ToTH1();
  hncantinue->SetTitle("Proportion of NC events originating from #bar{#nu}_{e}");
  hncantinue->Draw("hist");
  new TCanvas;
  TH1* hncnumu        = extrap.NCNumuProportion().ToTH1();
  hncnumu->SetTitle("Proportion of NC events originating from #nu_{#mu}");
  hncnumu->Draw("hist");
  new TCanvas;
  TH1* hncantinumu    = extrap.NCAntiNumuProportion().ToTH1();
  hncantinumu->SetTitle("Proportion of NC events originating from #bar{#nu}_{#mu}");
  hncantinumu->Draw("hist");

/* This plot is a check that all components add to 1 */
/*  new TCanvas;
  TH1* hTotal = (TH1*)hncnue->Clone();
  hTotal->Add(hncantinue);
  hTotal->Add(hncnumu);
  hTotal->Add(hncantinumu);
  hTotal->SetTitle("Sum of all Proportions");
  hTotal->Draw("hist");
*/

/* These plots show each extrapolated component */
  new TCanvas;
  extrap.OscNCComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.NueSurvComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.AntiNueSurvComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.NumuSurvComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.AntiNumuSurvComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.NueAppComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.AntiNueAppComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.NumuAppComponent().ToTH2(scale)->Draw("colz");
  new TCanvas;
  extrap.AntiNumuAppComponent().ToTH2(scale)->Draw("colz");

/* These plots are ratios of extrapolated components to FD MC
   The extrapolated components are collapsed to the extrapolated variable,
   i.e., calorimetric energy
   The FD MC are spectra that select a given component by truth
   The same cuts are applied to both spectra, so the ratio should be 1 */
/*  new TCanvas;
  TH1* hNCSurv = extrap.OscNCComponent().Unoscillated().ToTH1(scale);
  Spectrum sNCSurv = base.PredictComponent(&noosc, Flavors::kAll, Current::kNC, Sign::kBoth);
  hNCSurv->Divide(sNCSurv.ToTH1(scale));
  hNCSurv->SetTitle("Ratio of Extrapolation to FD MC: NC Surv");
  hNCSurv->Draw("hist");
  new TCanvas;
  TH1* hNueSurv = extrap.NueSurvComponent().Unoscillated().ToTH1(scale);
  Spectrum sNueSurv = base.PredictComponent(&noosc, Flavors::kNuEToNuE, Current::kCC, Sign::kNu);
  hNueSurv->Divide(sNueSurv.ToTH1(scale));
  hNueSurv->SetTitle("Ratio of Extrapolation to FD MC: #nu_{e} Surv");
  hNueSurv->Draw("hist");
  new TCanvas;
  TH1* hAntiNueSurv = extrap.AntiNueSurvComponent().Unoscillated().ToTH1(scale);
  Spectrum sAntiNueSurv = base.PredictComponent(&noosc, Flavors::kNuEToNuE, Current::kCC, Sign::kAntiNu);
  hAntiNueSurv->Divide(sAntiNueSurv.ToTH1(scale));
  hAntiNueSurv->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{e} Surv");
  hAntiNueSurv->Draw("hist");
  new TCanvas;
  TH1* hNumuSurv = extrap.NumuSurvComponent().Unoscillated().ToTH1(scale);
  Spectrum sNumuSurv = base.PredictComponent(&noosc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kNu);
  hNumuSurv->Divide(sNumuSurv.ToTH1(scale));
  hNumuSurv->SetTitle("Ratio of Extrapolation to FD MC: #nu_{#mu} Surv");
  hNumuSurv->Draw("hist");
  new TCanvas;
  TH1* hAntiNumuSurv = extrap.AntiNumuSurvComponent().Unoscillated().ToTH1(scale);
  Spectrum sAntiNumuSurv = base.PredictComponent(&noosc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kAntiNu);
  hAntiNumuSurv->Divide(sAntiNumuSurv.ToTH1(scale));
  hAntiNumuSurv->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{#mu} Surv");
  hAntiNumuSurv->Draw("hist");
  new TCanvas;
  TH1* hNueApp = extrap.NueAppComponent().Unoscillated().ToTH1(scale);
  hNueApp->Divide(sNueApp.ToTH1(scale));
  hNueApp->SetTitle("Ratio of Extrapolation to FD MC: #nu_{e} App");
  hNueApp->Draw("hist");
  new TCanvas;
  TH1* hAntiNueApp = extrap.AntiNueAppComponent().Unoscillated().ToTH1(scale);
  hAntiNueApp->Divide(sAntiNueApp.ToTH1(scale));
  hAntiNueApp->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{e} App");
  hAntiNueApp->Draw("hist");
  new TCanvas;
  TH1* hNumuApp = extrap.NumuAppComponent().Unoscillated().ToTH1(scale);
  hNumuApp->Divide(sNumuApp.ToTH1(scale));
  hNumuApp->SetTitle("Ratio of Extrapolation to FD MC: #nu_{#mu} App");
  hNumuApp->Draw("hist");
  new TCanvas;
  TH1* hAntiNumuApp = extrap.AntiNumuAppComponent().Unoscillated().ToTH1(scale);
  hAntiNumuApp->Divide(sAntiNumuApp.ToTH1(scale));
  hAntiNumuApp->SetTitle("Ratio of Extrapolation to FD MC: #bar{#nu}_{#mu} App");
  hAntiNumuApp->Draw("hist");
*/

/* This plot takes the extrapolation and base MC,
   applies null oscillations, and forms a ratio.
   If the ND Data and MC inputs to the extrapolation are identical,
   or if the extrapolation is the TrivialExtrapNC,
   the ratio should be 1. */
  new TCanvas;
  Spectrum predNoOsc = pred.PredictUnoscillated();
  Spectrum baseNoOsc = base.PredictUnoscillated();
  TH1* hratioNull = predNoOsc.ToTH1(scale);
  hratioNull->Divide(baseNoOsc.ToTH1(scale));
  hratioNull->SetTitle("Ratio of Extrapolation to FD MC (No Oscillations)");
  hratioNull->Draw("hist");

/* This plot takes the sterile prediction and standard 3 flavor prediction,
   applies the same oscillations, and forms a ratio.
   Since both predictions use the same extrapolation, the ratio should be 1. */
  new TCanvas;
  Spectrum predOsc = pred.Predict(calcst);
  Spectrum pred3FlavOsc = base.Predict(calc);
  TH1* hratioOsc = predOsc.ToTH1(scale);
  hratioOsc->Divide(pred3FlavOsc.ToTH1(scale));
  hratioOsc->SetTitle("Ratio of Sterile to 3 Flavor Prediction (Same Oscillations)");
  hratioOsc->Draw("hist");

  rootOut->Close();
}