Rd2dLDFFitter.cc

Go to the documentation of this file.

/* ################################################################################################

      Removed support for this module from framework and adst with revision 34470 (01.11.21)

#################################################################################################*/

#include "Rd2dLDFFitter.h"
#include <utl/ErrorLogger.h>

#include <cmath>

#include <fwk/CentralConfig.h>

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/BasicVector.h>
#include <utl/GeometryUtilities.h>
#include <utl/RadioGeometryUtilities.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/StringCompare.h>

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/RadioSimulation.h>


#include <revt/REvent.h>
#include <revt/Station.h>
#include <revt/Header.h>
#include <revt/StationRecData.h>
#include <revt/StationRRecDataQuantities.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

// include root stuff
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TFitResultPtr.h"
#include "TRandom3.h"

#include <cstddef>
#include <string>
#include <cstdlib>
#include <iostream>

#include <sstream>
#include <fstream>
#include <iomanip>
#include <stdlib.h>
#include <sys/stat.h>
#include <math.h>
#include <algorithm>
#include <vector>

using namespace evt;
using namespace fwk;
using namespace std;
using namespace utl;
using namespace det;


namespace Rd2dLDFFitter {

Rd2dLDFFitter::Rd2dLDFFitter() :
  fUseSilentStations(false)
{ }


VModule::ResultFlag
Rd2dLDFFitter::Init()
{
  INFODebug("Rd2dLDFFitter::Init()");

  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("Rd2dLDFFitter");

  // Read in the configurations of the xml file
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);

  topBranch.GetChild("FixCore").GetData(fFixCore);
  topBranch.GetChild("useSilentStations").GetData(fUseSilentStations);
  topBranch.GetChild("extrasigmacheck").GetData(fextrasigmacheck);
  topBranch.GetChild("randomseed").GetData(fRandomseed);

  topBranch.GetChild("A_scale").GetData(fA_scale);
  topBranch.GetChild("C3").GetData(fC3);
  topBranch.GetChild("C4").GetData(fC4);

  topBranch.GetChild("C1theta_0_10").GetData(fC1theta_0_10);
  topBranch.GetChild("C2theta_0_10").GetData(fC2theta_0_10);

  topBranch.GetChild("C1theta_10_20").GetData(fC1theta_10_20);
  topBranch.GetChild("C2theta_10_20").GetData(fC2theta_10_20);

  topBranch.GetChild("C1theta_20_30").GetData(fC1theta_20_30);
  topBranch.GetChild("C2theta_20_30").GetData(fC2theta_20_30);

  topBranch.GetChild("C1theta_30_40").GetData(fC1theta_30_40);
  topBranch.GetChild("C2theta_30_40").GetData(fC2theta_30_40);

  topBranch.GetChild("C1theta_40_50").GetData(fC1theta_40_50);
  topBranch.GetChild("C2theta_40_50").GetData(fC2theta_40_50);

  topBranch.GetChild("C1theta_50_55").GetData(fC1theta_50_55);
  topBranch.GetChild("C2theta_50_55").GetData(fC2theta_50_55);
  topBranch.GetChild("A_scale_50_55").GetData(fA_scale_50_55);

  topBranch.GetChild("C1theta_55").GetData(fC1theta_55);
  topBranch.GetChild("C2theta_55").GetData(fC2theta_55);
  topBranch.GetChild("A_scale_55").GetData(fA_scale_55);

  return eSuccess;
}


