GetSignal.cc

Go to the documentation of this file.

// OFFLINE
#include <fwk/CentralConfig.h>
#include <utl/Branch.h>
#include <utl/AugerUnits.h>
#include <det/Detector.h>

// LOCAL
#include <src/UniversalityConfig.h>
#include <src/AtmosphereUtilities.h>
#include <src/GetSignal.h>
#include <src/Functions.h>
#include <src/PowerConfig.h>


#include <utl/PhysicalFunctions.h>

using namespace std;


namespace un2 {


  vector<double>
  GetSignal(string model,
            double lgE,
            double zenith,
            double r,
            double psi,
            double height,
            double Xmax,
            double Rmu,
            string detector,
            int month,
            bool realAtmosphere,
            bool useThreshold) // this flag enables the threshold function needed for muonic components
  {
    const UniversalityConfig& config = PowerConfig::GetInstance();
    fwk::CentralConfig& cc = *fwk::CentralConfig::GetInstance(/*config.GetBootstrap()*/);

    const utl::Branch topBSig = cc.GetTopBranch("SignalModel").GetChild(model);
    const utl::Branch topBAtm = cc.GetTopBranch("AtmosphereParameters");

    // "un2" is the parameterised malargue atmosphere,
    // bariloche is equivalent to the CORSIKA atmosphere
    const string atmVersion = (realAtmosphere) ? "un2" : "bariloche";

    const map<unsigned int, string>& monthsDict = config.GetMonthsMapReverse();

    const double Xvert =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("Xvert").Get<double>();
    const double hs =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs").Get<double>();
    const double hs2 =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs2").Get<double>();

    const double DX = AnalyticDX(r, psi, zenith, Xmax, Xvert, hs, hs2, height);
    //cout << "DX " << DX << endl;

    // get reference signal
    //                       tot em mu emmu emhd
    vector<double> Signal = { 0., 0., 0., 0., 0. };

    for (const auto& comp : config.GetReducedComponents()) {


      // gamma
      double gamma = 0;

      {
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("gamma").Get<vector<double>>();

        for (unsigned int p = 0; p < a_pars.size(); ++p)
          gamma += pow(r, p) * a_pars[p];
      }


      double Sref;
      {
        const double n =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("n").Get<double>();
        const double rG =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("rG").Get<double>();
        const double s =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("s").Get<double>();

        Sref = NKG(r, n, rG, s);

      }

      if (useThreshold) {
        const double cThrsh =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("cThrsh").Get<double>();
        const double fThrsh =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("fThrsh").Get<double>();
        const double bg =
          topBSig.GetChild(detector).GetChild(comp).GetChild("Sref").GetChild("bg").Get<double>();

        Sref = ActivationFunctionBackground(Sref, fThrsh, cThrsh, bg);

      }

      // DeltaXmax
      double DeltaXmax = 0;
      {
        // only the constant and linear terms (a0, a1) are used, however this could be in need
        // for extension at some point
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("DeltaXmax").Get<vector<double>>();

        for (unsigned int p = 0; p < a_pars.size(); ++p)
          DeltaXmax += pow(r, p) * a_pars[p];  // TODO: horner
      }

      // DeltaX0
      double DeltaX0 = 0;
      {
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("DeltaX0").Get<vector<double>>();

        for (unsigned int p = 0; p < a_pars.size(); ++p)
          DeltaX0 += pow(r, p) * a_pars[p];
      }


      // DeltaXref
      // this is fixed, might be updated later!
      double DeltaXref;
      {
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("DeltaXref").Get<vector<double>>();
        // must be constant!
        DeltaXref = a_pars[0];
      }

      // lambda0
      double lambda0 = 0;
      {
        // only the constant and linear terms (a0, a1) are used, however this could be in need
        // for extension at some point
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("lambda0").Get<vector<double>>();

        for (unsigned int p = 0; p < a_pars.size(); ++p)
          lambda0 += pow(r, p) * a_pars[p];
      }

      // lambda1
      double lambda1 = 0;
      {
        // only the constant and linear terms (a0, a1) are used, however this could be in need
        // for extension at some point
        const vector<double> a_pars =
          topBSig.GetChild(detector).GetChild(comp).GetChild("lambda1").Get<vector<double>>();

        for (unsigned int p = 0; p < a_pars.size(); ++p)
          lambda1 += pow(r, p) * a_pars[p];
      }

      // final signal prediction for each component
      double Spred =
        ModifiedGH(DX, Sref, lgE, gamma, DeltaX0, DeltaXmax, DeltaXref, lambda0, 0/*lambda1*/, 0);

      const map<string, float> psiconsts = {{"em", 0.01}, {"mu", 0.05}, {"emmu", 0.05}, {"emhd", 0.02}};
      const double aziCorrection = exp(psiconsts.at(comp) * r/1000 * (DX-DeltaX0)/lambda0 * cos(psi) * sin(zenith));

      Spred *= aziCorrection;

      double hadrComponentFactor = 1.;
      if (detector=="SSD")
        hadrComponentFactor = 1.1;
      if (detector=="WCD")
        hadrComponentFactor = 1.;


      // the dependency of Rmu is implemented HERE:
      if (comp == "em") {
        Signal[0] += Spred;
        Signal[1]  = Spred;
      } else if (comp == "mu") {
        Signal[0] += Rmu * Spred;
        Signal[2]  = Rmu * Spred;
      } else if (comp == "emmu") {
        Signal[0] += Rmu * Spred;
        Signal[3]  = Rmu * Spred;
      } else if (comp == "emhd") {
        if ((1  + hadrComponentFactor*(Rmu - 1)) > 0) {
          Signal[0] += (1 + hadrComponentFactor*(Rmu -1)) * Spred;
          Signal[4]  = (1 + hadrComponentFactor*(Rmu -1)) * Spred;
        } else {
          Signal[0] += 0; // jut for completeness, adding "zero" if  1.2*Rmu is smaller than 0.2
          Signal[4] += 0;
        }
      }
    }

    return Signal;
  }


