UnivV2Rec.h

Go to the documentation of this file.

#ifndef __UnivV2Rec_h__
#define __UnivV2Rec_h__

#include <cmath>
#include <float.h>

#include <fwk/CentralConfig.h>
#include <utl/Branch.h>
#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerUnivRecData.h>
#include <evt/ShowerSimData.h>
#include <sevt/SEvent.h>

#include <utl/Minou.h>
#include <utl/Math.h>
#include <utl/GeometryUtilities.h>
#include <utl/PhysicalFunctions.h>
#include <utl/NormalDistribution.h>

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

namespace un2 {

  using namespace std;

  void SimpleReco(evt::Event& event);


  class un2Detector {

    private:

      /* */

    public:

      un2Detector(const int id,
          const sevt::Station& stat,
          const string detType,
          const double r,
          const double az,
          const double height,
          const double pftimeres,
          const double denseFlag) :
        Id(id),
        station(stat),
        type(detType),
        spRadius(r),
        spAzimuth(az),
        heightAboveSL(height),
        pfTimeResidual(pftimeres),
        isVirtual(denseFlag)
      {

        //INFO(boost::format("%s  r = %.2f m \t ψ = %.2f° \t h = %.1f m")%type %spRadius %Rad2Deg(spAzimuth) %heightAboveSL);

      };

      const unsigned int Id;
      const sevt::Station& station;

      const string type; // string
      const double spRadius; // meter
      const double spAzimuth; // rad
      const double heightAboveSL; // meter
      const double pfTimeResidual; // ns
      const bool isVirtual;

      double mcSpRadius = 0;
      double DX = 0;
      double signalUnit = 0;
      double avgAoP = 0;
      double nsBin = 0;
      double stationNSTime = 0;
      unsigned int signalStartSlot = 0;
      unsigned int signalStopSlot = 0;
      double t40RelToStart = 0;
      double Dt40 = 0;
      vector<double> referenceDt40; // {value, uncertainty} in ns
      vector<double> referenceDt0; // {value, uncertainty} in ns

      double weight = 1; /* artificially increases uncertainty, should be between 0 and 1 */

      double detBestFitXmax = 0;
      double detBestFitRmu = 0;
      double detBestFitT0 = 0;
      double detBestFitN = 0;
      double detChiSq = 0;

      double totalSignal;
      double referenceSignal;
      double uncertaintySignal;

      vector<double> integratedSignal = {0};
      vector<double> integratedNormalizedSignal = {0};

      double
      CalcUnit()
      {
        double unit = 0;

        if ("WCD" == type) {

          double tmp_AoP = 0;
          double tmp_unit = 0;
          int working_pmts = 0;
          for (unsigned int i=1; i<4; ++i) {
            try {
              tmp_AoP += station.GetPMT(i).GetRecData().GetAreaOverPeak();
              tmp_unit += station.GetPMT(i).GetRecData().GetVEMPeak() / station.GetPMT(i).GetRecData().GetVEMCharge();
              working_pmts += 1;
            } catch (utl::NonExistentComponentException& /*e*/) {
              // WARNING(boost::format("No signal in one WCD PMT"));
              continue;
            }

          }

          unit = tmp_unit / working_pmts;
          avgAoP = tmp_AoP / working_pmts;

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

          if (station.GetScintillatorPMT().HasCalibData() &&
              station.GetScintillatorPMT().HasRecData() &&
              station.GetScintillator().HasRecData() &&
              station.GetScintillator().HasMIPTrace()) {

            try {
              unit = station.GetScintillatorPMT().GetRecData().GetVEMPeak() /
                          station.GetScintillatorPMT().GetRecData().GetVEMCharge();
            } catch (utl::NonExistentComponentException& /*e*/) {
              // WARNING(boost::format("No signal in SSD PMT"));
              return 0;
            }

          } else {
            WARNING(boost::format("No SSD PMT data"));
            return 0;
          }

        }

        signalUnit = unit;

        return unit;
      }

      const utl::TraceD&
      GetTrace(sevt::StationConstants::SignalComponent traceType=sevt::StationConstants::eTotal)
      {

        if ("WCD" == type) {

          return station.GetVEMTrace(traceType);

        }

        if ("SSD" == type) {

          return station.GetScintillator().GetMIPTrace(traceType);

        } else {

          WARNING(type);
          WARNING("Detector type not found!");
          return station.GetVEMTrace(traceType);

        }

      }

      bool
      HasTrace(sevt::StationConstants::SignalComponent traceType)
      {

        if ("WCD" == type) {

          return station.HasVEMTrace(traceType);

        }

        if ("SSD" == type) {

          return station.GetScintillator().HasMIPTrace(traceType);

        } else {

          WARNING(type);
          WARNING("Detector type not found!");
          return station.HasVEMTrace(traceType);

        }

      }

    };



