CutOptimization.h

Go to the documentation of this file.

#pragma once
#include <map>

#include "CAFAna/XSec/SystematicDef.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/XSec/Multiverse.h"

#include "TString.h"
#include "TDirectory.h"
#include "TObjString.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TLegend.h"

namespace ana
{
  class CutOptimization
  {
  public:

    // static functions for plotting histograms that get saved to the
    // debug file
    static void PlotDebug(TH1 * hist, 
                          std::string draw_option, 
                          std::string title, 
                          std::string name);

    static void PlotDebug(std::vector<TH1 *> hists,
                          std::vector<std::string> labels,
                          std::string draw_option, 
                          std::string title, 
                          std::string name,
                          double ymax = -1);

    static void PlotDebug(const std::vector<TH1 *> universes, 
                          Spectrum * nominal,
                          std::string title, 
                          std::string name);

    static void PlotDebug(const std::vector<TH1 *> universes, 
                          TH1 * h_nominal,
                          std::string title, 
                          std::string name);

    enum OptimizeMethod_t {
      kBinByBin,         // bin by bin
      kIntegratedUp,     // integrate bin and above
      kIntegratedDown,   // integrate bin and below
      kIntegratedInsideBounds, // integrate between upper and lower bounds
    };    
    enum FOM_t {
      kdSigmaOverSigma
    };

    CutOptimization(const HistAxis * hist_axis,
                    const Cut * signal_cut,
                    const Cut * selection_cut);

    void DefineNominal(SystematicDef * nominal);
    void DefineSystematic(SystematicDef * syst);
    void DefineMultiverseSystematic(SystematicDef * syst, std::vector<Var> mv_weights);
    void DefineMultiverseSystematic(SystematicDef * syst, std::vector<SystShifts> mv_shifts);

    void SetSpillCuts(const SpillCut &);

    void Go();    
    void SaveTo(TDirectory * dir, const std::string& name) const;

    static std::unique_ptr<CutOptimization> LoadFrom(TDirectory * dir, const std::string& name);
    CutOptimization(const CutOptimization & rhs);
    CutOptimization(CutOptimization && rhs);

    CutOptimization & operator=(const CutOptimization & rhs);
    CutOptimization & operator=(CutOptimization && rhs);
    ~CutOptimization() = default;

    void Optimize(OptimizeMethod_t method, FOM_t fom, double data_pot,
                  TDirectory * save_dir, std::string plot_dump = ".",
                  bool debug = false);


    CutOptimization(Spectrum * nominal_signal,
                    Spectrum * nominal_bkgd,
                    Spectrum * nominal_selected_signal,
                    Spectrum * nominal_selected_bkgd,
                    Spectrum * nominal_selected,
                    std::vector<SystematicDef*> syst_defs,
                    std::map<SystematicDef*, UpDownPair<Spectrum > > syst_signal,
                    std::map<SystematicDef*, UpDownPair<Spectrum > > syst_bkgd,
                    std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_signal,
                    std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_bkgd,
                    std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected,
                    std::vector<SystematicDef*> mv_defs,
                    std::map<SystematicDef*, Multiverse * > mv_signal,
                    std::map<SystematicDef*, Multiverse * > mv_bkgd,
                    std::map<SystematicDef*, Multiverse * > mv_selected_signal,
                    std::map<SystematicDef*, Multiverse * > mv_selected_bkgd,
                    std::map<SystematicDef*, Multiverse * > mv_selected);
    static void Integrate(TH1 *& result, const TH1 * tmp, OptimizeMethod_t method);
    static void Integrate(Multiverse * mv, OptimizeMethod_t method);
    
    std::vector<SystematicDef*> GetSystDefs() { return fSystDefs; }
    std::vector<SystematicDef*> GetMVSystDefs() { return fMVSystDefs; }
    std::vector<SystematicDef*> GetAllSystDefs() { 
      std::vector<SystematicDef*> ret = fSystDefs;
      ret.insert(fSystDefs.end(), fMVSystDefs.begin(), fMVSystDefs.end());
      return ret;
    }

