PCAngleKinkFinder_module.cc File Reference

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include <memory>
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/GeometryCore.h"
#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "sbnobj/Common/Reco/PCAnglePlane.h"
#include "sbnobj/Common/Reco/PCAngleKink.h"
#include "TFitter.h"
#include "TMinuit.h"

Go to the source code of this file.

Classes
class	sbn::PCAngleKinkFinder
struct	sbn::PCAngleKinkFinder::FitResult

Namespaces
	sbn
	This module creates Common Analysis Files.

Functions
sbn::PCAngleKink	BuildKink (const sbn::PCAngle &max, const sbn::PCAngle &lo, const sbn::PCAngle &hi, const sbn::PCAngleKinkFinder::FitResult &fit, float est_angle)
std::tuple< std::vector < sbn::PCAngle >, unsigned >	FlattenBranch (const sbn::PCAnglePlane &plane, unsigned i_branch, bool allow_incomplete, float min_dist)
float	EstimateHeight (const std::vector< sbn::PCAngle > &angles, unsigned lo, unsigned hi)
int	FindMaxIndex (const std::vector< sbn::PCAngle > &angles, unsigned lo, unsigned hi)
int	FindFWHMHiIndex (const std::vector< sbn::PCAngle > &angles, unsigned start, unsigned end, float amp, float threshold)
int	FindFWHMLoIndex (const std::vector< sbn::PCAngle > &angles, unsigned start, unsigned end, float amp, float threshold)
double	ExpectedCosAngle (double dist_to_center, double radius, double mean_distance, double angle)
void	FitLikelihood (int &npar, double *g, double &result, double *par, int flag)

Variables
static std::vector< float >	fFitDataX
static std::vector< float >	fFitDataY
static unsigned	fFitCenter

Function Documentation

sbn::PCAngleKink BuildKink	(	const sbn::PCAngle &	max,
		const sbn::PCAngle &	lo,
		const sbn::PCAngle &	hi,
		const sbn::PCAngleKinkFinder::FitResult &	fit,
		float	est_angle
	)

Definition at line 91 of file PCAngleKinkFinder_module.cc.

                                                                                                                                                       {
  sbn::PCAngleKink kink;
  kink.maxWire = max.wire;
  
  kink.fwhm_distance = sqrt((hi.hit_time_cm - lo.hit_time_cm) * (hi.hit_time_cm - lo.hit_time_cm) 
      + (hi.hit_wire_cm - lo.hit_wire_cm) * (hi.hit_wire_cm - lo.hit_wire_cm)); 
  kink.max_angle = max.angle;
  kink.est_angle = est_angle;

  kink.position_max = {max.hit_time_cm, max.hit_wire_cm};
  kink.position_lo = {lo.hit_time_cm, lo.hit_wire_cm};
  kink.position_hi = {hi.hit_time_cm, hi.hit_wire_cm};
  
  kink.vec_lo_at_max = max.lo_vector;
  kink.vec_hi_at_max = max.hi_vector;
  
  kink.vec_lo_at_halfmax_lo = lo.lo_vector;
  kink.vec_hi_at_halfmax_lo = lo.hi_vector;
  
  kink.vec_lo_at_halfmax_hi = hi.lo_vector;
  kink.vec_hi_at_halfmax_hi = hi.hi_vector;

  kink.fit_chi2 = fit.chi2;
  kink.fit_angle = fit.angle;
  kink.fit_pitch = fit.pitch;

  return kink;
}

float EstimateHeight	(	const std::vector< sbn::PCAngle > &	angles,
		unsigned	lo,
		unsigned	hi
	)

Definition at line 189 of file PCAngleKinkFinder_module.cc.

                                                                                    {
  static const unsigned NAVG = 3;
  static const unsigned IAVG = (NAVG - 1) / 2;

  // if the "peak is too small -- ignore
  if (hi - lo + 1 < NAVG) return -1.;

  std::vector<float> avgd;
  // average out the vector
  for (unsigned i = IAVG + lo; i <= hi - IAVG; i++) {
    float sum = 0.;
    for (unsigned j = i - IAVG; j < i - IAVG + NAVG; j++) {
      sum += angles[j].angle;
    }

    avgd.push_back(sum / NAVG);
  }

  return *std::max_element(avgd.begin(), avgd.end());
}

double ExpectedCosAngle	(	double	dist_to_center,
		double	radius,
		double	mean_distance,
		double	angle
	)

