16 #include "NovaDAQConventions/DAQConventions.h" 39 std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
42 std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
52 std::vector<std::vector<art::Ptr<rb::CellHit> > >
56 std::vector<std::vector<art::Ptr<rb::CellHit> > >
93 produces< std::vector<rb::Cluster> >();
108 std::vector<art::Ptr<rb::CellHit > > hitlist;
109 for(
unsigned int i = 0;
i < hitcol->size(); ++
i){
111 hitlist.push_back(hit);
114 std::unique_ptr< std::vector<rb::Cluster> > clustercol(
new std::vector<rb::Cluster>);
119 evt.
put(std::move(clustercol));
125 std::unique_ptr< std::vector<rb::Cluster> >& clustercol)
131 int maxNeighbors = 1;
133 double lonerGap = 0.;
137 int minNumWeakNeighbors = 0;
151 minNumWeakNeighbors = 0;
162 minNumWeakNeighbors = 0;
170 minNumWeakNeighbors = 1;
182 const int nHits = hitlist.size();
187 std::vector<int> hit_chunkID;
188 std::vector<int> hit_numNeighbors;
195 for (
int i=0;
i<nHits;
i++) {
196 hit_numNeighbors.push_back(0);
198 for (
int i=0;
i<nHits;
i++) {
200 for (
int j=
i+1;
j<nHits;
j++) {
203 hit_numNeighbors[
i]++;
204 hit_numNeighbors[
j]++;
214 int currChunkID = -1;
219 for (
int i=0;
i<nHits;
i++) {
222 if (prevT+bigGap<
HitTime(hit_i)) {
225 hit_chunkID.push_back(currChunkID);
231 while (j<nHits&&hit_numNeighbors[j]<=maxNeighbors) j++;
239 if (!(prevT+bigGap<
HitTime(hit_j))) {
244 for (;
i<=j&&
i<nHits;
i++) {
249 hit_chunkID.push_back(currChunkID);
251 hit_i = hitlist[--
i];
256 const int nChunks = currChunkID+1;
273 bool savedAnyHits =
true;
276 unsigned int numHitsInChunk[nChunks];
278 while (savedAnyHits) {
280 savedAnyHits =
false;
283 for (
int iC=0;iC<nChunks;iC++) {
284 numHitsInChunk[iC] = 0;
286 for (
int i=0;
i<nHits;
i++) {
287 numHitsInChunk[hit_chunkID[
i]]++;
291 for (
int i=0;
i<nHits;
i++) {
295 if (numHitsInChunk[hit_chunkID[
i]]<
fMinHits) {
307 if (numHitsInChunk[hit_chunkID[
j]]<
fMinHits)
continue;
315 if (dt>lonerGap)
break;
320 hit_chunkID[
i] = hit_chunkID[
j];
346 std::vector<int> numWeakNeighbors;
347 numWeakNeighbors.resize(nHits);
349 for (
int i=0;
i<nHits-1;
i++) {
350 currChunkID = hit_chunkID[
i];
351 for (
int j=
i+1;
j<nHits;
j++) {
353 if(hit_chunkID[
j] != currChunkID)
continue;
356 numWeakNeighbors[
i]++;
357 numWeakNeighbors[
j]++;
365 for (
int i=0;
i<nHits;
i++)
367 if(numWeakNeighbors[
i] < minNumWeakNeighbors)
369 currChunkID = hit_chunkID[
i];
371 numHitsInChunk[currChunkID] -= 1;
373 hit_chunkID[
i] = nChunks;
380 std::vector<std::vector<art::Ptr<rb::CellHit> > >* tSlices
381 =
new std::vector<std::vector<art::Ptr<rb::CellHit> > >(nChunks+1);
382 std::vector<std::vector<art::Ptr<rb::CellHit> > >* slices;
385 for (
int i=0;
i<nHits;
i++) {
386 currChunkID = hit_chunkID[
i];
388 if (currChunkID>=0&&numHitsInChunk[currChunkID]>=
fMinHits) {
390 (*tSlices)[currChunkID].push_back(hit);
394 (*tSlices)[nChunks].push_back(hit);
401 slices =
new std::vector<std::vector<art::Ptr<rb::CellHit> > >;
403 for(
unsigned int ts = 0; ts < tSlices->size() - 1; ++ts)
405 int nHitsTime = (*tSlices)[ts].size();
408 std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
411 for(
unsigned int zs = 0; zs < zSlices.size(); ++zs)
416 slices->push_back(zSlices[zs]);
417 nHitsZ += zSlices[zs].size();
420 assert(nHitsTime == nHitsZ);
423 slices->push_back((*tSlices)[nChunks]);
438 std::vector<rb::Cluster> clusters(slices->size());
444 size_t ns = slices->size() - 1;
446 for(
size_t s = 0;
s <
ns; ++
s)
452 for(
size_t h = 0;
h < (*slices)[
s].size(); ++
h)
454 clusters[
s].Add((*slices)[
s][
h]);
460 for(
size_t h = 0;
h < (*slices)[
s].size(); ++
h)
463 (*slices)[
ns].push_back((*slices)[
s][
h]);
471 for(
size_t h = 0;
h < (*slices)[
ns].size(); ++
h)
473 clusters[
ns].Add((*slices)[ns][
h]);
477 clusters[
ns].SetNoise(
true);
479 if(slices)
delete slices;
481 int nHitsCalHit = hitlist.size();
489 for (
size_t iC=0;iC < clusters.size(); ++iC) {
490 if (iC==ns || clusters[iC].NCell()>=
fMinHits) {
491 nHitsSlicer += clusters[iC].NCell();
492 clustercol->push_back(clusters[iC]);
495 assert(clusters[iC].NCell() == 0 &&
" skipping a cluster with a non-zero size" );
498 assert(nHitsCalHit == nHitsSlicer);
508 const int allowedPlaneGap = weak ? 20 : 5;
509 const int allowedCellGap = weak ? 24 : 6;
513 TMath::Abs(hit_i->
Plane()-hit_j->
Plane()) <= allowedPlaneGap &&
514 TMath::Abs(hit_i->
Cell() -hit_j->
Cell() ) <= allowedCellGap );
521 std::unique_ptr<std::vector<rb::Cluster> >& slices)
528 if(
hits.size() == 0)
return;
534 slices->push_back(noiseClust);
538 std::vector<std::vector<art::Ptr<rb::CellHit> > > tSlices;
542 typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
544 it_t itStart =
hits.begin();
545 it_t itEnd =
hits.begin();
548 it_t itTest =
hits.begin();
555 if(itEnd ==
hits.end())
break;
558 if(itEnd ==
hits.end())
break;
562 noiseClust.
Add(*itStart);
576 std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
577 tSlices.push_back(slice);
584 if(itEnd ==
hits.end())
break;
589 std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
590 tSlices.push_back(slice);
593 for(it_t
it = itStart;
it != itEnd; ++
it) noiseClust.
Add(*
it);
597 for(
unsigned int s = 0;
s < tSlices.size(); ++
s){
598 std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
600 for(
unsigned int n = 0;
n < zSlices.size(); ++
n){
602 for(
unsigned int h = 0;
h < zSlices[
n].size(); ++
h){
603 clust.
Add(zSlices[
n][
h]);
605 slices->push_back(clust);
610 slices->push_back(noiseClust);
614 std::vector<std::vector<art::Ptr<rb::CellHit> > >
Slicer:: 617 std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices;
620 if(
hits.size() == 0)
return zSlices;
625 typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
626 it_t itBegin =
hits.begin();
627 it_t itEnd =
hits.begin();
630 const int endPlane = (*itEnd)->Plane();
633 if(itEnd ==
hits.end())
break;
635 const int nextPlane = (*itEnd)->Plane();
638 std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
639 if(slice.size() >
fMinHits || returnAllSlices) zSlices.push_back(slice);
646 if(itBegin !=
hits.end()){
647 std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
648 if(slice.size() >
fMinHits|| returnAllSlices) zSlices.push_back(slice);
651 std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
654 for(
unsigned int s = 0;
s < zSlices.size(); ++
s){
663 std::vector<std::vector<art::Ptr<rb::CellHit> > >
Slicer:: 674 std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
677 if(
hits.size() == 0)
return bSlices;
680 typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
681 it_t itBegin =
hits.begin();
682 it_t itEnd =
hits.begin();
691 if(itEnd ==
hits.end())
break;
696 if(block != nextBlock){
697 std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
698 bSlices.push_back(slice);
705 if(itBegin !=
hits.end()){
706 std::vector<art::Ptr<rb::CellHit> > slice(itBegin,itEnd);
707 bSlices.push_back(slice);
fvar< T > fabs(const fvar< T > &x)
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
unsigned short Plane() const
A collection of associated CellHits.
const daqchannelmap::DAQChannelMap * Map() const
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
DEFINE_ART_MODULE(TestTMapFile)
void NeighborSlice(std::vector< art::Ptr< rb::CellHit > > hitlist, std::unique_ptr< std::vector< rb::Cluster > > &clustercol)
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
std::vector< std::vector< art::Ptr< rb::CellHit > > > WindowSliceZ(std::vector< art::Ptr< rb::CellHit > > &hits, bool returnAllSlices=false) const
bool fSliceByBlock
If true, slices per block are written.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
bool fSeparateZGaps
If true, the neighbor algorithm will separate zGaps.
ProductID put(std::unique_ptr< PROD > &&product)
unsigned short Cell() const
Far Detector at Ash River, MN.
block
print "ROW IS " print row
bool AreNeighbors(art::Ptr< rb::CellHit > hit_i, art::Ptr< rb::CellHit > hit_j, bool weak=false)
Prototype Near Detector on the surface at FNAL.
unsigned int fMinHits
The minimum number of hits needed to create a slice (this is in any view)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
bool fDoNotSliceAnything
If true, one big slice is written out.
Near Detector in the NuMI cavern.
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
float HitTime(art::Ptr< rb::CellHit > hit)
void produce(art::Event &evt)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
std::string fAlgorithm
Neighbor/Window.
assert(nhit_max >=nhit_nbins)
Class to help Slicer4D manage the collection of hits.
std::string fCellHitInput
Read CellHits from Cal/[fCellHitFolderIn].
std::vector< std::vector< art::Ptr< rb::CellHit > > > BlockSlice(std::vector< art::Ptr< rb::CellHit > > &hits) const
virtual block_t computeBlock(lchan logicalchan) const =0
Which block is this lchan in?
void WindowSlice(std::vector< art::Ptr< rb::CellHit > > hitlist, std::unique_ptr< std::vector< rb::Cluster > > &clustercol)
Slicer(fhicl::ParameterSet const &pset)
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).