FdDoubleBumpFinder.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <evt/Event.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/ShowerFRecData.h>
#include <evt/ShowerRecData.h>
#include <evt/VGaisserHillasParameter.h>
#include <evt/GaisserHillas4Parameter.h>
#include <evt/MultipleGaisserHillasParameters.h>


#include <fdet/Camera.h>
#include <fdet/Eye.h>
#include <fdet/FDetector.h>
#include <fdet/Pixel.h>
#include <fdet/Telescope.h>

#include <fevt/Eye.h>
#include <fevt/EyeHeader.h>
#include <fevt/EyeRecData.h>
#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Pixel.h>
#include <fevt/PixelRecData.h>
#include <fevt/Telescope.h>
#include <fevt/TelescopeRecData.h>

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/VModule.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

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

#include <utl/AugerUnits.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/RandomEngine.h>
#include <utl/Reader.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/TimeInterval.h>
#include <utl/TimeStamp.h>
#include <utl/Transformation.h>

#include "FdDoubleBumpFinder.h"

#include "ProfileFitter.h"

#include <sstream>
#include <iostream>
#include <vector>

//do not use namespace evt, det
using namespace std;
using namespace fwk;
using namespace utl;
using namespace atm;
using namespace fevt;
using namespace sevt;
using namespace evt::gh;
using namespace fdDoubleBumpFinder;


FdDoubleBumpFinder::FdDoubleBumpFinder() : fProfileFitter()
{ }


VModule::ResultFlag
FdDoubleBumpFinder::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("FdDoubleBumpFinder");

  if (!topB) {
    ERROR("Could not find branch FdDoubleBumpFinder");
    return eFailure;
  }

  Branch preSelectEventsB = topB.GetChild("preSelectEvents");
  if (!preSelectEventsB) {
    ERROR("Could not find branch preSelectEvents");
    return eFailure;
  }
  preSelectEventsB.GetData(fPreSelectEvents);

  Branch maxZenithB = topB.GetChild("maxZenith");
  if (!maxZenithB) {
    ERROR("Could not find branch maxZenith");
    return eFailure;
  }
  maxZenithB.GetData(fMaxZenith);

  Branch minAxisPixelsB = topB.GetChild("minAxisPixels");
  if (!minAxisPixelsB) {
    ERROR("Could not find branch minAxisPixels");
    return eFailure;
  }
  minAxisPixelsB.GetData(fMinAxisPixels);

  Branch maxCoreTankDistB = topB.GetChild("maxCoreTankDist");
  if (!maxCoreTankDistB) {
    ERROR("Could not find branch maxCoreTankDist");
    return eFailure;
  }
  maxCoreTankDistB.GetData(fMaxCoreTankDist);

  Branch cutValueB =  topB.GetChild("cutValue");
  if (!cutValueB) {
    ERROR("Could not find branch cutValue");
    return eFailure;
  }
  cutValueB.GetData(fCutValue);

  Branch constraintXFirstB = topB.GetChild("constraintXFirst");
  if (!constraintXFirstB) {
    ERROR("Could not find branch constraintXFirst");
    return eFailure;
  }
  constraintXFirstB.GetData (fConstraintXFirst);

  Branch constraintLambdaB = topB.GetChild("constraintLambda");
  if (!constraintLambdaB) {
    ERROR("Could not find branch constraintLambda");
    return eFailure;
  }
  constraintLambdaB.GetData(fConstraintLambda);

  Branch constraintXFirstMinusXMaxB = topB.GetChild("constraintXFirstMinusXMax");
  if (!constraintXFirstMinusXMaxB) {
    ERROR("Could not find branch constraintXFirstMinusXMax");
    return eFailure;
  }
  constraintXFirstMinusXMaxB.GetData(fConstraintXFirstMinusXMax);

  Branch verbosityB = topB.GetChild("verbosity");
  if (!verbosityB) {
    ERROR("Could not find branch verbosity");
    return eFailure;
  }
  verbosityB.GetData(fVerbosity);

  const string version = GetVersionInfo(VModule::eRevisionNumber);

  std::ostringstream info;
  info << " Version: " << version
       << "\n\t Parameters:"
       << endl
       << "\n\t Preselection of events = " << fPreSelectEvents
       << endl
       << "\n\t Cut value              = " << fCutValue
       << "\n\t Max core tank dis      = " << fMaxCoreTankDist / m                      << " m"
       << "\n\t Min axis pixels        = " << fMinAxisPixels
       << "\n\t Max zenith angle       = " << fMaxZenith / degree                       << "°"
       << endl
       << "\n\t Constraint XFirst      = " << fConstraintXFirst / g * cm * cm           << " g/cm²"
       << "\n\t Constraint Lambda      = " << fConstraintLambda / g * cm * cm           << " g/cm²"
       << "\n\t Constraint XFirst-XMax = " << fConstraintXFirstMinusXMax / g * cm * cm  << " g/cm²"
       << endl;

  if (fVerbosity >- 1) {
    INFO(info);
  }

  return eSuccess;
}


