Risetime1000LLL.cc

Go to the documentation of this file.

// File "Risetime1000LLL.cc"
//
// This module calculates the risetime for an event.  Typically this
// is a value interpolated from a fit to the individual risetimes of
// each candidate station in the shower plane.  An option to
// recalculate the risetime of each station is provided.
//
// D.Barnhill (August 10, 2005)
// Taking M. Healy's original code and modifying a few things to make
// it more simple and reliable
//
// M.Healy (Sep 18, 2006): Further modifing the code to allow the
// error parameterization from leeds.  Also commenting out the linear
// fit for risetime to eliminate events that are on the cusp in terms
// of quality.
//
// M.Healy (Oct 9, 2006): Export the chi2 of the risetime fit for
// collection by subsiquent modules.
//
// M.Healy (Nov 21, 2006): Fixed a bug related to the inclusion of
// non-canidate stations from the event selection.
//
// M.Healy (Nov 29, 2006): Major upgrade to eliminate vestigial code
// and bring the module up to distribution quality standards.  Version
// to be submitted to the offline SVN shortly.
//
// M.Healy (Dec 1, 2006): Used the major upgrade to add functionality
// that will interface with the ADST produced by Karlsruhe.  This will
// allow browsing the fit results in the ADST viewer.
//
// M.Healy (Dec 14, 2006): Bug fixes, and changes to compile with
// version v2r2p3 of the offline.
//
// M.Healy (Jan 8, 2007): Strip out code to call RecShower->GetZenith()
// because it is not implimented.  Reverted to old method for getting
// the shower zenith angle.
//
// M.Healy (Jan 17, 2007): Eliminated the streaming of the TF1 object
// in the fit results.  This was not working with the ADST writer so
// resorted to the use of two doubles.
//
// M.Healy (May 15, 2007): The module is changed to use variables in
// the event now.  This is intended to replace the static variables
// that are currently in use.  For now I leave much of the old code as
// is so that the changes are incrimental.  I don't want to break
// anything just before a release.
//
// Author:    M.D.Healy, D.Barnhill
//
// C. Jarne  (2014/02) New asymmetry correction.
// C. Jarne  (2014/02) Modification in  uncertanty estimation according
// to gap2013-079 (now correlation term included).
//
// Modified: Feb 2014, C. Jarne
// Created:   Sep 24, 2004
// Modified:  May 15, 2007
//********************************************************************

// Standard c++ headers
#include <sstream>
#include <cmath>
#include <vector>

// Offline headers
#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/SEvent.h>
#include <sevt/Header.h>
#include <sevt/Station.h>
#include <sevt/StationRecData.h>
#include <sevt/StationConstants.h>
#include <sevt/PMT.h>
#include <sevt/PMTRecData.h>

#include <det/Detector.h>

#include <utl/Reader.h>
#include <utl/AugerUnits.h>
#include <utl/ErrorLogger.h>
#include <utl/StringCompare.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>

// ROOT headers
#include <TFormula.h>
#include <TMatrixD.h>
#include <TVirtualFitter.h>

// Headers for this module
#include "Risetime1000LLL.h"

using namespace std;
using namespace fwk;
using namespace utl;
using namespace Risetime1000;
using std::ostringstream;
using std::vector;


double Risetime1000LLL::fgRisetime = -1;
double Risetime1000LLL::fgRisetimeError = -1;
double Risetime1000LLL::fgRisetimeReducedChi2 = -1;
RisetimeResultsLLL* Risetime1000LLL::fgRisetimeResults = nullptr;


//#define DEBUG_Risetime1000LLL


Risetime1000LLL::~Risetime1000LLL()
{
  delete fgRisetimeResults;
  fgRisetimeResults = nullptr;
}


