cluster::ClusterParamsAlg Class Reference

#include "ClusterParamsAlg.h"

Public Member Functions
	ClusterParamsAlg ()
	ClusterParamsAlg (const std::vector< util::PxHit > &)
void	Initialize ()
void	SetMinNHits (size_t nhit)
size_t	MinNHits () const
int	SetHits (const std::vector< util::PxHit > &)
void	SetRefineDirectionQMin (double qmin)
void	SetVerbose (bool yes=true)
template<typename Stream >
void	TimeReport (Stream &stream) const
void	GetFANNVector (std::vector< float > &data)
void	PrintFANNVector ()
void	FillParams (util::GeometryUtilities const &gser, bool override_DoGetAverages=false, bool override_DoGetRoughAxis=false, bool override_DoGetProfileInfo=false, bool override_DoRefineStartPointsAndDirection=false, bool override_DoGetFinalSlope=false, bool override_DoTrackShowerSep=false, bool override_DoEndCharge=false)
const cluster_params &	GetParams () const
void	GetAverages (bool override=false)
void	GetRoughAxis (bool override=false)
void	GetProfileInfo (util::GeometryUtilities const &gser, bool override=false)
void	RefineStartPoints (util::GeometryUtilities const &gser)
void	GetFinalSlope (util::GeometryUtilities const &gser, bool override=false)
void	GetEndCharges (util::GeometryUtilities const &gser, bool override_=false)
void	RefineDirection (bool override=false)
void	RefineStartPointAndDirection (util::GeometryUtilities const &gser, bool override=false)
void	TrackShowerSeparation (bool override=false)
void	setNeuralNetPath (std::string s)
void	FillPolygon (util::GeometryUtilities const &gser)
void	GetOpeningAngle ()
const util::PxPoint &	RoughStartPoint ()
const util::PxPoint &	RoughEndPoint ()
double	RoughSlope ()
double	RoughIntercept ()
double	StartCharge (util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
	Returns the expected charge at the beginning of the cluster. More...
double	EndCharge (util::GeometryUtilities const &gser, float length=1., unsigned int nbins=10)
	Returns the expected charge at the end of the cluster. More...
float	MultipleHitWires ()
	Returns the number of multiple hits per wire. More...
float	MultipleHitDensity (util::GeometryUtilities const &gser)
	Returns the number of multiple hits per wire. More...
void	EnableFANN ()
void	DisableFANN ()
size_t	GetNHits () const
const std::vector< util::PxHit > &	GetHitVector () const
int	Plane () const
void	SetPlane (int p)

Public Attributes
cluster::cluster_params	fParams
std::string	fNeuralNetPath
std::vector< std::string >	fTimeRecord_ProcName
std::vector< double >	fTimeRecord_ProcTime

Protected Member Functions
double	IntegrateFitCharge (util::GeometryUtilities const &gser, double from_length, double to_length, unsigned int fit_first_bin, unsigned int fit_end_bin)
	Integrates the charge between two positions in the cluster axis. More...

Static Protected Member Functions
static double	LinearIntegral (double m, double q, double x1, double x2)
	Returns the integral of f(x) = mx + q defined in [x1, x2]. More...

Protected Attributes
size_t	fMinNHits
	Cut value for # hits: below this value clusters are not evaluated. More...
std::vector< util::PxHit >	fHitVector
bool	verbose
std::vector< double >	fChargeCutoffThreshold
int	fPlane
double	fQMinRefDir
std::vector< double >	fChargeProfile
std::vector< double >	fCoarseChargeProfile
std::vector< double >	fChargeProfileNew
int	fCoarseNbins
int	fProfileNbins
int	fProfileMaximumBin
double	fProfileIntegralForward
double	fProfileIntegralBackward
double	fProjectedLength
double	fBeginIntercept
double	fEndIntercept
double	fInterHigh_side
double	fInterLow_side
bool	fFinishedGetAverages
bool	fFinishedGetRoughAxis
bool	fFinishedGetProfileInfo
bool	fFinishedRefineStartPoints
bool	fFinishedRefineDirection
bool	fFinishedGetFinalSlope
bool	fFinishedRefineStartPointAndDirection
bool	fFinishedTrackShowerSep
bool	fFinishedGetEndCharges
double	fRough2DSlope
double	fRough2DIntercept
util::PxPoint	fRoughBeginPoint
util::PxPoint	fRoughEndPoint
bool	enableFANN

Detailed Description

Definition at line 24 of file ClusterParamsAlg.h.

Constructor & Destructor Documentation

cluster::ClusterParamsAlg::ClusterParamsAlg

(

)

Definition at line 32 of file ClusterParamsAlg.cxx.

References enableFANN, fMinNHits, Initialize(), and verbose.

  {
    fMinNHits = 10;
    enableFANN = false;
    verbose = true;
    Initialize();
  }

cluster::ClusterParamsAlg::ClusterParamsAlg

(

const std::vector< util::PxHit > &

inhitlist

)

Definition at line 40 of file ClusterParamsAlg.cxx.

References enableFANN, fMinNHits, SetHits(), and verbose.

  {
    fMinNHits = 10;
    enableFANN = false;
    verbose = true;
    SetHits(inhitlist);
  }

Member Function Documentation

void cluster::ClusterParamsAlg::DisableFANN ( )

inline

Definition at line 227 of file ClusterParamsAlg.h.

{ enableFANN = false; }

void cluster::ClusterParamsAlg::EnableFANN

(

)

Definition at line 172 of file ClusterParamsAlg.cxx.

References enableFANN.

  {
    enableFANN = true;
  }

double cluster::ClusterParamsAlg::EndCharge	(	util::GeometryUtilities const &	gser,
		float	length = 1.,
		unsigned int	nbins = 10
	)

Returns the expected charge at the end of the cluster.

Parameters

nbins	use at least this number of charge bins from charge profile
length	space before the end of cluster where to collect charge, in cm the expected charge at the end of the cluster

See also

StartCharge(), IntegrateFitCharge()

This method returns the charge under the last length cm of the cluster. See StartCharge() for a detailed explanation. For even more details, see IntegrateFitCharge().

Definition at line 1448 of file ClusterParamsAlg.cxx.

References fChargeProfile, fHitVector, fProjectedLength, GetProfileInfo(), and IntegrateFitCharge().

Referenced by cluster::StandardClusterParamsAlg::EndCharge(), and GetEndCharges().

  {
    switch (fHitVector.size()) {
    case 0: return 0.;
    case 1:
      return fHitVector.back().charge;
      // the "default" is the rest of the function
    } // switch

    // need the available number of bins and the axis length
    GetProfileInfo(gser);
    const unsigned int MaxBins = fChargeProfile.size();

    // this is the range of the fit:
    const unsigned int fit_first_bin = MaxBins > nbins ? MaxBins - nbins : 0,
                       fit_last_bin = MaxBins;

    // now determine the integration range, in bin units;
    // get to the end, and go length backward;
    // note that length can be pathologic (0, negative...); not our problem!
    const double from = fProjectedLength - length, to = fProjectedLength;

    return IntegrateFitCharge(gser, from, to, fit_first_bin, fit_last_bin);
  } // ClusterParamsAlg::EndCharge()

void cluster::ClusterParamsAlg::FillParams	(	util::GeometryUtilities const &	gser,
		bool	override_DoGetAverages = false,
		bool	override_DoGetRoughAxis = false,
		bool	override_DoGetProfileInfo = false,
		bool	override_DoRefineStartPointsAndDirection = false,
		bool	override_DoGetFinalSlope = false,
		bool	override_DoTrackShowerSep = false,
		bool	override_DoEndCharge = false
	)

Runs all the functions which calculate cluster params and stashes the results in the private ClusterParams struct.

Parameters

override_DoGetAverages	force re-execution of GetAverages()
override_DoGetRoughAxis	force re-execution of GetRoughAxis()
override_DoGetProfileInfo	force re-execution of GetProfileInfo()
override_DoRefineStartPoints	force re-execution of RefineStartPoints()
override_DoGetFinalSlope	force re-execution of GetFinalSlope()
override_DoEndCharge	force re-execution of GetEndCharges()

Definition at line 177 of file ClusterParamsAlg.cxx.

References GetAverages(), GetEndCharges(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), RefineStartPointAndDirection(), and TrackShowerSeparation().

  {
    GetAverages(override_DoGetAverages);
    GetRoughAxis(override_DoGetRoughAxis);
    GetProfileInfo(gser, override_DoGetProfileInfo);
    RefineStartPointAndDirection(gser, override_DoStartPointsAndDirection);
    GetFinalSlope(gser, override_DoGetFinalSlope);
    GetEndCharges(gser, override_DoEndCharge);
    TrackShowerSeparation(override_DoTrackShowerSep);
  }

void cluster::ClusterParamsAlg::FillPolygon

(

util::GeometryUtilities const &

gser

)

Definition at line 1236 of file ClusterParamsAlg.cxx.

References fHitVector, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::PolyObject, and util::GeometryUtilities::SelectPolygonHitList().

  {
    TStopwatch localWatch;
    localWatch.Start();

    if (fHitVector.size()) {
      std::vector<const util::PxHit*> container_polygon;
      gser.SelectPolygonHitList(fHitVector, container_polygon);
      //now making Polygon Object
      std::pair<float, float> tmpvertex;
      //make Polygon Object as in mac/PolyOverlap.cc
      std::vector<std::pair<float, float>> vertices;
      for (unsigned int i = 0; i < container_polygon.size(); i++) {
        tmpvertex = std::make_pair(container_polygon.at(i)->w, container_polygon.at(i)->t);
        vertices.push_back(tmpvertex);
      }
      fParams.PolyObject = Polygon2D(vertices);
    }

    fTimeRecord_ProcName.push_back("FillPolygon");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

void cluster::ClusterParamsAlg::GetAverages

(

bool

override = false

)

Calculates the following variables: mean_charge mean_x mean_y charge_wgt_x charge_wgt_y eigenvalue_principal eigenvalue_secondary multi_hit_wires N_Wires

Parameters

override

force recalculation of variables

Definition at line 195 of file ClusterParamsAlg.cxx.

References lar::util::StatCollector< T, W >::add(), lar::util::StatCollector< T, W >::Average(), cluster::cluster_params::charge_wgt_x, cluster::cluster_params::charge_wgt_y, cluster::cluster_params::eigenvalue_principal, cluster::cluster_params::eigenvalue_secondary, fFinishedGetAverages, fHitVector, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::mean_ADC, cluster::cluster_params::mean_charge, cluster::cluster_params::mean_x, cluster::cluster_params::mean_y, cluster::cluster_params::multi_hit_wires, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Wires, lar::util::StatCollector< T, W >::RMS(), cluster::cluster_params::rms_ADC, cluster::cluster_params::rms_charge, lar::util::StatCollector< T, W >::Sum(), cluster::cluster_params::sum_ADC, and cluster::cluster_params::sum_charge.

Referenced by FillParams(), GetRoughAxis(), cluster::StandardClusterParamsAlg::Integral(), cluster::StandardClusterParamsAlg::IntegralStdDev(), cluster::StandardClusterParamsAlg::MultipleHitDensity(), MultipleHitDensity(), MultipleHitWires(), cluster::StandardClusterParamsAlg::NHits(), cluster::StandardClusterParamsAlg::SummedADC(), and cluster::StandardClusterParamsAlg::SummedADCStdDev().

  {
    if (!override) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedGetAverages) return;
    }

    TStopwatch localWatch;
    localWatch.Start();

    TPrincipal fPrincipal(2, "D");

    fParams.N_Hits = fHitVector.size();

    std::map<double, int> wireMap;

    lar::util::StatCollector<double> charge, sumADC;

    int uniquewires = 0;
    int multi_hit_wires = 0;
    for (auto& hit : fHitVector) {
      double data[2];
      data[0] = hit.w;
      data[1] = hit.t;
      fPrincipal.AddRow(data);
      fParams.charge_wgt_x += hit.w * hit.charge;
      fParams.charge_wgt_y += hit.t * hit.charge;
      charge.add(hit.charge);
      sumADC.add(hit.sumADC);

      if (wireMap[hit.w] == 0) { uniquewires++; }
      if (wireMap[hit.w] == 1) { multi_hit_wires++; }
      wireMap[hit.w]++;
    }

    fParams.sum_charge = charge.Sum();
    fParams.mean_charge = charge.Average();
    fParams.rms_charge = charge.RMS();

    fParams.sum_ADC = sumADC.Sum();
    fParams.mean_ADC = sumADC.Average();
    fParams.rms_ADC = sumADC.RMS();

    fParams.N_Wires = uniquewires;
    fParams.multi_hit_wires = multi_hit_wires;

    if (fPrincipal.GetMeanValues()->GetNrows() < 2) { throw cluster::CRUException(); }

    fParams.mean_x = (*fPrincipal.GetMeanValues())[0];
    fParams.mean_y = (*fPrincipal.GetMeanValues())[1];
    fParams.mean_charge = fParams.sum_charge / fParams.N_Hits;

    if (fParams.sum_charge != 0.) {
      fParams.charge_wgt_x /= fParams.sum_charge;
      fParams.charge_wgt_y /= fParams.sum_charge;
    }
    else { // "SNAFU"; use the mean
      fParams.charge_wgt_x = fParams.mean_x;
      fParams.charge_wgt_y = fParams.mean_y;
    }

    fPrincipal.MakePrincipals();

    fParams.eigenvalue_principal = (*fPrincipal.GetEigenValues())[0];
    fParams.eigenvalue_secondary = (*fPrincipal.GetEigenValues())[1];

    fFinishedGetAverages = true;

    fTimeRecord_ProcName.push_back("GetAverages");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

void cluster::ClusterParamsAlg::GetEndCharges	(	util::GeometryUtilities const &	gser,
		bool	override_ = false
	)

Calculates the following variables: start_charge end_charge

Parameters

override_

force recompute the variables

See also

StartCharge(), EndCharge()

Definition at line 1344 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::end_charge, EndCharge(), fFinishedGetEndCharges, fParams, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::start_charge, and StartCharge().

Referenced by FillParams().

  {
    if (!override_) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedGetEndCharges) return;
    }

    TStopwatch localWatch;
    localWatch.Start();

    fParams.start_charge = StartCharge(gser);
    fParams.end_charge = EndCharge(gser);

    fFinishedGetEndCharges = true;

    fTimeRecord_ProcName.push_back("GetEndCharges");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  } // ClusterParamsAlg::GetEndCharges()

