FitVarsCPT.cxx

Go to the documentation of this file.

#include "CAFAna/Vars/FitVarsCPT.h"

#include "Utilities/func/MathUtil.h"

#include <cassert>

namespace ana
{
  // declared as 'extern' in FitVarsCPT.h
  const FitDmSq32CPT kFitDmSq32nu(osc::ENuSign::kNu);
  const FitDmSq32CPT kFitDmSq32bar(osc::ENuSign::kNuBar);
  const FitDmSq32ScaledCPT kFitDmSq32Scalednu(osc::ENuSign::kNu);
  const FitDmSq32ScaledCPT kFitDmSq32Scaledbar(osc::ENuSign::kNuBar);
  const FitDmSq32CPTHierarchy kFitDmSq32nuNorm( osc::ENuSign::kNu, false );
  const FitDmSq32CPTHierarchy kFitDmSq32barNorm(osc::ENuSign::kNuBar, false );
  const FitDmSq32CPTHierarchy kFitDmSq32nuInv(  osc::ENuSign::kNu, true );
  const FitDmSq32CPTHierarchy kFitDmSq32barInv( osc::ENuSign::kNuBar, true );
  const FitDmSq21CPT kFitDmSq21nu(osc::ENuSign::kNu);
  const FitDmSq21CPT kFitDmSq21bar(osc::ENuSign::kNuBar);
  const FitSinSq2Theta23CPT kFitSinSq2Theta23nu(osc::ENuSign::kNu);
  const FitSinSq2Theta23CPT kFitSinSq2Theta23bar(osc::ENuSign::kNuBar);
  const FitSinSqTheta23CPT kFitSinSqTheta23nu(osc::ENuSign::kNu);
  const FitSinSqTheta23CPT kFitSinSqTheta23bar(osc::ENuSign::kNuBar);
  const FitSinSqTheta23CPTOctant kFitSinSqTheta23nuFirst(osc::ENuSign::kNu,      false);
  const FitSinSqTheta23CPTOctant kFitSinSqTheta23barFirst(osc::ENuSign::kNuBar,  false);
  const FitSinSqTheta23CPTOctant kFitSinSqTheta23nuSecond(osc::ENuSign::kNu,     true);
  const FitSinSqTheta23CPTOctant kFitSinSqTheta23barSecond(osc::ENuSign::kNuBar, true);
  const FitSinSq2Theta13CPT kFitSinSq2Theta13nu(osc::ENuSign::kNu);
  const FitSinSq2Theta13CPT kFitSinSq2Theta13bar(osc::ENuSign::kNuBar);
  const FitSinSq2Theta12CPT kFitSinSq2Theta12nu(osc::ENuSign::kNu);
  const FitSinSq2Theta12CPT kFitSinSq2Theta12bar(osc::ENuSign::kNuBar);


  //---------------------------------------------------------------------------
  double FitDmSq21CPT::GetValue(const osc::IOscCalcAdjustable* osc) const
  { return osc::DowncastToCPT(osc)->GetDmsq21(fNuSign); }

  void FitDmSq21CPT::SetValue(osc::IOscCalcAdjustable* osc,
                              double val) const
  { osc::DowncastToCPT(osc)->SetDmsq21(Clamp(val), fNuSign); }


  //---------------------------------------------------------------------------
  double FitDmSq32CPT::GetValue(const osc::IOscCalcAdjustable* osc) const
    { return osc::DowncastToCPT(osc)->GetDmsq32(fNuSign); }

  void FitDmSq32CPT::SetValue(osc::IOscCalcAdjustable* osc,
                              double val) const
    { osc::DowncastToCPT(osc)->SetDmsq32(Clamp(val), fNuSign); }


  //---------------------------------------------------------------------------
  double FitDmSq32ScaledCPT::GetValue(const osc::IOscCalcAdjustable* osc) const
  { return osc::DowncastToCPT(osc)->GetDmsq32(fNuSign)*1000.0; }

  void FitDmSq32ScaledCPT::SetValue(osc::IOscCalcAdjustable* osc,
                              double val) const
  { osc::DowncastToCPT(osc)->SetDmsq32(Clamp(val)/1000.0, fNuSign); }

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

  double FitSinSq2Theta23CPT::GetValue(const osc::IOscCalcAdjustable* osc) const
    { 
      return util::sqr(sin(2*osc::DowncastToCPT(osc)->GetTh23(fNuSign) )); }

  void FitSinSq2Theta23CPT::SetValue(osc::IOscCalcAdjustable* osc,
                                     double val) const
    { 
      osc::DowncastToCPT(osc)->SetTh23(asin(sqrt( Clamp(val) ))/2, fNuSign); }

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

