SdCompositionParameters.cc

Go to the documentation of this file.

// File "SdCompositionParameters.cc"
//
// This module calculates composition sensitive variables for an event.
// In includes Risetime calculation from the former Module "Risetime1000LLL", which is obsolete now.
// The risetime 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.
//
//
// Author:    Philipp Papenbreer, M.D.Healy, D.Barnhill
// Modified: Feb 2014, P.Papenbreer
// Created:   May 24, 2016
// Modified:  May 24, 2016
//********************************************************************
//
//
//
//History of Risetime1000LLL.cc parts:
// 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.
//
//
// C. Jarne  (2014/02) New asymmetry correction.
// C. Jarne  (2014/02) Modification in  uncertanty estimation according
// to gap2013-079 (now correlation term included).

// 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>

#include <fwk/VModule.h>

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

// Headers for this module
#include "SdCompositionParameters.h"
#include <evt/ShowerSRecDataQuantities.h>

using namespace std;
using namespace fwk;
using namespace utl;
using namespace sevt;
using namespace SdCompParam;
using std::ostringstream;
using std::vector;


// These are static variables and will be refered to in this manner throughout
double SdCompositionParameters::fgRisetime = -1;
double SdCompositionParameters::fgRisetimeError = -1;
double SdCompositionParameters::fgRisetimeReducedChi2 = -1;
SdCompositionParameterResults* SdCompositionParameters::fgCompositionParameterResults = NULL;
DeltaResults* SdCompositionParameters::fgDeltaResults = NULL;
// End static variables

//#define DEBUG_Risetime1000LLL


SdCompositionParameterResults::SdCompositionParameterResults() :
  fFitPar0(0),
  fFitPar1(0),
  fRisetime1000(-1),
  fRisetime1000Error(-1),
  fRisetime1000Chi2(-1),
  fRisetime1000NDF(-1),
  fXmax(-1),
  fXmaxErrorUp(-1),
  fXmaxErrorDown(-1)
{ }


SdCompositionParameterResults::~SdCompositionParameterResults()
{
  for (std::vector<StationSdParameterData*>::iterator sIt = fStationData.begin();
       sIt != fStationData.end(); ++sIt)
    delete *sIt;
}


DeltaResults::DeltaResults()
{
  Delta = -1;
}


DeltaResults::~DeltaResults()
{
  for (std::vector<StationSdParameterData*>::iterator sIt = fStationData.begin();
       sIt != fStationData.end(); ++sIt)
    delete *sIt;
}


SdCompositionParameters::SdCompositionParameters()
{
  // Define the Risetime
  fRiseTimeStartPercent = 0.1;  //10%
  fRiseTimeStopPercent = 0.5;   //50%
  fForceRiseTimeRecalculation = false;  //False

  // Fit the Risetime
  fDoRiseTimeFit = true;                    //True
  fMinimumSignalForRiseTimeFit = 10;        //10.0 VEM
  fMinimumDistanceForRiseTimeFit = 0;       //0.0 meters
  fMaximumDistanceForRiseTimeFit = 1600;    //1600.0 meters
  fIncludeSaturatedForRiseTimeFit = false;  //False
  fDistanceToInterpolateFitResult = 1000;   //1000.0 meters
  fRTWeightingFunction = NULL;
  fRTWeights = NULL;
  fRejectedStations = 0;
  fRTRejectedStations = 0;                  //Stations only rejected for Risetime studies
  fLDFRejectedStations = 0;                  //Stations rejected for LDF studies
  //Calculate Leeds Delta
  fDoCalculateLeedsDelta = true;
  fDoCalculateLDFParameters = true;
}


SdCompositionParameters::~SdCompositionParameters()
{
  delete[] fRTWeightingFunction;
  fRTWeightingFunction = 0;
  delete fgCompositionParameterResults;
  fgCompositionParameterResults = 0;
}