  vector<vector<double>>
  GetTimeQuantile(double zenith,
                  double r,
                  double psi,
                  double height,
                  double Xmax,
                  string detector,
                  string quantile,
                  int month,
                  bool realAtmosphere)
  // returns station start time !!! RELATIVE TO PLANE SHOWER FRONT !!! dependent on DX/Xmax
  {
    const UniversalityConfig& config = PowerConfig::GetInstance();
    fwk::CentralConfig& cc = *fwk::CentralConfig::GetInstance(/*config.GetBootstrap()*/);

    const utl::Branch topBTimeModel = cc.GetTopBranch("TimeQuantileModel");
    const utl::Branch topBAtm = cc.GetTopBranch("AtmosphereParameters");

    const string atmVersion = (realAtmosphere) ? "un2" : "bariloche";

    const map<unsigned int, string>& monthsDict = config.GetMonthsMapReverse();
    const double Xvert =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("Xvert").Get<double>();
    const double hs =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs").Get<double>();
    const double hs2 =
      topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs2").Get<double>();

    const double DX = AnalyticDX(r, psi, zenith, Xmax, Xvert, hs, hs2, height);

    vector<vector<double>> vals = {};

    for (auto& comp : config.GetComponents()) {

      //                          c ddx t0 sigma
      const vector<string> tags = { "c", "ddx", "t0", "sigma" }; // k
      vector<double> function_parameters = { 0., 0., 0., 0. }; // k

      for (unsigned short int k = 0; k < function_parameters.size(); ++k) {

        const vector<string> a_pars = { "a0", "a1", "a2", "a3", "a4" }; // i
        vector<double> a_vals = { 0,    0,    0,    0,    0 }; // i

        for (short unsigned int i = 0; i < a_pars.size(); ++i) {

          vector<double> aa_pars; // j
          try {
            topBTimeModel.GetChild(comp).GetChild(detector).GetChild(quantile).GetChild(tags[k]).GetChild(a_pars[i]).GetData(aa_pars);
          } catch (utl::XMLParseException& /*e*/) {
            ERROR(boost::format("SPECIFIED TIME QUANTILE NOT VALID: %s\nUSING t00 instead!") % quantile);
            topBTimeModel.GetChild(detector).GetChild("t00").GetChild(tags[k]).GetChild(a_pars[i]).GetData(aa_pars);
          }

          for (short unsigned int j = 0; j < aa_pars.size(); ++j) {
            // ai = aa0 + aa1 sin(theta) + aa2 sin^2(theta)
            a_vals[i] += aa_pars[j] * pow(sin(zenith), j);  // TODO: horner
          }

          // c = a0 + a1 r + a2 r^2 + a3 r^3 + a4 r^4
          function_parameters[k] += a_vals[i] * pow(r, i);

        }

      }

      const double c = function_parameters[0];
      const double ddx = function_parameters[1];
      const double t0 = function_parameters[2];

      const double sigma = function_parameters[3];
      const double t = SignalStartTimePlaneFront(r, psi, zenith, DX + ddx, Xvert, hs, hs2) / c + t0;

      vals.push_back({ t, sigma });

    }

    return vals;
  }


