RdGeoCeLDFFitter.cc

Go to the documentation of this file.

#include "RdGeoCeLDFFitter.h"

#include <utl/ErrorLogger.h>

#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 root stuff
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TFitResultPtr.h"
#include "TRandom3.h"

#include <cstddef>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <sys/stat.h>

#include "Minuit2/Minuit2Minimizer.h"
#include "Minuit2/MnUserParameters.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnContours.h"
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MinimumParameters.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnPlot.h"
#include "Minuit2/MnUserCovariance.h"
#include <boost/accumulators/statistics/variance.hpp>

#include "LikelihoodFunction.h"

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


namespace RdGeoCeLDFFitter {

  // gets the most probable Xmax value for a given energy E and mass number A according to EPOS-LHC model
  double
  GetXmaxGumble(const double e, const double a)
  {
    const double lgE = log10(e / utl::eV) - 19;
    const double lnA = log(a);

    const double p0 = 775.589 - 7.047 * lnA - 2.427 * lnA * lnA;
    const double p1 = 57.589 - 0.743 * lnA + 0.214 * lnA * lnA;
    const double p2 = -0.820 - 0.169 * lnA - 0.027 * lnA * lnA;

    const double xmax = p0 + p1 * lgE + p2 * lgE * lgE;
    return xmax;
  }


  VModule::ResultFlag
  RdGeoCeLDFFitter::Init()
  {
    const Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdGeoCeLDFFitter");
    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("minNuberOfStationsForCoreFit").GetData(fMinNuberOfStationsForCoreFit);
    topBranch.GetChild("seedForCoreRandomization").GetData(fRandomSeed);
    topBranch.GetChild("fixCore").GetData(fFixCore);

    return eSuccess;
  }


  VModule::ResultFlag
  RdGeoCeLDFFitter::Run(evt::Event& event)
  {
    // 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();

    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);

    const int Nant = rEvent.GetNumberOfSignalStations();
    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();

    ostringstream info;
    info << "zenith, azimuth (rad) = " << zenith / utl::rad << ", " << azimuth / utl::rad << "\n"
            "zenith, azimuth (deg) = " << zenith / utl::deg << ", " << azimuth / utl::deg;
    INFODebug(info);

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

    for (const auto& station : rEvent.CandidateStationsRange()) {
      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 signal stations after rejection of saturated stations < 3 -> doing nothing!");
      return eSuccess;
    }

    const bool fitCore = !fFixCore && int(numberOfSignalStations) >= fMinNuberOfStationsForCoreFit;
    info.str("");
    info << "Core fit will ";
    if (!fitCore)
      info << "not ";
    info << "be performed. The number of signal stations is " << numberOfSignalStations;
    INFOIntermediate(info);

    /*
    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;
    } */

    const double sineAlpha =
      RadioGeometryUtilities::GetLorentzVector(
        showerrrec.GetReferenceAxis(event), showerrrec.GetMagneticFieldVector()
      ).GetMag();

    // set EventFitData
    EventFitData eventFitData;
    eventFitData.sinalpha = sineAlpha;
    eventFitData.zenith = zenith;
    info.str("");
    info << "setting zenith = " << zenith / utl::deg << "  sinalpha = " << sineAlpha;
    INFODebug(info);

    std::vector<StationFitData> vStationFitData;
    std::vector<StationFitData> vStationFitDataSignalStations;

    // we use the station iterator here to save the position in the vxB frame for all stations
    for (auto& station : rEvent.StationsRange()) {
      revt::StationRecData& sRec = station.GetRecData();

      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);

      if (station.IsSaturated())  // skip saturated stations
        continue;

      if (!fUserSubThresholdStations && !station.HasSignal())  // skip sub threshold stations
        continue;

      if (station.IsRejected())
        continue;

      StationFitData tmpStationFitData;
      tmpStationFitData.f = sRec.GetParameter(revt::eSignalEnergyFluenceMag);
      tmpStationFitData.ferror = sRec.GetParameterError(revt::eSignalEnergyFluenceMag);
      tmpStationFitData.f_vB = sRec.GetParameter(revt::eSignalEnergyFluenceVxB);
      tmpStationFitData.ferror_vB = sRec.GetParameterError(revt::eSignalEnergyFluenceVxB);
      tmpStationFitData.f_vvB = sRec.GetParameter(revt::eSignalEnergyFluenceVxVxB);
      tmpStationFitData.ferror_vvB = sRec.GetParameterError(revt::eSignalEnergyFluenceVxVxB);
      tmpStationFitData.x = X_VxB;
      tmpStationFitData.y = Y_VxB;
      tmpStationFitData.distance = sqrt(X_VxB * X_VxB + Y_VxB * Y_VxB);
      tmpStationFitData.hasSignal = station.HasSignal();

