NuMuTrigger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NuMuTrigger
3 // Module Type: filter
4 // File: NuMuTrigger_module.cc
5 //
6 // Generated at Wed May 8 14:35:56 2013 by Matthew Tamsett using artmod
7 // from art v1_02_06.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include "art/Framework/Core/EDFilter.h"
11 #include "art/Framework/Core/ModuleMacros.h"
12 #include "canvas/Persistency/Common/FindOneP.h"
13 #include "canvas/Persistency/Common/FindManyP.h"
14 #include "art/Framework/Principal/Event.h"
15 #include "canvas/Persistency/Common/Assns.h"
16 
17 #include "DDTBaseDataProducts/BaseProducts.h"
18 #include "DDTBaseDataProducts/Track3D.h"
19 #include "DDTBaseDataProducts/HitList.h"
20 #include "DDTBaseDataProducts/TriggerDecision.h"
21 
22 #include "DAQChannelMap/DAQChannelMap.h"
23 #include "DDTUtilities/AssociationUtil.h"
24 
25 #include "TVector3.h"
26 
27 //---------------------------------------------------------------------
28 namespace novaddt {
29  class NuMuTrigger;
30 }
31 //---------------------------------------------------------------------
32 class novaddt::NuMuTrigger : public art::EDFilter {
33  public:
34  explicit NuMuTrigger(fhicl::ParameterSet const & p);
35  virtual ~NuMuTrigger();
36  bool filter(art::Event & e) override;
37  void findBoundingBox(novaddt::HitList & xz_hits, novaddt::HitList & yz_hits,
38  TVector3 & start, TVector3 & end);
39  double n3DCellsToEdge(novaddt::Track3D & track, bool entering);
40  void endJob() override;
41  private:
42  unsigned int _prescale; ///< prescale factor (1 out of every this many passes are issued as real triggers)
43  std::string _trackModuleLabel; ///< label of module making the Tracks
44  std::string _trackInstanceLabel; ///< instance label making the Tracks
45  std::string _sliceModuleLabel; ///< label of module making the HitList
46  std::string _sliceInstanceLabel; ///< instance label making the HitList
47  std::string _trackToSliceInstanceLabel;///< instance of track to slice association produced by track module
48  unsigned int _minPlanesCrossed; ///< minimum planes crossed for a track to be used
49  unsigned int _minCellsPerView; ///< minimum number of cells per view for a good track
50  double _minCellsToEdge; ///< containment cut based on minimum number of cell projected to the edge of the detector
51  double _sliceBoxXMin; ///< lower x slice box fiducial cut
52  double _sliceBoxXMax; ///< upper x slice box fiducial cut
53  double _sliceBoxYMin; ///< lower y slice box fiducial cut
54  double _sliceBoxYMax; ///< upper y slice box fiducial cut
55  double _sliceBoxZMin; ///< lower z slice box fiducial cut
56  double _sliceBoxZMax; ///< upper z slice box fiducial cut
57  double _trackBoxXMin; ///< lower x track box fiducial cut
58  double _trackBoxXMax; ///< upper x track box fiducial cut
59  double _trackBoxYMin; ///< lower y track box fiducial cut
60  double _trackBoxYMax; ///< upper y track box fiducial cut
61  double _trackBoxZMin; ///< lower z track box fiducial cut
62  double _trackBoxZMax; ///< upper z track box fiducial cut
63  double _angularCut; ///< cosmic rejection cut on angle
64  double _enteringCut; ///< minimum plane crossed to be defined as entering
65  double _cosmicPIDCut; ///< cosmic pid (MinCellsToEdge/100)^2 + cos(beam)^2 - (dY_unit*4)^2
66  double _angleVarCut; ///< cosmic rejection cut on angle: cos(beam)^2 - (dY_unit)^2
67  bool _checkProjectionContainment;
68  bool _checkLowZContainment;
69  bool _checkBoxContainment;
70  bool _checkSliceContainment;
71  bool _checkAngle;
72  bool _checkEntering;
73  bool _checkPID;
74  bool _checkAngleVar;
75  std::string _detector; ///< Set detector specific configuration options, such as the location of the detector boundaries
76  ///< available options: NDOS, FD, FD_2DB - ND will be added later
77  bool _verbose;
78  unsigned int _nEvents = 0;
79  unsigned int _nSlices = 0;
80  unsigned int _nTracks = 0;
81  unsigned int _n3DTracks = 0;
82  unsigned int _nPass_slice_containment = 0;
83  unsigned int _nPass_3DTracks = 0;
84  unsigned int _nPass_dX = 0;
85  unsigned int _nPass_dY = 0;
86  unsigned int _nPass_dZ = 0;
87  unsigned int _nPass_projection_containment = 0;
88  unsigned int _nPass_box_containment = 0;
89  unsigned int _nPass_angle = 0;
90  unsigned int _nPass_entering = 0;
91  unsigned int _nPass_cosmic_PID = 0;
92  unsigned int _nPass_angle_var = 0;
93  unsigned int _trigger_counts = 0;
94  unsigned int _prescaled_trigger_counts = 0;
95 };
96 //---------------------------------------------------------------------
97 novaddt::NuMuTrigger::NuMuTrigger(fhicl::ParameterSet const & p)
98 : _prescale (p.get<unsigned int>("prescale" ))
99  , _trackModuleLabel (p.get<std::string >("TrackModuleLabel" ))
100  , _trackInstanceLabel(p.get<std::string >("TrackInstanceLabel", ""))
101  , _sliceModuleLabel (p.get<std::string >("SliceModuleLabel" ))
102  , _sliceInstanceLabel(p.get<std::string >("SliceInstanceLabel", ""))
103  , _trackToSliceInstanceLabel(p.get<std::string >("TrackToSliceInstanceLabel", ""))
104  , _minPlanesCrossed (p.get<unsigned int>("MinPlanesCrossed" ))
105  , _minCellsPerView (p.get<unsigned int>("MinCellsPerView" ))
106  , _minCellsToEdge (p.get<double> ("MinCellsToEdge" ))
107  , _sliceBoxXMin (p.get<double> ("SliceBoxXMin" ))
108  , _sliceBoxXMax (p.get<double> ("SliceBoxXMax" ))
109  , _sliceBoxYMin (p.get<double> ("SliceBoxYMin" ))
110  , _sliceBoxYMax (p.get<double> ("SliceBoxYMax" ))
111  , _sliceBoxZMin (p.get<double> ("SliceBoxZMin" ))
112  , _sliceBoxZMax (p.get<double> ("SliceBoxZMax" ))
113  , _trackBoxXMin (p.get<double> ("TrackBoxXMin" ))
114  , _trackBoxXMax (p.get<double> ("TrackBoxXMax" ))
115  , _trackBoxYMin (p.get<double> ("TrackBoxYMin" ))
116  , _trackBoxYMax (p.get<double> ("TrackBoxYMax" ))
117  , _trackBoxZMin (p.get<double> ("TrackBoxZMin" ))
118  , _trackBoxZMax (p.get<double> ("TrackBoxZMax" ))
119  , _angularCut (p.get<double> ("AngularCut" ))
120  , _enteringCut (p.get<double> ("EnteringCut" ))
121  , _cosmicPIDCut (p.get<double> ("CosmicPID" ))
122  , _angleVarCut (p.get<double> ("AngleVarCut" ))
123  , _checkProjectionContainment (p.get<bool> ("CheckProjectionContainment" ))
124  , _checkLowZContainment (p.get<bool> ("CheckLowZContainment" ))
125  , _checkBoxContainment (p.get<bool> ("CheckBoxContainment" ))
126  , _checkSliceContainment (p.get<bool> ("CheckSliceContainment" ))
127  , _checkAngle (p.get<bool> ("CheckAngle" ))
128  , _checkEntering (p.get<bool> ("CheckEntering" ))
129  , _checkPID (p.get<bool> ("CheckPID" ))
130  , _checkAngleVar (p.get<bool> ("CheckAngleVar" ))
131  , _detector (p.get<std::string> ("Detector" ))
132  , _verbose (p.get<bool> ("Verbose" ))
133 {
134  std::cout << "=== novaddt::NuMuTrigger instantiate" << std::endl;
135  std::cout << "\t prescale: " << _prescale << std::endl;
136  std::cout << "\t TrackModuleLabel: " << _trackModuleLabel << std::endl;
137  std::cout << "\t TrackInstanceLabel: " << _trackInstanceLabel << std::endl;
138  std::cout << "\t SliceModuleLabel: " << _sliceModuleLabel << std::endl;
139  std::cout << "\t SliceInstanceLabel: " << _sliceInstanceLabel << std::endl;
140  std::cout << "\t TrackToSliceInstanceLabel: " << _trackToSliceInstanceLabel << std::endl;
141  std::cout << "\t Cut on slice containment: " << _checkSliceContainment << std::endl;
142  std::cout << "\t - X slice cuts: " << _sliceBoxXMin << " to " << _sliceBoxXMax << std::endl;
143  std::cout << "\t - Y slice cuts: " << _sliceBoxYMin << " to " << _sliceBoxYMax << std::endl;
144  std::cout << "\t - Z slice cuts: " << _sliceBoxZMin << " to " << _sliceBoxZMax << std::endl;
145  std::cout << "\t MinCellsPerView: " << _minCellsPerView << std::endl;
146  std::cout << "\t MinPlanesCrossed: " << _minPlanesCrossed << std::endl;
147  std::cout << "\t Cut on projection containment: " << _checkProjectionContainment << std::endl;
148  std::cout << "\t - Detector: " << _detector << std::endl;
149  std::cout << "\t - MinCellsToEdge: " << _minCellsToEdge << std::endl;
150  std::cout << "\t - low Z containment: " << _checkLowZContainment << std::endl;
151  std::cout << "\t Cut on track box containment: " << _checkBoxContainment << std::endl;
152  std::cout << "\t - X box cuts: " << _trackBoxXMin << " to " << _trackBoxXMax << std::endl;
153  std::cout << "\t - Y box cuts: " << _trackBoxYMin << " to " << _trackBoxYMax << std::endl;
154  std::cout << "\t - Z box cuts: " << _trackBoxZMin << " to " << _trackBoxZMax << std::endl;
155  std::cout << "\t Cut on angle: " << _checkAngle << std::endl;
156  std::cout << "\t - Angular cut: " << _angularCut << std::endl;
157  std::cout << "\t Cut on entering: " << _checkEntering << std::endl;
158  std::cout << "\t - Entering cut: " << _enteringCut << std::endl;
159  std::cout << "\t Cut on pid: " << _checkPID << std::endl;
160  std::cout << "\t - Cosmic PID cut " << _cosmicPIDCut << std::endl;
161  std::cout << "\t Cut on angle variable: " << _checkAngleVar << std::endl;
162  std::cout << "\t - Angle Var cut: " << _angleVarCut << std::endl;
163  std::cout << "\t Verbose: " << _verbose << std::endl;
164 
165  // check that the requested detector is configured
166  assert( (_detector == "NDOS") ||
167  (_detector == "FD") ||
168  (_detector == "FD_2DB")
169  );
170 
171  produces< std::vector<novaddt::TriggerDecision>>();
172  produces< art::Assns<novaddt::TriggerDecision, novaddt::HitList> >();
173 }
174 //---------------------------------------------------------------------
175 novaddt::NuMuTrigger::~NuMuTrigger()
176 {
177  // Clean up dynamic memory and other resources here.
178 }
179 //---------------------------------------------------------------------
180 bool novaddt::NuMuTrigger::filter(art::Event & event)
181 {
182  // Implementation of required member function here.
183  if (_verbose) std::cout<< "=== novaddt::NuMuTrigger filter. Event: "
184  << event.id().event()
185  << std::endl;
186  // make the trigger decision and assignment vectors we are going to write out to the event
187  std::unique_ptr<std::vector<novaddt::TriggerDecision>> trigger_decisions(new std::vector<novaddt::TriggerDecision>);
188  std::unique_ptr<art::Assns<novaddt::TriggerDecision, novaddt::HitList> > assn(new art::Assns<novaddt::TriggerDecision, novaddt::HitList >);
189 
190  // get the slices
191  art::Handle< std::vector<novaddt::HitList> > slices;
192  event.getByLabel(_sliceModuleLabel, _sliceInstanceLabel, slices);
193  if (_verbose) std::cout<< "\tgot " << slices->size() << " slices" << std::endl;
194  _nSlices+=slices->size();
195 
196  // get the tracks for each hit list
197  const std::string the_slices = _trackModuleLabel+":"+_trackToSliceInstanceLabel;
198  art::FindManyP<novaddt::Track3D> find_tracks(slices, event, the_slices);
199  assert(find_tracks.isValid());
200 
201  // loop over slices
202  for(size_t i_slice = 0; i_slice < slices->size(); ++i_slice){
203  // get the hit list
204  novaddt::HitList hit_list = slices->at(i_slice);
205  if (_verbose){
206  std::cout<< "\tslice[" << i_slice << "]: " << hit_list.size() << " hits" << std::endl;
207  // loop over hits, useful for debugging only at the moment
208  for(size_t i_hit = 0; i_hit < hit_list.size(); ++i_hit){
209  novaddt::DAQHit hit(hit_list[i_hit]);
210  std::string view = "y";
211  if( hit.View().val == daqchannelmap::X_VIEW ) view = "x";
212  std::cout<< "\t\thit[" << i_hit
213  << "]: TDC: " << hit.TDC().val
214  << ", ADC: " << hit.ADC().val
215  << ", plane: " << hit.Plane().val
216  << ", " << view
217  << "-cell: " << hit.Cell().val
218  << std::endl;
219  } // end of loop on hits
220  }
221 
222  //---- check the trigger conditions
223 
224  // slice box containment
225  if (_checkSliceContainment){
226  novaddt::HitList xz_hits;
227  novaddt::HitList yz_hits;
228  // loop over hits in each slice
229  for(unsigned int j = 0; j < hit_list.size(); ++j){
230  novaddt::DAQHit hit(hit_list[j]);
231  if( hit.View().val == daqchannelmap::X_VIEW ){
232  xz_hits.push_back(hit);
233  } else {
234  yz_hits.push_back(hit);
235  }
236  } // end of loop on hits
237  TVector3 slice_start(0,0,0);
238  TVector3 slice_end(0,0,0);
239  findBoundingBox(xz_hits, yz_hits, slice_start, slice_end);
240  // x
241  if ((slice_start.X() < _sliceBoxXMin) || (slice_start.X() > _sliceBoxXMax) ||
242  (slice_end.X() < _sliceBoxXMin) || (slice_end.X() > _sliceBoxXMax)){
243  if (_verbose) std::cout << "\t\t - failed due to slice x-containment" << std::endl;
244  continue;
245  }
246  // y
247  if ((slice_start.Y() < _sliceBoxYMin) || (slice_start.Y() > _sliceBoxYMax) ||
248  (slice_end.Y() < _sliceBoxYMin) || (slice_end.Y() > _sliceBoxYMax)){
249  if (_verbose) std::cout << "\t\t - failed due to slice y-containment" << std::endl;
250  continue;
251  }
252  // z
253  if ((slice_start.Z() < _sliceBoxZMin) || (slice_start.Z() > _sliceBoxZMax) ||
254  (slice_end.Z() < _sliceBoxZMin) || (slice_end.Z() > _sliceBoxZMax)){
255  if (_verbose) std::cout << "\t\t - failed due to slice z-containment" << std::endl;
256  continue;
257  }
258  } // end of check on slice containment
259  _nPass_slice_containment++;
260 
261  // find the tracks for this hit list
262  std::vector< art::Ptr<novaddt::Track3D> > tracks = find_tracks.at(i_slice);
263  if (_verbose) std::cout<< "\t - got " << tracks.size() << " tracks" << std::endl;
264  if (!tracks.size() ){
265  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to lack of tracks" << std::endl;
266  continue;
267  }
268  _nTracks+=tracks.size();
269 
270  // get the hit list for each track
271  art::FindOneP<novaddt::HitList> fohl(tracks, event, _trackModuleLabel);
272  assert(fohl.isValid());
273  assert(fohl.size() == tracks.size());
274 
275  double longest_length = -9;
276  unsigned int i_longest = 0;
277  bool track_found = false;
278 
279  for(size_t i_track = 0; i_track < tracks.size(); ++i_track){
280  art::Ptr<novaddt::Track3D> track = tracks[i_track];
281  // get the hit list for this track
282  art::Ptr<novaddt::HitList> this_hit_list = fohl.at(i_track);
283 
284  TVector3 length = track->End() - track->Start();
285  double this_length = std::abs(length.Mag());
286 
287  if (_verbose) std::cout<< "\t\t3D-track[" << i_track
288  << "]\t view: " << track->View()
289  << ", hits: " << this_hit_list->size()
290  << ", X: " << track->Start().X() << " - " << track->End().X()
291  << ", Y: " << track->Start().Y() << " - " << track->End().Y()
292  << ", Z: " << track->Start().Z() << " - " << track->End().Z()
293  << ", length: " << this_length
294  << std::endl;
295 
296  //// loop over hits, useful for debugging only at the moment
297  //if (_verbose){
298  //for(size_t i_hit = 0; i_hit < this_hit_list->size(); ++i_hit){
299  //std::string h_view = "y";
300  //if( this_hit_list->at(i_hit).View().val == daqchannelmap::X_VIEW) h_view = "x";
301  //std::cout<< "\t\t\thit[" << i_hit
302  //<< "]: \tplane: " << this_hit_list->at(i_hit).Plane().val
303  //<< ",\tcell: " << this_hit_list->at(i_hit).Cell().val
304  //<< ",\t view: " << h_view
305  //<< ",\tTDC: " << this_hit_list->at(i_hit).TDC().val
306  //<< ",\tADC: " << this_hit_list->at(i_hit).ADC().val
307  //<< std::endl;
308  //} // end of loop on hits
309  //}
310 
311  // only select 3D tracks
312  if (!track->Is3D()) continue;
313  _n3DTracks++;
314  // select longest track
315  if (this_length > longest_length){
316  longest_length = this_length;
317  i_longest = i_track;
318  track_found = true;
319  }
320  } // end of loop on tracks
321 
322  if (!track_found){
323  if (_verbose) std::cout<< "\t\t - failed the trigger decision due lack of 3D tracks" << std::endl;
324  continue;
325  }
326  _nPass_3DTracks++;
327  novaddt::Track3D track = *tracks[i_longest];
328  // get the hit list for this track
329  art::Ptr<novaddt::HitList> hits_on_track = fohl.at(i_longest);
330 
331  if (_verbose) std::cout<< "\t\tlongest track[" << i_longest
332  << "], hits: " << hits_on_track->size()
333  << ", x: " << track.Start().X() << " - " << track.End().X()
334  << ", y: " << track.Start().Y() << " - " << track.End().Y()
335  << ", z: " << track.Start().Z() << " - " << track.End().Z()
336  << std::endl;
337 
338  // check the trigger conditions
339  // track exists in both views
340  unsigned int x_hits = 0;
341  unsigned int y_hits = 0;
342  for(size_t i_hit = 0; i_hit < hits_on_track->size(); ++i_hit){
343  if( hits_on_track->at(i_hit).View().val == daqchannelmap::X_VIEW){
344  x_hits++;
345  } else {
346  y_hits++;
347  }
348  } // end of loop on hits
349  // dX
350  if (x_hits < _minCellsPerView){
351  if (_verbose) std::cout<< "\t\t - failed the trigger decision due lack of x hits" << std::endl;
352  continue;
353  }
354  _nPass_dX++;
355  // dY
356  if (y_hits < _minCellsPerView){
357  if (_verbose) std::cout<< "\t\t - failed the trigger decision due lack of y hits" << std::endl;
358  continue;
359  }
360  _nPass_dY++;
361  // dZ
362  double dZ = std::fabs(track.End().Z() - track.Start().Z());
363  if (dZ < _minPlanesCrossed){
364  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to dZ" << std::endl;
365  continue;
366  }
367  _nPass_dZ++;
368  // projection containment
369  double minCellsToEdge = n3DCellsToEdge(track, ((!_checkLowZContainment)||_checkEntering));
370  if (_verbose) std::cout<< "\t\t - min distance to any edge: " << minCellsToEdge << std::endl;
371  if (_checkProjectionContainment && (minCellsToEdge < _minCellsToEdge)){
372  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to containment cut" << std::endl;
373  continue;
374  }
375  _nPass_projection_containment++;
376  // box containment
377  if (_checkBoxContainment && (std::min(track.Start().X(), track.End().X()) < _trackBoxXMin)){
378  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to X min containment cut" << std::endl;
379  continue;
380  }
381  if (_checkBoxContainment && (std::min(track.Start().Y(), track.End().Y()) < _trackBoxYMin)){
382  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to Y min containment cut" << std::endl;
383  continue;
384  }
385  if (_checkBoxContainment && (std::min(track.Start().Z(), track.End().Z()) < _trackBoxZMin)){
386  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to Z min containment cut" << std::endl;
387  continue;
388  }
389  if (_checkBoxContainment && (std::max(track.Start().X(), track.End().X()) > _trackBoxXMax)){
390  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to X max containment cut" << std::endl;
391  continue;
392  }
393  if (_checkBoxContainment && (std::max(track.Start().Y(), track.End().Y()) > _trackBoxYMax)){
394  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to Y max containment cut" << std::endl;
395  continue;
396  }
397  if (_checkBoxContainment && (std::max(track.Start().Z(), track.End().Z()) > _trackBoxZMax)){
398  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to Z max containment cut" << std::endl;
399  continue;
400  }
401  _nPass_box_containment++;
402  // angle
403  double scaled_dZ = (dZ * 6.7) / 4.; // cells are longer in z than in x or y
404  double dX = std::fabs( track.End().X() - track.Start().X() );
405  double dY = std::fabs( track.End().Y() - track.Start().Y() );
406  double CosBeam = ( scaled_dZ/ std::sqrt(dX*dX+dY*dY+scaled_dZ*scaled_dZ) );
407  if (_verbose) std::cout<< "\t\t - cos(beam): " << CosBeam << std::endl;
408  if (_checkAngle && (CosBeam < _angularCut)){
409  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to CosBeam cut" << std::endl;
410  continue;
411  }
412  _nPass_angle++;
413  // entering
414  double min_Z = std::min(track.Start().Z(), track.End().Z());
415  if (_verbose) std::cout<< "\t\t - min Z: " << min_Z << std::endl;
416  if(_checkEntering && (min_Z > _enteringCut)){
417  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to Entering cut" << std::endl;
418  continue;
419  }
420  _nPass_entering++;
421  // cosmic PID
422  TVector3 track_3d = TVector3();
423  track_3d.SetXYZ(dX,dY,scaled_dZ);
424  double dY_unit = track_3d.Unit().Y();
425  double cosmic_pid = ((minCellsToEdge/100.)*(minCellsToEdge/100.)) +
426  (CosBeam*CosBeam) -
427  (dY_unit*dY_unit);
428  if (_verbose) std::cout<< "\t\t - cosmic PID: " << cosmic_pid
429  << ", from nCells: " << minCellsToEdge
430  << ", cos(beam): " << CosBeam
431  << ", dY unit: " << dY_unit
432  << std::endl;
433  if (_checkPID && (cosmic_pid < _cosmicPIDCut)){
434  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to cosmic PID cut" << std::endl;
435  continue;
436  }
437  _nPass_cosmic_PID++;
438  // angle variable
439  // 2-D cut in angle var vs n hits space
440  //double hit_ratio = double(hits_on_track->size()) / double(hit_list.size());
441  //if (_checkAngleVar && (hit_ratio > 0.88)){
442  //double angle_var = (CosBeam*CosBeam) - (dY_unit*dY_unit);
443  //if( hit_list.size() >50 &&
444  //hit_list.size()<=135 &&
445  //(angle_var<((1.1/85.)*(hit_list.size()-1.65)))
446  //){
447  //if (_verbose) std::cout<< "\t\t - failed the trigger decision due to angle variable cut" << std::endl;
448  //continue;
449  //}
450  //if( hit_list.size() >135 &&
451  //hit_list.size()<=250 &&
452  //(angle_var<((0.6/115.)*(hit_list.size()-0.6)))
453  //){
454  //if (_verbose) std::cout<< "\t\t - failed the trigger decision due to angle variable cut" << std::endl;
455  //continue;
456  //}
457  //if( hit_list.size() >250 &&
458  //hit_list.size()<=500 &&
459  //(angle_var<((0.1/250.)*(hit_list.size()+0.6)))
460  //){
461  //if (_verbose) std::cout<< "\t\t - failed the trigger decision due to angle variable cut" << std::endl;
462  //continue;
463  //}
464  //if( hit_list.size() >500 &&
465  //(angle_var<0.8)
466  //){
467  //if (_verbose) std::cout<< "\t\t - failed the trigger decision due to angle variable cut" << std::endl;
468  //continue;
469  //}
470  //}
471  // 1-D cut on angle var
472  double angle_var = (CosBeam*CosBeam) - (dY_unit*dY_unit);
473  if (_checkAngleVar && (angle_var < _angleVarCut)){
474  if (_verbose) std::cout<< "\t\t - failed the trigger decision due to angle variable cut" << std::endl;
475  continue;
476  }
477  _nPass_angle_var++;
478  // if we reach here then we've passed the trigger
479  if (_verbose) std::cout<< "\t\t - passed the trigger" << std::endl;
480  _trigger_counts++;
481  // check the prescale
482  if (_trigger_counts%_prescale == _prescale-1) {
483  _prescaled_trigger_counts++;
484  if (_verbose) std::cout<< "\t\t - not prescaled" << std::endl;
485  // get the min and max TDC from the parent slice
486  // and use these to define the trigger window
487  // safest not to assume this is time sorted
488  uint64_t min_tdc = hit_list.front().TDC().val;
489  uint64_t max_tdc = hit_list.back().TDC().val;
490  for(size_t i_hit = 0; i_hit < hit_list.size(); ++i_hit){
491  if (hit_list[i_hit].TDC().val < min_tdc) min_tdc = hit_list[i_hit].TDC().val;
492  if (hit_list[i_hit].TDC().val > max_tdc) max_tdc = hit_list[i_hit].TDC().val;
493  } // end of loop on hits
494 
495  trigger_decisions->push_back(novaddt::TriggerDecision(
496  min_tdc,
497  max_tdc - min_tdc,
498  novaddt::TriggerMask(daqdataformats::trigID::TRIG_ID_DATA_NUMU),
499  _prescale)
500  );
501  // get the Ptr for the slice
502  art::Ptr<novaddt::HitList> slice_ptr(slices, i_slice);
503 
504  // get the Ptr for the trigger decision
505  art::ProductID td_id = this->getProductID<std::vector<novaddt::TriggerDecision> >();
506  art::Ptr<novaddt::TriggerDecision> td_ptr(td_id, trigger_decisions->size()-1, event.productGetter(td_id));
507 
508  // Association
509  assn->addSingle(td_ptr, slice_ptr);
510 
511  if (_verbose) std::cout<< "\t\t - passed the trigger decision, start time: "
512  << trigger_decisions->back().start()
513  << ", duration: " << trigger_decisions->back().duration()
514  << std::endl;
515  } // end of check on prescale
516 
517  } // end of loop on tracks
518 
519  bool goodTrigger = (trigger_decisions->size() > 0);
520  if (_verbose) std::cout << "=== end of novaddt::NuMuTrigger filter. Trigger decision: " << goodTrigger
521  << " from " << trigger_decisions->size()
522  << " passes (associations: " << assn->size()
523  << ")" << std::endl;
524 
525  event.put(std::move(trigger_decisions));
526  event.put(std::move(assn));
527 
528  return goodTrigger;
529 }
530 //---------------------------------------------------------------------
531 void novaddt::NuMuTrigger::findBoundingBox(novaddt::HitList & xz_hits, novaddt::HitList & yz_hits,
532  TVector3 & start, TVector3 & end)
533 {
534  //if (_verbose) std::cout << "finding bounding box, "
535  //<< xz_hits.size() << " xz hits, "
536  //<< yz_hits.size() << " yz hits"
537  //<< std::endl;
538  // find z boundaries
539  // default to the mid point
540  int xz_z_min = 448;
541  int xz_z_max = 448;
542  int xz_x_min = 192;
543  int xz_x_max = 192;
544  if (xz_hits.size() > 0){
545  xz_z_min = xz_hits[0].Plane().val;
546  xz_z_max = xz_hits[0].Plane().val;
547  xz_x_min = xz_hits[0].Cell().val;
548  xz_x_max = xz_hits[0].Cell().val;
549  }
550  for (auto const& hit : xz_hits){
551  if (hit.Plane().val < xz_z_min) xz_z_min = hit.Plane().val;
552  if (hit.Plane().val > xz_z_max) xz_z_max = hit.Plane().val;
553  if (hit.Cell().val < xz_x_min) xz_x_min = hit.Cell().val;
554  if (hit.Cell().val > xz_x_max) xz_x_max = hit.Cell().val;
555  }
556 
557  int yz_z_min = 448;
558  int yz_z_max = 448;
559  int yz_y_min = 192;
560  int yz_y_max = 192;
561  if (yz_hits.size() > 0){
562  yz_z_min = yz_hits[0].Plane().val;
563  yz_z_max = yz_hits[0].Plane().val;
564  yz_y_min = yz_hits[0].Cell().val;
565  yz_y_max = yz_hits[0].Cell().val;
566  }
567  for (auto const& hit : yz_hits){
568  if (hit.Plane().val < yz_z_min) yz_z_min = hit.Plane().val;
569  if (hit.Plane().val > yz_z_max) yz_z_max = hit.Plane().val;
570  if (hit.Cell().val < yz_y_min) yz_y_min = hit.Cell().val;
571  if (hit.Cell().val > yz_y_max) yz_y_max = hit.Cell().val;
572  }
573 
574  int z_min = std::min(xz_z_min,yz_z_min);
575  int z_max = std::max(xz_z_max,yz_z_max);
576 
577  start.SetXYZ( double (xz_x_min),
578  double (yz_y_min),
579  double (z_min));
580  end.SetXYZ( double (xz_x_max),
581  double (yz_y_max),
582  double (z_max));
583  return;
584 }
585 //---------------------------------------------------------------------
586 double novaddt::NuMuTrigger::n3DCellsToEdge(novaddt::Track3D & track, bool entering)
587 {
588  //std::cout<< "--- Find distance to Edge, entering: " << entering << std::endl;
589  //// Debug
590  //std::cout<< "--- track xz: " << xz_track.Start().X()
591  //<< " - " << xz_track.End().X()
592  //<< ", z: " << xz_track.Start().Z()
593  //<< " - " << xz_track.End().Z()
594  //<< std::endl;
595 
596  //std::cout<< "--- track yz: " << yz_track.Start().X()
597  //<< " - " << yz_track.End().X()
598  //<< ", z: " << yz_track.Start().Z()
599  //<< " - " << yz_track.End().Z()
600  //<< std::endl;
601 
602  // get the direction of the track
603  // x
604  double min_X = std::min(track.Start().X(), track.End().X());
605  double max_X = std::max(track.Start().X(), track.End().X());
606  // y
607  double min_Y = std::min(track.Start().Y(), track.End().Y());
608  double max_Y = std::max(track.Start().Y(), track.End().Y());
609  // Z
610  double min_Z = std::min(track.Start().Z(), track.End().Z());
611  double max_Z = std::max(track.Start().Z(), track.End().Z());
612  // direction
613  TVector3 trajectory = TVector3(max_X - min_X,
614  max_Y - min_Y,
615  max_Z - min_Z);
616  TVector3 dir = trajectory.Unit();
617 
618  double xlow, xhigh, ylow, yhigh, zlow, zhigh;
619  if (_detector == "FD"){
620  xlow = 0;
621  xhigh = 383;
622  ylow = 0;
623  yhigh = 383;
624  zlow = 0;
625  zhigh = 895;
626  } else if (_detector == "FD_2DB"){
627  xlow = 0;
628  xhigh = 383;
629  ylow = 0;
630  yhigh = 383;
631  zlow = 0;
632  zhigh = 127;
633  } else {
634  // default to NDOS
635  xlow = 0;
636  xhigh = 63;
637  ylow = 0;
638  yhigh = 95;
639  zlow = 0;
640  zhigh = 198;
641  }
642  // Make sure we're inside the box!
643  // if not set the track edge to match the detector boundary
644  if (min_X<=xlow) min_X = xlow;
645  if (max_X>=xhigh) max_X = xhigh;
646  if (min_Y<=ylow) min_Y = ylow;
647  if (max_Y>=yhigh) max_Y = yhigh;
648  if (min_Z<=zlow) min_Z = zlow;
649  if (max_Z>=zhigh) max_Z = zhigh;
650 
651  // If we're looking for entering tracks, then we don't want to check containment to the lower Z edge:
652  if (entering) zlow = -999;
653 
654  // Compute the distances to the x/y/z walls
655  double dx = 99.E99;
656  double dy = 99.E99;
657  double dz = 99.E99;
658  // we will assume that the track could be going in either direction
659  dx = std::min( ((xhigh - max_X)/dir.X()),
660  ((min_X - xlow) /dir.X())
661  );
662  dy = std::min( ((yhigh - max_Y)/dir.Y()),
663  ((min_Y - ylow) /dir.Y())
664  );
665  dz = std::min( ((zhigh - max_Z)/dir.Z()),
666  std::fabs((min_Z - zlow) /dir.Z()) // introduce fabs here so that we can offset the z_low for entering tracks
667  );
668  // Choose the shortest distance
669  double d = 0.0;
670  if (dx<dy && dx<dz) d = dx;
671  else if (dy<dz && dy<dx) d = dy;
672  else if (dz<dx && dz<dy) d = dz;
673 
674  // Make the step
675  TVector3 distance = TVector3(dir.X()*d, dir.Y()*d, dir.Z()*d);
676 
677  //// Debug
678  //if (_verbose) std::cout<< "--- combined track x: " << min_X
679  //<< " - " << max_X
680  //<< ", y: " << min_Y
681  //<< " - " << max_Y
682  //<< ", z: " << min_Z
683  //<< " - " << max_Z
684  //<< std::endl;
685 
686  //if (_verbose) std::cout<< "--- unit dir x: " << dir.X()
687  //<< ", y: " << dir.Y()
688  //<< ", z: " << dir.Z()
689  //<< std::endl;
690 
691  //if (_verbose) std::cout<< "--- dx: " << dx
692  //<< ", dy: " << dy
693  //<< ", dz: " << dz
694  //<< " => d: " << d
695  //<< std::endl;
696 
697  //if (_verbose) std::cout<< "--- distance x: " << distance.X()
698  //<< ", y: " << distance.Y()
699  //<< ", z: " << distance.Z()
700  //<< ", mag: " << distance.Mag()
701  //<< std::endl;
702  //assert(distance.Mag() > 1.);
703  return distance.Mag();
704 }
705 //---------------------------------------------------------------------
706 void novaddt::NuMuTrigger::endJob()
707 {
708  std::cout << "--- novaddt::NuMuTrigger endJob" << std::endl;
709  std::cout << "\tNumber of slices seen: " << _nSlices << std::endl;
710  std::cout << "\t - pass slice containment: " << _nPass_slice_containment << std::endl;
711  std::cout << "\tNumber of tracks seen: " << _nTracks << std::endl;
712  std::cout << "\tNumber of 3D tracks seen: " << _n3DTracks << std::endl;
713  std::cout << "\t - have a 3D track: " << _nPass_3DTracks << std::endl;
714  std::cout << "\t - pass dX: " << _nPass_dX << std::endl;
715  std::cout << "\t - pass dY: " << _nPass_dY << std::endl;
716  std::cout << "\t - pass dZ: " << _nPass_dZ << std::endl;
717  std::cout << "\t - pass projection containment: " << _nPass_projection_containment << std::endl;
718  std::cout << "\t - pass box containment: " << _nPass_box_containment << std::endl;
719  std::cout << "\t - pass angle: " << _nPass_angle << std::endl;
720  std::cout << "\t - pass entering: " << _nPass_entering << std::endl;
721  std::cout << "\t - pass cosmic PID: " << _nPass_cosmic_PID << std::endl;
722  std::cout << "\t - pass angle variable: " << _nPass_angle_var << std::endl;
723  std::cout << "\t - pass the trigger: " << _trigger_counts << std::endl;
724  std::cout << "\t after prescale: " << _prescaled_trigger_counts << std::endl;
725 }
726 //---------------------------------------------------------------------
727 DEFINE_ART_MODULE(novaddt::NuMuTrigger)
novaddt::NuMuTrigger::_trackBoxYMin
double _trackBoxYMin
lower y track box fiducial cut
Definition: NuMuTrigger_module.cc:59
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:415
novaddt::NuMuTrigger::_minCellsPerView
unsigned int _minCellsPerView
minimum number of cells per view for a good track
Definition: NuMuTrigger_module.cc:49
novaddt::NuMuTrigger::_sliceInstanceLabel
std::string _sliceInstanceLabel
instance label making the HitList
Definition: NuMuTrigger_module.cc:46
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
novaddt::TDC::val
value_type val
Definition: BaseProducts.h:34
novaddt::NuMuTrigger::_sliceBoxXMax
double _sliceBoxXMax
upper x slice box fiducial cut
Definition: NuMuTrigger_module.cc:52
febshutoff_auto.end
end
Definition: febshutoff_auto.py:54
novaddt::NuMuTrigger::NuMuTrigger
NuMuTrigger(fhicl::ParameterSet const &p)
Definition: NuMuTrigger_module.cc:97
novaddt::DAQHit::Plane
novaddt::Plane const & Plane() const
Definition: DAQHit.h:70
novaddt::DAQHit::TDC
novaddt::TDC const & TDC() const
Definition: DAQHit.h:74
novaddt::HitList
std::vector< DAQHit > HitList
Definition: HitList.h:15
novaddt::NuMuTrigger::_minCellsToEdge
double _minCellsToEdge
containment cut based on minimum number of cell projected to the edge of the detector ...
Definition: NuMuTrigger_module.cc:50
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
novaddt::Cell::val
value_type val
Definition: BaseProducts.h:109
track
Definition: event.h:19
novaddt::NuMuTrigger::_angularCut
double _angularCut
cosmic rejection cut on angle
Definition: NuMuTrigger_module.cc:63
site_stats_from_log.get
get
Definition: site_stats_from_log.py:189
novaddt::NuMuTrigger::n3DCellsToEdge
double n3DCellsToEdge(novaddt::Track3D &track, bool entering)
Definition: NuMuTrigger_module.cc:586
novaddt::NuMuTrigger::_sliceModuleLabel
std::string _sliceModuleLabel
label of module making the HitList
Definition: NuMuTrigger_module.cc:45
novaddt::NuMuTrigger::_trackBoxXMin
double _trackBoxXMin
lower x track box fiducial cut
Definition: NuMuTrigger_module.cc:57
DEFINE_ART_MODULE
DEFINE_ART_MODULE(TestTMapFile)
distance
unsigned distance(const T &t1, const T &t2)
Definition: Multiplet_module.cc:33
novaddt::NuMuTrigger::_trackBoxZMax
double _trackBoxZMax
upper z track box fiducial cut
Definition: NuMuTrigger_module.cc:62
std::abs
float abs(float number)
Definition: d0nt_math.hpp:39
novaddt::Track3D::End
TVector3 const & End() const
Definition: Track3D.h:46
novaddt::NuMuTrigger::_nPass_cosmic_PID
unsigned int _nPass_cosmic_PID
Definition: NuMuTrigger_module.cc:91
febshutoff_auto.start
start
Definition: febshutoff_auto.py:53
novaddt::NuMuTrigger::filter
bool filter(art::Event &e) override
Definition: NuMuTrigger_module.cc:180
novaddt::NuMuTrigger::findBoundingBox
void findBoundingBox(novaddt::HitList &xz_hits, novaddt::HitList &yz_hits, TVector3 &start, TVector3 &end)
Definition: NuMuTrigger_module.cc:531
dy
double dy[NP][NC]
Definition: checkMeasured_Target.C:19
demo0.length
length
Definition: demo0.py:21
dx
double dx[NP][NC]
Definition: checkMeasured_Target.C:18
dz
double dz[NP][NC]
Definition: checkMeasured_Target.C:20
novaddt::NuMuTrigger::_sliceBoxXMin
double _sliceBoxXMin
lower x slice box fiducial cut
Definition: NuMuTrigger_module.cc:51
novaddt::DAQHit::ADC
novaddt::ADC const & ADC() const
Definition: DAQHit.h:73
novaddt::DAQHit::View
novaddt::View const & View() const
Definition: DAQHit.h:72
daqchannelmap::X_VIEW
Identifier for the X measuring view of the detector (top)
Definition: DAQChannelMapConstants.h:118
d
Float_t d
Definition: plot.C:236
novaddt::NuMuTrigger::_prescale
unsigned int _prescale
prescale factor (1 out of every this many passes are issued as real triggers)
Definition: NuMuTrigger_module.cc:42
novaddt::Track3D::Start
TVector3 const & Start() const
Definition: Track3D.h:45
calib::j
const double j
Definition: BetheBloch.cxx:29
novaddt::Plane::val
value_type val
Definition: BaseProducts.h:84
novaddt::Track3D::View
int const & View() const
Definition: Track3D.h:44
novaddt::Track3D::Is3D
bool const & Is3D() const
Definition: Track3D.h:43
novaddt::NuMuTrigger::_trackInstanceLabel
std::string _trackInstanceLabel
instance label making the Tracks
Definition: NuMuTrigger_module.cc:44
novaddt::NuMuTrigger::_nPass_box_containment
unsigned int _nPass_box_containment
Definition: NuMuTrigger_module.cc:88
om::cout
OStream cout
Definition: OStream.cxx:6
novaddt::NuMuTrigger::_trackModuleLabel
std::string _trackModuleLabel
label of module making the Tracks
Definition: NuMuTrigger_module.cc:43
novaddt::NuMuTrigger::_trackBoxZMin
double _trackBoxZMin
lower z track box fiducial cut
Definition: NuMuTrigger_module.cc:61
Event.h
novaddt::NuMuTrigger::_prescaled_trigger_counts
unsigned int _prescaled_trigger_counts
Definition: NuMuTrigger_module.cc:94
novaddt::NuMuTrigger::_nPass_slice_containment
unsigned int _nPass_slice_containment
Definition: NuMuTrigger_module.cc:82
dir
TDirectory * dir
Definition: macro.C:5
hit
Definition: structs.h:12
novaddt::NuMuTrigger::_trackBoxYMax
double _trackBoxYMax
upper y track box fiducial cut
Definition: NuMuTrigger_module.cc:60
novaddt::NuMuTrigger::_trackToSliceInstanceLabel
std::string _trackToSliceInstanceLabel
instance of track to slice association produced by track module
Definition: NuMuTrigger_module.cc:47
novaddt::View::val
value_type val
Definition: BaseProducts.h:137
novaddt::NuMuTrigger::_cosmicPIDCut
double _cosmicPIDCut
cosmic pid (MinCellsToEdge/100)^2 + cos(beam)^2 - (dY_unit*4)^2
Definition: NuMuTrigger_module.cc:65
novaddt::NuMuTrigger::_nPass_projection_containment
unsigned int _nPass_projection_containment
Definition: NuMuTrigger_module.cc:87
novaddt::DAQHit::Cell
novaddt::Cell const & Cell() const
Definition: DAQHit.h:71
novaddt::NuMuTrigger::_enteringCut
double _enteringCut
minimum plane crossed to be defined as entering
Definition: NuMuTrigger_module.cc:64
ana::assert
assert(nhit_max >=nhit_nbins)
novaddt::ADC::val
value_type val
Definition: BaseProducts.h:65
novaddt::NuMuTrigger::_sliceBoxZMin
double _sliceBoxZMin
lower z slice box fiducial cut
Definition: NuMuTrigger_module.cc:55
novaddt::NuMuTrigger::_sliceBoxYMax
double _sliceBoxYMax
upper y slice box fiducial cut
Definition: NuMuTrigger_module.cc:54
std::min
T min(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:405
e
Float_t e
Definition: plot.C:35
novaddt::NuMuTrigger::_angleVarCut
double _angleVarCut
cosmic rejection cut on angle: cos(beam)^2 - (dY_unit)^2
Definition: NuMuTrigger_module.cc:66
novaddt::NuMuTrigger::_minPlanesCrossed
unsigned int _minPlanesCrossed
minimum planes crossed for a track to be used
Definition: NuMuTrigger_module.cc:48
novaddt::NuMuTrigger::_sliceBoxYMin
double _sliceBoxYMin
lower y slice box fiducial cut
Definition: NuMuTrigger_module.cc:53
art::Ptr
Definition: fwd.h:28
novaddt::NuMuTrigger::_trackBoxXMax
double _trackBoxXMax
upper x track box fiducial cut
Definition: NuMuTrigger_module.cc:58
novaddt::NuMuTrigger::_sliceBoxZMax
double _sliceBoxZMax
upper z slice box fiducial cut
Definition: NuMuTrigger_module.cc:56
string
enum BeamMode string
Generated on Mon Jan 25 2021 14:12:18 for NOvA by   doxygen 1.8.11