T0RecoSCECalibrations_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       T0RecoSCECalibrations
// Module Type: analyzer
// File:        T0RecoSCECalibrations_module.cc
//
//
// Hannah Rogers   - hannah.rogers@colostate.edu
// based on uboonecode modules by David Caratelli and Chris Barnes
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.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"

// data-products
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/OpFlash.h"
#include "lardataobj/AnalysisBase/T0.h"
#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/AnalysisBase/CosmicTag.h"
#include "lardataobj/MCBase/MCTrack.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "lardata/Utilities/AssociationUtil.h"

// ROOT
#include "TVector3.h"
#include <TTree.h>
#include "TTimeStamp.h"

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

class T0RecoSCECalibrations;

class T0RecoSCECalibrations : public art::EDAnalyzer {
public:
  explicit T0RecoSCECalibrations(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for calsses
  // without bare pointers or othe resource use.
        
  // Plugins should not be copied or assigned.
  T0RecoSCECalibrations(T0RecoSCECalibrations const &) = delete;
  T0RecoSCECalibrations(T0RecoSCECalibrations &&) = delete;
  T0RecoSCECalibrations & operator = (T0RecoSCECalibrations const &) = delete;
  T0RecoSCECalibrations & operator = (T0RecoSCECalibrations &&) = delete;
        
  void beginJob() override;
        
  //Required functions.
  void analyze(art::Event const & e) override;
        
private:

  // Delcare member data here.
  std::string fTrackProducer;
  std::string fT0Producer;
  std::string fFlashProducer;
  std::string fHitProducer;
  std::string fTriggerProducer;
  std::string fCaloProducer;
        
  bool      fUseMC;
  double      fTPCResolution; // [cm]
  double      fDriftVelocity; // [cm/us]
        
  bool _debug;
  bool fCathode;
  bool fData;
        
  double TOP, BOTTOM, FRONT, BACK, det_width; // [cm]
        
  std::vector<double> flash_times;
  std::vector<size_t> flash_idx_v;
        
  double fTimeRes;
        
  double fPEmin;
        
  double fRecoT0TimeOffset;
        
  //functions to be used throughout module
  bool   TrackEntersTop     (const std::vector<TVector3>& sorted_trk);
  bool   TrackEntersFront   (const std::vector<TVector3>& sorted_trk);
  bool   TrackEntersBack    (const std::vector<TVector3>& sorted_trk);
  bool   TrackEntersAnode   (const std::vector<TVector3>& sorted_trk, const int driftDir);
  bool   TrackEntersSide    (const std::vector<TVector3>& sorted_trk);
  bool   TrackExitsBottom   (const std::vector<TVector3>& sorted_trk);
  bool   TrackExitsFront    (const std::vector<TVector3>& sorted_trk);
  bool   TrackExitsBack     (const std::vector<TVector3>& sorted_trk);
  bool   TrackExitsAnode    (const std::vector<TVector3>& sorted_trk, const int driftDir);
  bool   TrackExitsSide     (const std::vector<TVector3>& sorted_trk);
  
  void   SortTrackPoints      (const recob::Track& track, std::vector<TVector3>& sorted_trk);
  void   SplitTrack(const recob::Track& track, std::vector<TVector3>& sorted_trk);
        
  double GetEnteringTimeCoord (const std::vector<TVector3>& sorted_trk);
  double GetExitingTimeCoord  (const std::vector<TVector3>& sorted_trk);
        
  std::pair<double,size_t> FlashMatch(const double reco_time);
        
  double MatchTracks(std::vector<TVector3>& sorted_trk, std::vector<TLorentzVector> mcpart);
  std::vector<std::vector< TLorentzVector>> BuildMCParticleList(const art::Handle<std::vector<simb::MCParticle>>& mcpart_h, const double& fTPCResolution, const double& width);
        
  // Tree parameters
  TTree *tree;
  double mc_time;
  double rc_time;
  double t_match;
  double pe_flash;
  double dt_flash;
  double dt_mc;
        
  double evttime;
  int year_month_day;
  int hour_min_sec;
        
  double length;
  double driftDir;
  bool TPC_edge;
  double rc_xs, rc_xe, rc_xs_corr, rc_xe_corr;
  double rc_ys, rc_ye;
  double rc_zs, rc_ze;
        
  int anode;
  int cathode;
  int sister_track;
        
  int run, subrun, event, trackNum;
        
  TTree *evTree;
  int trk_ctr;
  int ev_ctr;
        
  TTree *caloTree;
  //TTree *caloSCETree;
  int plane;
  double x, y, z;
  double dQdx, dEdx, dx, resRange;
  double trkLength, kinE;
  double vt0;
         
        
};

T0RecoSCECalibrations::T0RecoSCECalibrations(fhicl::ParameterSet const & p)
  :
  EDAnalyzer(p)
{

  fTrackProducer     = p.get<std::string>("TrackProducer"    );
  fHitProducer       = p.get<std::string>("HitProducer"      );
  fFlashProducer     = p.get<std::string>("FlashProducer"    );
  fT0Producer        = p.get<std::string>("T0Producer"       );
  fTriggerProducer   = p.get<std::string>("TriggerProducer"  );
  fCaloProducer    = p.get<std::string>("CaloProducer"     );
  fUseMC             = p.get<bool>       ("UseMC"            );
  fTPCResolution     = p.get<double>     ("Resolution"       );
  fTimeRes           = p.get<double>     ("TimeRes"          );
  fRecoT0TimeOffset  = p.get<double>     ("RecoT0TimeOffset" );
  fPEmin             = p.get<double>     ("PEmin"            );
  fCathode           = p.get<bool>       ("CathodeOnly"      );
  _debug             = p.get<bool>       ("debug"            );
  fData              = p.get<bool>       ("Data"             );
        
        
  // get boundaries based on detector bounds
  auto const* geom = lar::providerFrom<geo::Geometry>();
  
  TOP = fTPCResolution;
  BOTTOM = fTPCResolution;
  FRONT = fTPCResolution;
  BACK = fTPCResolution;

  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 top = center[1] + TPC.HalfHeight() - fTPCResolution;
    double bottom = center[1] - TPC.HalfHeight() + fTPCResolution;
    double front = center[2] - TPC.HalfLength() + fTPCResolution;
    double back = center[2] + TPC.HalfLength() - fTPCResolution;
   
    if (top > TOP) TOP = top;
    if (bottom < BOTTOM) BOTTOM = bottom;
    if (front < FRONT) FRONT = front;
    if (back  > BACK) BACK = 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();
  double temp   = detProp.Temperature();
  // Determine the drift velocity from 'efield' and 'temp'
  fDriftVelocity = detProp.DriftVelocity(efield,temp);
}  
  
void T0RecoSCECalibrations::beginJob(){

  art::ServiceHandle<art::TFileService> tfs;
  tree = tfs->make<TTree>("tree","SCE calibrations variables");
  tree->Branch("mc_time",&mc_time,"mc_time/D");
  tree->Branch("rc_time",&rc_time,"rc_time/D");
  tree->Branch("t_match",&t_match,"t_match/D");
  tree->Branch("dt_flash",&dt_flash,"dt_flash/D");
  tree->Branch("dt_mc",&dt_flash,"dt_mc/D");
  tree->Branch("pe_flash",&pe_flash,"pe_flash/D");
  tree->Branch("length", &length, "length/D");
  tree->Branch("TPC_edge",&TPC_edge,"TPC_edge/B");
  tree->Branch("driftDir",&driftDir,"driftDir/I");
  //Add branches for time of the track
  tree->Branch("evttime", &evttime, "evttime/D");
  tree->Branch("year_month_day", &year_month_day, "year_month_day/I");
  tree->Branch("hour_min_sec", &hour_min_sec, "hour_min_sec/I");
  // Add branches for the first and last x, y, and z coordinates of the rc tracks and the mc tracks
  tree->Branch("rc_xs",&rc_xs,"rc_xs/D");
  tree->Branch("rc_xs_corr",&rc_xs_corr,"rc_xs/D");
  tree->Branch("rc_ys",&rc_ys,"rc_ys/D");
  tree->Branch("rc_zs",&rc_zs,"rc_zs/D");
  tree->Branch("rc_xe",&rc_xe,"rc_xe/D");
  tree->Branch("rc_xe_corr",&rc_xe_corr,"rc_xe_corr/D");
  tree->Branch("rc_ye",&rc_ye,"rc_ye/D");
  tree->Branch("rc_ze",&rc_ze,"rc_ze/D");
  // information on whether track enters/exits which sides
  tree->Branch("anode",&anode,"anode/I"  );
  tree->Branch("cathode",&cathode,"cathode/I");
  tree->Branch("sister_track",&sister_track,"sister_track/I");
  tree->Branch("run",&run,"run/I");
  tree->Branch("subrun",&subrun,"subrun/I");
  tree->Branch("event",&event,"event/I");
  tree->Branch("trackNum",&trackNum,"trackNum/I");
        
  evTree = tfs->make<TTree>("evTree","Event information");
  evTree->Branch("trk_ctr",&trk_ctr,"trk_ctr/I");
  evTree->Branch("ev_ctr",&ev_ctr,"ev_ctr/I");
  ev_ctr = 0;
  trk_ctr = 0;
        
  caloTree = tfs->make<TTree>("caloTree","Uncorrected calorimetry information");
  caloTree->Branch("anode",&anode,"anode/I"  );
  caloTree->Branch("cathode",&cathode,"cathode/I");
  caloTree->Branch("run",&run,"run/I");
  caloTree->Branch("subrun",&subrun,"subrun/I");
  caloTree->Branch("event",&event,"event/I");
  caloTree->Branch("trackNum",&trackNum,"trackNum/I");
  caloTree->Branch("plane",&plane,"plane/I");
  caloTree->Branch("x",&x,"x/D");
  caloTree->Branch("y",&y,"y/D");
  caloTree->Branch("z",&z,"z/D");
  caloTree->Branch("vt0",&vt0,"vt0/D");
  caloTree->Branch("dQdx",&dQdx,"dQdx/D");
  caloTree->Branch("dEdx",&dEdx,"dEdx/D");
  caloTree->Branch("dx",&dx,"dx/D");
  caloTree->Branch("resRange",&resRange,"resRange/D");
  caloTree->Branch("trkLength",&trkLength,"trkLength/D");
  caloTree->Branch("kinE",&kinE,"kinE/D");
  caloTree->Branch("dt_flash",&dt_flash,"dt_flash/D");
  caloTree->Branch("TPC_edge",&TPC_edge,"TPC_edge/B");
        
  /*caloSCETree = tfs->make<TTree>("caloTree","SCE corrected calorimetry information");
    caloSCETree->Branch("run",&run,"run/I");
    caloSCETree->Branch("subrun",&subrun,"subrun/I");
    caloSCETree->Branch("event",&event,"event/I");
    caloSCETree->Branch("trackNum",&trackNum,"trackNum/I");
    caloSCETree->Branch("x",&x,"x/D");
    caloSCETree->Branch("y",&y,"y/D");
    caloSCETree->Branch("z",&z,"z/D");
    caloSCETree->Branch("dQdx",&dQdx,"dQdx/D");
    caloSCETree->Branch("dEdx",&dEdx,"dEdx/D");
    caloSCETree->Branch("dx",&dx,"dx/D");
    caloSCETree->Branch("resRange",&resRange,"resRange/D");
    caloSCETree->Branch("trkLength",&trkLength,"trkLength/D");
    caloSCETree->Branch("kinE",&kinE,"kinE/D");*/
        
}
  
void T0RecoSCECalibrations::analyze(art::Event const & e){


  event = e.event();
  subrun = e.subRun();
  run = e.run();
  trackNum = 0;
  ev_ctr++;
        
  std::vector<std::vector<TLorentzVector>> mcpart_list;
        
  if(_debug) std::cout << "NEW EVENT" << std::endl;
  if (_debug) std::cout << "top: " << TOP << "\nbottom: " << BOTTOM << "\nfront: " << FRONT << "\nback: " << BACK << std::endl;  
  
  
  flash_times.clear();
  flash_idx_v.clear();
        
  // load Flash
  if (_debug) { std::cout << "loading flash from producer " << fFlashProducer << std::endl; }
  art::Handle<std::vector<recob::OpFlash> > flash_h;
  if(!fCathode){
    flash_h = e.getHandle<std::vector<recob::OpFlash> >(fFlashProducer);

    // make sure flash looks good
    if(!flash_h.isValid()) {
      std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate Flash!"<<std::endl;
      throw std::exception();
    }
  }
        
  if(_debug&&fData) std::cout << "loading trigger time from producer " << fTriggerProducer << std::endl;
  art::Handle<std::vector<recob::OpFlash> > trigger_h;
  double trigger_time = 0;

  if(fData){
    trigger_h = e.getHandle<std::vector<recob::OpFlash> >(fTriggerProducer);
    trigger_time = trigger_h->at(0).Time();
  }

  // load tracks previously created for which T0 reconstruction should occur
  if (_debug) { std::cout << "loading track from producer " << fTrackProducer << std::endl; }
  auto track_h = e.getHandle<std::vector<recob::Track> >(fTrackProducer);

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

  std::vector<art::Ptr<recob::Track> > TrkVec;
  art::fill_ptr_vector(TrkVec, track_h);

  // grab 2d hits associated with tracks
  art::FindMany<recob::Hit> trk_hit_assn_v(track_h, e, fHitProducer);
        
  //load calorimetry data products associated with tracks
  art::FindMany<anab::Calorimetry> trk_calo_assn_v(track_h,e,fCaloProducer);
        
  if(_debug){
    std::cout << "Number of tracks: " << TrkVec.size() << std::endl;
    std::cout << "Number of hits: " << trk_hit_assn_v.size() << std::endl;
    std::cout << "Number of calorimetry products: " << trk_calo_assn_v.size() << "\n" << std::endl;
  }

  // load MCParticles
  auto mcpart_h = e.getHandle<std::vector<simb::MCParticle> >("largeant");

  // if we should use MCParticle
  if (fUseMC){
    // make sure particles exist
    if(!mcpart_h.isValid()) {
      std::cerr<<"\033[93m[ERROR]\033[00m ... could not locate MCParticle!"<<std::endl;
      throw std::exception();
    }

    mcpart_list = BuildMCParticleList(mcpart_h,fTPCResolution,det_width); 
  }// if use MCParticle

  // grab T0 objects associated with tracks
  art::FindMany<anab::T0> trk_t0_assn_v(track_h, e, fT0Producer );      

  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
  // prepare a vector of optical flash times, if flash above some PE cut value
  if(!fCathode){
    size_t flash_ctr = 0;
    for (auto const& flash : *flash_h){
      if (flash.TotalPE() > fPEmin){
        flash_times.push_back( flash.Time() - trigger_time);
        flash_idx_v.push_back(flash_ctr);
        if (_debug) std::cout << "\t flash time : " << flash.Time() - trigger_time << ", PE : " << flash.TotalPE() << std::endl;
      }
      flash_ctr += 1;
    }// for all flashes
  
    if (_debug) { std::cout << "Selected a total of " << flash_times.size() << " OpFlashes" << std::endl; }
  }
  
  // loop through reconstructed tracks
  size_t trk_ctr2 = -1;

  for (auto& track : TrkVec){ 
    trk_ctr2 ++;
    trk_ctr ++;
    if (_debug) std::cout << "Looping through reco track " << trk_ctr2 << " or " << trk_ctr << std::endl;  
  
    const std::vector<const recob::Hit*>& Hit_v = trk_hit_assn_v.at(trk_ctr2);
    if(_debug) std::cout << "     Hits loaded" << std::endl;
    const std::vector<const anab::Calorimetry*>& Calo_v = trk_calo_assn_v.at(trk_ctr2);
    if(_debug) std::cout << "     Calorimetry loaded" << std::endl;
                
    mc_time = 0.;
    rc_time = 0.;
    t_match = 0.;
    pe_flash = 0.;
    dt_flash = -9999.;
    dt_mc = -9999.;
    length = 0.;
    driftDir = 0;
    rc_xs = 0.;
    rc_xe = 0.; 
    rc_xs_corr = 0.;
    rc_xe_corr = 0.;
    rc_ys = 0.;
    rc_ye =0.; 
    rc_zs = 0.; 
    rc_ze = 0.;
    anode = 0;
    cathode = 0;
    sister_track = 0;
    TPC_edge = false;
                
    plane = -9999;
    trkLength = 0.;
    kinE = 0.;
    x = 0.;
    y = 0.;
    z = 0.;
    dx = 0.;
    dEdx = 0.;
    dQdx = 0.;
    resRange = 0.;
    vt0 = 0.;
        
    evttime = -9999.;
    year_month_day = -9999;
    hour_min_sec = -9999;
  
    // get sorted points for the track object [assuming downwards going]
    std::vector<TVector3> sorted_trk;
    SortTrackPoints(*track,sorted_trk);
    if(_debug) std::cout << "\tTrack goes from (" << sorted_trk.at(0).X() << ", " << sorted_trk.at(0).Y() << ", " << sorted_trk.at(0).Z() << ") --> (" << sorted_trk.at(sorted_trk.size()-1).X() << ", " << sorted_trk.at(sorted_trk.size()-1).Y() << ", " << sorted_trk.at(sorted_trk.size()-1).Z() << ")" << std::endl;
      
    if( sqrt(pow(sorted_trk.at(0).X() - sorted_trk.at(sorted_trk.size()-1).X(),2.0) + pow(sorted_trk.at(0).Y() - sorted_trk.at(sorted_trk.size()-1).Y(),2.0) + pow(sorted_trk.at(0).Z() - sorted_trk.at(sorted_trk.size()-1).Z(),2.0)) < 50 ){
      if(_debug) std::cout << "\tTrack too short. Skipping." << std::endl;
      continue;
    }

    // Determine if the track crosses the cathode 
    auto const* geom = lar::providerFrom<geo::Geometry>();   
    auto const* hit = Hit_v.at(0);
    const geo::WireID wireID = hit->WireID();
    const auto TPCGeoObject = geom->TPC(wireID.TPC,wireID.Cryostat);
    short int driftDir1 = TPCGeoObject.DetectDriftDirection();
    bool cross_cathode = false;
    for (size_t ii = 1; ii < Hit_v.size(); ii++) {
      const geo::WireID wireID2 = Hit_v.at(ii)->WireID();
      const auto TPCGeoObject2 = geom->TPC(wireID2.TPC,wireID2.Cryostat);
      short int driftDir2 = TPCGeoObject2.DetectDriftDirection(); 
                
      if(driftDir1 + driftDir2 == 0){
        cross_cathode = true;
        if(_debug) std::cout << "\tCrosses cathode!" << std::endl;
        //continue;
        ii = Hit_v.size();
      }
    }
    if(_debug) std::cout << "\tCross cathode = " << cross_cathode << std::endl;
                
        
    // -------------------------------------------------------------------------------
    //CATHODE CROSSERS
    if(cross_cathode){ 
      anode = 0;
      cathode = 1;
                
      //const std::vector<const anab::T0*>& T0_v = trk_t0_assn_v.at(trk_ctr);
      //auto t0 = T0_v.at(0);
      //rc_time = t0->Time();
                
        
      //Get Time and date
      art::Timestamp ts = e.time();
      TTimeStamp tts(ts.timeHigh(), ts.timeLow());
      evttime = tts.AsDouble();

      UInt_t year = 0;
      UInt_t month = 0;
      UInt_t day = 0;
      year_month_day = tts.GetDate(kTRUE,0,&year,&month,&day);

      UInt_t hour = 0;
      UInt_t min = 0;
      UInt_t sec = 0;
      hour_min_sec = tts.GetTime(kTRUE,0,&hour,&min,&sec);
        
      SplitTrack(*track,sorted_trk);
      std::vector<TVector3> top_trk = {sorted_trk.at(0), sorted_trk.at(1)};
      std::vector<TVector3> bottom_trk = {sorted_trk.at(2), sorted_trk.at(3)};
        
      if(_debug)std::cout << "\tCathode-crossing track from (" << top_trk.at(0).X() << ", " << top_trk.at(0).Y() << ", " << top_trk.at(0).Z() << ") --> (" << top_trk.at(1).X() << ", " << top_trk.at(1).Y() << ", " << top_trk.at(1).Z() << ") --> (" << bottom_trk.at(0).X() << ", " << bottom_trk.at(0).Y() << ", " << bottom_trk.at(0).Z() << ") --> (" << bottom_trk.at(1).X() << ", " << bottom_trk.at(1).Y() << ", " << bottom_trk.at(1).Z() << ")" << std::endl;
        
      // Top Track!
      if(_debug) std::cout << "\tTop track" << std::endl;
        
      sister_track = trackNum + 1;
      auto const &top = top_trk.at(0);
      auto const &bottom = top_trk.at(1);
                
      rc_xs = top.X();
      rc_xs_corr = rc_xs;
      rc_ys = top.Y();
      rc_zs = top.Z();
      rc_xe = bottom.X();
      rc_xe_corr = rc_xe;
      rc_ye = bottom.Y();
      rc_ze = bottom.Z();
        
      length = sqrt(pow(rc_xs - rc_xe,2.0) + pow(rc_ys - rc_ye,2.0) + pow(rc_zs - rc_ze,2.0));
      if(rc_xs<0.0) driftDir = -1;
      else driftDir = 1;
                
      tree->Fill();
      trackNum++;
                
      // Bottom Track
      if(_debug) std::cout << "\tBottom track" << std::endl;
      sister_track = trackNum - 1;
                
      auto const &top2 = bottom_trk.at(0);
      auto const &bottom2 = bottom_trk.at(1);
                
      rc_xs = top2.X();
      rc_xs_corr = rc_xs;
      rc_ys = top2.Y();
      rc_zs = top2.Z();
      rc_xe = bottom2.X();
      rc_xe_corr = rc_xe;
      rc_ye = bottom2.Y();
      rc_ze = bottom2.Z();
        
      length = sqrt(pow(rc_xs - rc_xe,2.0) + pow(rc_ys - rc_ye,2.0) + pow(rc_zs - rc_ze,2.0));
      if(rc_xs<0.0) driftDir = -1;
      else driftDir = 1;
                
      tree->Fill();
                
                
      //Save calorimetry information
      for (size_t ii = 0; ii < Calo_v.size(); ii++) {
        plane = ii;
        const anab::Calorimetry* calo_obj = Calo_v.at(ii);
        trkLength = calo_obj->Range();
        kinE = calo_obj->KineticEnergy();
                        
        for (size_t jj = 0; jj < calo_obj->dEdx().size(); jj++){
          x = calo_obj->XYZ().at(jj).X();
          y = calo_obj->XYZ().at(jj).Y();
          z = calo_obj->XYZ().at(jj).Z();
                                
          dx = calo_obj->TrkPitchVec().at(jj);
          dEdx = calo_obj->dEdx().at(jj);
          dQdx = calo_obj->dQdx().at(jj);
          resRange = calo_obj->ResidualRange().at(jj);
                                
          caloTree->Fill();
                        
        }
                
      }
      trackNum++;
      // ------------------------------------------------------------------------------------
      // ANODE PIERCERS 
    } else {
      if(fCathode) {std::cout << "\tSKIPPING ANODE-PEIRCING TRACKS" << std::endl; continue;}
        
      if(_debug) std::cout << "\t\tThis track starts in TPC " << wireID.TPC << " which has a drift direction of " << driftDir1 << std::endl; 
      driftDir = driftDir1;
        
      // create root trees variables
      auto const &top = sorted_trk.at(0);
      auto const &bottom = sorted_trk.at(1);
    
      rc_xs = top.X();
      rc_ys = top.Y();
      rc_zs = top.Z();
      rc_xe = bottom.X();
      rc_ye = bottom.Y();
      rc_ze = bottom.Z();
      length = track->Length();

      //Get Time and date
      art::Timestamp ts = e.time();
      TTimeStamp tts(ts.timeHigh(), ts.timeLow());
      evttime = tts.AsDouble();
    

      UInt_t year = 0;
      UInt_t month = 0;
      UInt_t day = 0;
      year_month_day = tts.GetDate(kTRUE,0,&year,&month,&day);

      UInt_t hour = 0;
      UInt_t min = 0;
      UInt_t sec = 0;
      hour_min_sec = tts.GetTime(kTRUE,0,&hour,&min,&sec);
    

      // keep track of whether it goes thorugh the anode or cathode
      anode = 0, cathode = 0;

      // 1st category: tracks which ENTER SIDE
      if ( TrackEntersSide(sorted_trk) == true ) {

        if (_debug) std::cout << "\t track enters side" << std::endl;

        // we are not done. We need to check that the track either: 1) Exits the bottom. 2) exits the front or 3) exits the back of the TPC.
        bool tagged = false;

        // tracks that exit the bottom
        if ( (TrackExitsBottom(sorted_trk) == true) ) {
          tagged = true;
          if (_debug) std::cout << "\t track exits bottom" << std::endl;
        }
        // tracks that exit the front
        if ( (TrackExitsFront(sorted_trk) == true) and (TrackEntersFront(sorted_trk) == false)) {
          tagged = true;
          if (_debug) std::cout << "\t track exits front" << std::endl;
        }
        // tracks that exit the back
        if ( (TrackExitsBack(sorted_trk) == true) and (TrackEntersBack(sorted_trk) == false) ) {
          tagged = true;
          if (_debug) std::cout << "\t track exits back" << std::endl;
        }

        // has either of these 3 conditions been met? if no, skip this track
        if (tagged == false) continue;
      
        // figure out if it enters the anode or cathode
        bool enters_anode = TrackEntersAnode(sorted_trk, driftDir);
      
        // get the X coordinate of the point piercing the anode/cathode (upon ENTERING) 
        double trkX = GetEnteringTimeCoord(sorted_trk);
      
        // reconstruct track T0 w.r.t. trigger time

        // The 'trkX' enters on the anode, the side of the TPC with a lower x value than the cathode
        if (enters_anode){
          rc_time = (det_width - (double)driftDir * trkX) / fDriftVelocity + fRecoT0TimeOffset;
          if(_debug) std::cout << "\tTrack enters anode. " << rc_time << "  = ( " << det_width << " - " << driftDir*trkX << ") / " << fDriftVelocity << std::endl;
          anode = 1;
        }

      }// if the track enters the side

      // case in which the track exits the side
      if (TrackExitsSide(sorted_trk) == true) {

        if (_debug) std::cout << "\t track exits side" << std::endl;

        // we are not done. We need to check that the track either: 1) Enters the bottom. 2) enters the front or 3) enters the back of the TPC.
        bool tagged = false;

        // track enters the top
        if ( (TrackEntersTop(sorted_trk) == true) ) {       
          tagged = true;
          if (_debug) std::cout << "\t track enters the top" << std::endl;
        }

        if ( (TrackEntersFront(sorted_trk) == true) and (TrackExitsFront(sorted_trk) == false)) {
          tagged = true;
          if (_debug) std::cout << "\t track enters front" << std::endl;
        }

        if ( (TrackEntersBack(sorted_trk) == true) and (TrackExitsBack(sorted_trk) == false)  ) {
          tagged = true;
          if (_debug) std::cout << "\t track enters back" << std::endl;
        }

        // has either of these 3 conditions been met? if no, skip this track
        if (tagged == false) continue;
        // figure out if it enters the anode or cathode
        bool exits_anode = TrackExitsAnode(sorted_trk, driftDir);
      
        // get the X coordinate of the point piercing the anode/cathode (upon ENTERING) 
        double trkX = GetExitingTimeCoord(sorted_trk);
      
        // reconstruct track T0 w.r.t. trigger time

        // The 'trkX' enters on the anode, the side of the TPC with a lower x value than the cathode
        if (exits_anode){
          rc_time = (det_width - (double)driftDir * trkX) / fDriftVelocity + fRecoT0TimeOffset;
          if(_debug) std::cout << "\tTrack exits anode. " << rc_time << "  = ( " << det_width << " - " << driftDir*trkX << ") / " << fDriftVelocity << std::endl;
          anode = 1;
        }

      }// if the track exits the side
        
      if(!anode) continue; 
    

      if (_debug)
        std::cout << "\t this track has a reconstructed time = " << rc_time << std::endl;
    
      rc_xs_corr = rc_xs + driftDir*rc_time*fDriftVelocity;
      rc_xe_corr = rc_xe + driftDir*rc_time*fDriftVelocity;
    
      // Determine if track hits edge of readout window
      //bool TPC_edge = false;
      for (auto& hits : Hit_v){
        auto peakHit = hits->PeakTime();
        //if(_debug) std::cout << "\t\tHit time track " << trk_ctr << ": " << peakHit << " in TPC " << hits->WireID().TPC << " plane " << hits->WireID().Plane << " and wire " << hits->WireID().Wire << std::endl;
        if( (peakHit > detProp.ReadOutWindowSize() - 50.0) || (peakHit < 50.0) ){
          TPC_edge = true;
          if(_debug) std::cout << "\t\tHit time out of range (0 - " << detProp.ReadOutWindowSize() << "): " << peakHit << std::endl;
          continue;
        }
      }
    
      /*if (TPC_edge) {
        if(_debug) std::cout << "\tHit time too close to edge" << std::endl;
        continue;
        }*/
        
      // flash matching
      auto const& flash_match_result = FlashMatch(rc_time);
      const art::Ptr<recob::OpFlash> flash_ptr(flash_h, flash_match_result.second );
      if (_debug)
        std::cout << "\t matched to flash w/ index " << flash_match_result.second << " w/ PE " << flash_ptr->TotalPE() << " and time " << flash_ptr->Time() - trigger_time << " vs reco time " << rc_time << std::endl;

    
      pe_flash = flash_ptr->TotalPE();
      t_match = flash_ptr->Time() - trigger_time;
      dt_flash = t_match - rc_time;
    
      // if we should use MC info -> continue w/ MC validation
      if (fUseMC == true){
        // loop through MCParticles to find the one that matches.
        mc_time = 0;
        double displacement = 1.0;
          
        for (size_t j=0; j < mcpart_list.size(); j++){
          
          std::vector<TLorentzVector> mcpart = mcpart_list.at(j);
          
          //if(_debug) std::cout <<"\t\tMCParticle: (" << mcpart.at(0).X() << ", " << mcpart.at(0).Y() << ", " << mcpart.at(0).Z() << ") -> (" << mcpart.at(1).X() << ", " << mcpart.at(1).Y() << ", " << mcpart.at(1).Z() << ")" << std::endl;
          
          // try matching to MC
          double tmp_disp = MatchTracks(sorted_trk, mcpart);
          if ((tmp_disp > displacement)||isnan(tmp_disp))
            continue;
          
          // matched -> get MCTrack time and reconstructed track reconstructed T0
          displacement = tmp_disp;
          mc_time = mcpart.at(0).T() / 1000.;

          if(_debug) std::cout << "\tParticle number: " << j << " with sigma comparison value of " << displacement << "\n\t\tTrack matched to MC value with t0 of " << mc_time << ", flash t0 of " << t_match << " and reconstructed time of " << rc_time << std::endl;
          
        } // for all MCParticles
      }// if we should use MCParticles
      
      tree->Fill();
      
      //Save calorimetry information
      for (size_t ii = 0; ii < Calo_v.size(); ii++) {
        plane = ii;
        const anab::Calorimetry* calo_obj = Calo_v.at(ii);
        trkLength = calo_obj->Range();
        kinE = calo_obj->KineticEnergy();
        vt0 = driftDir*rc_time*fDriftVelocity;
                        
        for (size_t jj = 0; jj < calo_obj->dEdx().size(); jj++){
          x = calo_obj->XYZ().at(jj).X();
          y = calo_obj->XYZ().at(jj).Y();
          z = calo_obj->XYZ().at(jj).Z();
                                
          dx = calo_obj->TrkPitchVec().at(jj);
          dEdx = calo_obj->dEdx().at(jj);
          dQdx = calo_obj->dQdx().at(jj);
          resRange = calo_obj->ResidualRange().at(jj);
                                
          caloTree->Fill();
                        
        }
                
      }
      
      trackNum++;
        
    }
    
  }
    
  evTree->Fill();
}  
  
double T0RecoSCECalibrations::MatchTracks(std::vector<TVector3>& sorted_trk, std::vector<TLorentzVector> mcpart)
{
  //assumes both the reco track
  //and mctrack are downwards going
  
  //auto const& mctrk_s = mcpart.at(0);
  //auto const& mctrk_e = mcpart.at(1);
  //auto const& track_s = sorted_trk.at(0);
  //auto const& track_e = sorted_trk.at(sorted_trk.size()-1);
  
  double mcX = mcpart.at(0).X(), mcY = mcpart.at(0).Y(), mcZ = mcpart.at(0).Z();
  double mcX_end = mcpart.at(1).X(), mcY_end = mcpart.at(1).Y(), mcZ_end = mcpart.at(1).Z();
  double trkX = sorted_trk.at(0).X(), trkY = sorted_trk.at(0).Y(), trkZ = sorted_trk.at(0).Z();
  double trkX_end = sorted_trk.at(sorted_trk.size()-1).X(), trkY_end = sorted_trk.at(sorted_trk.size()-1).Y(), trkZ_end = sorted_trk.at(sorted_trk.size()-1).Z();
  
  if(mcY<mcY_end){
    double tmpX = mcX, tmpY = mcY, tmpZ = mcZ;
    mcX = mcX_end;
    mcY = mcY_end;
    mcZ = mcZ_end;
    mcX_end = tmpX;
    mcY_end = tmpY;
    mcZ_end = tmpZ;
  }
  
  double trkThXZ = atan((trkZ_end - trkZ)/(trkX_end - trkX));
  double trkThYZ = atan((trkZ_end - trkZ)/(trkY_end - trkY));
  double sigThXZ = sqrt(2.0)*fTPCResolution/sqrt(pow(trkZ_end-trkZ,2.0)+pow(trkX_end-trkX,2.0));
  double sigThYZ = sqrt(2.0)*fTPCResolution/sqrt(pow(trkZ_end-trkZ,2.0)+pow(trkY_end-trkY,2.0));
  double mcThXZ = atan((mcZ_end - mcZ)/(mcX_end - mcX));
  double mcThYZ = atan((mcZ_end - mcZ)/(mcY_end - mcY));
  
  return sqrt(pow((trkThXZ -  mcThXZ)/sigThXZ,2.0) + pow((trkThYZ -  mcThYZ)/sigThYZ,2.0) + pow((trkY - mcY)/fTPCResolution,2.0) + pow((trkZ - mcZ)/fTPCResolution,2.0))/4.0;

  
  
  /*
  // if track start is above and mctrk start is above
  if ( ( track_s.Y() > track_e.Y() ) and ( mctrk_s.Y() > mctrk_e.Y() ) ){
  if ( (fabs(mctrk_s.Y()-track_s.Y()) < res) and (fabs(mctrk_s.Z()-track_s.Z()) < res) and (fabs(mctrk_e.Y()-track_e.Y()) < res) and (fabs(mctrk_e.Z()-track_e.Z()) < res) )
  return true;
  }
  // if track start is above and mctrk start is below
  if ( ( track_s.Y() > track_e.Y() ) and ( mctrk_s.Y() < mctrk_e.Y() ) ){
  if ( (fabs(mctrk_e.Y()-track_s.Y()) < res) and (fabs(mctrk_e.Z()-track_s.Z()) < res) and (fabs(mctrk_s.Y()-track_e.Y()) < res) and (fabs(mctrk_s.Z() - track_e.Z()) < res) )
  return true;
  }
  // if track start is below and mctrk start is above
  if ( ( track_s.Y() < track_e.Y() ) and ( mctrk_s.Y() > mctrk_e.Y() ) ){
  if ( (fabs(mctrk_s.Y()-track_e.Y()) < res) and (fabs(mctrk_s.Z()-track_e.Z()) < res) and (fabs(mctrk_e.Y()-track_s.Y()) < res) and (fabs(mctrk_e.Z()-track_s.Z()) < res) )
  return true;
  }
  // if track start is below and mctrk start is below
  if ( ( track_s.Y() < track_e.Y() ) and ( mctrk_s.Y() < mctrk_e.Y() ) ){
  if ( (fabs(mctrk_e.Y()-track_e.Y()) < res) and (fabs(mctrk_e.Z()-track_e.Z()) < res) and (fabs(mctrk_s.Y() - track_s.Y()) < res) and (fabs(mctrk_s.Z()-track_s.Z()) < res) )
  return true;
  }
  
  return false;
  */
}
  
  
std::pair<double,size_t> T0RecoSCECalibrations::FlashMatch(const double reco_time){
  
  // loop through all reco'd flash times and see if one matches
  // the reco time from the track
  double dt_min = 8000.; // us
  size_t idx_min = flash_times.size();

  for (size_t i=0; i < flash_times.size(); i++){
    auto const& time = flash_times[i];
    double dt = fabs(time - reco_time);
    if (dt < dt_min){
      dt_min  = dt;
      idx_min = flash_idx_v[i];
    }
  }

  std::pair<double,size_t> ret(dt_min,idx_min);
  return ret;
}


bool   T0RecoSCECalibrations::TrackEntersTop(const std::vector<TVector3>& sorted_trk)
{
  // check that the first point in the track
  // pierces the top boundary of the TPC
  // This track either will pierce the top of the TPC or is just about to (the '_TOP' variable is just below the actual coordinate position of the top in Y)

  if (sorted_trk.at(0).Y() > TOP)
    return true;

  return false;
}


bool T0RecoSCECalibrations::TrackEntersFront(const std::vector<TVector3>& sorted_trk)
{

  // Determine if the track enters the
  // front of the TPC based on if the position
  // of its initial Z-coordinate is less than
  // the location of the front of the TPC in Z
  
  // First define 'top_pt' to mean the point at the start of the track
  auto const& top_pt = sorted_trk.at(0);

  if (top_pt.Z() < FRONT)
    return true;

  // I may include the case in which I check
  // the y-coordinates as well, but I will not
  // implement that at this time
  
  // If this condition is not satisfied, then return 'false' (the track was not determined
  // within resolution to enter the front of the TPC)
  return false;
}


bool T0RecoSCECalibrations::TrackEntersBack(const std::vector<TVector3>& sorted_trk)
{

  // Determines if the track enters the
  // back of the TPC based on if the position
  // of its initial Z-coordinate is greater
  // than the location of the back of the
  // TPC in Z
  
  // First define 'top_pt' to mean the point at the start of the track
  auto const& top_pt = sorted_trk.at(0);

  if (top_pt.Z() > BACK)
    return true;

  // If this condition is not satisfied, then return 'false' (the track was not determined
  // within resolution to enter the back of the TPC)
  return false;
}


bool   T0RecoSCECalibrations::TrackEntersAnode(const std::vector<TVector3>& sorted_trk, const int driftDir)
{

  // we know the track enters either the
  // anode or cathode
  // at this point figure out
  // if it ENTERS the ANODE or CATHODE
  // ANODE: top point must be at lower X-coord
  // than bottom point
  // CATHODE: top point must be at larger X-coord
  // than bottom point
  // assume track has already been sorted
  // such that the 1st point is the most elevated in Y coord.
  // return TRUE if passes the ANODE
  
  auto const& top    = sorted_trk.at(0);
  auto const& bottom = sorted_trk.at( sorted_trk.size() - 1 );

  if ( ( (top.X() < bottom.X()) && driftDir<0 ) || ( (top.X() > bottom.X()) && driftDir>0 ) )
    return true;

  return false;
}


bool   T0RecoSCECalibrations::TrackEntersSide(const std::vector<TVector3>& sorted_trk)
{
  
  // check that the top-most point
  // is not on the top of the TPC
  // nor on the front & back of the TPC
  
  auto const& top_pt = sorted_trk.at(0);

  // if highest point above the TOP -> false
  if (top_pt.Y() > TOP)
    return false;

  // if highest point in Z close to front or back
  // -> FALSE
  if ( (top_pt.Z() < FRONT) or (top_pt.Z() > BACK) )
    return false;


  // If the function makes it this far, then it will enter through one of the sides of the TPC
  return true;
}


bool   T0RecoSCECalibrations::TrackExitsBottom(const std::vector<TVector3>& sorted_trk)
{

  // check that the last point in the track
  // pierces the bottom boundary of the TPC
  if ( sorted_trk.at( sorted_trk.size() - 1).Y() < BOTTOM )
    return true;

  return false;
}


bool   T0RecoSCECalibrations::TrackExitsFront(const std::vector<TVector3>& sorted_trk)
{

  // Determine if the track exits the
  // front of the TPC based on if the position
  // of its final Z-coordinate is less than
  // the location of the front of the TPC in Z
  
  // First define 'bottom_pt' to mean the point at the end of the track
  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);

