T0RecoAnodePiercers_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       T0RecoAnodePiercers
// Module Type: producer
// File:        T0RecoAnodePiercers.cc
//
// Joshua Thompson - joshualt@fnal.gov
// developed from work by Hannah Rogers   - hannah.rogers@colostate.edu
// based on uboonecode modules by David Caratelli and Chris Barnes
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.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 "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art_root_io/TFileService.h"

// services etc...
#include "larcore/Geometry/Geometry.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"

// data-products
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/OpFlash.h"
#include "lardataobj/AnalysisBase/T0.h"
#include "lardata/Utilities/AssociationUtil.h"
//#include "duneprototypes/Protodune/singlephase/DataUtils/ProtoDUNEPFParticleUtils.h"
//#include "duneprototypes/Protodune/singlephase/DataUtils/ProtoDUNETrackUtils.h"

// ROOT
#include "TVector3.h"

// C++
#include <memory>
#include <iostream>
#include <utility>

class T0RecoAnodePiercers;

class T0RecoAnodePiercers : public art::EDProducer {
        public:
        explicit T0RecoAnodePiercers(fhicl::ParameterSet const & p);
        // The destructor generated by the compiler is fine for classes
        // without bare pointers or other resource use.

        // Plugins should not be copied or assigned.
        T0RecoAnodePiercers(T0RecoAnodePiercers const &) = delete;
        T0RecoAnodePiercers(T0RecoAnodePiercers &&) = delete;
        T0RecoAnodePiercers & operator = (T0RecoAnodePiercers const &) = delete;
        T0RecoAnodePiercers & operator = (T0RecoAnodePiercers &&) = delete;

        //Required functions.

        void produce(art::Event& event);

        private:

        // Functions to be used throughout module

        void   SortTrackPoints  (const recob::Track& track, std::vector<TVector3>& sorted_trk);

        size_t FlashMatch               (const double reco_time, std::vector<double> op_times_v);

        // Declare member data here

        std::string             fPFPProducer;
        std::string     fTrackProducer;
        std::string     fHitProducer;
        std::string             fTriggerProducer;

        std::string     fFlashProducer;
        std::string     fFlashProducerMC;
        std::string     fFlashProducerData;

        double          fEdgeWidth; // [cm]
                int                     fReadoutEdgeTicks;
        
        int                     fMinPE;
        double          fMinTrackLength;

        double          fFlashScaleFactor;
        double          fFlashTPCOffset;

        bool            fDebug;
        
        double          fMinDtData;
        double          fMaxDtData;

        double          fMinDtMC;
        double          fMaxDtMC;

        double          DriftVelocity; // [cm/us]
        unsigned int    ReadoutWindow;

        double          det_top;
        double          det_bottom;
        double          det_front;
        double          det_back;
        double          det_width; // [cm]
        
        std::vector<double>     op_times;

        std::vector<size_t>     flash_id_v;
        art::Handle<std::vector<recob::OpFlash> >       flash_h;

        bool            MC;

        // Track parameters
        double          anode_rc_time;
        double          length;
        double          rc_xs, rc_xe;
                //double                rc_xs_corr, rc_xe_corr;
        double          rc_ys, rc_ye;
        double          rc_zs, rc_ze;

        // Flash parameters
        double          matched_flash_time;
        double          corrected_matched_flash_time;
        //double                matched_flash_time_width;

        double          matched_flash_pe;       
        //double                matched_flash_centre_y;
        //double                matched_flash_centre_z;
        //double                matched_flash_max_pe_det_x;
        //double                matched_flash_max_pe_det_y;
        //double                matched_flash_max_pe_det_z;
        //double                matched_flash_width_y;
        //double                matched_flash_width_z;
        
        // Reco results
        double          dt_flash_reco;

        // Track info
        bool            readout_edge;
        bool            TPC_entering_candidate;
        bool            TPC_exiting_candidate;
        bool            anode_piercing_candidate;
        //bool          cathode_crossing_track;
        
                double          dtMin;
        double          dtMax;
        };

