UpMuRecoAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//// Class:       UpMuRecoAna
//// Module Type: analyzer
//// File:        UpMuRecoAna_module.cc
////
//////////////////////////////////////////////////////////////////////////


#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include <art/Framework/Services/Registry/ServiceHandle.h>
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"

#include <iostream>
#include <set>
#include <cmath>
#include <algorithm>
#include <boost/math/special_functions.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include "TNtuple.h"

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "MEFinder/MEClusters.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Vertex.h"
#include "Calibrator/Calibrator.h"
#include "Utilities/func/MathUtil.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include <messagefacility/MessageLogger/MessageLogger.h>

namespace upmuana {
  class UpMuRecoAna;

  inline double getErr(double PE) {
    return 231594.0/(2267.16+pow(PE,2.01011)) + 9.55689;
  }

  inline double getDist(double x1, double y1, double z1, 
                        double x2, double y2, double z2) {
    return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2));
  }

  bool recoHitComp (rb::RecoHit lhs, rb::RecoHit rhs) {
    if (lhs.Y() == rhs.Y()) {
      if (lhs.Z() == rhs.Z()) {
        return lhs.T() < rhs.T();
      }
      else return lhs.Z() < rhs.Z();
    }
    else return lhs.Y() < rhs.Y();
  }
  
  bool pairHitComp (std::pair<rb::CellHit, rb::RecoHit> lhp, 
                    std::pair<rb::CellHit, rb::RecoHit> rhp) {
    return recoHitComp(lhp.second, rhp.second);
  }
  
  inline int getDiblock(int plane) { return (plane/64)+1; } // start at 1
}

class upmuana::UpMuRecoAna : public art::EDAnalyzer {
public:
  explicit UpMuRecoAna(fhicl::ParameterSet const & p);
  virtual ~UpMuRecoAna();

  void analyze(art::Event const & e) override;

  void beginJob() override;

private:
  std::string fHitsLabel;
  std::string fHitsInstance;
  std::string fSliceModuleLabel;
  std::string fSliceInstance;
  std::string fTrackModuleLabel;
  std::string fTrackToSliceInstance;
  std::string fFuzzyKVertexLabel;
  std::string fFuzzyKVertexInstance;
  std::string fMichelELabel;

  double fMinY, fMinX, fMinZ, fMaxY, fMaxX, fMaxZ;

  double      _TrackLen; 
  unsigned    _TrackHitsXY;
  unsigned    _TrackHitsX;
  unsigned    _TrackHitsY; 
  int         _dX; 
  int         _dY; 
  int         _dZ; 
  double      _Chi2;
  double      _LLR; 
  double      _R2X; 
  double      _R2Y; 
  double      _MinSlope;
  double      _MaxSlope;
  bool        fIsSim;
  std::string fParticleModuleLabel;
  std::string fParticleInstance;

  std::string fMCTruthLabel;

  double cw, pw;

  double getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>, 
                          bool(*)(std::pair<rb::CellHit,rb::RecoHit>, 
                                  std::pair<rb::CellHit,rb::RecoHit>)> ,
                std::vector<rb::RecoHit> &outliers,
                double &P_up, double &P_dn, double &chi2, double &slope);
  void LLR(std::vector<double>& eT,
           std::vector<double>& mT,
           std::vector<double>& mTErr, double& slope, double& chi2,
           double& P_up, double& P_dn, std::vector<rb::RecoHit> &in,
           std::vector<rb::RecoHit> &outliers);

  void LinFit(const std::vector<double>& x, 
              const std::vector<double>& y, double *fitpar);
  void LinFit(const std::vector<double>& x, 
              const std::vector<double>& y,
              const std::vector<double>& ye, double *fitpar);
  
 
  struct hit_info {
    float Run, SubRun, Event, TrackID, HitID, X, Y, Z, T,
      deltaT, eT, PE, PECorr, ADC, Plane, Cell, View, GoodTiming, Reco,
      Outlier, dirX, dirY, dirZ, eleAngle, GeV, absT, Diblock;
  };
  
  struct time_window {
    float startT;
    float endT;
  };

  float containmentType(rb::Track);
  float containmentType(sim::Particle);

  TNtuple* fTrackNtuple;
  TNtuple* fHitNtuple;
  TNtuple* fVertexNtuple;
  TNtuple* fEventNtuple;
  TNtuple* fParticleNtuple;
  TNtuple* fSliceNtuple;
  TNtuple* fFLSNtuple;
  unsigned int fNTracks = 0;
};

upmuana::UpMuRecoAna::UpMuRecoAna(fhicl::ParameterSet const & p)
  : EDAnalyzer           (p)
  , fHitsLabel           (p.get<std::string>("HitModuleLabel"  ))
  , fHitsInstance        (p.get<std::string>("HitInstanceLabel"))
  , fSliceModuleLabel    (p.get<std::string>("SliceModuleLabel"))
  , fSliceInstance       (p.get<std::string>("SliceInstanceLabel"))
  , fTrackModuleLabel    (p.get<std::string>("TrackModuleLabel"))
  , fTrackToSliceInstance(p.get<std::string>("TrackToSliceInstance"))
  , fFuzzyKVertexLabel   (p.get<std::string>("FuzzyKVertexLabel"))
  , fFuzzyKVertexInstance(p.get<std::string>("FuzzyKVertexInstance"))
  , fMichelELabel        (p.get<std::string>("MichelELabel"))
  , fMinY                (p.get<double>("MinY"))
  , fMinX                (p.get<double>("MinX"))
  , fMinZ                (p.get<double>("MinZ"))
  , fMaxY                (p.get<double>("MaxY"))
  , fMaxX                (p.get<double>("MaxX"))
  , fMaxZ                (p.get<double>("MaxZ"))
  , _TrackLen            (p.get<double>("TrackLen"))
  , _TrackHitsXY         (p.get<unsigned>("TrackHitsXY"))
  , _TrackHitsX          (p.get<unsigned>("TrackHitsX"))
  , _TrackHitsY          (p.get<unsigned>("TrackHitsY"))
  , _dX                  (p.get<int>("dX"))
  , _dY                  (p.get<int>("dY"))
  , _dZ                  (p.get<int>("dZ"))
  , _Chi2                (p.get<double>("Chi2"))
  , _LLR                 (p.get<double>("LLR"))
  , _R2X                 (p.get<double>("R2X"))
  , _R2Y                 (p.get<double>("R2Y"))
  , _MinSlope            (p.get<double>("MinSlope"))
  , _MaxSlope            (p.get<double>("MaxSlope"))
  , fIsSim               (p.get<bool>("IsSim"))  
  , fParticleModuleLabel (p.get<std::string>("ParticleModuleLabel"))
  , fParticleInstance    (p.get<std::string>("ParticleInstance"))
  , fMCTruthLabel        (p.get<std::string>("MCTruthLabel"))
{
  cw = 3.97;   //cell width in cm
  pw = 6.654;  //plane width in cm
}

