#include "SmallClusterFinderAlg.h"

Public Member Functions
	SmallClusterFinderAlg (fhicl::ParameterSet const &pset)

void	SelectLocalHitlist (util::GeometryUtilities const &gser, std::vector< art::Ptr< recob::Hit >> hitlist, std::vector< art::Ptr< recob::Hit >> &hitlistlocal, double wire_start, double time_start, double radlimit) const

void	SelectLocalHitlist (util::GeometryUtilities const &gser, std::vector< art::Ptr< recob::Hit >> hitlist, std::vector< art::Ptr< recob::Hit >> &hitlistlocal, double wire_start, double time_start, double radlimit, std::vector< int > &) const

int	GetPlaneAndTPC (art::Ptr< recob::Hit > a, unsigned int &plane, unsigned int &cryostat, unsigned int &time, unsigned int &wire) const

std::vector< art::Ptr< recob::Hit > >	CreateHighHitlist (util::GeometryUtilities const &gser, std::vector< art::Ptr< recob::Hit >> const &hitlist, std::vector< art::Ptr< recob::Hit >> &hitlistleftover) const

void	FindSmallClusters (util::GeometryUtilities const &gser, detinfo::DetectorClocksData const &dataClocks, detinfo::DetectorPropertiesData const &detProp, std::vector< art::Ptr< recob::Hit >> allHits)

void	sortLeftoverHits ()

std::vector< std::vector< art::Ptr< recob::Hit > > >	GetSmallClustersByPlane (unsigned int iPlane)

std::vector< art::Ptr< recob::Hit > >	GetLeftoversByPlane (unsigned int iPlane)

Private Member Functions
void	ClearandResizeVectors ()

Private Attributes
double	fWiretoCm

double	fTimetoCm

double	fWireTimetoCmCm

float	fTimeTick

float	fDriftVelocity

double	fWirePitch

std::vector< unsigned int >	fNWires

unsigned int	fNPlanes

double	fNTimes

geo::WireReadoutGeom const *	wireReadoutGeom = &art::ServiceHandle<geo::WireReadout>()->Get()

std::vector< std::vector< std::vector< art::Ptr< recob::Hit > > > >	smallClustList

std::vector< std::vector< art::Ptr< recob::Hit > > >	hitlistrefined

std::vector< std::vector< art::Ptr< recob::Hit > > >	hitlistbyplane

std::vector< std::vector< art::Ptr< recob::Hit > > >	hitlistleftover

bool	verbose

double	fRadiusSizePar

double	fNHitsInClust

Detailed Description

Definition at line 58 of file SmallClusterFinderAlg.h.

Constructor & Destructor Documentation

cluster::SmallClusterFinderAlg::SmallClusterFinderAlg ( fhicl::ParameterSet const & pset )

explicit

Definition at line 56 of file SmallClusterFinderAlg.cxx.

References ClearandResizeVectors(), fNHitsInClust, fNPlanes, fRadiusSizePar, fhicl::ParameterSet::get(), geo::WireReadoutGeom::Nplanes(), verbose, and wireReadoutGeom.

 {
   fNPlanes = wireReadoutGeom->Nplanes();
   fRadiusSizePar = pset.get<double>("RadiusSizePar");
   fNHitsInClust = pset.get<double>("NHitsInClust");
   verbose = pset.get<bool>("Verbose");
   ClearandResizeVectors();
 }

Member Function Documentation

void cluster::SmallClusterFinderAlg::ClearandResizeVectors ( )

private

Definition at line 187 of file SmallClusterFinderAlg.cxx.

References fNPlanes, hitlistbyplane, hitlistleftover, hitlistrefined, and smallClustList.

Referenced by FindSmallClusters(), and SmallClusterFinderAlg().

 {
   smallClustList.clear();
   hitlistbyplane.clear();
   hitlistrefined.clear();
   hitlistleftover.clear();
   smallClustList.resize(fNPlanes);
   hitlistbyplane.resize(fNPlanes);
   hitlistrefined.resize(fNPlanes);
   hitlistleftover.resize(fNPlanes);
 }