  if (bottom_pt.Z() < FRONT)
    return true;
  
  return false;
}


bool   T0RecoSCECalibrations::TrackExitsBack(const std::vector<TVector3>& sorted_trk)
{

  // Determine if the track exits the
  // front of the TPC based on if the position
  // of its final Z-coordinate is less than
  // the location of the front of the TPC in Z
  
  // First define 'bottom_pt' to mean the point at the end of the track
  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);

  if (bottom_pt.Z() > BACK)
    return true;

  return false;
}


bool   T0RecoSCECalibrations::TrackExitsAnode(const std::vector<TVector3>& sorted_trk, const int driftDir)
{

  // Check, once it's known that the track doesn't exit out of the bottom, whether it's the anode or
  // the cathode that it exits out of
  // This can be done by direct analogy with the 'Anode' function (shown in this file as the 'TrackEntersAnode') function written by D. Caratelli
  // Define 'top' as the point at the start of the track, and 'bottom' as the point at the end of the track

  auto const& top    = sorted_trk.at(0);
  auto const& bottom = sorted_trk.at(sorted_trk.size() - 1);

  // Check to see which point has a lower x coordinate
  // If the bottom does, then it exits out of the anode
  // If the top does, then it exits out of the cathode
  if ( ( (bottom.X() < top.X()) && driftDir<0 ) || ( (bottom.X() > top.X()) && driftDir>0 ) ) 
    return true;
        
  return false; // Otherwise, the top is less than the bottom, so the track ended closer to the cathode and exited there
}


