.. _program_listing_file_larflow_Reco_ShowerLikelihoodBuilder.cxx:

Program Listing for File ShowerLikelihoodBuilder.cxx
====================================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_Reco_ShowerLikelihoodBuilder.cxx>` (``larflow/Reco/ShowerLikelihoodBuilder.cxx``)

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

.. code-block:: cpp

   #include "ShowerLikelihoodBuilder.h"
   
   #include "larcv/core/DataFormat/DataFormatTypes.h"
   #include "larcv/core/DataFormat/EventImage2D.h"
   #include "DataFormat/larflow3dhit.h"
   #include "DataFormat/mcshower.h"
   #include "LArUtil/SpaceChargeMicroBooNE.h"
   #include "ublarcvapp/MCTools/MCPixelPGraph.h"
   #include "cluster_functions.h"
   
   namespace larflow {
   namespace reco {
   
     ShowerLikelihoodBuilder::ShowerLikelihoodBuilder()
     {
       _hll          = new TH2F("lfshower_ll","",2000, -10, 190, 1000, 0, 100 );
       _hll_weighted = new TH2F("lfshower_ll_weighted","",2000, -10, 190, 1000, 0, 100 );
       _psce = new larutil::SpaceChargeMicroBooNE();
     }
   
     ShowerLikelihoodBuilder::~ShowerLikelihoodBuilder()
     {
       // if ( _hll ) delete _hll;
       // _hll = nullptr;
     }
   
     void ShowerLikelihoodBuilder::process( larcv::IOManager& iolcv, larlite::storage_manager& ioll )
     {
   
   
   
       // get input data
       larcv::EventImage2D* ev_adc      =
         (larcv::EventImage2D*)iolcv.get_data( larcv::kProductImage2D, "wiremc" );
       larcv::EventImage2D* ev_seg      =
         (larcv::EventImage2D*)iolcv.get_data( larcv::kProductImage2D, "segment" );
       larcv::EventImage2D* ev_trueflow =
         (larcv::EventImage2D*)iolcv.get_data( larcv::kProductImage2D, "larflow" );
   
       larlite::event_mcshower* ev_mcshower
         = (larlite::event_mcshower*)ioll.get_data( larlite::data::kMCShower, "mcreco" );
   
       std::vector<larcv::Image2D> masked_v;
       std::vector<larcv::Image2D> badch_v;
       for ( size_t p=0; p<3; p++ ) {
         auto const& adc = ev_adc->Image2DArray()[p];
         auto const& seg = ev_seg->Image2DArray()[p];
   
         // image with shower pixels only
         larcv::Image2D masked(adc.meta());
         masked.paint(0.0);
   
         // dummy badch, where every channel is a dummy channel
         larcv::Image2D badch(adc.meta());
         badch.paint(0.0);
         
         // find shower pixels
         for ( size_t r=0; r<adc.meta().rows(); r++) {
           for ( size_t c=0; c<adc.meta().cols(); c++) {
   
             int pid = (int)seg.pixel(r,c);
             if (pid>0) badch.set_pixel(r,c,50.0);
             
             if ( pid !=larcv::kROIEminus && pid!=larcv::kROIGamma )
               continue;
             
             // copy pixel
             if ( adc.pixel(r,c)>10 )
               masked.set_pixel( r, c, adc.pixel(r,c) );
             else
               masked.set_pixel( r, c, 15.0 );
           }
         }
         masked_v.emplace_back( std::move(masked) );
         badch_v.emplace_back( std::move(badch) );
       }
   
       // images made, now we make the triplets
       tripletalgo.process( masked_v, badch_v, 1.0, true );
       tripletalgo.make_truth_vector( ev_trueflow->Image2DArray() );
   
       // save the true triplets as larflow3dhits    
       std::vector<larlite::larflow3dhit> truehit_v;
       
       for ( size_t i=0; i<tripletalgo._triplet_v.size(); i++ ) {
   
         if ( tripletalgo._truth_v[i]!=1 )
           continue;
         
         larlite::larflow3dhit hit;
         // ==========================================
         // DEFINITION OF DATA STORED IN LARFLOW3DHIT
         // * [0-2]:   x,y,z
         // * [3-9]:   6 flow direction scores + 1 max scire (deprecated based on 2-flow paradigm. for triplet, [8] is the only score stored 
         // * [10-12]: 3 ssnet scores, (bg,track,shower)
         // * [13]:    1 keypoint label score
         // ==========================================
   
         hit.resize( 14, 0 );
         for (size_t v=0; v<3; v++ ) hit[v] = tripletalgo._pos_v[i][v];
         hit[12] = 1.0;
         hit[13] = 0.0;
         hit.tick = tripletalgo._triplet_v[i][3];
         hit.srcwire = tripletalgo._triplet_v[i][2];
         hit.targetwire.resize(3);
         
         
         for (size_t p=0; p<3; p++)
           hit.targetwire[p] = tripletalgo._sparseimg_vv[p][ tripletalgo._triplet_v[i][p] ].col;
         hit.idxhit = i;
   
         truehit_v.emplace_back( std::move(hit) );
           
       }// end of triplet loop
       std::cout << "Stored " << truehit_v.size() << " true shower hits" << std::endl;
   
       // cluster the truehits
       _make_truehit_clusters( truehit_v );
   
       // now we parse the truth, by building the particle graph
       ublarcvapp::mctools::MCPixelPGraph mcpg;
       mcpg.set_adc_treename( "wire" );
       mcpg.buildgraph( iolcv, ioll );
   
       // Get neutrino particles
       std::vector<ublarcvapp::mctools::MCPixelPGraph::Node_t*> nodes_v
         = mcpg.getNeutrinoParticles();
   
       // we loop over the different MC showers
       // we try to associate a truehit cluster to the mc shower object
       struct ShowerInfo_t {
         int idx;
         int pid;
         int origin;
         int priority;
         int matched_cluster;
         std::vector<float> shower_dir;
         std::vector<float> shower_dir_sce;
         std::vector<float> shower_vtx;
         std::vector<float> shower_vtx_sce;      
         bool operator<(const ShowerInfo_t& rhs ) {
           if ( priority<rhs.priority ) return true;
           else if ( priority==rhs.priority && idx<rhs.idx ) return true;
           return false;
         };
       };
       std::vector<ShowerInfo_t> shower_info_v;
         
       for ( size_t idx=0; idx<ev_mcshower->size(); idx++ ) {
         auto const& mcsh = ev_mcshower->at( idx );
   
         ShowerInfo_t info;
         info.idx = idx;
         info.pid = mcsh.PdgCode();
         info.origin = mcsh.Origin();
         if ( info.origin==1 && abs(info.pid)==11 )
           info.priority = 0;
         else if ( info.origin==1 && abs(info.pid)==22 )
           info.priority = 1;
         else if ( info.origin==1 )
           info.priority = 2;
         else if ( info.origin==2 && abs(info.pid)==11 )
           info.priority = 3;
         else
           info.priority = 4;
   
         info.shower_dir.resize(3,0);
         info.shower_vtx.resize(3,0);
         info.shower_dir_sce.resize(3,0);
         info.shower_vtx_sce.resize(3,0);
         auto const& detprofdir = mcsh.DetProfile().Momentum().Vect();
         for (int i=0; i<3; i++) {
           info.shower_dir[i] = detprofdir[i]/detprofdir.Mag();
           info.shower_vtx[i] = mcsh.DetProfile().Position()[i];
         }
   
         float dirnorm = 0.;
         for (int v=0; v<3; v++)
           dirnorm += info.shower_dir[v]*info.shower_dir[v];
         dirnorm = sqrt(dirnorm);
         for (int v=0; v<3; v++)
           info.shower_dir[v] /= dirnorm;
   
         // adjust shower dir due to SCE
         float cos_sce = 1.0;      
         if ( abs(info.pid)==11 ) {
           
           // SCE corrected direction
           std::vector<float> shower_end(3,0);
           for (int p=0; p<3; p++ ) shower_end[p] = info.shower_vtx[p] + 3.0*info.shower_dir[p];
           
           std::vector<double> offset = _psce->GetPosOffsets( info.shower_vtx[0], info.shower_vtx[1], info.shower_vtx[2] );
           info.shower_vtx_sce  = { info.shower_vtx[0] - (float)offset[0] + (float)0.6,
                                    info.shower_vtx[1] + (float)offset[1],
                                    info.shower_vtx[2] + (float)offset[2] };
           offset = _psce->GetPosOffsets( shower_end[0], shower_end[1], shower_end[2] );
           std::vector<float> shower_end_sce  = { shower_end[0] - (float)offset[0] + (float)0.6,
                                                  shower_end[1] + (float)offset[1],
                                                  shower_end[2] + (float)offset[2] };
           float scenorm = 0.;
           for (int i=0; i<3; i++ ) {
             info.shower_dir_sce[i] = shower_end_sce[i]-info.shower_vtx_sce[i];
             scenorm += info.shower_dir_sce[i]*info.shower_dir_sce[i];
             cos_sce += info.shower_dir[i]*info.shower_dir_sce[i];
           }
           
           scenorm = sqrt(scenorm);
           for (int i=0; i<3; i++ ) info.shower_dir_sce[i] /= scenorm;
           cos_sce /= scenorm;
         }
         else {
           for (int i=0; i<3; i++ ) {
             info.shower_dir_sce[i] = info.shower_dir[i];
             info.shower_vtx_sce[i] = info.shower_vtx[i];
           }
         }
   
         std::cout << "[ShowerLikelihoodBuilder] true shower " << std::endl;
         std::cout << " pid=" << info.pid << std::endl;
         std::cout << " dir-truth=(" << info.shower_dir[0] << "," << info.shower_dir[1] << "," << info.shower_dir[2] << ")" << std::endl;
         std::cout << " dir-sce=(" << info.shower_dir_sce[0] << "," << info.shower_dir_sce[1] << "," << info.shower_dir_sce[2] << ")" << std::endl;
         std::cout << " cos(truth*sce)=" << cos_sce << std::endl;
         std::cout << " vertex-truth=(" << info.shower_vtx[0] << "," << info.shower_vtx[1] << "," << info.shower_vtx[2] << ")" << std::endl;
         std::cout << " vertex-sce=(" << info.shower_vtx_sce[0] << "," << info.shower_vtx_sce[1] << "," << info.shower_vtx_sce[2] << ")" << std::endl;
   
         shower_info_v.push_back( info );
       }
   
       std::sort( shower_info_v.begin(), shower_info_v.end() );    
       std::cout << "[ShowerLikelihoodBuilder] saved " << shower_info_v.size() << " showers" << std::endl;
   
       std::vector<int> claimed_cluster_v( cluster_v.size(), 0 );
       std::vector<int> cluster_match_v( shower_info_v.size(), -1 );
       int iidx = 0;
       for ( auto& info : shower_info_v ) {
         std::cout << "[" << iidx << "] pid=" << info.pid
                   << " origin=" << info.origin
                   << " vtx=(" << info.shower_vtx_sce[0] << "," << info.shower_vtx_sce[1] << "," << info.shower_vtx_sce[2] << ")"
                   << std::endl;
         int match_cluster_idx = _find_closest_cluster( claimed_cluster_v, info.shower_vtx_sce, info.shower_dir_sce );
         shower_info_v[iidx].matched_cluster = match_cluster_idx;
         std::cout << "   matched cluster index: " << match_cluster_idx << std::endl;
         iidx++;        
       }
   
       /*
   
       // cluster hits into cluster_t objects, using dbscan
       _analyze_clusters( truehit_v, shower_dir_sce, shower_vtx_sce );
       */
   
       std::vector< larlite::larflow3dhit > clustered_truehits_v;
       for ( auto& cluster : cluster_v ) {
         for ( auto& hitidx : cluster.hitidx_v ) {
           clustered_truehits_v.push_back( truehit_v[hitidx] );
         }
       }
       
       // fill profile histogram
   
       // break into clusters
       // truth ID the trunk cluster
       // save vars for the trunk verus non-trunk clutsters
       // 
       // _fillProfileHist( truehit_v, shower_dir_sce, shower_vtx_sce );
   
       // mcpg.printGraph();
       
       // save the larflow hits
       larlite::event_larflow3dhit* evout = (larlite::event_larflow3dhit*)ioll.get_data(larlite::data::kLArFlow3DHit, "trueshowerhits" );
       for (auto& hit : clustered_truehits_v )
         evout->emplace_back( std::move(hit) );
   
       // save the shower object we are basing the info on
       larlite::event_mcshower* mcshowerout = (larlite::event_mcshower*)ioll.get_data(larlite::data::kMCShower, "truthshower" );
       for ( auto& info : shower_info_v ) {
         mcshowerout->push_back( ev_mcshower->at( info.idx ) );
       }
   
       // save the pca's of the clusters
       larlite::event_pcaxis*  evout_pcaxis = (larlite::event_pcaxis*)ioll.get_data(larlite::data::kPCAxis, "truthshower");
       for ( size_t cidx=0; cidx<cluster_v.size(); cidx++ ) {
         auto const& c = cluster_v[cidx];
         larlite::pcaxis llpca = cluster_make_pcaxis( c, cidx );
         evout_pcaxis->push_back( llpca );
       }
       
       // saved masked shower image
       larcv::EventImage2D* evimgout = (larcv::EventImage2D*)iolcv.get_data(larcv::kProductImage2D, "trueshoweradc" );
       evimgout->Emplace( std::move(masked_v) );
   
     }
   
     void ShowerLikelihoodBuilder::_fillProfileHist( const std::vector<larlite::larflow3dhit>& truehit_v,
                                                     std::vector<float>& shower_dir,
                                                     std::vector<float>& shower_vtx )
     {
   
       // get distance of point from pca-axis
       // http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
   
       std::cout << "Fill hits for shower[ "
                 << "dir=(" << shower_dir[0] << "," << shower_dir[1] << "," << shower_dir[2] << ") "
                 << "vtx=(" << shower_vtx[0] << "," << shower_vtx[1] << "," << shower_vtx[2] << ") "
                 << "] with nhits=" << truehit_v.size()
                 << std::endl;
   
       std::vector<float> shower_end(3);
       std::vector<float> d3(3);
       float len3 = 1.0;
       for (int i=0; i<3; i++ ) {
         shower_end[i] = shower_vtx[i] + shower_dir[i];
         d3[i] = shower_dir[i];
       }
         
       for ( auto const& hit : truehit_v ) {
         
         std::vector<float> pt = { hit[0], hit[1], hit[2] };
   
         std::vector<float> d1(3);
         std::vector<float> d2(3);
   
         float len1 = 0.;
         for (int i=0; i<3; i++ ) {
           d1[i] = pt[i] - shower_vtx[i];
           d2[i] = pt[i] - shower_end[i];
           len1 += d1[i]*d1[i];
         }
         len1 = sqrt(len1);
   
         // cross-product
         std::vector<float> d1xd2(3);
         d1xd2[0] =  d1[1]*d2[2] - d1[2]*d2[1];
         d1xd2[1] = -d1[0]*d2[2] + d1[2]*d2[0];
         d1xd2[2] =  d1[0]*d2[1] - d1[1]*d2[0];
         float len1x2 = 0.;
         for ( int i=0; i<3; i++ ) {
           len1x2 += d1xd2[i]*d1xd2[i];
         }
         len1x2 = sqrt(len1x2);
         float rad  = len1x2/len3; // distance of point from PCA-axis
   
         float proj = 0.;
         for ( int i=0; i<3; i++ )
           proj += shower_dir[i]*d1[i];
   
         // std::cout << "  hit: (" << pt[0] <<  "," << pt[1] << "," << pt[2] << ") "
         //           << " dist=" << len1
         //           << " proj=" << proj
         //           << " rad=" << rad
         //           << std::endl;
         float w=1.;
         if ( rad>0 ) w = 1.0/rad;
         _hll_weighted->Fill( proj, rad, w );
         _hll->Fill( proj, rad );
       }
       
     }
   
     void ShowerLikelihoodBuilder::_make_truehit_clusters( std::vector< larlite::larflow3dhit >& truehit_v )
     {
   
       cluster_v.clear();
       float maxdist = 1.0;
       int minsize = 50;
       int maxkd = 5;
       cluster_larflow3dhits( truehit_v, cluster_v, maxdist, minsize, maxkd );
       for ( auto& cluster : cluster_v )
         cluster_pca( cluster );
     }
   
     int ShowerLikelihoodBuilder::_find_closest_cluster( std::vector<int>& claimed_cluster_v,
                                                         std::vector<float>& shower_vtx,
                                                         std::vector<float>& shower_dir )
     {
   
       // we define the trunk of the cluster as the one closest to the shower start
   
       float min_dist2vtx = 1.0e9;
       std::vector<float> trunk_endpt(3,0);
       int best_matched_cluster = -1;
       for ( size_t idx=0; idx<cluster_v.size(); idx++ ) {
   
         auto& cluster = cluster_v[idx];
         float dist2vtx[2] = {0.};
         for (int e=0; e<2; e++) {
           dist2vtx[e] = 0.;
           for (int i=0; i<3; i++) {
             dist2vtx[e] += (cluster.pca_ends_v[e][i]-shower_vtx[i])*(cluster.pca_ends_v[e][i]-shower_vtx[i]);
           }
           
           if ( dist2vtx[e]<min_dist2vtx ) {
             best_matched_cluster = idx;
             min_dist2vtx = dist2vtx[e];
             trunk_endpt  = cluster.pca_ends_v[e];
           }
         }
         float mindist2vtx = (dist2vtx[0]<dist2vtx[1]) ? dist2vtx[0] : dist2vtx[1];
         // std::cout << "[cluster " << idx << "] mindist2vtx=" << mindist2vtx
         //           << " start=(" << cluster.pca_ends_v[0][0] << "," << cluster.pca_ends_v[0][1] << "," << cluster.pca_ends_v[0][2] << ") "
         //           << " end=("   << cluster.pca_ends_v[1][0] << "," << cluster.pca_ends_v[1][1] << "," << cluster.pca_ends_v[1][2] << ") "        
         //           << std::endl;
       }
   
       std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] "
                 << "trunk cluster index=" << best_matched_cluster << " "
                 << "dist2vertex=" << min_dist2vtx
                 << std::endl;
   
       if ( min_dist2vtx>3.0 )
         return -1;
   
       claimed_cluster_v[ best_matched_cluster ] += 1;
       return best_matched_cluster;
     }  
     
     void ShowerLikelihoodBuilder::_analyze_clusters( std::vector< larlite::larflow3dhit >& truehit_v,
                                                      std::vector<float>& shower_dir,
                                                      std::vector<float>& shower_vtx )
     {
   
       cluster_v.clear();
       _trunk_cluster = -1;
       
       float maxdist = 1.0;
       int minsize = 50;
       int maxkd = 5;
       cluster_larflow3dhits( truehit_v, cluster_v, maxdist, minsize, maxkd );
   
       std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] number of clusters: " << cluster_v.size() << std::endl;
   
       // we define the trunk of the cluster as the one closest to the shower start
   
       float min_dist2vtx = 1.0e9;
       std::vector<float> trunk_endpt(3,0);
       for ( size_t idx=0; idx<cluster_v.size(); idx++ ) {
         auto& cluster = cluster_v[idx];
         std::cout << " [" << idx << "] analyze cluster size=" << cluster.points_v.size() << std::endl;      
         cluster_pca( cluster );
         float dist2vtx[2] = {0.};
         for (int e=0; e<2; e++) {
           for (int i=0; i<3; i++) {
             dist2vtx[e] += (cluster.pca_ends_v[e][i]-shower_vtx[i])*(cluster.pca_ends_v[e][i]-shower_vtx[i]);
           }
   
           if ( dist2vtx[e]<min_dist2vtx ) {
             _trunk_cluster = idx;
             min_dist2vtx = dist2vtx[e];
             trunk_endpt = cluster.pca_ends_v[e];
           }
         }
       }
   
       std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] "
                 << "trunk cluster index=" << _trunk_cluster << " "
                 << "dist2vertex=" << min_dist2vtx
                 << std::endl;
   
       if ( min_dist2vtx>3.0 )
         return;
       
       // get points near the vertex of the trunk cluster, within 5 cm
       cluster_t near_vtx;
       for ( int idx=0; idx<(int)cluster_v[_trunk_cluster].points_v.size(); idx++ ) {
         float dist = 0.;
         for (int i=0; i<3; i++) {
           dist += (cluster_v[_trunk_cluster].points_v[idx][i]-shower_vtx[i])*(cluster_v[_trunk_cluster].points_v[idx][i]-shower_vtx[i]);
         }
         dist = sqrt(dist);
         if ( dist<5.0 ) {
           near_vtx.points_v.push_back( cluster_v[_trunk_cluster].points_v[idx] );
         }
       }
       cluster_pca( near_vtx );
       std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] "
                 << "size of trunk portion: " << near_vtx.points_v.size() << " spacepoints" << std::endl;
   
       // calculate cosine between trunk pca and trunk start direction
       float trunkpcacos = 0.;
       for (int i=0; i<3; i++ ) {
         trunkpcacos += shower_dir[i]*near_vtx.pca_axis_v[0][i];
       }
       std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] "
                 << " cos(shower_dir * firstpca)=" << trunkpcacos
                 << std::endl;
   
       // now we measure relationships to the main cluster
       // (1) distance of cluster-endpoint to trunk pca-axis of nearest endpt
       // (2) cosine of pca-axes
       // (3) impact parameter: smallest distance between trunk and cluster pca-axis
   
       // loop over each cluster
       int idx=-1;
       for ( auto& clust : cluster_v ) {
         idx++;
         // skip trunk cluster
         if ( idx==_trunk_cluster ) continue;
         std::cout << "[ ShowerLikelihoodBuilder::_analyze_clusters ] compare trunk with cluster[" << idx << "]" << std::endl;
   
         // [2: cos(pca)]
         float cospca = 0.;
         for (int i=0; i<3; i++) {
           cospca += near_vtx.pca_axis_v[0][i]*clust.pca_axis_v[0][i];
         }
         std::cout << "  cluster-trunk pca-cosine: " << cospca << std::endl;
   
         // [3: impact param]
         float impact_dist, proj_trunk, proj_clust;
         _impactdist( trunk_endpt, near_vtx.pca_axis_v[0],
                      clust.pca_center, clust.pca_axis_v[0],
                      impact_dist, proj_trunk, proj_clust );
         std::cout << "  impact dist: " << impact_dist << std::endl;
         std::cout << "  impact pt trunk proj: " << proj_trunk << std::endl;
         std::vector<float> impactpt(3,0);
         for (int i=0; i<3; i++) {
           impactpt[i] = trunk_endpt[i] + proj_trunk*near_vtx.pca_axis_v[0][i];
         }
   
         // [1: distance]
         // we use distance between impact point on trunk and endpoints
         float enddist[2] = { 0, 0 };
         for (int e=0; e<2; e++) {      
           for(int i=0; i<3; i++) {
             enddist[e] += (impactpt[i]-clust.pca_ends_v[e][i])*(impactpt[i]-clust.pca_ends_v[e][i]);
           }
           enddist[e] = sqrt(enddist[e]);
         }
         float clust_dist = (enddist[0]<enddist[1]) ? enddist[0] : enddist[1];
         std::cout << "  dist of cluster to trunk: " << clust_dist << std::endl;
         
       }
       
     }
   
     void ShowerLikelihoodBuilder::_dist2line( const std::vector<float>& ray_start,
                                               const std::vector<float>& ray_dir,
                                               const std::vector<float>& pt,
                                               float& radial_dist, float& projection )
     {
   
       std::vector<float> d1(3);
       std::vector<float> d2(3);
   
       float len3 = 0.;
       std::vector<float> end(3,0);
       for (int i=0; i<3; i++) {
         end[i] = ray_start[i] + ray_dir[i];
         len3 += ray_dir[i]*ray_dir[i];
       }
       len3 = sqrt(len3);
   
       float len1 = 0.;
       for (int i=0; i<3; i++ ) {
         d1[i] = pt[i] - ray_start[i];
         d2[i] = pt[i] - end[i];
         len1 += d1[i]*d1[i];
       }
       len1 = sqrt(len1);
   
       // cross-product
       std::vector<float> d1xd2(3);
       d1xd2[0] =  d1[1]*d2[2] - d1[2]*d2[1];
       d1xd2[1] = -d1[0]*d2[2] + d1[2]*d2[0];
       d1xd2[2] =  d1[0]*d2[1] - d1[1]*d2[0];
       float len1x2 = 0.;
       for ( int i=0; i<3; i++ ) {
         len1x2 += d1xd2[i]*d1xd2[i];
       }
       len1x2 = sqrt(len1x2);
       radial_dist  = len1x2/len3; // distance of point from PCA-axis
       
       projection = 0.;
       for ( int i=0; i<3; i++ )
         projection += ray_dir[i]*d1[i]/len3;
       
     }
   
     void ShowerLikelihoodBuilder::_impactdist( const std::vector<float>& l_start,
                                                const std::vector<float>& l_dir,
                                                const std::vector<float>& m_start,
                                                const std::vector<float>& m_dir,
                                                float& impact_dist,
                                                float& proj_l,
                                                float& proj_m )
     {
   
       impact_dist = -1.0;
       proj_l  = 0.;
       proj_m  = 0.;
       
       // L = a + bs where a,b in R3, s in R
       // M = c + dt where c,d in R3, t in R
       
       float b2 = 0.; // l dir squared
       float d2 = 0.; // m dir squred
       float bd = 0.; // innerproduct l and m
       std::vector<float> e(3,0.); // line segment of start points m->l
       float be = 0.; // inner product l_dir and e
       float de = 0.; // inner product m_dir and e
       for (int i=0; i<3; i++) {
         b2 += l_dir[i]*l_dir[i];
         d2 += m_dir[i]*m_dir[i];
         bd += l_dir[i]*m_dir[i];
         e[i] = l_start[i]-m_start[i];
         be += e[i]*l_dir[i];
         de += e[i]*m_dir[i];
       }
   
       float A = -1.0*(b2*d2 - bd*bd);
       if (A==0) { return; } // parallel lines
   
       float s = (-1.0*b2*de + be*bd)/A;
       float t = (d2*be - be*bd)/A;
   
       float D = 0.; // distance
       for (int i=0; i<3; i++ ) {
         float dd = e[i] + l_dir[i]*t - m_dir[i]*s;
         D += dd*dd;
       }
       D = sqrt(D);
       impact_dist = D;
   
       proj_l = s;
       proj_m = t;
       
     }
   
   }
   }