std::vector< art::Ptr< recob::Hit > > cluster::SmallClusterFinderAlg::CreateHighHitlist	(	util::GeometryUtilities const &	gser,
		std::vector< art::Ptr< recob::Hit >> const &	hitlist,
		std::vector< art::Ptr< recob::Hit >> &	hitlistleftover
	)		const

Definition at line 262 of file SmallClusterFinderAlg.cxx.

References fNHitsInClust, fRadiusSizePar, GetPlaneAndTPC(), recob::Hit::PeakTime(), SelectLocalHitlist(), and verbose.

Referenced by FindSmallClusters().

 {
 
   std::vector<art::Ptr<recob::Hit>>
     hitlist_total; //This is the final result, a list of hits that are small clusters
 
   std::vector<art::Ptr<recob::Hit>> hitlistlocal;
 
   for (unsigned int ix = 0; ix < hitlist.size(); ix++) {
 
     art::Ptr<recob::Hit> const& theHit = hitlist[ix]; //grab a hit from the list
 
     double time = theHit->PeakTime();
     unsigned int plane(0), cstat(0), tpc(0), wire(0);
     GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
     //use the wire and time of this hit as a seed.
     // ^^^^^ This could probably be optimized?
 
     //get ALL of the hits from hitlist that are within the distance fRadiusSizePar of the seed hit.
     SelectLocalHitlist(gser, hitlist, hitlistlocal, (double)wire, time, fRadiusSizePar);
 
     if (hitlistlocal.size() < fNHitsInClust) {
       hitlist_total.push_back(theHit); //Add this hit if there are less than fNHitsInClust nearby.
       if (verbose)
         std::cout << " adding hit @ w,t " << wire << " " << time << " on plane " << plane
                   << std::endl;
     }
     else {
       //Add this hit to the leftover pile
       leftovers.push_back(theHit);
     }
 
     hitlistlocal.clear(); //clear the local hit list, and look at the next hit.
   }
 
   /*
     This method could definitely be optimized.  It creates a local hit list for each particle,
     while if there is a local hit list that is sufficiently separated from all others it should be OK to
     add them all at once.  This is a task for future coders!
   */
 
   return hitlist_total;
 }

void cluster::SmallClusterFinderAlg::FindSmallClusters	(	util::GeometryUtilities const &	gser,
		detinfo::DetectorClocksData const &	dataClocks,
		detinfo::DetectorPropertiesData const &	detProp,
		std::vector< art::Ptr< recob::Hit >>	allHits
	)

These lines determine the conversion factors to take wires and times to CMs

Definition at line 67 of file SmallClusterFinderAlg.cxx.

References ClearandResizeVectors(), CreateHighHitlist(), detinfo::DetectorPropertiesData::DriftVelocity(), detinfo::DetectorPropertiesData::Efield(), fDriftVelocity, fNPlanes, fRadiusSizePar, fTimeTick, fTimetoCm, fWirePitch, fWireTimetoCmCm, fWiretoCm, GetPlaneAndTPC(), hitlistbyplane, hitlistleftover, hitlistrefined, recob::Hit::PeakTime(), geo::WireReadoutGeom::Plane(), detinfo::sampling_rate(), SelectLocalHitlist(), util::size(), smallClustList, detinfo::DetectorPropertiesData::Temperature(), verbose, w, and wireReadoutGeom.