upmuana::UpMuRecoAna::~UpMuRecoAna()
{}

void upmuana::UpMuRecoAna::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;

  fTrackNtuple = tfs->make<TNtuple>("track_ntuple", "Track Ntuple",
                                    "Run:SubRun:Event:SliceID:TrackID:Nhits:NRecohits:NOutliers:ProbUp:ProbDn:LLR:Chi2:Slope:LLRX:Chi2X:SlopeX:LLRY:Chi2Y:SlopeY:R2X:R2Y:StartX:StartY:StartZ:StartT:EndX:EndY:EndZ:EndT:TrackHitsX:TrackHitsY:Length:dirX:dirY:dirZ:eleAngle:totalE:containment:avgTX:avgTY:meRecoGeV:meRecoNHits:meRecoX:meRecoY:meRecoZ:meRecoDist:meRecoDeltaT:meRecoAtTrackEnd");
  fHitNtuple  = tfs->make<TNtuple>("hit_ntuple", "Hit Ntuple",
                                   "Run:SubRun:Event:TrackID:HitID:X:Y:Z:T:deltaT:eT:PE:PECorr:ADC:Plane:Cell:View:GoodTiming:Reco:Outlier:dirX:dirY:dirZ:eleAngle:GeV:Diblock");

  fVertexNtuple = tfs->make<TNtuple>("vertex_ntuple", "Vertex Ntuple",
                                     "Run:SubRun:Event:Ntracks:totalHits:totalE:containment:X:Y:Z:T");

  fEventNtuple = tfs->make<TNtuple>("event_ntuple", "Event Ntuple",
                                    "Run:SubRun:Event:Ntracks:NcontainedT:Nvertices:NcontainedV:avgNhits:avgNRecohits:avgLen:hasUpMu:hasFCUpMu:hasTGUpMu:nDiblocks");
  
  fSliceNtuple = tfs->make<TNtuple>("slice_ntuple", "Slice Ntuple",
                                    "Run:SubRun:Event:SliceID:Ntracks:containment");

  if (fIsSim){
    fParticleNtuple = tfs->make<TNtuple>("part_ntuple", "Particle Ntuple",
                                         "Run:SubRun:Event:ParticleID:PDG:Length:StartX:StartY:StartZ:StartT:StartE:EndX:EndY:EndZ:EndT:EndE:StartPX:StartPY:StartPZ:StartEleAngle:containment:MotherID:TrackID:EDep:nFLS");
    /*
    fFLSNtuple      = tfs->make<TNtuple>("fls_ntuple", "FLS Ntuple",
    "Run:SubRun:Event:TrackID:PDG:Cell:Plane:EntryX:EntryY:EntryZ:EntryT:ExitX:ExitY:ExitZ:ExitT:EDep:XAverage:YAverage:ZAverage");*/
  } 
}

