Multiplet_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: Multiplet
3 // Module Type: filter
4 // File: Multiplet_module.cc
5 //
6 // This is a slicer designed to keep only very small slices, with finding SN
7 // nu hits in mind. It doesn't try to do 3D matching, just simply looks for
8 // things that are near in plane number and time. Single hits are probably
9 // best rejected: noise is predominantly singletons, even though so are almost
10 // half of your SN nu hits :(
11 //
12 // Generated at Tue Sep 9 04:40:50 2014 by Alec Habig using artmod
13 // from cetpkgsupport v1_07_00.
14 ////////////////////////////////////////////////////////////////////////
15 
20 
26 
27 #include <memory>
28 #include <list>
29 #include <algorithm>
30 
31 //a small function to get distance for unsigned values
32 template<class T>
33 unsigned distance(const T& t1, const T& t2){
34  return (t2>t1)?(t2-t1):(t1-t2);
35 }
36 
37 template <class T>
38 unsigned delta(const novaddt::DAQHit& h1, const novaddt::DAQHit& h2){
39  return distance(static_cast<T>(h1).val, static_cast<T>(h2).val);
40 }
41 
42 double deltaTDC(const novaddt::TDC& h1, const novaddt::TDC& h2){
43  double DT = (h2>h1)?(double(h2.val-h1.val)):(-double(h1.val-h2.val));
44  return DT+h2.fraction()-h1.fraction();
45 }
46 
47 namespace novaddt {
48  class Multiplet;
49 }
50 
52 public:
53  explicit Multiplet(fhicl::ParameterSet const & p);
54  // The destructor generated by the compiler is fine for classes
55  // without bare pointers or other resource use.
56 
57  // Plugins should not be copied or assigned.
58  Multiplet(Multiplet const &) = delete;
59  Multiplet(Multiplet &&) = delete;
60  Multiplet & operator = (Multiplet const &) = delete;
61  Multiplet & operator = (Multiplet &&) = delete;
62 
63  // Required functions.
64  bool filter(art::Event & event) override;
65 
66  // Selected optional functions.
67  void beginJob() override;
68  void endJob() override;
69 
70  //calculate the corrected time difference
71  double deltaTCorr(const novaddt::DAQHit& h1, const novaddt::DAQHit& h2);
72  // Check if hit belongs to group
73  bool isHitCloseToSlice(const DAQHit& h, const HitList& frag){
74  //require that candidate hit is close to any hit in slice
75  for(const auto& h2: frag){
76  if(isHitCloseToHit(h,h2))return true;
77  }
78  return false;
79  }
80 
81  bool isHitCloseToHit( const DAQHit& h1, const DAQHit& h2){
82  if (delta<Plane>(h1,h2)>fMaxPlaneSep) return false;
83  if (h1.View()==h2.View())
84  return (delta<Cell>(h1,h2)<=fMaxCellSep)&&
85  (std::abs(deltaTDC(h1,h2))<=fMaxDTxx);
86 
87  else
88  return(!fDoApplyDistCorrection||deltaTCorr(h1,h2)<=fMaxDTxy);
89  }
90 
91 private:
92  template<class HitCollection>
93  void findslices(const HitCollection& inhits, std::unique_ptr< std::vector< novaddt::HitList > >& product);
94 
95  template<class HitCollection>
96  HitCollection ApplyCorrection(const HitCollection& hits);
98 
99 
100 
101 private:
102  // Declare member data here.
104 
106  unsigned fMaxPlaneSep;
107  unsigned fMaxCellSep;
108  unsigned fMaxTDCSep;
109  double fMaxDTxx;
110  double fMaxDTxy;
111  unsigned fMinSliceSize;
112  unsigned fMaxSliceSize;
113  bool fVerbose;
116 
117  // Counters
118  int _nEvents = 0;
119  int _nHitsIn = 0;
120  int _nHitsOut = 0;
121  int _nSlicesOut = 0;
122 };
123 
124 
126  fDetUtils("fd"),
127  fInputTag (p.get< std::string >("InputTag" )),
128  fMaxPlaneSep (p.get< unsigned >("slicer.MaxPlaneSep" )),
129  fMaxCellSep (p.get< unsigned >("slicer.MaxCellSep" )),
130  fMaxTDCSep (p.get< unsigned >("slicer.MaxTDCSep" )),
131  fMaxDTxx (p.get< double >("slicer.MaxDT.xx" )),
132  fMaxDTxy (p.get< double >("slicer.MaxDT.xy" )),
133  fMinSliceSize (p.get< unsigned >("slicer.MinSliceSize" )),
134  fMaxSliceSize (p.get< unsigned >("slicer.MaxSliceSize" )),
135  fVerbose (p.get< bool >("Verbose" )),
136  fDoApplyTimingCorrection (p.get< bool >("ApplyTimingCorrection", false)),
137  fDoApplyDistCorrection (p.get< bool >("ApplyDistCorrection", false))
138  // Initialize member data
139 {
140  std::cout << "--- novaddt::Multiplet begin" << std::endl;
141  std::cout << "\t Input hitlists tag: " << fInputTag << std::endl;
142 
143  std::cout << "\t Max Z separation: " << fMaxPlaneSep << std::endl;
144  std::cout << "\t Max XY separation: " << fMaxCellSep << std::endl;
145  std::cout << "\t Max time separation: " << fMaxTDCSep << std::endl;
146  std::cout << "\t Min slice size: " << fMinSliceSize << std::endl;
147  std::cout << "\t Max slice size: " << fMaxSliceSize << std::endl;
148  std::cout << "\t Apply timing correction: " << fDoApplyTimingCorrection << std::endl;
149  std::cout << "\t Apply distance correction: " << fDoApplyDistCorrection << std::endl;
151  std::cout << "\t \t Max DT for hits in same view: " << fMaxDTxx<< std::endl;
152  std::cout << "\t \t Max DT for hits in XY views: " << fMaxDTxy<< std::endl;
153  }
154  std::cout << "\t Verbose: " << fVerbose << std::endl;
155 
156 
157  // Call appropriate Produces<>() functions
158  produces<std::vector<HitList>>();
159 }
160 //////////////////////////////////////////////////////////////////////
161 
163 {
164  _nEvents++;
165 
166  // Hitlist from event (input)
167  auto hitsHandle=event.getValidHandle<novaddt::HitList>(fInputTag);
168  //make a set from a list
169  novaddt::HitSet allHits(*hitsHandle);
170  //apply correction if needed
172  allHits = ApplyCorrection(allHits);
173 
174  _nHitsIn+=allHits.size();
175 
176  //make a vector of hitlists out of hits after slicing
177  std::unique_ptr< std::vector<HitList> >product (new std::vector<HitList>);
178 
179  // Find slices
180  findslices(allHits,product);
181 
182  _nSlicesOut+=product->size();
183  for (auto slice: *product) _nHitsOut+=slice.size();
184 
185  // Add filtered hitlists to event
186  event.put(std::move(product ));
187  return true;
188 }
189  //make a set from a list
191 {
192  // Implementation of optional member function here.
193 }
194 
195 //////////////////////////////////////////////////////////////////////
197 {
198  if (!fVerbose) return;
199  std::cout << "--- novaddt::Multiplet end " << std::endl;
200  std::cout << "\t# of events: " << _nEvents << std::endl;
201  std::cout << "\t# of hits input: " << _nHitsIn << std::endl;
202  std::cout << "\t# of hits output: " << _nHitsOut << std::endl;
203  std::cout << "\t# of slices output: " << _nSlicesOut << std::endl;
204 }
205 
206 std::ostream& operator << (std::ostream& s, novaddt::DAQHit& hit){
207  s<<"adc="<<hit.ADC().val
208  <<" p="<<hit.Plane().val
209  <<" c="<<hit.Cell().val<<"("<<hit.View().val<<")"
210  <<" t="<<hit.TDC().val;
211  return s;
212 }
213 
214 //////////////////////////////////////////////////////////////////////
216 {
217  double dist1 = fDetUtils.GetReadoutDistance(h1,h2.Cell());
218  double dist2 = fDetUtils.GetReadoutDistance(h2,h1.Cell());
219  double dtc1 = fDetUtils.DistTimeOffset(dist1)+fDetUtils.PigTimeOffset(h1);
220  double dtc2 = fDetUtils.DistTimeOffset(dist2)+fDetUtils.PigTimeOffset(h2);
221  double dtc = (dtc2-dtc1)*fDetUtils.ns_to_TDC();
222  double DT = deltaTDC(h1,h2);
223  /* std::cout<<"dist = "<<dist1<<" , "<<dist2
224  <<"\t dtc = "<<dtc<<";\tDT = "<<DT
225  <<"("<<delta<TDC>(h1,h2)<<")"
226  <<"\t -> "<<std::abs(DT-dtc)<<std::endl;
227  */
228  return std::abs(DT-dtc);
229 }
230 
231 //////////////////////////////////////////////////////////////////////
233  //applying the timing correction
234  double dt = hit.TDC().fraction()-fDetUtils.CellTimeOffset(hit)*fDetUtils.ns_to_TDC();
235  //get int and frac part
236  long dtI = std::floor(dt);
237  double dtF = dt-dtI;
238  //make a TDC object
239  novaddt::TDC tdc(hit.TDC().val + dtI, dtF*100);
240  /* std::cout<<"dt = "<<dt<<" ( "<<dtI<<" + "<<dtF<<" ) "<<std::endl;
241  std::cout<<"tdc0 = "<<hit.TDC().val<<" + "<<hit.TDC().fraction()<<std::endl;
242  std::cout<<"tdc1 = "<<tdc.val<<" + "<<tdc.fraction()<<std::endl;
243  */
244  return DAQHit(hit.DetID(),hit.Chan(), hit.ADC(), tdc);
245 }
246 
247 //////////////////////////////////////////////////////////////////////
248 template<class HitCollection>
249 HitCollection novaddt::Multiplet::ApplyCorrection(const HitCollection& hits)
250 {
251  HitCollection res;
252  //create a callable object that produces new hits
253  using namespace std::placeholders;
254  auto corrector = std::bind(&novaddt::Multiplet::HitCorrection, this, _1);
255  //apply correction to all hits, populating the res collection
256  std::transform( hits.begin(), hits.end(),
257  std::inserter(res,res.end()),
258  corrector);
259 
260  return res;
261 }
262 
263 //////////////////////////////////////////////////////////////////////
264 template<class HitCollection>
265 void novaddt::Multiplet::findslices(const HitCollection & inhits,
266  std::unique_ptr<std::vector<HitList> > & product)
267 {
268  if (!inhits.empty()) {
269 
270  // copy input hitlist to a list that we can start slicing up
271  std::list<DAQHit> hits(inhits.begin(),inhits.end());
272 
273  if (fVerbose) std::cout << "findslices got " << inhits.size()
274  << " hits, working with " << hits.size()<<std::endl;
275  // Loop over all hits
276  while (!hits.empty()) {
277  // this will be the first hit in the slice, put it there
278  HitList frag;
279  frag.emplace_back(hits.front());
280  // remove it from the front of the list, since it's now in a slice
281  hits.pop_front();
282  // point to the first hit in the slice
283  auto lastInSlice = frag.back();
284  auto testhit = hits.begin();
285  while (testhit != hits.end()) {
286  // since hits are time ordered, march along until we get out of time,
287  // with the last hit already in the slice, adding hits that are
288  // also in space with the last hit already in the slice
289  if(std::abs(deltaTDC(*testhit,lastInSlice)) > fMaxTDCSep) break;
290  //we're out of time, this fragment ends
291 
292  if(isHitCloseToSlice(*testhit,frag))
293  {
294  // add this hit to slice
295  frag.emplace_back(*testhit);
296  lastInSlice = frag.back();
297  // remove it from list of potential hits, since a hit can only
298  // belong to one slice (after moving to next hit)
299  auto erasehit = testhit;
300  testhit++;
301  hits.erase(erasehit);
302  }
303  else testhit++;
304 
305  } // end loop over rest of hits
306 
307  // we've finished building a slice, save it to product
308  if (fVerbose) std::cout << "Sliced: " << frag.size() << " hits from "
309  << frag.front().TDC().val << " to "
310  << frag.back().TDC().val <<std::endl;
311  if ((frag.size() <= fMaxSliceSize)
312  && (frag.size() >= fMinSliceSize)){
313  product->emplace_back(frag);
314  }
315  } // loop over full hitlist
316 
317  }
318  // Last slices seem to be noise slices.
319  // Save the slices that don't meet the cuts in here
320  // TODO: right now we're just dropping "noise" hits, more efficient.
321  // Make an empty slice just in case
322  //HitList noise;
323  //product->emplace_back(noise);
324 }
325 //////////////////////////////////////////////////////////////////////
novaddt::utils::DetectorUtils fDetUtils
Multiplet & operator=(Multiplet const &)=delete
HitCollection ApplyCorrection(const HitCollection &hits)
unsigned delta(const novaddt::DAQHit &h1, const novaddt::DAQHit &h2)
Multiplet(fhicl::ParameterSet const &p)
value_type val
Definition: BaseProducts.h:34
void findslices(const HitCollection &inhits, std::unique_ptr< std::vector< novaddt::HitList > > &product)
DAQHit HitCorrection(const DAQHit &hit)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
double deltaTCorr(const novaddt::DAQHit &h1, const novaddt::DAQHit &h2)
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
value_type val
Definition: BaseProducts.h:109
bool isHitCloseToHit(const DAQHit &h1, const DAQHit &h2)
void beginJob() override
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
double DistTimeOffset(double dist) const
double CellTimeOffset(const DAQHit &hit) const
float abs(float number)
Definition: d0nt_math.hpp:39
std::ostream & operator<<(std::ostream &s, novaddt::DAQHit &hit)
bool isHitCloseToSlice(const DAQHit &h, const HitList &frag)
const XML_Char * s
Definition: expat.h:262
Definition: Cand.cxx:23
void hits()
Definition: readHits.C:15
daqchannelmap::dchan const & Chan() const
Definition: DAQHit.h:57
novaddt::ADC const & ADC() const
Definition: DAQHit.h:73
double DT(const novaddt::TDC &t1, const novaddt::TDC &t2)
novaddt::View const & View() const
Definition: DAQHit.h:72
value_type val
Definition: BaseProducts.h:84
double t2
double deltaTDC(const novaddt::TDC &h1, const novaddt::TDC &h2)
double PigTimeOffset(const DAQHit &hit) const
Definition: DetectorUtils.h:69
double GetReadoutDistance(const DAQHit &hit, const novaddt::Cell &cell) const
double fraction() const
Definition: BaseProducts.h:37
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T product(std::vector< T > dims)
int const & DetID() const
Definition: DAQHit.h:75
Definition: event.h:1
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
value_type val
Definition: BaseProducts.h:137
bool filter(art::Event &event) override
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
void endJob() override
value_type val
Definition: BaseProducts.h:65
art::InputTag fInputTag
double T
Definition: Xdiff_gwt.C:5