VModule::ResultFlag
FdDoubleBumpFinder::Run(evt::Event& event)
{
  if (fVerbosity > -1)
    INFO("---Starting FdDoubleBumpFinder");

  if (!event.HasFEvent()) {
    if (fVerbosity > -1)
      WARNING("---No FD Event");
    return eContinueLoop;
  }

  FEvent& fdEvent = event.GetFEvent();
  for (FEvent::EyeIterator eyeIter = fdEvent.EyesBegin(ComponentSelector::eHasData);
       eyeIter != fdEvent.EyesEnd(ComponentSelector::eHasData); ++eyeIter) {

    if (fVerbosity > -1) {
      ostringstream ostr;
      ostr << "---Processing eye " << eyeIter->GetId();
      INFO(ostr);
    }

    if (fPreSelectEvents) {
      if (!IsPreselected(event, *eyeIter)) {
        if (fVerbosity > -1)
          INFO("---Not Preselected");
        continue;
      }
    }

    if (!Fit(*eyeIter)) {
      if (fVerbosity > -1)
        INFO("---Not Fitted");
      continue;
    }

    if (!FillRecData(event, *eyeIter)) {
      if (fVerbosity > -1)
        WARNING("!!!Error writing to event");
      continue;
    }
  }

  return eSuccess;
}


bool
FdDoubleBumpFinder::Scan(const vector<double>& depth,
                         const vector<double>& dEdX, const vector<double>& dEdXError,
                         double& xMax1, double& xMax2, double& dEdXMax1, double& dEdXMax2)
{
  double trackMax, trackMin;
  if (depth.size() > 1) {
    trackMax = depth.back();
    trackMin = depth.front();
  } else {
    if (fVerbosity > 0)
      WARNING("Depth.size() < 2; aborting scan");
    return false;
  }

  vector<double>::const_iterator ite = dEdX.begin();
  double energyMax = *ite;
  while (ite < dEdX.end()) {
    if (*ite > energyMax)
      energyMax = *ite;
    ite++;
  }

  const double energyMin = energyMax/10.;

  double chi2Min = 1e99;

  vector<double>::const_iterator itd;
  vector<double>::const_iterator itee;
  double xfirst1         = fConstraintXFirst / g * cm * cm;
  double lambda          = fConstraintLambda / g * cm * cm;
  double xfirstminusxmax = fConstraintXFirstMinusXMax / g * cm * cm;

  for (double ixMax1 = trackMin; ixMax1 <= trackMax; ixMax1 += trackMax/10) {
    for (double ixMax2 = ixMax1; ixMax2 <= trackMax; ixMax2 += trackMax/10) {
      for (double idEdXMax1 = energyMin; idEdXMax1 < energyMax; idEdXMax1 += energyMax/10) {
        for (double idEdXMax2 = energyMin; idEdXMax2 < energyMax; idEdXMax2 += energyMax/10) {
          double chi2=0;
          itd  = depth.begin();
          ite  = dEdX.begin();
          itee = dEdXError.begin();

          while (itd < depth.end()) {
            if (*itee != 0)
              chi2 += pow(profileFit::DoubleGHFunction(*itd, xfirst1, lambda, ixMax1, idEdXMax1, ixMax2 + xfirstminusxmax, lambda, ixMax2, idEdXMax2) - *ite,2)/ *itee / *itee;
            ++itd;
            ++ite;
            ++itee;
          }
          if (chi2 < chi2Min) {
            chi2Min  = chi2;
            xMax1    = ixMax1;
            xMax2    = ixMax2;
            dEdXMax1 = idEdXMax1;
            dEdXMax2 = idEdXMax2;
          }
        }
      }
    }
  }

  if (chi2Min == 1e99) {
    if (fVerbosity > 0)
      WARNING("No Min Found; aborting scan");
    return false;
  }

  if (fVerbosity > 0) {
    ostringstream ostr;
    ostr << "---Scan(): ";
    ostr << setprecision(2) << "Chi² = " << chi2Min
         << "\t" << "XMax1 = " << xMax1
         << "\t" << "XMax2 = " << xMax2
         << "\t" << "dEdXMax1 = " << dEdXMax1
         << "\t" << "dEdXMax2 = " << dEdXMax2;
    INFO(ostr);
  }

  return true;
}


