ParametersNumu.cxx

Go to the documentation of this file.

//
//  ParametersNumu.cpp
//  nova-mrb
//
//  Created by Brian Rebel on 10/29/15.
//  Copyright © 2015 Brian Rebel. All rights reserved.
//
#include <bitset>

#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/CellGeo.h"
#include "AnalysisSkimmer/parameters/ParametersNumu.h"
#include "RunHistory/service/RunHistoryService.h"

#include "StandardRecord/CVNLabels.h"

//------------------------------------------------------------------------------
std::pair<int, int> calcFirstLastLivePlane(int             plane,
                                           std::bitset<14> binary){
  //                                       novadaq::cnv::DetId det){
  //if not the FD, return full ND size                                                                                                                                                
  //if (det != novadaq::cnv::kFARDET) return std::pair<int, int>(0,213);
  int testDB = (plane / 64);
  int minblock = testDB;
  int maxblock = testDB;
  //find block boundaries;                                                                                                                                                            
  for (int i = testDB - 1; i >=0 ; --i){
    if (binary[i]) minblock = i;
    else break;
  }
  for (int i = testDB + 1; i < 14; ++i){
    if (binary[i]) maxblock = i;
    else break;
  }
  return std::pair<int, int>(64 * (minblock), 64 * (maxblock + 1) - 1);
}

//------------------------------------------------------------------------------
skim::ParametersNumu::ParametersNumu()
: fEdgeXPlusXZ (0.)
, fEdgeXMinusXZ(0.)
, fEdgeFrontXZ (0.)
, fEdgeBackXZ  (0.)
, fEdgeTopYZ   (0.)
, fEdgeBottomYZ(0.)
, fEdgeFrontYZ (0.)
, fEdgeBackYZ  (0.)
{
  art::ServiceHandle<geo::Geometry> geom;
  
  fNumCells    = geom->Plane(0)->Ncells() - 1;
  fNumDiblocks = std::ceil(geom->NPlanes() / 64);
}

//------------------------------------------------------------------------------
skim::ParametersNumu::ParametersNumu(size_t            const& bestPIDTrack,
                                     rb::Energy        const& energy,
                                     remid::ReMId      const& remid,
                                     cosrej::CosRejObj const& cosrej,
                                     cosrej::NueCosRej const& nuecosrej,
                                     rb::Vertex        const& vertex,
                                     std::vector<rb::Track  const*>           & tracks,
                                     std::vector<rb::Shower const*>           & showers,
                                     std::vector<rb::Prong  const*>           & prongs,
                                     std::vector<slid::ShowerLID const*>      & slids,
                                     numue::NumuE      const& numue,
                                     rb::Cluster       const& slice,
                                     rb::Track         const& track,
                                     cvn::Result       const& cvn,
                                     cvn::Result       const& cvn2017,
                                     bool              const& isData)
