ExtrapSterile.h

Go to the documentation of this file.

#pragma once

#include <memory>

#include "CAFAna/Core/Cut.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Extrap/ModularExtrapComponent.h"
#include "CAFAna/Core/OscillatableSpectrum.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/Cuts/TruthCuts.h"

class TDirectory;

namespace ana
{
  // This file implements a sterile neutrino compatible extrapolation.
  // If multiple sterile extrapolations are written,
  // This file will need to contain an generic sterile extrapolation,
  // class IExtrapSterile: public IExtrap
  // The individual sterile extrapolations will inherit from this

  /// A ModularExtrap compatible with sterile oscillations.
  /// Only extrapolates one sample.
  class ModularExtrapSterile: public ModularExtrap
  {
    /// Sets up all components to use FD MC--internal use only.
    /** Use a named constructor to create a ModularExtrapSterile. */
    ModularExtrapSterile(
      SpectrumLoaderBase& farMCswapLoader,
      SpectrumLoaderBase& farMCnonswapLoader,
      SpectrumLoaderBase& farMCtauswapLoader,
      const HistAxis& axis,
      const Cut& fdcut,
      const SystShifts& shiftMC,
      const Var& weight
    );

  public:
    /// Creates a NC disappearance extrapolation.
    /** NCSurvDecomp and NCNDcut are for survival components selected as NCs.
        NumuOscDecomp and NumuNDCut are for components oscillated from numus,
        which then get selected as NCs. */
    static ModularExtrapSterile NCDisappearance(
      Loaders& loaders,
      const IDecomp& NCSurvDecomp,
      const IDecomp& NumuOscDecomp,
      const HistAxis& axis,
      const HistAxis& axisNumuND,
      const Cut& fdcut,
      const Cut& NCNDcut,
      const Cut& NumuNDcut,
      const SystShifts& shiftMC = kNoShift,
      const Var& weight = kUnweighted
    );

    /// Creates a NC disappearance extrapolation with individual spectrum loaders.
    /** NCSurvDecomp and NCNDcut are for survival components selected as NCs.
        NumuOscDecomp and NumuNDCut are for components oscillated from numus,
        which then get selected as NCs. */
    static ModularExtrapSterile NCDisappearance(
      SpectrumLoaderBase& nearMCLoader,
      SpectrumLoaderBase& farMCswapLoader,
      SpectrumLoaderBase& farMCnonswapLoader,
      SpectrumLoaderBase& farMCtauswapLoader,
      const IDecomp& NCSurvDecomp,
      const IDecomp& NumuOscDecomp,
      const HistAxis& axis,
      const HistAxis& axisNumuND,
      const Cut& fdcut,
      const Cut& NCNDcut,
      const Cut& NumuNDcut,
      const SystShifts& shiftMC = kNoShift,
      const Var& weight = kUnweighted
    );

    /// Creates a NC disappearance "extrapolation" by returning FD MC.
    static ModularExtrapSterile TrivialExtrapNC(
      Loaders& loaders,
      const HistAxis& axis,
      const Cut& fdcut,
      const SystShifts& shiftMC = kNoShift,
      const Var& weight = kUnweighted
    );

    /// Creates a NC disappearance "extrapolation" by returning FD MC,
    /// using individual SpectrumLoaders
    static ModularExtrapSterile TrivialExtrapNC(
      SpectrumLoaderBase& farMCswapLoader,
      SpectrumLoaderBase& farMCnonswapLoader,
      SpectrumLoaderBase& farMCtauswapLoader,
      const HistAxis& axis,
      const Cut& fdcut,
      const SystShifts& shiftMC = kNoShift,
      const Var& weight = kUnweighted
    );

    // Prevent copying because we own objects on the free store.
    ModularExtrapSterile(const ModularExtrapSterile&) = delete;
    ModularExtrapSterile& operator=(const ModularExtrapSterile&) = delete;
    ModularExtrapSterile(ModularExtrapSterile&&) = default;
    ModularExtrapSterile& operator=(ModularExtrapSterile&&) = default;
    virtual ~ModularExtrapSterile() = default;

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

    /// Return the oscillatable NC component
    OscillatableSpectrum OscNCComponent() const;

    /// Override the ModularExtrap method
    Spectrum NCTotalComponent() override;

    /// Return the proportion of NCs that originate as a given neutrino flavor/sign
    Ratio NCNueProportion() const;
    Ratio NCAntiNueProportion() const;
    Ratio NCNumuProportion() const;
    Ratio NCAntiNumuProportion() const;

  protected:
    /// Constructor necessary for loading an extrapolation from file
    ModularExtrapSterile(ModularExtrap&& load);

    std::unique_ptr<Spectrum> fNCNueNumerator;
    std::unique_ptr<Spectrum> fNCAntiNueNumerator;
    std::unique_ptr<Spectrum> fNCNumuNumerator;
    std::unique_ptr<Spectrum> fNCAntiNumuNumerator;
    std::unique_ptr<Spectrum> fNCDenominator;
  };

  /// Adapted from TruthCuts.h,
  /// these cuts select NC events discerning neutrino flavor.
  class NCFlavSel
  {
  public:
    NCFlavSel(int pdg, int pdgorig) : fPdg(pdg), fPdgOrig(pdgorig)
    {
    }

    bool operator()(const caf::SRProxy* sr) const;
  protected:
    int fPdg, fPdgOrig;
  };

  /// Select NC \f$ \nu_e\to\nu_e \f$
  const Cut kNCBeamNue(NCFlavSel(12, 12));
  /// Select NC \f$ \nu_\mu\to\nu_e \f$
  const Cut kNCNueApp(NCFlavSel(12, 14));
  /// Select NC \f$ \nu_\mu\to\nu_\mu \f$
  const Cut kNCBeamNumu(NCFlavSel(14, 14));
  /// Select NC \f$ \nu_e\to\nu_\mu \f$
  const Cut kNCNumuApp(NCFlavSel(14, 12));
  /// Select NC \f$ \nu_\mu\to\nu_\tau \f$
  const Cut kNCTauFromMu(NCFlavSel(16, 14));
  /// Select NC \f$ \nu_e\to\nu_\tau \f$
  const Cut kNCTauFromE(NCFlavSel(16, 12));
}