bool
FdDoubleBumpFinder::Fit(fevt::Eye& eye)
{
  // Check for eye data
  if (eye.HasRecData()) {
    const EyeRecData& recData = eye.GetRecData();

    // Check for shower reconstruction
    if (recData.HasFRecShower()) {
      const evt::ShowerFRecData& showerData = recData.GetFRecShower();

      // Check for dE/dX profile
      if (showerData.HasEnergyDeposit()) {
        const TabulatedFunctionErrors& profile = showerData.GetEnergyDeposit();

        if (profile.GetNPoints() < 2) {
          if (fVerbosity > 0)
            WARNING("Need >= 2 points for profile fit; aborting fit");
          return false;
        }

        vector<double> depth;
        transform(profile.XBegin(), profile.XEnd(), back_inserter(depth),
                  [](const auto& elem) { return elem / (g/cm2); });

        vector<double> dEdX;
        transform(profile.YBegin(), profile.YEnd(), back_inserter(dEdX),
                  [](const auto& elem) { return elem / (PeV / (g/cm2)); });

        vector<double> dEdXError;
        transform(profile.YErrBegin(), profile.YErrEnd(), back_inserter(dEdXError),
                  [](const auto& elem) { return elem / (PeV / (g/cm2)); });

        fProfileFitter.SetVerbosity(fVerbosity);
        fProfileFitter.SetThreshold(-1e99);
        profileFit::EFitFunctionType theType = profileFit::eDoubleGaisserHillas;
        fProfileFitter.SetFunctionType(theType);

        vector<double> fitParameters;
        vector<double> fitSteps;
        vector<double> fitLimitsDown;
        vector<double> fitLimitsUp;

        double xMax1start;
        double xMax2start;
        double dEdXMax1start;
        double dEdXMax2start;

        // Search for a local minimum in GH function
        if (Scan(depth, dEdX, dEdXError,
                 xMax1start, xMax2start, dEdXMax1start, dEdXMax2start))
          {
            //Set starting values for double GH fit
            fitParameters.resize(profileFit::eNDGHpars);
            fitSteps.resize(profileFit::eNDGHpars);
            fitLimitsDown.resize(profileFit::eNDGHpars);
            fitLimitsUp.resize(profileFit::eNDGHpars);

            fitParameters[profileFit::eDGHXmax1]   = xMax1start;
            fitSteps[profileFit::eDGHXmax1]        = 10.;
            fitParameters[profileFit::eDGHnMax1]   = dEdXMax1start;
            fitSteps[profileFit::eDGHnMax1]        = dEdXMax1start*0.05;
            fitParameters[profileFit::eDGHLambda1 ]= fConstraintLambda / g * cm * cm;
            fitSteps[profileFit::eDGHLambda1]      = 5.;
            fitParameters[profileFit::eDGHX01]     = fConstraintXFirst / g * cm * cm;
            fitSteps[profileFit::eDGHX01]          = 10;

            fitParameters[profileFit::eDGHXmax2]   = xMax2start;
            fitSteps[profileFit::eDGHXmax2]        = 10;
            fitParameters[profileFit::eDGHnMax2]   = dEdXMax2start;
            fitSteps[profileFit::eDGHnMax2]        = dEdXMax2start*0.05;
            fitParameters[profileFit::eDGHLambda2] = fConstraintLambda / g * cm * cm;
            fitSteps[profileFit::eDGHLambda2]      = 5.;
            fitParameters[profileFit::eDGHX02]     = xMax2start + fConstraintXFirstMinusXMax / g * cm * cm;
            fitSteps[profileFit::eDGHX02]          = 10;

            fitLimitsDown[profileFit::eDGHXmax1]   = 0;
            fitLimitsUp[profileFit::eDGHXmax1]     = 2000;
            fitLimitsDown[profileFit::eDGHnMax1]   = 0;
            fitLimitsUp[profileFit::eDGHnMax1]     = 0;
            fitLimitsDown[profileFit::eDGHLambda1] = 10;
            fitLimitsUp[profileFit::eDGHLambda1]   = 150;
            fitLimitsDown[profileFit::eDGHX01]     = -1000;
            fitLimitsUp[profileFit::eDGHX01]       = 500;

            fitLimitsDown[profileFit::eDGHXmax2]   = 0;
            fitLimitsUp[profileFit::eDGHXmax2]     = 2000;
            fitLimitsDown[profileFit::eDGHnMax2]   = 0;
            fitLimitsUp[profileFit::eDGHnMax2]     = 0;
            fitLimitsDown[profileFit::eDGHLambda2] = 10;
            fitLimitsUp[profileFit::eDGHLambda2]   = 150;
            fitLimitsDown[profileFit::eDGHX02]     = -1000;
            fitLimitsUp[profileFit::eDGHX02]       = 1000;
          }
        fProfileFitter.SetStartParameters(fitParameters, fitSteps, fitLimitsDown, fitLimitsUp);
        fProfileFitter.SetData(depth.size(), &depth[0], &dEdX[0], &dEdXError[0]);

        // Double GH fit
        if (fProfileFitter.Fit() && fProfileFitter.GetStatus() > 0) {
          fDGHChi2 = fProfileFitter.GetChi2();
          fDGHNdf = fProfileFitter.GetNdf();
          const vector<double>& fitPar = fProfileFitter.GetFitParameters();
          fXMax1 = fitPar[profileFit::eDGHXmax1] * (g / cm2);
          fXMax2 = fitPar[profileFit::eDGHXmax2] * (g / cm2);
          fLambda1 = fitPar[profileFit::eDGHLambda1] * (g / cm2);
          fLambda2 = fitPar[profileFit::eDGHLambda2] * (g / cm2);
          fXFirst1 = fitPar[profileFit::eDGHX01] * (g / cm2);
          fXFirst2 = fitPar[profileFit::eDGHX02] * (g / cm2);
          fdEdXMax1 = fitPar[profileFit::eDGHnMax1] * PeV / (g / cm2);
          fdEdXMax2 = fitPar[profileFit::eDGHnMax2] * PeV / (g / cm2);

          if (!showerData.HasGHParameters()) {
            if (fVerbosity > 0)
              WARNING("No GH Paramaters; aborting fit");
            return false;
          }

          const evt::GaisserHillas4Parameter* gh4Pars =
            dynamic_cast<evt::GaisserHillas4Parameter*>(
              &eye.GetRecData().GetFRecShower().GetGHParameters()
              );

          if (!gh4Pars) {
            if (fVerbosity > 0)
              WARNING("Casting of GH Parameters failed; aborting fit");
            return false;
          }

          //Set starting values for single GH fit
          theType = profileFit::eGaisserHillas;
          fProfileFitter.SetFunctionType(theType);
          fitParameters.resize(profileFit::eNGHpars);
          fitSteps.resize(profileFit::eNGHpars);
          fitLimitsDown.resize(profileFit::eNGHpars);
          fitLimitsUp.resize(profileFit::eNGHpars);

          fitParameters[profileFit::eGHXmax]   = gh4Pars->GetXMax() / g * cm2;
          fitSteps[profileFit::eGHXmax]        = 10.;
          fitParameters[profileFit::eGHnMax]   = gh4Pars->GetNMax() / PeV * g / cm2;
          fitSteps[profileFit::eGHnMax]        = 0.05*gh4Pars->GetNMax() / PeV * g / cm2;
          fitParameters[profileFit::eGHLambda] = gh4Pars->GetShapeParameter(eLambda) / g*cm2;
          fitSteps[profileFit::eGHLambda]      = 5.;
          fitParameters[profileFit::eGHX0]     = gh4Pars->GetShapeParameter(eX0) / g * cm2;
          fitSteps[profileFit::eGHX0]          = 10;
          fitLimitsDown[profileFit::eGHXmax]   = 0;
          fitLimitsUp[profileFit::eGHXmax]     = 0;
          fitLimitsDown[profileFit::eGHnMax]   = 0;
          fitLimitsUp[profileFit::eGHnMax]     = 0;
          fitLimitsDown[profileFit::eGHLambda] = 10;
          fitLimitsUp[profileFit::eGHLambda]   = 150;
          fitLimitsDown[profileFit::eGHX0]     = -1000;
          fitLimitsUp[profileFit::eGHX0]       = 500;

          fProfileFitter.SetStartParameters(fitParameters, fitSteps, fitLimitsDown, fitLimitsUp);
          fProfileFitter.SetData(depth.size(), &depth[0], &dEdX[0], &dEdXError[0]);

          //Single GH fit (done with Chi² to be comparable to double GH fit)
          if (fProfileFitter.Fit() && fProfileFitter.GetStatus() > 0) {
            fGHChi2 = fProfileFitter.GetChi2();
            return true;
          }
          else {
            return false;
          }
        }
      }
      else {
        if (fVerbosity > 0)
          WARNING("No dE/dX profile available: aborting fit");
      }
    }
    else {
      if (fVerbosity > 0)
        WARNING("No shower rec data available; aborting fit");
    }
  }
  else {
    if (fVerbosity > 0)
      WARNING("No rec data available; aborting fit");
  }

  return false;
}