void cluster::ClusterParamsAlg::GetFANNVector

(

std::vector< float > &

data

)

This function returns a feature vector suitable for a neural net This function uses the data from cluster_params but packages it up in a different way, and so is inappropriate to include in clusterParams.hh. That's why it's here.

Parameters

data	takes a reference to a vector< float>

Definition at line 84 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::closing_angle, cluster::cluster_params::closing_angle_charge_wgt, cluster::cluster_params::eigenvalue_principal, cluster::cluster_params::eigenvalue_secondary, fParams, cluster::cluster_params::hit_density_1D, cluster::cluster_params::length, cluster::cluster_params::mean_charge, cluster::cluster_params::modified_hit_density, cluster::cluster_params::multi_hit_wires, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, cluster::cluster_params::N_Wires, cluster::cluster_params::opening_angle, cluster::cluster_params::opening_angle_charge_wgt, PI, cluster::cluster_params::RMS_charge, cluster::cluster_params::sum_charge, and cluster::cluster_params::width.

Referenced by PrintFANNVector(), and TrackShowerSeparation().

  {
    unsigned int length = 13;
    if (data.size() != length) data.resize(length);
    data[0] = fParams.opening_angle / PI;
    data[1] = fParams.opening_angle_charge_wgt / PI;
    data[2] = fParams.closing_angle / PI;
    data[3] = fParams.closing_angle_charge_wgt / PI;
    data[4] = fParams.eigenvalue_principal;
    data[5] = fParams.eigenvalue_secondary;
    data[6] = fParams.width / fParams.length;
    data[7] = fParams.hit_density_1D / fParams.modified_hit_density;
    data[8] = fParams.multi_hit_wires / fParams.N_Wires;
    data[9] = fParams.N_Hits_HC / fParams.N_Hits;
    data[10] = fParams.modified_hit_density;
    data[11] = fParams.RMS_charge / fParams.mean_charge;
    data[12] = log(fParams.sum_charge / fParams.length);
    return;
  }

void cluster::ClusterParamsAlg::GetFinalSlope	(	util::GeometryUtilities const &	gser,
		bool	override = false
	)

Calculates the following variables: hit_density_1D hit_density_2D angle_2d direction

Parameters

override

[description]

Calculates the following variables: angle_2d modified_hit_density

end testing

Definition at line 870 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::angle_2d, fChargeProfile, fChargeProfileNew, fCoarseChargeProfile, fCoarseNbins, fFinishedGetFinalSlope, fFinishedRefineStartPoints, fHitVector, fParams, fProfileNbins, fProjectedLength, fRough2DIntercept, fRough2DSlope, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2Dangle(), util::GeometryUtilities::Get2DDistance(), util::GeometryUtilities::GetPointOnLine(), cluster::cluster_params::hit_density_1D, cluster::cluster_params::length, cluster::cluster_params::mean_charge, cluster::cluster_params::modified_hit_density, cluster::cluster_params::modmeancharge, PI, RefineStartPoints(), cluster::cluster_params::start_point, and verbose.

Referenced by FillParams().

  {
    if (!override) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedGetFinalSlope) return;
      //Try to run the previous function if not yet done.
      if (!fFinishedRefineStartPoints) RefineStartPoints(gser);
    }
    else {
      //Try to run the previous function if not yet done.
      if (!fFinishedRefineStartPoints) RefineStartPoints(gser);
    }

    TStopwatch localWatch;
    localWatch.Start();
    const int NBINS = 720;
    std::vector<int> fh_omega_single(NBINS, 0); //720,-180., 180.

    double current_maximum = 0;
    double curr_max_bin = -1;

    for (auto hit : fHitVector) {

      double omx = gser.Get2Dangle((util::PxPoint*)&hit,
                                   &fParams.start_point); // in rad and assuming cm/cm space.
      int nbin = (omx + TMath::Pi()) * (NBINS - 1) / (2 * TMath::Pi());
      if (nbin >= NBINS) nbin = NBINS - 1;
      if (nbin < 0) nbin = 0;
      fh_omega_single[nbin] += hit.charge;

      if (fh_omega_single[nbin] > current_maximum) {
        current_maximum = fh_omega_single[nbin];
        curr_max_bin = nbin;
      }
    }
    // FIXME: using two different definitions of PI in the same calculation?
    //        2022-04-18 CHG
    fParams.angle_2d = (curr_max_bin / 720 * (2 * TMath::Pi())) - TMath::Pi();
    fParams.angle_2d *= 180 / PI;
    if (verbose) std::cout << " Final 2D angle: " << fParams.angle_2d << " degrees " << std::endl;

    double mod_angle =
      (fabs(fParams.angle_2d) <= 90) ?
        fabs(fParams.angle_2d) :
        180 - fabs(fParams.angle_2d); //want to transfer angle to radians and from 0 to 90.

    if (
      cos(
        mod_angle * PI /
        180.)) { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).
      if (mod_angle <= 75.)
        fParams.modified_hit_density =
          fParams.hit_density_1D / (2.45 * cos(mod_angle * PI / 180.) + 0.2595);
      else
        fParams.modified_hit_density =
          fParams.hit_density_1D / (1.036 * mod_angle * PI / 180. - 0.2561);

      //calculate modified mean_charge:
      fParams.modmeancharge = fParams.mean_charge / (1264 - 780 * cos(mod_angle * PI / 180.));
    }
    else
      fParams.modified_hit_density = fParams.hit_density_1D;

    // do not use for now.
    bool drawProfileHisto = false;
    if (drawProfileHisto) {

      fProfileNbins = (fParams.length);
      double corr_factor = 1;
      if (
        cos(
          mod_angle * PI /
          180.)) { //values here based on fit of hit_density_1D vs. mod_angle defined as above (in ArgoNeuT).

        if (mod_angle <= 75.)
          corr_factor *= (2.45 * cos(mod_angle * PI / 180.) + 0.2595);
        else
          corr_factor *= (1.036 * mod_angle * PI / 180. - 0.2561);
      }

      fProfileNbins =
        (int)(fProfileNbins / 2 * corr_factor + 0.5); // +0.5 to round to nearest sensible value

      if (verbose) std::cout << " number of final profile bins " << fProfileNbins << std::endl;

      fChargeProfile.clear();
      fCoarseChargeProfile.clear();
      fChargeProfile.resize(fProfileNbins, 0);
      fCoarseChargeProfile.resize(fCoarseNbins, 0);

      fChargeProfileNew.clear();
      fChargeProfileNew.resize(fProfileNbins, 0);

      util::PxPoint BeginOnlinePoint;
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &fParams.start_point, BeginOnlinePoint);

      for (auto hit : fHitVector) {

        util::PxPoint OnlinePoint;
        // get coordinates of point on axis.
        gser.GetPointOnLine(fRough2DSlope, &BeginOnlinePoint, &hit, OnlinePoint);

        double linedist = gser.Get2DDistance(&OnlinePoint, &BeginOnlinePoint);
        int fine_bin = (int)(linedist / fProjectedLength * (fProfileNbins - 1));

        if (fine_bin < fProfileNbins) //only fill if bin number is in range
        {
          fChargeProfileNew.at(fine_bin) += hit.charge;
        }
      }

      TH1F* charge_histo = new TH1F("charge dist", "charge dist", fProfileNbins, 0, fProfileNbins);

      TH1F* charge_diff = new TH1F("charge diff", "charge diff", fProfileNbins, 0, fProfileNbins);

      for (int ix = 0; ix < fProfileNbins - 1; ix++) {
        charge_histo->SetBinContent(ix, fChargeProfileNew[ix]);

        if (ix > 2 && ix < fProfileNbins - 3) {
          double diff =
            (1. / 12. * fChargeProfileNew[ix - 2] - 2. / 3. * fChargeProfileNew[ix - 1] +
             2. / 3. * fChargeProfileNew[ix + 1] - 1. / 12. * fChargeProfileNew[ix + 2]) /
            (double)fProfileNbins;
          charge_diff->SetBinContent(ix, diff);
        }
      }

      TCanvas* chCanv = new TCanvas("chCanv", "chCanv", 600, 600);
      chCanv->cd();
      charge_histo->SetLineColor(3);
      charge_histo->Draw("");
      charge_diff->SetLineColor(2);
      charge_diff->Draw("same");
      chCanv->Update();
    }


    fTimeRecord_ProcName.push_back("GetFinalSlope");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());

    fFinishedGetFinalSlope = true;
    return;
  }