: ParametersNumu()
{

  art::ServiceHandle<geo::Geometry>                geom;
  art::ServiceHandle<nova::dbi::RunHistoryService> rh;
  
  fParameters.fBestPIDTrack = 1.*bestPIDTrack;

  // Energy Variables
  fParameters.fTrackCCE  = numue.TrkCCE();
  fParameters.fNeutrinoE = fParameters.fTrackCCE;
  
  fParameters.fNDHadronicCal = numue.NDHadCalTran() + numue.NDHadCalCat();
  
  // ReMId variable
  fParameters.fReMIdVal = remid.Value();
  
  // Slice variables
  fParameters.fSliceHits         = slice.NCell();
  fParameters.fSliceContigPlanes = slice.MostContiguousPlanes(geo::kXorY);
  fParameters.fSliceMaxPlane     = slice.MaxPlane();
  fParameters.fSliceMinPlane     = slice.MinPlane();
  for(int i = 0;i < 10; ++i)fParameters.fSliceNHitsNearEdge[i] = 0;

  // check to see the minimum number of cells from edge of the detector
  fParameters.fSliceCellsFromEdge = 1e6;
  for(auto ch : slice.AllCells() ){
    size_t cellsFromEdge = std::min((unsigned int)ch->Cell(), geom->Plane(ch->Plane())->Ncells() - 1 - ch->Cell());
    LOG_DEBUG("EvaluatorNumu")
    << ch->Cell()
    << " "
    << geom->Plane(ch->Plane())->Ncells()-1-ch->Cell()
    << " "
    << fParameters.fSliceCellsFromEdge;
    
    if( cellsFromEdge < fParameters.fSliceCellsFromEdge ){
      fParameters.fSliceCellsFromEdge = cellsFromEdge;
      if     (ch->View() == geo::kX && cellsFromEdge == ch->Cell())                                          fParameters.fSliceClosestEdge = 0;
      else if(ch->View() == geo::kX && cellsFromEdge == geom->Plane(ch->Plane())->Ncells() - 1 - ch->Cell()) fParameters.fSliceClosestEdge = 1;
      else if(ch->View() == geo::kY && cellsFromEdge == ch->Cell())                                          fParameters.fSliceClosestEdge = 2;
      else if(ch->View() == geo::kY && cellsFromEdge == geom->Plane(ch->Plane())->Ncells() - 1 - ch->Cell()) fParameters.fSliceClosestEdge = 3;
      else fParameters.fSliceClosestEdge = -5;
    }
    
    if(cellsFromEdge < 10) {
      for(unsigned int i = 9; ( i >= cellsFromEdge && i < 10); i -= 1) ++fParameters.fSliceNHitsNearEdge[i];
    }
  }
  
  // CosRej variables
  fParameters.fCosRejKalFwdCell     = cosrej.KalFwdCell();
  fParameters.fCosRejKalBakCell     = cosrej.KalBakCell();
  fParameters.fCosRejAngleKal       = cosrej.AngleKal();
  fParameters.fCosRejCosFwdCell     = cosrej.CosFwdCell();
  fParameters.fCosRejCosBakCell     = cosrej.CosBakCell();
  fParameters.fCosRejConCosPID      = cosrej.ConCosPID();
  fParameters.fCosRejKalBakCellND   = cosrej.KalBakCellND();
  fParameters.fCosRejKalFwdCellND   = cosrej.KalFwdCellND();
  fParameters.fCosRejKalYPosAtTrans = cosrej.KalYPosAtTrans();
  fParameters.fCosRejKalSpeed       = cosrej.KFitSpeed();
  fParameters.fCosRejCosSpeed       = cosrej.CFitSpeed();
  fParameters.fCosRejCosDirScore    = cosrej.CDirScore();
  
  // Track variables
  fParameters.fTrackStartZ    = track.Start().Z();
  fParameters.fTrackStopZ     = track.Stop().Z();
  fParameters.fTrackStart     = track.Start();
  fParameters.fTrackStop      = track.Stop();
  fParameters.fTrackStartTime = track.Cell(0)->TNS();
  fParameters.fTrackStopTime  = track.Cell(track.NCell()-1)->TNS();
  fParameters.fTrackDir       = track.Dir();
  for(int i = 0; i < 10; ++i) fParameters.fTrackNHitsNearEdge[i] = 0;
  fParameters.fTrackEarliestHitPos = EarliestLatestHitPos(track, false).first;
  fParameters.fTrackLatestHitPos   = EarliestLatestHitPos(track, false).second;

  // Vetex variables (2017)
  fVertex                     = vertex.GetXYZ();
  fParameters.fVertexX        = fVertex.X();
  fParameters.fVertexY        = fVertex.Y();
  fParameters.fVertexZ        = fVertex.Z();
  fParameters.fVertexMaxW     = std::max(std::abs(fVertex.X()),
                                         std::abs(fVertex.Y()));

  // Shower variables (2017)
  fParameters.fNumShowers     = 1.*showers.size();
  if(showers.size()>0){
    fParameters.fShowerCalE     = showers[0]->CalorimetricEnergy();
    fParameters.fShowerTotalE   = showers[0]->TotalGeV();
    fShowerDir                  = showers[0]->Dir();
    fParameters.fShowerDirX     = fShowerDir.X();
    fParameters.fShowerDirY     = fShowerDir.Y();
    fParameters.fShowerDirZ     = fShowerDir.Z();
    fShowerVertex               = showers[0]->Start();
    fParameters.fShowerVertexX  = fShowerVertex.X();
    fParameters.fShowerVertexY  = fShowerVertex.Y();
    fParameters.fShowerVertexZ  = fShowerVertex.Z();
    fShowerEnd                  = showers[0]->Stop();
    fParameters.fShowerEndX     = fShowerEnd.X();
    fParameters.fShowerEndY     = fShowerEnd.Y();
    fParameters.fShowerEndZ     = fShowerEnd.Z();

    fParameters.fShowerNumXCell = showers[0]->NXCell();
    fParameters.fShowerNumYCell = showers[0]->NYCell();    
  } // showers.size()

  // Shower LID variables (2017)
  fParameters.fNumShowersLID  = 1.*slids.size();
  if(slids.size()>0){
    fParameters.fShowerVtxDist  = slids[0]->Gap();
  }

  fParameters.fShowerMinX     = std::numeric_limits<double>::max();
  fParameters.fShowerMinY     = std::numeric_limits<double>::max();
  fParameters.fShowerMaxX     = std::numeric_limits<double>::min();
  fParameters.fShowerMaxY     = std::numeric_limits<double>::min();
  fParameters.fShowerMinVtxZ  = std::numeric_limits<double>::max();
  fParameters.fShowerMinEndZ  = std::numeric_limits<double>::max();
  fParameters.fShowerMaxVtxZ  = std::numeric_limits<double>::min();
  fParameters.fShowerMaxEndZ  = std::numeric_limits<double>::min();
  fParameters.fProngLength    = std::numeric_limits<double>::min();

  fParameters.fOtherTrackMaxEndZ = std::numeric_limits<double>::min();
  fParameters.fOtherTrackMaxVtxZ = std::numeric_limits<double>::min();

  fParameters.fCosShowers     = 2.; // default value will pass any cuts


  //////////////////////// CVN ///////////////////////////
  // Differnet versions have different number of outputs
  //
  // Explanation for Production4 is in :
  //     https://cdcvs.fnal.gov/redmine/projects/nova_reconstruction/wiki/Prod4_tag

  
  fParameters.fCVN2017Val      = cvn2017.fOutput[0]+cvn2017.fOutput[1]+cvn2017.fOutput[2]+cvn2017.fOutput[3];

  // check Output Vector Size
  int cvnOutputSize = cvn.fOutput.size();

  if( cvnOutputSize == 16 ) {         /// 2017 and before

    // CVN Value formula from CAFMaker: FillPIDs.cxx -->
    // srcvn.numuid = result.fOutput[0]+result.fOutput[1]+result.fOutput[2]+result.fOutput[3];)
    fParameters.fCVNVal         = cvn.fOutput[0]+cvn.fOutput[1]+cvn.fOutput[2]+cvn.fOutput[3];

  } else if( cvnOutputSize == 392 ) {  /// 2018 ...
                                       /// See StandardRecord/CVNLabels.h  for constants

    float numuid = 0;
    //std::vector<float>::const_iterator maxIt = std::max_element(cvn.fOutput.cbegin(), cvn.fOutput.cend());

    for(int idx = 0; idx < caf::kNumCVNFinalStates; ++idx){
      const caf::CVNFinalState state = caf::cvnStates[idx];
      const double prob = cvn.fOutput[idx];
      if(state.nuPdg == 14) numuid += prob;
    }

    ///  srcvn.nueid += cvn.fOutput[caf::kCVN_nE_Other];
    numuid += cvn.fOutput[caf::kCVN_nM_Other];

    fParameters.fCVNVal         =   numuid;
  }

  /////////////////////////////////////////////////////

  float nShowerHit = 0.;
  if(showers.size() > 1)
    fParameters.fCosShowers = ((fShowerDir * showers[1]->Dir())/
                               (fShowerDir.Mag() * showers[1]->Dir().Mag()));

  // The ProngLength cut applies for the longest prong, so find that in
  // the following loop

  for(auto prg : prongs){
    if(prg->TotalLength() > fParameters.fProngLength)
      fParameters.fProngLength = prg->TotalLength();
  }

  for(auto shw : showers){
    nShowerHit += 1.*shw->NCell();
    auto vtx    = shw->Start();
    auto end    = shw->Stop();
    auto maxX   = std::max(vtx.X(), end.X());
    auto maxY   = std::max(vtx.Y(), end.Y());
    auto minX   = std::min(vtx.X(), end.X());
    auto minY   = std::min(vtx.Y(), end.Y());
    if(fParameters.fShowerMinX    > minX   ) fParameters.fShowerMinX    = minX;
    if(fParameters.fShowerMinY    > minY   ) fParameters.fShowerMinY    = minY;
    if(fParameters.fShowerMaxX    < maxX   ) fParameters.fShowerMaxX    = maxX;
    if(fParameters.fShowerMaxY    < maxY   ) fParameters.fShowerMaxY    = maxY;
    if(fParameters.fShowerMinVtxZ > vtx.Z()) fParameters.fShowerMinVtxZ = vtx.Z();
    if(fParameters.fShowerMinEndZ > end.Z()) fParameters.fShowerMinEndZ = end.Z();
    if(fParameters.fShowerMaxVtxZ < vtx.Z()) fParameters.fShowerMaxVtxZ = vtx.Z();
    if(fParameters.fShowerMaxEndZ < end.Z()) fParameters.fShowerMaxEndZ = end.Z();
  }


  // Check start and end Z for tracks other then best RemId one
  for(size_t itrk=0; itrk<tracks.size(); ++itrk){

    if( itrk == bestPIDTrack )             continue;  // Skip best pid track
    if(tracks[itrk]->View() != geo::kXorY) continue;  // 3D tracks only as in CAFMaker
    auto vtx    = tracks[itrk]->Start();
    auto end    = tracks[itrk]->Stop();
    if( fParameters.fOtherTrackMaxEndZ < end.Z() ) fParameters.fOtherTrackMaxEndZ = end.Z(); 
    if( fParameters.fOtherTrackMaxVtxZ < vtx.Z() ) fParameters.fOtherTrackMaxVtxZ = vtx.Z(); 
  }

  ///----------------------------------------------------------------------------

  // now check the run history and the diblocks in this run
  // the lines below are copied from CAFMaker_module.cc
  unsigned int dibmask = (rh->GoodDiBlockMask()&rh->GetConfiguration());
  
  unsigned int dibfirst     = 0;
  //unsigned int diblast      = 0;
  unsigned int firstDBPlane = 0;
  //unsigned int lastDBPlane  = geom->NPlanes();
  
  if (dibmask) {
    int iD;
    iD=0; while (!((dibmask>>iD)&1)) iD++; dibfirst=iD+1;
    //iD=0; while (dibmask>>iD)        iD++; diblast=iD;
    
    firstDBPlane = (dibfirst - 1) * 64;
    //lastDBPlane  = (diblast * 64) - 1;
  }
  
  std::bitset<14> binary(dibmask);
  
  std::pair<int,int> planesA = calcFirstLastLivePlane(slice.MinPlane(), binary);
  std::pair<int,int> planesB = calcFirstLastLivePlane(slice.MaxPlane(), binary);
  
  if ((planesA.first != planesB.first) || (planesA.second != planesB.second))fParameters.fDiffLastFirstLiveDiblock = -1;
  else fParameters.fDiffLastFirstLiveDiblock = (planesA.second - planesA.first + 1) * 0.015625;

  fParameters.fPlanesToFront = slice.MinPlane() - planesA.first;
  fParameters.fPlanesToBack  = planesA.second - slice.MaxPlane();
  
  LOG_DEBUG("ParametersNumu")
    << "Inside: ParametersNumu.cxx--"
    << " FirstPlane: "
    << fParameters.fSliceMinPlane
    << " firstDBPlane "
    << firstDBPlane
    << " planesA.first "
    << planesA.first;

  // give cosmic track related parameters obviously stupid values
  for(int i = 0; i < 10; ++i) fParameters.fCosTrackNHitsNearEdge[i] = std::numeric_limits<int>::min();
  
  fParameters.fCosmicTrackStart                = TVector3(std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min());
  fParameters.fCosmicTrackStop                 = TVector3(std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min());
  fParameters.fCosmicTrackEarliestHitPos       = TVector3(std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min());
  fParameters.fCosmicTrackLatestHitPos         = TVector3(std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min());
  fParameters.fCosmicTrackDir                  = TVector3(std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min(),
                                                          std::numeric_limits<double>::min());
  fParameters.fCosTrackStartTime               = std::numeric_limits<float>::min();
  fParameters.fCosTrackStopTime                = std::numeric_limits<float>::min();
  fParameters.fCosTrackXViewHits               = std::numeric_limits<int  >::min();
  fParameters.fCosTrackYViewHits               = std::numeric_limits<int  >::min();
  fParameters.fCosTrackTotalADC                = std::numeric_limits<float>::min();
  fParameters.fCosTrackTotalGeV                = std::numeric_limits<float>::min();
  fParameters.fCosTrackTotalCalorimetricEnergy = std::numeric_limits<float>::min();
  fParameters.fCosTrackMinTNS                  = std::numeric_limits<float>::min();
  fParameters.fCosTrackMinXYZ                  = TVector3(std::numeric_limits<float>::min(),
                                                          std::numeric_limits<float>::min(),
                                                          std::numeric_limits<float>::min());
  fParameters.fCosTrackMinPlane                = std::numeric_limits<int  >::min();
  fParameters.fCosTrackMaxTNS                  = std::numeric_limits<float>::min();
  fParameters.fCosTrackMaxXYZ                  = TVector3(std::numeric_limits<float>::min(),
                                                          std::numeric_limits<float>::min(),
                                                          std::numeric_limits<float>::min());
  fParameters.fCosTrackMaxPlane                = std::numeric_limits<int>::min();
  

  //Added for 2017 Analysis ...
  fParameters.fDistAllTop     = std::min(nuecosrej.StartMinDistToTop(),
                                         nuecosrej.StopMinDistToTop());
  fParameters.fDistAllBottom  = std::min(nuecosrej.StartMinDistToBottom(),
                                         nuecosrej.StopMinDistToBottom());
  fParameters.fDistAllWest    = std::min(nuecosrej.StartMinDistToWest(),
                                         nuecosrej.StopMinDistToWest());
  fParameters.fDistAllEast    = std::min(nuecosrej.StartMinDistToEast(),
                                         nuecosrej.StopMinDistToEast());
  fParameters.fDistAllFront   = std::min(nuecosrej.StartMinDistToFront(),
                                         nuecosrej.StopMinDistToFront());
  fParameters.fDistAllBack    = std::min(nuecosrej.StartMinDistToBack(),
                                         nuecosrej.StopMinDistToBack());

  fParameters.fNueCosRejPngPtp =  nuecosrej.ProngTransMom();

  return;
}

