#include "SmallClusterFinderAlg.h"

Public Member Functions
	SmallClusterFinderAlg ()

	SmallClusterFinderAlg (fhicl::ParameterSet const &pset)

void	reconfigure (fhicl::ParameterSet const &pset)

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

void	SelectLocalHitlist (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 > &)

int	GetPlaneAndTPC (art::Ptr< recob::Hit > a, unsigned int &p, unsigned int &cs, unsigned int &t, unsigned int &w)

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

void	FindSmallClusters (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

art::ServiceHandle< geo::Geometry >	geom

util::GeometryUtilities	gser

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 46 of file SmallClusterFinderAlg.h.

Constructor & Destructor Documentation

cluster::SmallClusterFinderAlg::SmallClusterFinderAlg ( )

METHODS global

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

Definition at line 64 of file SmallClusterFinderAlg.cxx.

References ClearandResizeVectors(), fNPlanes, and geom.

                                                    {
         
         fNPlanes = geom -> Nplanes();
         
         ClearandResizeVectors();
 
 }

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

Definition at line 75 of file SmallClusterFinderAlg.cxx.

References ClearandResizeVectors(), fNPlanes, geom, and reconfigure().

                                                                                 {
         
         fNPlanes = geom -> Nplanes();
         //call reconfigure to grab all the good 
         this -> reconfigure(pset);
         ClearandResizeVectors();
 
 }

Member Function Documentation

void cluster::SmallClusterFinderAlg::ClearandResizeVectors ( )

private

Definition at line 255 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);
         return;
 }

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