const std::vector<util::PxHit>& cluster::ClusterParamsAlg::GetHitVector ( ) const

inline

Definition at line 230 of file ClusterParamsAlg.h.

{ return fHitVector; }

size_t cluster::ClusterParamsAlg::GetNHits ( ) const

inline

Definition at line 229 of file ClusterParamsAlg.h.

Referenced by cluster::StandardClusterParamsAlg::NInputHits(), and cmtool::CPriorityAlgoBase::Priority().

{ return fHitVector.size(); }

void cluster::ClusterParamsAlg::GetOpeningAngle

(

)

const cluster_params& cluster::ClusterParamsAlg::GetParams ( ) const

inline

Definition at line 82 of file ClusterParamsAlg.h.

Referenced by cluster::StandardClusterParamsAlg::EndOpeningAngle(), cluster::StandardClusterParamsAlg::Integral(), cluster::StandardClusterParamsAlg::IntegralStdDev(), cluster::StandardClusterParamsAlg::MultipleHitDensity(), cluster::StandardClusterParamsAlg::NHits(), cmtool::CFAlgoShowerCompat::PrintClusterInfo(), cluster::StandardClusterParamsAlg::StartAngle(), cluster::StandardClusterParamsAlg::StartOpeningAngle(), cluster::StandardClusterParamsAlg::SummedADC(), cluster::StandardClusterParamsAlg::SummedADCStdDev(), and cluster::StandardClusterParamsAlg::Width().

{ return fParams; }

void cluster::ClusterParamsAlg::GetProfileInfo	(	util::GeometryUtilities const &	gser,
		bool	override = false
	)

Calculates the following variables: opening_angle opening_angle_highcharge closing_angle closing_angle_highcharge offaxis_hits

Parameters

override

[description]

Definition at line 337 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::end_point, fBeginIntercept, fChargeCutoffThreshold, fChargeProfile, fCoarseChargeProfile, fCoarseNbins, fEndIntercept, fFinishedGetProfileInfo, fFinishedGetRoughAxis, fHitVector, fInterHigh_side, fInterLow_side, fParams, fPlane, fProfileIntegralBackward, fProfileIntegralForward, fProfileMaximumBin, fProfileNbins, fProjectedLength, fRough2DIntercept, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2DDistance(), util::GeometryUtilities::GetPointOnLine(), util::GeometryUtilities::GetPointOnLineWSlopes(), GetRoughAxis(), leg, cluster::cluster_params::N_Hits, util::PxPoint::plane, cluster::cluster_params::start_point, cluster::cluster_params::sum_charge, util::PxPoint::t, verbose, util::PxPoint::w, weight, and cluster::cluster_params::width.

Referenced by EndCharge(), FillParams(), IntegrateFitCharge(), RefineStartPointAndDirection(), and cluster::StandardClusterParamsAlg::Width().

  {
    if (!override) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedGetProfileInfo) return;
      //Try to run the previous function if not yet done.
      if (!fFinishedGetRoughAxis) GetRoughAxis(true);
    }
    else {
      if (!fFinishedGetRoughAxis) GetRoughAxis(true);
    }

    TStopwatch localWatch;
    localWatch.Start();

    bool drawOrtHistos = false;

    //these variables need to be initialized to other values?
    if (fRough2DSlope == -999.999 || fRough2DIntercept == -999.999) GetRoughAxis(true);

    //get slope of lines orthogonal to those found crossing the shower.
    double inv_2d_slope = 0;
    if (fabs(fRough2DSlope) > 0.001)
      inv_2d_slope = -1. / fRough2DSlope;
    else
      inv_2d_slope = -999999.;

    double InterHigh = -999999;
    double InterLow = 999999;
    double WireHigh = -999999;
    double WireLow = 999999;
    fInterHigh_side = -999999;
    fInterLow_side = 999999;
    double min_ortdist(999), max_ortdist(-999); // needed to calculate width
    //loop over all hits. Create coarse and fine profiles
    // of the charge weight to help refine the start/end
    // points and have a first guess of direction

    for (auto const& hit : fHitVector) {

      //calculate intercepts along
      double intercept = hit.t - inv_2d_slope * (double)(hit.w);
      double side_intercept = hit.t - fRough2DSlope * (double)(hit.w);

      util::PxPoint OnlinePoint;
      // get coordinates of point on axis.
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &hit, OnlinePoint);

      double ortdist = gser.Get2DDistance(&OnlinePoint, &hit);
      if (ortdist < min_ortdist) min_ortdist = ortdist;
      if (ortdist > max_ortdist) max_ortdist = ortdist;

      if (inv_2d_slope != -999999) //almost all cases
      {
        if (intercept > InterHigh) { InterHigh = intercept; }

        if (intercept < InterLow) { InterLow = intercept; }
      }
      else //slope is practically horizontal. Care only about wires.
      {
        if (hit.w > WireHigh) { WireHigh = hit.w; }
        if (hit.w < WireLow) { WireLow = hit.w; }
      }

      if (side_intercept > fInterHigh_side) { fInterHigh_side = side_intercept; }

      if (side_intercept < fInterLow_side) { fInterLow_side = side_intercept; }
    }
    // end of first HitIter loop, at this point we should
    // have the extreme intercepts

    // Second loop. Fill profiles.

    // get projected points at the beginning and end of the axis.

    util::PxPoint HighOnlinePoint, LowOnlinePoint, BeginOnlinePoint, EndOnlinePoint;

    if (inv_2d_slope != -999999) //almost all cases
    {
      gser.GetPointOnLineWSlopes(fRough2DSlope, fRough2DIntercept, InterHigh, HighOnlinePoint);
      gser.GetPointOnLineWSlopes(fRough2DSlope, fRough2DIntercept, InterLow, LowOnlinePoint);
    }
    else //need better treatment of horizontal events.
    {
      util::PxPoint ptemphigh(fPlane, WireHigh, (fInterHigh_side + fInterLow_side) / 2);
      util::PxPoint ptemplow(fPlane, WireLow, (fInterHigh_side + fInterLow_side) / 2);
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemphigh, HighOnlinePoint);
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemplow, LowOnlinePoint);
    }
    if (verbose) {
      std::cout << " Low online point: " << LowOnlinePoint.w << ", " << LowOnlinePoint.t
                << " High: " << HighOnlinePoint.w << ", " << HighOnlinePoint.t << std::endl;
      //define BeginOnlinePoint as the one with lower wire number (for now), adjust intercepts accordingly
    }
    if (HighOnlinePoint.w >= LowOnlinePoint.w) {
      BeginOnlinePoint = LowOnlinePoint;
      fBeginIntercept = InterLow;
      EndOnlinePoint = HighOnlinePoint;
      fEndIntercept = InterHigh;
    }
    else {
      BeginOnlinePoint = HighOnlinePoint;
      fBeginIntercept = InterHigh;
      EndOnlinePoint = LowOnlinePoint;
      fEndIntercept = InterLow;
    }

    fProjectedLength = gser.Get2DDistance(&BeginOnlinePoint, &EndOnlinePoint);

    fCoarseNbins = 2;

    fProfileNbins = (fProjectedLength > 100) ? 100 : fProjectedLength;
    if (fProfileNbins < 10) fProfileNbins = 10;

    fChargeProfile.clear();
    fCoarseChargeProfile.clear();
    fChargeProfile.resize(fProfileNbins, 0);
    fCoarseChargeProfile.resize(fCoarseNbins, 0);

    // Some fitting variables to make a histogram:

    // TODO this is nonsense for small clusters
    std::vector<double> ort_profile;
    const int NBINS = 100;
    ort_profile.resize(NBINS);

    std::vector<double> ort_dist_vect;
    ort_dist_vect.reserve(fHitVector.size());

    double current_maximum = 0;
    for (auto& hit : fHitVector) {

      util::PxPoint OnlinePoint;
      // get coordinates of point on axis.
      gser.GetPointOnLine(fRough2DSlope, &BeginOnlinePoint, &hit, OnlinePoint);
      double ortdist = gser.Get2DDistance(&OnlinePoint, &hit);

      double linedist = gser.Get2DDistance(&OnlinePoint, &BeginOnlinePoint);
      ort_dist_vect.push_back(ortdist);
      int ortbin;
      if (ortdist == 0)
        ortbin = 0;
      else if (max_ortdist == min_ortdist)
        ortbin = 0;
      else
        ortbin = (int)(ortdist - min_ortdist) / (max_ortdist - min_ortdist) * (NBINS - 1);

      ort_profile.at(ortbin) += hit.charge;

      //calculate the weight along the axis, this formula is based on rough guessology.
      // there is no physics motivation behind the particular numbers, A.S.
      // \todo: switch to something that is motivated and easier to
      //        spell than guessology.  C.A.
      double weight = (ortdist < 1.) ? 10 * (hit.charge) : 5 * (hit.charge) / ortdist;

      int fine_bin =
        fProjectedLength ? (int)(linedist / fProjectedLength * (fProfileNbins - 1)) : 0;
      int coarse_bin =
        fProjectedLength ? (int)(linedist / fProjectedLength * (fCoarseNbins - 1)) : 0;

      if (fine_bin < fProfileNbins) //only fill if bin number is in range
      {
        fChargeProfile.at(fine_bin) += weight;

        //find maximum bin on the fly:
        if (fChargeProfile.at(fine_bin) > current_maximum && fine_bin != 0 &&
            fine_bin != fProfileNbins - 1) {
          current_maximum = fChargeProfile.at(fine_bin);
          fProfileMaximumBin = fine_bin;
        }
      }

      if (coarse_bin < fCoarseNbins) //only fill if bin number is in range
        fCoarseChargeProfile.at(coarse_bin) += weight;

    } // end second loop on hits. Now should have filled profile vectors.

    if (verbose) std::cout << "end second loop " << std::endl;

    double integral = 0;
    int ix = 0;
    for (ix = 0; ix < NBINS; ix++) {
      integral += ort_profile.at(ix);
      if (integral >= 0.95 * fParams.sum_charge) { break; }
    }

    fParams.width =
      2 * (double)ix / (double)(NBINS - 1) *
      (double)(max_ortdist -
               min_ortdist); // multiply by two because ortdist is folding in both sides.

    if (verbose) std::cout << " after width  " << std::endl;

    if (drawOrtHistos) {
      TH1F* ortDistHist = new TH1F(
        "ortDistHist", "Orthogonal Distance to axis;Distance (cm)", 100, min_ortdist, max_ortdist);
      TH1F* ortDistHistCharge =
        new TH1F("ortDistHistCharge",
                 "Orthogonal Distance to axis (Charge Weighted);Distance (cm)",
                 100,
                 min_ortdist,
                 max_ortdist);
      TH1F* ortDistHistAndrzej = new TH1F("ortDistHistAndrzej",
                                          "Orthogonal Distance Andrzej weighting",
                                          100,
                                          min_ortdist,
                                          max_ortdist);

      for (int index = 0; index < fParams.N_Hits; index++) {
        ortDistHist->Fill(ort_dist_vect.at(index));
        ortDistHistCharge->Fill(ort_dist_vect.at(index), fHitVector.at(index).charge);
        double weight = (ort_dist_vect.at(index) < 1.) ?
                          10 * (fHitVector.at(index).charge) :
                          (5 * (fHitVector.at(index).charge) / ort_dist_vect.at(index));
        ortDistHistAndrzej->Fill(ort_dist_vect.at(index), weight);
      }
      ortDistHist->Scale(1.0 / ortDistHist->Integral());
      ortDistHistCharge->Scale(1.0 / ortDistHistCharge->Integral());
      ortDistHistAndrzej->Scale(1.0 / ortDistHistAndrzej->Integral());

      TCanvas* ortCanv = new TCanvas("ortCanv", "ortCanv", 600, 600);
      ortCanv->cd();
      ortDistHistAndrzej->SetLineColor(3);
      ortDistHistAndrzej->Draw("");
      ortDistHistCharge->Draw("same");
      ortDistHist->SetLineColor(2);
      ortDistHist->Draw("same");

      TLegend* leg = new TLegend(.4, .6, .8, .9);
      leg->SetHeader("Charge distribution");
      leg->AddEntry(ortDistHist, "Unweighted Hits");
      leg->AddEntry(ortDistHistCharge, "Charge Weighted Hits");
      leg->AddEntry(ortDistHistAndrzej, "Charge Weighted by Guessology");
      leg->Draw();

      ortCanv->Update();
    }

    fProfileIntegralForward = 0;
    fProfileIntegralBackward = 0;

    //calculate the forward and backward integrals counting int the maximum bin.

    for (int ibin = 0; ibin < fProfileNbins; ibin++) {
      if (ibin <= fProfileMaximumBin) fProfileIntegralForward += fChargeProfile.at(ibin);
      if (ibin >= fProfileMaximumBin) fProfileIntegralBackward += fChargeProfile.at(ibin);
    }

    // now, we have the forward and backward integral so start
    // stepping forward and backward to find the trunk of the
    // shower. This is done making sure that the running
    // integral until less than 1% is left and the bin is above
    // a set theshold many empty bins.

    //forward loop
    double running_integral = fProfileIntegralForward;
    int startbin, endbin;
    for (startbin = fProfileMaximumBin; startbin > 1 && startbin < (int)(fChargeProfile.size());
         startbin--) {
      running_integral -= fChargeProfile.at(startbin);
      if (fChargeProfile.at(startbin) < fChargeCutoffThreshold.at(fPlane) &&
          running_integral / fProfileIntegralForward < 0.01)
        break;
      else if (fChargeProfile.at(startbin) < fChargeCutoffThreshold.at(fPlane) &&
               startbin - 1 > 0 &&
               fChargeProfile.at(startbin - 1) < fChargeCutoffThreshold.at(fPlane) &&
               startbin - 2 > 0 &&
               fChargeProfile.at(startbin - 2) < fChargeCutoffThreshold.at(fPlane))
        break;
    }

    //backward loop
    running_integral = fProfileIntegralBackward;
    for (endbin = fProfileMaximumBin; endbin > 0 && endbin < fProfileNbins - 1; endbin++) {
      running_integral -= fChargeProfile.at(endbin);
      if (fChargeProfile.at(endbin) < fChargeCutoffThreshold.at(fPlane) &&
          running_integral / fProfileIntegralBackward < 0.01)
        break;
      else if (fChargeProfile.at(endbin) < fChargeCutoffThreshold.at(fPlane) &&
               endbin + 1 < fProfileNbins - 1 && endbin + 2 < fProfileNbins - 1 &&
               fChargeProfile.at(endbin + 1) < fChargeCutoffThreshold.at(fPlane) &&
               fChargeProfile.at(endbin + 2) < fChargeCutoffThreshold.at(fPlane))
        break;
    }

    // now have profile start and endpoints. Now translate to wire/time.
    // Will use wire/time that are on the rough axis.
    // fProjectedLength is the length from fInterLow to interhigh a
    // long the rough_2d_axis on bin distance is:
    // util::PxPoint OnlinePoint;

    if (inv_2d_slope != -999999) //almost all cases
    {
      double ort_intercept_begin =
        fBeginIntercept + (fEndIntercept - fBeginIntercept) / fProfileNbins * startbin;
      gser.GetPointOnLineWSlopes(
        fRough2DSlope, fRough2DIntercept, ort_intercept_begin, fRoughBeginPoint);

      double ort_intercept_end =
        fBeginIntercept + (fEndIntercept - fBeginIntercept) / fProfileNbins * endbin;
      fRoughBeginPoint.plane = fPlane;

      gser.GetPointOnLineWSlopes(
        fRough2DSlope, fRough2DIntercept, ort_intercept_end, fRoughEndPoint);
      fRoughEndPoint.plane = fPlane;
    }
    else {
      double wire_begin = WireLow + (WireHigh - WireLow) / fProfileNbins * startbin;

      util::PxPoint ptemphigh(fPlane, wire_begin, (fInterHigh_side + fInterLow_side) / 2);
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemphigh, fRoughBeginPoint);
      fRoughBeginPoint.plane = fPlane;

      double wire_end = WireLow + (WireHigh - WireLow) / fProfileNbins * endbin;

      util::PxPoint ptemplow(fPlane, wire_end, (fInterHigh_side + fInterLow_side) / 2);
      gser.GetPointOnLine(fRough2DSlope, fRough2DIntercept, &ptemplow, fRoughEndPoint);
      fRoughEndPoint.plane = fPlane;
    }

    if (verbose)
      std::cout << "  rough start points " << fRoughBeginPoint.w << ", " << fRoughBeginPoint.t
                << " end: " << fRoughEndPoint.w << " " << fRoughEndPoint.t << std::endl;

    // ok. now have fRoughBeginPoint and fRoughEndPoint. No decision about direction has been made yet.
    fParams.start_point = fRoughBeginPoint;
    fParams.end_point = fRoughEndPoint;

    fFinishedGetProfileInfo = true;

    fTimeRecord_ProcName.push_back("GetProfileInfo");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