bool   T0RecoSCECalibrations::TrackExitsSide(const std::vector<TVector3>& sorted_trk)
{

  // check that the bottom-most point
  // is not on the bottom of the TPC
  // nor on the front & back of the TPC

  auto const& bottom_pt = sorted_trk.at(sorted_trk.size() - 1);

  // if lowest point below the BOTTOM -> false
  // Within this resolution, this means that it's likely that the track exited out of the bottom (at a point earlier on in the process than the last point) OR is just about to

  if (bottom_pt.Y() <  BOTTOM)
    return false;

  // if lowest point in Z close to front or back
  // -> FALSE
  // If the the bottom point is less than the front, then the track has already pierced the front of the TPC and exited that way OR is likely just about to
  // If the bottom point is greater than the back, then the track has already pierced the back of the TPC and exited that way OR is likely just about to
  if ( (bottom_pt.Z() < FRONT) or (bottom_pt.Z() > BACK) )
    return false;

  return true;
}

void T0RecoSCECalibrations::SplitTrack(const recob::Track& track, std::vector<TVector3>& sorted_trk){

  sorted_trk.clear();
        
  TVector3 track_neg, track_neg_c, track_pos_c,track_pos;
  double neg_x = 0.0;
  double pos_x = 0.0;
  double neg_c = -2.0*fTPCResolution;
  double pos_c = 2.0*fTPCResolution;
        
  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.X() < neg_x){
      neg_x = trk_loc.X();
      track_neg = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
    }
    if (trk_loc.X() > pos_x){
      pos_x = trk_loc.X();
      track_pos = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
    }
    if ((trk_loc.X() < 0.0) && (trk_loc.X() > neg_c)){
      neg_c = trk_loc.X();
      track_neg_c = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
    }
    if ((trk_loc.X() > 0.0) && (trk_loc.X() < pos_c)){
      pos_c = trk_loc.X();
      track_pos_c = {trk_loc.X(), trk_loc.Y(), trk_loc.Z()};
    }
  }
        
  if( track_neg.Y() > track_pos.Y()){
    sorted_trk.push_back(track_neg);
    sorted_trk.push_back(track_neg_c);
    sorted_trk.push_back(track_pos_c);
    sorted_trk.push_back(track_pos);
  } else {
    sorted_trk.push_back(track_pos);
    sorted_trk.push_back(track_pos_c);
    sorted_trk.push_back(track_neg_c);
    sorted_trk.push_back(track_neg);
  }
        
}

