MichelDecomp.cxx

Go to the documentation of this file.

#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/LoadFromRegistry.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Utilities.h"
#include "CAFAna/Vars/Vars.h"

#include "Utilities/func/MathUtil.h"

#include <cassert>

#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TCanvas.h"
#include "TProfile.h"
#include "TAxis.h"
#include "TMatrixD.h"
#include "TMath.h"
#include "TVectorD.h"
#include "TDirectory.h"
#include "TObjString.h"
#include "TRandom.h"

namespace ana
{
  REGISTER_LOADFROM("MichelDecomp", IDecomp, MichelDecomp);

  //----------------------------------------------------------------------


  const Binning MEBinning = Binning::Simple(3,0,3);
  const HistAxis MEAxis("N_{Michels}",MEBinning,kNMichels);

  //----------------------------------------------------------------------
  //----------------------------------------------------------------------
  MichelDecomp::MichelDecomp(SpectrumLoaderBase& NDMCloader,
                             SpectrumLoaderBase& NDDataloader,
                             const HistAxis&     EAxis,
                             const Cut&          cut,
                             const IDecomp*      nueEstimate,
                             const SystShifts&   shiftMC,
                             const SystShifts&   shiftData,
                             const Var&          wei)
    : fNC         (NDMCloader, EAxis,
                   cut && kIsNC && !kIsAntiNu, shiftMC, wei),
      fAntiNC     (NDMCloader, EAxis,
                   cut && kIsNC && kIsAntiNu, shiftMC, wei),
      fNumu       (NDMCloader, EAxis,
                   cut && kIsNumuCC && !kIsAntiNu, shiftMC, wei),
      fAntiNumu   (NDMCloader, EAxis,
                   cut && kIsNumuCC && kIsAntiNu, shiftMC, wei),
      fNue        (NDMCloader, EAxis,
                   cut && kIsBeamNue && !kIsAntiNu, shiftMC, wei),
      fAntiNue    (NDMCloader, EAxis,
                   cut && kIsBeamNue && kIsAntiNu, shiftMC, wei),
      fData       (NDDataloader, EAxis, cut, shiftData, wei),
      fMichelNC   (NDMCloader, EAxis, MEAxis,
                   cut && kIsNC, shiftMC, wei),
      fMichelNumu (NDMCloader, EAxis, MEAxis,
                   cut && kIsNumuCC, shiftMC, wei),
      fMichelNue  (NDMCloader, EAxis, MEAxis,
                   cut && kIsBeamNue, shiftMC, wei),
      fMichelData (NDDataloader, EAxis, MEAxis,
                   cut, shiftData, wei),
      fNueEstimate(nueEstimate),
      fOwnNueEstimate(false),
      isDecomposed(false)
  {}

  MichelDecomp::MichelDecomp(SpectrumLoaderBase& NDMCloader,
                             SpectrumLoaderBase& NDDataloader,
                             std::string         label,
                             const Binning&      bins,
                             const Var&          var,
                             const Cut&          cut,
                             const IDecomp*      nueEstimate,
                             const SystShifts&   shiftMC,
                             const SystShifts&   shiftData,
                             const Var&          wei)
    : MichelDecomp(NDMCloader, NDDataloader,
                   HistAxis(label, bins, var),
                   cut, nueEstimate, shiftMC, shiftData, wei)
  {}

  MichelDecomp::MichelDecomp(Loaders&            loaders,
                             const HistAxis&     axis,
                             const Cut&          cut,
                             const IDecomp*      nueEstimate,
                             const SystShifts&   shiftMC,
                             const SystShifts&   shiftData,
                             const Var&          wei)
    : MichelDecomp(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
                   loaders.GetLoader(caf::kNEARDET, Loaders::kData),
                   axis, cut, nueEstimate, shiftMC, shiftData, wei)
  {}

  MichelDecomp::~MichelDecomp()
  {
    if(fOwnNueEstimate) delete fNueEstimate;
  }