void cluster::ClusterParamsAlg::GetRoughAxis

(

bool

override = false

)

Calculates the following variables: verticalness fRough2DSlope fRough2DIntercept

Parameters

override

[description]

Definition at line 268 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::charge_wgt_x, cluster::cluster_params::charge_wgt_y, cluster::cluster_params::cluster_angle_2d, fFinishedGetAverages, fFinishedGetRoughAxis, fHitVector, fParams, fRough2DIntercept, fRough2DSlope, fTimeRecord_ProcName, fTimeRecord_ProcTime, GetAverages(), geo::vect::isfinite(), cluster::cluster_params::mean_charge, cluster::cluster_params::mean_x, cluster::cluster_params::mean_y, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, PI, cluster::cluster_params::RMS_charge, cluster::cluster_params::rms_x, and cluster::cluster_params::rms_y.

Referenced by FillParams(), GetProfileInfo(), and cluster::StandardClusterParamsAlg::StartAngle().

  {
    if (!override) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedGetRoughAxis) return;
      //Try to run the previous function if not yet done.
      if (!fFinishedGetAverages) GetAverages(true);
    }
    else {
      //Try to run the previous function if not yet done.
      if (!fFinishedGetAverages) GetAverages(true);
    }

    TStopwatch localWatch;
    localWatch.Start();

    double rmsx = 0.0;
    double rmsy = 0.0;
    double rmsq = 0.0;
    //using the charge weighted coordinates find the axis from slope
    double ncw = 0;
    double sumtime = 0;     //from sum averages
    double sumwire = 0;     //from sum averages
    double sumwiretime = 0; //sum over (wire*time)
    double sumwirewire = 0; //sum over (wire*wire)
    //next loop over all hits again

    for (auto& hit : fHitVector) {
      // First, abuse this loop to calculate rms in x and y
      rmsx += pow(fParams.mean_x - hit.w, 2) / fParams.N_Hits;
      rmsy += pow(fParams.mean_y - hit.t, 2) / fParams.N_Hits;
      rmsq += pow(fParams.mean_charge - hit.charge, 2) / fParams.N_Hits;
      //if charge is above avg_charge

      if (hit.charge > fParams.mean_charge) {
        ncw += 1;
        sumwire += hit.w;
        sumtime += hit.t;
        sumwiretime += hit.w * hit.t;
        sumwirewire += pow(hit.w, 2);
      } //for high charge
    }   //For hh loop

    fParams.rms_x = sqrt(rmsx);
    fParams.rms_y = sqrt(rmsy);
    fParams.RMS_charge = sqrt(rmsq);

    fParams.N_Hits_HC = ncw;
    //Looking for the slope and intercept of the line above avg_charge hits

    if ((ncw * sumwirewire - sumwire * sumwire) > 0.00001)
      fRough2DSlope = (ncw * sumwiretime - sumwire * sumtime) /
                      (ncw * sumwirewire - sumwire * sumwire); //slope for cw
    else
      fRough2DSlope = 999;
    fRough2DIntercept =
      (std::isfinite(fParams.charge_wgt_x) && std::isfinite(fParams.charge_wgt_y)) ?
        fParams.charge_wgt_y - fRough2DSlope * (fParams.charge_wgt_x) : //intercept for cw
        0.;

    //Getthe 2D_angle
    fParams.cluster_angle_2d = atan(fRough2DSlope) * 180 / PI;

    fFinishedGetRoughAxis = true;

    fTimeRecord_ProcName.push_back("GetRoughAxis");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

void cluster::ClusterParamsAlg::Initialize

(

)

Definition at line 126 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::Clear(), fChargeCutoffThreshold, fFinishedGetAverages, fFinishedGetEndCharges, fFinishedGetFinalSlope, fFinishedGetProfileInfo, fFinishedGetRoughAxis, fFinishedRefineDirection, fFinishedRefineStartPointAndDirection, fFinishedRefineStartPoints, fFinishedTrackShowerSep, fHitVector, fParams, fProfileIntegralBackward, fProfileIntegralForward, fProfileMaximumBin, fRough2DIntercept, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::PxPoint::t, and util::PxPoint::w.

