.. _program_listing_file_larflow_PrepFlowMatchData_WireOverlap.cxx:

Program Listing for File WireOverlap.cxx
========================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_PrepFlowMatchData_WireOverlap.cxx>` (``larflow/PrepFlowMatchData/WireOverlap.cxx``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #include "WireOverlap.h"
   
   #include <iostream>
   #include <ctime>
   
   #include "larlite/core/LArUtil/Geometry.h"
   #include "larflow/LArFlowConstants/LArFlowConstants.h"
   
   namespace larflow {
   namespace prep {
   
     bool WireOverlap::_isbuilt = false;
     std::vector< std::vector<int> >   WireOverlap::_wire_targetoverlap[6];
     std::vector< std::vector<int> >   WireOverlap::_wire_otheroverlap[6];
     std::map< std::pair<int,int>, int > WireOverlap::_planeflow2mapindex;
   
     /*
      * constructor
      */
     WireOverlap::WireOverlap() {};
   
     void WireOverlap::_build() {
   
       if (_isbuilt) return;
      
       std::clock_t start = std::clock();
       
       const larutil::Geometry* geo = larutil::Geometry::GetME();
   
       float dt_intersect = 0.;
       
       // build wire overlap data
       for ( int flowdir=0; flowdir<(int)larflow::kNumFlows; flowdir++ ) {
   
         int srcplane, tarplane;
         larflow::LArFlowConstants::getFlowPlanes((larflow::FlowDir_t)flowdir,srcplane,tarplane);
         int otherplane = larflow::LArFlowConstants::getOtherPlane(srcplane,tarplane);
   
         std::cout << "[WireOverlap::build] constructing overlap maps for flow[" << flowdir << "]: "
                   << "planes [" << srcplane << "] -> [" << tarplane << ", " << otherplane << "] "
                   << std::endl;
           
                   
         _wire_targetoverlap[flowdir].resize( geo->Nwires(srcplane) );        
         _wire_otheroverlap[flowdir].resize(  geo->Nwires(srcplane) );
           
         for (int isrc=0; isrc<(int)geo->Nwires(srcplane); isrc++) {
   
           Double_t xyzstart[3];
           Double_t xyzend[3];
           geo->WireEndPoints( (UChar_t)srcplane, (UInt_t)isrc, xyzstart, xyzend );
   
           float u1 = geo->WireCoordinate( xyzstart, (UInt_t)tarplane );
           float u2 = geo->WireCoordinate( xyzend,   (UInt_t)tarplane );
   
           float umin = (u1<u2) ? u1 : u2;
           float umax = (u1>u2) ? u1 : u2;
   
           umin -= 5.0;
           umax += 5.0;
   
           if ( umin<0 ) umin = 0;
           if ( umax<0 ) umax = 0;
   
           if ( (int)umin>=geo->Nwires(tarplane) ) umin = (float)geo->Nwires(tarplane)-1;
           if ( (int)umax>=geo->Nwires(tarplane) ) umax = (float)geo->Nwires(tarplane)-1;
           
           int nwires = umax-umin+1;
           _wire_targetoverlap[flowdir][isrc].reserve( nwires );
           _wire_otheroverlap[flowdir][isrc].reserve( nwires );
   
           for (int i=0; i<nwires; i++) {
             int tarwire = (int)umin+i;
             UInt_t src_ch = geo->PlaneWireToChannel( (UInt_t)srcplane, (UInt_t)isrc );
             UInt_t tar_ch = geo->PlaneWireToChannel( (UInt_t)tarplane, (UInt_t)tarwire );
             
             
             Double_t y,z;
             //std::clock_t xs_start = std::clock();
             bool crosses  = geo->ChannelsIntersect( src_ch, tar_ch, y, z );
             //dt_intersect += float(std::clock()-xs_start)/float(CLOCKS_PER_SEC);
             
             if ( crosses ) {
               Double_t pos[3] = { 0, y, z };
               int otherwire = geo->WireCoordinate( pos, otherplane );
   
               // for debug
               // if ( srcplane==1 && tarplane==2 )
               //   std::cout << "[" << srcplane << "," << tarplane << "," << otherplane << "] "
               //             << "wire=(" << xyzstart[1] << "," << xyzstart[2] << ")->(" << xyzend[1] << "," << xyzend[2] << ") "
               //             << "channels intersecting: (" << isrc << "," << tarwire << ") "
               //             << "@ (y,z)=(" << y << "," << z << ") -> otherwire="
               //             << otherwire << std::endl;
               
               if ( otherwire>=0 && otherwire<geo->Nwires(otherplane) ) {
                 _wire_targetoverlap[flowdir][isrc].push_back( tarwire );                    
                 _wire_otheroverlap[flowdir][isrc].push_back(  otherwire );
               }
               
   
             }
   
           }//loop over scan range
           
         }//loop over source wires
   
         _isbuilt = true;
       }//end of loop over source planes
       std::clock_t end = std::clock();
       std::cout << "[wireOverlap] elapsed " << float(end-start)/float(CLOCKS_PER_SEC) << " secs "
                 << " (intersect=" << dt_intersect << " secs)"
                 << std::endl;
       //std::cin.get();
     }
   
     std::vector< std::vector<int> > WireOverlap::getOverlappingWires( int sourceplane, int targetplane, int source_wire ) {
       
       if ( !_isbuilt ) _build();
       
       //int iflow = _planeflow2mapindex[ std::pair<int,int>(sourceplane,targetplane) ];
       int iflow = (int)larflow::LArFlowConstants::getFlowDirection(sourceplane,targetplane);
       std::vector< std::vector<int> > overlap(2);
       overlap[0] = _wire_targetoverlap[iflow][source_wire];
       overlap[1] = _wire_otheroverlap[iflow][source_wire];
       return overlap;
       
     }
   
     
   }
   }