MichelETrigger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MichelETrigger
3 // Module Type: filter
4 // File: MichelETrigger_module.cc
5 //
6 // Generated at Wed Jul 19 10:03:20 2017 by Matthew Judah using artmod
7 // from cetpkgsupport v1_10_02.
8 ////////////////////////////////////////////////////////////////////////
9 
10 
18 #include "fhiclcpp/ParameterSet.h"
26 
28 
29 //DDT Includes
37 
38 //C++ Includes
39 #include <memory>
40 #include <sstream>
41 #include <string>
42 #include <TMath.h>
43 #include <algorithm>
44 #include <numeric>
45 #include <iomanip>
46 #include <iostream>
47 #include <list>
48 #include <vector>
49 #include <cmath>
50 
51 //ROOT Includes
52 #include "TH1F.h"
53 #include "TH2F.h"
54 #include "TVector3.h"
55 
56 
57 #include <memory>
58 
59 namespace novaddt {
60  class MichelETrigger;
61 }
62 
64 public:
65 
66  int nEvents = 0;
67  int nTracks = 0;
68  int nStoppedTracks = 0;
69  int nClusters1 = 0;
70  int nMichels = 0;
71  int nTrig = 0;
72 
73  explicit MichelETrigger(fhicl::ParameterSet const & p);
74  // The destructor generated by the compiler is fine for classes
75  // without bare pointers or other resource use.
76 
77  // Plugins should not be copied or assigned.
78  MichelETrigger(MichelETrigger const &) = delete;
79  MichelETrigger(MichelETrigger &&) = delete;
80  MichelETrigger & operator = (MichelETrigger const &) = delete;
82 
83  // Required functions.
84  bool filter(art::Event & e) override;
85 
86  float GetSmartPrescale(int xcell, int ycell);
87  bool GetSmartPrescaleDecision(double prescale);
88  // Selected optional functions.
89 
90  void endJob() override;
91  void reconfigure(fhicl::ParameterSet const & p);
92 
93 private:
94  const int kTDC_PER_USEC = 64;
95  const double cw = 3.97;
96  const double pw = 6.654;
97 
98  // Declare member data here.
99  unsigned _prescale;
106  double fDistMuon;
107  double fDistMichel;
108  double fMaxE;
118 };
119 
120 
122 :
123  _prescale (p.get<unsigned>("prescale") ),
124  _tracklabel (p.get<std::string>("TrackLabel") ),
125  _hitslabel (p.get<std::string>("HitsTag") ),
126  fTimeUpperTrack (p.get<double>("time1_u")*kTDC_PER_USEC),
127  fTimeLowerTrack (p.get<double>("time1_l")*kTDC_PER_USEC),
128  fTimeUpperCluster (p.get<double>("time2_u")*kTDC_PER_USEC),
129  fTimeLowerCluster (p.get<double>("time2_l")*kTDC_PER_USEC),
130  fDistMuon (p.get<double>("MaxDistance") ),
131  fDistMichel (p.get<double>("DistanceMichel") ),
132  fMaxE (p.get<double>("EnergyCutOff") ),
133  fFiducial_Z_min (p.get<int>("fiducial.ZetRange.min") ),
134  fFiducial_Z_max (p.get<int>("fiducial.ZetRange.max") ),
135  fFiducial_X_min (p.get<int>("fiducial.XRange.min") ),
136  fFiducial_X_max (p.get<int>("fiducial.XRange.max") ),
137  fFiducial_Y_min (p.get<int>("fiducial.YRange.min") ),
138  fFiducial_Y_max (p.get<int>("fiducial.YRange.max") ),
139  fapplySmartPrescale(p.get<bool>("applySmartPrescale") ),
140  fadditionalPrescale(p.get<double>("additionalPrescale") )
141 {
142  produces<std::vector<TriggerDecision>>();
143 }
144 
145 bool cmp_tdc(const novaddt::DAQHit& h, const novaddt::TDC& tdc){
146  return h.TDC() < tdc;
147 }
148 
149 float novaddt::MichelETrigger::GetSmartPrescale(int xcell, int ycell)
150 {
151  //m and b come from the conversion of 383 cells (X/Y dimension of FD)
152  //to 95 cells
153  float m = 62./246.;
154  float b = 14 - 54*m;
155  float xcell_convert = m*xcell + b;
156  float ycell_convert = m*ycell + b;
157  if(xcell_convert > 94) xcell_convert = 94;
158  if(ycell_convert > 94) ycell_convert = 94;
159  return prescale_table[floor(xcell_convert)][floor(ycell_convert)]
161 }
162 
163 
165  {
166  // Is it good enough for small prescale factors ~1/100?
167  double r = (double)rand() / (double)RAND_MAX;
168  return r < prescale;
169  }
170 
171 
173 {
174  nEvents++;
175 
176  std::unique_ptr<std::vector<TriggerDecision>>
177  trigger_decisions(new std::vector<TriggerDecision>);
178  bool isTriggering = false;
179 
181  e.getByLabel(_tracklabel, Tracks);
182 
183  //get the hit list for each track, use this to calculate mean track time
185  assert(fohl.isValid());
186  assert(fohl.size() == Tracks->size());
187 
188  //Grab all time sorted hits
189  auto cellHits = e.getValidHandle<novaddt::HitList>(_hitslabel);
190 
191  //This will store the Michel Clusters id'd at the end of a track
192  std::vector<std::vector<novaddt::DAQHit>> result;
193 
194  //Loop over all the tracks in this event
195  for(uint iii = 0; iii < Tracks->size(); ++iii){
196  nTracks++;
197  //Grab track
198  novaddt::Track3D track = Tracks->at(iii);
199 
200  //Make Sure track is 3D
201  if(!track.Is3D()) continue;
202 
203  TVector3 trk_end[2];
204  trk_end[0] = track.Start();
205  trk_end[1] = track.End();
206 
207  for(int which_end = 0; which_end < 2; which_end++){
208  if( (trk_end[which_end].X() < fFiducial_X_min) ||
209  (trk_end[which_end].X() > fFiducial_X_max)) continue;
210  if( (trk_end[which_end].Y() < fFiducial_Y_min) ||
211  (trk_end[which_end].Y() > fFiducial_Y_max)) continue;
212  if( (trk_end[which_end].Z() < fFiducial_Z_min) ||
213  (trk_end[which_end].Z() > fFiducial_Z_max)) continue;
214 
215  nStoppedTracks++;
216 
217  //Track stops in the detector. Start looking for the Michel Cluster
218  std::vector<novaddt::DAQHit> meCluster;
219 
220  //Get Average time of the stopped track
221  //get the hit list for this track
222  art::Ptr<novaddt::HitList> this_hit_list = fohl.at(iii);
223 
224  double avg_t = 0;
225  double totW = 0;
226  double min_TDC = 0;
227  double max_TDC = 0;
228 
229  for(uint trackHits = 0; trackHits < this_hit_list->size(); trackHits++){
230  if(trackHits == 0)
231  min_TDC = this_hit_list->at(trackHits).TDC().val;
232  if( this_hit_list->at(trackHits).ADC().val < 50) continue;
233  if( this_hit_list->at(trackHits).TDC().val < min_TDC)
234  min_TDC = this_hit_list->at(trackHits).TDC().val;
235  double w = this_hit_list->at(trackHits).ADC().val;
236  avg_t += this_hit_list->at(trackHits).TDC().val * w;
237  totW += w;
238  }// loop over track hits
239 
240  avg_t = avg_t/totW;
241 
242 
243  //Get start and end cell hits that fall in the time window
244  const double firstT = (avg_t+fTimeLowerTrack);
245  const double lastT = (avg_t+fTimeUpperTrack);
246 
247  std::vector<novaddt::DAQHit> daqHits = *(cellHits);
248 
249 
250  const auto cellHitsBegin = std::lower_bound(daqHits.begin(),
251  daqHits.end(), firstT,
252  cmp_tdc);
253  const auto cellHitsEnd = std::lower_bound(cellHitsBegin, daqHits.end(),
254  lastT, cmp_tdc);
255 
256  //std::cout << avg_t << "\t" << firstT << "\t" << lastT << std::endl;
257 
258  //Hold onto potential Michel Cluster hits
259  std::vector<novaddt::DAQHit> meList;
260  double michel_pos[3] = {0,0,0};
261  int seed_plane_pos = -5;
262  double seed_time = 0;
263  double maxE = 0;
264  double distance_min = 1000;
265  for(auto cellHitsIter = cellHitsBegin; cellHitsIter != cellHitsEnd; ++cellHitsIter)
266  {
267  novaddt::DAQHit c = *cellHitsIter;
268  if(c.ADC() < 50) continue;
269 
270  // std::cout << c.TDC().val << std::endl;
271 
272  double cellpos[3] = {0,0,0};
273  cellpos[0] = c.Cell().val;
274  cellpos[1] = c.Cell().val;
275  cellpos[2] = c.Plane().val;
276 
277  if(c.View().val == daqchannelmap::X_VIEW)
278  cellpos[1] = trk_end[which_end].Y();
279  else
280  cellpos[0] = trk_end[which_end].X();
281 
282  TVector3 distance;
283  distance.SetXYZ((cellpos[0] - trk_end[which_end].X())*cw,
284  (cellpos[1] - trk_end[which_end].Y())*cw,
285  (cellpos[2] - trk_end[which_end].Z())*pw);
286 
287  double dist = distance.Mag();
288 
289  if(dist < fDistMuon && c.ADC().val < fMaxE){
291  meList.push_back(hit);
292  //std::cout << (hit.TDC().val - avg_t)/kTDC_PER_USEC << "\t" <<
293  // dist << "\t" << hit.ADC().val << std::endl;
294  if(hit.ADC() > maxE && dist < distance_min){
295  seed_time = hit.TDC().val;
296  seed_plane_pos = cellpos[2];
297  maxE = hit.ADC().val;
298  michel_pos[0] = cellpos[0];
299  michel_pos[1] = cellpos[1];
300  michel_pos[2] = cellpos[2];
301  }
302  }
303  }// First loop over cell hits
304 
305  //Didn't find any me candidates don't bother looking any further
306  if(seed_plane_pos < 0) continue;
307  nClusters1++;
308 
309  max_TDC = seed_time;// fTimeUpperCluster;
310  double dist2 = 1000;
311  auto iter = meList.begin();
312  while( iter != meList.end()){
313  novaddt::DAQHit c = *iter;
314 
315  double plane = c.Plane().val;
316 
317  dist2 = abs( plane - seed_plane_pos)*pw;
318 
319  //Make second timing cut
320  //Select only cell hits that happen
321  //within a certain time of the michel e candidate hit
322 
323  if( ((seed_time - fTimeLowerCluster) < (*iter).TDC().val ) &&
324  ((seed_time + fTimeUpperCluster) > (*iter).TDC().val ) &&
325  (dist2 < fDistMichel)){
326  meCluster.push_back(*iter);
327  if ((*iter).TDC().val > max_TDC)
328  max_TDC = (*iter).TDC().val;
329  }
330  iter++;
331  }//end second loop on michel hits
332 
333  //Number Cut: Only Keep Michel Cluster at least one cell
334  //hit and has cell hits in both views
335  int x_hits = 0;
336  int y_hits = 0;
337  double TotalADC = 0;
338  if(meCluster.size() > 0){
339  for(uint ihit = 0; ihit < meCluster.size(); ihit++){
340  if(meCluster[ihit].View().val == daqchannelmap::X_VIEW)
341  x_hits++;
342  else
343  y_hits++;
344  TotalADC += meCluster[ihit].ADC().val;
345  }
346  }
347 
348  if( x_hits > 0 && y_hits > 0){
349  result.push_back(meCluster);
350  _trigger_counts++;
351  nMichels++;
352  //std::cout << michel_pos[0] << "\t" << michel_pos[1] << "\t" <<
353  //michel_pos[2] << "\t" << std::endl;
354  }
355  if((_trigger_counts%_prescale == _prescale-1) && (x_hits > 0) && (y_hits > 0)) {
357  double prescaleValue = GetSmartPrescale(michel_pos[0],michel_pos[1]);
358  bool shouldTrigger = GetSmartPrescaleDecision(prescaleValue);
359  if(shouldTrigger){
360  trigger_decisions->emplace_back(min_TDC, max_TDC - min_TDC,
362  _prescale);
363  nTrig++;
364  }
365  }
366  else{
367  trigger_decisions->emplace_back(min_TDC, max_TDC - min_TDC, daqdataformats::trigID::TRIG_ID_DATA_MICHEL, _prescale);
368  nTrig++;
369  }
370  }
371  }// loop over ends of track
372  }//end loop over tracks
373 
374  if(result.size() > 0 && trigger_decisions->size() > 0) isTriggering = true;
375  e.put(std::move(trigger_decisions));
376  return isTriggering;
377 } //end Filter()
378 
379 
380 
382 {
383  std::cout << "=== novaddt::MichelETrigger endJob" << std::endl;
384  std::cout << "\tNumber of events: " << nEvents << std::endl;
385  std::cout << "\tNumber of tracks: " << nTracks << std::endl;
386  std::cout << "\tNumber of stopped tracks: " << nStoppedTracks
387  << std::endl;
388  std::cout << "\tNumber of cluster seeds: " << nClusters1 << std::endl;
389  std::cout << "\tNumber of Michels tagged: " << nMichels << std::endl;
390  std::cout << "\tNumber of Michels Triggered: " << nTrig << std::endl;
391 }
392 
393 
394 
396 {
397  // Implementation of optional member function here.
398 }
399 
value_type val
Definition: BaseProducts.h:34
void reconfigure(fhicl::ParameterSet const &p)
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
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 GetSmartPrescaleDecision(double prescale)
value_type val
Definition: BaseProducts.h:109
Definition: event.h:19
void abs(TH1 *hist)
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
TVector3 const & End() const
Definition: Track3D.h:46
Float_t Y
Definition: plot.C:38
double dist
Definition: runWimpSim.h:113
Float_t Z
Definition: plot.C:38
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
bool filter(art::Event &e) override
bool cmp_tdc(const novaddt::DAQHit &h, const novaddt::TDC &tdc)
novaddt::ADC const & ADC() const
Definition: DAQHit.h:73
novaddt::View const & View() const
Definition: DAQHit.h:72
Identifier for the X measuring view of the detector (top)
float GetSmartPrescale(int xcell, int ycell)
TVector3 const & Start() const
Definition: Track3D.h:45
MichelETrigger & operator=(MichelETrigger const &)=delete
value_type val
Definition: BaseProducts.h:84
bool const & Is3D() const
Definition: Track3D.h:43
Definition: View.py:1
OStream cout
Definition: OStream.cxx:6
const std::vector< std::vector< float > > prescale_table
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
const hit & b
Definition: hits.cxx:21
value_type val
Definition: BaseProducts.h:137
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
value_type val
Definition: BaseProducts.h:65
MichelETrigger(fhicl::ParameterSet const &p)
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
double maxE
Definition: plot_hist.C:8
Float_t e
Definition: plot.C:35
Float_t X
Definition: plot.C:38
Float_t w
Definition: plot.C:20
Definition: fwd.h:28
unsigned int uint