void upmuana::UpMuRecoAna::analyze(art::Event const & e)
{
  // getting all hits for backtracker
  
  art::Handle< std::vector<rb::Cluster > > slicecol;
  e.getByLabel( fSliceModuleLabel, slicecol);
  art::PtrVector<rb::Cluster > slces;
  
  for(int i = 0; i < (int)slicecol->size(); i++){
    art::Ptr<rb::Cluster> sli(slicecol, i);
    slces.push_back(sli);
  }
  
  
  // A vector of all hits in the event is needed by backtracker
  
  int nSlices = slces.size();
  
  art::PtrVector< rb::CellHit>  allHits;
  for(int s = 0; s < nSlices; s++){
    for(int c = 0; c<(int)slces[s]->NCell(); c++)
      allHits.push_back( slces[s]->Cell(c) );
  }
  //seems redundent to get all hits twice in two different forms,
  //need to improve this
  

  art::ServiceHandle<cheat::BackTracker>   bt;
  art::Handle<std::vector<rb::CellHit> > hits;
  //art::PtrVector< rb::CellHit>  allHits;
  e.getByLabel(fHitsLabel, fHitsInstance, hits);

  //art::Handle<std::vector<rb::Cluster> > slices;
  //e.getByLabel(fSliceModuleLabel, fSliceInstance, slices);

  art::Handle<std::vector<rb::Track> > tracks;
  e.getByLabel(fTrackModuleLabel, "", tracks);
  art::FindOneP<rb::Cluster> fmSlice(tracks, e,
                                      art::InputTag(fTrackModuleLabel,""));

  //  art::FindOneP<rb::Cluster> fohl(tracks, e, fTrackModuleLabel);
  //  assert(fohl.size() == tracks->size());

  art::Handle<std::vector<rb::Prong> > fuzzyKprongs;
  e.getByLabel(fFuzzyKVertexLabel, fFuzzyKVertexInstance, fuzzyKprongs);
  art::FindManyP<rb::Vertex> fmVert(fuzzyKprongs, e,
                                    art::InputTag(fFuzzyKVertexLabel,fFuzzyKVertexInstance));
  std::vector<art::Ptr<rb::Vertex> > vertices;


  float hasUpMu=0, hasFCUpMu=0, hasTGUpMu=0;

  size_t nContainedT = 0, nVertices = 0, nContainedV = 0;
  float avgNhits = 0, avgNRecohits = 0, avgLen = 0;
  
  for (size_t i_prong=0; i_prong<fuzzyKprongs->size(); ++i_prong) {
    auto theVertexVec = fmVert.at(i_prong);

    for(size_t i_vecvert=0; i_vecvert < theVertexVec.size(); ++i_vecvert) {

      art::Ptr<rb::Vertex> theVertex = theVertexVec.at(i_vecvert);

      bool seen = false;
      // check if this vertex has been seen before
      for (size_t i_vert=0; i_vert<vertices.size(); ++i_vert) {
        if (&(*vertices.at(i_vert)) == &(*theVertex)) {
          seen = true; break;
        }
        else
          vertices.push_back(theVertex);
      }
      if (seen) continue;

      ++nVertices;
      /* float nTracks = theVertex->NTracks(),*/
      float totalHits=0, 
        vertTotalE=0, containment=4; // assume fully contained
      /*
      //std::cout << "Clusters: " << theVertex->NClusters() << std::endl;

      for (size_t i_vprong=0; i_vprong<nTracks; ++i_vprong) {
        rb::Track theProng = *theVertex->Track(i_vprong);
        //if (containmentType(theProng) != 4)
        //  containment = 3; // in-produced
      
        totalHits += theProng.NCell();

        // get total E
        for (size_t i_hit=0; i_hit<theProng.NCell(); ++i_hit) {
          rb::RecoHit theRecoHit = theProng.RecoHit(i_hit);

          if (!theRecoHit.IsCalibrated()) continue;

          vertTotalE += theRecoHit.GeV();
        } // i_hit
      } // i_vprong
      */

      // simplified containment criteria: is the vertex inside the detector?
      if (theVertex->GetX() < fMaxX && theVertex->GetX() > fMinX &&
          theVertex->GetY() < fMaxY && theVertex->GetY() > fMinY &&
          theVertex->GetZ() < fMaxZ && theVertex->GetZ() > fMinZ)
        containment = 4;
      else containment = 1;

      if (containment == 4) ++nContainedV;

      float vertex_entries[11] = {(float)e.id().run(), (float)e.id().subRun(),
                                  (float)e.id().event(), 0/*nTracks*/, totalHits, vertTotalE, containment, 
                                  (float)theVertex->GetX(), (float)theVertex->GetY(), 
                                  (float)theVertex->GetZ(), (float)theVertex->GetT()};
      fVertexNtuple->Fill(vertex_entries);
    } // i_vecvert
  } // i_prong

  unsigned int nGoodTiming = 0;
  for (size_t i_hit=0; i_hit < hits->size(); ++i_hit) {
    if (hits->at(i_hit).GoodTiming()) ++nGoodTiming;
  }

  //std::cout << "Number of hits: " << hits->size() << std::endl;
  //std::cout << "Number of hits with fine timing: " << nGoodTiming << std::endl;
  //std::cout << "Number of tracks: " << tracks->size() << std::endl;
  //std::cout << "Number of vertices: " << nVertices << std::endl;

  std::vector<const rb::Cluster*> slices;
  std::vector<unsigned> tracks_per_slice;
  std::vector<float> slice_containment;

  std::vector<hit_info> allHitInfo;  //ntuple info for all hits
  std::vector<time_window> goodTimeWindows;
 
  // in case we want all hit info
  //time_window full_event_window;
  //full_event_window.startT = -5e6;
  //full_event_window.endT = 5e6;
  //goodTimeWindows.push_back(full_event_window);
  art::FindManyP<me::TrkME> fo(tracks, e, fMichelELabel);
  std::vector<float> activeDiblocks; //keeps track of diblocks with a hit
  for (size_t i_track=0; i_track < tracks->size(); ++i_track) {
    bool keepHits = true;
    rb::Track theTrack = tracks->at(i_track);
    float containment = containmentType(theTrack);
    // find the slice number, if any
    float slice=-1;
    bool found_slice = false;
    const rb::Cluster *theSlice = &(*fmSlice.at(i_track));
    for (size_t i_slice=0; i_slice<slices.size(); ++i_slice) {
      //std::cout << "i_slice " << i_slice << std::endl;
      if (theSlice == slices.at(i_slice))
        { slice = (float)i_slice; found_slice = true; 
          ++tracks_per_slice.at(i_slice);
          if(containment != 4) slice_containment.at(i_slice) = 1;
          break; }
    }
    if (!found_slice) {
      slice = (float)slices.size();
      slices.push_back(theSlice);
      tracks_per_slice.push_back(1);
      if(containment!=4)
        slice_containment.push_back(1);
      else
        slice_containment.push_back(4);
    }
    
    if (!theTrack.Is3D()) continue;

    //std::cout << "3D track id: " << theTrack.ID() << std::endl;

    // get hits in track
    art::PtrVector<rb::CellHit> trackHits = theTrack.AllCells();
    std::set< std::pair<rb::CellHit,rb::RecoHit>, 
              bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
                      std::pair<rb::CellHit,rb::RecoHit>)>
              sortedTrackHits (pairHitComp);
    std::set< std::pair<rb::CellHit,rb::RecoHit>, 
              bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
                      std::pair<rb::CellHit,rb::RecoHit>)>
              sortedTrackHitsX (pairHitComp);
    std::set< std::pair<rb::CellHit,rb::RecoHit>, 
              bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
                      std::pair<rb::CellHit,rb::RecoHit>)>
              sortedTrackHitsY (pairHitComp);
    
    std::vector<double> x_hit; //maybe should change to short? 
    std::vector<double> y_hit;
    std::vector<double> zy_hit;
    std::vector<double> zx_hit;
    //Get Michel Cluster associated with track
    std::vector<art::Ptr<me::TrkME> > michelE  = fo.at(i_track); 
    //initialize to garbage values
    //will change if ME is found
    float  meRecoGeV    = -5;
    float  meRecoHits   = -5;
    float  meRecoPos[3] = {-5, -5, -5};
    float  meRecoDist   = -5;
    float  meRecoDeltaT = -5;
    float  meAtEnd      = 1;
    if(!michelE.empty()){
      meRecoGeV = michelE[0]->TotalGeV();
      meRecoHits = michelE[0]->NCell();
      meRecoPos[0] = michelE[0]->MeanX();
      meRecoPos[1] = michelE[0]->MeanY();
      meRecoPos[2] = michelE[0]->MeanZ();
      meRecoDist = michelE[0]->DistToTrk();
      meRecoDeltaT = michelE[0]->DeltaT();
    } //end Michel if 
    
   
    for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
      art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
      if (!theHit->GoodTiming()) continue;
      rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
      if (!theRecoHit.IsCalibrated()) continue;
      sortedTrackHits.insert(std::make_pair(*theHit,theRecoHit));
      unsigned short iC  = theHit->Cell();
      unsigned short iP  = theHit->Plane();
      geo::View_t   view = theHit->View();
      if(view == geo::kY){
        x_hit.push_back(iC);
        zx_hit.push_back(iP);
        sortedTrackHitsY.insert(std::make_pair(*theHit,theRecoHit));
      }
      else if(view == geo::kX){
        y_hit.push_back(iC);
        zy_hit.push_back(iP);
        sortedTrackHitsX.insert(std::make_pair(*theHit,theRecoHit));
      }
    } // i_hit loop
    
    //    std::cout << "***Sorted RecoHits***" << std::endl;
    //    for (auto i_hit=sortedTrackHits.begin(); 
    //              i_hit != sortedTrackHits.end(); ++i_hit)
    //      std::cout << "Y = " << i_hit->Y() << std::endl;
    //    std::cout << "*********************" << std::endl;
    
    double fitpar[3];
    LinFit(x_hit, zx_hit, fitpar);
    double r2x = fitpar[2];
    LinFit(y_hit, zy_hit, fitpar);
    double r2y = fitpar[2];
    unsigned nxhit = x_hit.size();
    unsigned nyhit = y_hit.size();

    std::vector<rb::RecoHit> outliers, outliersX, outliersY;
    double P_up, P_dn, chi2, slope;
    
    double llr = getLLR(sortedTrackHits, outliers, P_up, P_dn, chi2, 
                        slope);
    
    double chi2X, slopeX, P_upX, P_dnX, chi2Y, slopeY, P_upY, P_dnY,
      llrX = getLLR(sortedTrackHitsX, outliersX, P_upX, P_dnX, 
                    chi2X, slopeX),
      llrY = getLLR(sortedTrackHitsY, outliersY, P_upY, P_dnY,
                    chi2Y, slopeY);
    
    rb::RecoHit start;
    rb::RecoHit end;
    if (sortedTrackHits.size() >= 1) start = (sortedTrackHits.begin()->second);
    else { start = rb::RecoHit(); end = rb::RecoHit(); }
    if (sortedTrackHits.size() >= 2) end = (--sortedTrackHits.end())->second;
    else end = start;
    float dist = theTrack.TotalLength();
    //getDist(start.X(), start.Y(), start.Z(),
    //end.X(),   end.Y(),   end.Z());
    //std::cout << "P_up: " << P_up << ", P_dn: " << P_dn << ", LLR: " << llr
    //          << std::endl;
    
    // averages of all hit directions
    float tDirX=0, tDirY=0, tDirZ=0, tEleAngle, trackTotalE=0;
    
    //Checks to see if hits should be kept
    if(chi2 > _Chi2)                             keepHits = false;
    if(dist < _TrackLen)                         keepHits = false;
    if(TMath::Abs(start.X() - end.X()) < _dX*cw) keepHits = false;
    if(TMath::Abs(start.Y() - end.Y()) < _dY*cw) keepHits = false;
    if(TMath::Abs(start.Z() - end.Z()) < _dZ*pw) keepHits = false;
    if(r2x < _R2X)                               keepHits = false;
    if(r2y < _R2Y)                               keepHits = false;  
    if(nxhit < _TrackHitsX)                      keepHits = false;
    if(nyhit < _TrackHitsY)                      keepHits = false;
    if((nxhit+nyhit) < _TrackHitsXY)             keepHits = false;
    
    if(keepHits){
      time_window window;
      window.startT = start.T() - 2000;
      window.endT   = end.T() + 2000;
      goodTimeWindows.push_back(window);
    }

    float avgTX=0, avgTY=0;

    // now fill hit ntuple
    for (size_t i_hit=0; i_hit < trackHits.size(); ++i_hit) {
      art::Ptr<rb::CellHit> theHit = trackHits.at(i_hit);
      float goodtime = theHit->GoodTiming() ? 1 : 0;
      rb::RecoHit theRecoHit = theTrack.RecoHit(theHit);
      float PE      = theHit->PE(),
            reco    = theRecoHit.IsCalibrated() ? 1 : 0,
            PECorr  = theRecoHit.IsCalibrated() ? theRecoHit.PECorr() : 
        theHit->PE(),
        outlier = 0;
      TVector3 theDir = theTrack.InterpolateDir(theRecoHit.Z());
      float eleAngle = atan(theDir.y()/sqrt(theDir.x()*theDir.x() + 
                                            theDir.z()*theDir.z()));

      tDirX += theDir.x()/trackHits.size();
      tDirY += theDir.y()/trackHits.size();
      tDirZ += theDir.z()/trackHits.size();

      float GeV = 0;

      // only calibrated hits are used for the fit
      if (reco) {
        GeV = theRecoHit.GeV();
        for (size_t i_out=0; i_out < outliers.size(); ++i_out) { 
          if (!outliers.at(i_out).IsCalibrated()) continue;
          if (outliers.at(i_out).PECorr() == theRecoHit.PECorr() &&
              outliers.at(i_out).X() == theRecoHit.X()           &&
              outliers.at(i_out).Y() == theRecoHit.Y()           &&
              outliers.at(i_out).Z() == theRecoHit.Z()           &&
              outliers.at(i_out).T() == theRecoHit.T()) 
            { outlier = 1; break; }
        }
      }

      trackTotalE += GeV;
      if (theHit->View() == geo::kX)
        avgTX += theRecoHit.T()/sortedTrackHitsX.size();
      else if (theHit->View() == geo::kY)
        avgTY += theRecoHit.T()/sortedTrackHitsY.size();

      float deltaT = (float) getErr(PE);
      float eT = (float)(getDist(theRecoHit.X(), theRecoHit.Y(), theRecoHit.Z(),
                                 start.X(), start.Y(), start.Z())/29.97),
            ADC = theHit->ADC();
      //store hit info for later
      hit_info info;
      info.Run = (float)e.id().run();
      info.SubRun = (float)e.id().subRun();
      info.Event = (float)e.id().event();
      info.TrackID = (float)i_track;
      info.HitID = (float)i_hit;
      info.X = (float)theRecoHit.X();
      info.Y = (float)theRecoHit.Y();
      info.Z = (float)theRecoHit.Z();
      info.T = (float)theRecoHit.T() - sortedTrackHits.begin()->second.T();
      info.deltaT = deltaT;
      info.eT =     eT;
      info.PE = PE;
      info.PECorr = PECorr;
      info.ADC = ADC;
      info.Plane = (float)theHit->Plane();
      info.Cell  =  (float)theHit->Cell();
      info.View  = (float)theHit->View();
      info.GoodTiming = goodtime;
      info.Reco = reco;
      info.Outlier = outlier;
      info.dirX = (float)theDir.x();
      info.dirY = (float)theDir.y();
      info.dirZ = (float)theDir.z();
      info.eleAngle = eleAngle;
      info.GeV = GeV;
      info.absT = theRecoHit.T();
      info.Diblock = getDiblock(theHit->Plane());
      allHitInfo.push_back(info);

      if(std::find(activeDiblocks.begin(),activeDiblocks.end(), 
                   info.Diblock)==activeDiblocks.end())
        activeDiblocks.push_back(info.Diblock);
    } // i_hit loop

    
    tEleAngle = atan(tDirY/sqrt(tDirX*tDirX + tDirZ*tDirZ));

    // update event-wide variables
    if (containment == 4) ++nContainedT;
    avgNhits += trackHits.size()/tracks->size();
    avgNRecohits += sortedTrackHits.size()/tracks->size();
    avgLen += dist/tracks->size();

    float track_entries[48]= 
      { (float)e.id().run(), (float)e.id().subRun(),
        (float)e.id().event(), slice, (float)i_track, 
        (float)trackHits.size(), (float)sortedTrackHits.size(), 
        (float)outliers.size(), (float)P_up, (float)P_dn, (float)llr, 
        (float)chi2, (float)slope, (float)llrX, (float)chi2X, (float)slopeX,
        (float)llrY, (float)chi2Y, (float)slopeY, (float)r2x, (float)r2y, 
        start.X(), start.Y(), start.Z(), start.T(), end.X(), end.Y(), 
        end.Z(), end.T(), (float)nxhit,(float)nyhit, dist, tDirX, tDirY, 
        tDirZ, tEleAngle, trackTotalE, containment, avgTX, avgTY,
        meRecoGeV, meRecoHits, meRecoPos[0], meRecoPos[1],
        meRecoPos[2], meRecoDist, meRecoDeltaT, meAtEnd
      };
    fTrackNtuple->Fill(track_entries);

  } // i_track loop

  // fill slice ntuple
  for (size_t i_slice = 0; i_slice < slices.size(); ++i_slice) {
    float slice_entries[6] = { (float)e.id().run(), (float)e.id().subRun(),
                               (float)e.id().event(), (float)i_slice,
                               (float)tracks_per_slice.at(i_slice),
                               (float)slice_containment.at(i_slice)};
    fSliceNtuple->Fill(slice_entries);
  }

  for(size_t i_hitInfo = 0; i_hitInfo < allHitInfo.size(); ++i_hitInfo){
    hit_info info = allHitInfo.at(i_hitInfo);
    bool added = false;
    for(size_t i_window = 0; i_window < goodTimeWindows.size(); i_window++){
      if (added) break;
      time_window window = goodTimeWindows.at(i_window);
      if(info.absT > window.startT && info.absT < window.endT){
        float hit_entries[26] = {info.Run, info.SubRun, info.Event,
                                 info.TrackID, info.HitID, info.X,
                                 info.Y, info.Z, info.T, info.deltaT,
                                 info.eT, info.PE, info.PECorr, info.ADC,
                                 info.Plane, info.Cell,
                                 info.View, info.GoodTiming, info.Reco,
                                 info.Outlier, info.dirX, info.dirY, 
                                 info.dirZ, info.eleAngle, info.GeV, info.Diblock};
        fHitNtuple->Fill(hit_entries);
        added = true;
      } 
    }// i_window
  }// i_hitInfo

 art::Handle<std::vector<sim::Particle> > particles;
 if (fIsSim) {
   e.getByLabel(fParticleModuleLabel, fParticleInstance, particles);
   std::vector<sim::Particle> michelElectrons;
   for (size_t i_part=0; i_part<particles->size(); ++i_part) {
     sim::Particle thePart = particles->at(i_part);
     float startEleAngle = atan(thePart.Py() / sqrt(thePart.Px()*thePart.Px() +
                                                    thePart.Pz()*thePart.Pz())),
       containment = containmentType(thePart);

     if ((thePart.PdgCode()==13 || thePart.PdgCode()==-13) &&
         thePart.Py() > 0 && fabs(startEleAngle)*180/M_PI > 10.0 &&
         thePart.E() > 0.5) {
       hasUpMu = 1;
       if (containment == 4) hasFCUpMu=1;
        if (containment == 1) hasTGUpMu=1;
     }
     
     float length =
       sqrt((thePart.Vx()-thePart.EndX())*(thePart.Vx()-thePart.EndX())
             +(thePart.Vy()-thePart.EndY())*(thePart.Vy()-thePart.EndY())
             +(thePart.Vz()-thePart.EndZ())*(thePart.Vz()-thePart.EndZ()));
     
     //std::vector<sim::FLSHit> flsHits= bt->ParticleToFLSHit(thePart.TrackId(), thePart.PdgCode());
     std::vector<sim::FLSHit> flsHits= bt->ParticleToFLSHit(thePart.TrackId());
     
     //float detLength = -1;
     //float entryT = 9000000;
     //float exitT = 0;
     //sim::FLSHit first;
     //sim::FLSHit last;     
     float eDep = 0;
     for(size_t i_fls = 0; i_fls < flsHits.size(); i_fls++){
       sim::FLSHit h = flsHits[i_fls];
       eDep += h.GetEdep(); 
       /*
       float T = (h.GetEntryT() + h.GetExitT()) / 2; //avg T of flshit
       
       if(T < entryT){
         entryT = T;
         first = h;
       }
       if(T > exitT){
         exitT = T;
         last = h;
       }
       
       float fls_entries[19] = 
         {(float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
          (float)h.GetTrackID(), (float)h.GetPDG(), (float)h.GetCellID(),
          (float)h.GetPlaneID(), (float)h.GetEntryX(), (float)h.GetEntryY(),
          (float)h.GetEntryZ(), (float)h.GetEntryT(), (float)h.GetExitX(),
          (float)h.GetExitY(), (float)h.GetExitZ(), (float)h.GetExitT(),
          (float)h.GetEdep(), (float)h.GetXAverage(), (float)h.GetYAverage(),
          (float)h.GetZAverage()
         };
       */
       //       fFLSNtuple->Fill(fls_entries);
     
     } //end loop on FLS hits
     
     /*
     if(flsHits.size() > 1){ //particles often have no hits
       detLength = sqrt((entryX - exitX)*(entryX - exitX) +
                        (entryY - exitY)*(entryY - exitY) +
                        (entryZ - exitZ)*(entryZ - exitZ));
                        }
     
     if(flsHits.size() > 5){
       UInt_t firstPlane = first.GetPlaneID();
       UInt_t lastPlane = last.GetPlaneID();
       UInt_t firstCell = first.GetCellID();
       UInt_t lastCell = last.GetCellID();
       art::ServiceHandle<geo::Geometry> geom;
       Double_t firstPos[4];
       //double   firstW;
       Double_t lastPos[4];
       //double   lastW;
       const geo::PlaneGeo* firstP = geom->Plane(firstPlane);
       const geo::PlaneGeo* lastP  = geom->Plane(lastPlane);
       const geo::CellGeo*  firstC = firstP->Cell(firstCell);
       const geo::CellGeo*  lastC  = lastP->Cell(lastCell);

       firstC->GetCenter(firstPos);
       lastC->GetCenter(lastPos);
       
       TVector3 traj = {thePart.Vx()-thePart.EndX(),
                        thePart.Vy()-thePart.EndY(),
                        thePart.Vz()-thePart.EndZ()};
       
       double scale = (firstPos[2]-lastPos[2])/traj.Z();
       TVector3 path = traj*scale;
       detLength = path.Mag();              
       }
       }*/
     float particle_entries[25] = 
       { (float)e.id().run(), (float)e.id().subRun(), (float)e.id().event(),
         (float)i_part, (float)thePart.PdgCode(), length,
          (float)thePart.Vx(), (float)thePart.Vy(), (float)thePart.Vz(),
         (float)thePart.T(), (float)thePart.E(), (float)thePart.EndX(),
         (float)thePart.EndY(), (float)thePart.EndZ(),
          (float)thePart.EndT(),(float)thePart.EndE(),
         (float)thePart.Px(), (float)thePart.Py(), (float)thePart.Pz(),
         startEleAngle, containment,(float)thePart.Mother(), 
         (float)thePart.TrackId(), eDep, (float) flsHits.size()};
     fParticleNtuple->Fill(particle_entries);
   }
 }
 
 float event_entries[14] = { (float)e.id().run(), (float)e.id().subRun(),
             (float)e.id().event(), (float)tracks->size(), (float)nContainedT, 
             (float)nVertices, (float)nContainedV, avgNhits, avgNRecohits, 
             avgLen, hasUpMu, hasFCUpMu, hasTGUpMu, 
             (float) activeDiblocks.size() };
 fEventNtuple->Fill(event_entries);


} // analyze() function


