NCCoordinateConverter.cxx

Go to the documentation of this file.

#include "NCCoordinateConverter.h"

#include "MessageService/MsgService.h"

#include "NCUtils/NCUtility.h"

#include <stdexcept>

CVSID("$Id: NCCoordinateConverter.cxx,v 1.32 2009/11/26 12:03:29 bckhouse Exp $");

using namespace NCType;
using namespace NC::Fitter;
using NC::Utility::SQR;

ClassImp(NC::CoordinateConverter)


NC::SystPars::SystPars()
{
  for(int n = 0; n < kNumSystematicParameters; ++n) fShifts[n] = 0;
}


bool NC::SystPars::operator!=(const SystPars& s) const
{
  for(int n = 0; n < kNumSystematicParameters; ++n)
    if(s.fShifts[n] != fShifts[n]) return true;

  return false;
}


bool NC::SystPars::operator==(const SystPars& s) const
{
  return !(*this != s);
}


double NC::SystPars::ShowerScale(Detector::Detector_t det) const
{
  if(det == Detector::kNear) return 1+fShifts[kAbsoluteHadronicCalibration];

  return (1+fShifts[kRelativeHadronicCalibration])*
         (1+fShifts[kAbsoluteHadronicCalibration]);
}


double NC::SystPars::TrackScale() const
{
  return 1+fShifts[kTrackEnergy];
}


double NC::SystPars::NormScale() const
{
  return 1+fShifts[kNormalization];
}


double NC::SystPars::NDCleaningScale(NCType::EEventType nccc,
                                     double showerEnergy) const
{
  if(nccc == NCType::kCC) return 1;

  if(showerEnergy <= 0.5) return 1 + 0.152*fShifts[kNCNearClean];
  if(showerEnergy <= 1.0) return 1 + 0.029*fShifts[kNCNearClean];
  if(showerEnergy <= 1.5) return 1 + 0.004*fShifts[kNCNearClean];

  return 1;
}


double NC::SystPars::NCBkgScale() const
{
  return 1+fShifts[kNCBackground];
}


double NC::SystPars::CCBkgScale() const
{
  return 1+fShifts[kCCBackground];
}


double NC::SystPars::PenaltyForShifts() const
{
  double pen = 0;

  for(int n = 0; n < kNumSystematicParameters; ++n){
    const double shift = fShifts[n];
    if(shift) pen += shift/NCType::kParams[n].sigma;
  }

  return pen;
}


NC::CoordinateConverter::CoordinateConverter()
{
  for(int n = 0; n < kNumParameters; ++n){
    fLoc[n] = -1;
    fFixedVals[n] = -9999;
  }
}


const Registry& NC::CoordinateConverter::DefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  r.Set("PrecScale",         1.0);

  r.Set("DeltaCP3FlavorVal", 0.0);

  // Whether to fit particular parameters
  r.Set("FitUE3Sqr",            false);
  r.Set("FitTheta13",           false);
  r.Set("FitDelta1",            false);
  r.Set("FitDelta2",            false);

  r.Set("FitTheta14",           true);
  r.Set("FitTheta24",           true);
  r.Set("FitTheta23",           true);

  // Fixed values to use for the parameters if not fitting them
  r.Set("UE3SqrVal",            0.);
  r.Set("Theta13Val",           0.);
  r.Set("Delta1Val",            0.);
  r.Set("Delta2Val",            0.);

  r.Set("Theta14Val",           0.);
  r.Set("Theta24Val",           0.);
  r.Set("Theta23Val",           0.);


  r.LockValues();
  return r;
}


void NC::CoordinateConverter::TryAddSystematic(const Registry& r,
                                               NCType::EFitParam p)
{
  int tmpb;
  TString key = "Fit"+NCType::kParams[p].name;
  if(r.Get(key, tmpb) && tmpb) fLoc[p] = AddSystematicParameter(p);
}


void NC::CoordinateConverter::PrepareSystematics(const Registry& r)
{
  TryAddSystematic(r, kNormalization);
  TryAddSystematic(r, kRelativeHadronicCalibration);
  TryAddSystematic(r, kAbsoluteHadronicCalibration);
  TryAddSystematic(r, kTrackEnergy);
  TryAddSystematic(r, kNCBackground);
  TryAddSystematic(r, kCCBackground);
  TryAddSystematic(r, kNCNearClean);
}