bool
FdDoubleBumpFinder::IsPreselected(evt::Event& event, fevt::Eye& eye)
  const
{
  if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
    if (fVerbosity > 0)
      INFO ("No FRecShower; rejecting event");
    return false;
  }

  const fevt::EyeRecData& eyeData = eye.GetRecData();
  const evt::ShowerFRecData& frecShower = eyeData.GetFRecShower();

  // Can be shot by lidar (zenith <= 60°)
  const Point& core = frecShower.GetCorePosition();
  const CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(core);
  const Vector& showerAxis = frecShower.GetAxis();

  if (showerAxis.GetTheta(localCS) > fMaxZenith) {
    if (fVerbosity > 0) {
      ostringstream ostr;
      ostr << "Zenith angle = " << showerAxis.GetTheta(localCS) / degree
           << "°; rejecting event";
      INFO (ostr);
    }
    return false;
  }

  //Minumum Pixel (ADST FDSelection.cc:323)
  if (eyeData.GetTimeFitNDof() - frecShower.GetStationIds().size() + 3 <= fMinAxisPixels) {
    if (fVerbosity > 0) {
      ostringstream ostr;
      ostr << "N Axis Pixel = " << eyeData.GetNTimeFitPixels()
           << "; rejecting event";
      INFO(ostr);
    }
    return false;
  }


  //Has Energy Reconstruction
  if (!frecShower.HasGHParameters()) {
    if (fVerbosity > 0)
      INFO("No GH Parameters; rejecting event");
    return false;
  }


  //No (HighGain)Saturated Pixels (ADST FDSelection.cc:187)
  for (Eye::ConstTelescopeIterator telIter =
       eye.TelescopesBegin(ComponentSelector::eHasData);
       telIter != eye.TelescopesEnd(ComponentSelector::eHasData);
       ++telIter)
  {
    for (Telescope::ConstPixelIterator pixelIter =
         telIter->PixelsBegin(ComponentSelector::eHasData);
         pixelIter != telIter->PixelsEnd(ComponentSelector::eHasData);
         ++pixelIter)
    {
      if (pixelIter->IsHighGainSaturation()) {
        if (fVerbosity > 0)
          INFO("Saturated; rejecting event");
        return false;
      }
    }
  }

  //Distance to axis (from RecDataProvider.cc CheckHybridStations())
  if (!event.HasSEvent()) {
    if (fVerbosity > 0)
      INFO("No SEvent; rejecting event");
    return false;
  }
  const SEvent& sevent = event.GetSEvent();

  det::Detector& Det = det::Detector::GetInstance();
  const sdet::SDetector& theSDetector = Det.GetSDetector();
  CoordinateSystemPtr CS =
    Det.GetFDetector().GetEye(eye.GetId()).GetLocalCoordinateSystem();

  Vector vertical(0,0,1,CS,Vector::kCartesian);
  Vector SDP = eyeData.GetSDP();
  AxialVector horizontal = cross(SDP, vertical);
  horizontal.Normalize();

  double Chi0 = eyeData.GetChiZero();
  double Rp   = eyeData.GetRp();

  Vector eyeCoreVec = Rp / sin(kPi - Chi0) * horizontal;
  eyeCoreVec.TransformTo(CS) ;

  Transformation rot(Transformation::Rotation(-Chi0, SDP, CS));

  Vector axis(rot*horizontal);                // apply rotation
  if (axis.GetTheta(CS) > (kPi/2))
    axis *= -1; // reflect if necessary
  axis.Normalize();

  bool first = true;
  double hottestRho = 0;
  double hottestSignal = 0;
  unsigned int hottestId = 0;
  for (SEvent::ConstStationIterator sditer = sevent.StationsBegin();
       sditer != sevent.StationsEnd(); ++sditer)
  {
    if (!sditer->HasRecData())
      continue;

    const StationRecData& recdata = sditer->GetRecData();

    if (sditer->IsSilent())
      continue;

    if (!sditer->IsCandidate())
      continue;

    unsigned int thisStationId = (unsigned int)sditer->GetId();
    double signal = recdata.GetTotalSignal();

    const sdet::Station& currentSt = theSDetector.GetStation(*sditer);
    const Point& SdPos = currentSt.GetPosition();
    Vector eyeTankVec(SdPos.GetX(CS), SdPos.GetY(CS), SdPos.GetZ(CS), CS,
      Vector::kCartesian);
    Vector coreTankVec = eyeTankVec - eyeCoreVec ;

    double tankProj      = (coreTankVec * axis);
    Vector tankInAxisVec = eyeCoreVec + tankProj * axis;
    Vector tankAxisVec   = coreTankVec - (tankProj * axis);
    double rho           = tankAxisVec.GetMag();  // perp. distance to axis

    const vector<unsigned short int>& hybridStations =
      frecShower.GetStationIds();

    for (unsigned int k = 0; k < hybridStations.size(); k++) {
      if (hybridStations[k] == thisStationId) {
        if (first) {
          first = false;
          hottestRho = rho;
          hottestSignal = signal;
          hottestId = thisStationId;
        }
        else if (signal > hottestSignal) { //find hottest hybrid station
          hottestSignal = signal;
          hottestRho = rho;
          hottestId = thisStationId;
        }
        break; //after station with right ID exit loop
      }
    }
  }
  const double fdStationAxisDistance = hottestRho;

  //things from coreTankDistCut (ADST/Analysis/src/FDSelection.cc:696)
  if (hottestId < 1) { //hybrid stations != SD stations
    if (fVerbosity > 0)
      INFO("No hottest station; rejecting event");
    return false;
  }

  // hasfdevent, hastriggeredpixels, hasreconstructedpixels, hasfdtimefit,
  // hasfdsdpfit
  if (!(eyeData.GetSDPFitNDof() > 0) || !(eyeData.GetTimeFitNDof() > 0)) {
    if (fVerbosity > 0)
      INFO("No time fit pixels; rejecting event");
    return false;
  }

  if (fdStationAxisDistance >= fMaxCoreTankDist) {
    if (fVerbosity > 0) {
      ostringstream ostr;
      ostr << "CoreTankDist: " << fdStationAxisDistance
           << "; rejecting event";
      INFO(ostr);
    }
    return false;
  }

  return true;
}