VModule::ResultFlag
Rd2dLDFFitter::Run(evt::Event& event)
{
  INFODebug("Rd2dLDFFitter::Run()");

  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("No radio event found!");
    return eSuccess;
  }

  // initialize data structures
  revt::REvent& rEvent = event.GetREvent();
  if (!event.HasRecShower())
    event.MakeRecShower();
  if (!event.GetRecShower().HasRRecShower())
    event.GetRecShower().MakeRRecShower();
  ShowerRRecData& showerrrec = event.GetRecShower().GetRRecShower();

  TGraph2DErrors lateralDistribution2DvxB;

  //const det::Detector& detector = det::Detector::GetInstance();
  const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();


  // get reference coordinate system of detector (usually PampaAmarilla)
  const CoordinateSystemPtr referenceCS =
    det::Detector::GetInstance().GetReferenceCoordinateSystem();

  const Point core = showerrrec.GetReferenceCorePosition(event);
  utl::CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(core);
  utl::Vector showerAxis = showerrrec.GetReferenceAxis(event);
  const double zenith = showerAxis.GetTheta(localCS);
  //const double azimuth = showerAxis.GetPhi(localCS);
  //double zenith = acos(showerAxis.GetZ(CS) / showerAxis.GetMag());

  const int nant = showerrrec.GetParameter(revt::eNumberOfStationsWithPulseFound);

  if (nant < 3) {
    WARNING("Number of signalstations less than 3 -> doing nothing");
    return eSuccess;
  }

  Vector magneticField = showerrrec.GetMagneticFieldVector();
  utl::RadioGeometryUtilities transformation =
    utl::RadioGeometryUtilities(showerAxis, localCS, magneticField);

  magneticField /= magneticField.GetMag();

  unsigned int numberOfSignalStations = 0;
  unsigned int numberOfSilentStations = 0;

  for (revt::REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
       sIt != rEvent.CandidateStationsEnd(); ++sIt) {
    revt::Station& station = *sIt;
    //StationRecData& sRec = station.GetRecData();
    if (station.IsSaturated())   // skip saturated stations
      continue;
    if (!station.HasSignal())    // skip silent station
      ++numberOfSilentStations;
    else
      ++numberOfSignalStations;
  }
  //For now, satured stations are no included.
  if (numberOfSignalStations < 3) {
    WARNING("Number of signalstations after rejection of saturated stations < 3 -> doing nothing!");
    return eSuccess;
  }
  double ratio_sil2sig = 1;
  if (numberOfSilentStations > numberOfSignalStations) {  // limit influence of silent stations only if we have more silent than signal stations
    ratio_sil2sig = numberOfSilentStations / numberOfSignalStations;
  }
  double index_signal_stations[numberOfSignalStations][5];
  unsigned int i = 0;
  for (revt::REvent::CandidateStationIterator sIt = rEvent.CandidateStationsBegin();
       sIt != rEvent.CandidateStationsEnd(); ++sIt) {
    revt::Station& station = *sIt;
    revt::StationRecData& sRec = station.GetRecData();
    if (station.IsSaturated())   // skip saturated stations
      continue;
    if (!fUseSilentStations && !station.HasSignal())   // skip silent stations
      continue;

    const Point sPos = rDetector.GetStation(station).GetPosition();
    double x_VxB = 0;
    double y_VxB = 0;
    double z_VxB = 0;

    // due to the choice of coordinate system the station position will be evaluated relative to the core position
    transformation.GetVectorInShowerPlaneVxB(x_VxB, y_VxB, z_VxB, sPos);
    sRec.SetParameter(revt::eLDFFitStationPositionVxB, x_VxB);
    sRec.SetParameter(revt::eLDFFitStationPositionVxVxB, y_VxB);
    sRec.SetParameter(revt::eLDFFitStationPositionV, z_VxB);

    const double integral = sRec.GetParameter(revt::eSignalEnergyFluenceMag);
    const double uncertainty = sRec.GetParameterError(revt::eSignalEnergyFluenceMag);
    //const double noiseLevel = sRec.GetParameter(revt::eNoiseEnergyFluenceMag);

    stringstream sstr;
    sstr << "Signal: " << integral << "+/-" << uncertainty << "\n"
            "Station position (reference CS): " << sPos.GetX(referenceCS) / utl::m << ", "
         << sPos.GetY(referenceCS) / utl::m << ", " << sPos.GetZ(referenceCS) / utl::m << "\n"
            "Station Position (localCS): " << sPos.GetX(localCS) / utl::m << ", "
         << sPos.GetY(localCS) / utl::m << ", " << sPos.GetZ(localCS) / utl::m << "\n"
            "core position (reference CS): " << core.GetX(referenceCS) << ", "
         << core.GetY(referenceCS) << ", " << core.GetZ(referenceCS) << "\n"
            "showerAxis (local cs): " << showerAxis.GetX(localCS) << ", "
         << showerAxis.GetY(localCS) << ", " << showerAxis.GetZ(localCS) << "\n"
            "magnetic field (localCS): " << magneticField.GetX(localCS) << ",  "
         << magneticField.GetY(localCS) << ", " << magneticField.GetZ(localCS) << "\n"
            "magnetic field (referenceCS): " << magneticField.GetX(referenceCS) << ",  "
         << magneticField.GetY(referenceCS) << ", " << magneticField.GetZ(referenceCS) << "\n"
            "position in vxB, vxvxB, z: " << x_VxB << ", " << y_VxB << ", " << z_VxB << endl;
    INFODebug(sstr.str());

    // set subthreshold stations to 0 and increase their uncertainty by N_silent/N_signal
    if (!station.HasSignal()) {
      lateralDistribution2DvxB.SetPoint(lateralDistribution2DvxB.GetN(), x_VxB, y_VxB, integral);
      lateralDistribution2DvxB.SetPointError(lateralDistribution2DvxB.GetN() - 1, 0., 0.,
                                              uncertainty * ratio_sil2sig);
    } else {
      lateralDistribution2DvxB.SetPoint(lateralDistribution2DvxB.GetN(), x_VxB, y_VxB, integral);
      lateralDistribution2DvxB.SetPointError(lateralDistribution2DvxB.GetN() - 1, 0., 0.,
                                              uncertainty);
      index_signal_stations[i][0] = lateralDistribution2DvxB.GetN() - 1;
      index_signal_stations[i][1] = x_VxB;
      index_signal_stations[i][2] = y_VxB;
      index_signal_stations[i][3] = integral;
      index_signal_stations[i][4] = uncertainty;
      ++i;
    }
  }

  double chi2 = 1e20;
  double sigma_max = 500;  //Upper limit for the sigma parameter

  double covRoot[9][9];

  TF2 ldf_pre_fit(
    "LDF_pre_fit", "[0]*TMath::Exp(-((x-[1])*(x-[1])/([2]*[2])+(y-[3])*(y-[3])/([2]*[2])))",
    lateralDistribution2DvxB.GetXmin() - 100, lateralDistribution2DvxB.GetXmax() + 100,
    lateralDistribution2DvxB.GetYmin() - 100, lateralDistribution2DvxB.GetYmax() + 100
  );

  TF2 ldf2DFit(
    "LDF2DFit",
    "[0]*TMath::Exp(-((x-([1]+[5]))**2+(y-[2])**2)/[3]**2)-[6]*[0]*TMath::Exp(-((x-([1]+[4]))**2+(y-[2])**2)/TMath::Exp([7]+[8]*[3])**2)",
    lateralDistribution2DvxB.GetXmin() - 100, lateralDistribution2DvxB.GetXmax() + 100,
    lateralDistribution2DvxB.GetYmin() - 100, lateralDistribution2DvxB.GetYmax() + 100
  );

  ldf_pre_fit.SetParameter(0, lateralDistribution2DvxB.GetZmax());
  ldf_pre_fit.SetParLimits(0, 1e-20 * 6.24150934 * 1e18, 1e-14 * 6.24150934 * 1e18);
  ldf_pre_fit.SetParameter(1, 0);
  ldf_pre_fit.SetParameter(2, 100);
  ldf_pre_fit.SetParLimits(2, 0, 1000);
  ldf_pre_fit.SetParameter(3, 0);

  lateralDistribution2DvxB.Fit(&ldf_pre_fit);

  if (gMinuit->fCstatu == "CONVERGED ") {
    ldf2DFit.SetParameter(0, ldf_pre_fit.GetParameter(0));
    ldf2DFit.SetParLimits(0, 1e-20 * 6.24150934 * 1e18, 1e-14 * 6.24150934 * 1e18);
    ldf2DFit.SetParameter(1, ldf_pre_fit.GetParameter(1));
    ldf2DFit.SetParameter(2, ldf_pre_fit.GetParameter(3));
    ldf2DFit.SetParameter(3, ldf_pre_fit.GetParameter(2));
    ldf2DFit.SetParLimits(3, 0, sigma_max);
  } else {
    ldf2DFit.SetParameter(0, lateralDistribution2DvxB.GetZmax());
    ldf2DFit.SetParLimits(1, 1e-20 * 6.24150934 * 1e18, 1e-14 * 6.24150934 * 1e18);
    ldf2DFit.SetParameter(1, 0);
    ldf2DFit.SetParameter(2, 0);
    ldf2DFit.SetParameter(3, 50);
    ldf2DFit.SetParLimits(3, 0, sigma_max);
  }

  double a_scale = fA_scale;
  double c1 = 0;
  double c2 = 0;
  double c0 = 0.41;

  if (zenith < 10*utl::degree) {
    c1 = fC1theta_0_10;
    c2 = fC2theta_0_10;
  } else if (zenith < 20*utl::degree) {
    c1 = fC1theta_10_20;
    c2 = fC2theta_10_20;
  } else if (zenith < 30*utl::degree) {
    c1 = fC1theta_20_30;
    c2 = fC2theta_20_30;
  } else if (zenith < 40*utl::degree) {
    c1 = fC1theta_30_40;
    c2 = fC2theta_30_40;
  } else if (zenith < 50*utl::degree) {
    c1 = fC1theta_40_50;
    c2 = fC2theta_40_50;
  } else if (zenith < 55*utl::degree) {
    c1 = fC1theta_50_55;
    c2 = fC2theta_50_55;
    a_scale = fA_scale_50_55;
    c0 = 0.46;
  } else {
    c1 = fC1theta_55;
    c2 = fC2theta_55;
    a_scale = fA_scale_55;
    c0 = 0.71;
  }

  ldf2DFit.FixParameter(4, c2);
  ldf2DFit.FixParameter(5, c1);
  ldf2DFit.FixParameter(6, a_scale);
  ldf2DFit.FixParameter(7, fC3);
  ldf2DFit.FixParameter(8, fC4);

  const bool fitCore = !(numberOfSignalStations < 5 || fFixCore == true);
  if (!fitCore) {
    ldf2DFit.FixParameter(1, 0.);
    ldf2DFit.FixParameter(2, 0.);
  }

  TFitResultPtr result = lateralDistribution2DvxB.Fit(&ldf2DFit, "S");
  chi2 = ldf2DFit.GetChisquare();
  INFODebug("chi2 with prefit input " + to_string(chi2));
  bool fitSucessful = false;
  if (gMinuit->fCstatu == "CONVERGED " || gMinuit->fCstatu == "OK " || gMinuit->fCstatu == "SUCCESSFUL") {
    fitSucessful = true;
    gMinuit->mnemat(&covRoot[0][0], 9);
  }

  //vary the initial parameter for the width of the big Gaussian to avoid local minima
  for (int k = 1; k < 9; ++k) {
    ldf2DFit.SetParameter(3, k * 25);
    TFitResultPtr result_temp = lateralDistribution2DvxB.Fit(&ldf2DFit, "S");
    INFODebug("varying sigma= " + to_string(k*25) + "=> chi2 " + to_string(ldf2DFit.GetChisquare()));
    if (ldf2DFit.GetChisquare() < chi2) {
      if (gMinuit->fCstatu == "CONVERGED " || gMinuit->fCstatu == "OK " || gMinuit->fCstatu == "SUCCESSFUL") {
        fitSucessful = true;
        result = result_temp;
        chi2 = ldf2DFit.GetChisquare();
        gMinuit->mnemat(&covRoot[0][0], 9);
      }
    }
  }

  if (!fitSucessful) {
    INFO("None of the fit attempts delivered a successful fit -> SKIP");
    rEvent.SetRejected(revt::REvent::e2DLDFFitfailed);
    return eSuccess;
  }
  if (result->Parameter(3) > sigma_max - 1.) {
    rEvent.SetRejected(revt::REvent::e2DLDFsigmaatlimit);
    INFO("final reconstructed sigma at parameter limit -> SKIP");
    return eSuccess;
  }

  /* vary the energy densities 1000 times within their uncertainty and re-fit the LDF.
  The resulting distribution of the fit parameters can be used to deselect events with unstable fits. */
  if (fextrasigmacheck) {
    TRandom3 r(fRandomseed);
    //const evt::Header& header = event.GetHeader();
    ldf2DFit.SetFitResult(*result);
    vector<double> var_A;
    vector<double> var_sigma;
    for (int k = 0; k < 1000; ++k) {
      for (unsigned int l = 0; l < numberOfSignalStations; ++l) {
        lateralDistribution2DvxB.SetPoint(
            index_signal_stations[l][0], index_signal_stations[l][1], index_signal_stations[l][2],
            r.Gaus(index_signal_stations[l][3], index_signal_stations[l][4]));
      }
      if (fInfoLevel >= eInfoDebug) {
        lateralDistribution2DvxB.Fit(&ldf2DFit, "SV");
      } else {
        lateralDistribution2DvxB.Fit(&ldf2DFit, "SQ");
      }

      if (gMinuit->fCstatu == "CONVERGED " ||
          gMinuit->fCstatu == "OK " ||
          gMinuit->fCstatu == "SUCCESSFUL") {
        if (ldf2DFit.GetParameter(3) < sigma_max - 1.) {
          var_A.push_back(ldf2DFit.GetParameter(0));
          var_sigma.push_back(ldf2DFit.GetParameter(3));
        }
      }
    }

    if (var_A.size() < 500) {
      rEvent.SetRejected(revt::REvent::e2DLDFvariationfails);
      INFO("less than 500 of the 1000 additional fits were successful -> SKIP!");
      return eSuccess;
    }

    //deselect events based on the the distributions of fitted parameters
    const double q_A_15_9 = Quantile(var_A, 15.9);
    const double q_A_50 = Quantile(var_A, 50.);
    const double q_A_84_1 = Quantile(var_A, 84.1);
    INFODebug("amplitude quantiles" + to_string(q_A_15_9) + " " + to_string(q_A_50) +
              " " + to_string(q_A_84_1));

    const double q_sigma_15_9 = Quantile(var_sigma, 15.9);
    const double q_sigma_50 = Quantile(var_sigma, 50.);
    const double q_sigma_84_1 = Quantile(var_sigma, 84.1);
    INFODebug("sigma quantiles" + to_string(q_sigma_15_9) + " " + to_string(q_sigma_50) +
              " " + to_string(q_sigma_84_1));

    bool unstableFit = false;
    if (q_A_50 < (result->Parameter(0) - result->ParError(0)) ||
        q_A_50 > (result->Parameter(0) + result->ParError(0))) {
      rEvent.SetRejected(revt::REvent::e2DLDFmeanAdistoutside);
      INFO("2DLDF Fit rejected, difference between mean of A-distribution and A too big!");
      unstableFit = true;
    }
    if (q_sigma_50 < (result->Parameter(3) - result->ParError(3)) ||
        q_sigma_50 > (result->Parameter(3) + result->ParError(3))) {
      rEvent.SetRejected(revt::REvent::e2DLDFmeanSigmadistoutside);
      INFO("2DLDF Fit rejected, difference between mean of Sigma-distribution and Sigma too big!");
      unstableFit = true;
    }
    if ((q_A_50 - q_A_15_9) / (q_A_84_1 - q_A_50) < 0.5 ||
        (q_A_50 - q_A_15_9) / (q_A_84_1 - q_A_50) > 2) {
      rEvent.SetRejected(revt::REvent::e2DLDFasymmetryAdist);
      INFO("2DLDF Fit rejected, A-distribution too asymmetric!");
      unstableFit = true;
    }
    if ((q_sigma_50 - q_sigma_15_9) / (q_sigma_84_1 - q_sigma_50) < 0.5 ||
        (q_sigma_50 - q_sigma_15_9) / (q_sigma_84_1 - q_sigma_50) > 2) {
      rEvent.SetRejected(revt::REvent::e2DLDFasymmetrySigmadist);
      INFO("2DLDF Fit rejected, Sigma-distribution too asymmetric!");
      unstableFit = true;
    }
    if (unstableFit)
      return eSuccess;
  }

  // fit was successful and stable (if fextrasigmacheck is set)
  // save fit result, increase recStage and compute energy

  showerrrec.Set2dLDFFitResult(*result);
  if (fitCore){
    double x = result->Parameter(1);
    double y = result->Parameter(2);
    utl::Point radioCore = transformation.GetVectorFromShowerPlaneVxB(x, y);

    showerrrec.SetParameter(revt::eCoreX, radioCore.GetX(referenceCS));
    showerrrec.SetParameter(revt::eCoreY, radioCore.GetY(referenceCS));
    showerrrec.SetParameter(revt::eCoreZ, radioCore.GetZ(referenceCS));
  }

  double A = result->Parameter(0);
  double sigma = result->Parameter(3);
  double A_error = result->ParError(0);
  double sigma_error = result->ParError(3);
  double sigmaA_error = result->CovMatrix(0, 3);
  double radiationEnergy = A * TMath::Pi()
      * (pow(sigma, 2) - c0 * TMath::Exp(2 * fC3 + 2 * fC4 * sigma));
  double radiationEnergyError = pow(
      pow(TMath::Pi(), 2) * pow(c0 * TMath::Exp(2 * fC3 + 2 * fC4 * sigma) - pow(sigma, 2), 2)
          * pow(A_error, 2)
          + pow(A, 2) * pow(TMath::Pi(), 2)
              * pow((2 * c0 * fC4 * TMath::Exp(2 * fC4 * sigma + 2 * fC3) - 2 * sigma), 2)
              * pow(sigma_error, 2)
          - TMath::Pi() * (2 * c0 * fC4 * TMath::Exp(2 * fC4 * sigma + 2 * fC3) - 2 * sigma)
              * sigmaA_error,
      0.5);

  showerrrec.SetParameter(revt::eRadiationEnergy, radiationEnergy);
  showerrrec.SetParameterError(revt::eRadiationEnergy, radiationEnergyError);

  INFOIntermediate("Final result: radiation energy = " + to_string(radiationEnergy)
                    + " +/- " + to_string(radiationEnergyError));
  if (fInfoLevel >= eInfoIntermediate) {
    if (event.HasSimShower() && event.GetSimShower().HasRadioSimulation()) {
      const double radiationEnergySim =
          event.GetSimShower().GetRadioSimulation().GetRadiationEnergy();
      if (radiationEnergySim != 0.) {
        const double Erad_resolution = (radiationEnergy - radiationEnergySim) / radiationEnergySim * 100;
        INFOIntermediate(to_string(Erad_resolution) + " % wrt true/simulated radiation energy");
      }
    }
  }

  stringstream sstr;
  sstr << "Final result: sigma = " << sigma << " +/- " << sigma_error << "\n"
       << "Final result: x = " << result->Parameter(1) << " +/- " << result->ParError(1) << "\n"
       << "Final result: y = " << result->Parameter(2) << " +/- " << result->ParError(2) << "\n";
  INFODebug(sstr.str());

  // calculate cosmic ray energy
  const double sineAlpha =
    RadioGeometryUtilities::GetLorentzVector(
      showerrrec.GetReferenceAxis(event),
      showerrrec.GetMagneticFieldVector()
    ).GetMag();
  const double d = 2.3137 * 1e14;  // energy calibration parameters from Aab et al., PRD 93, 122005 (2016)
  const double c = 0.505;
  const double energy = d * std::pow(radiationEnergy / (sineAlpha * sineAlpha), c);
  const double energyError =
    c * d * std::pow(1 / (sineAlpha * sineAlpha), c) *
    std::pow(radiationEnergy, c - 1) * radiationEnergyError;
  showerrrec.SetParameter(revt::eCosmicRayEnergy, energy);
  showerrrec.SetParameterError(revt::eCosmicRayEnergy, energyError);

  double chi2signal = 0.;
  ldf2DFit.SetFitResult(*result);
  for (unsigned int u = 0; u < numberOfSignalStations; ++u) {
    chi2signal += pow(
        (index_signal_stations[u][3] -
         ldf2DFit.Eval(index_signal_stations[u][1], index_signal_stations[u][2]))
          / index_signal_stations[u][4],
        2);
  }

  int NDF = 0;
  if (!fitCore) {
    NDF = numberOfSignalStations - 2;
  } else {
    NDF = numberOfSignalStations - 4;
  }

  showerrrec.SetParameter(revt::e2DLDFfitSignalChi2, chi2signal);
  showerrrec.SetParameter(revt::e2DLDFfitSignalNDF, NDF);

  return eSuccess;
}


VModule::ResultFlag
Rd2dLDFFitter::Finish()
{
  INFODebug("Rd2dLDFFitter::Finish()");
  return eSuccess;
}


double
Rd2dLDFFitter::Quantile(vector<double> data, const double percent)
  const
{
  sort(data.begin(), data.end());
  const unsigned int length = data.size() - 1;
  const int index_above = ceil(length * percent / 100.);
  const int index_below = floor(length * percent / 100.);
  if (index_above == index_below) {
    return data[index_above];
  } else {
    const double m = data[index_above] - data[index_below];
    const double b = data[index_above] - m * index_above;
    return length * percent / 100. * m + b;
  }
}

}