VModule::ResultFlag
SdCompositionParameters::Init()
{
  CentralConfig* const cc = CentralConfig::GetInstance();
  std::ostringstream info;

  info << "Configuring the FADC Parameters module.\n"
          "***************************************\n";

  Branch topBranch = cc->GetTopBranch("SdCompositionParameters");

  //Risetime Stuff
  //Branch RiseTimeBranch = topBranch.GetChild("RecalculateRisetime");
  Branch riseTimeBranch = topBranch.GetFirstChild();
  riseTimeBranch.GetChild("RiseTimeStart").GetData(fRiseTimeStartPercent);
  riseTimeBranch.GetChild("RiseTimeStop").GetData(fRiseTimeStopPercent);
  AttributeMap attributes = riseTimeBranch.GetAttributes();
  if (StringEquivalent(attributes["ReCalc"], "yes")) {
    info << "Going to recalculate risetimes for stations.\n"
            "Pulse RiseTime defined as:\n"
         << fRiseTimeStartPercent*100 << "% to " << fRiseTimeStopPercent*100
         << "% of total signal.\n";
    fForceRiseTimeRecalculation = true;
  } else {
    info << "Using existing station risetimes (from SdCalibrator).\n";
    fForceRiseTimeRecalculation = false;
  }
  //End Risetime Stuff

  //Risetime Fit Stuff
  info << "\nConfiguring the risetime fit routines.\n";
  Branch risetimeFitBranch;
  //if(RisetimeFitBranch = topBranch.GetChild("DoRisetimeFit"))
  if ((risetimeFitBranch = riseTimeBranch.GetNextSibling())) {
    attributes = risetimeFitBranch.GetAttributes();
  }

  if (StringEquivalent(attributes["Fit"], "yes")) {
    info << "Fitting the risetime for the event is enabled.\n";

    fDoRiseTimeFit = true;
    risetimeFitBranch.GetChild("MinimumSignal").GetData(fMinimumSignalForRiseTimeFit);
    risetimeFitBranch.GetChild("MinimumDistance").GetData(fMinimumDistanceForRiseTimeFit);
    risetimeFitBranch.GetChild("MaximumDistance").GetData(fMaximumDistanceForRiseTimeFit);
    risetimeFitBranch.GetChild("RisetimeEvaluatedAt").GetData(fDistanceToInterpolateFitResult);
    info << "\tMinimum Signal " << fMinimumSignalForRiseTimeFit << " VEM\n"
         << "\tMinimum Distance " << fMinimumDistanceForRiseTimeFit/m << " meters\n"
         << "\tMaximum Distance " << fMaximumDistanceForRiseTimeFit/m << " meters\n"
         << "\tRisetime Distance " << fDistanceToInterpolateFitResult/m << " meters\n";

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

    std::string weightingFunction;
    risetimeFitBranch.GetChild("WeightAsFunctionDistance").GetData(weightingFunction);
    fRTWeightingFunction = new char[weightingFunction.size() + 1];
    std::strcpy(fRTWeightingFunction, weightingFunction.c_str());
    info << "\tWeighting Function: " << fRTWeightingFunction << '\n';
  } else
    info << "Fitting the risetime for the event is disabled.\n";

  fRTWeights = new TFormula("RiseTimeWeights", fRTWeightingFunction);
  if (fRTWeights->IsZombie()) {
    ostringstream error;
    error << "Weight formula defined in SdCompositionParameters.xml contains an error!\n"
             "Current formula is: weight=" << fRTWeightingFunction << "\n"
             "This formula does not conform to the ROOT formula class guidelines.";
    ERROR(error);
    //fRTWeights->Delete();
    // this was memory leak: fRTWeights = NULL;
    delete fRTWeights;
    fRTWeights = 0;
    return eFailure;
  }
  //End Risetime Fit Stuff

  topBranch.GetChild("UseDLECorrectedTrace").GetData(fUseDLECorrectedTrace);
  if (fUseDLECorrectedTrace)
    fComponent = sevt::StationConstants::eDirectLightSubtracted;
  else
    fComponent = sevt::StationConstants::eTotal;

  Branch CalculateLeedsDeltaBranch;
  if ((CalculateLeedsDeltaBranch = topBranch.GetChild("LeedsDelta"))) {
    fDoCalculateLeedsDelta = true;
    attributes = CalculateLeedsDeltaBranch.GetAttributes();
    CalculateLeedsDeltaBranch.GetChild("CalculateLeedsDelta").GetData(fDoCalculateLeedsDelta);
    if (fDoCalculateLeedsDelta) {
      std::string currentFunction;
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkFunction").GetData(currentFunction);
      fRisetimeBenchmarkFunction = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeBenchmarkFunction, currentFunction.c_str());
      fRTBenchmark = new TFormula("RiseTimeBenchmark", fRisetimeBenchmarkFunction);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterA").GetData(currentFunction);
      fRisetimeBenchmarkParameterA = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeBenchmarkParameterA, currentFunction.c_str());
      fRTBenchmarkParameterA = new TFormula("RisetimeBenchmarkParA", fRisetimeBenchmarkParameterA);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterA_0").GetData(fRisetimeBenchmarkParameterA_0);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterA_1").GetData(fRisetimeBenchmarkParameterA_1);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterA_2").GetData(fRisetimeBenchmarkParameterA_2);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterB").GetData(currentFunction);
      fRisetimeBenchmarkParameterB = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeBenchmarkParameterB, currentFunction.c_str());
      fRTBenchmarkParameterB = new TFormula("RisetimeBenchmarkParB", fRisetimeBenchmarkParameterB);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterB_0").GetData(fRisetimeBenchmarkParameterB_0);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterB_1").GetData(fRisetimeBenchmarkParameterB_1);
      CalculateLeedsDeltaBranch.GetChild("RisetimeBenchmarkParameterB_2").GetData(fRisetimeBenchmarkParameterB_2);
      CalculateLeedsDeltaBranch.GetChild("RisetimeAzimuthCorrectionFunction").GetData(currentFunction);
      fRisetimeAzimuthCorrectionFunction = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeAzimuthCorrectionFunction, currentFunction.c_str());
      fRTAzimuthCorrections = new TFormula("RiseTimeAzimuthCorrection", fRisetimeAzimuthCorrectionFunction);
      CalculateLeedsDeltaBranch.GetChild("RisetimeDistanceCorrectionFunction").GetData(currentFunction);
      fRisetimeDistanceCorrectionFunction = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeDistanceCorrectionFunction, currentFunction.c_str());
      fRTDistanceCorrections = new TFormula("RiseTimeDistanceCorrection", fRisetimeDistanceCorrectionFunction);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction1").GetData(currentFunction);
      fRisetimeZenithCorrectionFunction1 = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeZenithCorrectionFunction1, currentFunction.c_str());
      fRTZenithCorrections1 = new TFormula("RiseTimeZenithCorrection1", fRisetimeZenithCorrectionFunction1);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction1Parameterl0").GetData(fRisetimeZenithCorrectionFunction1Parameterl0);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction1Parameterl1").GetData(fRisetimeZenithCorrectionFunction1Parameterl1);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction1Parameterl2").GetData(fRisetimeZenithCorrectionFunction1Parameterl2);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction2").GetData(currentFunction);
      fRisetimeZenithCorrectionFunction2 = new char[currentFunction.size() + 1];
      std::strcpy(fRisetimeZenithCorrectionFunction2, currentFunction.c_str());
      fRTZenithCorrections2 = new TFormula("RiseTimeZenithCorrection2", fRisetimeZenithCorrectionFunction2);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction2Parameterm0").GetData(fRisetimeZenithCorrectionFunction2Parameterm0);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction2Parameterm1").GetData(fRisetimeZenithCorrectionFunction2Parameterm1);
      CalculateLeedsDeltaBranch.GetChild("RisetimeZenithCorrectionFunction2Parameterm2").GetData(fRisetimeZenithCorrectionFunction2Parameterm2);
    }
  }

  Branch CalculateLDFParametersBranch;
  if ((CalculateLDFParametersBranch = topBranch.GetChild("LDFParameters"))) {
    fDoCalculateLDFParameters = true;
    attributes = CalculateLDFParametersBranch.GetAttributes();
    CalculateLDFParametersBranch.GetChild("CalculateLDFParameters").GetData(fDoCalculateLDFParameters);
    if (fDoCalculateLDFParameters) {
      CalculateLDFParametersBranch.GetChild("LDFParametersOptimumDistance").GetData(fLDFParametersOptimumDistance);
      CalculateLDFParametersBranch.GetChild("LDFParametersScaleDistance").GetData(fLDFParametersScaleDistance);
      CalculateLDFParametersBranch.GetChild("Cut1000").GetData(fCut1000);
    }
  }

  Branch GeneralParametersBranch;
  if ((GeneralParametersBranch = topBranch.GetChild("GeneralParameters"))) {
    std::string errorFunction;
    GeneralParametersBranch.GetChild("StationRisetimeErrorFunction").GetData(errorFunction);
    fRTErrorFunction = new char[errorFunction.size() + 1];
    std::strcpy(fRTErrorFunction, errorFunction.c_str());
    fRTErrorFunc = new TFormula("RiseTimeErrorFunc", fRTErrorFunction);
    GeneralParametersBranch.GetChild("StationRisetimeErrorParameter_Ja0").GetData(fRTErrorFunctionParameterJa0);
    GeneralParametersBranch.GetChild("StationRisetimeErrorParameter_Ja1").GetData(fRTErrorFunctionParameterJa1);
    GeneralParametersBranch.GetChild("StationRisetimeErrorParameter_Jb0").GetData(fRTErrorFunctionParameterJb0);
    GeneralParametersBranch.GetChild("StationRisetimeErrorParameter_Jb1").GetData(fRTErrorFunctionParameterJb1);
  }
  info <<"\n***************************************\n"
         "FADC Parameters module initialized.\n";
  INFO(info);

  return eSuccess;
}