  class LDFFit : public utl::Minou::Base {

    private:

      vector<un2Detector> fun2Detectors;

      // meta parameters on event level
      const utl::Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("UniversalityFitterV2");
      string fModel = topB.GetChild("HadronicInteractionModel").GetDataString();
      double fXmax = 0;
      double fLogE = 0;
      double fTheta = 0;
      bool fIsReal = false;
      int fMonth = 3;

    public:

      LDFFit() :
                //name start  step   low   high  fixed
        Base({
              {"Xmax", 780.,   .1, 600., 1200., true},
              {"Rmu",    1.,  .01,   0.,    3., true},
              {"lgE",   19.,   .1,   18,   21., false}})

      {}

      void
      AddDetectorData(un2Detector d)
      {

        fun2Detectors.push_back(d);

      }

      void
      AddMetaData(double avgXmax, double logE, double theta, bool sim, int month)
      {

        fXmax     = avgXmax;
        fLogE     = logE;
        fTheta    = theta;
        fIsReal   = !sim;
        fMonth    = month;

      }

      double
      operator()(const vector<double>& p)
      const
      {

        double logLikelihood = 0;

        for(unsigned int d = 0; d < fun2Detectors.size(); ++d) {

          const bool useThrs = true;

          double ll = 0;

          const double predSignal = GetSignal(fModel, p[2], fTheta,
                fun2Detectors[d].spRadius, fun2Detectors[d].spAzimuth, fun2Detectors[d].heightAboveSL,
                p[0], p[1],
                fun2Detectors[d].type, fMonth, fIsReal, useThrs)[0];

          //cout << fun2Detectors[d].type << " " << fun2Detectors[d].totalSignal << " " << predSignal << endl;

          if (fun2Detectors[d].totalSignal > 5) {

            ll += -utl::Sqr((fun2Detectors[d].totalSignal-predSignal)/fun2Detectors[d].uncertaintySignal);

          } else {

            ll += LogarithmOfTruncatedGaussianPDF(predSignal, fun2Detectors[d].totalSignal, fun2Detectors[d].uncertaintySignal);

          }

          logLikelihood += ll * fun2Detectors[d].weight;

        }

        return (!std::isnan(logLikelihood)) ? -logLikelihood : std::numeric_limits<double>::max();

      }

  };


  class SignalRatioFit : public utl::Minou::Base {

    private:

      vector<vector<un2Detector>> fun2DetectorPairs;
      vector<un2Detector> fun2Detectors;

      const utl::Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("UniversalityFitterV2");
      string fModel = topB.GetChild("HadronicInteractionModel").GetDataString();
      double fXmax = 0;
      double fLogE = 0;
      double fTheta = 0;
      bool fIsReal = false;
      int fMonth = 3;

    public:

      SignalRatioFit() :
                //name start  step   low   high  fixed
        Base({
              {"Xmax", 780.,    .1, 600., 1200., false},
              {"Rmu",   1.0,   .01,   0.,    3., false},
              {"lgE",   19.,    .1,   18,   21., false}})

      {}

      void
      AddDetector(un2Detector d)
      {

        fun2Detectors.push_back(d);

      }

      int
      NumPairs()
      {

        return fun2DetectorPairs.size();

      }

      void
      GeneratePairs()
      {

        for (const auto& det : fun2Detectors) {

          if ("WCD" == det.type) {
            continue;
          }

          const unsigned int id = det.Id;

          for (const auto& matchDet : fun2Detectors) {

            if (id == matchDet.Id && "WCD" == matchDet.type) {

              fun2DetectorPairs.push_back({det, matchDet});

            }

          }

        }

        INFO(boost::format("Generated %i detector pairs for AugerPrime fit.") % fun2DetectorPairs.size());
        for (const auto& dets : fun2DetectorPairs) {

          INFO(boost::format("%i (%i) %.2f MIP / %.2f VEM") % dets.at(0).Id % dets.at(1).Id % dets.at(0).totalSignal % dets.at(1).totalSignal);

        }

      }

      void
      AddMetaData(double avgXmax, double logE, double theta, bool sim, int month)
      {

        fXmax     = avgXmax;
        fLogE     = logE;
        fTheta    = theta;
        fIsReal   = !sim;
        fMonth    = month;

      }