Referenced by cluster::StandardClusterParamsAlg::Clear(), ClusterParamsAlg(), and SetHits().

  {
    fTimeRecord_ProcName.clear();
    fTimeRecord_ProcTime.clear();

    TStopwatch localWatch;
    localWatch.Start();

    //--- Initilize attributes values ---//
    fFinishedGetAverages = false;
    fFinishedGetRoughAxis = false;
    fFinishedGetProfileInfo = false;
    fFinishedRefineStartPoints = false;
    fFinishedRefineDirection = false;
    fFinishedGetFinalSlope = false;
    fFinishedRefineStartPointAndDirection = false;
    fFinishedTrackShowerSep = false;
    fFinishedGetEndCharges = false;

    fRough2DSlope = -999.999;     // slope
    fRough2DIntercept = -999.999; // slope

    fRoughBeginPoint.w = -999.999;
    fRoughEndPoint.w = -999.999;

    fRoughBeginPoint.t = -999.999;
    fRoughEndPoint.t = -999.999;

    fProfileIntegralForward = 999.999;
    fProfileIntegralBackward = 999.999;
    fProfileMaximumBin = -999;

    fChargeCutoffThreshold.clear();
    fChargeCutoffThreshold.reserve(3);
    fChargeCutoffThreshold.push_back(500);
    fChargeCutoffThreshold.push_back(500);
    fChargeCutoffThreshold.push_back(1000);

    fHitVector.clear();

    fParams.Clear();

    fTimeRecord_ProcName.push_back("Initialize");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

double cluster::ClusterParamsAlg::IntegrateFitCharge ( util::GeometryUtilities const &  gser, double  from_length, double  to_length, unsigned int  fit_first_bin, unsigned int  fit_end_bin  )

protected

Integrates the charge between two positions in the cluster axis.

Parameters

from_length	position on the axis to start integration from, in cm
to_length	position on the axis to end integration at, in cm
fit_first_bin	first bin for the charge fit
fit_end_bin	next-to-last bin for the charge fit

Returns

the charged fit integrated in the specified range, in ADC counts

See also

StartCharge(), EndCharge()

This function provides an almost-punctual charge at a position in the axis. Since the effective punctual charge is 0 ADC counts by definition, the charge can be integrated for some length. The procedure is made of two steps:

1. the charge profile is parametrized with a linear fit within the specified region
2. an integration of that fit is performed along the segment specified.

The region at point 1. is from fit_first_bin to fit_end_bin. These are specified in bin units. The binning is the one of the charge profile. It is suggested that a few bins are always kept, say 5 to 10, to reduce statistical fluctuations but maintaining a decent hypothesis of linearity along the range. The linear fit weighs all the bins in the profile the same.

The region at point to is from from_length to to_length, and it is measured in cm along the cluster axis, starting at the start of the cluster.

Definition at line 1371 of file ClusterParamsAlg.cxx.

References lar::util::details::SimplePolyFitterDataBase< T, D >::add(), fChargeProfile, lar::util::details::SimplePolyFitterBase< T, 1U >::FitParameters(), fProjectedLength, GetProfileInfo(), and LinearIntegral().

Referenced by EndCharge(), and StartCharge().

  {
    // first compute the information on the charge profile
    GetProfileInfo(gser);

    // first things first
    if (fit_first_bin > fit_end_bin) std::swap(fit_first_bin, fit_end_bin);

    // how many bins can we use?
    const unsigned int nbins =
      std::min(fit_end_bin - fit_first_bin, (unsigned int)fChargeProfile.size());
    if (nbins == 0) return 0;

    // move the specified range within reason
    if (fit_end_bin > fChargeProfile.size()) {
      fit_end_bin = fChargeProfile.size();
      fit_first_bin = fit_end_bin - nbins;
    }

    // if we have to use only one bin, so be it
    if (nbins < 2) return fChargeProfile[fit_first_bin];

    // fit the charge profile vs. bin number;
    // we assume that the first bin (#0) starts from the very beginning of the
    // cluster, that is at axis coordinate 0.
    lar::util::LinearFit<double> fitter;
    for (unsigned int iBin = fit_first_bin; iBin < fit_end_bin; ++iBin) {
      // should we use a Poisson error instead of no error?
      fitter.add((double)iBin, fChargeProfile[iBin]);
    } // for

    // get the linear fit parameters; [0] intercept [1] slope
    lar::util::LinearFit<double>::FitParameters_t fit;
    try {
      fit = fitter.FitParameters();
    }
    catch (std::range_error const&) { // oops...
      // this is actually unexpected, since the only reason for the fit to fail
      // would be a singular fit matrix, that should be pretty much impossible
      // given that the bin coordinates are well behaved and there are at least
      // two of them
      std::cerr << "IntegrateFitCharge(): linear fit failed!" << std::endl;
      return 0.;
    }

    // now determine the bins corresponding to the length to integrate;
    // note that length can be pathologic (0, negative...); not our problem!
    const double length_to_bin = (double(fChargeProfile.size() - 1) / fProjectedLength);
    const double from_bin = from_length * length_to_bin, to_bin = to_length * length_to_bin;

    // we want the integral between from_bin and to_bin now:
    return LinearIntegral(fit[1] /* m */, fit[0] /* q */, from_bin, to_bin);
  } // ClusterParamsAlg::IntegrateFitCharge()

double cluster::ClusterParamsAlg::LinearIntegral ( double  m, double  q, double  x1, double  x2  )

staticprotected

Returns the integral of f(x) = mx + q defined in [x1, x2].

Definition at line 1365 of file ClusterParamsAlg.cxx.

Referenced by IntegrateFitCharge().

  {
    return m / 2. * (cet::square(x2) - cet::square(x1)) + q * (x2 - x1);
  }

size_t cluster::ClusterParamsAlg::MinNHits ( ) const

inline

Definition at line 33 of file ClusterParamsAlg.h.

{ return fMinNHits; }

float cluster::ClusterParamsAlg::MultipleHitDensity

(

util::GeometryUtilities const &

gser

)

Returns the number of multiple hits per wire.

Returns

the number of multiple hits per wire

This returns the number of wires with mmore than one hit belonging to this cluster, divided by the cluster length in cm.

Definition at line 1488 of file ClusterParamsAlg.cxx.

References fHitVector, fParams, GetAverages(), cluster::cluster_params::length, cluster::cluster_params::multi_hit_wires, and RefineStartPoints().

  {
    if (fHitVector.size() < 2) return 0.0F;

    // compute all the averages
    GetAverages();
    RefineStartPoints(gser); // fParams.length

    // return the relevant information
    return std::isnormal(fParams.length) ? fParams.multi_hit_wires / fParams.length : 0.;
  } // ClusterParamsAlg::MultipleHitDensity()

float cluster::ClusterParamsAlg::MultipleHitWires

(

)

Returns the number of multiple hits per wire.

Returns

the number of multiple hits per wire

This returns the fraction of wires that have more than one hit belonging to this cluster.

Definition at line 1476 of file ClusterParamsAlg.cxx.

References fHitVector, fParams, GetAverages(), cluster::cluster_params::multi_hit_wires, and cluster::cluster_params::N_Wires.

  {
    if (fHitVector.size() < 2) return 0.0F;

    // compute all the averages
    GetAverages();

    // return the relevant information
    return std::isnormal(fParams.N_Wires) ? fParams.multi_hit_wires / fParams.N_Wires : 0.;
  } // ClusterParamsAlg::MultipleHitWires()

int cluster::ClusterParamsAlg::Plane ( ) const

inline

Definition at line 231 of file ClusterParamsAlg.h.

{ return fPlane; }

void cluster::ClusterParamsAlg::PrintFANNVector

(

)

For debugging purposes, prints the result of GetFANNVector in a nicely formatted form.

Returns

[description]

Definition at line 104 of file ClusterParamsAlg.cxx.

References GetFANNVector(), and verbose.

  {
    std::vector<float> data;
    GetFANNVector(data);
    if (verbose) {
      std::cout << "Printing FANN input vector:\n"
                << "   Opening Angle (normalized)  ... : " << data[0] << "\n"
                << "   Opening Angle charge weight  .. : " << data[1] << "\n"
                << "   Closing Angle (normalized)  ... : " << data[2] << "\n"
                << "   Closing Angle charge weight  .. : " << data[3] << "\n"
                << "   Principal Eigenvalue  ......... : " << data[4] << "\n"
                << "   Secondary Eigenvalue  ......... : " << data[5] << "\n"
                << "   Width / Length  ............... : " << data[6] << "\n"
                << "   Hit Density / M.H.D.  ......... : " << data[7] << "\n"
                << "   Percent MultiHit Wires  ....... : " << data[8] << "\n"
                << "   Percent High Charge Hits  ..... : " << data[9] << "\n"
                << "   Modified Hit Density  ......... : " << data[10] << "\n"
                << "   Charge RMS / Mean Charge ...... : " << data[11] << "\n"
                << "   log(Sum Charge / Length) ...... : " << data[12] << "\n";
    }
  }

void cluster::ClusterParamsAlg::RefineDirection

(

bool

override = false

)

This function calculates the opening/closing angles It also makes a decision on whether or not to flip the direction the cluster. Then the start point is refined later.

Parameters

override

[description]

Definition at line 1028 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::closing_angle, cluster::cluster_params::closing_angle_charge_wgt, cluster::cluster_params::end_point, fCoarseChargeProfile, fFinishedRefineDirection, fFinishedRefineStartPoints, fHitVector, fParams, fQMinRefDir, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, cluster::cluster_params::mean_charge, cluster::cluster_params::N_Hits, cluster::cluster_params::N_Hits_HC, cluster::cluster_params::opening_angle, cluster::cluster_params::opening_angle_charge_wgt, PI, cluster::cluster_params::start_point, cluster::cluster_params::sum_charge, util::PxPoint::t, verbose, and util::PxPoint::w.

Referenced by cluster::StandardClusterParamsAlg::EndOpeningAngle(), RefineStartPointAndDirection(), and cluster::StandardClusterParamsAlg::StartOpeningAngle().

  {
    //
    // We don't use "override"? Should we remove? 05/01/14
    //
    if (!override) override = true;

    TStopwatch localWatch;
    localWatch.Start();

    // if(!override) { //Override being set, we skip all this logic.
    //   //OK, no override. Stop if we're already finshed.
    //   if (fFinishedRefineDirection) return;
    //   //Try to run the previous function if not yet done.
    //   if (!fFinishedGetProfileInfo) GetProfileInfo(true);
    // } else {
    //   //Try to run the previous function if not yet done.
    //   if (!fFinishedGetProfileInfo) GetProfileInfo(true);
    // }

    // double wire_2_cm = gser.WireToCm();
    // double time_2_cm = gser.TimeToCm();

    // Removing the conversion since these points are set above in cm-cm space
    //   -- Corey

    util::PxPoint this_startPoint, this_endPoint;

    if (fFinishedRefineStartPoints) {
      this_startPoint = fParams.start_point;
      this_endPoint = fParams.end_point;
    }
    else {
      this_startPoint = fRoughBeginPoint;
      this_endPoint = fRoughEndPoint;
    }
    if (verbose) {
      std::cout << "Angle: Start point: (" << this_startPoint.w << ", " << this_startPoint.t
                << ")\n";
      std::cout << "Angle: End point  : (" << this_endPoint.w << ", " << this_endPoint.t << ")\n";
    }
    double endStartDiff_x = (this_endPoint.w - this_startPoint.w);
    double endStartDiff_y = (this_endPoint.t - this_startPoint.t);
    double hit_counter_forward = 0;
    double hit_counter_backward = 0;

    if (verbose && endStartDiff_y == 0 && endStartDiff_x == 0) {
      std::cerr << "Error:  end point and start point are the same!\n";

      fTimeRecord_ProcName.push_back("RefineDirection");
      fTimeRecord_ProcTime.push_back(localWatch.RealTime());
      return;
    }

    double percentage = 0.90;
    double percentage_HC = 0.90 * fParams.N_Hits_HC / fParams.N_Hits;
    const int NBINS = 200;
    const double wgt = 1.0 / fParams.N_Hits;

    // Containers for the angle histograms
    std::vector<float> opening_angle_bin(NBINS, 0.0);
    std::vector<float> closing_angle_bin(NBINS, 0.0);
    std::vector<float> opening_angle_highcharge_bin(NBINS, 0.0);
    std::vector<float> closing_angle_highcharge_bin(NBINS, 0.0);
    std::vector<float> opening_angle_chargeWgt_bin(NBINS, 0.0);
    std::vector<float> closing_angle_chargeWgt_bin(NBINS, 0.0);
    //hard coding this for now, should use SetRefineDirectionQMin function
    fQMinRefDir = 25;

    int index = -1;
    for (auto& hit : fHitVector) {
      index++;
      //skip this hit if below minimum cutoff param
      if (hit.charge < fQMinRefDir) continue;

      //weight_total = hit.charge;

      // Compute forward mean
      double hitStartDiff_x = (hit.w - this_startPoint.w);
      double hitStartDiff_y = (hit.t - this_startPoint.t);

      if (hitStartDiff_x == 0 && hitStartDiff_y == 0) continue;

      double cosangle_start = (endStartDiff_x * hitStartDiff_x + endStartDiff_y * hitStartDiff_y);
      cosangle_start /= (pow(pow(endStartDiff_x, 2) + pow(endStartDiff_y, 2), 0.5) *
                         pow(pow(hitStartDiff_x, 2) + pow(hitStartDiff_y, 2), 0.5));

      if (cosangle_start > 0) {
        // Only take into account for hits that are in "front"
        //no weighted average, works better as flat average w/ min charge cut
        hit_counter_forward++;
      }

      // Compute backward mean
      double hitEndDiff_x = (hit.w - this_endPoint.w);
      double hitEndDiff_y = (hit.t - this_endPoint.t);
      if (hitEndDiff_x == 0 && hitEndDiff_y == 0) continue;

      double cosangle_end = (endStartDiff_x * hitEndDiff_x + endStartDiff_y * hitEndDiff_y) * (-1.);
      cosangle_end /= (pow(pow(endStartDiff_x, 2) + pow(endStartDiff_y, 2), 0.5) *
                       pow(pow(hitEndDiff_x, 2) + pow(hitEndDiff_y, 2), 0.5));

      if (cosangle_end > 0) {
        //no weighted average, works better as flat average w/ min charge cut
        hit_counter_backward++;
      }

      int N_bins_OPEN = (NBINS - 1) * acos(cosangle_start) / PI;
      int N_bins_CLOSE = (NBINS - 1) * acos(cosangle_end) / PI;
      if (N_bins_OPEN < 0) N_bins_OPEN = 0;
      if (N_bins_CLOSE < 0) N_bins_CLOSE = 0;

      opening_angle_chargeWgt_bin[N_bins_OPEN] += hit.charge / fParams.sum_charge;
      closing_angle_chargeWgt_bin[N_bins_CLOSE] += hit.charge / fParams.sum_charge;
      opening_angle_bin[N_bins_OPEN] += wgt;
      closing_angle_bin[N_bins_CLOSE] += wgt;

      //Also fill bins for high charge hits
      if (hit.charge > fParams.mean_charge) {
        opening_angle_highcharge_bin[N_bins_OPEN] += wgt;
        closing_angle_highcharge_bin[N_bins_CLOSE] += wgt;
      }
    }

    //initialize iterators for the angles
    int iBin(0), jBin(0), kBin(0), lBin(0), mBin(0), nBin(0);
    double percentage_OPEN(0.0), percentage_CLOSE(0.0), percentage_OPEN_HC(0.0),
      percentage_CLOSE_HC(0.0), //The last 2 are for High Charge (HC)
      percentage_OPEN_CHARGEWGT(0.0), percentage_CLOSE_CHARGEWGT(0.0);

    for (iBin = 0; percentage_OPEN <= percentage && iBin < NBINS; iBin++) {
      percentage_OPEN += opening_angle_bin[iBin];
    }

    for (jBin = 0; percentage_CLOSE <= percentage && jBin < NBINS; jBin++) {
      percentage_CLOSE += closing_angle_bin[jBin];
    }

    for (kBin = 0; percentage_OPEN_HC <= percentage_HC && kBin < NBINS; kBin++) {
      percentage_OPEN_HC += opening_angle_highcharge_bin[kBin];
    }

    for (lBin = 0; percentage_CLOSE_HC <= percentage_HC && lBin < NBINS; lBin++) {
      percentage_CLOSE_HC += closing_angle_highcharge_bin[lBin];
    }

    for (mBin = 0; percentage_OPEN_CHARGEWGT <= percentage && mBin < NBINS; mBin++) {
      percentage_OPEN_CHARGEWGT += opening_angle_chargeWgt_bin[mBin];
    }

    for (nBin = 0; percentage_CLOSE_CHARGEWGT <= percentage && nBin < NBINS; nBin++) {
      percentage_CLOSE_CHARGEWGT += closing_angle_chargeWgt_bin[nBin];
    }

    double opening_angle = iBin * PI / NBINS;
    double closing_angle = jBin * PI / NBINS;
    double opening_angle_highcharge = kBin * PI / NBINS;
    double closing_angle_highcharge = lBin * PI / NBINS;
    double opening_angle_charge_wgt = mBin * PI / NBINS;
    double closing_angle_charge_wgt = nBin * PI / NBINS;

    double value_1 = closing_angle / opening_angle - 1;
    double value_2 = (fCoarseChargeProfile[0] / fCoarseChargeProfile[1]);
    if (value_2 < 100.0 && value_2 > 0.01)
      value_2 = log(value_2);
    else
      value_2 = 0.0;
    double value_3 = closing_angle_charge_wgt / opening_angle_charge_wgt - 1;

    // Using a sigmoid function to determine flipping.
    // I am going to weigh all of the values above (1, 2, 3) equally.
    // But, since value 2 in particular, I'm going to cut things off at 5.

    if (value_1 > 3) value_1 = 3.0;
    if (value_1 < -3) value_1 = -3.0;
    if (value_2 > 3) value_2 = 3.0;
    if (value_2 < -3) value_2 = -3.0;
    if (value_3 > 3) value_3 = 3.0;
    if (value_3 < -3) value_3 = -3.0;

    double Exp = exp(value_1 + value_2 + value_3);
    double sigmoid = (Exp - 1) / (Exp + 1);

    bool flip = false;
    if (sigmoid < 0) flip = true;
    if (flip) {
      if (verbose) std::cout << "Flipping!" << std::endl;
      std::swap(opening_angle, closing_angle);
      std::swap(opening_angle_highcharge, closing_angle_highcharge);
      std::swap(opening_angle_charge_wgt, closing_angle_charge_wgt);
      std::swap(fParams.start_point, fParams.end_point);
      std::swap(fRoughBeginPoint, fRoughEndPoint);
    }
    else if (verbose) {
      std::cout << "Not Flipping!\n";
    }

    fParams.opening_angle = opening_angle;
    fParams.opening_angle_charge_wgt = opening_angle_charge_wgt;
    fParams.closing_angle = closing_angle;
    fParams.closing_angle_charge_wgt = closing_angle_charge_wgt;

    fFinishedRefineDirection = true;

    fTimeRecord_ProcName.push_back("RefineDirection");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  } //end RefineDirection

void cluster::ClusterParamsAlg::RefineStartPointAndDirection	(	util::GeometryUtilities const &	gser,
		bool	override = false
	)

Definition at line 1261 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::direction, cluster::cluster_params::end_point, fFinishedGetProfileInfo, fFinishedRefineStartPointAndDirection, fParams, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, GetProfileInfo(), RefineDirection(), RefineStartPoints(), cluster::cluster_params::start_point, util::PxPoint::t, verbose, and util::PxPoint::w.

Referenced by FillParams().

  {
    // This function is meant to pick the direction.
    // It refines both the start and end point, and then asks
    // if it should flip.

    TStopwatch localWatch;
    localWatch.Start();

    if (verbose) std::cout << " here!!! " << std::endl;

    if (!override) { //Override being set, we skip all this logic.
      //OK, no override. Stop if we're already finshed.
      if (fFinishedRefineStartPointAndDirection) {
        fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
        fTimeRecord_ProcTime.push_back(localWatch.RealTime());
        return;
      }
      //Try to run the previous function if not yet done.
      if (!fFinishedGetProfileInfo) GetProfileInfo(gser, true);
    }
    else {
      //Try to run the previous function if not yet done.
      if (!fFinishedGetProfileInfo) GetProfileInfo(gser, true);
    }
    if (verbose) {
      std::cout << "REFINING .... " << std::endl;
      std::cout << "  Rough start and end point: " << std::endl;
      std::cout << "    s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
                << std::endl;
      std::cout << "    e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
                << std::endl;
    }
    RefineStartPoints(gser);
    if (verbose) {
      std::cout << "  Once Refined start and end point: " << std::endl;
      std::cout << "    s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
                << std::endl;
      std::cout << "    e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
                << std::endl;
      std::swap(fParams.start_point, fParams.end_point);
      std::swap(fRoughBeginPoint, fRoughEndPoint);
    }
    RefineStartPoints(gser);
    if (verbose) {
      std::cout << "  Twice Refined start and end point: " << std::endl;
      std::cout << "    s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
                << std::endl;
      std::cout << "    e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
                << std::endl;
      std::swap(fParams.start_point, fParams.end_point);
      std::swap(fRoughBeginPoint, fRoughEndPoint);
    }
    RefineDirection();
    if (verbose) {
      std::cout << "  Final start and end point: " << std::endl;
      std::cout << "    s: (" << fParams.start_point.w << ", " << fParams.start_point.t << ")"
                << std::endl;
      std::cout << "    e: (" << fParams.end_point.w << ", " << fParams.end_point.t << ")"
                << std::endl;
    }
    fParams.direction = (fParams.start_point.w < fParams.end_point.w) ? 1 : -1;

    fTimeRecord_ProcName.push_back("RefineStartPointAndDirection");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

void cluster::ClusterParamsAlg::RefineStartPoints

(

util::GeometryUtilities const &

gser

)

Calculates the following variables: length width

Calculates the following variables: length width hit_density_1D hit_density_2D direction

Parameters

override

[description]

Definition at line 686 of file ClusterParamsAlg.cxx.

References util::PxHit::charge, d, cluster::cluster_params::end_point, fFinishedRefineStartPoints, fHitVector, fParams, fRough2DSlope, fRoughBeginPoint, fRoughEndPoint, fTimeRecord_ProcName, fTimeRecord_ProcTime, util::GeometryUtilities::Get2DDistance(), cluster::cluster_params::hit_density_1D, cluster::cluster_params::hit_density_2D, cluster::cluster_params::length, cluster::cluster_params::N_Hits, util::PxPoint::plane, util::GeometryUtilities::SelectLocalHitlist(), cluster::cluster_params::start_point, util::PxPoint::t, verbose, util::PxPoint::w, cluster::cluster_params::width, and z.

Referenced by GetFinalSlope(), MultipleHitDensity(), and RefineStartPointAndDirection().

  {
    TStopwatch localWatch;
    localWatch.Start();

    // need to define physical direction with openind angles and pass that to Ryan's line finder.

    // Direction decision has been moved entirely to Refine Direction()
    // opening and closing angles are already set
    // they are associated with the start and endpoints.
    // If refine direction opted to switch the shower direction,
    // it also switched opening and closing angles.

    // fParams.direction=  ...

    // Direction is decided by RefineDirection()
    util::PxPoint startHit, endHit;
    startHit = fRoughBeginPoint;
    endHit = fRoughEndPoint;

    //Ryan's Shower Strip finder work here.
    //First we need to define the strip width that we want
    double d = 0.6; //this is the width of the strip.... this needs to be tuned to something.
    //===============================================================================================================
    // Will need to feed in the set of hits that we want.
    //  const std::vector<util::PxHit*> whole;
    std::vector<const util::PxHit*> subhit;
    double linearlimit = 8;
    double ortlimit = 12;
    double lineslopetest(fRough2DSlope);
    util::PxHit averageHit;
    //also are we sure this is actually doing what it is supposed to???
    //are we sure this works?
    //is anybody even listening?  Were they eaten by a grue?
    gser.SelectLocalHitlist(
      fHitVector, subhit, startHit, linearlimit, ortlimit, lineslopetest, averageHit);

    if (!(subhit.size()) || subhit.size() < 3) {
      if (verbose)
        std::cout << "Subhit list is empty or too small. Using rough start/end points..."
                  << std::endl;
      fParams.start_point = fRoughBeginPoint;
      fParams.end_point = fRoughEndPoint;

      fFinishedRefineStartPoints = true;

      fTimeRecord_ProcName.push_back("RefineStartPoints");
      fTimeRecord_ProcTime.push_back(localWatch.RealTime());

      return;
    }

    double avgwire = averageHit.w;
    double avgtime = averageHit.t;
    //vertex in tilda-space pair(x-til,y-til)
    std::vector<std::pair<double, double>> vertil;
    //vector of the sum of the distance of a vector to every vertex in tilda-space
    std::vector<double> vs;
    // $$This needs to be corrected//this is the good hits that are between strip
    std::vector<const util::PxHit*> ghits;
    ghits.reserve(subhit.size());
    double fardistcurrent = 0;
    util::PxHit startpoint;
    //KAZU!!! I NEED this too//this will need to come from somewhere...
    //This is some hit that is from the way far side of the entire shower cluster...
    util::PxPoint farhit = fRoughEndPoint;

    //=============building the local vector===================
    //this is clumsy... but just want to get something to work for now.
    //THIS is REAL Horrible and CLUMSY... We can make things into helper functions soon.
    std::vector<util::PxHit> returnhits;
    std::vector<double> radiusofhit;
    std::vector<int> closehits;
    //==============================================================================
    //Now we need to do the transformation into "tilda-space"
    for (unsigned int a = 0; a < subhit.size(); a++) {
      for (unsigned int b = a + 1; b < subhit.size(); b++) {
        if (subhit.at(a)->w != subhit.at(b)->w) {
          double xtil = ((subhit.at(a)->t - avgtime) - (subhit.at(b)->t - avgtime));
          xtil /= ((subhit.at(a)->w - avgwire) - (subhit.at(b)->w - avgwire));
          double ytil = (subhit.at(a)->w - avgwire) * xtil - (subhit.at(a)->t - avgtime);
          //push back the tilda vertex point on the pair
          std::pair<double, double> tv(xtil, ytil);
          vertil.push_back(tv);
        } //if Wires are not the same
      }   //for over b
    }     //for over a
    // look at the distance from a tilda-vertex to all other tilda-verticies
    for (unsigned int z = 0; z < vertil.size(); z++) {
      double vd = 0; //the distance for vertex... just needs to be something 0
      for (unsigned int c = 0; c < vertil.size(); c++)
        vd += std::hypot(vertil.at(z).first - vertil.at(c).first,
                         vertil.at(z).second - vertil.at(c).second);
      vs.push_back(vd);
      vd = 0.0;
    } //for z loop

    if (vs.size() == 0) //al hits on same wire?!
    {
      if (verbose)
        std::cout << "vertil list is empty. all subhits are on the same wire?" << std::endl;
      fParams.start_point = fRoughBeginPoint;
      fParams.end_point = fRoughEndPoint;

      fFinishedRefineStartPoints = true;

      fTimeRecord_ProcName.push_back("RefineStartPoints");
      fTimeRecord_ProcTime.push_back(localWatch.RealTime());
      return;
    }
    //need to find the min of the distance of vertex in tilda-space
    //this will get me the area where things are most linear
    int minvs = std::min_element(vs.begin(), vs.end()) - vs.begin();
    // now use the min position to find the vertex in tilda-space
    //now need to look a which hits are between the tilda lines from the points
    //in the tilda space everything in wire time is based on the new origin which is at the average wire/time
    double tilwire = vertil.at(minvs).first; //slope
    double tiltimet = -vertil.at(minvs).second +
                      d * sqrt(1 + pow(tilwire, 2)); //negative cept is accounted for top strip
    double tiltimeb = -vertil.at(minvs).second -
                      d * sqrt(1 + pow(tilwire, 2)); //negative cept is accounted for bottom strip
    // look over the subhit list and ask for which are inside of the strip
    for (unsigned int a = 0; a < subhit.size(); a++) {
      double dtstrip =
        (-tilwire * (subhit.at(a)->w - avgwire) + (subhit.at(a)->t - avgtime) - tiltimet) /
        sqrt(tilwire * tilwire + 1);
      double dbstrip =
        (-tilwire * (subhit.at(a)->w - avgwire) + (subhit.at(a)->t - avgtime) - tiltimeb) /
        sqrt(tilwire * tilwire + 1);

      if ((dtstrip < 0.0 && dbstrip > 0.0) || (dbstrip < 0.0 && dtstrip > 0.0)) {
        ghits.push_back(subhit.at(a));
      } //if between strip
    }   //for a loop
    //=======================Do stuff with the good hits============================

    //of these small set of hits just fit a simple line.
    //then we will need to see which of these hits is farthest away

    for (unsigned int g = 0; g < ghits.size(); g++) {
      // should call the helper funtion to do the fit
      //but for now since there is no helper function I will just write it here......again
      // now work on calculating the distance in wire time space from the far point
      //farhit needs to be a hit that is given to me
      double fardist = std::hypot(ghits.at(g)->w - farhit.w, ghits.at(g)->t - farhit.t);
      //come back to this... there is a better way to do this probably in the loop
      //there should also be a check that the hit that is farthest away has subsequent hits after it on a few wires
      if (fardist > fardistcurrent) {
        fardistcurrent = fardist;
        //if fardist... this is the point to use for the start point
        startpoint.t = ghits.at(g)->t;
        startpoint.w = ghits.at(g)->w;
        startpoint.plane = ghits.at(g)->plane;
        startpoint.charge = ghits.at(g)->charge;
      }
    } //for ghits loop
    //This can be the new start point

    //should this be here? Id argue this might be done once outside.
    fParams.length = gser.Get2DDistance((util::PxPoint*)&startpoint, &endHit);
    if (fParams.length)
      fParams.hit_density_1D = fParams.N_Hits / fParams.length;
    else
      fParams.hit_density_1D = 0;

    if (fParams.length && fParams.width)
      fParams.hit_density_2D = fParams.N_Hits / fParams.length / fParams.width;
    else
      fParams.hit_density_2D = 0;

    fParams.start_point = startpoint;

    fFinishedRefineStartPoints = true;

    fTimeRecord_ProcName.push_back("RefineStartPoints");
    fTimeRecord_ProcTime.push_back(localWatch.RealTime());
  }

const util::PxPoint& cluster::ClusterParamsAlg::RoughEndPoint ( )

inline

Definition at line 158 of file ClusterParamsAlg.h.

{ return fRoughEndPoint; }

double cluster::ClusterParamsAlg::RoughIntercept ( )

inline

Definition at line 161 of file ClusterParamsAlg.h.

{ return fRough2DIntercept; }

double cluster::ClusterParamsAlg::RoughSlope ( )

inline

Definition at line 160 of file ClusterParamsAlg.h.

{ return fRough2DSlope; }

const util::PxPoint& cluster::ClusterParamsAlg::RoughStartPoint ( )

inline

Definition at line 157 of file ClusterParamsAlg.h.

{ return fRoughBeginPoint; }

int cluster::ClusterParamsAlg::SetHits

(

const std::vector< util::PxHit > &

inhitlist

)

Definition at line 48 of file ClusterParamsAlg.cxx.

References fHitVector, fMinNHits, fPlane, Initialize(), and verbose.

Referenced by ClusterParamsAlg(), and cluster::StandardClusterParamsAlg::SetHitsFromPointers().

  {
    Initialize();

    // Make default values
    // Is done by the struct
    if (!(inhitlist.size())) {
      throw CRUException("Provided empty hit list!");
      return -1;
    }

    fHitVector.reserve(inhitlist.size());
    for (auto h : inhitlist)
      fHitVector.push_back(h);

    fPlane = fHitVector[0].plane;

    if (fHitVector.size() < fMinNHits) {
      if (verbose)
        std::cout << " the hitlist is too small. Continuing to run may result in crash!!! "
                  << std::endl;
      return -1;
    }

    return fHitVector.size();
  }

void cluster::ClusterParamsAlg::SetMinNHits ( size_t  nhit )

inline

Definition at line 31 of file ClusterParamsAlg.h.

Referenced by cmtool::CMManagerBase::SetClusters().

{ fMinNHits = nhit; }

void cluster::ClusterParamsAlg::setNeuralNetPath ( std::string  s )

inline

Definition at line 151 of file ClusterParamsAlg.h.

{ fNeuralNetPath = s; }

void cluster::ClusterParamsAlg::SetPlane

(

int

p

)

Definition at line 75 of file ClusterParamsAlg.cxx.

References cluster::cluster_params::end_point, fHitVector, fParams, fPlane, fRoughBeginPoint, fRoughEndPoint, util::PxPoint::plane, and cluster::cluster_params::start_point.

  {
    fPlane = p;
    for (auto& h : fHitVector)
      h.plane = p;
    fRoughBeginPoint.plane = fRoughEndPoint.plane = p;
    fParams.start_point.plane = fParams.end_point.plane = p;
  }

void cluster::ClusterParamsAlg::SetRefineDirectionQMin ( double  qmin )

inline

Definition at line 37 of file ClusterParamsAlg.h.

{ fQMinRefDir = qmin; }

void cluster::ClusterParamsAlg::SetVerbose ( bool  yes = true )

inline

Definition at line 39 of file ClusterParamsAlg.h.

Referenced by cmtool::CMManagerBase::SetClusters(), and cluster::StandardClusterParamsAlg::SetVerbose().

{ verbose = yes; }

double cluster::ClusterParamsAlg::StartCharge	(	util::GeometryUtilities const &	gser,
		float	length = 1.,
		unsigned int	nbins = 10
	)

Returns the expected charge at the beginning of the cluster.

Parameters

nbins	use at least this number of charge bins from charge profile
length	space at the start of cluster where to collect charge, in cm the expected charge at the beginning of the cluster

See also

EndCharge(), IntegrateFitCharge()

ClusterParamsAlg extracts a binned charge profile, parametrized versus the distance from the start of the cluster. All the charge on the plane orthogonal to cluster axis is collapsed into the point where that plane intersects the axis. The resulting 1D distribution is then binned.

This method returns the charge under the first length cm of the cluster.

This method considers the first nbins of this charge distribution and through a linear fit determines the expected charge at the first bin. Then, it scales the result to reflect how much charge would be deposited in a space of length centimetres, according to this linear fit.

Note that length may be 0 (charge will be 0) or negative (sort of extrapolation ahead of the cluster start).

For more details, see IntegrateFitCharge().

Definition at line 1430 of file ClusterParamsAlg.cxx.

References fHitVector, and IntegrateFitCharge().

Referenced by GetEndCharges(), and cluster::StandardClusterParamsAlg::StartCharge().

  {
    switch (fHitVector.size()) {
    case 0: return 0.;
    case 1:
      return fHitVector.back().charge;
      // the "default" is the rest of the function
    } // switch

    // this is the range of the fit:
    const unsigned int fit_first_bin = 0, fit_last_bin = nbins;

    return IntegrateFitCharge(gser, 0., length, fit_first_bin, fit_last_bin);
  } // ClusterParamsAlg::StartCharge()

template<typename Stream >

void cluster::ClusterParamsAlg::TimeReport

(

Stream &

stream

)

const

Definition at line 345 of file ClusterParamsAlg.h.

  {

    stream << "  <<ClusterParamsAlg::TimeReport>> starts...\n";
    for (size_t i = 0; i < fTimeRecord_ProcName.size(); ++i) {

      stream << "    Function: " << fTimeRecord_ProcName[i].c_str()
             << " ... Time = " << fTimeRecord_ProcTime[i] << " [s]\n";
    }
    stream << "  <<ClusterParamsAlg::TimeReport>> ends...\n";
  }

void cluster::ClusterParamsAlg::TrackShowerSeparation

(

bool

override = false

)

Definition at line 1329 of file ClusterParamsAlg.cxx.

References enableFANN, GetFANNVector(), and verbose.

Referenced by FillParams().

  {
    if (!override) return;
    if (verbose) std::cout << " ---- Trying T/S sep. ------ \n";
    if (enableFANN) {
      if (verbose) std::cout << " ---- Doing T/S sep. ------- \n";
      std::vector<float> FeatureVector, outputVector;
      GetFANNVector(FeatureVector);
    }
    else {
      if (verbose) std::cout << " ---- Failed T/S sep. ------ \n";
    }
  }

Member Data Documentation

bool cluster::ClusterParamsAlg::enableFANN

protected

Definition at line 287 of file ClusterParamsAlg.h.

Referenced by ClusterParamsAlg(), EnableFANN(), and TrackShowerSeparation().

double cluster::ClusterParamsAlg::fBeginIntercept

protected

Definition at line 267 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

std::vector<double> cluster::ClusterParamsAlg::fChargeCutoffThreshold

protected

Definition at line 248 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

std::vector<double> cluster::ClusterParamsAlg::fChargeProfile

protected

Definition at line 254 of file ClusterParamsAlg.h.

Referenced by EndCharge(), GetFinalSlope(), GetProfileInfo(), and IntegrateFitCharge().

std::vector<double> cluster::ClusterParamsAlg::fChargeProfileNew

protected

Definition at line 257 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope().

std::vector<double> cluster::ClusterParamsAlg::fCoarseChargeProfile

protected

Definition at line 255 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), GetProfileInfo(), and RefineDirection().

int cluster::ClusterParamsAlg::fCoarseNbins

protected

Definition at line 259 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), and GetProfileInfo().