void NC::CoordinateConverter::PrepareOscillations(const Registry& r)
{
  const int    kNumUMu3SqrBins     = 100;
  const double kUMu3SqrStart       = 1e-3;
  const double kUMu3SqrEnd         = 1;
  const double kDeltaUMu3Sqr       = (kUMu3SqrEnd-kUMu3SqrStart)/kNumUMu3SqrBins;

  const int    kNumDeltaMSqrBins   = 300;
  const double kDeltaMSqrStart     = 1.5;
  const double kDeltaMSqrEnd       = 4.5;
  const double kDeltaDeltaMSqr     = (kDeltaMSqrEnd-kDeltaMSqrStart)/kNumDeltaMSqrBins;

  int tmpb;

  switch(fOscillationModel){
  case kThreeFlavor:
    fLoc[kDeltaMSqr32] = AddParameter("dmsq_32", "#Deltam^{2}_{32} (10^{-3} eV^{2})",
                                      kDeltaMSqrStart,
                                      kDeltaMSqrEnd,
                                      kDeltaDeltaMSqr, true);

    if(r.Get("FitUE3Sqr", tmpb) && tmpb){
      fLoc[kTheta13] = AddParameter("theta13", "#theta_{13}",
                                    0, TMath::Pi()*2, -1, true);
    }

    if(r.Get("FitTheta23", tmpb) && tmpb){
      fLoc[kTheta23] = AddParameter("theta23", "#theta_{23}",
                                    0, TMath::Pi()*2, -1, true);
    }

    return;

  case kNoOscillations:
    return;

  case kFourFlavorGeneral:
  case kFourFlavorDelta43IsBig:
  case kFourFlavorDelta43Is0:
  case kFourFlavorDelta41Is0:

    if(r.Get("FitTheta13", tmpb) && tmpb){
      fLoc[kTheta13] = AddParameter("theta13", "#theta_{13}",
                                    0, TMath::Pi()*2, -1, true);
    }

    if(r.Get("FitTheta23", tmpb) && tmpb){
      fLoc[kTheta23] = AddParameter("theta23", "#theta_{23}",
                                    0, TMath::Pi()/2, -1, true);
    }

    if(fOscillationModel != kFourFlavorDelta41Is0){

      if(r.Get("FitTheta14", tmpb) && tmpb){
        fLoc[kTheta14] = AddParameter("theta14", "#theta_{14}",
                                      0, TMath::Pi()*2, -1, true);
      }

      if(r.Get("FitTheta24", tmpb) && tmpb){
        fLoc[kTheta24] = AddParameter("theta24", "#theta_{24}",
                                      0, TMath::Pi()/3, -1, true);
      }

      if(r.Get("FitDelta1", tmpb) && tmpb){
        fLoc[kDelta1] = AddParameter("delta1", "#delta_{1}",
                                     0, TMath::Pi()*2, -1, true);
      }
      
      if(fOscillationModel != kFourFlavorDelta43Is0){
        if(r.Get("FitDelta2", tmpb) && tmpb){
          fLoc[kDelta2] = AddParameter("delta2", "#delta_{2}",
                                       0, TMath::Pi()*2, -1, true);
        }
      }
    } // end if not Delta41=0

    if(fOscillationModel != kFourFlavorDelta43Is0){
      fLoc[kTheta34] = AddParameter("theta34", "#theta_{34}",
                                    0, TMath::Pi()/3, -1, true);
    }


    // TODO - How big should the ranges on all these mass scales be?
    fLoc[kDeltaMSqr31] = AddParameter("dmsq_31", "#Deltam^{2}_{31} (10^{-3} eV^{2})",
                                      kDeltaMSqrStart,
                                      kDeltaMSqrEnd,
                                      kDeltaDeltaMSqr, true);

    if(fOscillationModel != kFourFlavorDelta43IsBig &&
       fOscillationModel != kFourFlavorDelta43Is0 &&
       fOscillationModel != kFourFlavorDelta41Is0){

      fLoc[kDeltaMSqr41] = AddParameter("dmsq_41", "#Deltam^{2}_{41} (10^{-3} eV^{2})",
                                        kDeltaMSqrStart,
                                        kDeltaMSqrEnd,
                                        kDeltaDeltaMSqr, true);

      fLoc[kDeltaMSqr43] = AddParameter("dmsq_43", "#Deltam^{2}_{43} (10^{-3} eV^{2})",
                                        kDeltaMSqrStart,
                                        kDeltaMSqrEnd,
                                        kDeltaDeltaMSqr, true);
    }

    return;

  case kSterileFraction:
    fLoc[kFs] = AddParameter("fs", "f_{s}", 0, 1, -1, true);

    fLoc[kUMu3Sqr] = AddParameter("umu3", "|U_{#mu3}|^{2}",
                                  kUMu3SqrStart,
                                  kUMu3SqrEnd,
                                  kDeltaUMu3Sqr, true);

    fLoc[kDeltaMSqr32] = AddParameter("dmsq",
                                      "#Deltam^{2}_{32} (10^{-3} eV^{2})",
                                      kDeltaMSqrStart,
                                      kDeltaMSqrEnd,
                                      kDeltaDeltaMSqr, true);

    if(r.Get("FitUE3Sqr", tmpb) && tmpb){
      fLoc[kUE3Sqr] = AddParameter("ue3", "|U_{e3}|^{2}", 0, 1, -1, true);
    }

    return;

  case kSterileFractionTauNorm:
    //fLoc[kFs] = AddParameter("fs", "f_{s}", 0, 1, -1, true);

    fLoc[kUMu3Sqr] = AddParameter("umu3", "|U_{#mu3}|^{2}",
                                  kUMu3SqrStart,
                                  kUMu3SqrEnd,
                                  kDeltaUMu3Sqr, true);
    fLoc[kTauScale] = AddParameter("tauscale", "A_{#tau}",
                                  -0.1,
                                  2,
                                  0.01, true);
    fLoc[kDeltaMSqr32] = AddParameter("dmsq",
                                      "#Deltam^{2}_{32} (10^{-3} eV^{2})",
                                      kDeltaMSqrStart,
                                      kDeltaMSqrEnd,
                                      kDeltaDeltaMSqr, true);
    
    if(r.Get("FitUE3Sqr", tmpb) && tmpb){
      fLoc[kUE3Sqr] = AddParameter("ue3", "|U_{e3}|^{2}", 0, 1, -1, true);
    }
    return;

  case kDecay:
    fLoc[kTheta] = AddParameter("theta", "#theta", 0, TMath::Pi()/2, -1, true);

    fLoc[kAlpha] = AddParameter("alpha", "#alpha (GeV/km)", 0, 8e-3, -1, true);

    // Using a lower limit of zero gets us the pure-decay dchisq
    fLoc[kDeltaMSqr] = AddParameter("dmsq_32", "#Deltam^{2} (10^{-3} eV^{2})",0, 4.5, -1, true);

    return;

  case kDecoherence:
    fLoc[kTheta] = AddParameter("theta", "#theta", 0, TMath::Pi()/2, -1, true);

    // TODO - range?
    fLoc[kMu] = AddParameter("mu", "#mu", 0, 0.2, -1, true);

    return;
  }

  assert(0 && "Unknown oscillation model");
}