  /// Calls ComputeScaleFactor for every process / E Bin, saves scale factors
  void MichelDecomp::Decompose() const{
    std::vector<double> numuScales;
    std::vector<double> ncScales;
    std::vector<double> nueScales;

    TH1* EIdxLimitHist = fNumu.ToTH1(1e20);
    for (int EIdx = 1; EIdx <= EIdxLimitHist->GetNbinsX(); EIdx++){
      // Get appropriate nuescale from input decomp

      // Stores sums of our histograms, to reduce computing in functions
      MDCMPHelper hSums;
      hSums.fSumData   = GetSum(fMichelData, EIdx);
      hSums.fSumNumu   = GetSum(fMichelNumu, EIdx);
      hSums.fSumNC     = GetSum(fMichelNC,   EIdx);
      hSums.fSumNue    = GetSum(fMichelNue,  EIdx);
      hSums.fNueScale  = GetNueScale(EIdx);
      hSums.fNumuScale = GetNumuScale(EIdx);
      hSums.fNCScale   = GetNCScale(EIdx);

      // Typical nue slices have about 1/8 SlcME's per slice
      // So, if there are more than about 80 events in a bin, we have 
      // at least 10 SlcME's in that bin, which is enough to run MDCMP
      // If we have fewer, proportionally split NC's and CC's

      // Instead of cheacking total number of events in a bin, look at 
      // number of data events seen, less the nue component predicted
      /*if (EIdx == 1 || EIdx == 2 || EIdx == 8 || EIdx == 9 || EIdx == 10 ||
          EIdx ==11 || EIdx ==16 || EIdx ==17 || EIdx ==18 || EIdx == 19 ||
          EIdx ==20 || EIdx ==21 || EIdx > 23)*/
      if ((GetNumuMCContent(EIdx) < .2) || (GetNCMCContent(EIdx) == .0))
        //SplitByProp(hSums, numuScales, ncScales, nueScales);
        SplitByNue(hSums, numuScales, ncScales, nueScales);
      else
        SplitByMichels(EIdx, hSums, numuScales, ncScales, nueScales);
    }
    delete EIdxLimitHist;

    fFitResult.numuScales = numuScales;
    fFitResult.ncScales   = ncScales;
    fFitResult.nueScales  = nueScales;
    isDecomposed = true;
  }

  void MichelDecomp::SplitByMichels(int EIdx, MDCMPHelper hSums,
                                    std::vector<double> &numuScales,
                                    std::vector<double> &ncScales,
                                    std::vector<double> &nueScales) const
  {
      // We have enough data in this bin to do MichelDecomp
      double numuScale = ComputeScaleFactor(EIdx, hSums);
      double ncScale   =
        (hSums.fSumData - numuScale*hSums.fSumNumu -
         hSums.fNueScale*hSums.fSumNue) / hSums.fSumNC;
      numuScales.push_back(numuScale);
      ncScales.push_back(ncScale);
      nueScales.push_back(hSums.fNueScale);
  }

  void MichelDecomp::SplitByProp(MDCMPHelper hSums,
                                 std::vector<double> &numuScales,
                                 std::vector<double> &ncScales,
                                 std::vector<double> &nueScales) const
  {
    // Return MC if there are no events in the MC
    if (hSums.fSumNumu + hSums.fSumNC + hSums.fSumNue == 0){
      numuScales.push_back(1);
      ncScales.push_back(1);
      nueScales.push_back(1);
      return;
    }
    double DataMC = hSums.fSumData /
                    (hSums.fSumNue + hSums.fSumNumu + hSums.fSumNC);
      numuScales.push_back(DataMC);
      ncScales.push_back(DataMC);
      nueScales.push_back(DataMC);
  }
  
  void MichelDecomp::SplitByNue(MDCMPHelper hSums,
                                std::vector<double> &numuScales,
                                std::vector<double> &ncScales,
                                std::vector<double> &nueScales) const
  {
      numuScales.push_back(hSums.fNumuScale);
      ncScales.push_back(hSums.fNCScale);
      nueScales.push_back(hSums.fNueScale);
  }

  /// Nue scale from input decomp
  double MichelDecomp::GetNueScale(int EIdx) const {
    Spectrum nue = fNueEstimate->NueComponent() +
                   fNueEstimate->AntiNueComponent();
    Ratio nuerat = nue / (fNue+fAntiNue);
    TH1* hnuerat = nuerat.ToTH1();
    double scale = hnuerat->GetBinContent(EIdx);
    delete hnuerat;
    return scale;
  }
  