Referenced by cluster::SmallClusterFinder::produce().

 {
   fDriftVelocity = detProp.DriftVelocity(detProp.Efield(), detProp.Temperature());
   fWirePitch = wireReadoutGeom->Plane({0, 0, 0}).WirePitch();
   fTimeTick = sampling_rate(clockData) / 1000.;
   fWiretoCm = fWirePitch;
   fTimetoCm = fTimeTick * fDriftVelocity;
   fWireTimetoCmCm = (fTimeTick * fDriftVelocity) / fWirePitch;
 
   ClearandResizeVectors();
 
   // Catch the case were there are no hits in the event:
   if (allHits.size() == 0) {
     if (verbose) std::cout << " no hits received! exiting " << std::endl;
     return;
   }
 
   //sort the hits into hits by plane
   for (auto const& theHit : allHits) {
     unsigned int p(0), w(0), t(0), cs(0); //c=channel, p=plane, w=wire, but what is t?
     GetPlaneAndTPC(theHit, p, cs, t, w);  //Find out what plane this hit is on.
     //add this hit to the list specific to this plane
     hitlistbyplane[p].push_back(theHit);
   } // End loop on hits.
 
   // Want to check that the wires are OK for each event, and by each wire.
   // This could be done in HitFinder, but I'm doing it here because I want to do
   // some things specifically for Small Clusters... Flags to look for:  Does
   // this
   // wire have more hits that the other wires that also have hits? Significantly
   // more?       Does the RMS of this hit lie well below the peak of the hit?
   // Significantly below?        More?   If a bad wire is found, remove all of
   // it's from the hitlists so they aren't clustered.
 
   //make the refined hit list for each plane.
   for (unsigned int ip = 0; ip < fNPlanes; ip++) {
     hitlistrefined[ip] = CreateHighHitlist(gser, hitlistbyplane[ip], hitlistleftover[ip]);
 
     //Check that the lists are populated correctly:
     if (verbose)
       std::cout << "At plane " << ip << ", found " << hitlistrefined[ip].size() << " hits with "
                 << hitlistleftover[ip].size() << " leftover" << std::endl;
     //add the leftover hits to the correct object:
   }
 
   //Now we have all of the gammas.  This is good!
   //Want to split all of the gammas up into individual clusters.
 
   //I was going to use lists instead of vectors to do this, because lists are more efficient
   //at insertion and removal at arbitrary places.  But, the number of gammas is so small that
   //it just doesn't seem worth it
 
   // The method to split the gammas into individual clusters that I want is:
   // Use the existing method below to find all the hits within a certain distance from a given hit
   // Take ALL of those hits and write them to a cluster.
   // Then, remove those hits from the list of gammas.
 
   // Repeat, until the list of gammas is empty.
 
   // There is an issue of knowing which hits to remove from the total gamma hit list.
   // To solve this, I'm going to overload SelectLocalHitList to take a reference to a vector
   // of ints as an argument.  When a hit is added to the local hit list,
   // it's index will be added to a vector that is returned by reference.
   // It's important to remove these hits in reverse order.  This is because removing a hit
   // changes the index of all of the hits after it
 
   //Now we need to take the lists of small clusters and sort it into the individual spots
   //going to end up with smallClustList[plane][iClust][Hit]
   //loop over planes of hits:
 
   for (unsigned int iplane = 0; iplane < fNPlanes; iplane++) {
 
     if (hitlistrefined[iplane].size() == 0) continue;
 
     //write the rest of the gammas one by one, in clusters:
     int i = 1;
 
     std::vector<art::Ptr<recob::Hit>> splittingVector = hitlistrefined[iplane];
 
     while (splittingVector.size() != 0) {
 
       //find the first small cluster of gammas:
       std::vector<int> index;
       std::vector<art::Ptr<recob::Hit>> thiscluster;
       thiscluster.clear();
       index.clear();
 
       //Just use the first hit in the list of gammas:
 
       art::Ptr<recob::Hit> theHit = splittingVector.front(); //grab a hit from the list
 
       double time = theHit->PeakTime();
       unsigned int plane(0), cstat(0), tpc(0), wire(0);
       GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
 
       SelectLocalHitlist(gser, splittingVector, thiscluster, wire, time, fRadiusSizePar, index);
 
       if (verbose)
         std::cout << "Done writing " << thiscluster.size() << " hits to cluster with ID "
                   << plane * 100 + i << std::endl;
       //make sure to add these hits to the object that stores them:
       smallClustList[plane].push_back(thiscluster);
 
       //Lastly, remove the gammas just clustered from the refinded hit list
       while (index.size() != 0) {
         splittingVector.erase(splittingVector.begin() + (index.back()));
         index.pop_back();
       }
       i++;
     }
   }
 }

