// OutNtupleProc.cc
// Contact person: Kevin Labe <klabe@uchicago.edu>
// See OutNtupleProc.hh for more details.
//—————————————————————---------------------------------------------——//

#include <RAT/OutNtupleProc.hh>
#include <RAT/Log.hh>
#include <RAT/Processor.hh>
#include <RAT/DB.hh>
#include <RAT/Version.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/EV.hh>
#include <RAT/MetaInformation.hh>

#include <TFile.h>
#include <TTree.h>
#include <TString.h>

#include <utility>
#include <stdexcept>

namespace RAT {

OutNtupleProc::OutNtupleProc() : Processor("OutNtupleProc")
{
    fDefaultFileName = "";
    fFileName = "";
    fFileNameSuffix = "";
    file = NULL;
    output = NULL;
    branchDS = NULL;
    branchRun = NULL;
    fAutoFlush = 1000000;
    fAutoSave = 10000000;
    branchDS = new DS::Entry();
    McIndex = 0;
}

void OutNtupleProc::BeginOfRun(DS::Run&)
{
// Extract default filename from database, to use if none specified by user
// with /rat/procset file
    DB *db = DB::Get();
    DBLinkPtr io = db->GetLink("IO");
    try{
        fDefaultFileName = io->GetS("default_output_filename");
        if( fDefaultFileName.find( "." ) == std::string::npos )
          fDefaultFileName += ".ntuple.root";
    } catch (DBNotFoundError & e){
        Log::Die("OutNtupleProc::BeginOfRun: Default file name not found"
                 " in IO[] ratdb.");
    }

    // -- Create the output file at this point. No need to wait for the event loop
    // At this point we have everything we need
    if (!file) {
      OpenFile();
    }
}

void OutNtupleProc::EndOfRun(DS::Run& run){
    fMetaHelper.AddRun(run);
}

OutNtupleProc::~OutNtupleProc()
{
    if (file) {
        file->cd();
        if(!output->Write()) Log::Die("Could not write ntuple tree!");

        // write the meta information
        fMetaHelper.SetStoredEvents(output -> GetEntries());
        MetaInformation::Get()->Finalise(fMetaHelper);

        // Write the MetaInformation (RAT and DB state)
        if(!MetaInformation::Get()->GetDSMeta().Write( "meta" ))
            Log::Die("Could not write meta to ROOT file");

        file->Close();
        delete file;
    }
}

void OutNtupleProc::SetS(const std::string& param, const std::string& value)
{
    if (param == "file"){
      fFileName = value;
//        if (!OpenFile(value))
//            Log::Die("OutNtupleProc::SetS: Cannot open file " + value);
    }
    else if ( param == "suffix" ) {
      fFileNameSuffix = value;
    }
    else{
        throw ParamUnknown(param);
    }
}

void OutNtupleProc::OpenFile()
{
  string theFilename = fDefaultFileName;
  if( !fFileName.empty() )
    theFilename = fFileName;
  if (!fFileNameSuffix.empty()) {
    info << "OutNtupleProc: Appending [" << fFileNameSuffix << "] to the output NTuple file name. " << newline;
    string tmp_file = "";
    if( theFilename.find( ".root" ) == string::npos ) {
      // -- there is no .root extension
      tmp_file = theFilename;
    } else {
      size_t pos = theFilename.rfind(".ntuple.root");
      // If there is no 'ntuple.root', look just for '.root'
      if (pos == std::string::npos) {
        pos = theFilename.rfind(".root");
      }
      tmp_file = theFilename.substr(0,pos);
    }
    theFilename = tmp_file;
    theFilename += "_";
    theFilename += fFileNameSuffix;
    theFilename += ".ntuple.root";
  } else {
    debug << "OutNtupleProc: There is **NO** appendix to be added [" << fFileNameSuffix << "] to the file name. " << newline;
  }

    info << "OutNtupleProc: Writing to " << theFilename << newline;
    file = TFile::Open(theFilename.c_str(), "RECREATE");
    if (!file)
      Log::Die("OutNtupleProc::DSEvent: Cannot open output file " + theFilename);

    output = new TTree("output","output");

    // runID uniquely identifies a run.
    // runID comes from GetRunID();
    output->Branch("runID",&mcdata.runID);

    // ratversion tells which version of rat was used.
    // ratversion comes from GetRatVersion();
    output->Branch("ratversion",&mcdata.ratversion);

    // mcIndex uniquely identifies each MC event.
    // It counts the events as they are fed into the processor.
    output->Branch("mcIndex",&McIndex);

    // evIndex uniquely identifies a trigger within a MC event.
    // It counts from 0 to GetEVCount()-1 within each MC event.
    output->Branch("evIndex",&EvIndex);

    // mc tells whether there is MC data associated to a trigger.
    // mc comes from ExistMC();
    output->Branch("mc",&mcdata.mcvalid);

    // eventID uniquely identifies each trigger.
    // eventID comes from GetGTID();
    output->Branch("eventID",&data.eventID);

// The following branches store data on the two most energetic particles in an
// event for which MC data is present.  If there is only one particle, the
// second set of data branches is filled with zeros.  The function
// SearchEnergies() is used to determine which particles are most energetic.

    // pdg1 is the PDG code of the most energetic particle in an event.
    // pdg1 comes from GetPDGCode();
    output->Branch("pdg1",&mcdata.pdg1);
    // pdg2 is the same for the second-most energetic particle.
    output->Branch("pdg2",&mcdata.pdg2);

    // parentpdg1 is the PDG code of the parent of the most energetic particle,
    // if one exists.  It also comes from GetPDGCode();
    output->Branch("parentpdg1",&mcdata.parentpdg1);
    // parentpdg2 is the same, for the second-most energetic particle.
    output->Branch("parentpdg2",&mcdata.parentpdg2);

    // parentKE1 is the kinetic energy of the parent of the most energetic particle,
    // if one exists.
    output->Branch("parentKE1",&mcdata.parentKE1);
    // parentKE2 is the same, for the second-most energetic particle.
    output->Branch("parentKE2",&mcdata.parentKE2);

    // parentMeta1 is the meta info of the parent of the most energetic particle,
    // if one exists.
    output->Branch("parentMeta1",&mcdata.parentMetaInfo1);
    // parentMeta2 is the same, for the second-most energetic particle.
    output->Branch("parentMeta2",&mcdata.parentMetaInfo2);

    // mcke1 is the true kinetic energy of the most energetic particle.
    // mcke1 comes from GetKE();
    output->Branch("mcke1",&mcdata.ke1);
    // mcke2 is the same for the second-most energetic particle.
    output->Branch("mcke2",&mcdata.ke2);

    // mcMomx1 is the true x momentum of the most energetic particle.
    // mcMomx1 comes from GetMomentum().X();
    // mcMomy1 and mcMomz1 are the same for the y and z directions.
    output->Branch("mcmom1x", &mcdata.mom1x);
    output->Branch("mcmom1y", &mcdata.mom1y);
    output->Branch("mcmom1z", &mcdata.mom1z);
    // mcMom*2 is the same for the second-most energetic particle.
    output->Branch("mcmom2x", &mcdata.mom2x);
    output->Branch("mcmom2y", &mcdata.mom2y);
    output->Branch("mcmom2z", &mcdata.mom2z);

    // mctime1 is the true time of the most energetic particle
    output->Branch("mctime1",&mcdata.time1);
    // mctime2 is the same for the second-most energetic particle
    output->Branch("mctime2",&mcdata.time2);
// End of two most energetic particles data.

    // mct is the true time of the event, counting from the start of the run.
    // mct comes from GetMCTime();
    output->Branch("mct",&mcdata.mct);

    // mcEdep is the true total energy deposited in the scintillator.
    // mcEdep comes from GetScintEnergyDeposit();
    output->Branch("mcEdep",&mcdata.Edep);

    // mcEdepQuenched is the true energy deposited with quenched energy removed.
    // mcEdepQuenched comes from GetTotScintEdepQuenched();
    output->Branch("mcEdepQuenched",&mcdata.EdepQuenched);

    // mcNCherPhotons is the number of Cherenkov photons produced
    // mcNCherPhotons comes from GetNCherPhotons();
    output->Branch("mcNCherPhotons",&mcdata.NCherPhotons);

    // mcPosx is the true x position of the event.  (In case of possible
    // ambiguity, it uses the vertex of the most energetic particle).
    // mcPosx comes from GetPosition().X();
    output->Branch("mcPosx",&mcdata.mcposx);
    // mcPosy is the same, for the y direction.
    output->Branch("mcPosy",&mcdata.mcposy);
    // mcPosz is the same, for the z direction.
    output->Branch("mcPosz",&mcdata.mcposz);
    // mcPosr is the radius
    output->Branch("mcPosr",&mcdata.mcposr);

    // nhits is the number of normal and HQE PMTs hit at least once.
    // nhits comes from GetNhits();
    output->Branch("nhits",&data.nhits);

    // owlnhits is the number of OWL PMTs hit at least once.
    // owlnhits comes from GetUncalPMTs().GetOWLCount()
    output->Branch("owlnhits",&data.owlnhits);

    // necknhits is the number of Neck PMTs hit at least once.
    // necknhits comes from GetUncalPMTs().GetNeckCount()
    output->Branch("necknhits",&data.necknhits);

    // nhits cleaned is the number of normal and HQE PMTs hit at least once
    // and pass the channel and hit cleaning cuts
    output->Branch("nhitsCleaned",&data.nhitsCleaned);

    // inTimeHits100 is the number of in-time hits with triggers enabled
    // occuring in a 89 ns time window
    output->Branch("inTimeHits100",&data.inTimeHits100);

    // inTimeHits20 is the number of in-time hits with triggers enabled
    // occuring in a 46 ns time window
    output->Branch("inTimeHits20",&data.inTimeHits20);

    // isCal is whether the event was tagged by a calibration souce
    // returns true if it was a tagged event
    output->Branch("isCal", &data.isCal);

    // fecdHits is a 32-bit integer that contains the FECD hits
    // in the event
    output->Branch("fecdHits", &data.fecdHits);

    // q is the sum of all charge samples in the event.
    // q comes from GetTotalCharge();
    output->Branch("q",&data.Q);

    // trigger is the trigger word.
    // trigger comes from GetTrigType();
    output->Branch("triggerWord",&data.trigger);

    // tubiiWord is the TUBii trigger word.
    // tubiiWord comes from GetTubiiTrig();
    output->Branch("tubiiWord",&data.tubii);

    // dcApplied is the applied data cleaning cuts
    // dcApplied comes from rEV->GetDataCleaningFlags().GetApplied().GetULong64_t(0)
    output->Branch("dcApplied", &data.dcApplied);

    // dcFlagged is the data cleaning word.
    // dcFlagged comes from rEV->GetDataCleaningFlags().GetFlags().GetULong64_t(0)
    output->Branch("dcFlagged",&data.dcFlagged);

    // clockCount10 is the tick count of the trigger on the 10MHz clock.
    // clockCount10 comes from GetClockCount10();
    output->Branch("clockCount10",&data.ClockCount10);

    // clockCount50 is the tick count of the trigger on the 50MHz clock.
    // clockCount50 comes from GetClockCount50();
    output->Branch("clockCount50",&data.ClockCount50);

    // uTDays is part of the time stamp of a trigger, counting number of days
    // since January 1, 2010.
    // uTDays comes from GetUTDays();
    output->Branch("uTDays",&data.UTDays);

    // uTSecs is part of the time stamp of a trigger, counting number of seconds
    // since the beginning of the day.
    // uTSecs comes from GetUTSecs();
    output->Branch("uTSecs",&data.UTSecs);

    // uTNSecs is part of the time stamp of a trigger, giving ns precision.
    // uTNSecs comes from GetUTNSecs();
    output->Branch("uTNSecs",&data.UTNSecs);

    // scintFit tells whether scintFitter results are present in the data.
    // scintFit comes from GetFitResult("scintFitter");
    output->Branch("scintFit",&data.scintFit);

    // waterFit tells whether waterFitter results are present in the data.
    // waterFit comes from GetFitResult("waterFitter");
    output->Branch("waterFit",&data.waterFit);

    // partialWaterFit tells whether partialWaterFitter results are present in the data.
    // partialWaterFit comes from GetFitResult("partialWaterFitter");
    output->Branch("partialWaterFit",&data.partialWaterFit);

    // partialFit tells whether partialFitter results are present in the data.
    // partialFit comes from GetFitResult("partialFitter");
    output->Branch("partialFit",&data.partialFit);

    // Also store the MPF outputs
    output->Branch("posx_MultiPath",&data.posx_MultiPath);
    output->Branch("posy_MultiPath",&data.posy_MultiPath);
    output->Branch("posz_MultiPath",&data.posz_MultiPath);
    output->Branch("posr_MultiPath",&data.posr_MultiPath);
    output->Branch("time_MultiPath",&data.time_MultiPath);

    output->Branch("posz_fastZ",&data.posz_fastZ);
    output->Branch("time_fastZ",&data.time_fastZ);

    // set fitValid bool - is the default fitter valid overall - and fitValidityMask
    // for each component for the fit result.
    output->Branch("fitValid",&data.fitValid);
    output->Branch("fitValidityMask",&data.fitValidityMask);

    // Add validities of additional fitters
    output->Branch("fitValid_MultiPath",&data.fitValid_MultiPath);
    output->Branch("fitValid_fastZ",&data.fitValid_fastZ);

    // posx is the reconstructed position of the event.  It can come from
    // scintFitter, waterFitter, partialWaterFitter or partialFitter.
    // posx comes from GetVertex().GetPosition().X();
    // The fitted position has asymmetric errors, recorded as posxErrorPos and
    // posxErrorNeg.  These come from GetVertex().GetPositivePositionError().X() and
    // GetNegativePositionError().X().
    output->Branch("posx",&data.posx);
    output->Branch("posxPosError",&data.posxPosError);
    output->Branch("posxNegError",&data.posxNegError);
    // posy is the same for the y direction.
    output->Branch("posy",&data.posy);
    output->Branch("posyPosError",&data.posyPosError);
    output->Branch("posyNegError",&data.posyNegError);
    // posz is the same for the z direction.
    output->Branch("posz",&data.posz);
    output->Branch("poszPosError",&data.poszPosError);
    output->Branch("poszNegError",&data.poszNegError);
    // posr is the radius
    output->Branch("posr",&data.posr);


    // posFOM is the figure of merit of the fitted position.
    // posFOM comes from GetFOM("PositionLogL");
    output->Branch("posFOM",&data.posFOM);
    // posFOM2 provides complementary information to the LogL figure of merit above.
    // posFOM2 comes from GetFOM("PositionSelectedNHit");
    output->Branch("posFOM2",&data.posFOM2);

    // dirx is the reconstructed direction of the event (waterFitter only).
    // dirx comes from GetVertex().GetDirection().X();
    // dirxPosError and dirxNegError are the associated error.
    // dirxPosError comes from GetVertex().GetPositiveDirectionError().X();
    output->Branch("dirx",&data.dirx);
    output->Branch("dirxPosError",&data.dirxPosError);
    output->Branch("dirxNegError",&data.dirxNegError);
    // diry is the same for the y direction.
    output->Branch("diry",&data.diry);
    output->Branch("diryPosError",&data.diryPosError);
    output->Branch("diryNegError",&data.diryNegError);
    // dirz is the same for the z direction.
    output->Branch("dirz",&data.dirz);
    output->Branch("dirzPosError",&data.dirzPosError);
    output->Branch("dirzNegError",&data.dirzNegError);

    // dirFOM is the figure of merit of the fitted direction.
    // dirFOM comes from GetFOM("DirectionLogL");
    output->Branch("dirFOM",&data.dirFOM);
    // dirFOM2 provides complementary information to the LogL figure of merit above.
    // dirFOM2 comes from GetFOM("DirectionSelectedNHit");
    output->Branch("dirFOM2",&data.dirFOM2);

    // energy is the reconstructed energy.  It can come from scintFitter,
    // waterFitter, partialWaterFitter or partialFitter. energy comes from GetVertex().GetEnergy();
    // energyPosError and energyNegError are the associated errors.
    // energyPosError comes from GetVertex().GetPositiveEnergyError();
    output->Branch("energy",&data.energy);
    output->Branch("energyPosError",&data.energyPosError);
    output->Branch("energyNegError",&data.energyNegError);

    // energyFOM* are figures of merit for EnergyRSP reconstructed energy
    output->Branch("energyFOMmedProb",&data.energyFOMmedProb);
    output->Branch("energyFOMmedProbHit",&data.energyFOMmedProbHit);
    output->Branch("energyFOMmedDev",&data.energyFOMmedDev);
    output->Branch("energyFOMmedDevHit",&data.energyFOMmedDevHit);
    output->Branch("energyFOMGtest",&data.energyFOMGtest);
    output->Branch("energyFOMUtest",&data.energyFOMUtest);

    // penergy energy is the reconstructed energy from penergy
    output->Branch("penergyEnergy",&data.penergyEnergy);
    output->Branch("penergyError",&data.penergyError);

    output->Branch("penergyAEnergy",&data.penergyAEnergy);
    output->Branch("penergyAError",&data.penergyAError);

    // time is the reconstructed time.  It can come from scintFitter,
    // waterFitter, partialWaterFitter or partialFitter.  time comes from GetVertex().GetTime();
    // timePosError and timeNegError are the associated errors.
    // timePosError comes from GetVertex().GetPositiveTimeError();
    output->Branch("time",&data.time);
    output->Branch("timePosError",&data.timePosError);
    output->Branch("timeNegError",&data.timeNegError);

    // itr is a classifier giving the ratio of hits falling within a time window to the total hits
    // itr comes from GetClassification("ITR");
    output->Branch("itr",&data.itr);

    //XSite is a classifier to distinguish gamma/electron
    //XSite comes from GetClassification("XSite");
    output->Branch("xsite",&data.xsite);

    // timingPeaks is a classifier counting the number of peaks in the PMT hits'
    // timeseries.  timingPeaks comes from GetClassification("timingPeaks");
    output->Branch("timingPeaks",&data.timingPeaks);

    // berkeleyAlphaBeta is a classifier giving the alpha-ness of an event on a scale from
    // 0 to 1.  berkeleyAlphaBeta comes from GetClassification("BerkeleyAlphaBeta");
    output->Branch("berkeleyAlphaBeta",&data.berkeleyAlphaBeta);

    // theta_ij is a classifier describing the anisotropy of an event.
    // theta_ij comes from GetClassification("isotropy");
    output->Branch("thetaij",&data.thetaij);

    // beta14 is a classifier describing the anisotropy of an event.
    // beta14 comes from GetClassification("isotropy");
    output->Branch("beta14",&data.beta14);

    // qpdt is a classifier describing the probability of seeing early charge hits
    // qpdt and nhitsEarly comes from GetClassification("QPDT");
    output->Branch("qpdtProbability",&data.qpdtProbability);
    output->Branch("qpdtNhitsEarly",&data.qpdtNhitsEarly);
    output->Branch("qpdtQMax",&data.qpdtQMax);

    // meanTime is a classifier calculating the mean reconstructed emission time.
    // meanTime comes from GetClassification("meanTime")
    output->Branch("meanTime",&data.meanTime);

    // earlyTime is a classifier counting the number of PMT hits not being in
    // causal contact. earlyTime comes from GetClassification("earlyFraction");
    output->Branch("earlyTime",&data.earlyTime);

    // preTriggerHits is a classifier counting the number of PMT hits before
    // the event trigger. Comes from GetClassification("preTriggerHits");
    output->Branch("preTriggerHits",&data.preTriggerHits);

    // isoRegions is a classifier describing the anisotropy of an event.
    // isotropyRegion comes from GetClassification("fullChi");
    output->Branch("isotropyRegions",&data.isotropyRegions);

    // latIsoRegions is a classifier describing the anisotropy of an event.
    // latIsotropyRegions comes from GetClassification("fullLatChi");
    output->Branch("latIsotropyRegions",&data.latIsotropyRegions);

    // biPoCumulTimeResid is a classifier describing the likelihood of an event being
    // caused by in-window bipo
    // biPoCumul comes from GetClassification("BiPoCumulTimeResid")
    output->Branch("biPoCumul",&data.biPoCumul);

    // biPoLikelihoodDiff is a classifier describing the likelihood of an event being
    // caused by in-window bipo
    // biPoLikelihood212 and biPoLikelihood214 come from GetClassification("BiPoLikelihoodDiff")
    // using PDFs for BiPo212 and BiPo214 respectively.
    output->Branch("biPoLikelihood212",&data.biPoLikelihood212);
    output->Branch("biPoLikelihood214",&data.biPoLikelihood214);

    // alphaBetaClassifier is a classifier describing the ratio of likelihoods between an
    // event being caused by in-window bipo or in-window two betas
    // alphaBeta212 and alphaBeta214 come from GetClassification("AlphaBetaClassifier")
    // using PDFs for BiPo212 and BiPo214 respectively.
    output->Branch("alphaBeta212",&data.alphaBeta212);
    output->Branch("alphaBeta214",&data.alphaBeta214);

    // positionANN is a fit method, however for the ntuple we take only the radius of the event
    output->Branch("annR",&data.annR);

    // nearAVAngular is a classifier that finds the describes the fraction of the total PMTs hit
    // that are within an angle from the event.  Values are in the range 0 - 1.
    output->Branch("nearAV",&data.nearAV);

    // skyshine classifier output
    output->Branch("skyShine",&data.skyShine);
    output->Branch("topOWLCount",&data.topOWLCount);

    // ext0Nu are classifiers offering pulse shape discrimination between 0nu and externals
    // based on event topology and hit time residuals.
    output->Branch("ext0NuAngleTl208AV",&data.ext0NuAngleTl208AV);
    output->Branch("ext0NuAngleTl208HDR",&data.ext0NuAngleTl208HDR);
    output->Branch("ext0NuAngleTl208Exwater",&data.ext0NuAngleTl208ExWater);
    output->Branch("ext0NuAnglePMTbg",&data.ext0NuAnglePMTbg);
    output->Branch("ext0NuTimeTl208AV",&data.ext0NuTimeTl208AV);
    output->Branch("ext0NuTimeTl208HDR",&data.ext0NuTimeTl208HDR);
    output->Branch("ext0NuTimeTl208ExWater",&data.ext0NuTimeTl208ExWater);
    output->Branch("ext0NuTimePMTbg",&data.ext0NuTimePMTbg);
    output->Branch("ext0NuTimeTl208AVNaive",&data.ext0NuTimeTl208AVNaive);
    output->Branch("ext0NuTimeTl208HDRNaive",&data.ext0NuTimeTl208HDRNaive);
    output->Branch("ext0NuTimeTl208ExWaterNaive",&data.ext0NuTimeTl208ExWaterNaive);
    output->Branch("ext0NuTimePMTbgNaive",&data.ext0NuTimePMTbgNaive);
    output->Branch("ext0NuTimeTl208AVOffset",&data.ext0NuTimeTl208AVOffset);
    output->Branch("ext0NuTimeTl208HDROffset",&data.ext0NuTimeTl208HDROffset);
    output->Branch("ext0NuTimeTl208ExWaterOffset",&data.ext0NuTimeTl208ExWaterOffset);
    output->Branch("ext0NuTimePMTbgOffset",&data.ext0NuTimePMTbgOffset);

    output->SetAutoSave(fAutoSave);
}


Processor::Result OutNtupleProc::DSEvent(DS::Run&, DS::Entry& ds)
{
    std::string fittype;
    if (!file) {
      OpenFile();
    }

    memset(&mcdata,0,sizeof(mcDS));
    *branchDS = ds;

// Handle run-level data
    mcdata.runID = ds.GetRunID();
    mcdata.ratversion = RAT::GetRATVersion();

// Handle MC-level data
    if (ds.MCExists()){
        mcdata.mcvalid = 1;
        DS::MC& rmc = ds.GetMC();
        mcdata.Edep = rmc.GetScintEnergyDeposit();
        mcdata.EdepQuenched = rmc.GetScintQuenchedEnergyDeposit();
        mcdata.NCherPhotons = rmc.GetNCherPhotons();
        mcdata.mct = rmc.GetMCTime();

// Collect data on two highest-energy particles
        const size_t mcpc = rmc.GetMCParticleCount();
        std::pair<int,int> index(-1,-1);
        if(mcpc == 1)
            index = std::pair<int,int>(0, -1);
        else if(mcpc > 1)
            index = SearchEnergies(rmc);
        GetParticleData(index, rmc);
    }

// Handle Trigger-level data
    size_t evc = ds.GetEVCount();
    if (evc != 0){
        for (size_t j = 0; j < evc; j++){
            memset(&data,0,sizeof(TriggerDS));
            fittype.clear();
            DS::EV& rEV = ds.GetEV(j);
            EvIndex = j;
            GetEventData(&rEV);

// Limited check of fit results
            fittype = CheckFits(&rEV);
            if(!fittype.empty())
                GetFitData(fittype, &rEV);
            GetPenergyData(&rEV);
            GetMultiPathData(&rEV);
            GetFastZData(&rEV);

// Limited check of classifier results
            try{ data.itr = rEV.
                          GetClassifierResult("ITR:"+fittype)
                                           .GetClassification("ITR");}
            catch(DS::DataNotFound & error)
                {data.itr = DS::INVALID;}
           try{ data.xsite = rEV.
                 GetClassifierResult("XSite:"+fittype)
                         .GetClassification("XSite");}
           catch(DS::DataNotFound & error)
                {data.xsite = DS::INVALID;}
            try{ data.timingPeaks = rEV.
                          GetClassifierResult("timingPeaks:"+fittype)
                                           .GetClassification("timingPeaks");}
            catch(DS::DataNotFound & error)
                {data.timingPeaks = DS::INVALID;}
            try{ data.thetaij = rEV.
                          GetClassifierResult("isotropy:"+fittype)
                                   .GetClassification("thetaij");}
            catch(DS::DataNotFound & error)
                {data.thetaij = DS::INVALID;}
            try{ data.beta14 = rEV.
                          GetClassifierResult("isotropy:"+fittype)
                                   .GetClassification("snobeta14");}
            catch(DS::DataNotFound & error)
                {data.beta14 = DS::INVALID;}
            try{ data.qpdtProbability = rEV.
                          GetClassifierResult("QPDT:"+fittype)
                                   .GetClassification("QPDT");
                 data.qpdtNhitsEarly = rEV.
                          GetClassifierResult("QPDT:"+fittype)
                                   .GetClassification("NHits_early");
                 data.qpdtQMax = rEV.
                          GetClassifierResult("QPDT:"+fittype)
                                   .GetClassification("QMax");
            }
            catch(DS::DataNotFound & error){
                 data.qpdtProbability = DS::INVALID;
                 data.qpdtNhitsEarly = DS::INVALID;
                 data.qpdtQMax = DS::INVALID;
            }
            try{ data.berkeleyAlphaBeta = rEV.
                          GetClassifierResult("BerkeleyAlphaBeta:"+fittype)
                                           .GetClassification("likelihood");}
            catch(DS::DataNotFound & error)
                {
                  try{ data.berkeleyAlphaBeta = rEV.
                          GetClassifierResult("BerkeleyAlphaBeta:default")
                                           .GetClassification("likelihood");}
                  catch(DS::DataNotFound & error)
                    {
                      data.berkeleyAlphaBeta = DS::INVALID;
                    }
                }
            try{ data.meanTime = rEV.
                          GetClassifierResult("meanTime:"+fittype)
                                           .GetClassification("meanTime");}
            catch(DS::DataNotFound & error)
                {data.meanTime = DS::INVALID;}
            try{ data.preTriggerHits = rEV.
                          GetClassifierResult("preTriggerHits:"+fittype)
                                           .GetClassification("preTriggerHits");}
            catch(DS::DataNotFound & error)
                {data.preTriggerHits = DS::INVALID;}
            try{ data.isotropyRegions = rEV.
                          GetClassifierResult("isoRegions:"+fittype)
                                           .GetClassification("fullChi");}
            catch(DS::DataNotFound & error)
                {data.isotropyRegions = DS::INVALID;}
            try{ data.latIsotropyRegions = rEV.
                          GetClassifierResult("isoRegions:"+fittype)
                                           .GetClassification("fullLatChi");}
            catch(DS::DataNotFound & error)
                {data.latIsotropyRegions = DS::INVALID;}
            try{ data.earlyTime = rEV.
                          GetClassifierResult("earlyTime:"+fittype)
                                           .GetClassification("earlyFraction");}
            catch(DS::DataNotFound & error)
                {data.earlyTime = DS::INVALID;}

            try{data.biPoCumul = rEV.
                          GetClassifierResult("BiPoCumulTimeResid:"+fittype)
                                           .GetClassification("gamma");}
            catch(DS::DataNotFound & error)
                {data.biPoCumul = DS::INVALID;}
            try{data.biPoLikelihood212 = rEV.
                          GetClassifierResult("BiPoLikelihoodDiff:212:"+fittype)
                                           .GetClassification("lglklhdDiff");}
            catch(DS::DataNotFound & error)
                {data.biPoLikelihood212 = DS::INVALID;}
            try{data.biPoLikelihood214 = rEV.
                          GetClassifierResult("BiPoLikelihoodDiff:214:"+fittype)
                                           .GetClassification("lglklhdDiff");}
            catch(DS::DataNotFound & error)
                {data.biPoLikelihood214 = DS::INVALID;}
            try{data.alphaBeta212 = rEV.
                          GetClassifierResult("AlphaBetaClassifier:212:"+fittype)
                                           .GetClassification("AlphaBetaClassifier");}
            catch(DS::DataNotFound & error)
                {data.alphaBeta212 = DS::INVALID;}
            try{data.alphaBeta214 = rEV.
                          GetClassifierResult("AlphaBetaClassifier:214:"+fittype)
                                           .GetClassification("AlphaBetaClassifier");}
            catch(DS::DataNotFound & error)
                {data.alphaBeta214 = DS::INVALID;}
            try{data.annR = rEV.
                GetFitResult("positionANN:"+fittype).GetVertex(0).GetPosition().Mag();}
            catch(DS::DataNotFound & error)
                {data.annR = DS::INVALID;}
            try{data.nearAV = rEV.
                GetClassifierResult("nearAVAngular:"+fittype).GetClassification("ratio");}
            catch(DS::DataNotFound & error)
                {data.nearAV = DS::INVALID;}

            try {
              data.skyShine = rEV.GetClassifierResult("skyshine")
                                 .GetClassification("sideFraction");
              data.topOWLCount = rEV.GetClassifierResult("skyshine")
                                    .GetClassification("topOWLCount");
            } catch (DS::DataNotFound& error) {
              data.skyShine = DS::INVALID;
              data.topOWLCount = DS::INVALID;
            }
            try{data.ext0NuAngleTl208AV = rEV.
                          GetClassifierResult("Ext0NuCosTheta:Tl208AV:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuAngleTl208AV = DS::INVALID;}
            try{data.ext0NuAngleTl208HDR = rEV.
                          GetClassifierResult("Ext0NuCosTheta:Tl208Hdropes:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuAngleTl208HDR = DS::INVALID;}
            try{data.ext0NuAngleTl208ExWater = rEV.
                          GetClassifierResult("Ext0NuCosTheta:Tl208Exwater:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuAngleTl208ExWater = DS::INVALID;}
            try{data.ext0NuAnglePMTbg = rEV.
                          GetClassifierResult("Ext0NuCosTheta:PMTbetagamma:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuAnglePMTbg = DS::INVALID;}
            try{data.ext0NuTimeTl208AV = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208AV:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208AV = DS::INVALID;}
            try{data.ext0NuTimeTl208HDR = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Hdropes:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208HDR = DS::INVALID;}
            try{data.ext0NuTimeTl208ExWater = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Exwater:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208ExWater = DS::INVALID;}
            try{data.ext0NuTimePMTbg = rEV.
                          GetClassifierResult("Ext0NuTimeRes:PMTbetagamma:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimePMTbg = DS::INVALID;}
            try{data.ext0NuTimeTl208AVNaive = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208AV_NAIVE:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208AVNaive = DS::INVALID;}
            try{data.ext0NuTimeTl208HDRNaive = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Hdropes_NAIVE:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208HDRNaive = DS::INVALID;}
            try{data.ext0NuTimeTl208ExWaterNaive = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Exwater_NAIVE:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208ExWaterNaive = DS::INVALID;}
            try{data.ext0NuTimePMTbgNaive = rEV.
                          GetClassifierResult("Ext0NuTimeRes:PMTbetagamma_NAIVE:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimePMTbgNaive = DS::INVALID;}
            try{data.ext0NuTimeTl208AVOffset = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208AV_OFFSET10:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208AVOffset = DS::INVALID;}
            try{data.ext0NuTimeTl208HDROffset = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Hdropes_OFFSET10:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208HDROffset = DS::INVALID;}
            try{data.ext0NuTimeTl208ExWaterOffset = rEV.
                          GetClassifierResult("Ext0NuTimeRes:Tl208Exwater_OFFSET10:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimeTl208ExWaterOffset = DS::INVALID;}
            try{data.ext0NuTimePMTbgOffset = rEV.
                          GetClassifierResult("Ext0NuTimeRes:PMTbetagamma_OFFSET10:"+fittype)
                                           .GetClassification("discriminant");}
            catch(DS::DataNotFound & error)
                {data.ext0NuTimePMTbgOffset = DS::INVALID;}

            // Record the collected data to output file
            if(-1 == output->Fill()) Log::Die("Could not fill ntuple file!");
        }
    }
    else{
        EvIndex = -1;
        memset(&data,0,sizeof(TriggerDS));
        // Record the collected data to output file
        if(-1 == output->Fill()) Log::Die("Could not fill ntuple file!");
    }
    McIndex++;
    return Processor::OK;
}


// Collects Event-level data.
void OutNtupleProc::GetEventData(DS::EV *rEV)
{
    data.eventID = rEV->GetGTID();
    data.nhits = rEV->GetNhits();
    data.owlnhits = rEV->GetUncalPMTs().GetOWLCount();
    data.necknhits = rEV->GetUncalPMTs().GetNeckCount();
    data.nhitsCleaned = rEV->GetNhitsCleaned();
    data.inTimeHits100 = rEV->GetInTimeHits100();
    data.inTimeHits20 = rEV->GetInTimeHits20();
    data.isCal = rEV->GetCalibrationEvent();
    data.Q = rEV->GetTotalCharge();
    data.trigger = rEV->GetTrigType();
    data.tubii = rEV->GetTubiiTrig();
    Int_t fPass = rEV->GetDataCleaningFlags().GetLatestPass();
    UInt_t fecd = 0;
    for(size_t iPMT = 0; iPMT < rEV->GetUncalPMTs().GetFECDCount(); iPMT++)
    {
      int channel = rEV->GetUncalPMTs().GetFECDPMT(iPMT).GetChannel();
      fecd += (1<<channel);
    }
    data.fecdHits = fecd;
    if(fPass >= 0)
    {
      data.dcApplied = rEV->GetDataCleaningFlags().GetApplied(fPass).GetULong64_t(0);
      data.dcFlagged = rEV->GetDataCleaningFlags().GetFlags(fPass).GetULong64_t(0);
    }
    data.ClockCount10 = rEV->GetClockCount10();
    data.ClockCount50 = rEV->GetClockCount50();
    data.UTDays = rEV->GetUniversalTime().GetDays();
    data.UTSecs = rEV->GetUniversalTime().GetSeconds();
    data.UTNSecs = rEV->GetUniversalTime().GetNanoSeconds();
}

// Sets fittype and checks fit-specific data.  Note that only one fit result
// is used even if several are available.  The order of preference is:
// scintFitter, then waterFitter, then partialFitter and partialWaterFitter,
// then quad fitter.
std::string OutNtupleProc::CheckFits(DS::EV *rEV)
{
    std::string fittype;
    try{
        data.fitValid = rEV->GetFitResult("scintFitter").GetValid();
        data.scintFit = true;
        fittype = "scintFitter";
    }
    catch( DS::DataNotFound & error ){}
    if (fittype.empty()){
        try{
            data.fitValid = rEV->GetFitResult("waterFitter").GetValid();
            data.waterFit = true;
            fittype = "waterFitter";
            GetWaterFitterData(fittype, rEV);
        }
        catch( DS::DataNotFound & error ){}
    }
    if (fittype.empty()){
        try{
            data.fitValid = rEV->GetFitResult("partialFitter").GetValid();
            data.partialFit = true;
            fittype = "partialFitter";
        }
        catch( DS::DataNotFound & error ){}
    }
    if (fittype.empty()){
        try{
            data.fitValid = rEV->GetFitResult("partialWaterFitter").GetValid();
            data.partialWaterFit = true;
            fittype = "partialWaterFitter";
            GetWaterFitterData(fittype, rEV);
        }
        catch( DS::DataNotFound & error ){}
    }
    if (fittype.empty()){
        try{
            data.fitValid = rEV->GetFitResult("quad").GetValid();
            fittype = "quad";
        }
        catch( DS::DataNotFound & error ){}
    }
    return fittype;
}

// Looks through Monte Carlo data to find the two particles of highest energies.
std::pair<int,int> OutNtupleProc::SearchEnergies(DS::MC& rmc)
{
    int mcpc = rmc.GetMCParticleCount();
    std::pair<int,int> index(0,0);
    std::pair<float,float> energy(0,0);
    DS::MCParticle& rparticle0 = rmc.GetMCParticle(0);
    DS::MCParticle& rparticle1 = rmc.GetMCParticle(1);
    float en0 = rparticle0.GetKineticEnergy();
    float en1 = rparticle1.GetKineticEnergy();
    if( en0 > en1 ){
        energy.first = en0;
        energy.second = en1;
        index.second = 1;
    }
    else{
        energy.first = en1;
        energy.second = en0;
        index.first = 1;
    }
    for (int pIndex = 2; pIndex < mcpc; pIndex++){
        DS::MCParticle& rparticle = rmc.GetMCParticle(pIndex);
        if (rparticle.GetKineticEnergy() > energy.first){
            index.second = index.first;
            index.first = pIndex;
        }
        else{
            if (rparticle.GetKineticEnergy() > energy.second)
                index.second = pIndex;
        }
    }
    return index;
}

// Collects some standardizable fit results
void OutNtupleProc::GetFitData(const std::string fittype, DS::EV *rEV)
{

    DS::FitResult fResult = rEV->GetFitResult(fittype);
    DS::FitVertex rvertex = fResult.GetVertex(0);

    // Position information
    if( rvertex.ContainsPosition() && rvertex.ValidPosition() ) {
        data.posx = rvertex.GetPosition().X();
        data.posy = rvertex.GetPosition().Y();
        data.posz = rvertex.GetPosition().Z();
        data.posr = rvertex.GetPosition().Mag();
    }
    else {
        data.posx = DS::INVALID;
        data.posy = DS::INVALID;
        data.posz = DS::INVALID;
        data.posr = DS::INVALID;
    }
    if( rvertex.ContainsPositivePositionError() && rvertex.ValidPositivePositionError() ) {
        data.posxPosError = rvertex.GetPositivePositionError().X();
        data.posyPosError = rvertex.GetPositivePositionError().Y();
        data.poszPosError = rvertex.GetPositivePositionError().Z();
    }
    else {
        data.posxPosError = DS::INVALID;
        data.posyPosError = DS::INVALID;
        data.poszPosError = DS::INVALID;
    }
    if( rvertex.ContainsNegativePositionError() && rvertex.ValidNegativePositionError() ) {
        data.posxNegError = rvertex.GetNegativePositionError().X();
        data.posyNegError = rvertex.GetNegativePositionError().Y();
        data.poszNegError = rvertex.GetNegativePositionError().Z();
    }
    else {
        data.posxNegError = DS::INVALID;
        data.posyNegError = DS::INVALID;
        data.poszNegError = DS::INVALID;
    }

    try {
      data.posFOM  = fResult.GetFOM("PositionLogL");
      data.posFOM2 = (UInt_t)fResult.GetFOM("PositionSelectedNHit");
    } catch(const std::out_of_range& oor) {
      data.posFOM  = DS::INVALID;
      data.posFOM2 = DS::INVALID;
    }

    // Energy information
    if( rvertex.ContainsEnergy() && rvertex.ValidEnergy() )
        data.energy = rvertex.GetEnergy();
    else
        data.energy = DS::INVALID;
    if( rvertex.ContainsPositiveEnergyError() && rvertex.ValidPositiveEnergyError() )
        data.energyPosError = rvertex.GetPositiveEnergyError();
    else
        data.energyPosError = DS::INVALID;
    if( rvertex.ContainsNegativeEnergyError() && rvertex.ValidNegativeEnergyError() )
        data.energyNegError = rvertex.GetNegativeEnergyError();
    else
        data.energyNegError = DS::INVALID;

    try{ data.energyFOMmedProb = fResult.GetFOM("EnergyMedianProb");}
    catch(const std::out_of_range& oor) { data.energyFOMmedProb = DS::INVALID; }
    try{ data.energyFOMmedProbHit = fResult.GetFOM("EnergyMedianProbHit");}
    catch(const std::out_of_range& oor) { data.energyFOMmedProbHit = DS::INVALID; }
    try{ data.energyFOMmedDev = fResult.GetFOM("EnergyMedianDev");}
    catch(const std::out_of_range& oor) { data.energyFOMmedDev = DS::INVALID; }
    try{ data.energyFOMmedDevHit = fResult.GetFOM("EnergyMedianDevHit");}
    catch(const std::out_of_range& oor) { data.energyFOMmedDevHit = DS::INVALID; }
    try{ data.energyFOMGtest = fResult.GetFOM("EnergyGtest");}
    catch(const std::out_of_range& oor) { data.energyFOMGtest = DS::INVALID; }
    try{ data.energyFOMUtest = fResult.GetFOM("EnergyUtest");}
    catch(const std::out_of_range& oor) { data.energyFOMUtest = DS::INVALID; }

    // Time information
    if( rvertex.ContainsTime() && rvertex.ValidTime() )
        data.time = rvertex.GetTime();
    else
        data.time = DS::INVALID;
    if( rvertex.ContainsPositiveTimeError() && rvertex.ValidPositiveTimeError() )
        data.timePosError = rvertex.GetPositiveTimeError();
    else
        data.timePosError = DS::INVALID;
    if( rvertex.ContainsNegativeTimeError() && rvertex.ValidNegativeTimeError() )
        data.timeNegError = rvertex.GetNegativeTimeError();
    else
        data.timeNegError = DS::INVALID;

    // Overall validity
    data.fitValidityMask = rvertex.GetValidityMask();
}

void OutNtupleProc::GetPenergyData(DS::EV *rEV){
    try{
        DS::FitResult fitResult = rEV->GetFitResult("penergy");
        data.penergyEnergy = fitResult.GetVertex(0).GetEnergy();
        data.penergyError  = fitResult.GetVertex(0).GetPositiveEnergyError();
    } catch ( DS::DataNotFound & error ){
        data.penergyEnergy = DS::INVALID;
    }

    try{
        DS::FitResult fitResult = rEV->GetFitResult("penergyA");
        data.penergyAEnergy = fitResult.GetVertex(0).GetEnergy();
        data.penergyAError  = fitResult.GetVertex(0).GetPositiveEnergyError();
    } catch ( DS::DataNotFound & error ){
        data.penergyAEnergy = DS::INVALID;
    }
}

void OutNtupleProc::GetMultiPathData(DS::EV *rEV){

    data.posx_MultiPath = DS::INVALID;
    data.posy_MultiPath = DS::INVALID;
    data.posz_MultiPath = DS::INVALID;
    data.posr_MultiPath = DS::INVALID;
    data.time_MultiPath = DS::INVALID;
    data.fitValid_MultiPath = false;

    try{
        DS::FitVertex rvertex = rEV->GetFitResult("multipath").GetVertex(0);
        if( rvertex.ContainsPosition() && rvertex.ValidPosition() ) {
            data.posx_MultiPath = rvertex.GetPosition().X();
            data.posy_MultiPath = rvertex.GetPosition().Y();
            data.posz_MultiPath = rvertex.GetPosition().Z();
            data.posr_MultiPath = rvertex.GetPosition().Mag();
            data.time_MultiPath = rvertex.GetTime();
            data.fitValid_MultiPath = rvertex.ValidPosition();
        }
    }
    catch ( DS::DataNotFound & error ){ /*Catch this case and pass on initialized (invalid) values*/ }
}

void OutNtupleProc::GetFastZData(DS::EV *rEV){

    data.posz_fastZ = DS::INVALID;
    data.time_fastZ = DS::INVALID;
    data.fitValid_fastZ = false;

    try{
        DS::FitVertex rvertex = rEV->GetFitResult("fastZ").GetVertex(0);
        if( rvertex.ContainsPosition() && rvertex.ValidPosition() ) {
            data.posz_fastZ = rvertex.GetPosition().Z();
            data.time_fastZ = rvertex.GetTime();
            data.fitValid_fastZ = rvertex.ValidPosition();
        }
    }
    catch ( DS::DataNotFound & error ){ /*Catch this case and pass on initialized (invalid) values*/ }
}


// Collects fit results specific to waterFitter
void OutNtupleProc::GetWaterFitterData(const std::string fittype, DS::EV *rEV)
{

    DS::FitResult fResult = rEV->GetFitResult(fittype);
    DS::FitVertex rvertex = fResult.GetVertex(0);

    if( rvertex.ContainsDirection() && rvertex.ValidDirection() ) {
        data.dirx = rvertex.GetDirection().X();
        data.diry = rvertex.GetDirection().Y();
        data.dirz = rvertex.GetDirection().Z();
    }
    else {
        data.dirx = DS::INVALID;
        data.diry = DS::INVALID;
        data.dirz = DS::INVALID;
    }
    if( rvertex.ContainsPositiveDirectionError() && rvertex.ValidPositiveDirectionError() ) {
        data.dirxPosError = rvertex.GetPositiveDirectionError().X();
        data.diryPosError = rvertex.GetPositiveDirectionError().Y();
        data.dirzPosError = rvertex.GetPositiveDirectionError().Z();
    }
    else {
        data.dirxPosError = DS::INVALID;
        data.diryPosError = DS::INVALID;
        data.dirzPosError = DS::INVALID;
    }
    if( rvertex.ContainsNegativeDirectionError() && rvertex.ValidNegativeDirectionError() ) {
        data.dirxNegError = rvertex.GetNegativeDirectionError().X();
        data.diryNegError = rvertex.GetNegativeDirectionError().Y();
        data.dirzNegError = rvertex.GetNegativeDirectionError().Z();
    }
    else {
        data.dirxNegError = DS::INVALID;
        data.diryNegError = DS::INVALID;
        data.dirzNegError = DS::INVALID;
    }

    try {
      data.dirFOM  = fResult.GetFOM("DirectionLogL");
      data.dirFOM2 = (UInt_t)fResult.GetFOM("DirectionSelectedNHit");
    } catch(const std::out_of_range& oor) {
      data.dirFOM  = DS::INVALID;
      data.dirFOM2 = DS::INVALID;
    }
}

// Collects the desired data on the two particles of highest energy.
void OutNtupleProc::GetParticleData(const std::pair<int,int> index, DS::MC& rmc)
{

  if(index.first == -1) return;
  const DS::MCParticle& part1 = rmc.GetMCParticle(index.first);

  mcdata.pdg1 = part1.GetPDGCode();
  mcdata.ke1 = part1.GetKineticEnergy();
  mcdata.mcposx = part1.GetPosition().X();
  mcdata.mcposy = part1.GetPosition().Y();
  mcdata.mcposz = part1.GetPosition().Z();
  mcdata.mcposr = part1.GetPosition().Mag();
  mcdata.mom1x = part1.GetMomentum().X();
  mcdata.mom1y = part1.GetMomentum().Y();
  mcdata.mom1z = part1.GetMomentum().Z();
  mcdata.time1 = part1.GetTime();

  // -- if there is a second particle grab it as well
  bool has_second = (index.second != -1);

  if(has_second){
    const DS::MCParticle& part2 = rmc.GetMCParticle(index.second);
    mcdata.pdg2 = part2.GetPDGCode();
    mcdata.ke2 = part2.GetKineticEnergy();
    mcdata.mom2x = part2.GetMomentum().X();
    mcdata.mom2y = part2.GetMomentum().Y();
    mcdata.mom2z = part2.GetMomentum().Z();
    mcdata.time2 = part2.GetTime();
  }

  // -- How grab the parents and loop for the information
  size_t mcpar = rmc.GetMCParentCount();
  // The previous logic was flawed. If, for some reason the most energetic particle
  // was not the first in terms of index, this would not give the right result.
  // This could, for example, happen in a BiPo decay with a delayed electron
  for (size_t ipar = 0; ipar < mcpar; ++ipar) {
    const DS::MCParticle& rparent = rmc.GetMCParent(ipar);

    if (rparent.GetVertexID() == part1.GetVertexID()) {
        mcdata.parentpdg1 = rparent.GetPDGCode();
        mcdata.parentKE1 = rparent.GetKineticEnergy();
        mcdata.parentMetaInfo1 = rparent.GetMetaInfo();
    }
    if (has_second) {
      if (rparent.GetVertexID() == rmc.GetMCParticle(index.second).GetVertexID()) {
          mcdata.parentpdg2 = rparent.GetPDGCode();
          mcdata.parentKE2 = rparent.GetKineticEnergy();
          mcdata.parentMetaInfo2 = rparent.GetMetaInfo();
      }
    }
  }
}

} // namespace RAT