double cluster::ClusterParamsAlg::fEndIntercept

protected

Definition at line 268 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

bool cluster::ClusterParamsAlg::fFinishedGetAverages

protected

Definition at line 273 of file ClusterParamsAlg.h.

Referenced by GetAverages(), GetRoughAxis(), and Initialize().

bool cluster::ClusterParamsAlg::fFinishedGetEndCharges

protected

Definition at line 281 of file ClusterParamsAlg.h.

Referenced by GetEndCharges(), and Initialize().

bool cluster::ClusterParamsAlg::fFinishedGetFinalSlope

protected

Definition at line 278 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), and Initialize().

bool cluster::ClusterParamsAlg::fFinishedGetProfileInfo

protected

Definition at line 275 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), Initialize(), and RefineStartPointAndDirection().

bool cluster::ClusterParamsAlg::fFinishedGetRoughAxis

protected

Definition at line 274 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), GetRoughAxis(), and Initialize().

bool cluster::ClusterParamsAlg::fFinishedRefineDirection

protected

Definition at line 277 of file ClusterParamsAlg.h.

Referenced by Initialize(), and RefineDirection().

bool cluster::ClusterParamsAlg::fFinishedRefineStartPointAndDirection

protected

Definition at line 279 of file ClusterParamsAlg.h.