T0RecoAnodePiercers::T0RecoAnodePiercers(fhicl::ParameterSet const & fcl)
        :
    EDProducer(fcl) {

        fDebug                          = fcl.get<bool>         ("Debug"        );

        fPFPProducer            = fcl.get<std::string>  ("PFPProducer"                  );
        fTrackProducer          = fcl.get<std::string>  ("TrackProducer"        );
        fHitProducer            = fcl.get<std::string>  ("HitProducer"          );

        fTriggerProducer        = fcl.get<std::string>  ("TriggerProducer"      );

        fFlashProducerData      = fcl.get<std::string>  ("DataFlashProducer"    );
        fFlashProducerMC        = fcl.get<std::string>  ("MCFlashProducer"      );

        fEdgeWidth                      = fcl.get<double>       ("EdgeWidth"                    );
        fReadoutEdgeTicks       = fcl.get<int>          ("ReadoutEdgeTicks"             );

        fMinPE                          = fcl.get<double>       ("MinPE"                        );
        fMinTrackLength         = fcl.get<double>       ("MinTrackLength"       );

        fMinDtData                      = fcl.get<double>       ("MinDtFlashRecoData"   );
        fMaxDtData                      = fcl.get<double>       ("MaxDtFlashRecoData"   );

        fMinDtMC                        = fcl.get<double>       ("MinDtFlashRecoMC"     );
        fMaxDtMC                        = fcl.get<double>       ("MaxDtFlashRecoMC"     );

        fFlashScaleFactor       = fcl.get<double>       ("FlashScaleFactor"     );
        fFlashTPCOffset         = fcl.get<double>       ("FlashTPCOffset"       );

        // get boundaries based on detector bounds
        auto const* geom = lar::providerFrom<geo::Geometry>();

        det_top = fEdgeWidth;
        det_bottom = fEdgeWidth;
        det_front = fEdgeWidth;
        det_back = fEdgeWidth;

        for (geo::TPCID const& tID: geom->IterateTPCIDs()) {
                geo::TPCGeo const& TPC = geom->TPC(tID);

                if(TPC.DriftDistance() < 25.0) continue;

                double origin[3] = {0.};
                double center[3] = {0.};
                TPC.LocalToWorld(origin, center);

                double tpc_top = center[1] + TPC.HalfHeight();
                double tpc_bottom = center[1] - TPC.HalfHeight();
                double tpc_front = center[2] - TPC.HalfLength();
                double tpc_back = center[2] + TPC.HalfLength();

                if (tpc_top     > det_top)              det_top = tpc_top;
                if (tpc_bottom  < det_bottom)   det_bottom = tpc_bottom;
                if (tpc_front   < det_front)    det_front = tpc_front;
                if (tpc_back    > det_back)     det_back = tpc_back;  

                det_width = TPC.DriftDistance();
                }

        // Use 'detp' to find 'efield' and 'temp'
        auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob();
        double efield = detProp.Efield();
        if(fDebug) std::cout << "Nominal electric field is: " << efield << " V/cm" << std::endl;

        double temp   = detProp.Temperature();
        if(fDebug) std::cout << "LAr temperature is: " << temp << " K" << std::endl;

        // Determine the drift velocity from 'efield' and 'temp'
        DriftVelocity = detProp.DriftVelocity(efield,temp);
        if(fDebug) std::cout << "Drift velocity is: " << DriftVelocity << " cm/us" << std::endl;

        // Get Readout window length

        ReadoutWindow = detProp.ReadOutWindowSize();
        if(fDebug) std::cout << "Readout window is: " << ReadoutWindow << " us" << std::endl;

        // Declare what the module produces

        produces < std::vector<anab::T0> >();
        produces < art::Assns<recob::PFParticle, anab::T0> >();
        }

