ReactorExperiment.cxx

Go to the documentation of this file.

#include "CAFAna/Experiment/ReactorExperiment.h"

#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Vars/FitVarsCPT.h"

#include "Utilities/func/MathUtil.h"

#include "TDirectory.h"
#include "TH1.h"
#include "TObjString.h"

#include <cassert>

namespace ana
{

  //----------------------------------------------------------------------
  double ReactorExperiment::ChiSq(osc::IOscCalcAdjustable* osc,
                                  const SystShifts& /*syst*/) const
  {
    return util::sqr((SSTh13(osc)-fBestFit)/fSigma);
  }

  //----------------------------------------------------------------------
  double ReactorExperiment::SSTh13(osc::IOscCalcAdjustable* osc) const
  {
    return kFitSinSq2Theta13.GetValue(osc);
  }

  //----------------------------------------------------------------------
  double ReactorExperimentNu::SSTh13(osc::IOscCalcAdjustable* osc) const
  {
    return kFitSinSq2Theta13nu.GetValue(osc);
  }

  //----------------------------------------------------------------------
  double ReactorExperimentNubar::SSTh13(osc::IOscCalcAdjustable* osc) const
  {
    return kFitSinSq2Theta13bar.GetValue(osc);
  }

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

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

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

    TH1D params("", "", 2, 0, 2);
    params.SetBinContent(1, fBestFit);
    params.SetBinContent(2, fSigma);
    params.Write("params");

    dir->Write();
    delete dir;

    tmp->cd();
  }

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

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "ReactorExperiment");
    delete tag;

    TH1* params = (TH1*)dir->Get("params");
    assert(params);

    const double bestFit = params->GetBinContent(1);
    const double sigma   = params->GetBinContent(2);

    delete dir;

    return std::unique_ptr<ReactorExperiment>(new ReactorExperiment(bestFit,
                                                                    sigma));
  }

  //----------------------------------------------------------------------
  const ReactorExperiment* DayaBayConstraint2014()
  {
    // Daya Bay, Neutrino 2014:
    // https://indico.fnal.gov/getFile.py/access?contribId=256&sessionId=15&resId=0&materialId=slides&confId=8022
    return new ReactorExperiment(.084, .005);
  }

  //----------------------------------------------------------------------
  const ReactorExperiment* WorldReactorConstraint2015()
  {
    // Weighted average of
    // Daya Bay: arXiv:1505.03456
    // RENO: arXiv:1410.7987
    // Double Chooz: 1406.7763
    return new ReactorExperiment(.086, .005);
  }

  //----------------------------------------------------------------------
  const ReactorExperiment* WorldReactorConstraint2016()
  {
    // PDG website as of Jun 2016
    return new ReactorExperiment(.085, .005);
  }

  //----------------------------------------------------------------------
  const ReactorExperiment* WorldReactorConstraint2017()
  {
    // http://pdg.lbl.gov/2017/tables/rpp2017-sum-leptons.pdf
    // ssth13=0.021+/-0.0011 -> ss2th13=0.082+/-0.004
    return new ReactorExperiment(.082, .004);
  }

  //----------------------------------------------------------------------
  // http://www.nu-fit.org/ https://arxiv.org/abs/1811.05487 NuFit 4.0 (2018)
  // assuming same constraints for neutrinos and antineutrinos
  // without SK-atm, assign conservative errors
  // ssth13=0.022+/-0.00065 -> ss2th13=0.086+/-0.002
  const ReactorExperiment* ReactorConstraintNuFit2018nu()
  {
    return new ReactorExperimentNu(.086, .002);
  }

  //----------------------------------------------------------------------
  const ReactorExperiment* ReactorConstraintNuFit2018bar()
  {
    return new ReactorExperimentNubar(.086, .002);
  }

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

  const ReactorExperiment* WorldReactorConstraint2019()
  {
    // http://pdg.lbl.gov/2019/tables/rpp2019-sum-leptons.pdf
    // ssth13=0.0218+/-0.0007 -> ss2th13=0.085+/-0.003
    return new ReactorExperiment(.085, .003);
  }


}