//------------------------------------------------------------------------------
skim::ParametersNumu::ParametersNumu(size_t            const& bestPIDTrack,
                                     rb::Energy        const& energy,
                                     remid::ReMId      const& remid,
                                     cosrej::CosRejObj const& cosrej,
                                     cosrej::NueCosRej const& nuecosrej,
                                     rb::Vertex        const& vertex,
                                     std::vector<rb::Track  const*>           & tracks,
                                     std::vector<rb::Shower const*>           & showers,
                                     std::vector<rb::Prong const*>            & prongs,
                                     std::vector<slid::ShowerLID const*>      & slids,
                                     rb::Track         const& cosmictrack,
                                     numue::NumuE      const& numue,
                                     rb::Cluster       const& slice,
                                     rb::Track         const& track,
                                     cvn::Result       const& cvn,
                                     cvn::Result       const& cvn2017,
                                     bool              const& isData)
: ParametersNumu(bestPIDTrack,
                 energy,
                 remid,
                 cosrej,
                 nuecosrej,
                 vertex,
                 tracks,
                 showers,
                 prongs,
                 slids,
                 numue,
                 slice,
                 track,
                 cvn,
                 cvn2017,
                 isData)
{
  if(fEdgeXPlusXZ == 0 && fEdgeXMinusXZ == 0){
    GetDetectorEdges();
  }

  for(int i = 0; i < 10; ++i) fParameters.fCosTrackNHitsNearEdge[i] = 0;

  fParameters.fCosmicTrackStart                = cosmictrack.Start();
  fParameters.fCosmicTrackStop                 = cosmictrack.Stop();
  fParameters.fCosmicTrackEarliestHitPos       = EarliestLatestHitPos(cosmictrack, true).first;
  fParameters.fCosmicTrackLatestHitPos         = EarliestLatestHitPos(cosmictrack, true).second;
  fParameters.fCosmicTrackDir                  = cosmictrack.Dir();
  fParameters.fCosTrackStartTime               = cosmictrack.Cell(0)->TNS();
  fParameters.fCosTrackStopTime                = cosmictrack.Cell(cosmictrack.NCell()-1)->TNS();
  fParameters.fCosTrackXViewHits               = cosmictrack.NXCell();
  fParameters.fCosTrackYViewHits               = cosmictrack.NYCell();
  fParameters.fCosTrackTotalADC                = cosmictrack.TotalADC();
  fParameters.fCosTrackTotalGeV                = cosmictrack.TotalGeV();
  fParameters.fCosTrackTotalCalorimetricEnergy = cosmictrack.CalorimetricEnergy();
  fParameters.fCosTrackMinTNS                  = cosmictrack.MinTNS();
  fParameters.fCosTrackMinXYZ                  = cosmictrack.MinXYZ();
  fParameters.fCosTrackMinPlane                = cosmictrack.MinPlane();
  fParameters.fCosTrackMaxTNS                  = cosmictrack.MaxTNS();
  fParameters.fCosTrackMaxXYZ                  = cosmictrack.MaxXYZ();
  fParameters.fCosTrackMaxPlane                = cosmictrack.MaxPlane();
  
  return;
}