double upmuana::UpMuRecoAna::getLLR(std::set< std::pair<rb::CellHit,rb::RecoHit>,
                                   bool(*)(std::pair<rb::CellHit,rb::RecoHit>,
                                           std::pair<rb::CellHit,rb::RecoHit>)> 
                                        sortedHits,
                                    std::vector<rb::RecoHit> &outliers,
                                    double &P_up, double &P_dn, double &chi2,
                                    double &slope)
{
  if (sortedHits.size() < 2) {
    P_up = 1e-30;
    P_dn = 1e-30;
    chi2 = 1e6;
    slope = -1e6;
    return 0;
  }

  std::vector<double> measuredTimes;
  std::vector<double> expectedTimes;
  std::vector<double> wgts;

  // first hit
  measuredTimes.push_back(0.0);
  expectedTimes.push_back(0.0);
  double err = getErr(sortedHits.begin()->first.PE());
  wgts.push_back(1.0/err/err);

  double  startX = sortedHits.begin()->second.X(),
    startY = sortedHits.begin()->second.Y(),
    startZ = sortedHits.begin()->second.Z(),
    startT = sortedHits.begin()->second.T();

  std::vector<rb::RecoHit> in;

  for (auto i_hit =  ++sortedHits.begin();
       i_hit != sortedHits.end();
       ++i_hit) {
    in.push_back(i_hit->second);
    measuredTimes.push_back(i_hit->second.T() - startT);
    double dist = getDist(i_hit->second.X(), i_hit->second.Y(), 
                          i_hit->second.Z(), startX, startY, startZ);
    expectedTimes.push_back(dist/29.97);
    err = getErr(i_hit->first.PE());
    wgts.push_back(1.0/err/err);
  }

  LLR(expectedTimes, measuredTimes, wgts, slope, chi2, P_up, P_dn, 
      in, outliers);

  return log(P_up/P_dn);
}