  /// Numu scale from input decomp
  double MichelDecomp::GetNumuScale(int EIdx) const {
    Spectrum numu = fNueEstimate->NumuComponent() +
                    fNueEstimate->AntiNumuComponent();
    Ratio numurat = numu / (fNumu+fAntiNumu);
    TH1* hnumurat = numurat.ToTH1();
    double scale = hnumurat->GetBinContent(EIdx);
    delete hnumurat;
    return scale;
  }
  
  /// NC scale from input decomp
  double MichelDecomp::GetNCScale(int EIdx) const {
    Spectrum nc = fNueEstimate->NCTotalComponent();
    Ratio ncrat = nc / (fNC+fAntiNC);
    TH1* hncrat = ncrat.ToTH1();
    double scale = hncrat->GetBinContent(EIdx);
    delete hncrat;
    return scale;
  }

  /// Getting numu relative MC content
  double MichelDecomp::GetNumuMCContent(int EIdx) const {
    Spectrum numu = fNumu+fAntiNumu;
    Ratio numurat = numu / (numu+fNC+fAntiNC+fNue+fAntiNue);
    TH1* hnumurat = numurat.ToTH1();
    double numucont = hnumurat->GetBinContent(EIdx);
    delete hnumurat;
    return numucont;
  }

  double MichelDecomp::GetNCMCContent(int EIdx) const {
    Spectrum nc = fNC+fAntiNC;
    Ratio ncrat = nc / (nc+fNumu+fAntiNumu+fNue+fAntiNue);
    TH1* hncrat = ncrat.ToTH1();
    double nccont = hncrat->GetBinContent(EIdx);
    delete hncrat;
    return nccont;
  }

  ///////////////////////////////////////////////////////////////////////////
  /// Starting code to run through decomp algebra
  ///////////////////////////////////////////////////////////////////////////
  /// Compute scale factors
  double MichelDecomp::ComputeScaleFactor(int EIdx, MDCMPHelper hSums) const {
    double dx = 0.01;

    // Gives the bounds for the scale factor after each iteration
    double MinScale = 0;
    double MaxScale =  ComputeMaxScale(hSums);
    double scale = (MaxScale-MinScale) / 2;

    // Continue for 10 iterations
    for (double N = 0; N < 10; N++){
      // Adjust dx to be within the current allowed scale factors
      if ((MaxScale - MinScale) < 2*dx)
        dx = (MaxScale - MinScale)/4;
      double likeplus  = MDCMPLogLikelihood(scale+dx, EIdx, hSums);
      double likeminus = MDCMPLogLikelihood(scale-dx, EIdx, hSums);
      MaxScale = (likeplus > likeminus) ? MaxScale : scale;
      MinScale = (likeplus > likeminus) ? scale : MinScale;
      scale = (MaxScale+MinScale)/2;
    }
    return scale;
  }


  double MichelDecomp::ComputeMaxScale(MDCMPHelper hSums) const{
    assert(hSums.fSumNumu > 0 && "Don't have enough stats.  Bailing.");
    double MaxScale = (hSums.fSumData - hSums.fNueScale*hSums.fSumNue) /
                       hSums.fSumNumu;
    return MaxScale;
  }

  /// Calculates the expected MC events, for a given scale factor
  double MichelDecomp::CalculateMCMean(double h,
                                       double Numu,
                                       double NC,
                                       double Nue,
                                       MDCMPHelper hSums) const{
    assert(hSums.fSumNC > 0 && "Don't have enough stats.  Bailing.");

    double mean =
      h*(Numu - hSums.fSumNumu / hSums.fSumNC * NC) +
      hSums.fSumData / hSums.fSumNC * NC -
      hSums.fNueScale*hSums.fSumNue / hSums.fSumNC * NC + hSums.fNueScale*Nue;
    return mean;
  }

  /// Calculates the LL for a given scale, based on the data
  double MichelDecomp::MDCMPLogLikelihood(double h, int EIdx,
                                          MDCMPHelper hSums) const{
    double LL = 0;

    TH1* mIdxLimitHist = fMichelNumu.ToTH2(1e20);
    for (int mIdx = 1; mIdx <= mIdxLimitHist->GetNbinsY(); mIdx++){
      double Data = GetTemplateContent(fMichelData, EIdx, mIdx);
      double Numu = GetTemplateContent(fMichelNumu, EIdx, mIdx);
      double NC   = GetTemplateContent(fMichelNC,   EIdx, mIdx);
      double Nue  = GetTemplateContent(fMichelNue,  EIdx, mIdx);
      double mean = CalculateMCMean(h, Numu, NC, Nue, hSums);
      if (mean == 0 && Data == 0)
        continue;
      if (Data == 0){
        LL += Data - mean;
        continue;
      }
      if (mean == 0)
        mean = 1e-10;
      LL += Data - Data * TMath::Log(Data / mean) - mean;
    }
    delete mIdxLimitHist;

    return LL;
  }