//------------------------------------------------------------------------------
skim::ParametersNumu::ParametersNumu(size_t            const& bestPIDTrack,
                                     rb::Energy        const& energy,
                                     remid::ReMId      const& remid,
                                     cosrej::CosRejObj const& cosrej,
                                     cosrej::NueCosRej const& nuecosrej,
                                     rb::Vertex        const& vertex,
                                     std::vector<rb::Track  const*>           & tracks,
                                     std::vector<rb::Shower const*>           & showers,
                                     std::vector<rb::Prong const*>            & prongs,
                                     std::vector<slid::ShowerLID const*>      & slids,
                                     qeef::QePId       const& qepid,
                                     numue::NumuE      const& numue,
                                     rb::Cluster       const& slice,
                                     rb::Track         const& track,
                                     cvn::Result       const& cvn,
                                     cvn::Result       const& cvn2017,
                                     bool              const& isData)
: ParametersNumu(bestPIDTrack,
                 energy,
                 remid,
                 cosrej,
                 nuecosrej,
                 vertex,
                 tracks,
                 showers,
                 prongs,
                 slids,
                 numue,
                 slice,
                 track,
                 cvn,
                 cvn2017,
                 isData)
{
  // QePID variables
  fParameters.fQePIDVal     = qepid.Value();
  fParameters.fQePIDNTracks = qepid.Ntrk();

  return;
}


//------------------------------------------------------------------------------
//skim::ParametersNumu::ParametersNumu(caf::StandardRecord const& sr,
//                                     bool                const& isData)
//: ParametersNumu()
//{
//  
//    //fParameters.fBestPIDTrack         = sr.sel.remid.bestidx;
//    //fParameters.fNeutrinoE            = sr.energy.numusimp.trkccE + (isData)*0.14*(sr.energy.numusimp.trkccE - sr.energy.mutrkE[sr.sel.remid.bestidx].E);
//    //fParameters.fNDHadronicCal        = sr.energy.numusimp.ndhadcalcatE + sr.energy.numusimp.ndhadcaltranE;
//  fParameters.fReMIdVal             = sr.sel.remid.pid;
//  fParameters.fQePIDVal             = sr.sel.qepid.pid;
//  fParameters.fQePIDNTracks         = sr.sel.qepid.ntrk;
//  fParameters.fSliceHits            = sr.slc.nhit;
//  fParameters.fSliceContigPlanes    = sr.slc.ncontplanes;
//  fParameters.fSliceMaxPlane        = sr.slc.lastplane;
//  fParameters.fSliceMinPlane        = sr.slc.firstplane;
//  fParameters.fSliceCellsFromEdge   = sr.slc.ncellsfromedge;
//  fParameters.fPlanesToFront        = sr.sel.contain.nplanestofront;
//  fParameters.fPlanesToBack         = sr.sel.contain.nplanestoback;
//  fParameters.fCosRejKalFwdCell     = sr.sel.contain.kalfwdcell;
//  fParameters.fCosRejKalBakCell     = sr.sel.contain.kalbakcell;
//  fParameters.fCosRejAngleKal       = sr.sel.cosrej.anglekal;
//  fParameters.fCosRejCosBakCell     = sr.sel.contain.cosbakcell;
//  fParameters.fCosRejCosFwdCell     = sr.sel.contain.cosfwdcell;
//  fParameters.fCosRejConCosPID      = sr.sel.cosrej.numucontpid;
//  fParameters.fCosRejKalBakCellND   = sr.sel.contain.kalbakcellnd;
//  fParameters.fCosRejKalFwdCellND   = sr.sel.contain.kalfwdcellnd;
//  fParameters.fCosRejKalYPosAtTrans = sr.sel.contain.kalyposattrans;
//  fParameters.fCosRejKalSpeed       = sr.sel.cosrej.kfitspeed;
//  fParameters.fCosRejCosSpeed       = sr.sel.cosrej.cfitspeed;
//  fParameters.fCosRejCosDirScore    = sr.sel.cosrej.cdirscore;
//    //fParameters.fTrackStartZ          = sr.trk.kalman[sr.sel.remid.bestidx].start.Z();
//    //fParameters.fTrackStopZ           = sr.trk.kalman[sr.sel.remid.bestidx].stop.Z();
//    //fParameters.fTrackStart           = sr.trk.kalman[sr.sel.remid.bestidx].start;
//    //fParameters.fTrackStop            = sr.trk.kalman[sr.sel.remid.bestidx].stop;
//  fParameters.fTrackStartTime       = 0;//Medbh: can't find a way to access this through standard record!! srtrack has no info on its own hits!!
//  fParameters.fTrackStopTime        = 0;
//    //fParameters.fTrackDir             = sr.trk.kalman[sr.sel.remid.bestidx].dir;
//    //fParameters.fTrackCCE             = sr.energy.numusimp.trkccE;
//  
//  art::ServiceHandle<nova::dbi::RunHistoryService> rh;
//  unsigned int dibmask  = (rh->GoodDiBlockMask()&rh->GetConfiguration());
//  std::bitset<14> binary(dibmask);
//
//  std::pair<int,int> planesA = calcFirstLastLivePlane(sr.slc.firstplane, binary);
//  std::pair<int,int> planesB = calcFirstLastLivePlane(sr.slc.lastplane, binary);
//
//  if ((planesA.first != planesB.first) || (planesA.second != planesB.second))fParameters.fDiffLastFirstLiveDiblock = -1;
//  else fParameters.fDiffLastFirstLiveDiblock = (planesA.second - planesA.first + 1) * 0.015625;
//
//  return;
//}