void T0RecoAnodePiercers::produce(art::Event& event){

        MC = !(event.isRealData());

        if(MC) {
                dtMin = fMinDtMC; //us
                dtMax = fMaxDtMC; //us
                }
        else {
                dtMin = fMinDtData; //us
                dtMax = fMaxDtData; //us
                }
 
        if(fDebug) std::cout << "\tAnode piercing tracks selected by cuts.\n\t\t Min length: " 
                << fMinTrackLength << ", Min PE: " << fMinPE << 
                " and flash-reco time difference between " << dtMin << " and " 
                << dtMax << " us." << std::endl;

        auto t0_v    = std::make_unique<std::vector<anab::T0>>();
        auto pfp_t0  = std::make_unique<art::Assns<recob::PFParticle, anab::T0>>();

        int track_number = 0;

        // Load detector clocks for later
        auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);;

        double TPC_trigger_offset = 0.0;

        //std::cout << "Event number: " << event.event() << std::endl;
        if(fDebug) std::cout << "Set event number: " << event.event() << "\nTop: " 
                << det_top << "\nBottom: " << det_bottom << "\nFront: " << det_front 
                << "\nBack: " << det_back << "\nEdge width: " << fEdgeWidth << std::endl;  

        op_times.clear();
        flash_id_v.clear();
        flash_h.clear();

        //Set flash producer to MC or data

        if(!MC) fFlashProducer = fFlashProducerData;
        if(MC) fFlashProducer = fFlashProducerMC;

        // load Flashes
        if (fDebug) std::cout << "Loading flash from producer " << fFlashProducer << std::endl;

        flash_h = event.getHandle<std::vector<recob::OpFlash> >(fFlashProducer);

        // make sure flashes look good
        if(!flash_h.isValid()) {
        std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate Flash!"<<std::endl;
        throw std::exception();
                }

        double trigger_time = 0;

        if(!MC){
                auto trigger_h = event.getHandle<std::vector<recob::OpFlash> >(fTriggerProducer);

                if( (!trigger_h) || trigger_h->empty()) {
                if(fDebug) std::cout << "\tTrigger not found. Skipping." << std::endl;
                        event.put(std::move(t0_v));
                        event.put(std::move(pfp_t0));
                return;
                }

                if(fDebug) std::cout << "Loading trigger time from producer " 
                        << fTriggerProducer << std::endl;

                trigger_time = trigger_h->at(0).Time();
                }


        // load PFParticles

        auto reco_particles_h = event.getValidHandle<std::vector<recob::PFParticle>>(fPFPProducer);
        const art::FindManyP<recob::Track> findTracks(reco_particles_h,event,fTrackProducer);
        auto reco_tracks_h = event.getValidHandle<std::vector<recob::Track>>(fTrackProducer);
        art::FindManyP<recob::Hit> findHits(reco_tracks_h, event, fTrackProducer);

        if(!reco_particles_h.isValid()) {
                std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate PFParticles!"<<std::endl;
                throw std::exception(); 
                }

        // Utilities for PFParticles and tracks
        //protoana::ProtoDUNEPFParticleUtils pfpUtil;
        //protoana::ProtoDUNETrackUtils trackUtil;

        // Get trigger to TPC Offset

        TPC_trigger_offset = clockData.TriggerOffsetTPC();
        if(fDebug) std::cout << "TPC time offset from trigger: " 
                << TPC_trigger_offset << " us" << std::endl;

        // Prepare a vector of optical flash times, if flash above some PE cut value

        size_t flash_ctr = 0;
        for (auto const& flash : *flash_h){
                if (flash.TotalPE() > fMinPE){
                        double op_flash_time;
                        if(!MC) op_flash_time = flash.Time() - trigger_time;
                        if(MC) op_flash_time = flash.Time() - trigger_time - TPC_trigger_offset;
                op_times.push_back(op_flash_time);
                        flash_id_v.push_back(flash_ctr);
                        if (fDebug) std::cout << "\t Flash: " << flash_ctr << " has time : " 
                        << op_flash_time << ", PE : " << flash.TotalPE() << std::endl;
                        }
                flash_ctr++;
                } // for all flashes

        if(fDebug) std::cout << "Selected a total of " << op_times.size() << " OpFlashes" << std::endl;

        // LOOP THROUGH RECONSTRUCTED PFPARTICLES

        size_t ev_particle_ctr = 0;

        for(unsigned int particle = 0; particle < reco_particles_h->size(); ++particle){

                const recob::PFParticle &pfparticle = (*reco_particles_h)[particle];

        // Only consider primary particles
        if(!pfparticle.IsPrimary()) continue;

                ev_particle_ctr++;


                //Replacing this with the hardcoded method to remove util dependency
                //const recob::Track* track = pfpUtil.GetPFParticleTrack(pfparticle,event,fPFPProducer,fTrackProducer);
                const recob::Track* track = 0x0;
                const std::vector<art::Ptr<recob::Track>> pfpTracks = findTracks.at(pfparticle.Self());
                if( pfpTracks.size() != 0 ){
                  track = (pfpTracks.at(0)).get(); 
                }
                if(track == 0x0) { 
                        if(fDebug) std::cout << "\tPFParticle " << ev_particle_ctr << " is not track like" << std::endl;
                        continue; 
                } 

                if (fDebug) std::cout << "\tLooping through reco PFParticle " << ev_particle_ctr << std::endl;  

                //Replacing this with the hardcoded method to remove util dependency
                //const std::vector<const recob::Hit*>& hit_v = trackUtil.GetRecoTrackHits(*track,event,fTrackProducer);
                std::vector<art::Ptr<recob::Hit>> inputHits = findHits.at(track->ID());
                std::vector<const recob::Hit*> hit_v;
                for(const art::Ptr<recob::Hit> hit : inputHits){
                  hit_v.push_back(hit.get());
                }


                anode_rc_time = ReadoutWindow;

                length = 0.;
                rc_xs = ReadoutWindow/10.;
                rc_xe = -ReadoutWindow/10.; 
                //rc_xs_corr = 399.;
                //rc_xe_corr = -399.;
                rc_ys = -99.;
                rc_ye = -99.; 
                rc_zs = -99.; 
                rc_ze = -99.;

                matched_flash_time = -ReadoutWindow;
                corrected_matched_flash_time = -ReadoutWindow;
                //matched_flash_time_width = -1.;

                matched_flash_pe = 0.;
                //matched_flash_centre_y = -99.;
                //matched_flash_centre_z = -99.;
                //matched_flash_width_y = -99.;
                //matched_flash_width_z = -99.;
                //matched_flash_max_pe_det_x = 0.;
                //matched_flash_max_pe_det_y = -99.;
                //matched_flash_max_pe_det_z = -99.;

                dt_flash_reco = 999;

                anode_piercing_candidate = false;
                //cathode_crossing_track = false;

                // Get sorted points for the track object [assuming downwards going]

                std::vector<TVector3> sorted_trk;
                SortTrackPoints(*track,sorted_trk);

                TVector3 track_start = sorted_trk.at(0);
                TVector3 track_end = sorted_trk.at(sorted_trk.size() - 1);

                if(fDebug) std::cout << "\t\tTrack goes from (" << track_start.X() << ", " 
                        << track_start.Y() << ", " << track_start.Z() << ") --> (" << track_end.X() << ", " 
                        << track_end.Y() << ", " << track_end.Z() << ")" << std::endl;

                if(track->Length() < fMinTrackLength )  {
                                if(fDebug) std::cout << "\t\t\tParticle track too short. Skipping." << std::endl;
                                continue;
                        }

                auto const* geom = lar::providerFrom<geo::Geometry>();   
                auto const* first_hit = hit_v.at(0);
                const geo::WireID wireID_start = first_hit->WireID();
                const auto TPCGeoObject_start = geom->TPC(wireID_start.TPC,wireID_start.Cryostat);
                short int driftDir_start = TPCGeoObject_start.DetectDriftDirection();

                short int driftDir_end = 0;
                //cathode_crossing_track = false;

            for (size_t ii = 1; ii < hit_v.size(); ii++) {
                const geo::WireID wireID_end = hit_v.at(ii)->WireID();
                        const auto TPCGeoObject_end = geom->TPC(wireID_end.TPC,wireID_end.Cryostat);
                        driftDir_end = TPCGeoObject_end.DetectDriftDirection(); 
                
                        if(driftDir_end + driftDir_start == 0){
                                //cathode_crossing_track = true;
                                ii = hit_v.size();
                                }
                        }
                // ------------------------------------------------------------------------------------
                // ANODE PIERCERS 

                if(fDebug) std::cout << "\t\tThis track starts in TPC " << wireID_start.TPC 
                << " which has a drift direction of " << driftDir_start << std::endl; 

                // create root trees variables

                rc_xs = track_start.X();
                rc_ys = track_start.Y();
                rc_zs = track_start.Z();
                rc_xe = track_end.X();
                rc_ye = track_end.Y();
                rc_ze = track_end.Z();
                length = track->Length();

                // Determine if track hits edge of readout window
                // NECESSARY TO REMOVE TRACKS THAT AREN'T FINISHED WHEN WINDOW CLOSES
                // AS WELL AS TRACKS THAT WERE BEING COLLECTED BEFORE THE WINDOW OPENED

                readout_edge = false;
                for (auto& hits : hit_v) {
                        auto hit_tick = hits->PeakTime();
                        if(hit_tick < fReadoutEdgeTicks || hit_tick > (ReadoutWindow - fReadoutEdgeTicks)){
                                readout_edge = true;
                                if(fDebug) std::cout << "\tTrack hits edge of readout window. "
                                        "Skipping." << std::endl;
                                continue;
                                }

                        double hit_time = clockData.TPCTick2TrigTime(hit_tick);

                        // If track within window, get reco time from earliest hit time
                        if (hit_time < anode_rc_time) anode_rc_time = hit_time;
                        }

                if(readout_edge) continue;

                TPC_entering_candidate = false;
                TPC_exiting_candidate = false;

                // Tracks which may enter TPC through an APA
                if (rc_ys < (det_top - fEdgeWidth) && rc_zs > (det_front + fEdgeWidth) 
                        && rc_zs < (det_back - fEdgeWidth)) {

                        // reconstruct track T0 w.r.t. trigger time
                        if( ( rc_xs > rc_xe && driftDir_start>0 ) || 
                                ( rc_xs < rc_xe && driftDir_start<0 ) ) {
                                TPC_entering_candidate = true;
                                if(fDebug) std::cout << "\t\tTrack may enter TPC through "
                                "anode. Reco t0: " << anode_rc_time << " us" << std::endl;
                                }
                        }

                // Tracks which may exit TPC through an APA
                if (rc_ye > (det_bottom + fEdgeWidth) && rc_ze > (det_front + fEdgeWidth) 
                        && rc_ze < (det_back - fEdgeWidth)) {

                        // reconstruct track T0 w.r.t. trigger time     
                        if( ( rc_xe > rc_xs && driftDir_end>0 ) || 
                        ( rc_xe < rc_xs && driftDir_end<0 ) ) { 
                                TPC_exiting_candidate = true;
                                if(fDebug) std::cout << "\t\tTrack may exit TPC through "
                                        "anode. Reco t0:" << anode_rc_time << " us" << std::endl;
                                }
                        }

                if(TPC_entering_candidate||TPC_exiting_candidate) 
                        anode_piercing_candidate = true;

                if(!anode_piercing_candidate) {
                        if(fDebug) std::cout << "\t\tTrack does not pierce anode." << std::endl;
                        continue; 
                        }

                if(TPC_entering_candidate&&TPC_exiting_candidate) {
                        if(fDebug) std::cout << "\t\tTrack neither enters nor exits"
                                " through non-anode TPC face. No useful end point for SCE" 
                                " measurement." << std::endl;
                        continue; 
                }

                //rc_xs_corr = rc_xs + driftDir_start*anode_rc_time*DriftVelocity;
                //rc_xe_corr = rc_xe + driftDir_end*anode_rc_time*DriftVelocity;

                // FLASH MATCHING

                size_t op_match_result = FlashMatch((anode_rc_time),op_times);

                if(op_match_result==99999) {
                        if(fDebug) std::cout << "Unable to match flash to track." << std::endl;
                        continue;
                        }

                const art::Ptr<recob::OpFlash> flash_ptr(flash_h, op_match_result);

                matched_flash_time = flash_ptr->Time() - trigger_time;
                if(!MC) corrected_matched_flash_time = fFlashScaleFactor*matched_flash_time + fFlashTPCOffset;
                if(MC) corrected_matched_flash_time = fFlashScaleFactor*matched_flash_time + fFlashTPCOffset - TPC_trigger_offset;
                //matched_flash_time_width = flash_ptr->TimeWidth();

                dt_flash_reco = corrected_matched_flash_time - anode_rc_time;

                matched_flash_pe = flash_ptr->TotalPE();

                //matched_flash_centre_y = flash_ptr->YCenter();
                //matched_flash_centre_z = flash_ptr->ZCenter();
                //matched_flash_width_y = flash_ptr->YWidth();
                //matched_flash_width_z = flash_ptr->ZWidth();

                /*unsigned int max_pe_channel = 9999;
                double max_pe = 0;
                unsigned int pd_ch;
                for(pd_ch = 0; pd_ch <= geom->MaxOpChannel(); pd_ch++) {
                        double channel_i_pe = flash_ptr->PE(pd_ch);
                        if(channel_i_pe > max_pe)  {
                                max_pe = channel_i_pe;
                                max_pe_channel = pd_ch;
                                }
                        }

                if(max_pe_channel==9999||max_pe==0) continue;                   

                matched_flash_max_pe_det_x = -det_width;
                if(max_pe_channel>143) matched_flash_max_pe_det_x = det_width;
                */
                /*double max_pe_det_v[3];
                        geom->OpDetGeoFromOpChannel(max_pe_channel).GetCenter(max_pe_det_v);
                        matched_flash_max_pe_det_x = max_pe_det_v[0];
                        matched_flash_max_pe_det_y = max_pe_det_v[1]; 
                        matched_flash_max_pe_det_z = max_pe_det_v[2];
                */

                /*if(fDebug) std::cout << "\t\tOpChannel " << max_pe_channel << 
                        " has maximum PE, and is located at: (" << 
                        matched_flash_max_pe_det_x << ", " << 
                        matched_flash_max_pe_det_y << ", " << 
                        matched_flash_max_pe_det_z << ")" << std::endl; */

                if(fDebug) std::cout << "\t\t Matched to flash w/ index " << op_match_result 
                        << " w/ PE " << matched_flash_pe << ", corrected time " << 
                        corrected_matched_flash_time << " us vs corrected reco time " << 
                        anode_rc_time << " us" << std::endl;

                // ---------------------------------------------------------------
                // CREATE T0 OBJECT AND ASSIGN TO PFPARTICLE

                if(length>fMinTrackLength&&matched_flash_pe>fMinPE&&
                        dt_flash_reco>dtMin&&dt_flash_reco<dtMax) {
                        if(fDebug) std::cout << "\t\tPFParticle: " << particle << 
                        " has a matched flash and passes cuts. Assigning T0: " 
                        << anode_rc_time << " us." << std::endl;

                        anab::T0 anode_t0 = anab::T0((anode_rc_time*1000), 0, 1, flash_ptr.key());
                t0_v->push_back(anode_t0);

                // Use the association utility
                // Need to convert our recob::PFParticle to an art::Ptr<recob::PFParticle
                art::Ptr<recob::PFParticle> pfp_ptr(reco_particles_h,particle);

                util::CreateAssn(*this, event, *t0_v, pfp_ptr, *pfp_t0);
                        }


                track_number++;

                }

        event.put(std::move(t0_v));
        event.put(std::move(pfp_t0));

        }

