BENDecomp.cxx

Go to the documentation of this file.

#include "3FlavorAna/Decomp/BENDecomp.h"

#include "StandardRecord/StandardRecord.h"
#include "StandardRecord/Proxy/SRProxy.h"

#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Cuts/BeamCuts.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Core/FileReducer.h"
#include "CAFAna/Core/LoadFromRegistry.h"

#include "3FlavorAna/Cuts/NueCuts2017.h"
#include "3FlavorAna/Cuts/NueCuts2018.h"
#include "3FlavorAna/Cuts/NueCuts2020.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/NumuCuts2017.h"
#include "3FlavorAna/Cuts/NumuCuts2018.h"
#include "3FlavorAna/Cuts/NumuCuts2020.h"
#include "3FlavorAna/Cuts/BeamNueCuts.h"
#include "3FlavorAna/Vars/HistAxes.h"

#include "Utilities/rootlogon.C"

#include "TDirectory.h"
#include "TObjString.h"
#include "TH3.h"
#include "TH2.h"
#include "TGraph.h"
#include "TLegend.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TVectorD.h"
#include "TLine.h"

#include <cassert>

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

  BENDecomp::BENDecomp(SpectrumLoaderBase& loaderMC,
                       SpectrumLoaderBase& loaderData,
                       const HistAxis& axis,
                       const Cut& cut,
                       const SystShifts& shiftMC,
                       const SystShifts& shiftData,
                       const Var& wei,
                       const HistAxis & axisNumuPi,
                       const Cut& cutNumuPi,
                       const double kaonNormalization
                       )
    : fData    (loaderData, axis, cut,                              shiftData, wei),
      fNCTot   (loaderMC,   axis, cut && kIsNC,                       shiftMC, wei),
      fNue     (loaderMC,   axis, cut && kIsBeamNue&&!kIsAntiNu,      shiftMC, wei), 
      fAntiNue (loaderMC,   axis, cut && kIsBeamNue&& kIsAntiNu,      shiftMC, wei),
      fNumu    (loaderMC,   axis, cut && kIsNumuCC &&!kIsAntiNu,      shiftMC, wei),
      fAntiNumu(loaderMC,   axis, cut && kIsNumuCC && kIsAntiNu,      shiftMC, wei),
      fTotal   (loaderMC,   axis, cut && kHasNu,                      shiftMC, wei),
      fNotNue  (loaderMC,   axis, cut && kHasNu && !(kIsBeamNue && !kIsAntiNu), shiftMC, wei),

      fNueCCFromPi    (loaderMC, axis, cut && kIsNueCCFromPiDk,                           shiftMC, wei),
      fNueCCFromKa    (loaderMC, axis, cut && kIsNueCCFromKaDk,                           shiftMC, wei),
      fNueCCFromOther (loaderMC, axis, cut && kIsNueCCNotFromPiDk && kIsNueCCNotFromKaDk, shiftMC, wei),
      fNueCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axis, cut && kIsNueCCFromPiDk,       shiftMC, wei, Spectrum::kSparse),

      fNumuSelData      (loaderData, axisNumuPi, cutNumuPi,                                              shiftData, wei),
      fNumuSelBkg         (loaderMC, axisNumuPi, cutNumuPi && (!kIsNumuCC || kIsAntiNu),         shiftMC, wei),
      fNumuSelCCFromPi    (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk,                           shiftMC, wei),
      fNumuSelCCFromKa    (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromKaDk,                           shiftMC, wei),
      fNumuSelCCFromOther (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk,shiftMC, wei),
      fNumuSelCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis,  axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk,      shiftMC, wei, Spectrum::kSparse),

      fNumuUncontainData(Spectrum::Uninitialized()),
      fNumuUncontainBkg(Spectrum::Uninitialized()),
      fNumuUncontainCCFromPi(Spectrum::Uninitialized()),
      fNumuUncontainCCFromKa(Spectrum::Uninitialized()),
      fNumuUncontainCCFromOther(Spectrum::Uninitialized()),

      fNueCCFromPiEstim(Spectrum::Uninitialized()),
      fNueCCFromKaEstim(Spectrum::Uninitialized()),

      fIsDecomposed(false),
      fNumuPiWeights(0),
      fIsFixedKaScale(true),   
      fKaonNormalization(kaonNormalization)
  {}

  BENDecomp::BENDecomp(Loaders& loaders,
                       const HistAxis& axis,
                       const Cut& cut,
                       const SystShifts& shiftMC,
                       const SystShifts& shiftData,
                       const Var& wei,
                       const HistAxis & axisNumuPi,
                       const Cut& cutNumuPi,
                       const double kaonNormalization
                       )
    : BENDecomp(loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
               loaders.GetLoader(caf::kNEARDET, Loaders::kData),
               axis, cut, shiftMC, shiftData, wei,
               axisNumuPi, cutNumuPi, kaonNormalization)
  {}

  BENDecomp::BENDecomp (EBENCalcKa kaopt, 
                        SpectrumLoaderBase& loaderMC,
                        SpectrumLoaderBase& loaderData,
                        const HistAxis& axis,
                        const Cut& cut,
                        const SystShifts& shiftMC,
                        const SystShifts& shiftData,
                        const Var& wei,
                        const HistAxis& axisNumuPi,
                        const Cut& cutNumuPi,
                        const HistAxis& axisNumuKa,
                        const Cut& cutNumuKa
                        )
    : fData    (loaderData, axis, cut,                              shiftData, wei),
      fNCTot   (loaderMC,   axis, cut && kIsNC,                       shiftMC, wei),
      fNue     (loaderMC,   axis, cut && kIsBeamNue&&!kIsAntiNu,      shiftMC, wei), 
      fAntiNue (loaderMC,   axis, cut && kIsBeamNue&& kIsAntiNu,      shiftMC, wei),
      fNumu    (loaderMC,   axis, cut && kIsNumuCC &&!kIsAntiNu,      shiftMC, wei),
      fAntiNumu(loaderMC,   axis, cut && kIsNumuCC && kIsAntiNu,      shiftMC, wei),
      fTotal   (loaderMC,   axis, cut && kHasNu,                      shiftMC, wei),
      fNotNue  (loaderMC,   axis, cut && !(kIsBeamNue && !kIsAntiNu) && kHasNu, shiftMC, wei),
      
      fNueCCFromPi    (loaderMC, axis, cut && kIsNueCCFromPiDk,                           shiftMC, wei),
      fNueCCFromKa    (loaderMC, axis, cut && kIsNueCCFromKaDk,                           shiftMC, wei),
      fNueCCFromOther (loaderMC, axis, cut && kIsNueCCNotFromPiDk && kIsNueCCNotFromKaDk, shiftMC, wei),
      fNueCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis, axis, cut && kIsNueCCFromPiDk,       shiftMC, wei, Spectrum::kSparse),
      
      fNumuSelData      (loaderData, axisNumuPi, cutNumuPi,                                              shiftData, wei),
      fNumuSelBkg         (loaderMC, axisNumuPi, cutNumuPi && (!kIsNumuCC || kIsAntiNu),         shiftMC, wei),
      fNumuSelCCFromPi    (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk,                           shiftMC, wei),
      fNumuSelCCFromKa    (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCFromKaDk,                           shiftMC, wei),
      fNumuSelCCFromOther (loaderMC, axisNumuPi, cutNumuPi && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk,shiftMC, wei),
      fNumuSelCCFromPiPtPz(loaderMC, ktpzAxis, ktpTAxis,  axisNumuPi, cutNumuPi && kIsNumuCCFromPiDk,      shiftMC, wei, Spectrum::kSparse),
      
      fNumuUncontainData(loaderData, axisNumuKa, cutNumuKa, shiftData, wei),
      fNumuUncontainBkg (loaderMC, axisNumuKa, cutNumuKa && (!kIsNumuCC || kIsAntiNu),           shiftMC, wei),
      fNumuUncontainCCFromPi(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCFromPiDk,               shiftMC, wei),
      fNumuUncontainCCFromKa(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCFromKaDk,               shiftMC, wei),
      fNumuUncontainCCFromOther(loaderMC, axisNumuKa, cutNumuKa && kIsNumuCCNotFromPiDk && kIsNumuCCNotFromKaDk, shiftMC, wei),
      
      fNueCCFromPiEstim(Spectrum::Uninitialized()),
      fNueCCFromKaEstim(Spectrum::Uninitialized()),
      
      fIsDecomposed(false),
      fNumuPiWeights(0),
      fIsFixedKaScale(false),   
      fKaonNormalization(1)
  {}

  BENDecomp::BENDecomp (EBENCalcKa kaopt,
                        Loaders& loaders,
                        const HistAxis& axis,
                        const Cut& cut,
                        const SystShifts& shiftMC,
                        const SystShifts& shiftData,
                        const Var& wei,
                        const HistAxis& axisNumuPi,
                        const Cut& cutNumuPi,
                        const HistAxis& axisNumuKa,
                        const Cut& cutNumuKa
                        )
    : BENDecomp(kaopt ,
                loaders.GetLoader(caf::kNEARDET, Loaders::kMC),
                loaders.GetLoader(caf::kNEARDET, Loaders::kData),
                axis, cut, shiftMC, shiftData, wei,
                axisNumuPi, cutNumuPi, axisNumuKa, cutNumuKa)
  {}