bool
FdDoubleBumpFinder::IsSelected(fevt::Eye& eye)
  const
{
  if (!eye.HasRecData() || !eye.GetRecData().HasFRecShower()) {
    if (fVerbosity > 0)
      INFO("No FRecShower; rejecting event");
    return false;
  }

  // Chi²-Cut
  if ((fDGHChi2 + 2) - fGHChi2 > fCutValue) {
    if (fVerbosity > 0) {
      ostringstream ostr;
      ostr << "Cut Not Fulfilled: (dghChi2 + 2) - chi2 = " << (fDGHChi2 + 2) - fGHChi2
           << "; rejecting event";
      INFO(ostr);
    }
    return false;
  }

  // Both XMax in FOV
  const TabulatedFunctionErrors& profile =
    eye.GetRecData().GetFRecShower().GetLongitudinalProfile();

  const unsigned int size = profile.GetNPoints();
  if (size <= 0) {
    if (fVerbosity > 0)
      INFO("Size of depth = 0; rejecting event");
    return false;
  }

  const double trackMax = *(--profile.XEnd());
  const double trackMin = *profile.XBegin();

  if (fXMax1 < trackMin || fXMax1 > trackMax ||
      fXMax2 < trackMin || fXMax2 > trackMax) {
    if (fVerbosity > 0) {
      ostringstream ostr;
      ostr << "XMax not in FoV; rejecting event: not between "
           << trackMin / (g/cm2)
           << "<>" << trackMax / (g/cm2);
      INFO(ostr);
    }
    return false;
  }
  return true;
}