      double
      operator()(const vector<double>& p)
      const
      {

        if (2>fun2DetectorPairs.size()) {
          WARNING("Too few detector pairs generated for ratio fit!");
          return 0;
        }

        double logLikelihood = 0;
        for(unsigned int d = 0; d < fun2DetectorPairs.size(); ++d) {

          const bool useThrs = true;
          const auto& ssd = fun2DetectorPairs[d].at(0);
          const auto& wcd = fun2DetectorPairs[d].at(1);
          const double recSignalSSD = ssd.totalSignal;
          const double recSignalWCD = wcd.totalSignal;

          const double predSignalSSD = GetSignal(fModel, p[2], fTheta,
                ssd.spRadius, ssd.spAzimuth, ssd.heightAboveSL,
                p[0], p[1],
                ssd.type, fMonth, fIsReal, useThrs)[0];
          const double predSignalWCD = GetSignal(fModel, p[2], fTheta,
                wcd.spRadius, wcd.spAzimuth, wcd.heightAboveSL,
                p[0], p[1],
                wcd.type, fMonth, fIsReal, useThrs)[0];

          const double z = recSignalSSD / recSignalWCD;
          const double zPred = predSignalSSD / predSignalWCD;
          // const double uncRatio = ssd.uncertaintySignal/wcd.uncertaintySignal;

          // const double corr = 0.6; // rough estimate for signal residual correlation
          //
          // const double alpha = uncRatio * corr;
          // const double beta = uncRatio * std::sqrt(1-utl::Sqr(corr));
          //
          // const double prob = 1./M_PI * beta/(utl::Sqr(z-alpha) + utl::Sqr(beta));

          // actually it is not perfectly gaussian, but deriving the PDF is quite tricky
          // and it is almost gaussian
          // sigma approximated from uncorrelated central normal ratio distribution
          const double mu = z-zPred;
          const double var = (utl::Sqr(ssd.uncertaintySignal/recSignalSSD) +
                              utl::Sqr(wcd.uncertaintySignal/recSignalWCD));
          const double sigma = std::sqrt(var);
          const double logprob = -0.5 * utl::Sqr(mu/sigma) - log(sigma*utl::kSqrt2Pi);
          // cout << p[1] << "\t" << z << "\t" << prob << endl;

          logLikelihood += logprob;
          //
          // REALLY needed, unfortunately
          if (predSignalWCD > 5) {

            logLikelihood += - 0.05 * utl::Sqr((predSignalWCD-recSignalWCD)/wcd.uncertaintySignal) -
                               log(wcd.uncertaintySignal*utl::kSqrt2Pi);
            logLikelihood += - 0.05 * utl::Sqr((predSignalSSD-recSignalSSD)/ssd.uncertaintySignal) -
                               log(ssd.uncertaintySignal*utl::kSqrt2Pi);

          }

        }

        return (!std::isnan(logLikelihood)) ? -logLikelihood : std::numeric_limits<double>::max();

      }

  };


  class RmuXmaxFit : public utl::Minou::Base {

    private:

      vector<un2Detector> fun2Detectors;

      // meta parameters on event level
      const utl::Branch topB = fwk::CentralConfig::GetInstance()->GetTopBranch("UniversalityFitterV2");
      string fModel = topB.GetChild("HadronicInteractionModel").GetDataString();
      double fXmax = 0;
      double fLogE = 0;
      double fTheta = 0;
      bool fIsReal = false;
      int fMonth = 3;
      bool fNorm = false;
      double fPenalty = 0;

    public:

      RmuXmaxFit() :
                //name start  step   low   high  fixed
        Base({
              {"Xmax", 780.,    .1, 600., 1200., false},
              {"Rmu",   1.0,   .01,   0.,    3., false},
              {"lgE",   19.,    .1,   18,   21., false}})

      {}


      void
      AddDetectorData(un2Detector d)
      {

        fun2Detectors.push_back(d);

      }


      void
      AddMetaData(double avgXmax, double logE, double theta,
                  bool sim, int month, bool normalized, double penalty=0)
      {

        fXmax     = avgXmax;
        fLogE     = logE;
        fTheta    = theta;
        fIsReal   = !sim;
        fMonth    = month;
        fNorm     = normalized;
        fPenalty  = penalty;

      }


