307 std::vector<std::vector<TLorentzVector>> mcpart_list;
324 std::cerr<<
"\033[93m[ERROR]\033[00m ... could not locate Flash!"<<
std::endl;
335 trigger_time = trigger_h->at(0).Time();
340 auto track_h =
e.getHandle<std::vector<recob::Track> >(
fTrackProducer);
343 if(!track_h.isValid()) {
344 std::cerr<<
"\033[93m[ERROR]\033[00m ... could not locate Track!"<<
std::endl;
348 std::vector<art::Ptr<recob::Track> > TrkVec;
352 art::FindMany<recob::Hit> trk_hit_assn_v(track_h,
e,
fHitProducer);
355 art::FindMany<anab::Calorimetry> trk_calo_assn_v(track_h,
e,
fCaloProducer);
358 std::cout <<
"Number of tracks: " << TrkVec.size() <<
std::endl;
359 std::cout <<
"Number of hits: " << trk_hit_assn_v.size() <<
std::endl;
360 std::cout <<
"Number of calorimetry products: " << trk_calo_assn_v.size() <<
"\n" <<
std::endl;
364 auto mcpart_h =
e.getHandle<std::vector<simb::MCParticle> >(
"largeant");
369 if(!mcpart_h.isValid()) {
370 std::cerr<<
"\033[93m[ERROR]\033[00m ... could not locate MCParticle!"<<
std::endl;
378 art::FindMany<anab::T0> trk_t0_assn_v(track_h,
e,
fT0Producer );
383 size_t flash_ctr = 0;
384 for (
auto const& flash : *flash_h){
385 if (flash.TotalPE() >
fPEmin){
388 if (
_debug) std::cout <<
"\t flash time : " << flash.Time() - trigger_time <<
", PE : " << flash.TotalPE() <<
std::endl;
397 size_t trk_ctr2 = -1;
399 for (
auto& track : TrkVec){
404 const std::vector<const recob::Hit*>& Hit_v = trk_hit_assn_v.at(trk_ctr2);
406 const std::vector<const anab::Calorimetry*>& Calo_v = trk_calo_assn_v.at(trk_ctr2);
447 std::vector<TVector3> sorted_trk;
449 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;
451 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 ){
457 auto const* geom = lar::providerFrom<geo::Geometry>();
458 auto const*
hit = Hit_v.at(0);
460 const auto TPCGeoObject = geom->TPC(wireID.
TPC,wireID.
Cryostat);
461 short int driftDir1 = TPCGeoObject.DetectDriftDirection();
462 bool cross_cathode =
false;
463 for (
size_t ii = 1; ii < Hit_v.size(); ii++) {
464 const geo::WireID wireID2 = Hit_v.at(ii)->WireID();
465 const auto TPCGeoObject2 = geom->TPC(wireID2.
TPC,wireID2.
Cryostat);
466 short int driftDir2 = TPCGeoObject2.DetectDriftDirection();
468 if(driftDir1 + driftDir2 == 0){
469 cross_cathode =
true;
475 if(
_debug) std::cout <<
"\tCross cathode = " << cross_cathode <<
std::endl;
505 std::vector<TVector3> top_trk = {sorted_trk.at(0), sorted_trk.at(1)};
506 std::vector<TVector3> bottom_trk = {sorted_trk.at(2), sorted_trk.at(3)};
508 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;
514 auto const &
top = top_trk.at(0);
515 auto const &
bottom = top_trk.at(1);
537 auto const &top2 = bottom_trk.at(0);
538 auto const &bottom2 = bottom_trk.at(1);
557 for (
size_t ii = 0; ii < Calo_v.size(); ii++) {
563 for (
size_t jj = 0; jj < calo_obj->
dEdx().size(); jj++){
564 x = calo_obj->
XYZ().at(jj).X();
565 y = calo_obj->
XYZ().at(jj).Y();
566 z = calo_obj->
XYZ().at(jj).Z();
582 if(
fCathode) {std::cout <<
"\tSKIPPING ANODE-PEIRCING TRACKS" <<
std::endl;
continue;}
584 if(
_debug) std::cout <<
"\t\tThis track starts in TPC " << wireID.
TPC <<
" which has a drift direction of " << driftDir1 <<
std::endl;
588 auto const &
top = sorted_trk.at(0);
589 auto const &
bottom = sorted_trk.at(1);
602 evttime = tts.AsDouble();
644 if (tagged ==
false)
continue;
688 if (tagged ==
false)
continue;
710 std::cout <<
"\t this track has a reconstructed time = " <<
rc_time <<
std::endl;
717 for (
auto& hits : Hit_v){
718 auto peakHit = hits->PeakTime();
720 if( (peakHit > detProp.ReadOutWindowSize() - 50.0) || (peakHit < 50.0) ){
722 if(
_debug) std::cout <<
"\t\tHit time out of range (0 - " << detProp.ReadOutWindowSize() <<
"): " << peakHit <<
std::endl;
736 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;
747 double displacement = 1.0;
749 for (
size_t j=0; j < mcpart_list.size(); j++){
751 std::vector<TLorentzVector> mcpart = mcpart_list.at(j);
757 if ((tmp_disp > displacement)||isnan(tmp_disp))
761 displacement = tmp_disp;
762 mc_time = mcpart.at(0).T() / 1000.;
764 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;
772 for (
size_t ii = 0; ii < Calo_v.size(); ii++) {
779 for (
size_t jj = 0; jj < calo_obj->
dEdx().size(); jj++){
780 x = calo_obj->
XYZ().at(jj).X();
781 y = calo_obj->
XYZ().at(jj).Y();
782 z = calo_obj->
XYZ().at(jj).Z();
constexpr std::uint32_t timeLow() const
double GetEnteringTimeCoord(const std::vector< TVector3 > &sorted_trk)
void SplitTrack(const recob::Track &track, std::vector< TVector3 > &sorted_trk)
constexpr std::uint32_t timeHigh() const
bool TrackExitsBack(const std::vector< TVector3 > &sorted_trk)
const std::vector< Point_t > & XYZ() const
CryostatID_t Cryostat
Index of cryostat.
std::vector< size_t > flash_idx_v
double MatchTracks(std::vector< TVector3 > &sorted_trk, std::vector< TLorentzVector > mcpart)
std::pair< double, size_t > FlashMatch(const double reco_time)
const std::vector< float > & ResidualRange() const
std::string fTriggerProducer
bool isValid() const noexcept
timescale_traits< TriggerTimeCategory >::time_point_t trigger_time
A point in time on the trigger time scale.
bool TrackExitsAnode(const std::vector< TVector3 > &sorted_trk, const int driftDir)
const std::vector< float > & dQdx() const
void SortTrackPoints(const recob::Track &track, std::vector< TVector3 > &sorted_trk)
bool TrackEntersSide(const std::vector< TVector3 > &sorted_trk)
bool TrackExitsFront(const std::vector< TVector3 > &sorted_trk)
std::string fFlashProducer
bool TrackExitsBottom(const std::vector< TVector3 > &sorted_trk)
bool TrackEntersAnode(const std::vector< TVector3 > &sorted_trk, const int driftDir)
std::vector< std::vector< TLorentzVector > > BuildMCParticleList(const art::Handle< std::vector< simb::MCParticle >> &mcpart_h, const double &fTPCResolution, const double &width)
bool TrackEntersTop(const std::vector< TVector3 > &sorted_trk)
Detector simulation of raw signals on wires.
bool TrackEntersBack(const std::vector< TVector3 > &sorted_trk)
const std::vector< float > & dEdx() const
std::string fCaloProducer
std::vector< double > flash_times
const std::vector< float > & TrkPitchVec() const
bool TrackEntersFront(const std::vector< TVector3 > &sorted_trk)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
bool TrackExitsSide(const std::vector< TVector3 > &sorted_trk)
const float & KineticEnergy() const
constexpr WireID()=default
Default constructor: an invalid TPC ID.
double GetExitingTimeCoord(const std::vector< TVector3 > &sorted_trk)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
TPCID_t TPC
Index of the TPC within its cryostat.
std::string fTrackProducer
const float & Range() const
cet::coded_exception< error, detail::translate > exception
QTextStream & endl(QTextStream &s)