MichelDecomp.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Decomp/IDecomp.h"

#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Vars/Vars.h"

#include <vector>

class TDirectory;

namespace ana
{
  class Loaders;

  /// This is the Binning scheme used for the # Michel Templates
  extern const Binning MEBinning;
  extern const HistAxis MEAxis;

  struct MDCMPFitResults
  {
    std::vector<double> numuScales;
    std::vector<double> ncScales;
    std::vector<double> nueScales;
  };

  /// Help reduce verbosity when passing hist values between MichelDecomp funcs
  struct MDCMPHelper
  {
    double fSumData;
    double fSumNumu;
    double fSumNC;
    double fSumNue;
    double fNueScale;
    
    // Other scales from input decomp
    double fNumuScale;
    double fNCScale;
  };

  class MichelDecomp: public IDecomp
  {
  public:
    MichelDecomp(SpectrumLoaderBase& NDMCloader,
                 SpectrumLoaderBase& NDdataloader,
                 const HistAxis&     axis,
                 const Cut&          cut,
                 const IDecomp*      nNueEstimate,
                 const SystShifts&   shiftMC = kNoShift,
                 const SystShifts&   shiftData = kNoShift,
                 const Var&          wei = kUnweighted);

    MichelDecomp(SpectrumLoaderBase& NDMCloader,
                 SpectrumLoaderBase& NDdataloader,
                 std::string         label,
                 const Binning&      bins,
                 const Var&          var,
                 const Cut&          cut,
                 const IDecomp*      nNueEstimate,
                 const SystShifts&   shiftMC = kNoShift,
                 const SystShifts&   shiftData = kNoShift,
                 const Var&          wei = kUnweighted);

    MichelDecomp(Loaders&            loaders,
                 const HistAxis&     axis,
                 const Cut&          cut,
                 const IDecomp*      nNueEstimate,
                 const SystShifts&   shiftMC = kNoShift,
                 const SystShifts&   shiftData = kNoShift,
                 const Var&          wei = kUnweighted);

    ~MichelDecomp();

    //nc
    virtual Spectrum NCTotalComponent() const override;
    virtual Spectrum NCComponent()      const override;
    virtual Spectrum NCAntiComponent()  const override;
    //end nc
    virtual Spectrum NumuComponent()     const override;
    virtual Spectrum AntiNumuComponent() const override;
    virtual Spectrum NueComponent()      const override;
    virtual Spectrum AntiNueComponent()  const override;

    virtual Spectrum Data_Component()       const override;

    virtual Spectrum MCTotal()           const{return NCAntiComponent() +NCComponent()+NumuComponent()+AntiNumuComponent()+NueComponent()+AntiNueComponent();}

    // Convenient for making decomposition blessed plots
    Spectrum MC_NCTotalComponent()  const override;
    Spectrum MC_NCComponent()       const override;
    Spectrum MC_NCAntiComponent()   const override;
    Spectrum MC_NumuComponent()     const override;
    Spectrum MC_AntiNumuComponent() const override;
    Spectrum MC_NueComponent()      const override;
    Spectrum MC_AntiNueComponent()  const override;

    virtual Spectrum MC_MCTotal()           const{return MC_NCTotalComponent() + MC_NumuComponent() + MC_AntiNumuComponent() + MC_NueComponent() + MC_AntiNueComponent();}

    void SavePlots(TDirectory* dir) const;


    virtual void SaveTo(TDirectory* dir, const std::string& name) const override;
    static std::unique_ptr<MichelDecomp> LoadFrom(TDirectory* dir, const std::string& name);

    const IDecomp* GetNueEstimateDecomp() const {return fNueEstimate;}
  protected:
    MichelDecomp(Spectrum nc, Spectrum antinc,
                 Spectrum numu, Spectrum antinumu,
                 Spectrum nue, Spectrum antinue, Spectrum data,
                 ReweightableSpectrum michelnc,
                 ReweightableSpectrum michelnumu,
                 ReweightableSpectrum michelnue,
                 ReweightableSpectrum micheldata,
                 const IDecomp* nueestimate)
      : fNC(nc),
        fAntiNC(antinc),
        fNumu(numu),
        fAntiNumu(antinumu),
        fNue(nue),
        fAntiNue(antinue),
        fData(data),
        fMichelNC(michelnc),
        fMichelNumu(michelnumu),
        fMichelNue(michelnue),
        fMichelData(micheldata),
        fNueEstimate(nueestimate),
        fOwnNueEstimate(false),
        isDecomposed(false) {}

    void   Decompose() const;


    // For bins with low stats PropDecomp
    void SplitByProp(MDCMPHelper hSums,
                     std::vector<double> &numuScales,
                     std::vector<double> &ncScales,
                     std::vector<double> &nueScales) const;
    // Or scales from input decomp                   
    void SplitByNue(MDCMPHelper hSums,
                    std::vector<double> &numuScales,
                    std::vector<double> &ncScales,
                    std::vector<double> &nueScales) const;
                    
    void SplitByMichels(int EIdx, MDCMPHelper hSums,
                        std::vector<double> &numuScales,
                        std::vector<double> &ncScales,
                        std::vector<double> &nueScales) const;


    // These functions cover the algebra --->
    double GetNueScale(int EIdx) const;
    double GetNumuScale(int EIdx) const;
    double GetNCScale(int EIdx) const;
    double GetNumuMCContent(int EIdx) const;
    double GetNCMCContent(int EIdx) const;
    double ComputeScaleFactor(int EIdx,
                              MDCMPHelper hSums) const;
    double CalculateMCMean(double h,
                           double Numu,
                           double NC,
                           double Nue,
                           MDCMPHelper hSums) const;
    double MDCMPLogLikelihood(double h,
                              int EIdx,
                              MDCMPHelper hSums) const;
    double ComputeMaxScale(MDCMPHelper hSums) const;
    double GetSum(ReweightableSpectrum spect, int EIdx) const;
    double GetTemplateContent (ReweightableSpectrum spect,
                               int EIdx, int mIdx) const;
    // <---

    // Apply scale factors from fit to MC spectra
    virtual Spectrum MCToDCMPComp(std::vector<double>, Spectrum) const;

    // Called by SavePlots
    void SaveCompPlots(TDirectory* dir, std::string comp,
                       Spectrum dcmp, Spectrum mc) const;
    void SaveTempPlots(TDirectory* dir) const;

    Spectrum fNC;
    Spectrum fAntiNC;
    Spectrum fNumu;
    Spectrum fAntiNumu;
    Spectrum fNue;
    Spectrum fAntiNue;
    Spectrum fData;
    ReweightableSpectrum fMichelNC;
    ReweightableSpectrum fMichelNumu;
    ReweightableSpectrum fMichelNue;
    ReweightableSpectrum fMichelData;

    const IDecomp* fNueEstimate; // Fix beam nue's to different decomp
    bool fOwnNueEstimate;

    mutable bool isDecomposed;
    mutable MDCMPFitResults fFitResult;
  };
}