Referenced by Initialize(), and RefineStartPointAndDirection().

bool cluster::ClusterParamsAlg::fFinishedRefineStartPoints

protected

Definition at line 276 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), Initialize(), RefineDirection(), and RefineStartPoints().

bool cluster::ClusterParamsAlg::fFinishedTrackShowerSep

protected

Definition at line 280 of file ClusterParamsAlg.h.

Referenced by Initialize().

std::vector<util::PxHit> cluster::ClusterParamsAlg::fHitVector

protected

This vector holds the pointer to hits. This should be used for computation for speed.

Definition at line 242 of file ClusterParamsAlg.h.

Referenced by EndCharge(), FillPolygon(), GetAverages(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), Initialize(), MultipleHitDensity(), MultipleHitWires(), RefineDirection(), RefineStartPoints(), SetHits(), SetPlane(), and StartCharge().

double cluster::ClusterParamsAlg::fInterHigh_side

protected

Definition at line 269 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

double cluster::ClusterParamsAlg::fInterLow_side

protected

Definition at line 270 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo().

size_t cluster::ClusterParamsAlg::fMinNHits

protected

Cut value for # hits: below this value clusters are not evaluated.

Definition at line 236 of file ClusterParamsAlg.h.

Referenced by ClusterParamsAlg(), and SetHits().