  double FitSinSqTheta23CPT::GetValue(
      const osc::IOscCalcAdjustable* osc) const
    { return util::sqr(sin( osc::DowncastToCPT(osc)->GetTh23(fNuSign) )); }

  void FitSinSqTheta23CPT::SetValue(osc::IOscCalcAdjustable* osc,
                                    double val) const
    { osc::DowncastToCPT(osc)->SetTh23(asin(sqrt( Clamp(val) )), fNuSign); }

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

  double FitSinSq2Theta13CPT::GetValue(const osc::IOscCalcAdjustable* osc) const
  { 
    return util::sqr(sin(2*osc::DowncastToCPT(osc)->GetTh13(fNuSign) )); }
  void FitSinSq2Theta13CPT::SetValue(osc::IOscCalcAdjustable* osc,
                                     double val) const
    { osc::DowncastToCPT(osc)->SetTh13(asin(sqrt( Clamp(val) ))/2, fNuSign); }

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

  double FitSinSq2Theta12CPT::GetValue(
      const osc::IOscCalcAdjustable* osc) const
    { return util::sqr(sin(2*osc::DowncastToCPT(osc)->GetTh12(fNuSign) )); }

  void FitSinSq2Theta12CPT::SetValue(osc::IOscCalcAdjustable* osc,
                                     double val) const
    { osc::DowncastToCPT(osc)->SetTh12(asin(sqrt( Clamp(val) ))/2, fNuSign); }

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

  SigmaDelta::SigmaDelta(osc::OscCalcCPT* osc,
                         const IConstrainedFitVar &x,
                         const IConstrainedFitVar &y)
    : fXvar(x), fYvar(y),
      fCenterDelta(fXvar.HighLimit() - fYvar.HighLimit()),
      fCenterSigma(fXvar.HighLimit() + fYvar.LowLimit())
  {
    // Check to make sure the allowed area is square. If not, the logic is
    // more complicated, but most pairs of CPT FitVars have matched ranges.
    assert ( fCenterDelta == (fXvar.LowLimit() - fYvar.LowLimit()) );
    assert ( fCenterSigma == (fXvar.LowLimit() + fYvar.HighLimit()) );
    osc->fSigDel[this] = {fXvar.GetValue(osc)+fYvar.GetValue(osc),
                          fXvar.GetValue(osc)-fYvar.GetValue(osc)};
  }

  double SigmaDelta::GetSigma(const osc::IOscCalcAdjustable* osc) const
    { return osc::DowncastToCPT(osc)->fSigDel[this].first; }
  double SigmaDelta::GetDelta(const osc::IOscCalcAdjustable* osc) const
    { return osc::DowncastToCPT(osc)->fSigDel[this].second; }

  void SigmaDelta::SetSigma(osc::IOscCalcAdjustable* osc, double val)
  {
    osc::OscCalcCPT* cpt = osc::DowncastToCPT(osc);
    cpt->fSigDel[this].first = val;
    SetFitVars(cpt);
  }

  void SigmaDelta::SetDelta(osc::IOscCalcAdjustable* osc, double val)
  {
    osc::OscCalcCPT* cpt = osc::DowncastToCPT(osc);
    cpt->fSigDel[this].second = val;
    SetFitVars(cpt);
  }

  double SigmaDelta::DeltaPenalty(double val) const
  {
    return 0; //have to check both vars, so just apply the penalty in sigma
  }

  double SigmaDelta::SigmaPenalty(double val,
                                  osc::IOscCalcAdjustable* osc) const
  {
    osc::OscCalcCPT* cpt = osc::DowncastToCPT(osc);
    double delta = cpt->fSigDel[this].second;

    return (   fXvar.Penalty( (val+delta)/2 , osc)
             + fYvar.Penalty( (val-delta)/2 , osc)
             + 1000*(delta > (fXvar.HighLimit()-fYvar.LowLimit()) )
             + 1000*(delta < (fXvar.LowLimit()-fYvar.HighLimit()) )
           ); //extra 1000 is in case contour hits impossible region
  }

  void SigmaDelta::SetFitVars(osc::OscCalcCPT* osc)
  {
    double sigma = osc->fSigDel[this].first;
    double delta = osc->fSigDel[this].second;
    //if requested value is impossible, try to at least maintain delta
    if( (delta>fCenterDelta)==(sigma>fCenterSigma) )
    {
      fXvar.SetValue(osc, (sigma+delta)/2);
      double xc = fXvar.GetValue(osc);
      fYvar.SetValue(osc, xc-delta);
    }
    else
    {
      fYvar.SetValue(osc, (sigma-delta)/2);
      double yc = fYvar.GetValue(osc);
      fXvar.SetValue(osc, yc+delta);
    }
  }

} //namespace