double NC::CoordinateConverter::
ChooseValue(NCType::EFitParam param, const NC::Fitter::CoordNDim& coords) const
{
  if(IsFit(param))
    return coords.at(fLoc[param]);
  else
    return fFixedVals[param];
}


void NC::CoordinateConverter::Prepare(const Registry& r, bool hush)
{
  int    tmpi;
  double tmpd;

  if(r.Get("DeltaCP3FlavorVal", tmpd)) fFixedVals[kDelta13]    = tmpd;
  if(r.Get("UE3SqrVal",         tmpd)) fFixedVals[kUE3Sqr]     = tmpd;
  if(r.Get("Theta13Val",        tmpd)) fFixedVals[kTheta13]    = tmpd;
  if(r.Get("Delta1Val",         tmpd)) fFixedVals[kDelta1]     = tmpd;
  if(r.Get("Delta2Val",         tmpd)) fFixedVals[kDelta2]     = tmpd;
  if(r.Get("Theta14Val",        tmpd)) fFixedVals[kTheta14]    = tmpd;
  if(r.Get("Theta24Val",        tmpd)) fFixedVals[kTheta24]    = tmpd;
  if(r.Get("Theta23Val",        tmpd)) fFixedVals[kTheta23]    = tmpd;


  if(r.Get("PrecScale",         tmpd)) fPrecScale         = tmpd;

  if(r.Get("OscillationModel",  tmpi)) fOscillationModel  = NCType::EOscModel(tmpi);

  PrepareOscillations(r);
  PrepareSystematics(r);

  if(!hush){
    MSG("NCCoordConv", Msg::kInfo) << "parameter names, labels, etc:" << endl;

    for(unsigned int i = 0; i < fFitParams.size(); ++i)
      MSG("NCCoordConv", Msg::kInfo) << fFitParams[i].ShortName() << "\t"
                                     << fFitParams[i].LatexName() << "\t"
                                     << fFitParams[i].Min() << " to "
                                     << fFitParams[i].Max() << " in "
                                     << fFitParams[i].Precision() << " with "
                                     << (fFitParams[i].UseLimits() ?
                                         "binding" : "non-binding")
                                     << " limits" << endl;
  }
}


