UpMuTestTrigger_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////
2 // Class: UpMuTestTrigger
3 // Module Type: filter
4 // File: UpMuTestTrigger_module.cc
5 //
6 // Generated at Tues July 17 2013 by Will Henderson
7 // Based on NuMuTestTrigger (Matthew Tamsett)
8 //////////////////////////////////////////////////////////////////////
9 
10 // NOvAART includes
16 
17 // NOvADDT includes
24 
25 // NOvADAQ includes
26 #include "DAQChannelMap/DAQChannelMap.h"
28 
29 // ROOT includes
30 #include "TVector2.h"
31 #include "TVector3.h"
32 #include "TH1F.h"
33 #include "TH2F.h"
34 #include "TFile.h"
35 
36 //////////////////////////////////////////////////////////////////////
37 namespace novaddt {
38  class UpMuTestTrigger;
39 }
40 //////////////////////////////////////////////////////////////////////
42 public:
43  explicit UpMuTestTrigger(fhicl::ParameterSet const & p);
44  virtual ~UpMuTestTrigger();
45  bool filter(art::Event & e) override;
46  void endJob() override;
47 private:
48  unsigned int _prescale; ///< prescale factor (1 out of every this many passes are issued as real triggers)
49  std::string _TrackModuleLabel; ///< label of module making the Tracks
50  std::string _TrackInstanceLabel; ///< instance label making the Tracks
51  std::string _SliceModuleLabel; ///< label of module making the HitList
52  std::string _SliceInstanceLabel; ///< instance label making the HitList
53 
54  // Minimum track hit cuts
55  unsigned int fMinPlanes; // minimum planes crossed (Z)
56  unsigned int fMinCellsX; // minimum number of cells in X view
57  unsigned int fMinCellsY; // minimum number of cells in Y view
58 
59  // Containment cuts
60  float fMinX; // lower x fiducial cut
61  float fMaxX; // upper x fiducial cut
62  float fMinY; // lower y fiducial cut
63  float fMaxY; // upper y fiducial cut
64  unsigned int fMinZ; // lower z fiducial cut
65  unsigned int fMaxZ; // upper z fiducial cut
66 
67  // Angle cut
68  double fAngleCut; // accept 30 deg vertical cone (default-- cos angle is angle from horizontal)
69 
70  // Direction cut
71  double fDirectionCut; // accept upward slope
72 
73  // Counters
74  unsigned int _nEvents =0;
75  unsigned int _nTracks =0;
76  unsigned int _nPass_dhorizontal =0;
77  unsigned int _nPass_dvertical =0;
78  unsigned int _nPass_through =0;
79  unsigned int _nPass_stopped =0;
80  unsigned int _nPass_fully =0;
81  unsigned int _nPass_angle =0;
82  unsigned int _nPass_direction =0;
83  unsigned int _trigger_counts = 0;
84  unsigned int _prescaled_trigger_counts = 0;
85 
86  // Switches
87  bool fcheckMinHits; // boolean switch for cut on minimum track hits
88  bool fcheckContainment; // boolean switch for cut on containment
89  bool fcheckAngle; // boolean switch for cut on angle
90  bool fcheckDirection; // boolean switch for cut on direction
91 
92  // Histograms
93  TH2F *HitsXZ = new TH2F("hitsXZ", ";z (plane); x (vertical cell)", 900, 0., 900., 400, 0., 400.); // hit locations in XZ view
94  TH2F *HitsYZ = new TH2F("hitsYZ", ";z (plane); y (horizontal cell)", 900, 0., 900., 400, 0., 400.); // hit locations in YZ view
95 
96  TH1F *fslope = new TH1F("slope",";Slope", 220, -55, 55); // slope of the fit for trigger
97  TH1F *fexpected = new TH1F("expected",";expTime (ns)", 80, -40, 40); // expected hit time difference
98  TH1F *fmeasured = new TH1F("measured", ";meaTime (ns)", 2000, -1000, 1000); // measured hit time difference
99 
100  TH2F *fdirslp = new TH2F("dirslp", ";angle (degrees); slope", 220, -55, 55, 100, 0, 180); // graph of slope vs. angle cos(Y) in the x-z plane
101  TH2F *fexpmea = new TH2F("expmea", ";expected time (ns); measured time (ns)", 200, 0, 120, 400, -500, 500); // graph of expected time vs. measured time
102 
103  /*
104  TH1D *ftrue; // true hit time
105  TH1D *fsmear; // difference between smeared time and true time
106  TH2D *fupmuper; // percent of upmus recognized vs. cut hits in the x-z plane
107  */
108 };
109 //////////////////////////////////////////////////////////////////////
111 : _prescale (p.get<unsigned int>("prescale" ))
112 , _TrackModuleLabel (p.get<std::string >("TrackModuleLabel" ))
113 , _TrackInstanceLabel (p.get<std::string >("TrackInstanceLabel", ""))
114 , _SliceModuleLabel (p.get<std::string >("SliceModuleLabel" ))
115 , _SliceInstanceLabel (p.get<std::string >("SliceInstanceLabel", ""))
116 , fMinPlanes (p.get<unsigned int>("MinPlanes" ))
117 , fMinCellsX (p.get<unsigned int>("MinCellsX" ))
118 , fMinCellsY (p.get<unsigned int>("MinCellsY" ))
119 , fMinX (p.get<float> ("Xmin" ))
120 , fMaxX (p.get<float> ("Xmax" ))
121 , fMinY (p.get<float> ("Ymin" ))
122 , fMaxY (p.get<float> ("Ymax" ))
123 , fMinZ (p.get<unsigned int>("Zmin" ))
124 , fMaxZ (p.get<unsigned int>("Zmax" ))
125 , fAngleCut (p.get<double> ("AngleCut" ))
126 , fDirectionCut (p.get<double> ("DirectionCut" ))
127 , fcheckMinHits (p.get<bool > ("CheckMinHits" ))
128 , fcheckContainment (p.get<bool > ("CheckContainment" ))
129 , fcheckAngle (p.get<bool > ("CheckAngle" ))
130 , fcheckDirection (p.get<bool > ("CheckDirection" ))
131 {
132  std::cout << "=== novaddt::UpMuTestTrigger begin ===" << std::endl;
133  std::cout << "\t prescale: " << _prescale << std::endl;
134  std::cout << "\t TrackModuleLabel: " << _TrackModuleLabel << std::endl;
135  std::cout << "\t TrackInstanceLabel: " << _TrackInstanceLabel << std::endl;
136  std::cout << "\t SliceModuleLabel: " << _SliceModuleLabel << std::endl;
137  std::cout << "\t SliceInstanceLabel: " << _SliceInstanceLabel << std::endl;
138  std::cout << "\t MinPlanes: " << fMinPlanes << std::endl;
139  std::cout << "\t MinCellsX: " << fMinCellsX << std::endl;
140  std::cout << "\t MinCellsY: " << fMinCellsY << std::endl;
141  std::cout << "\t X fiducial cuts: " << fMinX << " to " << fMaxX << std::endl;
142  std::cout << "\t Y fiducial cuts: " << fMinY << " to " << fMaxY << std::endl;
143  std::cout << "\t Z fiducial cuts: " << fMinZ << " to " << fMaxZ << std::endl;
144  std::cout << "\t Angle cut: " << fAngleCut << std::endl;
145  std::cout << "\t Direction cut: " << fDirectionCut << std::endl << std::endl;
146  std::cout << "\t CheckMinHits: " << fcheckMinHits << std::endl;
147  std::cout << "\t CheckContainment: " << fcheckContainment << std::endl;
148  std::cout << "\t CheckAngle: " << fcheckAngle << std::endl;
149  std::cout << "\t CheckDirection: " << fcheckDirection << std::endl;
150  produces< std::vector<novaddt::TriggerDecision>>();
151  produces< art::Assns<novaddt::TriggerDecision, novaddt::HitList> >();
152 }
153 //////////////////////////////////////////////////////////////////////
155 {
156  // Clean up dynamic memory and other resources
157 }
158 //////////////////////////////////////////////////////////////////////
160 {
161  //////////////////////////////
162  // Events //
163  //////////////////////////////
164 
165  // UpMuTestTrigger event status
166  LOG_DEBUG("UpMuTestTrigger") << "=== novaddt::UpMuTestTrigger filter === :Event: "
167  << event.id().event()
168  << std::endl;
169  _nEvents++;
170 
171  if(_nEvents % 25 == 0){
172  std::cout << "=== novaddt::UpMuTestTrigger === EVENT: " << _nEvents << std::endl;
173  }
174 
175  // Trigger decisions and assignment vectors written out to the event
176  std::unique_ptr<std::vector<novaddt::TriggerDecision>> trigger_decisions(new std::vector<novaddt::TriggerDecision>);
177  std::unique_ptr<art::Assns<novaddt::TriggerDecision, novaddt::HitList> > assn(new art::Assns<novaddt::TriggerDecision, novaddt::HitList >);
178 
179 
180 
181  //////////////////////////////
182  // Tracks //
183  //////////////////////////////
184 
185  // Obtaining tracks
187  event.getByLabel(_TrackModuleLabel, _TrackInstanceLabel, tracks);
188  LOG_DEBUG("UpMuTestTrigger") << "\tgot " << tracks->size() << " tracks" << std::endl;
189  // Obtaining hitlist for each track
191 
192  // Examining tracks
193  for(size_t i_track = 0; i_track < tracks->size(); ++++i_track){
194  _nTracks++;
195  novaddt::Track x_track = tracks->at(i_track);
196  novaddt::Track y_track = tracks->at(i_track+1);
197  // Obtain the hit list for this track
198  art::Ptr<novaddt::HitList> x_hit_list = fohl.at(i_track);
199  art::Ptr<novaddt::HitList> y_hit_list = fohl.at(i_track+1);
200 
201 
202 
203  //////////////////////////////
204  // Hits //
205  //////////////////////////////
206 
207  /*
208  // DEBUG:: Hits loop
209  // XZ hits
210  for(size_t i_hit = 0; i_hit < x_hit_list->size(); ++i_hit){
211  LOG_DEBUG("NuMuTestTrigger") << "\t\t\thit[" << i_hit
212  << "]: \tTDC: " << x_hit_list->at(i_hit).TDC().val
213  << ", \tADC: " << x_hit_list->at(i_hit).ADC().val
214  << ", \tplane: " << x_hit_list->at(i_hit).Plane().val
215  << ", \tx-cell: " << x_hit_list->at(i_hit).Cell().val
216  << std::endl;
217  } // End of loop on x-hits
218 
219  // YZ hits
220  for(size_t i_hit = 0; i_hit < y_hit_list->size(); ++i_hit){
221  LOG_DEBUG("NuMuTestTrigger") << "\t\t\thit[" << i_hit
222  << "]: \tTDC: " << y_hit_list->at(i_hit).TDC().val
223  << ", \tADC: " << y_hit_list->at(i_hit).ADC().val
224  << ", \tplane: " << y_hit_list->at(i_hit).Plane().val
225  << ", \ty-cell: " << y_hit_list->at(i_hit).Cell().val
226  << std::endl;
227  } // End of loop on Y-hits
228  */
229 
230 
231 
232  //////////////////////////////
233  // Trigger //
234  //////////////////////////////
235 
236  // Checking trigger conditions
237 
238  //*************************************************//
239  //******** Hits cut (minimal track length) ********//
240  //*************************************************//
241 
242  double x_dZ = std::fabs(x_track.EndZ() - x_track.StartZ());
243  double y_dZ = std::fabs(y_track.EndZ() - y_track.StartZ());
244  double dZ = std::max(x_dZ,y_dZ);
245 
246  if(fcheckMinHits){
247  //CUT::dX or dZ (horizonal projection minimum hits)
248  if (x_hit_list->size() < fMinCellsX || dZ < fMinPlanes){
249  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to lack horizontal plane hits (dX and dZ cut)" << std::endl;
250  continue;
251  }
253 
254  //CUT::dY (vertical projection minimum hits)
255  if (y_hit_list->size() < fMinCellsY){
256  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to lack of verticle plane hits (dY cut)" << std::endl;
257  continue;
258  }
260  }
261 
262  //***************************************************************************//
263  //******** Containment cut (through-going, stopped, fully-contained) ********//
264  //***************************************************************************//
265 
266  double min_Z = std::min(x_track.StartZ(),y_track.StartZ());
267  double max_Z = std::max(x_track.EndZ(), y_track.EndZ());
268 
269  if(fcheckContainment){
270  // Through-going muon veto (X and Z horizontal vetos)
271  //CUT::X containment
272  if ((x_track.StartV() < fMinX) || (x_track.StartV() > fMaxX) ||
273  (x_track.EndV() < fMinX) || (x_track.EndV() > fMaxX)){
274  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to x-containment" << std::endl;
275  continue;
276  }
277  //CUT::Z containment
278  if ((min_Z < fMinZ) || (max_Z> fMaxZ) ||
279  (max_Z < fMinZ) || (min_Z> fMaxZ)){
280  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to z-containment" << std::endl;
281  continue;
282  }
283  _nPass_through++;
284 
285  // Stopped muon veto (Top Y vertical veto)
286  //CUT::Y maximum containment
287  if ((y_track.StartV() > fMaxY) || (y_track.EndV() > fMaxY)){
288  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to maximum y-containment" << std::endl;
289  continue;
290  }
291  _nPass_stopped++;
292 
293  // Fully-contained muon veto (Bottom Y vertical veto)
294  //CUT::Y minimum containment
295  if ((y_track.StartV() < fMinY) || (y_track.EndV() < fMinY)){
296  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to minimum y-containment" << std::endl;
297  continue;
298  }
299  _nPass_fully++;
300  }
301 
302  //*******************************************//
303  //******** Angle cut (vertical cone) ********//
304  //*******************************************//
305 
306  double scaled_dZ = (dZ * 6.6) / 3.9; // Planes longer (z) than Cells (x or y)
307  double dX = std::fabs( x_track.EndV() - x_track.StartV() );
308  double dY = std::fabs( y_track.EndV() - y_track.StartV() );
309  double CosTrack = ( scaled_dZ/ std::sqrt(dX*dX+dY*dY+scaled_dZ*scaled_dZ) );
310 
311  //CUT::Upward facing 30 degree cone is default
312  if (fcheckAngle){
313  if (CosTrack > fAngleCut){
314  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to angle cut" << std::endl;
315  continue;
316  }
317  _nPass_angle++;
318  }
319 
320  /*
321  //// EXTRA::Alternative angular cut-- a simple loose cut on the inverse slope |dz|/|dy|
322  //// inspired by doc-db 8888
323  //// at the moment I'm not sure if the x, y & z coordinates provided are sensible
324  //// and in the same aspect ratio.
325  //// I doubt they are. This should be re-visited
326  //double dY = std::fabs( y_track.EndV() - y_track.StartV() );
327  ////double dZ = std::fabs( tracks->at(i_track).EndZ() - tracks->at(i_track).StartZ() );
328  //double inverse_slope = dZ / dY;
329  //std::cout << "\t\t\tdY: " << dY
330  //<< ", dZ: " << dZ
331  //<< ", inverse slope: " << inverse_slope
332  //<< std::endl;
333  //if (inverse_slope > _inverseSlopeCut){
334  //std::cout << "\t\t - failed the trigger decision due to inverse slope cut" << std::endl;
335  //continue;
336  //}
337  */
338 
339  //**************************************************//
340  //******** Direction cut (upward traveling) ********//
341  //**************************************************//
342 
343  if(fcheckDirection){
344  /*
345  // Uncertainty analysis
346  double min_tdc = std::min(x_hit_list->front().TDC().val,y_hit_list->front().TDC().val);
347  double max_tdc = std::max(x_hit_list->back().TDC().val, y_hit_list->back().TDC().val);
348  double dTime = max_tdc - min_tdc; // Time separation between track beginning and end
349  double dLength = std::sqrt(dX*dX+dY*dY+scaled_dZ*scaled_dZ); // Track length
350  */
351 
352  //################# ALGORITHM 1 #################//
353  // Slope correlation method ///
354  //###############################################//
355  /*
356  // Use T for Y cells
357 
358  // Sort hits by Y cell position
359  unsigned int nHits = y_hit_list->size();
360  //art::Handle<novaddt::HitList> hits;
361  HitList hit_list(*y_hit_list);
362  std::sort(hit_list.begin(), hit_list.end(),CompareDAQHit<Cell>());
363 
364  // Definitions
365  int loop = 0;
366  unsigned int nPairs = (nHits+1) * nHits / 2;
367  double tmte=0, tm=0, te=0, te2=0, slp, angle;
368  double mea[nPairs], exp[nPairs];
369 
370  for(size_t i_hit = 0; i_hit < nHits; i_hit++){
371  for (size_t j_hit = 0; j_hit < i_hit; j_hit++){
372  // Variables
373  //Ti = (hit_list.at(i_hit).TDC().val * 15.625) - ((-hit_list.at(i_hit).Cell().val + 192) * dC / 1.59); // Photon transport (fiber optical) corrected time
374  //Tj = (hit_list.at(j_hit).TDC().val * 15.625) - ((-hit_list.at(j_hit).Cell().val + 192) * dC / 1.59); // n = 1.59 index; (dx/2) = 192 cells;
375 
376  // Measured time
377  mea[loop] = hit_list.at(i_hit).TDC().val*15.625 - hit_list.at(j_hit).TDC().val*15.625; // TDC->ns conversion (64 MHz clock ticks)
378  fmeasured->Fill(mea[loop]);
379 
380  // Expected time
381  exp[loop] = hit_list.at(i_hit).Cell().val*0.13008999712728 - hit_list.at(j_hit).Cell().val*0.13008999712728; // Cell->ns converstion (3.9*1e7/299792458)
382  fexpected->Fill(exp[loop]);
383 
384  // Analysis
385  fexpmea->Fill(exp[loop], mea[loop]);
386 
387  // Counters
388  tmte += mea[loop] * exp[loop];
389  tm += mea[loop];
390  te += exp[loop];
391  te2 += exp[loop] * exp[loop];
392  loop++;
393 
394  // Check timing (j)
395  //std::cout << "Expected: " << exp[loop] << ", Measured: " << mea[loop] << std::endl;
396  } // End of x comp hits loops (j_hit)
397 
398  // Check timing (i)
399  //std::cout << Ti << std::endl;
400  } // End of x test hits loop (i_hit)
401 
402  // Metrics
403  slp = (loop*tmte-tm*te)/(loop*te2-te*te);
404  //angle = std::acos(dY/scaled_dZ) * 180 / 3.14159;
405  angle = std::acos(CosTrack) * 180 / 3.14159;
406  //std::cout << "Track time sigma angle: " << angle << ", Track time sigma slope: " << slp << std::endl;
407 
408  fslope->Fill(slp);
409  fdirslp->Fill(slp, angle);
410  */
411  //################# ALGORITHM 2 #################//
412  /// Adjacent cell analysis ///
413  //###############################################//
414 
415  // Get Y position for X cells-- use X cell's T for these for better resolution (add in Y cell's T too?) (use mergetracks/Track3D to get vertex?)
416 
417  unsigned int yHits = y_hit_list->size();
418  HitList check_hit_list(*y_hit_list); // Y hitlist
419  unsigned int xHits = x_hit_list->size();
420  HitList hit_list(*x_hit_list); // X hitlist
421 
422  // Sort X hits by Y cell position
423  std::sort(check_hit_list.begin(), check_hit_list.end(), CompareDAQHit<Plane>());
424  int sortY = check_hit_list.at(yHits-1).Cell().val - check_hit_list.at(0).Cell().val; // pos for +y, nega for -y:: .back() - .front()
425  if(sortY < 0){
426  std::sort(hit_list.rbegin(), hit_list.rend(), CompareDAQHit<Plane>()); // descending (reverse) order
427  //std::sort(hit_list.begin(), hit_list.end(), CompareDAQHit<Plane>());
428  //std::reverse(hit_list.begin(),hit_list.end());
429  } else{
430  std::sort(hit_list.begin(), hit_list.end(), CompareDAQHit<Plane>()); // ascending (normal) order
431  }
432 
433  // // Get Y position for X cells (relative to first cell)
434  // //std::sort(check_hit_list.begin(), check_hit_list.end(), CompareDAQHit<Cell>());
435  // double dYc, dYp, sum=0, m;
436  // for(size_t y_hit = 0; y_hit < (yHits - 1); y_hit++){
437  // dYc= check_hit_list.at(y_hit+1).Cell().val - check_hit_list.at(y_hit).Cell().val;
438  // dYp= check_hit_list.at(y_hit+1).Plane().val - check_hit_list.at(y_hit).Plane().val;
439  // sum += dYc/dYp;
440  // std::cout << sum << std::endl;
441  //
442  // }
443 
444  // Get Y position for X cells (relative to first cell) [need to improve slope determination with more points]
445  double dYc, dYp, m;
446  dYc= check_hit_list.at(yHits-1).Cell().val - check_hit_list.at(0).Cell().val;
447  dYp= check_hit_list.at(yHits-1).Plane().val - check_hit_list.at(0).Plane().val;
448  m = dYc / dYp; // mean slope of Y cell positions as a function of plane (take average of y and see if negative or positive-- can be used for sorting instead of current method)
449 
450  // Fill Y position for sorted X cells
451  double AltCell[xHits]; // Y position of X cells relative to first X cell
452  for(size_t x_hit = 0; x_hit < xHits; x_hit++){
453  AltCell[x_hit] = m * hit_list.at(x_hit).Plane().val + m * hit_list.at(0).Plane().val;
454  }
455 
456  // // Sort hits by Y cell position
457  // if(m > 0){
458  // std::sort(hit_list.begin(), hit_list.end(), CompareDAQHit<Cell>()); // ascending order (Plane?)
459  // } else{
460  // std::sort(hit_list.begin(), hit_list.end(), CompareDAQHit<Cell>());
461  // std::reverse(hit_list.begin(),hit_list.end());
462  // //std::sort(hit_list.rbegin(), hit_list.rend(), CompareDAQHit<Cell>()); // descending (reverse) order
463  // }
464 
465  //
466  // // Sort hits by Y cell position
467  // for(size_t i_hit = 0; i_hit < nHits; i_hit++){
468  // for (size_t j_hit = 0; j_hit < i_hit; j_hit++){
469  // if(AltCell[j_hit] < AltCell[i_hit]){
470  // DAQHit tmp = hit_list.at(i_hit);
471  // double tmpalt = AltCell[i_hit];
472  //
473  // hit_list.at(i_hit) = hit_list.at(j_hit);
474  // AltCell[i_hit] = AltCell[j_hit];
475  //
476  // hit_list.at(j_hit) = tmp;
477  // AltCell[j_hit] = tmpalt;
478  // }
479  // }
480  // }
481 
482  // Definitions
483  unsigned int nPairs = (xHits+1) * xHits / 2;
484  double tmte=0, tm=0, te=0, te2=0, slp, angle, dT, dC;
485  double mea[nPairs], exp[nPairs];
486  int loop = 0;
487 
488  for(size_t i_hit = 0; i_hit < xHits; i_hit++){
489  for (size_t j_hit = 0; j_hit < i_hit; j_hit++){
490  // Variables
491  //Ti = (hit_list.at(i_hit).TDC().val * 15.625) - ((-hit_list.at(i_hit).Cell().val + 192) * dC / 1.59); // Photon transport (fiber optical) corrected time
492  //Tj = (hit_list.at(j_hit).TDC().val * 15.625) - ((-hit_list.at(j_hit).Cell().val + 192) * dC / 1.59); // n = 1.59 index; (dx/2) = 192 cells;
493  dT = hit_list.at(i_hit).TDC().val * 15.625 - hit_list.at(j_hit).TDC().val * 15.625; // TDC->ns conversion (64 MHz clock ticks)
494  //dC = hit_list.at(i_hit).Cell().val * 0.13008999712728 - hit_list.at(j_hit).Cell().val * 0.13008999712728; // Cell->ns converstion (speed of light and cell scaling-- 3.9*1e7/299792458) // X position
495  dC = AltCell[i_hit] * 0.13008999712728 - AltCell[j_hit] * 0.13008999712728; // Y position
496 
497  // Measured time
498  mea[loop] = dT;
499  fmeasured->Fill(mea[loop]);
500 
501  // Expected time
502  exp[loop] = dC;
503  fexpected->Fill(exp[loop]);
504 
505  // Analysis
506  fexpmea->Fill(exp[loop], mea[loop]);
507 
508  // Counters
509  tmte += mea[loop] * exp[loop];
510  tm += mea[loop];
511  te += exp[loop];
512  te2 += exp[loop] * exp[loop];
513  loop++;
514 
515  // Check timing (j)
516  //std::cout << "Expected: " << exp[loop] << ", Measured: " << mea[loop] << std::endl;
517  } // End of x comp hits loops (j_hit)
518 
519  // Check timing (i)
520  //std::cout << Ti << std::endl;
521  } // End of x test hits loop (i_hit)
522 
523  // Metrics
524  slp = (loop*tmte-tm*te)/(loop*te2-te*te);
525  //angle = std::acos(dY/scaled_dZ) * 180 / 3.14159;
526  angle = std::acos(CosTrack) * 180 / 3.14159;
527  //std::cout << "Track time sigma angle: " << angle << ", Track time sigma slope: " << slp << std::endl;
528 
529  fslope->Fill(slp);
530  fdirslp->Fill(slp, angle);
531 
532  //---------------------------------------------------------------------------
533  // CUT::Slope must be positive to be upward going (assuming fDirectionCut==0)
534  if(slp < fDirectionCut){
535  LOG_DEBUG("UpMuTestTrigger") << "\t\t - failed the trigger decision due to direction cut" << std::endl;
536  continue;
537  }
539  }
540 
541 
542 
543  //////////////////////////////
544  // UpMu //
545  //////////////////////////////
546 
547  // Track has passed trigger decision and is now identified as upward traveling muon
548  _trigger_counts++;
549  // check the prescale
552  LOG_DEBUG("NuMuTrigger") << "\t\t - not prescaled" << std::endl;
553  //DEBUG::Analyze UpMu tracks
554  // XZ projection
555  for(size_t i_hit = 0; i_hit < x_hit_list->size(); ++i_hit){
556  HitsXZ -> Fill(x_hit_list->at(i_hit).Plane().val, x_hit_list->at(i_hit).Cell().val);
557  } // End of loop on x hits
558  // YZ projection
559  for(size_t i_hit = 0; i_hit < y_hit_list->size(); ++i_hit){
560  HitsYZ -> Fill(y_hit_list->at(i_hit).Plane().val, y_hit_list->at(i_hit).Cell().val);
561  } // End of loop on x hits
562 
563  uint64_t min_tdc_trig = std::min(x_hit_list->front().TDC().val,y_hit_list->front().TDC().val);
564  uint64_t max_tdc_trig = std::max(x_hit_list->back().TDC().val, y_hit_list->back().TDC().val);
565  trigger_decisions->push_back(novaddt::TriggerDecision(
566  min_tdc_trig,
567  max_tdc_trig - min_tdc_trig,
569  _prescale)
570  );
571  util::CreateAssn(*this, event, *trigger_decisions, x_hit_list, *assn);
572  util::CreateAssn(*this, event, *trigger_decisions, y_hit_list, *assn);
573  LOG_DEBUG("UpMuTestTrigger") << "\t\t - passed the trigger decision, start time: "
574  << trigger_decisions->back().start()
575  << ", duration: " << trigger_decisions->back().duration()
576  << std::endl;
577  } // End of check on prescale
578  } // End of loop on tracks
579 
580  bool goodTrigger = (trigger_decisions->size() > 0);
581  LOG_DEBUG("UpMuTestTrigger") << "=== end of novaddt::UpMuTestTrigger filter === :Trigger decision: " << goodTrigger
582  << " from " << trigger_decisions->size()
583  << " passes"<< std::endl;
584 
585  event.put(std::move(trigger_decisions));
586  event.put(std::move(assn));
587 
588  return goodTrigger;
589 } // End of loop on events
590 //////////////////////////////////////////////////////////////////////
592 {
593  std::cout << "=== novaddt::UpMuTestTrigger end === " << std::endl;
594  std::cout << "\tNumber of events: " << _nEvents << std::endl;
595  std::cout << "\tNumber of tracks: " << _nTracks << " / " << "Trigger Rejection" << std::endl;
596  std::cout << "\tNumber that pass min horizontal hits: " << _nPass_dhorizontal << " / " << 100 - (100*(float)(_nPass_dhorizontal)/(float)(_nTracks)) << "%" << std::endl;
597  std::cout << "\tNumber that pass min vertical hits: " << _nPass_dvertical << " / " << 100 - (100*(float)(_nPass_dvertical)/(float)(_nTracks)) << "%" << std::endl;
598  std::cout << "\tNumber that pass through containment (X/Z): " << _nPass_through << " / " << 100 - (100*(float)(_nPass_through)/(float)(_nTracks)) << "%" << std::endl;
599  std::cout << "\tNumber that pass stopped containment (MaxY): " << _nPass_stopped << " / " << 100 - (100*(float)(_nPass_stopped)/(float)(_nTracks)) << "%" << std::endl;
600  std::cout << "\tNumber that pass fully containment (MinY): " << _nPass_fully << " / " << 100 - (100*(float)(_nPass_fully)/(float)(_nTracks)) << "%" << std::endl;
601  std::cout << "\tNumber that pass angle: " << _nPass_angle << " / " << 100 - (100*(float)(_nPass_angle)/(float)(_nTracks)) << "%" << std::endl;
602  std::cout << "\tNumber that pass direction: " << _nPass_direction << " / " << 100 - (100*(float)(_nPass_direction)/(float)(_nTracks)) << "%" << std::endl;
603  std::cout << "\tNumber that pass the trigger: " << _trigger_counts << std::endl;
604  std::cout << "\t after prescale: " << _prescaled_trigger_counts << std::endl;
605  // Write analysis histograms
606  TFile f("upmu_hists.root","recreate");
607 
608  HitsXZ->Write();
609  HitsYZ->Write();
610 
611  fslope->Write();
612  fexpected->Write();
613  fmeasured->Write();
614  fdirslp->Write();
615  fexpmea->Write();
616  /*
617  ftrue->Write();
618  fsmear->Write();
619  fupmuper->Write();
620  */
621 }
622 
623 //////////////////////////////////////////////////////////////////////
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
Double_t angle
Definition: plot.C:86
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
float const & StartZ() const
Definition: Track.h:68
std::string _SliceInstanceLabel
instance label making the HitList
std::vector< DAQHit > HitList
Definition: HitList.h:15
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool filter(art::Event &e) override
unsigned int _prescale
prescale factor (1 out of every this many passes are issued as real triggers)
DEFINE_ART_MODULE(TestTMapFile)
T acos(T number)
Definition: d0nt_math.hpp:54
float const & EndV() const
Definition: Track.h:70
std::string _TrackModuleLabel
label of module making the Tracks
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
UpMuTestTrigger(fhicl::ParameterSet const &p)
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
float const & StartV() const
Definition: Track.h:67
float const & EndZ() const
Definition: Track.h:71
T min(const caf::Proxy< T > &a, T b)
Float_t e
Definition: plot.C:35
Definition: fwd.h:28
std::string _SliceModuleLabel
label of module making the HitList
std::string _TrackInstanceLabel
instance label making the Tracks