MichelEfinder_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MichelEfinder
3 // Module Type: producer
4 // File: MichelEfinder_module.cc
5 //
6 // Generated at Sun Jun 19 05:00:16 2016 by Andrey Sheshukov using artmod
7 // from cetpkgsupport v1_10_02.
8 //
9 // Michel electron finding by Matt Strait, separated in a module
10 // For tests
11 //
12 ////////////////////////////////////////////////////////////////////////
13 
21 #include "fhiclcpp/ParameterSet.h"
23 
28 
30 
31 #include <memory>
32 #include <algorithm>
33 
34 namespace novaddt {
35  class MichelEfinder;
36 }
37 
39 public:
40  explicit MichelEfinder(fhicl::ParameterSet const & p);
41  // The destructor generated by the compiler is fine for classes
42  // without bare pointers or other resource use.
43 
44  // Plugins should not be copied or assigned.
45  MichelEfinder(MichelEfinder const &) = delete;
46  MichelEfinder(MichelEfinder &&) = delete;
47  MichelEfinder & operator = (MichelEfinder const &) = delete;
49 
50  // Required functions.
51  void produce(art::Event & e) override;
52 private:
53  bool michelclose(const novaddt::DAQHit &trackend,
54  const novaddt::DAQHit &testcell);
55  void findmichelexclusion( const novaddt::HitList & ends,
56  const novaddt::HitList & allhits,
57  novaddt::HitList & answer);
58 
59 
60 private:
61  const int kTDC_PER_USEC = 64;
62 
65 
66  double fMaxDistance;
69 };
70 
71 
72 // Returns -1 if the testcell is later than the trackend by more than the
73 // time window. Otherwise returns 0 if the testcell is not within the
74 // michel exclusion zone and 1 if it is.
75 
76 
77 bool novaddt::MichelEfinder::michelclose(const novaddt::DAQHit & __restrict__ trackend,
78  const novaddt::DAQHit & __restrict__ testcell)
79 {
80  if(testcell.View()!=trackend.View())return false;
81 
82  const double dplane = 1.5*(trackend.Plane().val - testcell.Plane().val);
83  if(dplane > fMaxDistance) return false;
84 
85  const double dcell = trackend.Cell().val - testcell.Cell().val;
86  if(dcell > fMaxDistance) return false;
87 
88  double dist2 = pow(dplane, 2) + pow(dcell, 2);
89 
90  // Exclude a disc of hits up to N cells = N/1.5 planes away.
91  return dist2 <= fMaxDistance*fMaxDistance;
92 
93  // Now look, this is a very blunt cut, but we could get clever and notice
94  // that Michel hits ought to be in a line connected to the track end and
95  // < 5mus, whereas neutron hits are more random and can be much later, and
96  // B-12 hits, if we care, are 1-2 cells again connected to the track end.
97  // But let's see if this is sufficient before getting fancy.
98 }
99 
100 
101 //a comparison function for searching in time
102 bool cmp_tdc(const novaddt::DAQHit& h, const novaddt::TDC& tdc){
103  return h.TDC()<tdc;
104 }
105 
106 
107 
108 
109 
111  const novaddt::HitList & ends,
112  const novaddt::HitList & allhits,
113  novaddt::HitList & answer
114 )
115 {
116  // Didn't find any track ends, nothing to exclude
117  if(ends.empty()) return;
118 
119  for(const auto & end: ends){
120  //const int endview = end.View().val;
121 
122  //find the time range
123  // Only need to look at the time region around the end
124  auto tmin=end.TDC()+fDTime_min;
125  auto tmax=end.TDC()+fDTime_max;
126  const auto range_beg = std::lower_bound(allhits.begin(),allhits.end(),tmin, cmp_tdc);
127  const auto range_end = std::lower_bound(range_beg, allhits.end(),tmax, cmp_tdc);
128  for(auto hit = range_beg; hit != range_end; hit++){
129  bool res = michelclose(end, *hit);
130  if(res) answer.push_back(*hit);
131  }
132  }
133 }
134 
135 // 2D track: If in the Y view, assume downwards-going. Find the lowest, and
136 // therefore last, hit. If in the X view, can't tell which way the track is
137 // going. Return both ends.
138 static std::vector<novaddt::DAQHit>
139 findtrackends2d(const std::vector<novaddt::DAQHit> & trackhits)
140 {
141  std::vector<novaddt::DAQHit> answer;
142 
143  // Just in case, bail on no track hits
144  if(trackhits.empty()) return answer;
145 
146  // This is a track in the Y view
147  if(trackhits[0].View().val == daqchannelmap::Y_VIEW){
148  // According to doc-11570, the bottom y cell has number zero.
149  novaddt::DAQHit * lastycell = NULL;
150  {
151  int lowest = 1000;
152  for(const auto & hit: trackhits){
153  if(hit.Cell().val < lowest){
154  lowest = hit.Cell().val;
155  lastycell = (novaddt::DAQHit*)&hit;
156  }
157  }
158  }
159 
160  answer.push_back(*lastycell);
161  return answer;
162  }
163 
164  // This is a track in the X view. Find each end.
165  for(int increasingplane = 0; increasingplane < 2; increasingplane++){
166  novaddt::DAQHit * lastxcell = NULL;
167  {
168  int endiest = increasingplane? 0: 10000;
169  for(const auto & hit: trackhits){
170  if(increasingplane ^ (hit.Plane().val < endiest)){
171  endiest = hit.Plane().val;
172  lastxcell = (novaddt::DAQHit*)&hit;
173  }
174  }
175  }
176  answer.push_back(*lastxcell);
177  }
178 
179  return answer;
180 }
181 
182 
183 // 3D track: Assume track is downwards-going. Find lowest, and therefore last,
184 // Y-hit and from the Z-direction of the track, find the corresponding last
185 // X-hit
186 __attribute__((unused)) static std::vector<novaddt::DAQHit>
187 findtrackends3d(const std::vector<novaddt::DAQHit> & trackhits)
188 {
189  std::vector<novaddt::DAQHit> answer;
190 
191  // Just in case, bail on no track hits
192  if(trackhits.empty()) return answer;
193 
194  // According to doc-11570, the bottom y cell has number zero.
195  novaddt::DAQHit * lastycell = NULL;
196  {
197  int lowest = 1000;
198  for(const auto & hit: trackhits){
199  if(hit.View().val != daqchannelmap::Y_VIEW) continue;
200  if(hit.Cell().val < lowest){
201  lowest = hit.Cell().val;
202  lastycell = (novaddt::DAQHit*)&hit;
203  }
204  }
205  }
206 
207  // No y hits, return empty vector
208  if(lastycell == NULL) return answer;
209  answer.push_back(*lastycell);
210 
211  // Find out which way the track is going
212  novaddt::DAQHit * someotheryhit = NULL;
213  for(const auto & hit: trackhits){
214  if(hit.View().val != daqchannelmap::Y_VIEW) continue;
215  if(lastycell != (novaddt::DAQHit*)&hit)
216  someotheryhit = (novaddt::DAQHit*)&hit;
217  }
218 
219  const bool increasingplane = someotheryhit == NULL ||
220  lastycell->Plane().val > someotheryhit->Plane().val;
221 
222  // Find ends in x.
223  novaddt::DAQHit * lastxcell = NULL;
224  {
225  int endinest = increasingplane? 0: 10000;
226  for(const auto & hit: trackhits){
227  if(hit.View().val != daqchannelmap::X_VIEW) continue;
228 
229  if(increasingplane ^ (hit.Plane().val < endinest)){
230  endinest= hit.Plane().val;
231  lastxcell = (novaddt::DAQHit*)&hit;
232  }
233  }
234  }
235 
236  // y hits, but no x hits, return vector of size one.
237  if(lastxcell == NULL) return answer;
238  answer.push_back(*lastxcell);
239 
240  return answer;
241 }
242 
243 
244 
246 :fTracksTag(p.get<std::string>("TracksTag")),
247  fHitsTag(p.get<std::string>("HitsTag")),
248  fMaxDistance(p.get<double>("window.MaxDistance")),
249  fDTime_min(p.get<int>("window.TimeRange.min")*kTDC_PER_USEC),
250  fDTime_max(p.get<int>("window.TimeRange.max")*kTDC_PER_USEC)
251 // Initialize member data here.
252 {
253  produces<novaddt::HitList>();
254  mf::LogInfo("MichelEfinder")<<"config:"
255  <<"\n TracksTag:"<<fTracksTag
256  <<"\n HitsTag: "<<fHitsTag
257  <<"\n window: "
258  <<"\n\t MaxDistance: "<<fMaxDistance
259  <<"\n\t DTime limits: ["<<fDTime_min<<","<<fDTime_max<<"]";
260 }
261 
262 //-------------------------------------------
264 {
265  auto Tracks=e.getValidHandle<std::vector<novaddt::HitList>>(fTracksTag);
267 
268  std::vector<novaddt::DAQHit> trackEnds;
269  trackEnds.reserve(Tracks->size());
270 
271  mf::LogDebug("MichelEfinder")<<"find tracks ends for Ntrk="<<Tracks->size();
272  for(const auto& track: *Tracks){
273  auto ends=findtrackends2d(track);
274  trackEnds.insert(trackEnds.end(),ends.begin(),ends.end());
275  }
276  auto result=std::make_unique<novaddt::HitList>();
277  findmichelexclusion(trackEnds,*Hits,*result);
278 
279  mf::LogDebug("MichelEfinder")<<"found "<<result->size()<<" michel hits"
280  <<" around "<<trackEnds.size()<<" endpoints";
281  e.put(std::move(result));
282 }
283 
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
bool michelclose(const novaddt::DAQHit &trackend, const novaddt::DAQHit &testcell)
constexpr T pow(T x)
Definition: pow.h:75
Definition: event.h:19
void produce(art::Event &e) override
DEFINE_ART_MODULE(TestTMapFile)
bool cmp_tdc(const novaddt::DAQHit &h, const novaddt::TDC &tdc)
Identifier for the Y measuring view of the detector (side)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
__attribute__((unused)) static std
Identifier for the X measuring view of the detector (top)
MichelEfinder(fhicl::ParameterSet const &p)
MichelEfinder & operator=(MichelEfinder const &)=delete
value_type val
Definition: BaseProducts.h:84
static std::vector< novaddt::DAQHit > findtrackends2d(const std::vector< novaddt::DAQHit > &trackhits)
void findmichelexclusion(const novaddt::HitList &ends, const novaddt::HitList &allhits, novaddt::HitList &answer)
Definition: View.py:1
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: structs.h:12
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Float_t e
Definition: plot.C:35
enum BeamMode string