Definition at line 290 of file PCAngleKinkFinder_module.cc.

                                                                                                  {
  if (dist_to_center >= radius) return 1.;

  double theta = angle;
  double A = dist_to_center;
  double B = radius - dist_to_center;
  double R = radius;
  double l = mean_distance;

  // double M01 = (B*(B + lB)*(3*A*(A + lA)*lB + (B*B*lA + B*(4*A - lA)*lB + 2*A*lB*lB)*cos(theta))*sin(theta))/(12.*lB*(B*lA + A*lB));
  //double M00 = (A*(A + lA)*lB*(2*B*lA*(2*A + lA) + A*(A - lA)*lB) + B*(B + lB)*cos(theta)*(6*A*lA*(A + lA)*lB + lA*(B*B*lA + B*(4*A - lA)*lB + 2*A*lB*lB)*cos(theta)))/(12.*lA*lB*(B*lA + A*lB));
  //double M11 = (B*(B + lB)*(B*B*lA + B*(4*A - lA)*lB + 2*A*lB*lB)*sin(theta)*sin(theta))/(12.*lB*(B*lA + A*lB));

  double M01 = (B*(B + l)*(-6*A*(A + l)*(A + B - 2*R)*sin(theta) + (3*B*(A + B)*(B + l) - 6*B*(B + l)*R + 2*(2*B + l)*R*R)*sin(2*theta)))/(24.*l*R*R);
  double M00 = (A*(A + l)*(3*A*(A + B)*(A + l) - 6*A*(A + l)*R + 2*(2*A + l)*R*R) + B*(B + l)*cos(theta)*(-6*A*(A + l)*(A + B - 2*R) + (3*B*(A + B)*(B + l) - 6*B*(B + l)*R + 2*(2*B + l)*R*R)*cos(theta)))/(12.*l*R*R);
  double M11 = (B*(B + l)*(3*B*(A + B)*(B + l) - 6*B*(B + l)*R + 2*(2*B + l)*R*R)*sin(theta)*sin(theta))/(12.*l*R*R);

  double trace = M00 + M11;
  double det = M00 * M11 - M01 * M01;
  double eigen = 0.5 * (trace + sqrt(trace * trace - 4 *det));
  double vecA = M01;
  double vecB = eigen - M00;

  return vecA / sqrt(vecA * vecA + vecB * vecB);
}

int FindFWHMHiIndex	(	const std::vector< sbn::PCAngle > &	angles,
		unsigned	start,
		unsigned	end,
		float	amp,
		float	threshold
	)

Definition at line 238 of file PCAngleKinkFinder_module.cc.

                                                                                                                   {
  // handle bad amplitude
  if (amp < 0.) return -1; 

  // go below half-the amplitude
  int halfmax_ind = -1;

  // forward track to the end of this branch
  //
  // We always identify hits on the last part of any branch, so we don't need 
  // to worry about forward tracking to any other part of the branch
  for (unsigned i = end; i < angles.size(); i++) {
    if (angles[i].angle < amp/2.) {
      halfmax_ind = i;
      break;
    } 
  }
  return halfmax_ind;
} 

int FindFWHMLoIndex	(	const std::vector< sbn::PCAngle > &	angles,
		unsigned	start,
		unsigned	end,
		float	amp,
		float	threshold
	)

Definition at line 258 of file PCAngleKinkFinder_module.cc.

                                                                                                                   {
  // handle bad amplitude
  if (amp < 0.) return -1; 

  // backtrack
  int halfmax_ind = -1;

  for (int i = start; i >= 0; i--) {
    if (angles[i].angle < amp / 2.) {
      halfmax_ind = i;      
      break;
    }
  }

  return halfmax_ind;
}

int FindMaxIndex	(	const std::vector< sbn::PCAngle > &	angles,
		unsigned	lo,
		unsigned	hi
	)

Definition at line 211 of file PCAngleKinkFinder_module.cc.

                                                                                {
  static const unsigned NAVG = 3;
  static const unsigned IAVG = (NAVG - 1) / 2;

  // if the "peak is too small -- ignore
  if (hi - lo + 1 < (int)NAVG) return -1.;

  std::vector<float> avgd;
  // average out the vector
  for (unsigned i = lo + IAVG; i <= hi - IAVG; i++) {
    float sum = 0.;
    for (unsigned j = i - IAVG; j < i - IAVG + NAVG; j++) {
      sum += angles[j].angle;
    }

    avgd.push_back(sum / NAVG);
  }

  // index of the averaged
  unsigned avgd_index = std::distance(avgd.begin(), std::max_element(avgd.begin(), avgd.end()));

  // convert this to index in angles list
  unsigned angles_index = avgd_index + IAVG + lo;

  return angles_index;
} 

void FitLikelihood	(	int &	npar,
		double *	g,
		double &	result,
		double *	par,
		int	flag
	)