NC::SystPars NC::CoordinateConverter::
SystParsFromCoordNDim(const CoordNDim& coords) const
{
  NC::SystPars ret;

  for(int n = 0; n < kNumParameters; ++n){
    if(!IsSystematic(EFitParam(n))) continue;

    if(IsFit(EFitParam(n))) ret.fShifts[n] = coords[fLoc[n]];
  }

  return ret;
}


NC::OscProb::OscPars*
NC::CoordinateConverter::OscParsFromCoordNDim(const CoordNDim& coords) const
{
  using namespace NC::OscProb;

  // This can happen if eg we didn't do a fit. It's still the callers'
  // fault, but this is more polite than the alternatives (crashing).
  if(coords.empty() && fOscillationModel != kNoOscillations) return 0;

  try{
    switch(fOscillationModel){
    case kThreeFlavor:{
      ThreeFlavor* pars = new ThreeFlavor;

      pars->Theta23() = ChooseValue(kTheta23, coords);
      pars->DeltaMSqr32() = coords.at(fLoc[kDeltaMSqr32])*1.e-3;

      pars->Theta13() = ChooseValue(kTheta13, coords);

      pars->Theta12() = 0.61; //SNO+KamLAND
      pars->DeltaMSqr12() = 7.59e-5;//SNO+KamLAND
      pars->RockDensity() = 2.7;//g/cc standard rock
      pars->Delta13() = fFixedVals[kDelta13];//cp violating phase
      pars->Hierarchy()  = 1;

      return pars;
    }

    case kNoOscillations:{
      return new NoOscillations;
    }

    case kFourFlavorGeneral:{
      FourFlavorGeneral* pars = new FourFlavorGeneral;

      pars->Theta13() = ChooseValue(kTheta13, coords);
      pars->Delta1() =  ChooseValue(kDelta1,  coords);
      pars->Delta2() =  ChooseValue(kDelta2,  coords);
      pars->Theta14() = ChooseValue(kTheta14, coords);

      pars->Theta23() = ChooseValue(kTheta23, coords);
      pars->Theta24() = ChooseValue(kTheta24, coords);
      pars->Theta34() = coords.at(fLoc[kTheta34]);
      pars->DeltaMSqr31() = coords.at(fLoc[kDeltaMSqr31])*1e-3;
      pars->DeltaMSqr41() = coords.at(fLoc[kDeltaMSqr41])*1e-3;
      pars->DeltaMSqr43() = coords.at(fLoc[kDeltaMSqr43])*1e-3;

      return pars;
    }

    case kFourFlavorDelta43IsBig:{
      FourFlavorDelta43IsBig* pars = new FourFlavorDelta43IsBig;

      pars->Theta13() = ChooseValue(kTheta13, coords);
      pars->Delta1() =  ChooseValue(kDelta1,  coords);
      pars->Delta2() =  ChooseValue(kDelta2,  coords);
      pars->Theta14() = ChooseValue(kTheta14, coords);

      pars->Theta23() = ChooseValue(kTheta23, coords);
      pars->Theta24() = ChooseValue(kTheta24, coords);
      pars->Theta34() = coords.at(fLoc[kTheta34]);
      pars->DeltaMSqr31() = coords.at(fLoc[kDeltaMSqr31])*1e-3;

      return pars;
    }

    case kFourFlavorDelta41Is0:{
      FourFlavorDelta41Is0* pars = new FourFlavorDelta41Is0;

      pars->Theta13() = ChooseValue(kTheta13, coords);

      pars->Theta23() = ChooseValue(kTheta23, coords);
      pars->Theta34() = coords.at(fLoc[kTheta34]);
      pars->DeltaMSqr31() = coords.at(fLoc[kDeltaMSqr31])*1e-3;

      return pars;
    }

    case kFourFlavorDelta43Is0:{
      FourFlavorDelta43Is0* pars = new FourFlavorDelta43Is0;

      pars->Theta13() = ChooseValue(kTheta13, coords);
      pars->Delta1() =  ChooseValue(kDelta1,  coords);
      pars->Theta14() = ChooseValue(kTheta14, coords);

      pars->Theta23() = ChooseValue(kTheta23, coords);
      pars->Theta24() = ChooseValue(kTheta24, coords);
      pars->DeltaMSqr31() = coords.at(fLoc[kDeltaMSqr31])*1e-3;

      return pars;
    }

    case kSterileFraction:{
      SterileFraction* pars = new SterileFraction;

      pars->UE3Sqr() = ChooseValue(kUE3Sqr, coords);

      pars->UMu3Sqr() = coords.at(fLoc[kUMu3Sqr]);
      pars->DeltaMSqr32() = coords.at(fLoc[kDeltaMSqr32])*1.e-3;
      pars->Fs() = coords.at(fLoc[kFs]);


      return pars;
    }

    case kDecay:{
      Decay* pars = new Decay;

      pars->Theta() = coords.at(fLoc[kTheta]);
      pars->Alpha() = coords.at(fLoc[kAlpha]);
      pars->DeltaMSqr() = coords.at(fLoc[kDeltaMSqr])*1e-3;

      return pars;
    }

    case kSterileFractionTauNorm:{
      SterileFractionTauNorm* pars = new SterileFractionTauNorm;
      
      pars->UMu3Sqr() = coords.at(fLoc[kUMu3Sqr]);
      pars->DeltaMSqr32() = coords.at(fLoc[kDeltaMSqr32])*1.e-3;
      pars->Fs() = 0;//coords.at(fLoc[kFs]);
      pars->TauScale() = coords.at(fLoc[kTauScale]);
      
      pars->UE3Sqr() = ChooseValue(kUE3Sqr, coords);
      
      return pars;
    }
      
    case kDecoherence:{
      Decoherence* pars = new Decoherence;
      
      pars->Theta() = coords.at(fLoc[kTheta]);
      pars->Mu() = coords.at(fLoc[kMu]);

      return pars;
    }

    } // end switch

    assert(0 && "Unknown oscillation model");
  } // end try

  catch(std::out_of_range& e){
    MSG("NCCoordConv", Msg::kError) << "Out of range error in "
                                    << "OscParsForCoordNDim, with model "
                                    << fOscillationModel << " and coords size "
                                    << coords.size() << ". Aborting." << endl;
    abort();
  }
}


