Program Listing for File FlowTriples.cxx¶
↰ Return to documentation for file (larflow/PrepFlowMatchData/FlowTriples.cxx)
#include "FlowTriples.h"
#include "WireOverlap.h"
#include <ctime>
#include <sstream>
namespace larflow {
namespace prep {
FlowTriples::FlowTriples( int source, int target,
const std::vector<larcv::Image2D>& adc_v,
const std::vector<larcv::Image2D>& badch_v,
float threshold, bool save_index ) {
_sparseimg_vv = make_initial_sparse_image( adc_v, threshold );
_makeTriples( source, target, adc_v, badch_v, _sparseimg_vv, threshold, save_index );
}
FlowTriples::FlowTriples( int source, int target,
const std::vector<larcv::Image2D>& adc_v,
const std::vector<larcv::Image2D>& badch_v,
const std::vector< std::vector<PixData_t> >& sparseimg_vv,
float threshold, bool save_index ) {
_makeTriples( source, target, adc_v, badch_v, sparseimg_vv, threshold, save_index );
}
void FlowTriples::_makeTriples( int source, int target,
const std::vector<larcv::Image2D>& adc_v,
const std::vector<larcv::Image2D>& badch_v,
const std::vector< std::vector<PixData_t> >& sparseimg_vv,
float threshold, bool save_index ) {
std::clock_t start_ = std::clock();
_source_plane = source;
_target_plane = target;
_other_plane = larflow::LArFlowConstants::getOtherPlane( source, target );
// allocate space for triples
_triple_v.reserve( sparseimg_vv[_source_plane].size() );
int ndeadch_added = 0;
std::map< std::pair<int,int>, int > _added_dead_channel_pixels[3];
_deadch_to_add.resize( adc_v.size() );
// turn this on for a full dump
bool for_debug = false;
// make possible triples based on going from source pixel to target pixels
for ( auto& srcpix : sparseimg_vv[_source_plane] ) {
// activate conditional debug
// we puke for pixels above threshold, but sitting in dead channel.
// this is to trace effect of INFILL
// float badchval = badch_v[_source_plane].pixel( srcpix.row, srcpix.col );
// if ( badchval>0 )
// for_debug = true;
// else
// for_debug = false;
// we get all the 'target' and 'other' plane wires this pixel overlaps with
std::vector< std::vector<int> > overlap = larflow::prep::WireOverlap::getOverlappingWires( _source_plane, _target_plane, srcpix.col );
if ( for_debug ) {
std::cout << "FlowTriples[" << _source_plane << "," << _target_plane << "," << _other_plane << "] "
<< " srcwire=" << srcpix.col << " row=" << srcpix.row
<< " tar-overlap=[" << overlap[0].front() << "," << overlap[0].back() << "] "
<< " oth-overlap=[" << overlap[1].front() << "," << overlap[1].back() << "] "
<< std::endl;
}
if ( overlap[0].size()==0 ) {
if ( for_debug )
std::cout << " -- target no overlap" << std::endl;
continue;
}
// get iterator to overlap[0] vector: contains pixels in the target wire, the source wire overlaps with
auto it_overlap0 = overlap[0].begin();
// get the lowerbound pixel for (row,col) in target sparse image container (assumes sorted properly sort worked)
PixData_t lb( srcpix.row, (int)overlap[0][0], 0.0 );
auto it_target = std::lower_bound( sparseimg_vv[_target_plane].begin(), sparseimg_vv[_target_plane].end(), lb );
if ( for_debug ) {
if ( it_target!=sparseimg_vv[_target_plane].end() )
std::cout << " target_pix_lb(row,col)=(" << it_target->row << "," << it_target->col << ") ";
std::cout << " other=[" << overlap[1].front() << "," << overlap[1].back() << "]"
<< std::endl;
}
// iterator forward in target wire with charge, until we go to another row (different time)
if (for_debug) std::cout << " ... start of target pixel loop ..." << std::endl;
while ( it_target!=sparseimg_vv[_target_plane].end() && it_target->row==srcpix.row ) {
// break if we out of range over the target wire
if ( it_target->col>overlap[0].back() || it_target->row!=srcpix.row ) break;
// find position in overlap[0] in order to get overlap[1] element, i.e. the other wire the source+target wires intersect
it_overlap0 = std::lower_bound( it_overlap0, overlap[0].end(), it_target->col );
if ( it_overlap0==overlap[0].end() ) {
if ( for_debug ) std::cout << " target col not in list: break" << std::endl;
break; // didnt find the target column in the list
}
// scan up until matches
int ivec = -1;
while ( *it_overlap0<=it_target->col && it_overlap0!=overlap[0].end()) {
ivec = it_overlap0-overlap[0].begin();
it_overlap0++;
}
if ( for_debug )
std::cout << " .. target pixel: row=" << it_target->row << " col=" << it_target->col << " ivec=" << ivec << std::endl;
int otherwire = overlap[1][ivec];
if ( for_debug )
std::cout << " ... search for otherwire=" << otherwire << " ivec=" << ivec << std::endl;
// now find the other plane pixel in the sparse matrix.
// we allow for a little slop
// first search for pixel in other plane sparseimg vector (that is above threshold)
auto it_other = std::lower_bound( sparseimg_vv[_other_plane].begin(),
sparseimg_vv[_other_plane].end(),
PixData_t( srcpix.row, otherwire-2, 0.0) );
if ( for_debug ) {
if ( it_other!=sparseimg_vv[_other_plane].end() ) {
std::cout << " ... otherlb(r,c)=" << it_other->row << "," << it_other->col << ")" << std::endl;
}
}
// now we scan through the pixels in the other plane
bool found = false;
if ( it_other!=sparseimg_vv[_other_plane].end() && it_other->row==srcpix.row && it_other->col<=otherwire ) {
// pixels to search
while ( it_other!=sparseimg_vv[_other_plane].end() && it_other->row==srcpix.row ) {
// if a pixel in the other plane is found close to the one we searched for,
// we know it has charge, so we store the triple
//valid triple
if ( for_debug ) {
std::cout << " .... looking for row=" << srcpix.row << " otherwire=" << otherwire
<< " cols=(" << srcpix.col << "," << it_target->col << "," << it_other->col << ")"
<< std::endl;
}
if ( abs(it_other->col-otherwire)<=1 ) {
if (for_debug) {
std::cout << " ... found triple with charge!!!" << std::endl;
}
if (save_index) {
std::vector<int> trip = { srcpix.idx, it_target->idx, it_other->idx, srcpix.row }; // store positions in sparsematrix
_triple_v.push_back( trip );
}
else {
std::vector<int> trip = { srcpix.col, it_target->col, it_other->col, srcpix.row }; // store positions in sparsematrix
_triple_v.push_back( trip );
}
found = true;
break;
}
it_other++;
}
}
// if we did not find a pixel in the other plane, then we check if it is in a bad region
if ( !found && badch_v[ _other_plane ].pixel( srcpix.row, otherwire ) > 0 ) {
// badchannel, other wire lands in dead channel
// check we have this pixel
if ( for_debug ) {
std::cout << " ... looking for badch row=" << srcpix.row
<< " cols=(" << srcpix.col << "," << it_target->col << "," << otherwire << ")"
<< std::endl;
std::cout << " ... found triple with dead channel" << std::endl;
}
// store this dead pixel to add to the sparsematrix after we complete the search method
auto it_dead_other = _added_dead_channel_pixels[_other_plane].find( std::pair<int,int>(srcpix.row,otherwire) );
int otherplane_index = 0;
if ( it_dead_other==_added_dead_channel_pixels[_other_plane].end() ) {
// was not in the creation set, create this pixel in the other plane sparse image
otherplane_index = sparseimg_vv[_other_plane].size() + _deadch_to_add[_other_plane].size();
PixData_t badpix( srcpix.row, otherwire, 0.0 );
badpix.idx = otherplane_index;
// this screws up the search, but at the end, so won't matter?
_deadch_to_add[_other_plane].push_back( badpix );
_added_dead_channel_pixels[_other_plane][ std::pair<int,int>(srcpix.row,otherwire) ] = otherplane_index;
ndeadch_added++;
}
else {
otherplane_index = it_dead_other->second;
}
if (save_index) {
// store position in sparsematrix, omit last plane
std::vector<int> trip = { srcpix.idx, it_target->idx, otherplane_index, srcpix.row };
_triple_v.push_back( trip );
}
else {
// store position in image
std::vector<int> trip = { srcpix.col, it_target->col, otherwire, srcpix.row };
_triple_v.push_back( trip );
}
}
else {
if ( for_debug ) {
std::cout << "looking for row=" << srcpix.row
<< " cols=(" << srcpix.col << "," << it_target->col << "," << otherwire << ")"
<< std::endl;
if ( it_other==sparseimg_vv[_other_plane].end() )
std::cout << " ... no triple" << std::endl;
else
std::cout << " ... lowerbound found, but no match? lowerbound otherpix=(" <<it_other->row << "," << it_other->col << ")" << std::endl;
}
}
// iterate the target pixel
it_target++;
}//target while loop
}
// for ( auto& deadpix : _deadch_to_add[_other_plane] ) {
// sparseimg_vv[_other_plane].push_back( deadpix );
// }
_triple_v.shrink_to_fit();
std::clock_t end_ = std::clock();
std::cout << "[FlowTriples] for flow source[" << _source_plane << "] "
<< "to target[" << _target_plane << "] planes "
<< "found " << _triple_v.size() << " triples "
<< "ndeadch-added=" << ndeadch_added << " "
<< "elasped=" << float(end_-start_)/float(CLOCKS_PER_SEC)
<< std::endl;
}
std::vector<TH2D> FlowTriples::plot_triple_data( const std::vector<larcv::Image2D>& adc_v,
const std::vector< std::vector<PixData_t> >& sparseimg_vv,
std::string hist_stem_name ) {
std::vector<TH2D> out_v;
for ( int p=0; p<(int)adc_v.size(); p++ ) {
std::stringstream ss;
ss << "htriples_plane" << p << "_" << hist_stem_name;
auto const& meta = adc_v[p].meta();
TH2D hist( ss.str().c_str(), "",
meta.cols(), meta.min_x(), meta.max_x(),
meta.rows(), meta.min_y(), meta.max_y() );
out_v.emplace_back(std::move(hist));
}
std::vector<int> pl = { _source_plane, _target_plane, _other_plane };
for ( auto const& trip : _triple_v ) {
for (size_t i=0; i<3; i++ ) {
auto const& pix = sparseimg_vv[ pl[i] ][ trip[i] ];
out_v[ pl[i] ].SetBinContent( pix.col+1, pix.row+1, pix.val+1 );
}
}
return out_v;
}
std::vector<TH2D> FlowTriples::plot_sparse_data( const std::vector<larcv::Image2D>& adc_v,
const std::vector< std::vector<PixData_t> >& sparseimg_vv,
std::string hist_stem_name ) {
std::vector<TH2D> out_v;
for ( int p=0; p<(int)adc_v.size(); p++ ) {
std::stringstream ss;
ss << "htriples_plane" << p << "_" << hist_stem_name;
auto const& meta = adc_v[p].meta();
TH2D hist( ss.str().c_str(), "",
meta.cols(), meta.min_x(), meta.max_x(),
meta.rows(), meta.min_y(), meta.max_y() );
for ( auto const& pix : sparseimg_vv[p] ) {
if (pix.val>=10 )
hist.SetBinContent( pix.col+1, pix.row+1, pix.val );
}
out_v.emplace_back(std::move(hist));
}
return out_v;
}
std::vector< std::vector<FlowTriples::PixData_t> >
FlowTriples::make_initial_sparse_image( const std::vector<larcv::Image2D>& adc_v,
float threshold ) {
// sparsify planes: pixels must be above threshold
std::vector< std::vector<FlowTriples::PixData_t> > sparseimg_vv(adc_v.size());
for ( size_t p=0; p<adc_v.size(); p++ ) {
sparseimg_vv[p].reserve( (int)( 0.1 * adc_v[p].as_vector().size() ) );
for ( size_t r=0; r<adc_v[p].meta().rows(); r++ ) {
for ( size_t c=0; c<adc_v[p].meta().cols(); c++ ) {
float val = adc_v[p].pixel(r,c);
if ( val>=threshold ) {
sparseimg_vv[p].push_back( PixData_t((int)r,(int)c, val) );
}
}
}
// should be sorted in (r,c). do i pull the trigger and sort?
// std::sort( sparseimg_vv[p].begin(), sparseimg_vv[p].end() );
int idx=0;
for ( auto& pix : sparseimg_vv[p] ) {
pix.idx = idx;
idx++;
}
std::cout << "[FlowTriples] plane[" << p << "] has " << sparseimg_vv[p].size() << " (above threshold) pixels" << std::endl;
}
return sparseimg_vv;
}
}
}