void   T0RecoSCECalibrations::SortTrackPoints(const recob::Track& track, std::vector<TVector3>& sorted_trk)
{
  sorted_trk.clear();
                
  TVector3 track_start, track_end;      
  double start_y = BOTTOM - 2.0*fTPCResolution;
  double end_y = TOP + 2.0*fTPCResolution;
        
  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);
        
  /*  THIS METHOD ASSUMES THE TRACK IS SORTED AT ALL!!
  // vector to store 3D coordinates of
  // ordered track
  sorted_trk.clear();

  // take the reconstructed 3D track
  // and assuming it is downwards
  // going, sort points so that
  // the track starts at the top
  // which point is further up in Y coord?
  // start or end?
  auto const&N = track.NumberTrajectoryPoints();
  auto const&start = track.LocationAtPoint(0);
  auto const&end   = track.LocationAtPoint( N - 1 );

  // if points are ordered correctly
  if (start.Y() > end.Y()){
  for (size_t i=0; i < N; i++)
  sorted_trk.push_back( track.LocationAtPoint(i) );
  }
  
  // otherwise flip order
  else {
  for (size_t i=0; i < N; i++)
  sorted_trk.push_back( track.LocationAtPoint( N - i - 1) );
  }
  */
}


double T0RecoSCECalibrations::GetEnteringTimeCoord(const std::vector<TVector3>& sorted_trk)
{

  // get the drift-coordinate value
  // associated with the point
  // along the track piercing the anode / cathode
  // ** WHEN the track enters the anode / cathode
  return sorted_trk.at(0).X();
}