  /// Sums given spectrum over NME bins, for a given E Bin
  double MichelDecomp::GetSum(ReweightableSpectrum spect, int EIdx) const{
    TH2* hist = spect.ToTH2(fData.POT());
    double sum = hist->ProjectionX()->GetBinContent(EIdx);
    delete hist;
    return sum;
  }

  /// Gets bin content in template histogram, for a given E / NME bin.
  double MichelDecomp::GetTemplateContent(ReweightableSpectrum spec,
                                          int EIdx, int mIdx) const{
    TH2* hist = spec.ToTH2(fData.POT());
    double content = hist->GetBinContent(EIdx, mIdx);
    delete hist;
    return content;
  }
  ///////////////////////////////////////////////////////////////////////////
  /// Ending code to run through decomp algebra
  ///////////////////////////////////////////////////////////////////////////

  // Inputs a std::vector of scale factors, and a Comp spectrum,
  // And outputs decomp result with components scaled by MDCMP fit
  Spectrum MichelDecomp::MCToDCMPComp(std::vector<double> scales,
                                      Spectrum MCComp) const{
    Eigen::ArrayXd scaledMC = MCComp.GetEigen(fData.POT());
    for (int EIdx = 0; EIdx < scaledMC.size()-2; EIdx++){
      // Bin content in new histogram is old content * fit scale
      int    binidx  = EIdx+1; // Root counts bins starting at 1
      double bincont = scales[EIdx]*scaledMC[binidx];
      scaledMC[binidx] = bincont;
    }
    return Spectrum(std::move(scaledMC),
                    HistAxis(fMichelData.GetLabels(), fMichelData.GetBinnings()),
                    fData.POT(), fData.Livetime());
  }

  //----------------------------------------------------------------------
  void MichelDecomp::SaveCompPlots(TDirectory* dir, std::string comp,
                                   Spectrum dcmp, Spectrum mc) const
  {
    double pot = fData.POT();
    dir->cd();
    TH1* hdcmp = dcmp.ToTH1(pot);
    hdcmp->Write((comp+"_micheldecomp").c_str());
    TH1* hmc   = mc.ToTH1(pot);
    hmc->Write((comp+"_mc").c_str());

    TH1* hratio = (TH1*)hdcmp->Clone("hratio");
    hratio->Divide(hmc);
    hratio->Write((comp+"_ratio").c_str());

    delete hdcmp; delete hmc; delete hratio;
  }

  //----------------------------------------------------------------------
  void MichelDecomp::SaveTempPlots(TDirectory* dir) const
  {
    double pot = fData.POT();
    dir->cd();

    TH2* tempnumu = fMichelNumu.ToTH2(pot);
    TH2* tempnc   = fMichelNC.  ToTH2(pot);
    TH2* tempnue  = fMichelNue. ToTH2(pot);
    TH2* tempdata = fMichelData.ToTH2(pot);

    tempnumu->Write("NumuTemplate");
    tempnc  ->Write("NCTemplate");
    tempnue ->Write("NueTemplate");
    tempdata->Write("DataTemplate");

    delete tempnumu; delete tempnc; delete tempnue; delete tempdata;
  }