Definition at line 316 of file PCAngleKinkFinder_module.cc.

                                                                                {
  double angle = par[0];
  double res = par[1];
  double radius = par[2];
  double mean_distance = par[3];

  double chi2 = 0.;
  for (unsigned i = 0; i < fFitDataX.size(); i++) {
    double expectation = acos(ExpectedCosAngle(fFitDataX[i], radius, mean_distance, angle));
    chi2 += (expectation - fFitDataY[i]) * (expectation - fFitDataY[i]) / (fFitDataY[i] * fFitDataY[i] * res * res);
  }

  result = chi2 / fFitDataX.size();
}

std::tuple<std::vector<sbn::PCAngle>, unsigned> FlattenBranch	(	const sbn::PCAnglePlane &	plane,
		unsigned	i_branch,
		bool	allow_incomplete,
		float	min_dist
	)

Definition at line 121 of file PCAngleKinkFinder_module.cc.

                                                                                                                                                {
  std::vector<sbn::PCAngle> ret;

  std::cout << "Flattening angles for branch: " << i_branch << " plane: " << plane.plane.Plane << " hierarchy: ";
  for (int id: plane.branchHierarchy[i_branch]) std::cout << id << " ";
  std::cout << std::endl;

  // start at the end of the branch and count backwards
  for (unsigned i_branch_hierarchy = 0; i_branch_hierarchy < plane.branchHierarchy[i_branch].size(); i_branch_hierarchy++) {
    // look up the index of the current branch
    unsigned this_branch = std::distance(plane.branchIDs.begin(), std::find(plane.branchIDs.begin(), plane.branchIDs.end(), plane.branchHierarchy[i_branch][i_branch_hierarchy])); 
    if (this_branch == plane.branchIDs.size()) {
      std::cout << "BAD: UNFOUND BRANCH ID\n";
      continue;
    }
    std::cout << "BRANCH CHECK: " << plane.branchHierarchy[i_branch][i_branch_hierarchy] << " " << this_branch << " " << i_branch << std::endl;
    for (int i = plane.angles[this_branch].size()-1; i >= 0; i--) {
      // invalid if incomplete
      bool is_complete = (plane.angles[this_branch][i].complete || allow_incomplete) &&
           plane.angles[this_branch][i].dist_to_lo > min_dist && plane.angles[this_branch][i].dist_to_hi > min_dist &&
           plane.angles[this_branch][i].angle > -99.;
       
      if (!is_complete) continue;

      // Check if this hit is replicated in the following branch
      bool found_hit = false;
      if (i_branch_hierarchy > 0) {
        int i_check_branch_hierarchy = i_branch_hierarchy - 1;
        unsigned chk_branch = plane.branchIDs.size();
        int check_hitID = plane.angles[this_branch][i].hitID;
        while (i_check_branch_hierarchy >= 0 && chk_branch == plane.branchIDs.size()) {
          chk_branch = std::distance(plane.branchIDs.begin(), std::find(plane.branchIDs.begin(), plane.branchIDs.end(), plane.branchHierarchy[i_branch][i_check_branch_hierarchy]));
          i_check_branch_hierarchy --;
        }
        if (chk_branch < plane.branchIDs.size()) {
          for (unsigned j = 0; j < plane.angles[chk_branch].size(); j++) {
            if (plane.angles[chk_branch][j].hitID == check_hitID) {
              found_hit = true;
              break;
            }
          }
        }
      }
      // invalid if replicated
      if (found_hit) continue;

      // Valid!
      ret.push_back(plane.angles[this_branch][i]);
    }
  }

  // reverse the angles
  std::reverse(ret.begin(), ret.end());

  // get the number of valid hits in the main branch
  unsigned n_branch_hits = 0;
  for (unsigned i = 0; i < plane.angles[i_branch].size(); i++) {
    n_branch_hits += (plane.angles[i_branch][i].complete || allow_incomplete) &&
           plane.angles[i_branch][i].dist_to_lo > min_dist && plane.angles[i_branch][i].dist_to_hi > min_dist &&
           plane.angles[i_branch][i].angle > -99.;
  }

  // get the index of the start of this branch
  unsigned branch_start = ret.size() - n_branch_hits;
  return {ret, branch_start};
}

Variable Documentation

unsigned fFitCenter

static

Definition at line 88 of file PCAngleKinkFinder_module.cc.

std::vector<float> fFitDataX

static

Definition at line 86 of file PCAngleKinkFinder_module.cc.

std::vector<float> fFitDataY

static

Definition at line 87 of file PCAngleKinkFinder_module.cc.