  vector<vector<double>>
  GetTraceParameters(double zenith,
                     double r,
                     double psi,
                     double height,
                     double Xmax,
                     string detector,
                     int month,
                     bool realAtmosphere,
                     double planeFrontTimeResidual, // if parsed, quantile model is used!
                     double nsBin)
  {
    const UniversalityConfig& config = PowerConfig::GetInstance();
    fwk::CentralConfig& cc = *fwk::CentralConfig::GetInstance(/*config.GetBootstrap()*/);

    vector<vector<double>> traceParameters = {};

    const vector<vector<double>> t40RelToPF = \
         GetTimeQuantile(zenith, r, psi, height, Xmax, detector, "t40", month, realAtmosphere);

    vector<vector<double>> signalTraces(5);
    vector<vector<double>> signalCDFs(5);
    vector<vector<double>> traceMoments(5);

    const utl::Branch topBTraceModel = cc.GetTopBranch("TimeTraceModel");
    const utl::Branch topBAtm = cc.GetTopBranch("AtmosphereParameters");

    const string atmVersion = (realAtmosphere) ? "un2" : "bariloche";
    const map<unsigned int, string>& monthsDict = config.GetMonthsMapReverse();
    double Xvert;
    double hs;
    double hs2;
    topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("Xvert").GetData(Xvert);
    topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs").GetData(hs);
    topBAtm.GetChild(atmVersion).GetChild(monthsDict.at(month)).GetChild("hs2").GetData(hs2);

    const double DX = AnalyticDX(r, psi, zenith, Xmax, Xvert, hs, hs2, height);

    unsigned int componentNumber = 0;

    for (const auto& comp : config.GetComponents()) {

      vector<double> componentTrace;
      vector<double> componentCDF;
      vector<double> componentMoments;

      map<string, map<string, double>> readout = { {"exp", {{"a", 0}, {"b", 0}} },
                                                   {"std", {{"a", 0}, {"b", 0}} } };

      for (const auto& moment : {"exp", "std"}) {

        for (const auto& linPar : {"a", "b"}) {

          vector<double> ai_vals = {0, 0, 0, 0, 0};
          const vector<string> ai = {"a0", "a1", "a2", "a3", "a4"};

          for (unsigned int i=0; i < ai.size(); ++i) {

            vector<double> aii_vals;
            topBTraceModel.GetChild(detector).GetChild(comp).GetChild(moment).GetChild(linPar).GetChild(ai[i]).GetData(aii_vals);

            for (unsigned int j=0; j < aii_vals.size(); ++j) {

              ai_vals[i] += aii_vals[j] * pow(sin(zenith), j);

            }

          }

          for (unsigned int j=0; j < ai_vals.size(); ++j) {

            readout[moment][linPar] += ai_vals[j] * pow(r, j);

          }

        }

      }

      const double ExpA = readout["exp"]["a"];
      const double ExpB = readout["exp"]["b"];
      const double StdA = readout["std"]["a"];
      const double StdB = readout["std"]["b"];

      // 100ns is kind of a natural limit for these two quantities
      const double traceExpectationValue = (ExpA * DX + ExpB < 100) ? 100 : ExpA * DX + ExpB;
      const double traceStandardDeviation = (StdA * DX + StdB < 100) ? 100 : StdA * DX + StdB;

      traceMoments[componentNumber] = {traceExpectationValue, traceStandardDeviation};

      componentMoments = {traceExpectationValue, traceStandardDeviation};

      const double t40 = t40RelToPF[componentNumber][0] + planeFrontTimeResidual;

      const double sigmaLogNormalPar = LogNormalSigmaFromExpAndStd(traceExpectationValue, traceStandardDeviation, nsBin);
      const double muLogNormalParFromt40 = LogNormalMuFromTimeQuantileAndSigma(t40, 0.4, sigmaLogNormalPar, nsBin);

      // sometimes the planeFrontTimeResidual is far off the numerically allowed value, then we have to use
      // mu from the direct parameterization. this usually occurs only for distant stations
      const double muLogNormalPar = (-10 == planeFrontTimeResidual || muLogNormalParFromt40 != muLogNormalParFromt40) ?
        LogNormalMuFromExpAndStd(traceExpectationValue, traceStandardDeviation, nsBin) :
        muLogNormalParFromt40;

      traceParameters.push_back({muLogNormalPar, sigmaLogNormalPar});

      }

      return traceParameters;

    }