//--------------------- 38 North -------------------

  Spectrum BENDecomp::NueEstimate() const {
    if(!fIsDecomposed)
      Decompose();

    double thePOT = fNueCCFromPiEstim.POT();

    Eigen::ArrayXd aData = fData.GetEigen(thePOT);
    Eigen::ArrayXd aNue1D = (fNueCCFromPiEstim + fNueCCFromKaEstim + fNueCCFromOther).GetEigen(thePOT);
    for(int i = 0; i < aNue1D.size(); ++i){
      // Nue estimate always between 0 and data
      aNue1D[i] = std::min(std::max(0., aNue1D[i]), aData[i]);
    }

    return Spectrum(std::move(aNue1D),
                    HistAxis(fData.GetLabels(), fData.GetBinnings()),
                    thePOT, 0);
  }

  void BENDecomp::Decompose() const{
    MakeWeightsNumuFromPion();
    MakeWeightsNumuFromKaon();

    NueEstimateFromPi();
    NueEstimateFromKa();

    fIsDecomposed = true;
  }

  BENDecomp::~BENDecomp()
  {
    delete fNumuPiWeights;
  }

 void BENDecomp::NueEstimateFromPi() const {
    DontAddDirectory guard;

    auto thePOT = fNueCCFromPiPtPz.POT();
    auto hNue3D = fNueCCFromPiPtPz.ToTH3(thePOT); //pt pz var
    int nbinsx= hNue3D->GetNbinsX();
    int nbinsy= hNue3D->GetNbinsY();
    int nbinsz= hNue3D->GetNbinsZ();

    for (int ee = 1; ee < nbinsz+1; ee++){
      for (int pz = 1; pz < nbinsx+1; pz++){
        for (int pt = 1; pt < nbinsy+1; pt++){
          double ww = fNumuPiWeights->GetBinContent(pz,pt);
          double nnu = hNue3D->GetBinContent(pz,pt,ee);
          if(nnu>0 && ww>0)       hNue3D->SetBinContent(pz,pt,ee, ww*nnu);        
        }
      }
    }

    std::unique_ptr<TH1D> res(hNue3D->ProjectionZ(UniqueName().c_str()));

    fNueCCFromPiEstim = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(res->GetArray(), res->GetNbinsX()+2)),
                                 HistAxis(fData.GetLabels(), fData.GetBinnings()),
                                 thePOT, 0);
    delete hNue3D;
  }
  
  void BENDecomp::NueEstimateFromKa() const {
    fNueCCFromKaEstim = fNueCCFromKa;
    fNueCCFromKaEstim.Scale(fKaonNormalization);
  }

  void BENDecomp::MakeWeightsNumuFromKaon() const{
    if(!fIsFixedKaScale){
      double pot =  fNumuSelData.POT();
      double pionScale = 1.; double  kaonScale = 1.;
      double epsilon = 1;
      //hack hack hack hack
      //pion contained .75 to 3 GeV
      //kaon uncontained 4.5 to 10

      std::unique_ptr<TH1D> hpionContained(fNumuSelCCFromPi.ToTH1(pot));
      int bin_lo_pi = hpionContained->FindBin(0.75);
      int bin_hi_pi = hpionContained->FindBin(3);
      std::unique_ptr<TH1D> hkaonUncontain(fNumuUncontainCCFromKa.ToTH1(pot));
      int bin_lo_ka = hkaonUncontain->FindBin(4.5);
      int bin_hi_ka = hkaonUncontain->FindBin(10);

      while(epsilon>1E-6) {
        
        Spectrum kaonCorrect = fNumuSelCCFromKa;
        kaonCorrect.Scale (kaonScale);
        Spectrum rawContained = (fNumuSelData - kaonCorrect - fNumuSelBkg -  
                                 fNumuSelCCFromOther);
        
        std::unique_ptr<TH1D> hrawContained(rawContained.ToTH1(pot));
        
        pionScale = (hrawContained->Integral(bin_lo_pi,bin_hi_pi))/
          (hpionContained->Integral(bin_lo_pi,bin_hi_pi));
        
        Spectrum pionCorrect = fNumuUncontainCCFromPi;
        pionCorrect.Scale(pionScale);
        
        Spectrum rawUncontain = (fNumuUncontainData - pionCorrect - fNumuUncontainBkg -  
                                 fNumuUncontainCCFromOther);
        
    
        std::unique_ptr<TH1D> hrawUncontain(rawUncontain.ToTH1(pot));
        epsilon = kaonScale;
        kaonScale = (hrawUncontain->Integral(bin_lo_ka,bin_hi_ka))/(hkaonUncontain->Integral(bin_lo_ka,bin_hi_ka));
        epsilon = abs(epsilon - kaonScale)/kaonScale;
      }
      fKaonNormalization = kaonScale;
    }

}

  void BENDecomp::MakeWeightsNumuFromPion() const{
    double thePOT = fNumuSelData.POT();
    std::unique_ptr<TH1D> hratioNumu(((fNumuSelData - fNumuSelBkg - fNumuSelCCFromOther)/(fNumuSelCCFromPi + fNumuSelCCFromKa)).ToTH1());
    auto hNumu3D = fNumuSelCCFromPiPtPz.ToTH3(thePOT);
    TH2* tmp = (TH2*)hNumu3D->Project3D("yx");
    auto hWeight2D = (TH2*) tmp->Clone(); 
    delete tmp;
    tmp = NULL;

    int nbinsx= hNumu3D->GetNbinsX();
    int nbinsy= hNumu3D->GetNbinsY();
    int nbinsz= hNumu3D->GetNbinsZ();

    for (int pz = 1; pz < nbinsx+1; pz++){
      for (int pt = 1; pt < nbinsy+1; pt++){
        double sum=0,wsum=0/*,tmpn=0*/;
        for (int ee = 1; ee < nbinsz+1; ee++){
          double nnu = hNumu3D->GetBinContent(pz,pt,ee);
          double wei = hratioNumu->GetBinContent(ee);
          double ene = hratioNumu->GetBinCenter (ee);
          wsum+= nnu * wei*ene;
          sum+=  nnu * ene;
//        tmpn+= nnu;
        }
        if(wsum > 0)    hWeight2D->SetBinContent(pz,pt,wsum/sum);
/*
        if(wsum <= 0 && tmpn > 0) {
                        std::cout<< "WARNING: Check out MakeWeightsNumuFromPion() nnu = "
                         << tmpn <<" while weight not updated here: " <<hWeight2D->GetBinContent(pz, pt);       // seems to be ok so far...
        }
*/
      }
     }

    Simulation();
    CornerLabel("#nu-beam");
    //hWeight2D->SetTitle("Pi weights;pz (GeV/c);pT (GeV/c);weight");
    fNumuPiWeights = (TH2*) hWeight2D->Clone();
    delete hWeight2D;
    delete hNumu3D;
  }
  //----------------------------------------------------------------------
  Spectrum BENDecomp::NumuComponent() const
  {
    if(!fIsDecomposed) 
      Decompose();
    return (fData - NueEstimate())*(fNumu/fNotNue); 
  }

  //----------------------------------------------------------------------
  Spectrum BENDecomp::AntiNumuComponent() const
  {
    if(!fIsDecomposed) 
      Decompose();
    return (fData - NueEstimate())*(fAntiNumu/fNotNue); 
 }

  //-------------------------    WARNING    ------------------------------
  //----------------------------------------------------------------------
  Spectrum BENDecomp::NCTotalComponent() const
  {
    if(!fIsDecomposed) 
      Decompose();
    return (fData - NueEstimate())*(fNCTot/fNotNue);
  }

  //----------------------------------------------------------------------
  Spectrum BENDecomp::NCComponent() const 
  {
    std::cout << "BENDecomp::NCComponent() not implemented" << std::endl; 
    abort();
  }

  //----------------------------------------------------------------------
  Spectrum BENDecomp::NCAntiComponent() const 
  {
    std::cout << "BENDecomp::NCAntiComponent() not implemented" << std::endl; 
    abort();
  }

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

  //----------------------------------------------------------------------
  Spectrum BENDecomp::NueComponent() const
  {
    if(!fIsDecomposed) 
      Decompose();
    return NueEstimate();
  }
  //----------------------------------------------------------------------
  Spectrum BENDecomp::AntiNueComponent() const
  {
    if(!fIsDecomposed) 
      Decompose();
    return  (fData - NueEstimate())*(fAntiNue/fNotNue);
  }

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

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

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

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

    if(fIsFixedKaScale){
      TVectorD v(1);v[0]=fKaonNormalization;
      v.Write("KaonNormalization");
    }
    else{
      fNumuUncontainData.SaveTo(dir, "numu_uncontain_data");
      fNumuUncontainBkg.SaveTo(dir, "numu_uncontain_bkg");
      fNumuUncontainCCFromPi.SaveTo(dir, "numuCC_uncontain_from_pi");
      fNumuUncontainCCFromKa.SaveTo(dir, "numuCC_uncontain_from_ka");
      fNumuUncontainCCFromOther.SaveTo(dir, "numuCC_uncontain_from_other");
    }
    fData.SaveTo(dir, "data_comp");
    fNCTot.SaveTo(dir, "nc_tot_comp");
    fNue.SaveTo(dir, "nue_comp");
    fAntiNue.SaveTo(dir, "antinue_comp");
    fNumu.SaveTo(dir, "numu_comp");
    fAntiNumu.SaveTo(dir, "antinumu_comp");
    fTotal.SaveTo(dir, "total_comp");
    fNotNue.SaveTo(dir, "not_nue_comp");

    fNueCCFromPi.SaveTo(dir, "nueCC_from_pi");
    fNueCCFromKa.SaveTo(dir, "nueCC_from_ka");
    fNueCCFromOther.SaveTo(dir, "nueCC_from_other");
    fNueCCFromPiPtPz.SaveTo(dir, "nueCC_from_pi_pt_pz");

    fNumuSelData.SaveTo(dir, "numu_sel_data");
    fNumuSelBkg.SaveTo(dir, "numu_sel_bkg");
    fNumuSelCCFromPi.SaveTo(dir, "numuCC_from_pi");
    fNumuSelCCFromKa.SaveTo(dir, "numuCC_from_ka");
    fNumuSelCCFromOther.SaveTo(dir, "numuCC_from_other");
    fNumuSelCCFromPiPtPz.SaveTo(dir, "numuCC_from_pi_pt_pz");

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

    std::unique_ptr<BENDecomp> ret(new BENDecomp);
    
    ret->fData     = *Spectrum::LoadFrom(dir, "data_comp");
    ret->fNCTot    = *Spectrum::LoadFrom(dir, "nc_tot_comp");
    ret->fNue      = *Spectrum::LoadFrom(dir, "nue_comp");
    ret->fAntiNue  = *Spectrum::LoadFrom(dir, "antinue_comp");
    ret->fNumu     = *Spectrum::LoadFrom(dir, "numu_comp");
    ret->fAntiNumu = *Spectrum::LoadFrom(dir, "antinumu_comp");
    ret->fTotal    = *Spectrum::LoadFrom(dir, "total_comp");
    ret->fNotNue   = *Spectrum::LoadFrom(dir, "not_nue_comp");
    
    ret->fNueCCFromPi = *Spectrum::LoadFrom(dir, "nueCC_from_pi");
    ret->fNueCCFromKa = *Spectrum::LoadFrom(dir, "nueCC_from_ka");
    ret->fNueCCFromOther = *Spectrum::LoadFrom(dir, "nueCC_from_other");
    ret->fNueCCFromPiPtPz = *Spectrum::LoadFrom(dir, "nueCC_from_pi_pt_pz");

    ret->fNumuSelData = *Spectrum::LoadFrom(dir, "numu_sel_data");
    ret->fNumuSelBkg = *Spectrum::LoadFrom(dir, "numu_sel_bkg");
    ret->fNumuSelCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_from_pi");
    ret->fNumuSelCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_from_ka");
    ret->fNumuSelCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_from_other");
    ret->fNumuSelCCFromPiPtPz = *Spectrum::LoadFrom(dir, "numuCC_from_pi_pt_pz");

    if(dir->GetListOfKeys()->Contains("KaonNormalization")){
      auto v = *(TVectorD*)dir->Get("KaonNormalization");
      ret->fKaonNormalization = v[0];
      ret->fIsFixedKaScale=true;
    }
    else{
      ret->fNumuUncontainData = *Spectrum::LoadFrom(dir, "numu_uncontain_data");
      ret->fNumuUncontainBkg = *Spectrum::LoadFrom(dir, "numu_uncontain_bkg");
      ret->fNumuUncontainCCFromPi = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_pi");
      ret->fNumuUncontainCCFromKa = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_ka");
      ret->fNumuUncontainCCFromOther = *Spectrum::LoadFrom(dir, "numuCC_uncontain_from_other");
    }

    ret->fNumuPiWeights = 0;

    delete dir;

    return ret;
  }
  void BENDecomp::SavePlotsKa(TDirectory * dir){
    if(fIsFixedKaScale) {
      std::cerr << "Using a fixed kaon scale " <<  fKaonNormalization
                << ". No plots will be saved " << std::endl;
      return;
    }
    if(!fIsDecomposed) Decompose(); 
    rootlogon();
    auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500);
    dir->cd();
    auto pot = fNumuUncontainData.POT();
    auto hdata = fNumuUncontainData.ToTH1(pot);
    hdata->SetMarkerStyle(kFullCircle);
    hdata->Scale(.001);
    std::vector<std::pair<Spectrum, int> > sp_col = {{fNumuUncontainBkg,kGray+1},
                                                     {fNumuUncontainCCFromOther,kGreen+2},
                                                     {fNumuUncontainCCFromKa,kMagenta -9},
                                                     {fNumuUncontainCCFromPi,kAzure +4}};
    auto hs = new THStack();
    for (auto sp:sp_col){
      auto hist  = sp.first.ToTH1(pot);
      hist->Scale(.001);
      hist->SetFillColorAlpha(sp.second,0.3);
      hs->Add(hist, "hist");
    }

    hdata->Draw("e1");

    hdata->GetYaxis()->SetTitle("10^{3} Events/ 11.0 #times 10^{20} POT");
    hdata->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
    hs->Draw("same");
    hdata->Draw("e1 same");
    
    TLegend *leg = new TLegend(0.6, 0.55, 0.89,0.89);
    leg->SetHeader("Uncontained");
    TGraph * dummy = new TGraph;
    dummy->SetLineWidth(3);
    dummy->SetMarkerStyle(kFullCircle); leg->AddEntry(dummy->Clone(), "#nu_{#mu} data","lpe");
    dummy->SetMarkerStyle(0);
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kAzure+4); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from #pi^{#pm}","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta-9); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from K^{#pm}","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kGreen +2 ); leg->AddEntry(dummy->Clone(),"#nu_{#mu} from Other","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1 ); leg->AddEntry(dummy->Clone(),"Background","f");
    leg->Draw();
    CenterTitles(hdata);
    Preliminary();
    CornerLabel("#nu-beam");
    c1->Write("numuUncontained");

    auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
    hratioNueKa->GetYaxis()->SetTitle("#nu_{e} Ka Estimate / Ka MC");
    hratioNueKa->GetYaxis()->SetRangeUser(0, 0.2 + hratioNueKa->GetMaximum());
    hratioNueKa->Draw("hist");
    hratioNueKa->Write("nue1Dka");

  }
  
  void BENDecomp::SavePlotsPi(TDirectory * dir) {
    if(!fIsDecomposed) Decompose(); 
    rootlogon();
    auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,600);
    dir->cd();
    auto pot = fNumuSelData.POT();
    auto hdata2 = fNumuSelData.ToTH1(pot);
    hdata2->SetMarkerStyle(kFullCircle);
    hdata2->Scale(.001);
    std::vector<std::pair<Spectrum, int> > sp_col2 = {{fNumuSelBkg,kGray +1},
             {fNumuSelCCFromOther,kGreen +2},
             {fNumuSelCCFromKa,kMagenta -9 },
             {fNumuSelCCFromPi,kAzure +4}};
    auto hs2 = new THStack();
    for (auto sp:sp_col2){
      auto hist2  = sp.first.ToTH1(pot);
      hist2->SetFillColorAlpha(sp.second,0.3);
      hist2->Scale(.001);
      hs2->Add(hist2, "hist");
    }

    hdata2->Draw("e1");
    hdata2->GetYaxis()->SetTitle("10^{3} Events/ 11.0 #times 10^{20} POT");
    hdata2->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
    hs2->Draw("same");
    hdata2->Draw("e1 same");
    gPad->Update();


    TLegend *leg2 = new TLegend(0.6, 0.55, 0.76,0.89);
    leg2->SetHeader("Contained");
    TGraph * dummy = new TGraph;
    dummy->SetLineWidth(3);
    dummy->SetMarkerStyle(kFullCircle); leg2->AddEntry(dummy->Clone(), "#nu_{#mu} data","lpe");
    dummy->SetMarkerStyle(0);
    dummy->SetLineColor(kBlack);  dummy->SetFillColor(kAzure +4 ); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from #pi^{#pm}","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kMagenta -9 ); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from K^{#pm}","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kGreen +2); leg2->AddEntry(dummy->Clone(),"#nu_{#mu} from Other","f");
    dummy->SetLineColor(kBlack); dummy->SetFillColor(kGray +1); leg2->AddEntry(dummy->Clone(),"Background","f");

    leg2->Draw();
    CenterTitles(hdata2);
    Preliminary();
    CornerLabel("#nu-beam");
    c1->Write("numuContained");

    // Numu weights1D
    auto hratioNumu = ((fNumuSelData - fNumuSelBkg - fNumuSelCCFromOther)/(fNumuSelCCFromPi + fNumuSelCCFromKa)).ToTH1(); 
    hratioNumu->GetYaxis()->SetTitle("#nu_{#mu} (Data-Other) / (MC-Other)");

    TLine *line = new TLine(0.0,1.0,27.0,1.0);
    line->SetLineWidth(2);
    line->SetLineColor(kRed);
    hratioNumu->Draw("hist");
    CenterTitles(hratioNumu);
    Simulation();
    CornerLabel("#nu-beam");
    line->Draw("same");
    c1->Write("numu1Dwgt");
    // Numu pt pz wgt 2D
    fNumuPiWeights->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
    fNumuPiWeights->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
    fNumuPiWeights->SetTitle("Pion Weights");
    fNumuPiWeights->Draw("colz");
    CenterTitles(fNumuPiWeights);
    Simulation();
    CornerLabel("#nu-beam");
    c1->Write("numu2Dwgt");
    // Numu pt pz count 2D
    auto hnumu3d = fNumuSelCCFromPiPtPz.ToTH3(fNumuSelCCFromPiPtPz.POT());
    hnumu3d->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
    hnumu3d->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
    hnumu3d->GetZaxis()->SetTitle("Events");
    hnumu3d->Project3D("yx")->Draw("colz");
    Simulation();
    CornerLabel("#nu-beam");
    gPad->Update();
    c1->Write("numu2D");
    delete hnumu3d;
    // Nue pt pz count 2D
    auto hnue3d = fNueCCFromPiPtPz.ToTH3(fNueCCFromPiPtPz.POT());
    hnue3d->GetXaxis()->SetTitle("Ancestor forward momentum p_{z} (GeV/c)");
    hnue3d->GetYaxis()->SetTitle("Ancestor trans. momentum p_{t} (GeV/c)");
    hnue3d->GetZaxis()->SetTitle("Events");
    //hnue3d->SetTitle(";Ancestor forward momentum p_{z} (GeV/c);Ancestor trans. momentum p_{t} (GeV/c); Events");
    // hnue3d->GetZaxis()->Scale(1/1000);
    hnue3d->Project3D("yx")->Draw("colz");
    hnue3d->SetTitle("");
    Simulation();
    CornerLabel("#nu-beam");
    gPad->Update();
    c1->Write("nue2D");
    delete hnue3d;
    // Nue Pi estimate/mc ratio
    auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
    hratioNuePi->GetYaxis()->SetTitle("#nu_{e} Pi Estimate / Pi MC");
    hratioNuePi->GetYaxis()->SetRangeUser(0.94, 0.2+hratioNuePi->GetMaximum());
    hratioNuePi->SetTitle("");
    hratioNuePi->Draw("hist");
    line->Draw("same");
    Simulation();
    CornerLabel("#nu-beam");
    CenterTitles(hratioNuePi);
    c1->Write("nue1Dpi");

  }
  
  void  BENDecomp::SavePlots(TDirectory * dir) {
    if(!fIsDecomposed) Decompose();
    rootlogon();
    dir->cd();
    auto c1 = new TCanvas(UniqueName().c_str(),"c1",800,500); 
    double thePOT=fData.POT();
    
    auto hs = new THStack();
    auto hs2 = new THStack(); 

    std::vector<std::pair<Spectrum, int> > sp_col = {
      {fNueCCFromOther,kGray+1},
      {fNueCCFromKaEstim, kMagenta+2},
      {fNueCCFromPiEstim, kAzure +4}
    };
    
    std::vector<std::pair<Spectrum, int> > sp_col2 = {
      {fNueCCFromOther,kGray+1},
      {fNueCCFromKa, kMagenta+2},
      {fNueCCFromPi, kAzure+4}
    };

    for (auto sp:sp_col){
      auto hist  = sp.first.ToTH1(thePOT, sp.second);
      hs->Add(hist, "hist");
    }

    for (auto sp:sp_col2){
      auto hist  = sp.first.ToTH1(thePOT, sp.second);
      hist->SetLineStyle(2);
      hs2->Add(hist, "hist");
    }
//    fData.ToTH1(thePOT)->Draw("ep");
    NueEstimate()    .ToTH1(thePOT,kPink+9,1)->Draw("hist");
    hs->Draw("same");
    hs2->Draw("same");
    fNue             .ToTH1(thePOT,kPink+9,2)->Draw("hist same");
    NueEstimate()    .ToTH1(thePOT,kPink+9,1)->Draw("hist same");

    gStyle->SetTextSize(0.03);
    TLegend * leg = new TLegend(0.7, 0.6, 0.89,0.89);
    TGraph * dummy = new TGraph;
    dummy->SetLineWidth(3);
    dummy->SetLineColor(kPink+9); leg->AddEntry(dummy->Clone(),"Total Beam #nu_{e}","l");
    dummy->SetLineColor(kAzure+4); leg->AddEntry(dummy->Clone(),"#nu_{e} from Pi","l");
    dummy->SetLineColor(kMagenta+2); leg->AddEntry(dummy->Clone(),"#nu_{e} from Ka","l");
    dummy->SetLineColor(kGray+1); leg->AddEntry(dummy->Clone(),"#nu_{e} from Other","l");
    dummy->SetLineColor(kBlack  ); dummy->SetLineStyle(2); 
    leg->AddEntry(dummy->Clone(),"Uncorrected MC","l");
    leg->Draw();
    Simulation();
    CornerLabel("#nu-beam");
    c1->Write("BENresult");

    // Nue estimate/mc ratio
    auto hratioNue = (NueEstimate()/fNue).ToTH1();
    auto hratioNuePi = (fNueCCFromPiEstim/fNueCCFromPi).ToTH1();
    auto hratioNueKa = (fNueCCFromKaEstim/fNueCCFromKa).ToTH1();
    hratioNue->GetYaxis()->SetTitle("#nu_{e} Estimate / MC");
    hratioNue->GetYaxis()->SetRangeUser(0, hratioNue->GetMaximum()+0.2);
    hratioNue->SetLineColor(kPink+9);
    hratioNuePi->SetLineColor(kAzure+4);
    hratioNueKa->SetLineColor(kMagenta+2);
    hratioNue->Draw("hist");
    hratioNue->SetTitle("");
    hratioNuePi->Draw("hist same");
    hratioNueKa->Draw("hist same");
//    leg->GetListOfPrimitives()->RemoveLast();
//    leg->GetListOfPrimitives()->RemoveLast();
//    leg->SetY1(0.15);
//    leg->SetY2(0.45);
    leg->Draw();
    CenterTitles(hratioNue);
    Simulation();
    CornerLabel("#nu-beam");
    c1->Write("BENratios");
  }

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

  Spectrum BENDecomp::MC_NueComponent() const {
    return fNue;
  }
  Spectrum BENDecomp::MC_AntiNueComponent() const {
    return fAntiNue;
  }
  Spectrum BENDecomp::MC_NumuComponent() const {
    return fNumu;
  }
  Spectrum BENDecomp::MC_AntiNumuComponent() const {
    return fAntiNumu;
  }
  Spectrum BENDecomp::MC_NCTotalComponent() const {
    return fNCTot;
  }
  Spectrum BENDecomp::Data_Component() const {
    return fData;
  }

  void BENDecomp::ReduceForBEN2020Decaf(caf::StandardRecord* sr,
                                        const caf::SRProxy* srProxy)
  {
    
    ClearMultiNuInfo(sr);
    ClearSecondaryTrackInfo(sr);

    if(sr->hdr.det == caf::kFARDET) return;
   
    if(!ana::kNue2020NDDecafCut(srProxy) && (ana::kNumuUncontainNDDecafCut(srProxy) || ana::kNumuContainNDDecafCut(srProxy))){

          // remove these completely
          sr->parent = caf::SRParentBranch();
          sr->training = caf::SRTrainingBranch();
          sr->me = caf::SRMichelE();
          
          // slc branch
          // keeping some basic slice info for various cuts
          unsigned int tempncontplanes = sr->slc.ncontplanes;
          unsigned int tempnhit = sr->slc.nhit;
          unsigned int tempfirstplane = sr->slc.firstplane;
          unsigned int templastplane = sr->slc.lastplane;
          unsigned int tempncellsfromedge = sr->slc.ncellsfromedge;
          float tempcalE = sr->slc.calE;
          caf::SRVector3D tempmeanpos = sr->slc.meanpos;
          caf::SRVector3D tempboxmax = sr->slc.boxmax;
          caf::SRVector3D tempboxmin = sr->slc.boxmin;
          
          sr->slc = caf::SRSlice();
          
          sr->slc.ncontplanes = tempncontplanes;
          sr->slc.nhit = tempnhit;
          sr->slc.firstplane = tempfirstplane;
          sr->slc.lastplane = templastplane;
          sr->slc.ncellsfromedge = tempncellsfromedge;
          sr->slc.calE = tempcalE;
          sr->slc.meanpos = tempmeanpos;
          sr->slc.boxmin = tempboxmin;
          sr->slc.boxmax = tempboxmax;
          
          // energy branch
          // keeping only a few numu energy vars
          float temprecotrkcchadE = sr->energy.numu.recotrkcchadE;
          float temptrkccE = sr->energy.numu.trkccE;
          float temptrknonqeE = sr->energy.numu.trknonqeE;
          float temptrkqeE = sr->energy.numu.trkqeE;
          float tempndhadcaltranE = sr->energy.numu.ndhadcaltranE;
          float tempndhadcalcatE = sr->energy.numu.ndhadcalcatE;
          float tempE = sr->energy.numu.E;
          float temphadcalE = sr->energy.numu.hadcalE;
          float temphadtrkE = sr->energy.numu.hadtrkE;
          float tempndtrkcaltranE = sr->energy.numu.ndtrkcaltranE;
          float tempndtrkcalcatE = sr->energy.numu.ndtrkcalcatE;
          float tempndtrkcalactE = sr->energy.numu.ndtrkcalactE;
          float tempndtrklencat = sr->energy.numu.ndtrklencat;
          float tempndtrklenact = sr->energy.numu.ndtrklenact;

          sr->energy.numu = caf::SRNumuEnergy();

          sr->energy.numu.recotrkcchadE = temprecotrkcchadE;
          sr->energy.numu.trkccE = temptrkccE;
          sr->energy.numu.trknonqeE = temptrknonqeE;
          sr->energy.numu.trkqeE = temptrkqeE;
          sr->energy.numu.ndhadcaltranE = tempndhadcaltranE;
          sr->energy.numu.ndhadcalcatE = tempndhadcalcatE;
          sr->energy.numu.E = tempE;
          sr->energy.numu.hadcalE = temphadcalE;
          sr->energy.numu.hadtrkE = temphadtrkE;
          sr->energy.numu.ndtrkcaltranE = tempndtrkcaltranE;
          sr->energy.numu.ndtrkcalcatE = tempndtrkcalcatE;
          sr->energy.numu.ndtrkcalactE = tempndtrkcalactE;
          sr->energy.numu.ndtrklencat = tempndtrklencat;
          sr->energy.numu.ndtrklenact = tempndtrklenact;
          
          ResetBPFEnergy(sr);
         
          // sel branch 
          
          sr->sel.bpfid = caf::SRBpfId();
          sr->sel.xnuepid = caf::SRXnue();
          
          ResetLEMInfo(sr);
          ResetLIDInfo(sr);
          ResetRVPInfo(sr);
          ResetCosRejInfo(sr);
          
          //
          // track branch
          // need number of cosmic tracks, therefore we create a dummy vector 
          sr->trk.cosmic.tracks.clear(); 
          sr->trk.cosmic.tracks.resize(sr->trk.cosmic.ntracks);
          
          // need start and stop coordinates
          if(sr->trk.kalman.ntracks > 0){
            std::vector<caf::SRVector3D> temptrackstart;
            std::vector<caf::SRVector3D> temptrackstop;
            for(int i = 0; i < (int)sr->trk.kalman.tracks.size(); i++){
              temptrackstart.push_back(sr->trk.kalman.tracks[i].start);
              temptrackstop.push_back(sr->trk.kalman.tracks[i].stop);
            }
            float temptrackslen = sr->trk.kalman.tracks[0].len;

            sr->trk.kalman.tracks.clear();
            sr->trk.kalman.tracks.resize(sr->trk.kalman.ntracks);

            for(int i = 0; i < (int)sr->trk.kalman.tracks.size(); i++){
              sr->trk.kalman.tracks[i].start = temptrackstart[i];
              sr->trk.kalman.tracks[i].stop = temptrackstop[i];
            }
            sr->trk.kalman.tracks[0].len = temptrackslen;
          }
          
          // vertex branch
          // need vtx coordinates
          if(sr->vtx.elastic.IsValid && sr->vtx.elastic.fuzzyk.npng > 0){
            float tempvtxX = sr->vtx.elastic.vtx.X();
            float tempvtxY = sr->vtx.elastic.vtx.Y();
            float tempvtxZ = sr->vtx.elastic.vtx.Z();
            float tempnshwlid = sr->vtx.elastic.fuzzyk.nshwlid;
            float tempnpng = sr->vtx.elastic.fuzzyk.npng;
            float tempdedx0 = sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.dedx0;
            int templongestidx = sr->vtx.elastic.fuzzyk.longestidx;
            float templen = sr->vtx.elastic.fuzzyk.png[templongestidx].len; 
            
            std::vector<caf::SRVector3D> tempshwlidstart;
            std::vector<caf::SRVector3D> tempshwlidstop;
            std::vector<caf::SRVector3D> temppngstart;
            std::vector<caf::SRCVNParticleResult> tempcvnpart;
            std::vector<float> tempshwlidcale;

            for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.nshwlid; pngidx++){
             tempshwlidstart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.start);
             tempshwlidstop.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.stop);
            }
            for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.npng; pngidx++){
             temppngstart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].start);
             tempcvnpart.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].cvnpart);
             tempshwlidcale.push_back(sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.calE);
            }

                        
            sr->vtx.elastic.IsValid = true;
            sr->vtx.elastic.vtx.x = tempvtxX;
            sr->vtx.elastic.vtx.y = tempvtxY;
            sr->vtx.elastic.vtx.z = tempvtxZ;
           
            sr->vtx.elastic.fuzzyk.npng = tempnpng;
            sr->vtx.elastic.fuzzyk.nshwlid = tempnshwlid;
            sr->vtx.elastic.fuzzyk.png.resize(tempnpng);
            for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.nshwlid; pngidx++){
              sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.start = tempshwlidstart[pngidx];
              sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.stop = tempshwlidstop[pngidx];
            }
            for(int pngidx = 0; pngidx < (int)sr->vtx.elastic.fuzzyk.npng; pngidx++){
              sr->vtx.elastic.fuzzyk.png[pngidx].start = temppngstart[pngidx];
              sr->vtx.elastic.fuzzyk.png[pngidx].cvnpart = tempcvnpart[pngidx];
              sr->vtx.elastic.fuzzyk.png[pngidx].shwlid.calE = tempshwlidcale[pngidx];
            }
            sr->vtx.elastic.fuzzyk.png[0].shwlid.lid.dedx0 = tempdedx0;
            sr->vtx.elastic.fuzzyk.longestidx = templongestidx;
            sr->vtx.elastic.fuzzyk.png[templongestidx].len = templen;
            
            ClearHoughVertexInfo(sr);
          }
         
          // truth branch
          if(sr->hdr.ismc == 1 && sr->mc.nnu > 0) { 
                  
                  sr->mc.nu[0].rwgt.ppfx.vuniv.clear();
                  sr->mc.nu[0].rwgt.ppfx.nvuniv = 0;              
                  
                  ClearMichelTruthInfo(sr);
          }          

    }

  }

}