PredictionSterile.cxx

Go to the documentation of this file.

#include "NuXAna/Prediction/PredictionSterile.h"

#include <iostream>
#include <typeinfo>
#include <utility>

#include "OscLib/OscCalcSterile.h"

#include "NuXAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Core/LoadFromFile.h"

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

namespace ana
{
  //----------------------------------------------------------------------
  PredictionSterile::PredictionSterile(std::unique_ptr<ModularExtrapSterile> &&extrap)
    : PredictionExtrap(std::move(extrap))
      // This base class constructor initializes everything 
  {
  }

  //---------------------------------------------------------------------------
  PredictionSterile::PredictionSterile(PredictionExtrap&& load)
    : PredictionExtrap(std::move(load))
  {
  }

  //----------------------------------------------------------------------
  Spectrum PredictionSterile::PredictUnoscillated() const
  {
    // Use an osc::OscCalcSterile to avoid warnings/errors from PredictComponent
    // Just set to a 3 flavor calculator, with all oscillation parameters to 0
    osc::OscCalcSterile noosc;
    noosc.SetNFlavors(3);
    noosc.SetL(0);
    noosc.SetRho(0);
    noosc.SetAngle(1, 2, 0);
    noosc.SetAngle(1, 3, 0);
    noosc.SetAngle(2, 3, 0);
    noosc.SetDelta(1, 3, 0);
    noosc.SetDm(2, 0);
    noosc.SetDm(3, 0);
    return Predict(&noosc);
  }

  //----------------------------------------------------------------------
  Spectrum PredictionSterile::PredictComponent(osc::IOscCalc* calc,
                                               Flavors::Flavors_t flav,
                                               Current::Current_t curr,
                                               Sign::Sign_t sign) const
  {
    Spectrum ret = fExtrap->NCTotalComponent(); // Get binning
    ret.Clear();

    if(curr & Current::kCC){
      // Nothing new here
      ret += PredictionExtrap::PredictComponent(calc, flav, Current::kCC, sign);
    }
    if(curr & Current::kNC){
      osc::OscCalcSterile* calc_sterile = osc::DowncastToSterile(calc);

      if(calc_sterile) {
        // The PredictionSterile constructors all require a ModularExtrap Sterile,
        // so this should be safe unless someone else makes a change that breaks this
        ModularExtrapSterile* extrap = static_cast<ModularExtrapSterile*>(fExtrap.get());

        // Oscillate NCs as <start flavor> to any active flavor,
        // then multiply by the proportion of NCs that originate as <start flavor>
        if(flav & Flavors::kNuEToNuE   && sign & Sign::kNu)
          ret += extrap->OscNCComponent().Oscillated(calc_sterile, +12, 0)
                 *extrap->NCNueProportion();
        if(flav & Flavors::kNuEToNuE   && sign & Sign::kAntiNu)
          ret += extrap->OscNCComponent().Oscillated(calc_sterile, -12, 0)
                 *extrap->NCAntiNueProportion();
  
        if(flav & Flavors::kNuMuToNuMu && sign & Sign::kNu)
          ret += extrap->OscNCComponent().Oscillated(calc_sterile, +14, 0)
                 *extrap->NCNumuProportion();
        if(flav & Flavors::kNuMuToNuMu && sign & Sign::kAntiNu)
          ret += extrap->OscNCComponent().Oscillated(calc_sterile, -14, 0)
                 *extrap->NCAntiNumuProportion();
      }
      else {
        std::cerr << "Input calculator was not of type OscCalcSterile." << std::endl
                  << "Reverting to 3 Flavor oscillation with no NC oscillations." << std::endl;

        ret += PredictionExtrap::PredictComponent(calc, flav, Current::kNC, sign);
      }
    }

    return ret;
  }

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

    PredictionExtrap::SaveTo(dir, name);

    dir = (TDirectory*)dir->GetDirectory(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("PredictionSterile").Write("type", TObject::kOverwrite);

    delete dir;

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<PredictionSterile> PredictionSterile::LoadFrom(TDirectory* dir, const std::string& name)
  {
    TObjString* ptag = (TObjString*)dir->Get((name+"/type").c_str());
    assert(ptag);

    const TString tag = ptag->GetString();
    assert(tag.CompareTo("PredictionSterile") == 0);
    delete ptag;

    std::unique_ptr<PredictionSterile> ret(
                                           new PredictionSterile( std::move( *(PredictionExtrap::LoadFrom(dir, name).release()) ) )
    );

    return ret;
  }
}