void upmuana::UpMuRecoAna::LLR(std::vector<double>& eT, 
                               std::vector<double>& mT,
                               std::vector<double>& wgts, double& slope, 
                               double& chi2, double& P_up, double& P_dn, 
                               std::vector<rb::RecoHit>& in,
                               std::vector<rb::RecoHit>& outliers)
{
  // eT - x, mT - y
  size_t n = mT.size();

  // y = a + bx
  double a, b;
  if (eT.size() >= 2)
    chi2 = util::LinFit(eT, mT, wgts, b, a);
  else {
    chi2 = 1e6;
    slope = -1e6;
    P_up = 1e-30;
    P_dn = 1e-30;
    return;
  }
  //a = fitpar[0];
  //b = fitpar[1];

  size_t totAllowOutlier = n/10;
  size_t nOutliers = 0;
  std::vector<double> x_filt, y_filt, w_filt;
  for (size_t i=0; i < n; i++) {
    double y_fit = a + b*eT.at(i);
    if ( fabs(mT.at(i) - y_fit) < 5*(1.0/sqrt(wgts.at(i))) 
         || nOutliers>totAllowOutlier)
      {
        x_filt.push_back(eT.at(i));
        y_filt.push_back(mT.at(i));
        w_filt.push_back(wgts.at(i));
      }
    else {
      ++nOutliers;
      outliers.push_back(in[i]);
    }
  }

  if (x_filt.size() >= 2)
    chi2 = util::LinFit(x_filt, y_filt, w_filt, b, a);
  else {
    chi2 = 1e6;
    slope = -1e6;
    P_up = 1e-30;
    P_dn = 1e-30;
  }
  //a = fitpar[0];
  //b = fitpar[1];
  n = x_filt.size();
  if (n < 5) {
    slope = 0;
    chi2  = 999;
    P_up  = 1e-30;
    P_dn  = 1;
    return;
  }

  slope = b;
  //chi2 = 0;
  //for (size_t i=0; i<n; i++) {
  //  double y_exp = a + b*x_filt.at(i);
  //  chi2 += pow((y_filt.at(i)-y_exp) / ye_filt.at(i), 2.0);
  //}
  chi2 /= (double)(n-2);

  double  one_over_ee = 0.0,
    x_over_ee   = 0.0,
    y_over_ee   = 0.0;

  for (size_t i=0; i<n; ++i) {
    //double     e = ye_filt.at(i);
    one_over_ee += w_filt.at(i);
    x_over_ee   += x_filt.at(i)*w_filt.at(i);
    y_over_ee   += y_filt.at(i)*w_filt.at(i);
  }

  double up_inter = (y_over_ee-x_over_ee)/one_over_ee;
  double dn_inter = (y_over_ee+x_over_ee)/one_over_ee;

  double up_chi2 = 0.0, dn_chi2 = 0.0;
  for (size_t i=0; i<n; ++i) {
    double e = 1.0/sqrt(w_filt.at(i));
    double up_expected = up_inter + x_filt.at(i);
    double dn_expected = dn_inter - x_filt.at(i);
    up_chi2 += pow((y_filt.at(i)-up_expected) / e, 2.0);
    dn_chi2 += pow((y_filt.at(i)-dn_expected) / e, 2.0);
  }

  double prob_up = boost::math::gamma_q((double)(n-2)/2.0,up_chi2/2.0),
    prob_dn = boost::math::gamma_q((double)(n-2)/2.0,dn_chi2/2.0);

  if (prob_up < 1e-30) prob_up = 1e-30;
  if (prob_dn < 1e-30) prob_dn = 1e-30;

  P_up = prob_up;
  P_dn = prob_dn;

}