  vector<double>
  GetIntegratedTotalSignal(vector<double> t,
                           string model,
                           double lgE,
                           double zenith,
                           double r,
                           double psi,
                           double height,
                           double Xmax,
                           double Rmu,
                           string detector,
                           int month,
                           bool realAtmosphere,
                           bool useThreshold,
                           double planeFrontTimeResidual,
                           double nsBin,
                           string returnValue,
                           int component,
                           double normalized)
  {

    vector<double> Sig = {};
    vector<double> SigUncertainty = {};

    const double totalSignal = GetSignal(model, lgE, zenith, r, psi, height, Xmax, Rmu, detector, month, realAtmosphere, useThreshold)[component];
    const double norm = (normalized != 0) ? totalSignal/normalized : 1.;
    const vector<double> traceParameters = GetTraceParameters(zenith, r, psi, height, Xmax, detector, month, realAtmosphere, planeFrontTimeResidual, nsBin)[component];
    const double mu = traceParameters[0];
    const double sigma = traceParameters[1];

    for (unsigned int i = 0; i < t.size(); ++i) {

      Sig.push_back(LogNormalCDF(t[i]/nsBin, mu, sigma) * totalSignal / norm);

    }

    if (returnValue == "cdf" || returnValue == "integrated" || returnValue == "CDF" || returnValue == "total")
      return Sig;

    const double timeUncertaintyRaw = GetTimeQuantile(zenith, r, psi, height, Xmax, detector, "t40", month, realAtmosphere)[component][1];
    const double timeUncertainty = (timeUncertaintyRaw > 10) ? 10 : timeUncertaintyRaw; // dear darko, please trust me, dont change this number

    const double baseline = 1; // and some extra just to be sure

    for (unsigned int i = 0; i < t.size(); ++i) {

      // maybe replace with inline-option?
      double stdSigRaw = 0.;
      if (detector == "WCD") {

        utl::wcd::ESignalVarianceModel stdWCDmodel = utl::wcd::eGAP2014_035;
        //stdSigRaw = utl::wcd::SignalUncertainty(stdWCDmodel, cos(zenith), Sig[i]);
        stdSigRaw = LogNormalCDF(t[i]/nsBin, mu, sigma)*utl::wcd::SignalUncertainty(stdWCDmodel, cos(zenith), totalSignal);

      } else if (detector == "SSD") {

        utl::ssd::ESignalVarianceModel stdSSDmodel = utl::ssd::eGAP2019_000;
        //stdSigRaw = utl::ssd::SignalUncertainty(stdSSDmodel, cos(zenith), Sig[i]);
        stdSigRaw = LogNormalCDF(t[i]/nsBin, mu, sigma)*utl::ssd::SignalUncertainty(stdSSDmodel, cos(zenith), totalSignal);

      }

      const double stdSig = (normalized != 0) ? 1. : stdSigRaw; // if trace normalization is applied, no total signal estiamte should be used

      SigUncertainty.push_back(sqrt(utl::Sqr(LogNormalPDF(t[i]/nsBin, mu, sigma) * timeUncertainty * totalSignal /* Sig[i]*/) +  utl::Sqr(stdSig) + utl::Sqr(baseline)) / norm);

    }

    if (returnValue == "std" || returnValue == "unc" || returnValue == "STD" || returnValue == "error") {

      return SigUncertainty;

    } else {

      ERROR("GetIntegratedTotalSignal reached end of function without defined return value!");
      ERROR("Define returnValue! Was set to:");
      ERROR(returnValue);
      return {0};

    }

  }

  vector<double>
  GetIntegratedTotalSignalSum(vector<double> t,
                              string model,
                              double lgE,
                              double zenith,
                              double r,
                              double psi,
                              double height,
                              double Xmax,
                              double Rmu,
                              string detector,
                              int month,
                              bool realAtmosphere,
                              bool useThreshold,
                              double planeFrontTimeResidual,
                              double nsBin,
                              double normalized
                              )

  {
    // other than GetIntegratedTotalSignal, this calculated the integrated signals as the
    // sum of the components. This function hast other usecases than GetIntegratedSignal!

    vector<double> sumCompSignals;
    vector<vector<double>> compSignals;
    for (unsigned short int i=1; i<5; ++i /*the four components*/) {

      compSignals.push_back(GetIntegratedTotalSignal(t, model, lgE, zenith, r, psi, height, Xmax, Rmu,
            detector, month, realAtmosphere, useThreshold, planeFrontTimeResidual, nsBin, "cdf", i /*comp*/));

    }

    for (unsigned int i=0; i<t.size(); ++i) {

      double binSum = 0;

      for (unsigned int c=0; c < compSignals.size(); ++c) {

        binSum += compSignals[c][i];

      }

      sumCompSignals.push_back(binSum);

    }

    if (normalized != 0) {

      for (unsigned int i=0; i<sumCompSignals.size(); ++i) {

        sumCompSignals[i] = sumCompSignals[i] * normalized / sumCompSignals[sumCompSignals.size()-1];

      }

    }

    return sumCompSignals;

  }

}