CoordNDim NC::CoordinateConverter::
CoordNDimFromOscPars(const NC::OscProb::OscPars* pars) const
{
  assert(pars->OscillationModel() == fOscillationModel);

  CoordNDim ret;

  using namespace NCType;

  for(int n = 0; n < kNumParameters; ++n){
    const EFitParam par = EFitParam(n);
    if(IsFit(par) && !IsSystematic(par)){
      const unsigned int loc = FitterIndex(par);
      if(ret.size() <= loc) ret.resize(loc+1);
      double parval = pars->GetParameterValue(par);
      // Have to undo 1e3 factor in deltam values. TODO better
      if(par == kDeltaMSqr || par == kDeltaMSqr12 || par == kDeltaMSqr31 ||
         par == kDeltaMSqr32 || par == kDeltaMSqr41 || par == kDeltaMSqr43)
        parval *= 1e3;
      ret[loc] = parval;
    }
  }

  return ret;
}


int NC::CoordinateConverter::FitterIndex(EFitParam par) const
{
  assert(int(par) >= 0 && par < kNumParameters);
  assert(fLoc[par] != -1);
  return fLoc[par];
}


NCParameter NC::CoordinateConverter::ParameterForFitterIndex(int idx) const
{
  assert(idx >= 0 && idx < int(fFitParams.size()));

  return fFitParams[idx];
}