/*void upmuana::UpMuRecoAna::LinFit(std::vector<double>& x,
  std::vector<double>& y,
  std::vector<double>& y_err,
  double *fitpar) {
  const size_t n = x.size();
  double xsum=0, ysum=0, xxsum=0, yysum=0, xysum=0, eesum=0;

  for (size_t i=0; i<n; ++i) {
  xsum  += x[i]/y_err[i]/y_err[i];
  ysum  += y[i]/y_err[i]/y_err[i];
  xxsum += x[i]*x[i]/y_err[i]/y_err[i];
  yysum += y[i]*y[i]/y_err[i]/y_err[i];
  xysum += x[i]*y[i]/y_err[i]/y_err[i];
  eesum += 1./y_err[i]/y_err[i];
  }

  const double b = (xysum*eesum - xsum*ysum)/(xxsum*eesum - xsum*xsum);
  const double a = (xysum    - b*xxsum)/xsum;
  const double r = (xysum - xsum*ysum)/sqrt(xxsum - xsum*xsum)
  *(yysum - ysum*ysum);

  fitpar[0] = a;
  fitpar[1] = b;
  fitpar[2] = r*r;
  }*/

float upmuana::UpMuRecoAna::containmentType(sim::Particle thePart)
{
  TLorentzVector start = TLorentzVector(thePart.Position());
  TLorentzVector stop  = TLorentzVector(thePart.EndPosition());

  if (start.Y() > stop.Y()) // assume upward-going
  { TLorentzVector hold = start; start = stop; stop = hold; }

  if (start.Y() < fMinY || start.X() < fMinX || start.X() > fMaxX ||
      start.Z() < fMinZ || start.Z() > fMaxZ) { // entered from outside
    if (stop.Y() > fMaxY || 
        stop.X() < fMinX || 
        stop.X() > fMaxX ||
        stop.Z() < fMinZ ||
        stop.Z() > fMaxZ) { return 1; }//through-going
    return 2; // stopped
  } // started in detector
  if (stop.Y() > fMaxY || stop.X() < fMinX || stop.X() > fMaxX ||
      stop.Z() < fMinZ || stop.Z() > fMaxZ)
    { return 3; }// in-produced

  return 4; // fully contained
  
}