std::string cluster::ClusterParamsAlg::fNeuralNetPath

Definition at line 329 of file ClusterParamsAlg.h.

cluster::cluster_params cluster::ClusterParamsAlg::fParams

Definition at line 327 of file ClusterParamsAlg.h.

Referenced by FillPolygon(), GetAverages(), GetEndCharges(), GetFANNVector(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), Initialize(), MultipleHitDensity(), MultipleHitWires(), RefineDirection(), RefineStartPointAndDirection(), RefineStartPoints(), and SetPlane().

int cluster::ClusterParamsAlg::fPlane

protected

Definition at line 249 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), SetHits(), and SetPlane().

double cluster::ClusterParamsAlg::fProfileIntegralBackward

protected

Definition at line 263 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

double cluster::ClusterParamsAlg::fProfileIntegralForward

protected

Definition at line 262 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

int cluster::ClusterParamsAlg::fProfileMaximumBin

protected

Definition at line 261 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), and Initialize().

int cluster::ClusterParamsAlg::fProfileNbins

protected

Definition at line 260 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), and GetProfileInfo().

double cluster::ClusterParamsAlg::fProjectedLength

protected

Definition at line 264 of file ClusterParamsAlg.h.

Referenced by EndCharge(), GetFinalSlope(), GetProfileInfo(), and IntegrateFitCharge().

double cluster::ClusterParamsAlg::fQMinRefDir

protected

Definition at line 252 of file ClusterParamsAlg.h.

Referenced by RefineDirection().

double cluster::ClusterParamsAlg::fRough2DIntercept

protected

Definition at line 284 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), and Initialize().

double cluster::ClusterParamsAlg::fRough2DSlope

protected

Definition at line 283 of file ClusterParamsAlg.h.

Referenced by GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), Initialize(), and RefineStartPoints().

util::PxPoint cluster::ClusterParamsAlg::fRoughBeginPoint

protected

Definition at line 285 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), Initialize(), RefineDirection(), RefineStartPointAndDirection(), RefineStartPoints(), and SetPlane().

util::PxPoint cluster::ClusterParamsAlg::fRoughEndPoint

protected

Definition at line 286 of file ClusterParamsAlg.h.

Referenced by GetProfileInfo(), Initialize(), RefineDirection(), RefineStartPointAndDirection(), RefineStartPoints(), and SetPlane().

std::vector<std::string> cluster::ClusterParamsAlg::fTimeRecord_ProcName

Definition at line 331 of file ClusterParamsAlg.h.

Referenced by FillPolygon(), GetAverages(), GetEndCharges(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), Initialize(), RefineDirection(), RefineStartPointAndDirection(), and RefineStartPoints().

std::vector<double> cluster::ClusterParamsAlg::fTimeRecord_ProcTime

Definition at line 332 of file ClusterParamsAlg.h.

Referenced by FillPolygon(), GetAverages(), GetEndCharges(), GetFinalSlope(), GetProfileInfo(), GetRoughAxis(), Initialize(), RefineDirection(), RefineStartPointAndDirection(), and RefineStartPoints().

bool cluster::ClusterParamsAlg::verbose

protected

Definition at line 245 of file ClusterParamsAlg.h.

Referenced by ClusterParamsAlg(), GetFinalSlope(), GetProfileInfo(), PrintFANNVector(), RefineDirection(), RefineStartPointAndDirection(), RefineStartPoints(), SetHits(), and TrackShowerSeparation().

---

The documentation for this class was generated from the following files:

• larreco/v09_25_00/source/larreco/RecoAlg/ClusterRecoUtil/ClusterParamsAlg.h
• larreco/v09_25_00/source/larreco/RecoAlg/ClusterRecoUtil/ClusterParamsAlg.cxx