MMSlicer_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file MMSlicer.cxx
3 // \brief Based on slicer algorithm, but with a ADC cut
4 // Reproducing the trigger process
5 // \authors ricken
6 ///////////////////////////////////////////////////////////////////////////
13 
14 #include "Geometry/Geometry.h"
15 #include "NovaDAQConventions/DAQConventions.h"
16 #include "RecoBase/CellHit.h"
17 #include "RecoBase/Cluster.h"
18 
19 #include <cmath>
20 #include <string>
21 
22 // Slicer header
23 namespace slicer {
24 
25  class MMSlicer : public art::EDProducer {
26  public:
27  explicit MMSlicer(fhicl::ParameterSet const &pset);
28  virtual ~MMSlicer();
29 
30  void produce(art::Event& evt);
31 
32  protected:
33  // internal methods
34  void WindowSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
35  std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
36  float HitTime(art::Ptr<rb::CellHit> hit) { return hit->TNS(); }
37 
38  // Helper for WindowSlice
39  std::vector<std::vector<art::Ptr<rb::CellHit> > >
40  WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits) const;
41 
42 
43  // configuration settings
44  std::string fCellHitInput; ///< Read CellHits from Cal/[fCellHitFolderIn]
45  unsigned int fMinHits; ///< The minimum number of hits needed to create a slice (this is in any view)
46 
47  // Configuration specific to the Window slicer
48  double fTimeWindow; // ns
49  double fZGapAllowed; // planes
51 
52  };
53 
54 }
55 
56 namespace slicer
57 {
58  //----------------------------------------------------------------------
60  fCellHitInput (pset.get< std::string >("CellHitInput") ),
61  fMinHits (pset.get< int >("MinHits") ),
62  fTimeWindow (pset.get< double >("TimeWindow") ),
63  fZGapAllowed (pset.get< int >("ZGapAllowed") ),
64  fADCmin (pset.get< int >("ADCmin") ),
65  fADCmax (pset.get< int >("ADCmax") )
66  {
67  produces< std::vector<rb::Cluster> >();
68  }
69 
70  //----------------------------------------------------------------------
72  {
73  }
74 
75  //----------------------------------------------------------------------
77  {
78  // Load hit list (CellHits)
80  evt.getByLabel(fCellHitInput, hitcol);
81 
82  std::vector<art::Ptr<rb::CellHit > > hitlist;
83  for(unsigned int i = 0; i < hitcol->size(); ++i){
84  art::Ptr<rb::CellHit> hit(hitcol, i);
85  if (hit->ADC() > fADCmin && hit->ADC() < fADCmax)
86  hitlist.push_back(hit);
87  }
88 
89  std::unique_ptr< std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
90  WindowSlice(hitlist, clustercol);
91  evt.put(std::move(clustercol));
92  }
93 
94 
95  //----------------------------------------------------------------------
96  void MMSlicer::
98  std::unique_ptr<std::vector<rb::Cluster> >& slices)
99  {
100  // A timeslice is at most fTimeWindow long, and must contain at least
101  // fMinHits hits. Any fTimeWindow-sized subwindow drawn within it
102  // will also contain at least fMinHits hits.
103 
104  slices->clear();
105  if(hits.size() < fMinHits ) return;
106 
107  std::vector<std::vector<art::Ptr<rb::CellHit> > > tSlices;
108 
110 
111  typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
112  // The inclusive->exclusive range we're considering for the current slice
113  it_t itStart = hits.begin();
114  it_t itEnd = hits.begin();
115 
116  // We use this point to test if we can expand the slice
117  it_t itTest = hits.begin();
118 
119  while(true){
120  // Grow the slice by moving the end until it's fTimeWindow in size
121  while(HitTime(*itEnd) - HitTime(*itStart) < fTimeWindow){
122  ++itEnd;
123  if(itEnd == hits.end()) break; //out of hits
124  }
125  if(itEnd == hits.end()) break; // Break all the way out of the loop
126 
127  // The window doesn't have enough hits, drop first hit and try again
128  if(itEnd - itStart < fMinHits){
129  ++itStart;
130  // assert(itStart <= itEnd);
131  continue;
132  }
133 
134  // The window is long enough and has enough charge
135  // Catch itTest up until it defines the start of a fTimeWindow-sized
136  // window ending on itEnd
137  while (HitTime(*itEnd) - HitTime(*itTest) > fTimeWindow) ++itTest;
138 
139  // If growing the window more would cause the itTest-itEnd window
140  // not to have enough hits, then write out what we have.
141  if(itEnd+1 - itTest < fMinHits){
142  std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
143  tSlices.push_back(slice);
144  itStart = itEnd;
145  continue;
146  }
147 
148  // We can still legally grow the window. Do that and try again
149  ++itEnd;
150  if(itEnd == hits.end()) break;
151  }
152 
153  // If there are enough hits in the final slice, make it one too
154  if(itEnd - itStart > fMinHits){
155  std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
156  tSlices.push_back(slice);
157  }
158 
159  // std::cout<<"the number of time slices: "<<tSlices.size()<<std::endl;
160 
161  // Split every time slice up based on z if it has gaps in
162  for(unsigned int s = 0; s != tSlices.size(); ++s){
163  std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
164  WindowSliceZ(tSlices[s]);
165 
166  for(unsigned int n = 0; n != zSlices.size(); ++n){
167  rb::Cluster clust;
168  //Start adding hits:
169  for(unsigned int h = 0; h != zSlices[n].size(); ++h){
170  clust.Add(zSlices[n][h]);
171  }
172  if (clust.NXCell()>=3 && clust.NYCell()>=3)
173  slices->push_back(clust);
174  }
175  }
176  }
177 
178  ////////////////////////////////////////////////////////////////////////
179  std::vector<std::vector<art::Ptr<rb::CellHit> > > MMSlicer::
181  {
182  std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices;
183 
185 
186  // std::cout<<"total number hits in the time slice: "<<hits.size()<<std::endl;
187 
188  if (hits.size()<fMinHits+1) return zSlices;
189 
190  // The extremes, inclusive and exclusive of our candidate slice
191  typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
192  it_t itBegin = hits.begin();
193  it_t itEnd = hits.begin();
194 
195  while(itEnd != hits.end()-1){
196  const int endPlane = (*itEnd)->Plane();
197  ++itEnd;
198  const int nextPlane = (*itEnd)->Plane();
199 
200 
201  if(nextPlane - endPlane > fZGapAllowed){
202  std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
203  if (itEnd-itBegin > fMinHits) zSlices.push_back(slice);
204 
205  // if(slice.size() > fMinHits) zSlices.push_back(slice);
206  // Start a new slice
207  itBegin = itEnd;
208  }
209  } // end while
210 
211 
212  // Whatever hits are left, give them their own slice
213  // std::cout<<"remaing hits number: "<<itEnd-itBegin<<std::endl;
214  std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
215 
216  if (slice.size() > fMinHits) zSlices.push_back(slice);
217 
218  return zSlices;
219  }
220 
221  ////////////////////////////////////////////////////////////////////////
222 
223 } // end namespace slicer
224 /////////////////////////////////////////////////////////////////////////
225 
226 
227 
229 
float TNS() const
Definition: CellHit.h:46
std::string fCellHitInput
Read CellHits from Cal/[fCellHitFolderIn].
A collection of associated CellHits.
Definition: Cluster.h:47
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
DEFINE_ART_MODULE(TestTMapFile)
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
unsigned int fMinHits
The minimum number of hits needed to create a slice (this is in any view)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
void hits()
Definition: readHits.C:15
int evt
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void produce(art::Event &evt)
MMSlicer(fhicl::ParameterSet const &pset)
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: event.h:1
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
Class to help Slicer4D manage the collection of hits.
float HitTime(art::Ptr< rb::CellHit > hit)
std::vector< std::vector< art::Ptr< rb::CellHit > > > WindowSliceZ(std::vector< art::Ptr< rb::CellHit > > &hits) const
Encapsulate the geometry of one entire detector (near, far, ndos)
void WindowSlice(std::vector< art::Ptr< rb::CellHit > > hitlist, std::unique_ptr< std::vector< rb::Cluster > > &clustercol)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124