float upmuana::UpMuRecoAna::containmentType(rb::Track theTrack)
{
  TVector3 start = theTrack.Start();
  TVector3 stop  = theTrack.Stop();

  //std::cout << "Start: (" << start.X() <<", " << start.Y() << ", " << start.Z() << ")" << std::endl;
  //std::cout << "Stop:  (" << stop.X() << ", " << stop.Y() << ", " << stop.Z() << ")" << std::endl;

  if (start.y() > stop.y()) // assume upward-going, so switch
    { TVector3 hold = start; start = stop; stop = hold; }

  //std::cout << "Start: (" << start.X() <<", " << start.Y() << ", " << start.Z() << ")" << std::endl;
  //std::cout << "Stop:  (" << stop.X() << ", " << stop.Y() << ", " << stop.Z() << ")" << std::endl;

  if (start.y() < fMinY || start.x() < fMinX || start.x() > fMaxX ||
      start.z() < fMinZ || start.z() > fMaxZ) { // entered from outside
    if (stop.y() > fMaxY || 
        stop.x() < fMinX || 
        stop.x() > fMaxX ||
        stop.z() < fMinZ ||
        stop.z() > fMaxZ) { return 1; }//through-going
    return 2; // stopped
  } // started in detector
  if (stop.y() > fMaxY || stop.x() < fMinX || stop.x() > fMaxX ||
      stop.z() < fMinZ || stop.z() > fMaxZ)
    { return 3; }// in-produced

  return 4; // fully contained
  
}