NCParameter NC::CoordinateConverter::
ParameterForFitParam(NCType::EFitParam par) const
{
  return ParameterForFitterIndex(FitterIndex(par));
}


bool NC::CoordinateConverter::IsFit(NCType::EFitParam par) const
{
  return fLoc[par] != -1;
}


int NC::CoordinateConverter::NumParameters() const
{
  int fit = 0;
  for(int n = 0; n < NCType::kNumParameters; ++n)
    if(IsFit(NCType::EFitParam(n))) ++fit;
  return fit;
}


bool NC::CoordinateConverter::IsSystematic(NCType::EFitParam par) const
{
  assert(par < NCType::kNumParameters);
  return par < NCType::kNumSystematicParameters;
}


int NC::CoordinateConverter::AddParameter(string shortName,
                                          string latexName,
                                          double min,
                                          double max,
                                          double prec,
                                          bool limit)
{
  assert(fPrecScale > 0);

  const double defaultPrec = .01;
  if(prec < 0) prec = TMath::Abs((max - min)*defaultPrec);

  fFitParams.push_back(NCParameter(shortName, latexName, min, max,
                                   prec/fPrecScale, limit));
  return fFitParams.size()-1;
}


int NC::CoordinateConverter::AddSystematicParameter(EFitParam p)
{
  const double sigma = kParams[p].sigma;
  return AddParameter(kParams[p].name.Data(), kParams[p].latexName.Data(),
                      -sigma, sigma, -1, false);
}


TH1* AxesForNCParameter(TString name, TString title,
                        NCParameter xp, TString yLabel)
{
 return new TH1F((name != "") ? TString::Format("%s_%s",
                                                 name.Data(),
                                                 xp.ShortName().c_str())
                               : "",
                  TString::Format("%s;%s;%s",
                                  title.Data(),
                                  xp.LatexName().c_str(),
                                  yLabel.Data()),
                  100, xp.Min(), xp.Max());
}


TH1* NC::CoordinateConverter::AxesForParameter(TString name,
                                               TString title,
                                               NCType::EFitParam x,
                                               TString yLabel) const
{
  const NCParameter xp = ParameterForFitParam(x);

  return AxesForNCParameter(name, title, xp, yLabel);
}


TH1* NC::CoordinateConverter::AxesForParameter(TString name,
                                               TString title,
                                               int x,
                                               TString yLabel) const
{
  const NCParameter xp = ParameterForFitterIndex(x);

  return AxesForNCParameter(name, title, xp, yLabel);
}


TH2* AxesForNCParameters(TString name, TString title,
                         NCParameter xp, NCParameter yp)
{
  return new TH2F((name != "") ? TString::Format("%s_%s_%s",
                                                 name.Data(),
                                                 xp.ShortName().c_str(),
                                                 yp.ShortName().c_str())
                               : "",
                  TString::Format("%s;%s;%s",
                                  title.Data(),
                                  xp.LatexName().c_str(),
                                  yp.LatexName().c_str()),
                  100, xp.Min(), xp.Max(),
                  100, yp.Min(), yp.Max());
}


TH2* NC::CoordinateConverter::AxesForParameters(TString name,
                                                TString title,
                                                NCType::EFitParam x,
                                                NCType::EFitParam y) const
{
  const NCParameter xp = ParameterForFitParam(x);
  const NCParameter yp = ParameterForFitParam(y);

  return AxesForNCParameters(name, title, xp, yp);
}


TH2* NC::CoordinateConverter::AxesForParameters(TString name,
                                                TString title,
                                                int x, int y) const
{
  const NCParameter xp = ParameterForFitterIndex(x);
  const NCParameter yp = ParameterForFitterIndex(y);

  return AxesForNCParameters(name, title, xp, yp);
}


vector<double> NC::CoordinateConverter::
VectorFromSystPars(const NC::SystPars& s) const
{
  vector<double> shift;

  for(int n = 0; n < kNumParameters; ++n){
    if(!IsSystematic(EFitParam(n))) continue;

    if(IsFit(EFitParam(n))) shift.push_back(s.fShifts[n]/kParams[n].sigma);
  }

  return shift;
}

---