//------------------------------------------------------------------------------
void skim::ParametersNumu::Compare(NumuCutParameters const& compCuts)
{

  std::map<std::string, std::pair<double, double> > comps;
  
  comps["NeutrinoE"]              = std::make_pair(fParameters.fNeutrinoE,            compCuts.fNeutrinoE           );
  comps["NDHadronicCal"]          = std::make_pair(fParameters.fNDHadronicCal,        compCuts.fNDHadronicCal       );
  comps["ReMIdVal"]               = std::make_pair(fParameters.fReMIdVal,             compCuts.fReMIdVal            );
  comps["QePIDVal"]               = std::make_pair(fParameters.fQePIDVal,             compCuts.fQePIDVal            );
  comps["QePIDNTracks"]           = std::make_pair(fParameters.fQePIDNTracks,         compCuts.fQePIDNTracks        );
  comps["SliceHits"]              = std::make_pair(fParameters.fSliceHits,            compCuts.fSliceHits           );
  comps["SliceContigPlanes"]      = std::make_pair(fParameters.fSliceContigPlanes,    compCuts.fSliceContigPlanes   );
  comps["SliceMaxPlane"]          = std::make_pair(fParameters.fSliceMaxPlane,        compCuts.fSliceMaxPlane       );
  comps["SliceMinPlane"]          = std::make_pair(fParameters.fSliceMinPlane,        compCuts.fSliceMinPlane       );
  comps["SliceCellsFromEdge"]     = std::make_pair(fParameters.fSliceCellsFromEdge,   compCuts.fSliceCellsFromEdge  );
  comps["PlanesToFront"]          = std::make_pair(fParameters.fPlanesToFront,        compCuts.fPlanesToFront       );
  comps["PlanesToBack"]           = std::make_pair(fParameters.fPlanesToBack,         compCuts.fPlanesToBack        );
  comps["CosRejKalFwdCell"]       = std::make_pair(fParameters.fCosRejKalFwdCell,     compCuts.fCosRejKalFwdCell    );
  comps["CosRejKalBakCell"]       = std::make_pair(fParameters.fCosRejKalBakCell,     compCuts.fCosRejKalBakCell    );
  comps["CosRejAngleKal"]         = std::make_pair(fParameters.fCosRejAngleKal,       compCuts.fCosRejAngleKal      );
  comps["CosRejCosFwdCell"]       = std::make_pair(fParameters.fCosRejCosFwdCell,     compCuts.fCosRejCosFwdCell    );
  comps["CosRejCosBakCell"]       = std::make_pair(fParameters.fCosRejCosBakCell,     compCuts.fCosRejCosBakCell    );
  comps["CosRejConCosPID"]        = std::make_pair(fParameters.fCosRejConCosPID,      compCuts.fCosRejConCosPID     );
  comps["CosRejKalFwdCellND"]     = std::make_pair(fParameters.fCosRejKalFwdCellND,   compCuts.fCosRejKalFwdCellND  );
  comps["CosRejKalBakCellND"]     = std::make_pair(fParameters.fCosRejKalBakCellND,   compCuts.fCosRejKalBakCellND  );
  comps["CosRejKalYPosAtTrans"] = std::make_pair(fParameters.fCosRejKalYPosAtTrans,   compCuts.fCosRejKalYPosAtTrans);
  comps["TrackStartZ"]          = std::make_pair(fParameters.fTrackStart.Z(),         compCuts.fTrackStart.Z()      );
  comps["TrackStopZ"]             = std::make_pair(fParameters.fTrackStop.Z(),        compCuts.fTrackStop.Z()       );
  comps["TrackCCE"]               = std::make_pair(fParameters.fTrackCCE,             compCuts.fTrackCCE            );

  //Added 2017 ...
  comps["fDistAllTop"]      =  std::make_pair(fParameters.fDistAllTop,      compCuts.fDistAllTop);
  comps["fDistAllBottom"]   =  std::make_pair(fParameters.fDistAllBottom,   compCuts.fDistAllBottom);
  comps["fDistAllWest"]     =  std::make_pair(fParameters.fDistAllWest,     compCuts.fDistAllWest);
  comps["fDistAllEast"]     =  std::make_pair(fParameters.fDistAllEast,     compCuts.fDistAllEast);
  comps["fDistAllFront"]    =  std::make_pair(fParameters.fDistAllFront,    compCuts.fDistAllFront);
  comps["fDistAllBack"]     =  std::make_pair(fParameters.fDistAllBack,     compCuts.fDistAllBack);
  comps["fCVNVal"]          =  std::make_pair(fParameters.fCVNVal,          compCuts.fCVNVal);
  //Added 2018 ...
  comps["fCVN2017Val"]      =  std::make_pair(fParameters.fCVN2017Val,      compCuts.fCVN2017Val);

  comps["fNumShowers"]      =  std::make_pair(fParameters.fNumShowers,      compCuts.fNumShowers);
  comps["fShowerNumXCell"]  =  std::make_pair(fParameters.fShowerNumXCell,  compCuts.fShowerNumXCell);
  comps["fShowerNumYCell"]  =  std::make_pair(fParameters.fShowerNumYCell,  compCuts.fShowerNumYCell);
  comps["fCosShowers"]      =  std::make_pair(fParameters.fCosShowers,      compCuts.fCosShowers);
  comps["fShowerDirX"]      =  std::make_pair(fParameters.fShowerDirX,      compCuts.fShowerDirX);
  comps["fShowerDirY"]      =  std::make_pair(fParameters.fShowerDirY,      compCuts.fShowerDirY);
  comps["fShowerDirZ"]      =  std::make_pair(fParameters.fShowerDirZ,      compCuts.fShowerDirZ);
  comps["fShowerLength"]    =  std::make_pair(fParameters.fProngLength,     compCuts.fProngLength);
  comps["fShowerCalE"]      =  std::make_pair(fParameters.fShowerCalE,      compCuts.fShowerCalE);
  comps["fFracShowerHits"]  =  std::make_pair(fParameters.fFracShowerHits,  compCuts.fFracShowerHits);
  comps["fShowerVertexX"]   =  std::make_pair(fParameters.fShowerVertexX,   compCuts.fShowerVertexX);
  comps["fShowerVertexY"]   =  std::make_pair(fParameters.fShowerVertexY,   compCuts.fShowerVertexY);
  comps["fShowerVertexZ"]   =  std::make_pair(fParameters.fShowerVertexZ,   compCuts.fShowerVertexZ);
  comps["fShowerEndX"]      =  std::make_pair(fParameters.fShowerEndX,      compCuts.fShowerEndX);
  comps["fShowerEndY"]      =  std::make_pair(fParameters.fShowerEndY,      compCuts.fShowerEndY);
  comps["fShowerEndZ"]      =  std::make_pair(fParameters.fShowerEndZ,      compCuts.fShowerEndZ);
  comps["fShowerMinX"]      =  std::make_pair(fParameters.fShowerMinX,      compCuts.fShowerMinX);
  comps["fShowerMinY"]      =  std::make_pair(fParameters.fShowerMinY,      compCuts.fShowerMinY);
  comps["fShowerMaxX"]      =  std::make_pair(fParameters.fShowerMaxX,      compCuts.fShowerMaxX);
  comps["fShowerMaxY"]      =  std::make_pair(fParameters.fShowerMaxY,      compCuts.fShowerMaxY);

  comps["fShowerMinVtxZ"]   =  std::make_pair(fParameters.fShowerMinVtxZ,   compCuts.fShowerMinVtxZ);
  comps["fShowerMinEndZ"]   =  std::make_pair(fParameters.fShowerMinEndZ,   compCuts.fShowerMinEndZ);
  comps["fShowerMaxVtxZ"]   =  std::make_pair(fParameters.fShowerMaxVtxZ,   compCuts.fShowerMaxVtxZ);
  comps["fShowerMaxEndZ"]   =  std::make_pair(fParameters.fShowerMaxEndZ,   compCuts.fShowerMaxEndZ);
  comps["fShowerVtxDist"]   =  std::make_pair(fParameters.fShowerVtxDist,   compCuts.fShowerVtxDist);
  comps["fShowerPt"]        =  std::make_pair(fParameters.fShowerPt,        compCuts.fShowerPt);
  comps["fVertexX"]         =  std::make_pair(fParameters.fVertexX,         compCuts.fVertexX);
  comps["fVertexY"]         =  std::make_pair(fParameters.fVertexY,         compCuts.fVertexY);
  comps["fVertexZ"]         =  std::make_pair(fParameters.fVertexZ,         compCuts.fVertexZ);
  comps["fVertexMaxW"]      =  std::make_pair(fParameters.fVertexMaxW,      compCuts.fVertexMaxW);

  ///--
  comps["fOtherTrackMaxVtxZ"]   =  std::make_pair(fParameters.fOtherTrackMaxVtxZ,   compCuts.fOtherTrackMaxVtxZ);
  comps["fOtherTrackMaxEndZ"]   =  std::make_pair(fParameters.fOtherTrackMaxEndZ,   compCuts.fOtherTrackMaxEndZ);



  double fracDiff = 0.;
  
  for(auto itr : comps){
    
    auto skimSR = itr.second;
    
    fracDiff = std::abs((skimSR.first - skimSR.second)/skimSR.first);
    if( fracDiff > 1.e-3)
      LOG_DEBUG("ParametersNumu")
        << itr.first
        << ": "
        << skimSR.first
        << " "
        << skimSR.second
        << " "
        << fracDiff;
  }
  
  return;
}