void upmuana::UpMuRecoAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, double *fitpar){
  // http://en.wikipedia.org/wiki/Simple_linear_regression
  const auto n    = x_val.size();
  const auto x  = (std::accumulate(x_val.begin(), x_val.end(), 0.0))/n;                       // <x>
  const auto y  = (std::accumulate(y_val.begin(), y_val.end(), 0.0))/n;                       // <y>
  const auto xx = (std::inner_product(x_val.begin(), x_val.end(), x_val.begin(), 0.0))/n;     // <xx>
  const auto yy = (std::inner_product(y_val.begin(), y_val.end(), y_val.begin(), 0.0))/n;     // <yy>
  const auto xy = (std::inner_product(x_val.begin(), x_val.end(), y_val.begin(), 0.0))/n;     // <xy>

  const auto b    = (xy - x*y)/(xx - x*x);                                                    // slope
  const auto a    = y - b*x;                                                                  // intercept
  const auto r    = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y));                                   // Rxy - coeffcient of determination
  fitpar[0] = a;
  fitpar[1] = b;
  fitpar[2] = r*r;

}

void upmuana::UpMuRecoAna::LinFit(const std::vector<double>& x_val, const std::vector<double>& y_val, const std::vector<double>& y_err, double *fitpar){
  // http://en.wikipedia.org/wiki/Simple_linear_regression
  int n    = x_val.size();
  double x = 0; 
  double y = 0; 
  double xx = 0;
  double yy = 0;
  double xy = 0;
  double ee = 0;

  for ( int i=0; i<n; i++ ){

    x = x + x_val[i]/y_err[i]/y_err[i];
    y = y + y_val[i]/y_err[i]/y_err[i];

    xx = xx + x_val[i]*x_val[i]/y_err[i]/y_err[i];
    yy = yy + y_val[i]*y_val[i]/y_err[i]/y_err[i];
    xy = xy + x_val[i]*y_val[i]/y_err[i]/y_err[i];
    ee = ee + 1./y_err[i]/y_err[i];
  }

  const auto b    = (xy*ee - x*y)/(xx*ee - x*x);                                              // slope
  const auto a    = (xy - b*xx)/x;                                                            // intercept
  const auto r    = (xy - x*y)/sqrt((xx - x*x)*(yy - y*y));                                   // Rxy - coeffcient of determination
  fitpar[0] = a;
  fitpar[1] = b;
  fitpar[2] = r*r;
  
  
}

DEFINE_ART_MODULE(upmuana::UpMuRecoAna)