.. _program_listing_file_larflow_Ana_keypoint_recoana.cxx:

Program Listing for File keypoint_recoana.cxx
=============================================

|exhale_lsh| :ref:`Return to documentation for file <file_larflow_Ana_keypoint_recoana.cxx>` (``larflow/Ana/keypoint_recoana.cxx``)

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

.. code-block:: cpp

   #include <iostream>
   #include <string>
   
   #include "TTree.h"
   
   #include "larcv/core/DataFormat/IOManager.h"
   #include "larcv/core/DataFormat/EventPGraph.h"
   #include "DataFormat/storage_manager.h"
   #include "DataFormat/larflow3dhit.h"
   #include "LArUtil/LArProperties.h"
   #include "LArUtil/Geometry.h"
   #include "ublarcvapp/MCTools/LArbysMC.h"
   #include "larflow/KeyPoints/PrepKeypointData.h"
   #include "larflow/Reco/KPCluster.h"
   
   int main( int nargs, char** argv )
   {
   
     // keypoint reco analysis
     // output:
     //  1) ttree per event storing
     //     - mc event info
     //     - number of reco vertices near true neutrino vertex
     //     - number of reco vertices near true track keypoint
     //     - number of reco vertices that are false
     //  2) if wire-cell cosmic tagger info available, ttree per event storing
     //     - number of reco vertices near true neutrino vertex w/ wirecell mask
     //     - number of reco vertices near true track keypoint w/ wirecell mask
     //     - number of reco vertices that are false w/ wirecell mask
   
     // inputs
     // ------
   
     // 1) dlmerged file
     // 2) keypoint reco output
   
     std::string dlmerged_file = argv[1];
     std::string kpreco_file   = argv[2];
   
     // outputs
     // 1) ana tfile with ana ttree
     std::string outfilename   = argv[3];
     
     // Load inputs
     larcv::IOManager iocv( larcv::IOManager::kREAD, "larcv", larcv::IOManager::kTickBackward );
     iocv.add_in_file( dlmerged_file );
     iocv.reverse_all_products();
     iocv.initialize();
   
     larlite::storage_manager ioll( larlite::storage_manager::kREAD );
     ioll.add_in_filename( dlmerged_file );
     ioll.open();
   
     TFile kpreco_tfile(kpreco_file.c_str(),"open");
     TTree* kpreco_ttree = (TTree*)kpreco_tfile.Get("larflow_keypointreco");
     std::vector<larflow::reco::KPCluster>* ev_kpreco = 0;
     kpreco_ttree->SetBranchAddress( "kpcluster_v", &ev_kpreco );
   
     // load outputs
     TFile* out = new TFile(outfilename.c_str(),"new");
   
     // define output tree
     TTree* ev_ana = new TTree("kprecoana_event","Keypoint Reco Ana Event tree");
   
     int run,subrun,event;
     
     int n_reco_kp;
     int n_near_true_vtx;
     int n_near_true_trackend;
     int n_near_true_showerstart;
     int n_false_keypoint;
   
     int n_reco_kp_wct;  
     int n_near_true_vtx_wct;
     int n_near_true_trackend_wct;
     int n_near_true_showerstart_wct;
     int n_false_keypoint_wct;
   
     float vtx_qsum[3] = { 0., 0., 0. };
     float min_dist_to_vtx = 0.;
     float min_dist_to_vtx_dl = 0.;
     int has_good_dl_vertex = 0;
     
     ev_ana->Branch("run",&run,"run/I");
     ev_ana->Branch("subrun",&subrun,"subrun/I");
     ev_ana->Branch("event",&event,"event/I");  
   
     ev_ana->Branch("n_reco_kp",&n_reco_kp,"n_reco_kp/I");
     ev_ana->Branch("n_near_true_vtx",&n_near_true_vtx,"n_near_true_vtx/I");
     ev_ana->Branch("n_near_true_trackend",&n_near_true_trackend,"n_near_true_trackend/I");
     ev_ana->Branch("n_near_true_showerstart",&n_near_true_showerstart,"n_near_true_showerstart/I");
     ev_ana->Branch("n_false_keypoint",&n_false_keypoint,"n_false_keypoint/I");
   
     ev_ana->Branch("n_reco_kp_wct",&n_reco_kp_wct,"n_reco_kp_wct/I");  
     ev_ana->Branch("n_near_true_vtx_wct",&n_near_true_vtx_wct,"n_near_true_vtx_wct/I");
     ev_ana->Branch("n_near_true_trackend_wct",&n_near_true_trackend_wct,"n_near_true_trackend_wct/I");
     ev_ana->Branch("n_near_true_showerstart_wct",&n_near_true_showerstart_wct,"n_near_true_showerstart_wct/I");
     ev_ana->Branch("n_false_keypoint_wct",&n_false_keypoint_wct,"n_false_keypoint_wct/I");
   
     ev_ana->Branch("min_dist_to_vtx", &min_dist_to_vtx, "min_dist_to_vtx/F" );
     ev_ana->Branch("min_dist_to_vtx_dl", &min_dist_to_vtx_dl, "min_dist_to_vtx_dl/F" );
     ev_ana->Branch("has_good_dl_vertex", &has_good_dl_vertex, "has_good_dl_vertex/I" );
     ev_ana->Branch("vtx_qsum",vtx_qsum, "vtx_qsum[3]/F" );
   
     TTree* kp_ana = new TTree("kprecoana_keypt","Keypoint Reco Ana per True Keypoint tree");  
     float vtx_sce[3];
     float max_score;
     float max_score_dist;
     int is_nu_vtx;
     int n_nearby_5cm;
     kp_ana->Branch("vtx_sce",vtx_sce,"vtx_sce[3]/F");
     kp_ana->Branch("max_score",&max_score,"max_score/F");
     kp_ana->Branch("max_score_dist",&max_score_dist,"max_score_dist/F");
     kp_ana->Branch("n_nearby_5cm",&n_nearby_5cm,"n_nearby_5cm/I");
     kp_ana->Branch("is_nu_vtx",&is_nu_vtx,"is_nu_vtx/I");
   
     
     const float cut_off_dist = 10.0; // cm
     const float cut_off_distsq = cut_off_dist*cut_off_dist;
   
     // Keypoint Truth Data Maker
     larflow::keypoints::PrepKeypointData kpdata;
     ublarcvapp::mctools::LArbysMC lmc;
     lmc.bindAnaVariables( ev_ana );
     
     int nentries = iocv.get_n_entries();
     for (int ientry=0; ientry<nentries; ientry++) {
   
       // load trees
       iocv.read_entry(ientry);
       ioll.go_to(ientry);
       kpreco_ttree->GetEntry(ientry);
   
       // Get neutrino interaction truth
       lmc.process( ioll );
       float true_vtx[3] = { lmc._vtx_sce_x, lmc._vtx_sce_y, lmc._vtx_sce_z };
   
       // truth keypoints
       kpdata.process( iocv, ioll );
       
       // wirecell cosmic pixel image
       larcv::EventImage2D* ev_thrumu =
         (larcv::EventImage2D*)iocv.get_data(larcv::kProductImage2D,"thrumu");
       auto const& thrumu_v = ev_thrumu->Image2DArray();
   
       // ADC images
       larcv::EventImage2D* ev_adc =
         (larcv::EventImage2D*)iocv.get_data(larcv::kProductImage2D,"wire");
       auto const& adc_v = ev_adc->Image2DArray();
   
       // DL Vertex
       larcv::EventPGraph* ev_pgraph =
         (larcv::EventPGraph*)iocv.get_data(larcv::kProductPGraph,"test");
       
       std::cout << "[ENTRY " << ientry << "]" << std::endl;
       std::cout << "  number of truth keypoints: " << kpdata.getKPdata().size() << std::endl;
       std::cout << "  number of reco keypoints: " << ev_kpreco->size() << std::endl;
       std::cout << "  number of Wirecell images: " << thrumu_v.size() << std::endl;
       std::cout << "  true neutrino vertex: (" << true_vtx[0] << "," << true_vtx[1] << "," << true_vtx[2] << ")" << std::endl;
       std::cout << "  number of DL reco vertices: " << ev_pgraph->PGraphArray().size() << std::endl;
       lmc.printInteractionInfo();
       kpdata.printKeypoints();
       
       n_reco_kp = (int)ev_kpreco->size();
       
       // KP Reco loop
       n_near_true_vtx = 0;
       n_near_true_trackend = 0;
       n_near_true_showerstart = 0;
       n_false_keypoint = 0;
   
       n_reco_kp_wct = 0;
       n_near_true_vtx_wct = 0;
       n_near_true_trackend_wct = 0;
       n_near_true_showerstart_wct = 0;
       n_false_keypoint_wct = 0;
   
       // characterize vertex
       min_dist_to_vtx = 1.0e9;
   
       // get vertex activity
       for (int p=0; p<3; p++) vtx_qsum[p] = 0.0;
       if ( lmc._vtx_tick>adc_v[0].meta().min_y() && lmc._vtx_tick<adc_v[0].meta().max_y() ) {
         int row = adc_v[0].meta().row( lmc._vtx_tick );
   
         for (int p=0; p<3; p++ ) {
   
           if ( lmc._vtx_wire[p]>=adc_v[p].meta().min_x() && lmc._vtx_wire[p]<(int)adc_v[p].meta().max_x() ) {
             int col = adc_v[p].meta().col( lmc._vtx_wire[p] );
   
             for (int dr=-3; dr<=3; dr++ ) {
               int r=row+dr;
               if ( r<0 || r>=(int)adc_v[p].meta().rows() ) continue;
               for (int dc=-3; dc<=3; dc++) {
                 int c = col+dc;
                 if (c<0 || c>=(int)adc_v[p].meta().cols() ) continue;
                 if ( adc_v[p].pixel(r,c)>10.0 )
                   vtx_qsum[p] += adc_v[p].pixel(r,c);
               }//end of col loop
             }//end of row loop
           }//end of if valid wire
         }//end of plane loop
       }//end of if valid tick
   
       // get dl vertex info
       has_good_dl_vertex = 0;
       min_dist_to_vtx_dl = 1.0e9;
       for ( auto const& pgraph : ev_pgraph->PGraphArray() ) {
         for ( auto const& roi : pgraph.ParticleArray() ) {
           float dist = 0.;
           dist += ( roi.X()-true_vtx[0] )*( roi.X()-true_vtx[0] );
           dist += ( roi.Y()-true_vtx[1] )*( roi.Y()-true_vtx[1] );
           dist += ( roi.Z()-true_vtx[2] )*( roi.Z()-true_vtx[2] );
           dist = sqrt(dist);
           std::cout << "dl vtx: (" << roi.X() << "," << roi.Y() << "," << roi.Z() << ") dist2vtx=" << dist << " cm" << std::endl;
           if ( min_dist_to_vtx_dl>dist ) {
             min_dist_to_vtx_dl = dist;
           }
         }
       }
       if ( min_dist_to_vtx_dl<5.0 )
         has_good_dl_vertex = 1;
       
       for ( auto const& kpc : *ev_kpreco ) {
   
         // determine if on Wirecell-tagged cosmic
         // requirement: be on tagged pixel for 2/3 planes
         bool iscosmic = false;
   
         float tick = 3200 + kpc.max_pt_v[0]/larutil::LArProperties::GetME()->DriftVelocity()/0.5;
         if ( tick<=thrumu_v[0].meta().min_y() || tick>=thrumu_v[0].meta().max_y() ) {
           iscosmic = true;
         }
         else {
           int row = thrumu_v[0].meta().row( tick );        
           int nplanes_on_wctpixel = 0;
           for ( auto const& img : thrumu_v ) {
             // out of image?  we remove it implicitly by adding to counter
             std::vector<double> dpos = { (double)kpc.max_pt_v[0],
                                          (double)kpc.max_pt_v[1],
                                          (double)kpc.max_pt_v[2] };
             int nearestwire = larutil::Geometry::GetME()->NearestWire( dpos, img.meta().plane() );
             int col = img.meta().col((float)nearestwire);
   
             bool onpixel = false;
             for (int dr=-1; dr<=1; dr++) {
               int r = row+dr;
               if ( r<0 || r>=(int)img.meta().rows() ) continue;
               for (int dc=-1; dc<=1; dc++) {
                 int c = col+dc;
                 if ( c<0 || c>=(int)img.meta().cols() ) continue;
                 if ( img.pixel(r,c)>10.0 ) onpixel = true;
                 if ( onpixel ) break;
               }
               if ( onpixel ) break;            
             }
   
             if ( onpixel ) nplanes_on_wctpixel++;
           }
           if ( nplanes_on_wctpixel>=3 ) iscosmic = true;
         }
   
         // dist to vtx
         float reco_vtx_dist = 0.;
         for (int i=0; i<3; i++) {
           reco_vtx_dist += (kpc.max_pt_v[i]-true_vtx[i])*(kpc.max_pt_v[i]-true_vtx[i]);
         }
         reco_vtx_dist = sqrt(reco_vtx_dist);
         if ( reco_vtx_dist < min_dist_to_vtx )
           min_dist_to_vtx = reco_vtx_dist;
   
         
         // loop over truth keypoints and determine if near by
         int num_nearby = 0;
         float min_dist = 1.0e9;
         for ( auto const& kpd : kpdata.getKPdata() ) {
           
           float dist = 0.;
           for (int i=0; i<3; i++) {
             dist += (kpc.max_pt_v[i]-kpd.keypt[i])*(kpc.max_pt_v[i]-kpd.keypt[i]);
           }
           dist = sqrt(dist);
   
           if ( min_dist>dist )
             min_dist = dist;
   
           if ( dist>5.0 ) {
             continue;
           }
   
           // is nearby
           num_nearby++;
   
           // is neutrino vertex?
           float vtx_dist = 0.;
           for (int i=0; i<3; i++) {
             vtx_dist += (kpd.keypt[i]-true_vtx[i])*(kpd.keypt[i]-true_vtx[i]);
           }
           vtx_dist = sqrt(vtx_dist);
           if ( vtx_dist<1.0 ) {
             n_near_true_vtx++;
             if ( !iscosmic ) n_near_true_vtx_wct++;
             continue;
           }
   
           // is shower start?
           if ( kpd.is_shower==1 ) {
             n_near_true_showerstart++;
             if ( !iscosmic ) n_near_true_showerstart++;
             continue;
           }
   
           // is near track end
           n_near_true_trackend++;
           if ( !iscosmic ) n_near_true_trackend_wct++;
   
           
         }//end of loop over true keypoints
   
         //std::cout << "[reco keypoint]  closest true keypoint: " << min_dist << " cm" << std::endl;
         
         if ( num_nearby==0 ) {
           n_false_keypoint++;
           if ( !iscosmic ) n_false_keypoint_wct++;
         }
         
       }//end of loop over reco keypoints
   
       std::cout << "  nearby true nu vtx: " << n_near_true_vtx << " w/ wire-cell tagger: " << n_near_true_vtx_wct << std::endl;
       std::cout << "  nearby true shower start: " << n_near_true_showerstart << " w/ wire-cell tagger: " << n_near_true_showerstart_wct << std::endl;
       std::cout << "  nearby true track end: " << n_near_true_trackend << " w/ wire-cell tagger: " << n_near_true_trackend_wct << std::endl;
       std::cout << "  false keypoint: " << n_false_keypoint << " w/ wire-cell tagger: " << n_false_keypoint_wct << std::endl;
       
       // for ( auto const& kpd : kpdata.getKPdata() ) {
       //   // for each truth keypoint, we save max triplet within certain distance
       //   max_score_dist = 1.0e9;
       //   max_score = -1.0;
       //   float dist = 0.;
       //   for ( auto const& hit : *ev_lmhit ) {
       //     dist = 0.;
       //     for (int i=0; i<3; i++) {
       //       dist += (hit[i]-kpd.keypt[i])*(hit[i]-kpd.keypt[i]);
       //       vtx_sce[i] = kpd.keypt[i];
       //     }
       //     if ( dist<cut_off_distsq ) {
       //       if ( hit[13]>max_score ) {
       //         max_score = hit[13];
       //         max_score_dist = sqrt(dist);
       //       }
       //     }
       //   }
   
       ev_ana->Fill();
         
   
     }//end of event loop
   
     out->Write();
   
     return 0;
   }