//------------------------------------------------------------------------------
skim::ParametersNumu::~ParametersNumu()
{
  
}

namespace skim{
  //------------------------------------------------------------------------------
  std::ostream& operator<< (std::ostream& o, skim::ParametersNumu const& p)
  {
    o
    << "\nBestPIDTrack\t"        << p.ParametersStruct().fBestPIDTrack
    << "\nNeutrinoE\t"           << p.ParametersStruct().fNeutrinoE
    << "\nNDHadronicCal\t"       << p.ParametersStruct().fNDHadronicCal
    << "\nReMIdVal\t"            << p.ParametersStruct().fReMIdVal
    << "\nQePIDVal\t"            << p.ParametersStruct().fQePIDVal
    << "\nQePIDNTracks\t"        << p.ParametersStruct().fQePIDNTracks
    << "\nSliceHits\t"           << p.ParametersStruct().fSliceHits
    << "\nSliceContigPlanes\t"   << p.ParametersStruct().fSliceContigPlanes
    << "\nSliceMaxPlane\t"       << p.ParametersStruct().fSliceMaxPlane
    << "\nSliceMinPlane\t"       << p.ParametersStruct().fSliceMinPlane
    << "\nSliceCellsFromEdge\t"  << p.ParametersStruct().fSliceCellsFromEdge
    << "\nPlanesToFront\t"       << p.ParametersStruct().fPlanesToFront
    << "\nPlanesToBack\t"        << p.ParametersStruct().fPlanesToBack
    << "\nCosRejKalFwdCell\t"    << p.ParametersStruct().fCosRejKalFwdCell
    << "\nCosRejKalBakCell\t"    << p.ParametersStruct().fCosRejKalBakCell
    << "\nCosRejAngleKal\t"      << p.ParametersStruct().fCosRejAngleKal
    << "\nCosRejCosBakCell\t"    << p.ParametersStruct().fCosRejCosBakCell
    << "\nCosRejCosFwdCell\t"    << p.ParametersStruct().fCosRejCosFwdCell
    << "\nCosRejConCosPID\t"     << p.ParametersStruct().fCosRejConCosPID
    << "\nCosRejKalBakCellND\t"  << p.ParametersStruct().fCosRejKalBakCellND
    << "\nCosRejKalFwdCellND\t"  << p.ParametersStruct().fCosRejKalFwdCellND
    << "\nCosRejKalPosAtTrans\t" << p.ParametersStruct().fCosRejKalYPosAtTrans
    //<< "\nTrackStartZ\t"         << p.ParametersStruct().fTrackStartZ
    //<< "\nTrackStopZ\t"          << p.ParametersStruct().fTrackStopZ
    << "\nTrackStart\tx: "   << p.ParametersStruct().fTrackStart.X() << " y: " << p.ParametersStruct().fTrackStart.Y() << " z: " << p.ParametersStruct().fTrackStart.Z()
    << "\nTrackStop\tx: "   << p.ParametersStruct().fTrackStop.X() << " y: " << p.ParametersStruct().fTrackStop.Y() << " z: " << p.ParametersStruct().fTrackStop.Z()
    << "\nTrackDir\tx: "   << p.ParametersStruct().fTrackDir.X() << " y: " << p.ParametersStruct().fTrackDir.Y() << " z: " << p.ParametersStruct().fTrackDir.Z()
    << "\nTrackCCE\t"            << p.ParametersStruct().fTrackCCE  //;
    // Added 2017 
    << "\nfDistAllTop "     << p.ParametersStruct().fDistAllTop
    << "\nfDistAllBottom "  << p.ParametersStruct().fDistAllBottom
    << "\nfDistAllWest "    << p.ParametersStruct().fDistAllWest
    << "\nfDistAllEast "    << p.ParametersStruct().fDistAllEast
    << "\nfDistAllFront "   << p.ParametersStruct().fDistAllFront
    << "\nfDistAllBack "    << p.ParametersStruct().fDistAllBack
    << "\nfCVNVal "         << p.ParametersStruct().fCVNVal
    << "\nfNueCosRejPngPtp "<< p.ParametersStruct().fNueCosRejPngPtp
    // Added 2018
    << "\nfCVN2017Val "         << p.ParametersStruct().fCVN2017Val;
    return o;
  }
}