      // the noise level is added to reflect the uncertainty of the model
      //double uncertainty = sqrt(Sqr(integralUncertainty) + Sqr(noiseLevel));

      info.str("");
      info << "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 << endl << "\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);
      INFODebug(info);

      // 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;
      } */

      if (station.HasSignal()) {
        info.str("");
        info << "position in vxB, vxvxB, z: " << X_VxB << ", " << Y_VxB << ", " << Z_VxB << " "
                "f " << tmpStationFitData.f << " +- "
             << sRec.GetParameter(revt::eNoiseEnergyFluenceMag) << " +- "
             << sRec.GetParameterError(revt::eSignalEnergyFluenceMag) << " "
                "fvB " << tmpStationFitData.f_vB << " +- "
             << sRec.GetParameter(revt::eNoiseEnergyFluenceVxB) << " +- "
             << sRec.GetParameterError(revt::eSignalEnergyFluenceVxB) << " "
                "fvvB " << tmpStationFitData.f_vvB << " +- " << tmpStationFitData.ferror_vvB;
        INFODebug(info);
        vStationFitDataSignalStations.push_back(tmpStationFitData);
      }

      vStationFitData.push_back(tmpStationFitData);
    }

    // keep only a few non signal stations

    INFODebug("setting fit parameters");

    FitConfig fitConfig;

    LDFLikelihoodFunction fitFunction = LDFLikelihoodFunction(fitConfig, eventFitData, vStationFitData);

    double Erad_start = 10;
    double Erad_start_error = 10;
    double CRenergy = 1 * utl::EeV;

    // calculate start value for Erad from SD cosmic-ray energy
    if (event.GetRecShower().HasSRecShower()) {
      CRenergy = event.GetRecShower().GetSRecShower().GetEnergy();
      const double CRenergyError = event.GetRecShower().GetSRecShower().GetEnergyError() + 0.18 * CRenergy;
      Erad_start = 15.8 * utl::MeV * Sqr(sineAlpha * CRenergy / utl::EeV);
      Erad_start_error = 15.8 * utl::MeV * Sqr(sineAlpha) * 2 * CRenergy / utl::EeV * CRenergyError / utl::EeV;
    }

    const double XmaxEstimate = 0.5 * (GetXmaxGumble(CRenergy, 1) + GetXmaxGumble(CRenergy, 56));
    const double DXmaxEstimate = getAtmosphere(1564) / cos(zenith) - XmaxEstimate;

    info.str("");
    info << "start values: Erad = " << Erad_start / utl::MeV << " +- " << Erad_start_error / utl::MeV
         << " MeV, Dxmax = "  << DXmaxEstimate;
    INFOIntermediate(info);

    ROOT::Minuit2::MnUserParameters upar;
    upar.Add("Erad", Erad_start / utl::MeV, Erad_start_error / utl::MeV);
    upar.Add("dxmax", DXmaxEstimate, 50);
    upar.Add("coreX", 0., 100);
    upar.Add("coreY", 0., 100);

    upar.Fix("dxmax");
    upar.Fix("coreX");
    upar.Fix("coreY");

    if (fInfoLevel < eInfoDebug){
      // suppress printing the Minuit2 outout
      // Info in <Minuit2>: VariableMetricBuilder: no improvement in line search
      // Info in <Minuit2>: DavidonErrorUpdator: delgam < 0 : first derivatives increasing along search line
      // Info in <Minuit2>: VariableMetricBuilder: matrix not pos.def, gdel > 0
      gErrorIgnoreLevel = 1001;
    }

    ROOT::Minuit2::MnMigrad migrad_pre(fitFunction, upar);
    // minimize calls migrad algorithm. If migrad fails it calls simplex and then migrad again.
    ROOT::Minuit2::FunctionMinimum minimum_pre = migrad_pre();
    info.str("");
    info << "prefit: " << minimum_pre;
    INFOIntermediate(info);
    migrad_pre.Fix("Erad");

    if (fitCore) {
      migrad_pre.Release("coreX");
      migrad_pre.Release("coreY");

      minimum_pre = migrad_pre();
      info.str("");
      info << "prefit: " << minimum_pre;
      INFOIntermediate(info.str());
    }

    migrad_pre.Release("Erad");
    migrad_pre.Release("dxmax");
    migrad_pre.Fix("coreX");
    migrad_pre.Fix("coreY");
    minimum_pre = migrad_pre();
    info.str("");
    info << "Releasing Erad + dxmax, fixing corex, corey. prefit: " << minimum_pre;
    INFOIntermediate(info.str());

    if (fitCore) {
      migrad_pre.Release("coreX");
      migrad_pre.Release("coreY");
    }

    // minimize calls migrad algorithm. If migrad fails it calls simplex and then migrad again.
    ROOT::Minuit2::FunctionMinimum minimum = migrad_pre();
    info.str("");
    info << "minimum: " << minimum;
    INFOIntermediate(info);

    // abort if fit was not successful
    if (!minimum.IsValid()) {
      info.str("");
      info << "Geo ce LDF fit not converged. No result will be saved, returning eSuccess.";
      INFOFinal(info);
      return eSuccess;
    }

    ROOT::Minuit2::MnUserParameterState minUstate_tmp = minimum.UserState();
    ROOT::Minuit2::MnUserCovariance cov = minUstate_tmp.Covariance();

    // propagate uncertainty of SD core
    double eradRMS = 0;
    double dxmaxRMS = 0;
    if (!fitCore) {
      const int number_of_trials = 100;
      std::vector<double> radiation_energies;
      std::vector<double> dxmaxs;
      utl::Vector coreError = showerrrec.GetReferenceCoreError(event);
      Double_t mu[2] = { 0, 0 };
      Double_t sigma[2] = { 0, 0 };
      sigma[0] = coreError.GetX(localCS);
      sigma[1] = coreError.GetY(localCS);
      double rho = showerrrec.GetReferenceCoreErrorCorrelationXY(event);
      CorRand rnd(mu, sigma, rho, fRandomSeed);
      for (int i = 0; i < number_of_trials; ++i) {
        double dX;
        double dY;
        rnd.GetRnd(&dX, &dY);
        utl::Point tmpcore = core + utl::Vector(dX, dY, 0, localCS);
        utl::CoordinateSystemPtr tmpLocalCS = LocalCoordinateSystem::Create(tmpcore);
        utl::RadioGeometryUtilities tmptransformation = utl::RadioGeometryUtilities(showerAxis, tmpLocalCS, MagneticField);

        std::vector<StationFitData> tmpvStationFitData;
        for (const auto& station : rEvent.CandidateStationsRange()) {
          // we use the station iterator here to save the position in the vxB frame for all stations
          const auto& sRec = station.GetRecData();

          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
          tmptransformation.GetVectorInShowerPlaneVxB(X_VxB, Y_VxB, Z_VxB, sPos);

          if (station.IsSaturated())  // skip saturated stations
            continue;

          if (!fUserSubThresholdStations && !station.HasSignal())  // skip sub threshold stations
            continue;

          StationFitData tmpStationFitData;
          tmpStationFitData.f = sRec.GetParameter(revt::eSignalEnergyFluenceMag);
          tmpStationFitData.ferror = sRec.GetParameterError(revt::eSignalEnergyFluenceMag);
          tmpStationFitData.f_vB = sRec.GetParameter(revt::eSignalEnergyFluenceVxB);
          tmpStationFitData.ferror_vB = sRec.GetParameterError(revt::eSignalEnergyFluenceVxB);
          tmpStationFitData.f_vvB = sRec.GetParameter(revt::eSignalEnergyFluenceVxVxB);
          tmpStationFitData.ferror_vvB = sRec.GetParameterError(revt::eSignalEnergyFluenceVxVxB);
          tmpStationFitData.x = X_VxB;
          tmpStationFitData.y = Y_VxB;
          tmpStationFitData.distance = sqrt(X_VxB * X_VxB + Y_VxB * Y_VxB);
          tmpStationFitData.hasSignal = station.HasSignal();

          tmpvStationFitData.push_back(tmpStationFitData);
        }

        FitConfig fitConfig;

        LDFLikelihoodFunction tmpfitFunction = LDFLikelihoodFunction(fitConfig, eventFitData, tmpvStationFitData);
        ROOT::Minuit2::MnUserParameters tmpupar;
        tmpupar.Add("Erad", minUstate_tmp.Value("Erad"), 0.1 * minUstate_tmp.Value("Erad"));
        tmpupar.Add("dxmax", minUstate_tmp.Value("dxmax"), minUstate_tmp.Value("dxmax") * 0.1);
        tmpupar.Add("coreX", 0., 100);
        tmpupar.Add("coreY", 0., 100);

        tmpupar.Fix("coreX");
        tmpupar.Fix("coreY");

        ROOT::Minuit2::MnMigrad tmpmigrad(tmpfitFunction, tmpupar);
        ROOT::Minuit2::FunctionMinimum tmpminimum = tmpmigrad();
        ROOT::Minuit2::MnUserParameterState tmpminUstate = tmpminimum.UserState();
        radiation_energies.push_back(tmpminUstate.Value("Erad"));
        dxmaxs.push_back(tmpminUstate.Value("dxmax"));
      }

      for (unsigned int i = 0; i < radiation_energies.size(); ++i) {
        eradRMS += Sqr(radiation_energies[i] - minUstate_tmp.Value("Erad"));
        dxmaxRMS += Sqr(dxmaxs[i] - minUstate_tmp.Value("dxmax"));
      }

      eradRMS /= radiation_energies.size();
      dxmaxRMS /= dxmaxs.size();
      eradRMS = sqrt(eradRMS);
      dxmaxRMS = sqrt(dxmaxRMS);
    }

    const double radiationEnergy = minUstate_tmp.Value("Erad") * utl::MeV;
    showerrrec.SetParameter(revt::eGeoCeLDFErad, radiationEnergy);
    showerrrec.SetParameter(revt::eGeoCeLDFDxmax, minUstate_tmp.Value("dxmax"));

    if (fInfoLevel >= eInfoIntermediate && event.HasSimShower() && event.GetSimShower().HasRadioSimulation()) {
      const double radiationEnergySim = event.GetSimShower().GetRadioSimulation().GetRadiationEnergy();
      if (radiationEnergySim) {
        info.str("");
        info << " (" << (radiationEnergy - radiationEnergySim) / radiationEnergySim * 100
             << "% wrt true/simulated radiation energy";
        INFO(info);
      }
    }

    // save core position in auger CS
    utl::Point radioCore;
    if (fitCore) {
      double x = minUstate_tmp.Value("coreX");
      double y = minUstate_tmp.Value("coreY");
      showerrrec.SetParameter(revt::eGeoCeLDFCoreX, x);
      showerrrec.SetParameter(revt::eGeoCeLDFCoreY, y);

      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));
    } else {
      radioCore = core;  // if no core is fitted use the reference core
    }

    utl::RadioGeometryUtilities transformationRadioCore =
      utl::RadioGeometryUtilities(showerAxis, LocalCoordinateSystem::Create(radioCore), MagneticField);

    // we use the station iterator here to save the position in the vxB frame for all stations
    for (auto& station : rEvent.StationsRange()) {
      revt::StationRecData& sRec = station.GetRecData();

      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
      transformationRadioCore.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);
    }

    INFODebug("Calculating chi2...");
    // calculate chi2 only from signal stations
    LDFLikelihoodFunction fitFunctionChi2 = LDFLikelihoodFunction(fitConfig, eventFitData, vStationFitDataSignalStations);
    std::vector<double> pars = {
      minUstate_tmp.Value("Erad"), minUstate_tmp.Value("dxmax"),
      minUstate_tmp.Value("coreX"), minUstate_tmp.Value("coreY")
    };

    const double chi2 = fitFunctionChi2.GetChi2FullEnergyFluence(pars);
    int ndf = numberOfSignalStations - 4;
    if (!fitCore)
      ndf = numberOfSignalStations - 2;
    showerrrec.SetParameter(revt::eGeoCeLDFChi2, chi2);
    showerrrec.SetParameter(revt::eGeoCeLDFNDF, ndf);

    info.str("");
    info << "chi2/ndf = " << chi2 << "/" << ndf;
    INFOIntermediate(info);

    // fixme: check if any element is nan
    if (isnan(cov(0, 0))) {
      WARNING("Cov matrix contains nan. Can not estimate uncertainties");
      return eSuccess;
    }

    // save covariances
    const double radiationEnergyError = sqrt(Sqr(cov(0, 0)) +  Sqr(eradRMS)) * pow(utl::MeV, 2);

    showerrrec.SetParameterCovariance(revt::eGeoCeLDFErad, revt::eGeoCeLDFErad, radiationEnergyError);
    showerrrec.SetParameterCovariance(revt::eGeoCeLDFErad, revt::eGeoCeLDFDxmax, cov(0, 1) * utl::MeV);
    const double dxmaxError = sqrt(Sqr(cov(1, 1)) + Sqr(dxmaxRMS));
    showerrrec.SetParameterCovariance(revt::eGeoCeLDFDxmax, revt::eGeoCeLDFDxmax, dxmaxError);

    if (fitCore) {
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFErad, revt::eGeoCeLDFCoreX, cov(0, 2) * utl::MeV);
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFErad, revt::eGeoCeLDFCoreY, cov(0, 3) * utl::MeV);
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFDxmax, revt::eGeoCeLDFCoreX, cov(1, 2));
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFDxmax, revt::eGeoCeLDFCoreY, cov(1, 3));
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFCoreX, revt::eGeoCeLDFCoreX, cov(2, 2));
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFCoreY, revt::eGeoCeLDFCoreY, cov(3, 3));
      showerrrec.SetParameterCovariance(revt::eGeoCeLDFCoreX, revt::eGeoCeLDFCoreY, cov(2, 3));
    }

    return eSuccess;
  }


  VModule::ResultFlag
  RdGeoCeLDFFitter::Finish()
  {
    return eSuccess;
  }

}