VModule::ResultFlag
Risetime1000LLL::Init()
{
  std::ostringstream info;
  info << "Configuring the FADC Parameters module.\n"
          "***************************************\n";
  auto topB = CentralConfig::GetInstance()->GetTopBranch("Risetime1000");

  {
    auto recalculateRiseTimeB = topB.GetChild("RecalculateRisetime");
    recalculateRiseTimeB.GetChild("Recalculate").GetData(fForceRiseTimeRecalculation);
    recalculateRiseTimeB.GetChild("RiseTimeStartFraction").GetData(fRiseTimeStartFraction);
    recalculateRiseTimeB.GetChild("RiseTimeStopFraction").GetData(fRiseTimeStopFraction);
    if (fForceRiseTimeRecalculation) {
      info << "Going to recalculate risetimes for stations.\n"
              "Pulse RiseTime defined as:\n"
           << fRiseTimeStartFraction*100 << "% to " << fRiseTimeStopFraction*100
           << "% of total signal.\n";
    } else
      info << "Using existing station risetimes (from SdCalibrator).\n";
  }

  // Risetime Fit
  {
    info << "\nConfiguring the risetime fit routines.\n";
    auto risetimeFitB = topB.GetChild("RisetimeFit");
    risetimeFitB.GetChild("Fit").GetData(fDoRiseTimeFit);
    if (fDoRiseTimeFit) {
      info << "Fitting the risetime for the event is enabled.\n";
      risetimeFitB.GetChild("MinimumSignal").GetData(fMinimumSignalForRiseTimeFit);
      risetimeFitB.GetChild("MinimumDistance").GetData(fMinimumDistanceForRiseTimeFit);
      risetimeFitB.GetChild("MaximumDistance").GetData(fMaximumDistanceForRiseTimeFit);
      risetimeFitB.GetChild("RisetimeEvaluatedAt").GetData(fDistanceToInterpolateFitResult);
      info << "\tMinimum Signal " << fMinimumSignalForRiseTimeFit << " VEM\n"
           << "\tMinimum Distance " << fMinimumDistanceForRiseTimeFit/m << " m\n"
           << "\tMaximum Distance " << fMaximumDistanceForRiseTimeFit/m << " m\n"
           << "\tRisetime Distance " << fDistanceToInterpolateFitResult/m << " m\n";

      risetimeFitB.GetChild("IncludeSaturated").GetData(fIncludeSaturatedForRiseTimeFit);
      if (fIncludeSaturatedForRiseTimeFit)
        info << "\tIncluding saturated stations.\n";
      else
        info << "\tNot including saturated stations.\n";

      risetimeFitB.GetChild("WeightAsFunctionDistance").GetData(fRTWeightingFunction);
      info << "\tWeighting Function: " << fRTWeightingFunction << '\n';
    } else
      info << "Fitting the risetime for the event is disabled.\n";

    fRTWeights = new TFormula("RiseTimeWeights", fRTWeightingFunction.c_str());
    if (fRTWeights->IsZombie()) {
      ostringstream err;
      err << "Weight formula defined in Risetime1000LLL.xml contains an error!\n"
             "Current formula is: weight=" << fRTWeightingFunction << "\n"
             "This formula does not conform to the ROOT formula class guidelines.";
      ERROR(err);
      delete fRTWeights;
      fRTWeights = nullptr;
      return eFailure;
    }
  }

  topB.GetChild("UseDLECorrectedTrace").GetData(fUseDLECorrectedTrace);

  if (fUseDLECorrectedTrace)
    fComponent = sevt::StationConstants::eDirectLightSubtracted;
  else
    fComponent = sevt::StationConstants::eTotal;

  info <<"***************************************\n"
         "FADC Parameters module initialized.";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
Risetime1000LLL::Run(evt::Event& event)
{
  std::ostringstream info;

  fgRisetime = -1;
  fgRisetimeError = -1;
  fgRisetimeReducedChi2 = -1;
  delete fgRisetimeResults;
  fgRisetimeResults = nullptr;
  fRejectedStations = 0;

  if (!event.HasSEvent()) {
    info << "There is no SEvent";
    INFO(info);
    return eSuccess;
  }

  if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {

    const auto& sRecShower = event.GetRecShower().GetSRecShower();

    const auto& baryCS = sRecShower.GetBarycenterCoordinateSystem();
    const double sec = 1 / sRecShower.GetAxis().GetCosTheta(baryCS);

    // Added a correction to the MeanRiseTime for shower asymmetry as
    // suggested by the Leeds group in a private communication.
    // t1/2(correct) = t1/2 - g cos(zeta)
    // g = alpha + gamma r^2
    // alpha = -66.61 + 95.13 sec(theta) - 30.73 sec^2(theta)
    // gamma = -0.0009721 + 0.001993 sec(theta) - 0.001259 sec^2(theta) + 0.0002546 sec^3(theta)
    //
    // r is the distance to the shower plane (meters), theta is the zenith angle,
    // and zeta is the angle to a station with the zero direction defined
    // to be under the shower axis.
    // M.Healy August 8, 2006
    //const double alpha = -66.61 + secZenith*(95.13 - 30.73*secZenith);
    //const double gamma = -0.0009721 + secZenith*(0.001993 + secZenith*(-0.001259 + 0.0002546*secZenith));
    // New asymmetry correction (2014/01 C. Jarne)
    const double alpha = 96.73 + sec*(-282.40 + sec*(241.80 - sec*62.61));
    const double gamma = -0.0009572 + sec*(0.002068 + sec*(-0.001362 + sec*0.0002861));

    fgRisetimeResults = new RisetimeResultsLLL;

    auto& sEvent = event.GetSEvent();

    for (auto& currentStation : sEvent.StationsRange()) {
      if (!currentStation.HasVEMTrace(fComponent) || !currentStation.HasRecData())
        continue;

      unsigned short rejectCode = 0;

      // Setting the above variables
      auto& sRecData = currentStation.GetRecData();
      double stationRisetime = sRecData.GetRiseTime();

      if (fForceRiseTimeRecalculation || stationRisetime <= 0)
        stationRisetime = RecalculateRiseTime(currentStation);

      if (stationRisetime <= 0)
        continue;

      const double stationDistance = sRecData.GetSPDistance();
      const double stationDistanceError = sRecData.GetSPDistanceError();
      const double g = alpha + gamma * Sqr(stationDistance);
      // This gets the asymmetry angle as set by the reconstruction module.
      // Beware of custom reconstructions that do not set this information.
      const double zeta = sRecData.GetAzimuthShowerPlane();
      stationRisetime -= g * std::cos(zeta);
      const double signal =
        fUseDLECorrectedTrace ? sRecData.GetTotalSignalCleaned() : sRecData.GetTotalSignal();
      // slight mistake in error propagation here
      const double stationRisetimeError = fRTWeights->Eval(stationDistance/m, sec, signal);

      if (signal < fMinimumSignalForRiseTimeFit)
        rejectCode |= sevt::StationRecData::eLowSignal;
      if (!currentStation.IsCandidate())
        rejectCode |= sevt::StationRecData::eNotCandidate;
      if (!fIncludeSaturatedForRiseTimeFit && currentStation.IsLowGainSaturation())
        rejectCode |= sevt::StationRecData::eLowGainSaturated;
      if (stationDistance < fMinimumDistanceForRiseTimeFit || stationDistance > fMaximumDistanceForRiseTimeFit)
        rejectCode |= sevt::StationRecData::eNotInRing;

      if (rejectCode)
        ++fRejectedStations;

      const auto currentStationRisetimeData = new StationRisetimeDataLLL;
      currentStationRisetimeData->fStationId = currentStation.GetId();
      currentStationRisetimeData->fCorrRisetime = stationRisetime;
      currentStationRisetimeData->fCorrRisetimeError = stationRisetimeError;
      currentStationRisetimeData->fDistance = stationDistance;
      currentStationRisetimeData->fDistanceError = stationDistanceError;
      currentStationRisetimeData->fRejectCode = rejectCode;

#ifdef DEBUG_Risetime1000LLL
      ostringstream info;
      info << "Station Id: " << currentStationRisetimeData->fStationId << "\n"
              "Risetime: " << currentStationRisetimeData->fCorrRisetime << "\n"
              "Risetime Error: " << currentStationRisetimeData->fCorrRisetimeError << "\n"
              "Distance: " << currentStationRisetimeData->fDistance << "\n"
              "DistanceError: " << currentStationRisetimeData->fDistanceError << "\n"
              "RejectCode: " << currentStationRisetimeData->fRejectCode;
      INFO(info);
#endif

      fgRisetimeResults->fStationData.push_back(currentStationRisetimeData);
      sRecData.SetCorrectedRiseTime(currentStationRisetimeData->fCorrRisetime,
                                    currentStationRisetimeData->fCorrRisetimeError);
      sRecData.SetRiseTimeRejectionCode(currentStationRisetimeData->fRejectCode);
    }
  } else {
    info << "There is no surface reconstruction!\n"
            "Need a core position and zenith angle before this module.";
    INFO(info);
    return eSuccess;
  }

  if (fDoRiseTimeFit) {
    FitEventRiseTime();

    ostringstream info;
    info << "The RiseTime at " << fDistanceToInterpolateFitResult/km
         << " km is " << fgRisetime
         << " +- " << fgRisetimeError
         << " ns";
    INFO(info);

    if (event.HasRecShower() && event.GetRecShower().HasSRecShower()) {
      auto& sRecData = event.GetRecShower().GetSRecShower();
      if (!sRecData.HasRiseTime1000())
        sRecData.MakeRiseTime1000();

       auto& risetime1000Data = sRecData.GetRiseTime1000();
       risetime1000Data.SetRiseTime1000(fgRisetime, fgRisetimeError);
       risetime1000Data.SetChi2Ndof(fgRisetimeResults->fRisetime1000Chi2,
                                    fgRisetimeResults->fRisetime1000NDF);
       risetime1000Data.SetAlphaBeta(fgRisetimeResults->fFitPar0,
                                     fgRisetimeResults->fFitPar1);
       risetime1000Data.SetXmax(fgRisetimeResults->fXmax,
                                fgRisetimeResults->fXmaxErrorUp,
                                fgRisetimeResults->fXmaxErrorDown);
    }

  }

  return eSuccess;
}


double
Risetime1000LLL::FitEventRiseTime()
{
  const unsigned int pointsInGraph = fgRisetimeResults->fStationData.size() - fRejectedStations;

  // If only 2 stations pass the cuts then we will ignore the event
  // M. Healy (18 Sep 2006)
  if (pointsInGraph <= 2) {
    INFO("Event risetime uncalculated due to an insuffiecent number of acceptable stations.");
    return -1;
  }

  vector<double> distances(pointsInGraph);
  vector<double> distancesError(pointsInGraph);
  vector<double> risetimes(pointsInGraph);
  vector<double> risetimesError(pointsInGraph);
  unsigned int point = 0;

  for (const auto s : fgRisetimeResults->fStationData) {
    const auto& currentStation = *s;
    if (currentStation.fRejectCode)
      continue;

    distances[point] = currentStation.fDistance / km;
    distancesError[point] = currentStation.fDistanceError / km;
    risetimes[point] = currentStation.fCorrRisetime / ns;
    risetimesError[point] = currentStation.fCorrRisetimeError / ns;

    ostringstream info;
    info << "Passed Cut: station " << currentStation.fStationId << ", "
            "distance " << distances[point] << " km, "
            "risetime " << risetimes[point] << " ns";
    INFO(info);

    ++point;
  }

  ostringstream info;
  info << pointsInGraph << " stations passed the cuts.";
  INFO(info);

  TGraphErrors riseTimeGraph(pointsInGraph, &distances.front(), &risetimes.front(), nullptr, &risetimesError.front());
  riseTimeGraph.SetMarkerStyle(20);
  riseTimeGraph.SetMarkerColor(2);
  riseTimeGraph.SetLineColor(2);
  riseTimeGraph.SetLineWidth(2);
  TF1 risetimeFit("RisetimeFit", "40+[0]*x+[1]*x*x",
                  fMinimumDistanceForRiseTimeFit/km, fMaximumDistanceForRiseTimeFit/km);
  risetimeFit.SetParLimits(0, 0, 10000);
  risetimeFit.SetParLimits(1, 0, 10000);

  // Offline Modification in uncertanty estimation -2013 by C. Jarne (Before not correlation term included)
  // Added to calculate the correct uncertainties considering the correlation
  double risetime = -1;
  double risetimeErr = -1;
  double risetimeChi2 = -1;
  double risetimeNDF = -1;
  const int failed = riseTimeGraph.Fit(&risetimeFit, "Q", "", fMinimumDistanceForRiseTimeFit/km, fMaximumDistanceForRiseTimeFit/km);
  if (!failed) {
    risetime = risetimeFit.Eval(fDistanceToInterpolateFitResult/km);
    const auto fitter = TVirtualFitter::GetFitter();
    const int nPar = fitter->GetNumberTotalParameters();
    const TMatrixD covar(nPar, nPar, fitter->GetCovarianceMatrix());
    //covar.Print();
    const double rtcov = covar[0][1];
    const double rtcov2 = covar[1][1];
    const double rtcov1 = covar[0][0];

    //risetimeerr = sqrt(pow(risetimeFit.GetParError(0), 2) + pow(risetimeFit.GetParError(1), 2));
    // New estimation including the correlation term 2014/01 C. Jarne
    risetimeErr = sqrt(rtcov2 + rtcov1 + 2*rtcov);

    risetimeChi2 = risetimeFit.GetChisquare();
    risetimeNDF = risetimeFit.GetNDF();
  }

  fgRisetime = risetime*ns;
  fgRisetimeError = risetimeErr*ns;
  fgRisetimeReducedChi2 = risetimeChi2 / risetimeNDF;

  fgRisetimeResults->fFitPar0 = risetimeFit.GetParameter(0) * m/km;
  fgRisetimeResults->fFitPar1 = risetimeFit.GetParameter(1) * m/km * m/km;
  fgRisetimeResults->fRisetime1000 = risetime*ns;
  fgRisetimeResults->fRisetime1000Error = risetimeErr*ns;
  fgRisetimeResults->fRisetime1000Chi2 = risetimeChi2;
  fgRisetimeResults->fRisetime1000NDF = risetimeNDF;

  // This functionality not yet included
  fgRisetimeResults->fXmax = -1;
  fgRisetimeResults->fXmaxErrorUp = -1;
  fgRisetimeResults->fXmaxErrorDown = -1;

  return risetime*ns;
}


double
Risetime1000LLL::RecalculateRiseTime(const sevt::Station& station)
{
  const auto& recData = station.GetRecData();
  int startBin = recData.GetSignalStartSlot() - 4;
  const int stopBin = recData.GetSignalEndSlot();
  if (startBin >= stopBin || startBin < 0 || startBin > 5000)
    startBin = 0;

  int numPMTs = 0;
  double totalRisetime = 0;
  for (const auto& currentPMT : station.PMTsRange()) {
    if (!currentPMT.HasRecData() || !currentPMT.GetRecData().HasVEMTrace(fComponent))
      continue;

    ++numPMTs;
    const auto& pmtRec = currentPMT.GetRecData();
    const auto& vemTrace = pmtRec.GetVEMTrace(fComponent);
    const double risetimeStart =
      TraceAlgorithm::TimeAtRelativeSignalX(vemTrace, startBin, stopBin, 100*fRiseTimeStartFraction);
    const double risetimeStop =
      TraceAlgorithm::TimeAtRelativeSignalX(vemTrace, startBin, stopBin, 100*fRiseTimeStopFraction);

    const double diff = risetimeStop - risetimeStart;
    totalRisetime += diff;
  }

  if (!numPMTs) {
    ERROR("Valid station, but no PMTRecData.  Risetime recalculation failed!");
    return -1;
  } else
    return totalRisetime / numPMTs;
}