  private:
    TH1 * AbsUncertainty       (const UpDownPair<TH1> syst, const TH1 * nominal) const;
    TH1 * AbsUncertaintySquared(const UpDownPair<TH1> syst, const TH1 * nominal) const;

    static TH1 * ToHist(const Spectrum   * spec, OptimizeMethod_t method, double POT = -1);
    UpDownPair<TH1> ToUpDownHist(Multiverse * mv  , const TH1 * h_nominal);

    void OptimizedSigmaOverSigma(OptimizeMethod_t method, double data_pot,
                                 TDirectory * save_dir, std::string plot_dump = ".",
                                 bool debug = false);

    static TH1 * Efficiency(TH1 * num, TH1 * denom, OptimizeMethod_t method);    
    UpDownPair<TH1> Efficiency(UpDownPair<TH1> num, UpDownPair<TH1> denom, OptimizeMethod_t method);
    Multiverse * Efficiency(Multiverse * num, Multiverse * denom, OptimizeMethod_t method);
    
    SystematicDef * fNominalSystDef;
    Spectrum * fNominalSignal;
    Spectrum * fNominalBkgd;
    Spectrum * fNominalSelectedSignal;
    Spectrum * fNominalSelectedBkgd;
    Spectrum * fNominalSelected;

    std::vector<SystematicDef*> fSystDefs; // todo change to set
    std::map< SystematicDef*, UpDownPair<Spectrum> > fSignal;
    std::map< SystematicDef*, UpDownPair<Spectrum> > fBkgd;
    std::map< SystematicDef*, UpDownPair<Spectrum> > fSelectedSignal;
    std::map< SystematicDef*, UpDownPair<Spectrum> > fSelectedBkgd;
    std::map< SystematicDef*, UpDownPair<Spectrum> > fSelected;
    
    std::vector<SystematicDef*> fMVSystDefs;
    std::map< SystematicDef*, Multiverse *> fMVSignal;
    std::map< SystematicDef*, Multiverse *> fMVBkgd;
    std::map< SystematicDef*, Multiverse *> fMVSelectedSignal;
    std::map< SystematicDef*, Multiverse *> fMVSelectedBkgd;
    std::map< SystematicDef*, Multiverse *> fMVSelected;


    const Cut * fSignalCut = NULL;
    const Cut * fSelectionCut = NULL;
    const HistAxis * fOptiHistAxis = NULL;
  };

  inline UpDownPair<Spectrum> CopyUpDownSpectrum(const UpDownPair<Spectrum> & copy)
  {
    UpDownPair<Spectrum> ret;
    if(copy.up)   ret.up   = new Spectrum(*copy.up);
    if(copy.down) ret.down = new Spectrum(*copy.down);
    return ret;
  }

  inline UpDownPair<Spectrum> CopyUpDownSpectrum(UpDownPair<Spectrum> && copy)
  {
    UpDownPair<Spectrum> ret;
    if(copy.up)   ret.up   = std::move(copy.up);
    if(copy.down) ret.down = std::move(copy.down);
    return ret;
  }

  inline void SaveToUpDownSpectra(const UpDownPair<Spectrum> & save, 
                                  TDirectory * dir) 
  {
    TDirectory * tmp = gDirectory;
    dir->cd();
    TObjString("UpDownPair").Write("type");
    
    if(save.IsTwoSided()) {
      save.up->SaveTo(dir, "up");
      save.down->SaveTo(dir, "down");
    }
    else {
      save.up->SaveTo(dir, "onesided");
    }
    tmp->cd();
  }

  inline UpDownPair<Spectrum> LoadFromUpDownSpectra(TDirectory * dir) 
  {
    TObjString * ptag = (TObjString*) dir->Get("type");
    assert(ptag->GetString() == "UpDownPair");
    
    if(dir->GetDirectory("onesided")) {
      return UpDownPair<Spectrum>{Spectrum::LoadFrom(dir, "onesided").release()};
    }
    else {
      return UpDownPair<Spectrum>
        {Spectrum::LoadFrom(dir, "up").release(),
            Spectrum::LoadFrom(dir, "down").release()};
    }    
  }

}