Definition at line 339 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;
         //std::cout << "The hitlist size is " << hitlist.size() << std::endl;
         //for (unsigned int iHit =0; iHit < hitlist.size();iHit++) std::cout << (*hitlist[iHit]) << std::endl;
 
         
         for(unsigned int ix = 0; ix<  hitlist.size();  ix++){
                 
                 art::Ptr<recob::Hit> theHit = hitlist[ix]; //grab a hit from the list
         
         double time = theHit->PeakTime() ;  
         unsigned int plane(0),cstat(0),tpc(0),wire(0);
                 //std::cout << "The hit is " << (*theHit) << std::endl;
                 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(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 ( std::vector< art::Ptr< recob::Hit > > allHits )

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

Definition at line 98 of file SmallClusterFinderAlg.cxx.

References ClearandResizeVectors(), CreateHighHitlist(), detinfo::DetectorProperties::DriftVelocity(), detinfo::DetectorProperties::Efield(), fDriftVelocity, fNPlanes, fRadiusSizePar, fTimeTick, fTimetoCm, fWirePitch, fWireTimetoCmCm, fWiretoCm, geom, GetPlaneAndTPC(), hitlistbyplane, hitlistleftover, hitlistrefined, recob::Hit::PeakTime(), detinfo::DetectorProperties::SamplingRate(), SelectLocalHitlist(), smallClustList, detinfo::DetectorProperties::Temperature(), lar::dump::vector(), verbose, and w.

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

 { 
   const detinfo::DetectorProperties* detp = lar::providerFrom<detinfo::DetectorPropertiesService>();
   
         fDriftVelocity=detp->DriftVelocity(detp->Efield(),detp->Temperature());
         fWirePitch = geom -> WirePitch();
         fTimeTick=detp->SamplingRate()/1000.;
         fWiretoCm=fWirePitch;
         fTimetoCm=fTimeTick*fDriftVelocity;
         fWireTimetoCmCm=(fTimeTick*fDriftVelocity)/fWirePitch;
 
         ClearandResizeVectors();
 
         //A vector to hold hits, not yet filled:
         //std::vector< art::Ptr<recob::Hit> > hitlist;
 
         //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;
         }
 
         //std::cout << "Received " << allHits.size() << " hits from the art::producer." << std::endl;
   
         art::Ptr< recob::Hit> theHit;
   
         //sort the hits into hits by plane
         for(std::vector<art::Ptr<recob::Hit> >::iterator HitIter = allHits.begin();
                 HitIter != allHits.end(); HitIter ++){
                 theHit = *HitIter;
         unsigned int p(0),w(0), t(0),cs(0); //c=channel, p=plane, w=wire, but what is t?
         GetPlaneAndTPC(*HitIter,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);
         //std::cout << "Added a hit to plane " << p << std::endl;
     } // 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.
         
         //Start by looping over planes:
         for (unsigned int ip=0; ip< fNPlanes; ip++){
                 //Need to sort the hits into which are on the wires...
                 //Of course, larsoft doesn't allow the wires to know which hits are on them.
                 //So I guess we'll circumvent that by just making up our own struct for this.
         }
         
         
         //make the refined hit list for each plane.
         for(unsigned int ip=0;ip<fNPlanes;ip++) {
                 //std::cout << "Before refining, " << hitlistbyplane[ip].size() << " hits on plane " << ip << std::endl;
         hitlistrefined[ip]=CreateHighHitlist( 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
         //
         //Also, larsoft is so slow on its own that i think the difference is undetectable
         
         // 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
     
 
   
  
         // make an art::PtrVector of the clusters
         //std::auto_ptr<std::vector<recob::Cluster> > SmallClusterFinder(new std::vector<recob::Cluster>);
         //std::auto_ptr< art::Assns<recob::Cluster, recob::Hit> > assn(new art::Assns<recob::Cluster, recob::Hit>);
         
         //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){
                     
                    // std::cout << "\nThe hits remaining to be spilt are:" << std::endl;
                     //for (unsigned int j = 0; j < splittingVector.size();j++){
                     //  std::cout << *splittingVector[j] << std::endl;
                     //}
                     
                     //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);
                         //std::cout << "The hit is " << (*theHit) << std::endl;
                         GetPlaneAndTPC(theHit,plane,cstat,tpc,wire);
                     
                    // std::cout << "The time of the seed hit is " << time << std::endl;
                     
                     SelectLocalHitlist(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){
                         //for (int j = 0; j < index.size();j++) std::cout << index[j] << " ";
                         //std::cout << "\nSize of index is " << index.size() << std::endl;
                         //std::cout << "Size of hitlistrefined is " << splittingVector.size() << std::endl;
                         //std::cout << "Erasing element at " << index.back() << std::endl;
                         splittingVector.erase( splittingVector.begin() + (  index.back() ) );
                         index.pop_back();
                 }
                         i++;
                 }
     }
         return;
 }

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

Definition at line 416 of file SmallClusterFinderAlg.cxx.

References fNPlanes, and hitlistleftover.

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

                                                                                                   {
         //Check the plane number:
         if (/* iPlane >= 0 && */ iPlane < fNPlanes)
                 return hitlistleftover[iPlane];
         else {
                 std::vector <art::Ptr< recob::Hit> >  vec;
                 return vec;
         }
 }

int cluster::SmallClusterFinderAlg::GetPlaneAndTPC	(	art::Ptr< recob::Hit >	a,
		unsigned int &	p,
		unsigned int &	cs,
		unsigned int &	t,
		unsigned int &	w
	)

Definition at line 388 of file SmallClusterFinderAlg.cxx.

References recob::Hit::Channel(), geo::GeometryCore::ChannelToWire(), and geom.

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

 {
         art::ServiceHandle<geo::Geometry> geom;
         unsigned int channel = a->Channel(); 
         geom->ChannelToWire(channel);
         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 406 of file SmallClusterFinderAlg.cxx.

References fNPlanes, and smallClustList.

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

                                                                                                                     {
         //Check the plane number:
         if (/* iPlane >= 0 && */ iPlane < fNPlanes)
                 return smallClustList[iPlane];
         else {
                 std::vector< std::vector <art::Ptr< recob::Hit> > > vec;
                 return vec;
         }
 }

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

Constructor

Definition at line 85 of file SmallClusterFinderAlg.cxx.

References fNHitsInClust, fRadiusSizePar, fhicl::ParameterSet::get(), and verbose.

Referenced by cluster::SmallClusterFinder::reconfigure(), and SmallClusterFinderAlg().

 {
         fRadiusSizePar  =pset.get<  double      > ("RadiusSizePar");
         fNHitsInClust   =pset.get<  double      > ("NHitsInClust");
         verbose                 =pset.get<      bool    > ("Verbose");
         
         //fRadiusSizePar is used to exclude hits from a cluster outside of a certain size
         //fNHitsInClust ensures the clusters don't get too big
         //verbose determines the number of printouts - all or nothing, really.
  }

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

Definition at line 274 of file SmallClusterFinderAlg.cxx.

References util::GeometryUtilities::Get2DDistance(), GetPlaneAndTPC(), gser, recob::Hit::PeakTime(), and lar::dump::vector().

Referenced by CreateHighHitlist(), and FindSmallClusters().

 {
         //loop over the hits in "hitlist", which should contain the hits we're selecting from
         for(std::vector < art::Ptr < recob::Hit > >::const_iterator hitIter = hitlist.begin(); hitIter != hitlist.end();  hitIter++){
         art::Ptr<recob::Hit> theHit = (*hitIter);
         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);
     }
     return;
 }

void cluster::SmallClusterFinderAlg::SelectLocalHitlist	(	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
	)

Definition at line 300 of file SmallClusterFinderAlg.cxx.

References util::GeometryUtilities::Get2DDistance(), GetPlaneAndTPC(), gser, recob::Hit::PeakTime(), and lar::dump::vector().

 {
         //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(std::vector < art::Ptr < recob::Hit > >::const_iterator hitIter = hitlist.begin(); hitIter != hitlist.end();  hitIter++){
         art::Ptr<recob::Hit> theHit = (*hitIter);
         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++; 
     }
     //std::cout << "in select local hit list, index is:" << std::endl;
     //for (int j = 0; j < index.size();j++) std::cout << index[j] << " ";
 
     //need to make sure the index array is in order.  Sort it!
     std::sort (index.begin(), index.end());
     
     return;
 }