double T0RecoSCECalibrations::GetExitingTimeCoord(const std::vector<TVector3>& sorted_trk) 
{
  // get the drift-coordinate value
  // associated with the point
  // along the track piercing the anode / cathode
  // ** WHEN the track exits the anode / cathode
  return sorted_trk.at(sorted_trk.size() - 1).X();
}

std::vector<std::vector< TLorentzVector >> T0RecoSCECalibrations::BuildMCParticleList(const art::Handle<std::vector<simb::MCParticle> >& mcpart_h,const double& fTPCResolution, const double& width){

  std::vector<std::vector< TLorentzVector>> mcVec;

  for (size_t j=0; j < mcpart_h->size(); j++){
                        
    auto const& mcpart = mcpart_h->at(j);
                
    double x0 = mcpart.Position(0).X();
    double y0 = mcpart.Position(0).Y();
    double z0 = mcpart.Position(0).Z();
    double t0 = mcpart.T(0);
    double x1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).X();
    double y1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).Y();
    double z1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-1).Z();
    double t1 = mcpart.T(mcpart.NumberTrajectoryPoints()-1);
    std::cout << 3 << std::endl;
                
    //double xs = x0, ys = y0, zs = z0, ts = t0;
    //double xe = x1, ye = y1, ze = z1, te = t1;
                
    bool top_ok = false, bot_ok = false;
    std::cout << mcpart.NumberTrajectoryPoints() << std::endl;
    if(mcpart.NumberTrajectoryPoints()<2900){
      for(size_t k = 0; k<mcpart.NumberTrajectoryPoints(); k++){
        if ((x0 > width) || (x0 < -width) || (y0 > TOP + fTPCResolution) || (y0 < BOTTOM - fTPCResolution) || (z0 < FRONT - fTPCResolution) || (z0 > BACK + fTPCResolution) ){
          x0 = mcpart.Position(k+1).X();
          y0 = mcpart.Position(k+1).Y();
          z0 = mcpart.Position(k+1).Z();
          t0 = mcpart.T(k+1);
        } else top_ok = true;

        if ((x1 > width) || (x1 < -width) || (y1 > TOP + fTPCResolution) || (y1 < BOTTOM - fTPCResolution) || (z1 < FRONT - fTPCResolution) || (z1 > BACK + fTPCResolution) ){
          x1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).X();
          y1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).Y();
          z1 = mcpart.Position(mcpart.NumberTrajectoryPoints()-2-k).Z();
          t1 = mcpart.T(mcpart.NumberTrajectoryPoints()-2-k);
        } else bot_ok = true;

        if(top_ok&&bot_ok) k = mcpart.NumberTrajectoryPoints();
      }
    }
                        
    double xs = x0, ys = y0, zs = z0, ts = t0;
    double xe = x1, ye = y1, ze = z1, te = t1;
                
    TLorentzVector mc_top(TVector3(xs,ys,zs),ts);
    TLorentzVector mc_bot(TVector3(xe,ye,ze),te);
    //if (!top_ok||!bot_ok) continue;
                
    if (ys >  ye) mcVec.push_back({mc_top, mc_bot});
    else mcVec.push_back({mc_bot, mc_top});
    std::cout << 7 << std::endl;
    //if(_debug&&_particle) std::cout << " MCParticle: (" << x0 << ", " << y0 << ", " << z0 << ", " << t0/1000. << ") --> (" << xs << ", " << ys << ", " << zs << ", " << ts/1000.0 << ") --> (" << xe << ", " << ye << ", " << ze << ", " << te/1000.0 << ") --> (" << x1 << ", " << y1 << ", " << z1 << ", " << t1/1000.0 << ")" << std::endl;
                
  }
  return mcVec; 

}



DEFINE_ART_MODULE(T0RecoSCECalibrations)