//-------------------------------------------------------------------------------
//New way of projecting tracks to edge of detector

float skim::ParametersNumu::ProjectDistToEdge(TVector3 Pos, TVector3 Dir)
{

  //  float EdgeTop    =  780;
  //  float EdgeBottom = -780;
  //  float EdgeFront  =  7800;
  //  float EdgeBack   =  0;
  //  float EdgeWest   = -780;
  //  float EdgeEast   =  780;

  float MinDistToEdge = std::numeric_limits<float>::max();
  /*LOG_VERBATIM("ParametersNumu") << "Calculating projections";
  LOG_VERBATIM("ParametersNumu") << "X edges: " << fEdgeXMinusXZ << ", " << fEdgeXPlusXZ;
  LOG_VERBATIM("ParametersNumu")<< "Y edges: " << fEdgeBottomYZ << ", " << fEdgeTopYZ;
  LOG_VERBATIM("ParametersNumu")<< "Z edges: " << fEdgeBackYZ   << ", " << fEdgeFrontYZ;
  LOG_VERBATIM("ParametersNumu")<< "(not " << fEdgeBackXZ   << ", " << fEdgeFrontXZ << ")";
  LOG_VERBATIM("ParametersNumu") << "Dir(X,Y,Z): " << Dir.X() << ", " << Dir.Y() << ", " << Dir.Z() ;
  LOG_VERBATIM("ParametersNumu") << "Pos(X,Y,Z): " << Pos.X() << ", " << Pos.Y() << ", " << Pos.Z() ;*/

  if(Pos.X() > fEdgeXPlusXZ || Pos.X() < fEdgeXMinusXZ ||
     Pos.Y() > fEdgeTopYZ   || Pos.Y() < fEdgeBottomYZ ||
     Pos.Z() > fEdgeBackYZ || Pos.Z() < fEdgeFrontYZ)
    {
      LOG_VERBATIM("ParametersNumu") << "Pos out of range, return -99999";
      return -99999;
    }

  if(Dir.Z() != 0)
    {
      MinDistToEdge = ( fEdgeBackYZ - Pos.Z() ) / Dir.Z();         //distance to back of detector
      LOG_VERBATIM("ParametersNumu") << "Z dist to back = " << MinDistToEdge << "\t or \t";
      if(MinDistToEdge < 0)MinDistToEdge = ( fEdgeFrontYZ - Pos.Z() ) /Dir.Z();   //if that would involve moving backwards, use the distance to the front of the detector
      LOG_VERBATIM("ParametersNumu") << ( fEdgeFrontYZ - Pos.Z() ) /Dir.Z() ;
    }

  if(Dir.Y() != 0)
    {
      float YDist = ( fEdgeTopYZ - Pos.Y() ) / Dir.Y();
      LOG_VERBATIM("ParametersNumu") << "Y dist to top =" << YDist << "\t or \t";
      if(YDist < 0) YDist = ( fEdgeBottomYZ - Pos.Y() ) / Dir.Y();
      LOG_VERBATIM("ParametersNumu") << ( fEdgeBottomYZ - Pos.Y() ) / Dir.Y() ;
      if(YDist < MinDistToEdge) MinDistToEdge = YDist;
    }

  if(Dir.X() != 0)
    {
      float XDist = ( fEdgeXPlusXZ - Pos.X() ) / Dir.X();
      LOG_VERBATIM("ParametersNumu") << "X dist to East = " << XDist << "\t or \t";
      if(XDist < 0) XDist = ( fEdgeXMinusXZ - Pos.X() ) / Dir.X();
      LOG_VERBATIM("ParametersNumu") << ( fEdgeXMinusXZ - Pos.X() ) / Dir.X() ;
      if(XDist < MinDistToEdge) MinDistToEdge = XDist;
    }

  if(MinDistToEdge==99999 || MinDistToEdge < 0)MinDistToEdge = -99999;
  LOG_VERBATIM("ParametersNumu") << "MinDist = " << MinDistToEdge ;
  return MinDistToEdge;

  }