  //----------------------------------------------------------------------
  void MichelDecomp::SavePlots(TDirectory* dir) const
  {
    TDirectory* tmp = gDirectory;
    dir->cd();

    SaveCompPlots(dir->mkdir("numu"), "numu",
                             NumuComponent(), fNumu);
    SaveCompPlots(dir->mkdir("antinumu"), "antinumu",
                             AntiNumuComponent(), fAntiNumu);
    SaveCompPlots(dir->mkdir("nc"), "nc",
                             NCComponent(), fNC);
    SaveCompPlots(dir->mkdir("antinc"), "antinc",
                             NCAntiComponent(), fAntiNC);
    SaveCompPlots(dir->mkdir("nue"), "nue",
                             NueComponent(), fNue);
    SaveCompPlots(dir->mkdir("antinue"), "antinue",
                             AntiNueComponent(), fAntiNue);

    SaveTempPlots(dir->mkdir("templates"));

    tmp->cd();
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::NueComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.nueScales, fNue);
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::AntiNueComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.nueScales, fAntiNue);
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::NumuComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.numuScales, fNumu);
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::AntiNumuComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.numuScales, fAntiNumu);
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::NCTotalComponent() const {
    if (!isDecomposed)
      Decompose();
    return NCComponent()+NCAntiComponent();
  }

  Spectrum MichelDecomp::NCComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.ncScales, fNC);
  }

  Spectrum MichelDecomp::NCAntiComponent() const {
    if (!isDecomposed)
      Decompose();
    return MCToDCMPComp(fFitResult.ncScales, fAntiNC);
  }

  // And now make the MC distributions publicly accessible
  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_NueComponent() const {
    return fNue;
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_AntiNueComponent() const {
    return fAntiNue;
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_NumuComponent() const {
    return fNumu;
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_AntiNumuComponent() const {
    return fAntiNumu;
  }

  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_NCTotalComponent() const {
    return fNC+fAntiNC;
  }
  
  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_NCComponent() const {
    return fNC;
  }
  
  //----------------------------------------------------------------------
  Spectrum MichelDecomp::MC_NCAntiComponent() const {
    return fAntiNC;
  }
  
  //----------------------------------------------------------------------
  Spectrum MichelDecomp::Data_Component() const {
    return fData;
  }

  //----------------------------------------------------------------------
  void MichelDecomp::SaveTo(TDirectory* dir, const std::string& name) const{
    TDirectory* tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("MichelDecomp").Write("type");

    fNC.SaveTo(          dir, "NC");
    fAntiNC.SaveTo(      dir, "AntiNC");
    fNumu.SaveTo(        dir, "Numu");
    fAntiNumu.SaveTo(    dir, "AntiNumu");
    fNue.SaveTo(         dir, "Nue");
    fAntiNue.SaveTo(     dir, "AntiNue");
    fData.SaveTo(     dir, "Data");
    fMichelNC.SaveTo(    dir, "TempNC");
    fMichelNumu.SaveTo(  dir, "TempNumu");
    fMichelNue.SaveTo(   dir, "TempNue");
    fMichelData.SaveTo(  dir, "TempData");
    fNueEstimate->SaveTo(dir, "NueEstimate");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  std::unique_ptr<MichelDecomp> MichelDecomp::LoadFrom(TDirectory* dir, const std::string& name){
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    Spectrum* nc       = ana::LoadFrom<Spectrum>
      (dir, "NC").release();
    Spectrum* antinc   = ana::LoadFrom<Spectrum>
      (dir, "AntiNC").release();
    Spectrum* numu     = ana::LoadFrom<Spectrum>
      (dir, "Numu").release();
    Spectrum* antinumu = ana::LoadFrom<Spectrum>
      (dir, "AntiNumu").release();
    Spectrum* nue      = ana::LoadFrom<Spectrum>
      (dir, "Nue").release();
    Spectrum* antinue  = ana::LoadFrom<Spectrum>
      (dir, "AntiNue").release();
    Spectrum* data  = ana::LoadFrom<Spectrum>
      (dir, "Data").release();
    ReweightableSpectrum* tempnc   = ana::LoadFrom<ReweightableSpectrum>
      (dir, "TempNC").release();
    ReweightableSpectrum* tempnumu = ana::LoadFrom<ReweightableSpectrum>
      (dir, "TempNumu").release();
    ReweightableSpectrum* tempnue  = ana::LoadFrom<ReweightableSpectrum>
      (dir, "TempNue").release();
    ReweightableSpectrum* tempdata = ana::LoadFrom<ReweightableSpectrum>
      (dir, "TempData").release();
    const IDecomp* nueEst = ana::LoadFrom<IDecomp>
      (dir, "NueEstimate").release();

    MichelDecomp* ret =
      new MichelDecomp(*nc, *antinc,
                       *numu, *antinumu, *nue, *antinue, *data,
                       *tempnc, *tempnumu, *tempnue, *tempdata, nueEst);

    ret->fOwnNueEstimate = true;

    delete dir;

    return std::unique_ptr<MichelDecomp>(ret);
  }
}

//  LocalWords:  EIdx