      double
      operator()(const vector<double>& p)
      const
      {

        double logLikelihood = 0;

        /* time information used here */
        for (unsigned int d = 0; d < fun2Detectors.size(); ++d) {

          double ll = 0;
          vector<double> timeVec    = {};
          vector<double> timeVecDet = {};

          // time shift allowed to counter trigger issues
          const double t0 = p[3+d];
          // use parametrized start time as a functino of Xmax
          //const double t0 = GetTimeQuantile(fTheta, fun2Detectors[d].spRadius,
                                            //fun2Detectors[d].spAzimuth, fun2Detectors[d].heightAboveSL,
                                            //p[0], fun2Detectors[d].type, "t0", fMonth, fIsReal)[0].at(0);

          for (unsigned int t = 0; t < fun2Detectors[d].integratedSignal.size(); ++t) {
            timeVec.push_back(t*fun2Detectors[d].nsBin);
            timeVecDet.push_back(t*fun2Detectors[d].nsBin - t0);
          }

          // allow for normalization
          //const double norm = (fNorm) ? p[2+2*d+1] : 0;

          //const vector<double> st = GetIntegratedTotalSignal(timeVecDet, "EPOS-LHC", fLogE, fTheta,
          const vector<double> st = GetIntegratedTotalSignalSum(timeVecDet, fModel, p[2], fTheta,
              fun2Detectors[d].spRadius, fun2Detectors[d].spAzimuth, fun2Detectors[d].heightAboveSL,
              p[0] /*xmax*/, p[1] /*rmu*/, fun2Detectors[d].type, fMonth, fIsReal, true,
              //fun2Detectors[d].pfTimeResidual, fun2Detectors[d].nsBin, "cdf", 0 [>component integer<],
              fun2Detectors[d].pfTimeResidual, fun2Detectors[d].nsBin,
              fun2Detectors[d].integratedSignal.back());

          const vector<double> stunc = GetIntegratedTotalSignal(timeVec, fModel, p[2], fTheta,
              fun2Detectors[d].spRadius, fun2Detectors[d].spAzimuth, fun2Detectors[d].heightAboveSL,
              fXmax /*xmax*/, 1. /*rmu*/, fun2Detectors[d].type, fMonth, fIsReal, true,
              fun2Detectors[d].pfTimeResidual, fun2Detectors[d].nsBin, "unc", 0,
              fun2Detectors[d].totalSignal);

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

            //if (timeVec[t] < 50.) // first couple of bins are just too volatile, so skip 100 ns
                //continue;

            // model is gaussian
            ll += -1 * utl::Sqr((fun2Detectors[d].integratedSignal[t] - st[t]) /
                    (sqrt(utl::Sqr(stunc[t]) + fun2Detectors[d].totalSignal)) *
                    fun2Detectors[d].weight);
            //cout << boost::format("ll: %.3f \t sig: %.1f \t n*st: %.1f \t unc: %.1f\n")
            // %ll %fun2Detectors[d].integratedSignal[t] %(norm * st[t]) %stunc[t];

          }

          if (!std::isnan(ll)) {

            // check for -nan cases
            logLikelihood += ll;

          } else {

            logLikelihood += -1e5;

          }

          // cout << p[0] << "\t"<< ll << "\t" << d << "/" << fun2Detectors.size() << endl;

        }

        /* signal information used here */
        for(unsigned int d = 0; d < fun2Detectors.size(); ++d) {

          double ll = 0;

          const double predSignal = GetSignal(fModel, p[2], fTheta,
                fun2Detectors[d].spRadius, fun2Detectors[d].spAzimuth, fun2Detectors[d].heightAboveSL,
                p[0], p[1],
                fun2Detectors[d].type, fMonth, fIsReal, true)[0];

          if (fun2Detectors[d].totalSignal > 5) {

            ll += -utl::Sqr((fun2Detectors[d].totalSignal-predSignal)/fun2Detectors[d].uncertaintySignal);

          } else {

            ll += LogarithmOfTruncatedGaussianPDF(predSignal, fun2Detectors[d].totalSignal, fun2Detectors[d].uncertaintySignal);

          }

          /* correlation information is used here */
          // Rµ Xmax correlation penalty
          // CAVEAT:
          // this penalty constrains forces the reconstructed Xmax value to follow the implemented
          // distribution. this is good if you want to constrain your fit but might be misleading
          // if you want to check what the reconstruction is doing with the trace shape information of the
          // stations only.
          // for real data the Xmax/Rmu correlation is different than from simulation, therefore
          // this information must be parsed here, if the penalty is on. currently this is done using
          // the calibrated lnRµ0 estimated from EPOS-LHC protons.
          // Remove this for photon searches!
          if (!std::isnan(ll)) {

            const double penaltyContribution = log(XmaxDistribution(p[1], p[0], p[2], fIsReal));
            logLikelihood += fPenalty * penaltyContribution * 50 * pow(10, fLogE-19);
            logLikelihood += ll * fun2Detectors[d].weight * 0.1;

          } else {

            logLikelihood += -1e5;

          }

          // if you change anything at the Xmax/Rmu correlation you have to adjust the multiplicator
          // wrt to the station likelihood. this is a magic number, unfortunately!  the exact
          // value however does not change the Xmax reconstruction much.
          // check output for estimating fudge factor
          //cout << p[0] << " " << p[1] << "\t"<< logLikelihood << "\t" << penaltyContribution << "\tratio:" <<
            //logLikelihood / penaltyContribution / fun2Detectors.size() << endl;
        }

        // check if nan, if so return maximum float number
        return (!std::isnan(logLikelihood)) ? -logLikelihood : std::numeric_limits<double>::max();

      }

  };

};

#endif