//------------------------------------------------------------------------------
//Need detector edges for ProjectDistToEdge. Shamelessly stolen from livegeometry
void skim::ParametersNumu::GetDetectorEdges()
{
  art::ServiceHandle<geo::Geometry> geom;
  double xyz1[3], xyz2[3];

  if(!(geom->DetId() == novadaq::cnv::kFARDET || geom->DetId() == novadaq::cnv::kNEARDET)) {
    LOG_VERBATIM("ParametersNumu") << "UNKNOWN DETECTOR!!!";
    fEdgeXPlusXZ = 0;
    fEdgeXMinusXZ = 0;
    fEdgeFrontXZ = 0;
    fEdgeBackXZ = 0;
    fEdgeTopYZ = 0;
    fEdgeBottomYZ = 0;
    fEdgeFrontYZ = 0;
    fEdgeBackYZ = 0;
  }

  // start with XZ view, which is easiest since z position does NOT change with cell number; this is odd plane numbers
  geom->Plane(1)->Cell(0)->GetCenter(xyz1);
  fEdgeFrontXZ = xyz1[2]-3.5;  // the half width of the cell in Z is ~3.5 cm
  geom->Plane(63)->Cell(0)->GetCenter(xyz2);
  fEdgeXMinusXZ = (xyz1[0]+xyz2[0])/2.0 - 2.3;  // Plane 0 is 'higher' than plane 63; this repeats for each diblock exactly.  We don't want to worry about this though, so average it.  We pretend things are actual boxes with right angles.
  geom->Plane(1)->Cell(fNumCells)->GetCenter(xyz1);  // cell 383 is most +x, cell 0 is most -X
  geom->Plane(63)->Cell(fNumCells)->GetCenter(xyz2);
  fEdgeXPlusXZ = (xyz1[0]+xyz2[0])/2.0 + 2.3;  // 2.28 is approximate half width of cell here; again, round up slightly
  geom->Plane((fNumDiblocks*64)-1)->Cell(0)->GetCenter(xyz1);  // any cell in last plane will do for finding the back
  fEdgeBackXZ = xyz1[2]+3.5;

  // now YZ view, which has Y AND Z changing with cell number for the same plane.   The Z from cell 0 to cell 383 differs by about 2.5 cm.
  // I'll average the Z values over a DCM to find the DCM edges.  I never want a gap between a DCM edge and the detector edge, since that could unwittingly be considered 'good' even if the closest DCM is 'bad'.  So, I'll use the DCMs that stick out the least to define the edge.  This makes us 'lose' about 2 cm from parts of the edge at each end, but it's safest.

  geom->Plane(0)->Cell(0)->GetCenter(xyz1);  // front bottom most cell
  geom->Plane(0)->Cell(63)->GetCenter(xyz2); // front top cell of frontmost DCM
  fEdgeFrontYZ = (xyz1[2]+xyz2[2])/2.0 - 3.5;  // average the 2 to get the 'front'
  geom->Plane(62)->Cell(0)->GetCenter(xyz2); // back most bottom cell of frontmost bottom DCM
  fEdgeBottomYZ = (xyz1[1]+xyz2[1])/2.0 - 2.3; // reuse same point 1 for bottom calculation
  geom->Plane(0)->Cell(fNumCells)->GetCenter(xyz1); // topmost frontmost cell, in topmost frontmost DCM
  geom->Plane(62)->Cell(fNumCells)->GetCenter(xyz2); // topmost farther back cell in topmost frontmost DCM
  fEdgeTopYZ = (xyz1[1]+xyz2[1])/2.0 + 2.3;  // top edge
  geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells-63)->GetCenter(xyz1);  // frontmost topmost cell in backmost topmost DCM
  geom->Plane((fNumDiblocks*64)-2)->Cell(fNumCells)->GetCenter(xyz2);  // backmost topmost cell, in backmost topmost DCM
  fEdgeBackYZ = (xyz1[2]+xyz2[2])/2.0 + 3.5; // back edge

  // Muon catcher.  TODO: Make these more precise

  // Based on the gdml
  //fMCEdge = 199.65;
  //fMCTop = 75.3;
  //fMCFront = 1277.65;
  //fMCBack = 1577.54;

  // based on the event display
  //fMCEdge = 190;
  //fMCTop = 67.5;
  //fMCFront = 1278.4;
  //fMCBack = 1585;
  }


//------------------------------------------------------------------------------
//Returns the position of the earliest hit in time
std::pair<TVector3,TVector3> skim::ParametersNumu::EarliestLatestHitPos(rb::Track         const& track, bool cosmictrack)
{
  art::ServiceHandle<geo::Geometry> geom;

  float mintime = 1e6;
  float maxtime = 0.001;
  TVector3 mintimepos;//position of earliest hit in time
  TVector3 maxtimepos;//position of latest hit in time

  size_t trackCellsFromEdge = 1e6;
  int closestEdge = -5;

  for(unsigned int i = 0; i < track.NCell(); i++)
    {
      auto chit = track.Cell(i);
      if(chit->TNS() < mintime){
        mintime = chit->TNS();
        double xyz[3];
        geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
        mintimepos.SetXYZ(xyz[0], xyz[1], xyz[2]);
      }
      if (chit->TNS() > maxtime){
        maxtime = chit->TNS();
        double xyz[3];
        geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
        maxtimepos.SetXYZ(xyz[0], xyz[1], xyz[2]);
      }
      
      size_t cellsFromEdge = std::min((unsigned int)chit->Cell(), geom->Plane(chit->Plane())->Ncells() - 1 - chit->Cell());
      
      if( cellsFromEdge < trackCellsFromEdge ){
        trackCellsFromEdge = cellsFromEdge;
        if(chit->View() == geo::kX && cellsFromEdge == chit->Cell())closestEdge = 0;
        else if(chit->View() == geo::kX && cellsFromEdge == geom->Plane(chit->Plane())->Ncells() - 1 - chit->Cell()) closestEdge = 1;
        else if(chit->View() == geo::kY && cellsFromEdge == chit->Cell())closestEdge = 2;
        else if(chit->View() == geo::kY && cellsFromEdge == geom->Plane(chit->Plane())->Ncells() - 1 - chit->Cell()) closestEdge = 3;
        else closestEdge = -5;
      }
      
      if(cosmictrack){
        fParameters.fCosTrackNCellsFromEdge = trackCellsFromEdge;
        fParameters.fCosTrackClosestEdge = closestEdge;
        if(cellsFromEdge < 10) {
          for(unsigned int i = 9; ( i >= cellsFromEdge && i < 10 ); i-=1) fParameters.fCosTrackNHitsNearEdge[i]++;
        }
      } else {
        fParameters.fTrackNCellsFromEdge = trackCellsFromEdge;
        fParameters.fTrackClosestEdge = closestEdge;
        if(cellsFromEdge < 10) {
          for(unsigned int i=9; ( i >= cellsFromEdge && i < 10 ); i-=1) fParameters.fTrackNHitsNearEdge[i]++;
        }
      }
    }

  return std::pair<TVector3, TVector3> (mintimepos, maxtimepos);
}


//-------------------- Analysis 2017 Muon Energy

//  float skim::ParametersNumu::A2017NumuMuE  (float rawEM, float rawHA);
//  float skim::ParametersNumu::A2017NumuHadE (float rawEM, float rawHA);
//  float skim::ParametersNumu::A2017NumuRecoE(float rawEM, float rawHA);