void    T0RecoAnodePiercers::SortTrackPoints(const recob::Track& track, std::vector<TVector3>& sorted_trk)
{
        sorted_trk.clear();

        TVector3 track_start, track_end;        
        double start_y = det_bottom - fEdgeWidth;
        double end_y = det_top + fEdgeWidth;

        for (size_t ii = 0; ii < track.NumberTrajectoryPoints(); ii++){
                auto const& trk_loc = track.LocationAtPoint(ii);

                if ((trk_loc.X() < -998.)||(trk_loc.Y() < -998.)||(trk_loc.Z() < -998)) continue;

                if (trk_loc.Y() < end_y){
                        end_y = trk_loc.Y();
                        track_end = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
                        }

                if (trk_loc.Y() > start_y){
                        start_y = trk_loc.Y();
                        track_start = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
                        }
                }

        sorted_trk.push_back(track_start);
        sorted_trk.push_back(track_end);
        }

size_t  T0RecoAnodePiercers::FlashMatch(const double reco_time, std::vector<double> op_times_v)
{
        // loop through all reco'd flash times and see if one matches
        // the time from the track/particle

        double dt_min = 9999999.; 
        size_t matched_op_id = 99999;

        for (size_t i=0; i < op_times_v.size(); i++){
                double op_time_i = op_times_v[i];
                double corrected_op_time_i = op_time_i*fFlashScaleFactor + fFlashTPCOffset;
                double corrected_flash_reco_time_diff = corrected_op_time_i - reco_time;
                if (fabs(corrected_flash_reco_time_diff) < dt_min){
                        dt_min  = fabs(corrected_flash_reco_time_diff);
                        matched_op_id = flash_id_v[i];
                        }
                }

        return matched_op_id;
        }

DEFINE_ART_MODULE(T0RecoAnodePiercers)