std::vector< art::Ptr< recob::Hit > > cluster::SmallClusterFinderAlg::GetLeftoversByPlane ( unsigned int iPlane )

Definition at line 331 of file SmallClusterFinderAlg.cxx.

References fNPlanes, and hitlistleftover.

Referenced by cluster::SmallClusterFinder::produce().

 {
   if (iPlane < fNPlanes) return hitlistleftover[iPlane];
   return {};
 }

int cluster::SmallClusterFinderAlg::GetPlaneAndTPC	(	art::Ptr< recob::Hit >	a,
		unsigned int &	plane,
		unsigned int &	cryostat,
		unsigned int &	time,
		unsigned int &	wire
	)		const

Definition at line 310 of file SmallClusterFinderAlg.cxx.

References w.

Referenced by CreateHighHitlist(), FindSmallClusters(), and SelectLocalHitlist().

 {
   p = a->WireID().Plane;
   t = a->PeakTime();
   w = a->WireID().Wire;
 
   return 0;
 }

std::vector< std::vector< art::Ptr< recob::Hit > > > cluster::SmallClusterFinderAlg::GetSmallClustersByPlane ( unsigned int iPlane )

Definition at line 325 of file SmallClusterFinderAlg.cxx.

References fNPlanes, and smallClustList.

Referenced by cluster::SmallClusterFinder::produce().

 {
   if (iPlane < fNPlanes) return smallClustList[iPlane];
   return {};
 }

void cluster::SmallClusterFinderAlg::SelectLocalHitlist	(	util::GeometryUtilities const &	gser,
		std::vector< art::Ptr< recob::Hit >>	hitlist,
		std::vector< art::Ptr< recob::Hit >> &	hitlistlocal,
		double	wire_start,
		double	time_start,
		double	radlimit
	)		const

Definition at line 205 of file SmallClusterFinderAlg.cxx.

References util::GeometryUtilities::Get2DDistance(), and GetPlaneAndTPC().

Referenced by CreateHighHitlist(), and FindSmallClusters().

 {
   //loop over the hits in "hitlist", which should contain the hits we're selecting from
   for (auto const& theHit : hitlist) {
     double time = theHit->PeakTime();
     unsigned int plane, cstat, tpc, wire;
     GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
     //we now know which wire and what time this hit occurred on
 
     //calculate linear distance from start point and orthogonal distance from axis
     double linear_dist = gser.Get2DDistance(wire, time, wire_start, time_start);
 
     if (linear_dist < radlimit) hitlistlocal.push_back(theHit);
   }
 }

void cluster::SmallClusterFinderAlg::SelectLocalHitlist	(	util::GeometryUtilities const &	gser,
		std::vector< art::Ptr< recob::Hit >>	hitlist,
		std::vector< art::Ptr< recob::Hit >> &	hitlistlocal,
		double	wire_start,
		double	time_start,
		double	radlimit,
		std::vector< int > &	index
	)		const

Definition at line 230 of file SmallClusterFinderAlg.cxx.

References util::GeometryUtilities::Get2DDistance(), and GetPlaneAndTPC().

 {
   // loop over the hits in "hitlist", which should contain the hits we're selecting from
   int i = 0; //i keeps track of the index of the hit.
   for (auto const& theHit : hitlist) {
     double time = theHit->PeakTime();
     unsigned int plane, cstat, tpc, wire;
     GetPlaneAndTPC(theHit, plane, cstat, tpc, wire);
     //we now know which wire and what time this hit occurred on
 
     //calculate linear distance from start point and orthogonal distance from axis
     double linear_dist = gser.Get2DDistance(wire, time, wire_start, time_start);
 
     if (linear_dist < radlimit) {
       hitlistlocal.push_back(theHit);
       index.push_back(i);
     }
     i++;
   }
 
   //need to make sure the index array is in order.  Sort it!
   std::sort(index.begin(), index.end());
 }