VModule::ResultFlag
SdCompositionParameters::Run(evt::Event& theEvent)
{
  std::ostringstream info;
  fgRisetime = -1;
  fgRisetimeError = -1;
  fgRisetimeReducedChi2 = -1;
  delete fgCompositionParameterResults;
  fgCompositionParameterResults = nullptr;
  fRejectedStations = 0;
  fRTRejectedStations = 0;
  fLDFRejectedStations = 0;

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

  if (theEvent.HasRecShower()) {
    if (theEvent.GetRecShower().HasSRecShower()) {
      sevt::SEvent& theSEvent = theEvent.GetSEvent();
      fgCompositionParameterResults = new SdCompositionParameterResults;
      // Prefer the first method, but currently not implimented in v2r2 of the offline
      //const double secZenith = 1.0/std::cos(theEvent.GetRecShower().GetZenith()); // is this in radians?

      // Eliminate once RecShower-> GetAzimuth() and GetZenith() and GetCore() work
      // This line is using the current agreed method for Sd reconstuction where
      // GetX(SiteCoordinateSystem) returns the shower u and GetY(SiteCoordinateSystem)
      // returns the shower v.
      const CoordinateSystemPtr CS = det::Detector::GetInstance().GetSiteCoordinateSystem();
      const double this_u = theEvent.GetRecShower().GetSRecShower().GetAxis().GetX(CS);
      const double this_v = theEvent.GetRecShower().GetSRecShower().GetAxis().GetY(CS);
      const double EventThetaRec = std::asin(sqrt(this_u*this_u + this_v*this_v)); // in radians
      secZenith = 1 / std::cos(EventThetaRec);
      recEnergy = theEvent.GetRecShower().GetSRecShower().GetEnergy();
      // End hack for zenith angle

      for (sevt::SEvent::StationIterator sIt = theSEvent.StationsBegin(); sIt != theSEvent.StationsEnd(); ++sIt) {
        sevt::Station& currentStation = *sIt;
        if (theSEvent.HasStation(currentStation.GetId()) /*&& currentStation.HasVEMTrace(fComponent)*/ && currentStation.HasRecData()) {
          unsigned short rejectcode = 0;
          bool usestation = true;
          int saturationflag = -1;

          // Setting the above variables
          sevt::StationRecData& theStationRecData = currentStation.GetRecData();
          double stationrisetime = fForceRiseTimeRecalculation ?
            RecalculateRiseTime(currentStation) :
            fUseDLECorrectedTrace ? theStationRecData.GetRiseTimeCleaned() :
            theStationRecData.GetRiseTime();

          if (stationrisetime <= 0) {
            ostringstream info;
            info << "Station " << currentStation.GetId()
                 << " has no rise time information, calculating...";
            INFO(info);
            stationrisetime = RecalculateRiseTime(currentStation);
            if (stationrisetime <= 0)
              continue;
          }

          const double stationdistance = theStationRecData.GetSPDistance();
          const double stationdistanceerror = theStationRecData.GetSPDistanceError();
          // 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)
          /* Old risetime corrections.. these were replaced by adaptable ones in the xml
          const double alpha = 96.73 + secZenith*(-282.40 + secZenith*(241.80 - 62.61*secZenith));
          const double gamma = -0.0009572 + secZenith*(0.002068 + secZenith*(-0.001362 + 0.0002861*secZenith));
          const double g = alpha + gamma * stationdistance*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 = theStationRecData.GetAzimuthShowerPlane();
          stationrisetime = stationrisetime - g * std::cos(zeta);
          // M.Healy August 8, 2006
          // End additions by M.Healy
          */

          const double l = fRTZenithCorrections1->Eval(fRisetimeZenithCorrectionFunction1Parameterl0,fRisetimeZenithCorrectionFunction1Parameterl1,fRisetimeZenithCorrectionFunction1Parameterl2,secZenith);
          const double m = fRTZenithCorrections2->Eval(fRisetimeZenithCorrectionFunction2Parameterm0,fRisetimeZenithCorrectionFunction2Parameterm1,fRisetimeZenithCorrectionFunction2Parameterm2,secZenith);
          const double g = fRTDistanceCorrections->Eval(l,m,stationdistance);
          const double zeta = theStationRecData.GetAzimuthShowerPlane();
          stationrisetime = fRTAzimuthCorrections->Eval(stationrisetime, g, zeta);

          const double signal =
            fUseDLECorrectedTrace ?  theStationRecData.GetTotalSignalCleaned() : theStationRecData.GetTotalSignal();
          // slight mistake in error propogation here
          //const double stationrisetimeerror =
          //  fRTWeights->Eval(stationdistance/m, secZenith, signal);
          double pararray[7] = {stationdistance, signal, fRTErrorFunctionParameterJa0, fRTErrorFunctionParameterJa1, fRTErrorFunctionParameterJb0, fRTErrorFunctionParameterJb1, secZenith};

          const double stationrisetimeerror =
            fRTErrorFunc->EvalPar(&pararray[0]);//replace old error function by parametrization of Nicole Krohm
          const double photonbenchmark_A =
            fRTBenchmarkParameterA->Eval(fRisetimeBenchmarkParameterA_0, fRisetimeBenchmarkParameterA_1, fRisetimeBenchmarkParameterA_2, secZenith);
          const double photonbenchmark_B =
            fRTBenchmarkParameterB->Eval(fRisetimeBenchmarkParameterB_0, fRisetimeBenchmarkParameterB_1, fRisetimeBenchmarkParameterB_2, secZenith);
          const double photonbenchmark =
            fRTBenchmark->Eval(photonbenchmark_A, photonbenchmark_B, stationdistance);
          if (signal < fMinimumSignalForRiseTimeFit)
            rejectcode |= sevt::StationRecData::eLowSignal;
          if (!currentStation.IsCandidate())
            rejectcode |= sevt::StationRecData::eNotCandidate;
          if (currentStation.IsLowGainSaturation() && !fIncludeSaturatedForRiseTimeFit)
            rejectcode |= sevt::StationRecData::eLowGainSaturated;
          if (stationdistance < fMinimumDistanceForRiseTimeFit || stationdistance > fMaximumDistanceForRiseTimeFit)
            rejectcode |= sevt::StationRecData::eNotInRing;
          saturationflag = SaturationFlag(currentStation);
          bool stflag = RTCandidateFlag(theEvent, rejectcode, stationdistance, usestation, stationrisetime, signal);
          int num_pmts=0;
          for (sevt::Station::PMTIterator pmtIt = currentStation.PMTsBegin(); pmtIt != currentStation.PMTsEnd(); ++pmtIt) {
            sevt::PMT& currentPMT = *pmtIt;
            if (currentPMT.HasRecData() && currentPMT.GetRecData().HasVEMTrace(fComponent)){
              if(currentPMT.GetRecData().GetVEMTrace().GetSize()>0){
                ++num_pmts;
              }
              TraceD& aVEMTrace = currentPMT.GetRecData().GetVEMTrace(fComponent);

              int start_bin = currentStation.GetRecData().GetSignalStartSlot() - 4;
              const int stop_bin = currentStation.GetRecData().GetSignalEndSlot();

              if (start_bin >= stop_bin || start_bin < 0 || start_bin > 5000)
                start_bin = 0;

              const double risetime_start =
                TraceAlgorithm::TimeAtRelativeSignalX(aVEMTrace, start_bin, stop_bin, 100*fRiseTimeStartPercent);
              const double risetime_stop =
                TraceAlgorithm::TimeAtRelativeSignalX(aVEMTrace, start_bin, stop_bin, 100*fRiseTimeStopPercent);

              const double diff = risetime_stop - risetime_start;
              if ((stationrisetime - diff) / stationrisetime >= 0.0001)
                --num_pmts;
            }
          }

          usestation = num_pmts >= 1;

          // Done setting variables

          if (rejectcode)
            ++fRejectedStations;
          if (!rejectcode&&!stflag)
            ++fRTRejectedStations;
          if (rejectcode & sevt::StationRecData::eNotCandidate)
            ++fLDFRejectedStations;
          // Assign to the station risetime data
          StationSdParameterData* const currentStationSdParameterData = new StationSdParameterData;
          currentStationSdParameterData->fStationId = currentStation.GetId();
          currentStationSdParameterData->fCorrRisetime = stationrisetime;
          currentStationSdParameterData->fCorrRisetimeError = stationrisetimeerror;
          currentStationSdParameterData->fDistance = stationdistance;
          currentStationSdParameterData->fDistanceError = stationdistanceerror;
          currentStationSdParameterData->fSecZenith = secZenith;
          currentStationSdParameterData->fPhotonBenchmark = photonbenchmark;
          currentStationSdParameterData->fSignal = signal;
          currentStationSdParameterData->fRejectCode = rejectcode;
          currentStationSdParameterData->fUseStation = usestation;
          currentStationSdParameterData->fSaturationFlag = saturationflag;
          currentStationSdParameterData->fRTCandidateFlag = stflag;
          // End assign data

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

          fgCompositionParameterResults->fStationData.push_back(currentStationSdParameterData);
          theStationRecData.SetCorrectedRiseTime(currentStationSdParameterData->fCorrRisetime,
                                                 currentStationSdParameterData->fCorrRisetimeError);
          theStationRecData.SetRiseTimeRejectionCode(currentStationSdParameterData->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(theEvent);

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

    if (theEvent.HasRecShower()) {
      if (theEvent.GetRecShower().HasSRecShower()) {
        evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
        if (!theSRecData.HasRiseTime1000())
          theSRecData.MakeRiseTime1000();

         evt::RiseTime1000& theRisetime1000Data = theSRecData.GetRiseTime1000();
         theRisetime1000Data.SetRiseTime1000(fgRisetime, fgRisetimeError);
         theRisetime1000Data.SetChi2Ndof(fgCompositionParameterResults->fRisetime1000Chi2,
                                         fgCompositionParameterResults->fRisetime1000NDF);
         theRisetime1000Data.SetAlphaBeta(fgCompositionParameterResults->fFitPar0,
                                          fgCompositionParameterResults->fFitPar1);
         theRisetime1000Data.SetXmax(fgCompositionParameterResults->fXmax,
                                     fgCompositionParameterResults->fXmaxErrorUp,
                                     fgCompositionParameterResults->fXmaxErrorDown);
      }
    }
  }

  if (fDoCalculateLeedsDelta)
    LeedsDeltaCalculation(theEvent);

  if (fDoCalculateLDFParameters)
    LDFParametersCalculation(theEvent);

  PhotonEnergyCalculation(theEvent, recEnergy, secZenith);

  return eSuccess;
}


VModule::ResultFlag
SdCompositionParameters::Finish()
{
  INFO("SdCompositionParameters::Finish()");

  return eSuccess;
}


double
SdCompositionParameters::FitEventRiseTime(evt::Event& theEvent)
{
  evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
  const unsigned int pointsInGraph =
    fgCompositionParameterResults->fStationData.size() - fRejectedStations - fRTRejectedStations;

  // 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);
  vector<bool> stflag(pointsInGraph);
  unsigned int point = 0;
  int nstClose = 0;

  for (std::vector<StationSdParameterData*>::const_iterator stationrisetimeiterator = fgCompositionParameterResults->fStationData.begin();
       stationrisetimeiterator != fgCompositionParameterResults->fStationData.end(); ++stationrisetimeiterator) {
    const StationSdParameterData& currentStation = **stationrisetimeiterator;
    if (currentStation.fRejectCode)
      continue;
    if (!currentStation.fRTCandidateFlag)
      continue;
    distances[point] = currentStation.fDistance/km;
    distanceserror[point] = currentStation.fDistanceError/km;
    risetimes[point] = currentStation.fCorrRisetime/ns;
    risetimeserror[point] = currentStation.fCorrRisetimeError/ns;
    stflag[point] = currentStation.fRTCandidateFlag;

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

    if (currentStation.fDistance<=1000.){
      nstClose++;
    }
    ++point;
  }

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

  TGraphErrors riseTimeGraph(pointsInGraph, &distances.front(), &risetimes.front(), 0, &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;
  bool risetime1000Flag = false;
  bool risetime1000FlagAlt = false;

  const int failed = riseTimeGraph.Fit(&risetimeFit, "Q", "", fMinimumDistanceForRiseTimeFit/km, fMaximumDistanceForRiseTimeFit/km);
  if (!failed) {
    risetime = risetimeFit.Eval(fDistanceToInterpolateFitResult/km);
    TVirtualFitter* const 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();

    risetime1000Flag = SdCompositionParameters::SelectionRiseTime1000(pointsInGraph, nstClose, theEvent.GetRecShower().GetSRecShower().GetEnergy(), false);
    risetime1000FlagAlt = SdCompositionParameters::SelectionRiseTime1000(pointsInGraph, nstClose, theEvent.GetRecShower().GetSRecShower().GetEnergy(), true);
  }

  fgRisetime = risetime*ns;
  fgRisetimeError = risetimeerr*ns;
  fgRisetimeReducedChi2 = risetimechi2/risetimeNDF;

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

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

  theSRecData.SetParameter(eRiseTime1000, risetime*ns);
  theSRecData.SetParameterCovariance(eRiseTime1000, eRiseTime1000, risetimeerr*ns*risetimeerr*ns);
  theSRecData.SetParameter(eRiseTime1000Chi2, risetimechi2);
  theSRecData.SetParameter(eRiseTime1000NDF, risetimeNDF);
  theSRecData.SetParameter(eRiseTime1000NSt, pointsInGraph);
  theSRecData.SetParameter(eRiseTime1000NStClose, nstClose);
  theSRecData.SetParameter(eRiseTime1000Flag, risetime1000Flag);
  theSRecData.SetParameter(eRiseTime1000FlagAlt, risetime1000FlagAlt);
  return risetime*ns;
}


double
SdCompositionParameters::RecalculateRiseTime(sevt::Station& theStation)
{
  int num_pmts = 0;
  double total_risetime = 0;
  std::vector<double> risetime;
  for (sevt::Station::PMTIterator pmtIt = theStation.PMTsBegin(); pmtIt != theStation.PMTsEnd(); ++pmtIt) {
    sevt::PMT& currentPMT = *pmtIt;
    if (currentPMT.HasRecData() && currentPMT.GetRecData().HasVEMTrace(fComponent)) {
      ++num_pmts;
      sevt::PMTRecData& pmtRec = currentPMT.GetRecData();
      TraceD& aVEMTrace = pmtRec.GetVEMTrace(fComponent);

      int start_bin = theStation.GetRecData().GetSignalStartSlot() - 4;
      const int stop_bin = theStation.GetRecData().GetSignalEndSlot();

      if (start_bin >= stop_bin || start_bin < 0 || start_bin > 5000)
        start_bin = 0;

      const double risetime_start =
        TraceAlgorithm::TimeAtRelativeSignalX(aVEMTrace, start_bin, stop_bin, 100*fRiseTimeStartPercent);
      const double risetime_stop =
        TraceAlgorithm::TimeAtRelativeSignalX(aVEMTrace, start_bin, stop_bin, 100*fRiseTimeStopPercent);

      const double diff = risetime_stop - risetime_start;
      risetime.push_back(diff);
      total_risetime += diff;
    }
  }

  if (num_pmts == 0) {
    ERROR("Valid station, but no PMTRecData.  Risetime recalculation failed!");
    return -1;
  }
  const double meanRiseTime = total_risetime / num_pmts;
  /*double rmsRiseTime = 0;
  for (int i = 0; i < num_pmts; ++i)
    rmsRiseTime += pow(meanRiseTime - risetime[i], 2);
  rmsRiseTime = sqrt(rmsRiseTime / num_pmts);*/

  return meanRiseTime;
}


double
SdCompositionParameters::LeedsDeltaCalculation(evt::Event& theEvent)
{
  evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
  const unsigned int pointsInGraph =
    fgCompositionParameterResults->fStationData.size() - fRejectedStations - fRTRejectedStations;

  vector<double> distances(pointsInGraph);
  vector<double> risetimes(pointsInGraph);
  vector<double> risetimeserror(pointsInGraph);
  vector<double> benchmarks(pointsInGraph);
  vector<double> signals(pointsInGraph);
  vector<int> saturationflags(pointsInGraph);
  vector<bool> usestations(pointsInGraph);
  unsigned int point = 0;
  double deltaRiseTime = 0;

  for (std::vector<StationSdParameterData*>::const_iterator stationrisetimeiterator = fgCompositionParameterResults->fStationData.begin();
       stationrisetimeiterator != fgCompositionParameterResults->fStationData.end(); ++stationrisetimeiterator) {
    const StationSdParameterData& currentStation = **stationrisetimeiterator;
    if (currentStation.fRejectCode)
      continue;
    if (!currentStation.fRTCandidateFlag)
      continue;
    distances[point] = currentStation.fDistance;
    risetimes[point] = currentStation.fCorrRisetime/ns;
    risetimeserror[point] = currentStation.fCorrRisetimeError/ns;
    benchmarks[point] = currentStation.fPhotonBenchmark;
    signals[point] = currentStation.fSignal;
    saturationflags[point] = currentStation.fSaturationFlag;
    usestations[point] = currentStation.fUseStation;

    deltaRiseTime += (risetimes[point]-benchmarks[point])/risetimeserror[point];
    ++point;
  }
  deltaRiseTime = deltaRiseTime/double(point);
  bool deltaRiseTimeFlag = SelectionDeltaRiseTime(point, theEvent.GetRecShower().GetSRecShower().GetEnergy(), false);
  bool deltaRiseTimeFlagAlt = SelectionDeltaRiseTime(point, theEvent.GetRecShower().GetSRecShower().GetEnergy(), true);
  theSRecData.SetParameter(eLeedsDelta, deltaRiseTime);
  theSRecData.SetParameter(eDeltaRiseTimeFlag, deltaRiseTimeFlag);
  theSRecData.SetParameter(eDeltaRiseTimeFlagAlt, deltaRiseTimeFlagAlt);

  return deltaRiseTime;
}


double
SdCompositionParameters::LDFParametersCalculation(evt::Event& theEvent)
{
  if (theEvent.HasRecShower()) {
    if (theEvent.GetRecShower().HasSRecShower()) {
      evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
      double R_NKG = 0;
      double S_4 = 0;
      theSRecData.SetShowerSizeType(evt::ShowerSRecData::eS1000);
      double s1000 = theSRecData.GetShowerSize();
      double beta  = theSRecData.GetBeta();
      double gamma = theSRecData.GetGamma();
      const unsigned int pointsInGraph =
        fgCompositionParameterResults->fStationData.size() - fLDFRejectedStations;
      unsigned int point = 0;

      vector<double> signals(pointsInGraph);
      vector<double> distances(pointsInGraph);
      vector<double> saturationflags(pointsInGraph);

      for (std::vector<StationSdParameterData*>::const_iterator stationrisetimeiterator = fgCompositionParameterResults->fStationData.begin();
           stationrisetimeiterator != fgCompositionParameterResults->fStationData.end(); ++stationrisetimeiterator) {
        const StationSdParameterData& currentStation = **stationrisetimeiterator;
        if (currentStation.fRejectCode & sevt::StationRecData::eNotCandidate)
          continue;
        signals[point] = currentStation.fSignal;
        distances[point] = currentStation.fDistance/km;
        saturationflags[point] = currentStation.fSaturationFlag;
        bool sbflag = SbCandidateFlag(saturationflags[point], distances[point], fCut1000);
        if(sbflag){
          double signalexp = s1000 * pow(distances[point]/fLDFParametersOptimumDistance, beta) * pow((distances[point]+fLDFParametersScaleDistance)/(fLDFParametersOptimumDistance+fLDFParametersScaleDistance), (beta+gamma));
          R_NKG += signals[point] / signalexp;
          double signalexp_2 = s1000 * pow(distances[point]/fLDFParametersOptimumDistance, -4);
          S_4 += signals[point] / signalexp_2;
          ++point;
        }
      }
      if (point>0) {
        R_NKG = R_NKG/double(point);
        S_4 = S_4/double(point);
      }
      
      bool rnkg_event_sel = SelectionLDFPar(point, theEvent);
      theSRecData.SetParameter(eLDF_RNKG, R_NKG);
      theSRecData.SetParameter(eLDF_S4, S_4);
      theSRecData.SetParameter(eLDFFlag, rnkg_event_sel);
      theSRecData.SetParameter(eLDF_NSt, point);
    }
  }

  return 0;
}


double
SdCompositionParameters::PhotonEnergyCalculation(evt::Event& theEvent, double E_hadron, double secZenith)
{
  evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
  const double kElongRatePar1WGNew = 872.;
  const double kElongRatePar2WGNew = 64.;
  const double kElongRatePar3WGNew = 37.;
  const double kUniversalityPar1WGCorrGH = 3.03; // 3.04
  const double kUniversalityPar2WGCorrGH = 592.;
  const double kUniversalityPar3WGCorrGH = 152.;
  const double kUniversalityPar4WGCorrGH = 100.;

  TF1 fElongRate("fElongRate", "[0] + [1]*log10(x) + [2]*log10(x)*log10(x)", 0.5, 2.5);
  TF1 fUn("fUn", "[0] * (1.+(x-[3])/[1]) / (1. + pow((x-[3])/[2],2.))", -100., 2000.);
  TF1 fUnGH("fUnGH", "[0] * pow(1.+(x-[2])/[1], [1]/[3]) * exp(-(x-[2])/[3])", -100., 2000.);

  // Pars: [0] VEM/EeV,  [1] g/cm^2, [2] g/cm^2
  //       [3] is the shift between Xmax(axis) and Xmax(R=1000m) g/cm^2
  const double ElongRatePars[3] = {kElongRatePar1WGNew, kElongRatePar2WGNew, kElongRatePar3WGNew};
  const double UnPars[4] = {kUniversalityPar1WGCorrGH, kUniversalityPar2WGCorrGH, kUniversalityPar3WGCorrGH, kUniversalityPar4WGCorrGH};

  double _scale = 1;

  // Set parameters according to EERecOption
  for (int i = 0; i < fElongRate.GetNumberFreeParameters(); ++i)
    fElongRate.SetParameter(i,ElongRatePars[i]);

  for (int i = 0; i < fUnGH.GetNumberFreeParameters(); ++i)
    fUnGH.SetParameter(i, UnPars[i]);

  // Ground slant depth
  static const double xground = 880.0; //is this correct?
  double X = xground * secZenith;
  static const double kF0 = 2; // first guess is Ehadr*kF0
  static const double kERecTolerance = 1e-5;
  static const double kERecMinDeltaX = -50;
  // End hack for zenith angle
  theSRecData.SetShowerSizeType(evt::ShowerSRecData::eS1000);
  double s1000 = theSRecData.GetShowerSize();

  // Iterations start here (use EeV as internal energy unit)
  // initial guess
  double Eold_EeV = kF0 * E_hadron / 1e18;
  double epsilon = 1000;
  int niter = 0;
  double DxTolerance = kERecMinDeltaX * 2;
  double XM = 0;

  while ((X-XM) > DxTolerance && epsilon > kERecTolerance) {
    niter += 1;

    // find Xmax from elongation rate
    XM = fElongRate.Eval(Eold_EeV);

    // find energy from universal profile
    double E_EeV = 0;
    _scale = (0.83 + 2.3e-4 * (X - XM)); // with corr. factor
    E_EeV = pow(s1000/fUnGH.Eval(X-XM), 1./_scale);

    epsilon = std::abs(E_EeV - Eold_EeV);

    Eold_EeV = E_EeV;
  }

  int EgammaFlag = 1;
  // If no iteration or final DX <= -50 g/cm^2 (or 0 g/cm^2)
  if (niter == 0) {
    EgammaFlag = -1;
  } else if ((X - fElongRate.Eval(Eold_EeV)) <= kERecMinDeltaX) {
    EgammaFlag = 0;
  }
  /*
  recEv->SetEgammaRec(Eold_EeV*1e18, opt);
  recEv->SetEgammaIter(niter, opt);
  recEv->SetEgammaXmax(XM, opt);
  */

  theSRecData.SetParameter(eEgammaRec, Eold_EeV * 1e18);
  theSRecData.SetParameter(eEgammaIter, niter);
  theSRecData.SetParameter(eEgammaXmax, XM);
  theSRecData.SetParameter(eEgammaFlag, EgammaFlag);

  return 0;
}


int
SdCompositionParameters::SaturationFlag(sevt::Station& station)
{
  Int_t flag = -2;
  if (station.IsCandidate()) {
    flag = -1;
    const bool saturation = station.IsLowGainSaturation() || station.IsHighGainSaturation();
    if (!saturation) {
      flag = 1;
    } else if (!station.IsLowGainSaturation()) {
      flag = 0;
    }
  }
  return flag;
}


bool
SdCompositionParameters::RTCandidateFlag(evt::Event& theEvent, int saturationflag, double distance, bool usest, double rtcorr, double signal)
{
  evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();
  double beta = theSRecData.GetBeta();
  double gamma = theSRecData.GetGamma();
  theSRecData.SetShowerSizeType(evt::ShowerSRecData::eS1000);
  double s1000 = theSRecData.GetShowerSize();
  double ehadr = theSRecData.GetEnergy();

  double signalexp = s1000 * pow(distance/1000., beta) * pow((distance + 700.) / 1700., beta + gamma);

  double kMinDistBenchNicole = 600;
  double kMaxDistBenchNicole = 2000;
  // use only non low-gain saturated stations
  const bool flag =
    usest &&
    saturationflag >= 0 &&
    ehadr * 1e-18 > 1 &&
    rtcorr > 40 &&
    signal >= 6 && // 3
    signalexp >= 0 && // 10
    kMinDistBenchNicole < distance && distance < kMaxDistBenchNicole;

  return flag;
  //return 0;
}


bool
SdCompositionParameters::SbCandidateFlag(int saturationflag, double distance, bool cut1000)
{
  // Allow only HG-saturation or no saturation at all
  bool flag = (saturationflag >= 0);
  // If selected, use only stations > 1000m
  if (cut1000)
    flag = (flag && distance > 1);

  return flag;
}


bool
SdCompositionParameters::SelectionDeltaRiseTime(int nst, double hadrenergy, bool doAltSel)
{ 
  bool cut = (hadrenergy >= 1e18);
  if (doAltSel)
    cut = (cut && nst >= 4);
  else
    cut = (cut && nst >= 3);
  return cut;
}


bool
SdCompositionParameters::SelectionLDFPar(int nst, evt::Event& theEvent)
{
  bool cut = false;
  if (theEvent.HasRecShower()) {
    if (theEvent.GetRecShower().HasSRecShower()) {
      evt::ShowerSRecData& theSRecData = theEvent.GetRecShower().GetSRecShower();

      cut = (theSRecData.HasLDF() && nst >= 1);
    }
  }
  return cut;
}


bool
SdCompositionParameters::SelectionRiseTime1000(int nst, int nstClose, double hadrenergy, bool doAltSel)
{
  bool cut = hadrenergy * 1e-18 >= 1;
  if (doAltSel)
    cut = (cut && nst >= 3 && nstClose > 0);
  else
    cut = (cut && nst >= 4);

  return cut;
}