bool
FdDoubleBumpFinder::FillRecData(evt::Event& /*event*/, fevt::Eye& eye)
{
  using namespace evt;
  GaisserHillas4Parameter::ShapeParameter x01(eX0, fXFirst1);
  GaisserHillas4Parameter::ShapeParameter x02(eX0, fXFirst2);
  GaisserHillas4Parameter::ShapeParameter lambda1(eLambda, fLambda1);
  GaisserHillas4Parameter::ShapeParameter lambda2(eLambda, fLambda2);
  GaisserHillas4Parameter ghPar1(fdEdXMax1, fXMax1, x01, lambda1);
  GaisserHillas4Parameter ghPar2(fdEdXMax2, fXMax2, x02, lambda2);

  ghPar1.SetChiSquare(fDGHChi2, fDGHNdf);
  ghPar2.SetChiSquare(fDGHChi2, fDGHNdf);

  ShowerFRecData& shower = eye.GetRecData().GetFRecShower();
  MultipleGaisserHillasParameters& mghPars = shower.GetMultipleGHParameters();

  mghPars.Add(ghPar1);
  mghPars.Add(ghPar2);
  mghPars.SetChi2Improvement(fDGHChi2 - fGHChi2 - 2);

  if (mghPars.GetVirtualParameters().size() != 2) {
    if (fVerbosity > 0)
      WARNING("Error in writing to event; rejecting event");
    return false;
  }

  ostringstream ostr;
  ostr << "Data needed for Cut" << endl
       << "chi2improvement= " << mghPars. GetChi2Improvement()
       << " XMax1= " << (*(mghPars.GetVirtualParameters()[0])).GetXMax()/(g/cm2)
       << " XMax2= " << (*(mghPars.GetVirtualParameters()[1])).GetXMax()/(g/cm2)
       << endl;
  INFO(ostr); //always print

//  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//  ++++ if you want event IDs written to file, pls uncomment this +++++
//  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//   if (IsSelected(eye)) {
//     if (fVerbosity > -1)
//       INFO("---Double bump found");
//     ofstream pFile;
//     pFile.open("events.txt",ios::app);
//     if (!pFile.is_open()) {
//       if (fVerbosity > 0)
//         WARNING("Error in writing to file; rejecting event");
//       return false;
//     }
//     if (!event.HasSEvent()) {
//       if (fVerbosity > 0)
//         WARNING("Error in writing to file; rejecting event");
//       return false;
//     }
//     ostringstream id;
//     id << event.GetSEvent().GetHeader().GetId();
//     pFile << id.str() <<"\n";
//     pFile.close();
//   }

  return true;
}


VModule::ResultFlag
FdDoubleBumpFinder::Finish()
{
  if (fVerbosity > -1)
    INFO("---Closing FdDoubleBumpFinder");
  return eSuccess;
}