Leana_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////
2 /// \brief Module to analyze low energy events and noise
3 /// \author messier@indiana.edu
4 /// \date May 2011
5 ////////////////////////////////////////////////////////////////////
6 #include <vector>
7 
8 // ROOT includes
9 #include "TH1D.h"
10 #include "TH2D.h"
11 #include "TTree.h"
12 
13 // ART includes
18 
19 // NOvA includes
22 #include "Commissioning/LeanaNt.h"
23 #include "Geometry/Geometry.h"
24 #include "GeometryObjects/Geo.h"
25 #include "RawData/RawDigit.h"
26 #include "RawData/RawTrigger.h"
27 #include "TrackFit/RLFit.h"
28 #include "Utilities/func/MathUtil.h"
29 #include "Utilities/RawUtil.h"
30 
31 
32 namespace comi {
33  class LeanaNt;
34 
35  /// Summary information about previous muons
36  struct MuonCand {
37  unsigned char fId; ///< Stop? Through-going?
38  float fX1; ///< Start x position (cm)
39  float fY1; ///< Start y position (cm)
40  float fZ1; ///< Start z position (cm)
41  float fX2; ///< End x position (cm)
42  float fY2; ///< End y position (cm)
43  float fZ2; ///< End z position (cm)
44  double fT; ///< Time in spill (usec)
45  double fNOvAT; ///< "Global" NOvA time (usec)
46  };
47 
48  /// Analysis module for low energy events and noise
49  class Leana : public art::EDAnalyzer {
50  public:
51  explicit Leana(fhicl::ParameterSet const& pset);
52  virtual ~Leana();
53 
54  void beginJob();
55  void analyze(art::Event const& evt);
56  void reconfigure(const fhicl::ParameterSet& p);
57 
58  private:
59  void Ana1Hit(unsigned int iwin);
60  void AnaNHit(unsigned int iwin);
61  void MuonCandidate(unsigned int iwin,
62  unsigned int sz,
63  unsigned int nx,
64  unsigned int ny);
65  void FindClosestMuon(unsigned int i,
66  double x,
67  double y,
68  double z,
69  double t,
70  double* dstop,
71  double* tstop,
72  unsigned int* istop,
73  double* dperp,
74  double* tperp,
75  unsigned int* iperp,
76  bool assn=true);
77  void DumpMuon(unsigned int i);
78  void DumpMuonTable();
79  void DumpMichel(unsigned int i);
80  bool IsPointOnEdge(double x, double y, double z);
81 
82  private:
83  TH1D* fEventCount; ///< Event counter
84  TH1D* fSliceSz; ///< Histogram of Nhit per slice
85  TH2D* fSliceSzVsEvent; ///< Histogram of slice size vs. event
86 
87  TH1D* f1HitADC; ///< ADC spectrum for 1 hit slices
88  TH1D* f1HitTDC; ///< TDC spectrum for 1 hit slices
89  TH2D* f1HitADCvsTDC; ///< ADC vs TDC for 1 hit slices
90 
91  TH1D* f2HitDplane; ///< Distance from first to last plane
92  TH1D* f3HitDplane; ///< Distance from first to last plane
93 
94  TH1D* f3HitSamePlaneDcell; ///< Max-min cell in same plane
95  TH1D* f3HitDiffPlaneDcell; ///< Max-min cell in different planes
96 
97  LeanaNt* fNt; ///< Ntuple for hit clusters
98  TTree* fNtTree; ///< Ntuple for hit clusters
99 
100  private:
101  //
102  // Buffer size should be the calibration trigger rate * the size
103  // of fDeltaTnamu. For 20 seconds at 10 Hz, this suggests a
104  // minimum 2000
105  //
106  static const unsigned int MUON_BUFF_SZ = 4096;
107  static const char kNotMu = 0; ///< Not a good muon candidate
108  static const char kStopMu = 1; ///< A stopping muon candidate
109  static const char kThruMu = 2; ///< A through-going muon candidate
110 
111  private:
112  std::vector< art::Ptr<rawdata::RawDigit> > fDigi;
113 
114  std::vector<util::RawSlice> fWin; ///< The raw slices for this event
115  double fT0; ///< Event NOvA trigger time (usec)
116 
117  bool fDoFill; ///< Should we be filling ntuple / histogram?
118  unsigned int fMuIndx1; ///< Current muon list
119  unsigned int fMuIndx2; ///< Old non-associated muons
120  std::vector<MuonCand> fMuon[MUON_BUFF_SZ]; ///< Table of muon data
121 
122  private:
123  int fDumpMichel; ///< Dump table of Michel electrons?
124 
125  private:
126  double fTimeWindow; ///< [usec] Time window size
127  unsigned int fNhit; ///< Number of hits req'd in window
128  unsigned int fNhitX; ///< Number of hits req'd in x view
129  unsigned int fNhitY; ///< Number of hits req'd in y view
130 
131  int fDcellCut; ///< Allowed gap between hit cells
132  int fDplaneCut; ///< Allowed gap between hit planes
133 
134  unsigned int fMuonNhit; ///< Required # hits on muon track
135  unsigned int fMuonNhitX; ///< # required in x view
136  unsigned int fMuonNhitY; ///< # required in y view
137  double fEdgeBuffer; ///< [cm] Edge distance for stopping/entering decision
138  double fDeltaTnamu; ///< Time for non-association (s)
139  };
140 
141  // NOvA uses a 64 Mhz clock
142  static const double kUSEC_PER_TDC = (1.0E6/64.0E6);
143  static const int kTDC_PER_USEC = 64;
144 
145  //......................................................................
146 
148  EDAnalyzer(p),
149  fDoFill(true),
150  fMuIndx1(0),
151  fMuIndx2(0)
152  {
153  this->reconfigure(p);
154  }
155 
156  //......................................................................
157 
159  {
160  fTimeWindow = p.get<double>("TimeWindow");
161  fNhit = p.get<unsigned int>("Nhit");
162  fNhitX = p.get<unsigned int>("NhitX");
163  fNhitY = p.get<unsigned int>("NhitY");
164  fDcellCut = p.get<unsigned int>("DcellCut");
165  fDplaneCut = p.get<unsigned int>("DplaneCut");
166  fMuonNhit = p.get<unsigned int>("MuonNhit");
167  fMuonNhitX = p.get<unsigned int>("MuonNhitX");
168  fMuonNhitY = p.get<unsigned int>("MuonNhitY");
169  fEdgeBuffer = p.get<double>("EdgeBuffer");
170  fDeltaTnamu = p.get<double>("DeltaTnamu");
171  fDumpMichel = p.get<int>("DumpMichel");
172  }
173 
174  //......................................................................
175 
177  {
178  //
179  // Some macros to make booking histograms a little easier
180  //
181 #define H1(NAME,LABEL,N,X1,X2) \
182  NAME = f->make<TH1D>(#NAME,LABEL,N,X1,X2)
183 
184 #define H2(NAME,LABEL,NX,X1,X2,NY,Y1,Y2) \
185  NAME = f->make<TH2D>(#NAME,LABEL,NX,X1,X2,NY,Y1,Y2)
186 
188  //
189  // General event-level histograms
190  //
191  H1(fEventCount, ";;Events Seen", 1, -0.5,0.5);
192  H1(fSliceSz, ";N_{hit};N_{slices}", 101, -0.5,100.5);
193  H2(fSliceSzVsEvent,";Event;N_{hit}", 100, 0.0,1000000,101,-0.5,100.5);
194 
195  //
196  // 1-hit slice histograms
197  //
198  H1(f1HitADC,";ADC;N_{hit}",4096,-0.5,4095.5);
199  H1(f1HitTDC,";t [#musec]; N_{hit}",700,-100,600);
200  H2(f1HitADCvsTDC,";t [#musec]; ADC",
201  700,-100,600,4096,-0.5,4095.5);
202 
203  H1(f2HitDplane,";#Delta_{plane}",51,-0.5,50.5);
204  H1(f3HitDplane,";#Delta_{plane}",51,-0.5,50.5);
205 
206  H1(f3HitSamePlaneDcell,";#Delta_{cell}",51,-0.5,50.5);
207  H1(f3HitDiffPlaneDcell,";#Delta_{cell}",51,-0.5,50.5);
208 
209  fNt = new LeanaNt;
210  fNtTree = f->make<TTree>("le","Leana Analysis Tree");
211  fNtTree->Branch("fLeanaNt","comi::LeanaNt",&fNt,1600,99);
212  }
213 
214  //......................................................................
215 
217 
218  //......................................................................
219 
221  {
222  unsigned int i;
223  unsigned int fMaxSz = 10;
224 
225  fNt->run = evt.run();
226  fNt->subrun = evt.subRun();
227  fNt->evt = evt.id().event();
228 
229  //
230  // Get trigger and t0 data from the raw trigger information
231  //
233  evt.getByLabel("daq", trigv);
234  const rawdata::RawTrigger& trig = (*trigv)[0];
235  fT0 = trig.TDCT0()*kUSEC_PER_TDC;
236 
237  //
238  // Pull out the raw digits from the event
239  //
241  evt.getByLabel("daq", digidummy);
242 
243  //
244  // Shuffle the data into a vector we can manipulate
245  //
246  fDigi.clear();
247  for (i=0; i<digidummy->size(); ++i) {
248  fDigi.push_back(art::Ptr<rawdata::RawDigit>(digidummy, i));
249  }
251  badc->Apply(fDigi);
252  if (fDigi.size()==0) return;
253 
254  if (fDoFill) fEventCount->Fill(0.0);
255 
256  //
257  // Time sort and slice hits into windows
258  //
259  fWin.clear();
261  unsigned int tdcwin = rint(fTimeWindow*kTDC_PER_USEC);
263  if (fWin.size()<1) return;
264 
265  //
266  // Clear the muon table
267  //
268  fMuon[fMuIndx1].resize(fWin.size());
269  for (i=0; i<fWin.size(); ++i) {
270  fMuon[fMuIndx1][i].fId = kNotMu;
271  fMuon[fMuIndx1][i].fNOvAT = fT0;
272  fMuon[fMuIndx1][i].fT = 0.0;
273  }
274 
275  //
276  // Analyze the activity in each window in turn
277  //
278  for (i=0; i<fWin.size(); ++i) {
279  unsigned int sz = fWin[i].second-fWin[i].first+1;
280  unsigned int nx = 0;
281  unsigned int ny = 0;
282  util::CountXY(fDigi, fWin[i].first, fWin[i].second, &nx, &ny);
283 
284  fNt->iwin = i;
285  fNt->nx = nx;
286  fNt->ny = ny;
287 
288  if (fDoFill) fSliceSz->Fill(sz);
289  if (fDoFill) fSliceSzVsEvent->Fill(evt.id().event(), sz);
290 
291  // Analyze the small events
292  if (sz==1) this->Ana1Hit(i);
293  if (sz>1 && sz<=fMaxSz) this->AnaNHit(i);
294 
295  // Check if this activity is a muon
296  this->MuonCandidate(i, sz, nx, ny);
297  }
298 
299  // this->DumpMuonTable();
300 
301  //
302  // Compute the time between the muons at fMuIndx1 and fMuIndx2 in
303  // units of seconds
304  //
305  double tnassn = fabs(fMuon[fMuIndx2][0].fNOvAT -
306  fMuon[fMuIndx1][0].fNOvAT)/1.E6;
307  //
308  // Advance the muon indicies. The "current" entry (fMuIndx1) always
309  // advances. The secondary index (fMuIndx2) only advances if enough
310  // time have elapsed since the muon at the current index
311  //
313  if (tnassn>fDeltaTnamu) {
315  //
316  // Now that we have the data we need to compute backgrounds, we
317  // can start filling histograms.
318  //
319  fDoFill = true;
320  }
321  }
322 
323  //......................................................................
324 
325  void Leana::DumpMuon(unsigned int i)
326  {
327  std::cout << "[" << i << "] ";
328 
329  if (i>= fMuon[fMuIndx1].size()) {
330  std::cout << " not valid muon" << std::endl;
331  return;
332  }
333 
334  if (fMuon[fMuIndx1][i].fId == kNotMu) std::cout << " not-mu ";
335  if (fMuon[fMuIndx1][i].fId == kStopMu) std::cout << " stop-mu ";
336  if (fMuon[fMuIndx1][i].fId == kThruMu) std::cout << " thru-mu ";
337 
338  std::cout << " (x1,y1,z1)=("
339  << fMuon[fMuIndx1][i].fX1 << ","
340  << fMuon[fMuIndx1][i].fY1 << ","
341  << fMuon[fMuIndx1][i].fZ1 << ") ";
342 
343  std::cout << " (x2,y2,z2)=("
344  << fMuon[fMuIndx1][i].fX2 << ","
345  << fMuon[fMuIndx1][i].fY2 << ","
346  << fMuon[fMuIndx1][i].fZ2 << ") ";
347 
348  std::cout << " L=" <<
349  util::pythag((fMuon[fMuIndx1][i].fX1-fMuon[fMuIndx1][i].fX2),
350  (fMuon[fMuIndx1][i].fY1-fMuon[fMuIndx1][i].fY2),
351  (fMuon[fMuIndx1][i].fZ1-fMuon[fMuIndx1][i].fZ2));
352 
353  std::cout << " t=" << fMuon[fMuIndx1][i].fT << std::endl;
354  }
355 
356  //......................................................................
357 
359  {
360  std::cout << "*= indx1=" << fMuIndx1
361  << " has " << fMuon[fMuIndx1].size() << " entries\n";
362 
363  for (unsigned int i=0; i<fMuon[fMuIndx1].size(); ++i) this->DumpMuon(i);
364 
365  std::cout << " indx2=" << fMuIndx2
366  << " has " << fMuon[fMuIndx2].size() << " entries\n";
367  std::cout << " " << fMuon[fMuIndx2][0].fNOvAT << "\t"
368  << fMuon[fMuIndx2][0].fT << std::endl;
369 
370  }
371 
372  //......................................................................
373 
374  void Leana::DumpMichel(unsigned int imu)
375  {
376  std::cout << "||"
377  << "==================================="
378  << "==================================="
379  << std::endl;
380  std::cout << fNt->run << "/"
381  << fNt->subrun << "/"
382  << fNt->evt
383  << " Michel Candidate"
384  << std::endl;
385 
386  std::cout << "(nx,ny)=("
387  << fNt->nx << ","
388  << fNt->ny << ") "
389  << "(x,y,z,t)="
390  << fNt->x << ","
391  << fNt->y << ","
392  << fNt->z << ","
393  << fNt->t << ") "
394  << " (dstop,tstop)=("
395  << fNt->dstop << ","
396  << fNt->tstop << ")"
397  << std::endl;
398 
399  this->DumpMuon(imu);
400  }
401 
402  //......................................................................
403 
404  void Leana::Ana1Hit(unsigned int iwin)
405  {
406  unsigned int i = fWin[iwin].first;
407 
408  double t = fDigi[i]->TDC()*kUSEC_PER_TDC;
409  double q = fDigi[i]->ADC();
410 
411  if (fDoFill) f1HitADC->Fill(q);
412  if (fDoFill) f1HitTDC->Fill(t);
413  if (fDoFill) f1HitADCvsTDC->Fill(t,q);
414  }
415 
416  //......................................................................
417 
418  void Leana::AnaNHit(unsigned int iwin)
419  {
420  unsigned int i;
421 
424 
425  unsigned int i1 = fWin[iwin].first;
426  unsigned int i2 = fWin[iwin].second;
427  unsigned int j1, j2;
428  //
429  // Look for the smallest gap between each hit and its neighbors in
430  // terms of cell and plane numbers. Flag isolated hits as bad.
431  //
432  static const int kUNINIT = 999999;
433  std::vector<bool> goodhit(i2+1);
434  for (j1=i1; j1<=i2; ++j1) {
435  int dcell = kUNINIT; // distance to closest cell
436  int dplane = kUNINIT; // distance to closest plane
437 
438  int p1 = cmap->GetPlane(fDigi[j1].get());
439  int c1 = cmap->GetCell (fDigi[j1].get());
440  int v1 = geom->Plane(p1)->View();
441 
442  for (j2=i1; j2<=i2; ++j2) {
443  if (j1==j2) continue;
444 
445  int p2 = cmap->GetPlane(fDigi[j2].get());
446  int c2 = cmap->GetCell (fDigi[j2].get());
447  int v2 = geom->Plane(p2)->View();
448 
449  if (v1==v2) {
450  if (abs(c2-c1)<dcell) dcell = abs(c2-c1);
451  }
452  if (abs(p2-p1)<dplane) dplane = abs(p2-p1);
453  }
454  //
455  // If dcell and dplane are still uninitialized, its because there
456  // is only one hit in this view. Reset the values to 0.
457  //
458  if (dcell == kUNINIT) dcell = 0;
459  if (dplane == kUNINIT) dplane = 0;
460 
461  if (dcell<=fDcellCut && dplane<=fDplaneCut) goodhit[j1] = true;
462  else goodhit[j1] = false;
463  }
464 
465  int nx = 0;
466  int ny = 0;
467  int c1x = 999999, c2x = -999999;
468  int c1y = 999999, c2y = -999999;
469  int p1x = 999999, p2x = -999999;
470  int p1y = 999999, p2y = -999999;
471  double x=0, y=0, z=0, t=0;
472  double qx=0, qy=0;
473  double xyz[3];
474  for (i=i1; i<=i2; ++i) {
475  if (goodhit[i]==false) continue;
476 
477  int p = cmap->GetPlane(fDigi[i].get());
478  int c = cmap->GetCell (fDigi[i].get());
479 
480  geo::View_t v;
481  geom->CellInfo(p, c, &v, xyz);
482 
483  if (v==geo::kX) {
484  ++nx;
485  x += xyz[0];
486  qx += fDigi[i]->ADC();
487  if (c<c1x) c1x = c;
488  if (p<p1x) p1x = p;
489  if (c>c2x) c2x = c;
490  if (p>p2x) p2x = p;
491  }
492  if (v==geo::kY) {
493  ++ny;
494  y += xyz[1];
495  qy += fDigi[i]->ADC();
496  if (c<c1y) c1y = c;
497  if (p<p1y) p1y = p;
498  if (c>c2y) c2y = c;
499  if (p>p2y) p2y = p;
500  }
501  z += xyz[2];
502  t += fDigi[i]->TDC()*kUSEC_PER_TDC;
503  }
504  //
505  // Require at least one good hit in each view
506  //
507  if (nx<1 || ny<1) return;
508 
509  fNt->iwin = iwin;
510  fNt->nx = nx;
511  fNt->ny = ny;
512  fNt->x = x = x/(float)nx;
513  fNt->y = y = y/(float)ny;
514  fNt->z = z = z/(float)(nx+ny);
515  fNt->t = t = t/(float)(nx+ny);
516  fNt->qx = qx;
517  fNt->qy = qy;
518 
519  fNt->dcellx = c2x-c1x;
520  fNt->dcelly = c2y-c1y;
521  fNt->dplanex = p2x-p1x;
522  fNt->dplaney = p2y-p1y;
523 
524  fNt->compx = (float)nx/(float)((fNt->dcellx+1)*(fNt->dplanex+1));
525  fNt->compy = (float)ny/(float)((fNt->dcelly+1)*(fNt->dplaney+1));
526 
527  //
528  // Find the closest muon to this cluster
529  //
530  double dstop, tstop; // Distance and time to stop muon end point
531  double dperp, tperp; // Distance and time to muon track
532  unsigned int istop, iperp; // Indicies of matched muons
533  this->FindClosestMuon(iwin,
534  x, y, z, t,
535  &dstop, &tstop, &istop,
536  &dperp, &tperp, &iperp);
537  fNt->dstop = dstop;
538  fNt->tstop = tstop;
539  fNt->dperp = dperp;
540  fNt->tperp = tperp;
541 
542  //
543  // Find the closest muon to this cluster in an unassociated event
544  //
545  unsigned int istopna, iperpna;
546  this->FindClosestMuon(iwin,
547  x, y, z, t,
548  &dstop, &tstop, &istopna,
549  &dperp, &tperp, &iperpna,
550  false);
551  fNt->dstopna = dstop;
552  fNt->tstopna = tstop;
553  fNt->dperpna = dperp;
554  fNt->tperpna = tperp;
555 
556  if (fDoFill) fNtTree->Fill();
557 
558  //
559  // Check if this is a Michel candidate, if yes, print it to a log
560  // file
561  //
562  if (fDumpMichel==1) {
563  bool ismichel = true;
564  if ( fNt->dplanex>=(fNt->nx+fNt->ny) ) ismichel = false;
565  if ( fNt->dplaney>=(fNt->nx+fNt->ny) ) ismichel = false;
566  if ( fNt->dcellx >=(fNt->nx+fNt->ny) ) ismichel = false;
567  if ( fNt->dcelly >=(fNt->nx+fNt->ny) ) ismichel = false;
568  if ( 4*fNt->nx <= (fNt->nx+fNt->ny) ) ismichel = false;
569  if ( 4*fNt->ny <= (fNt->nx+fNt->ny) ) ismichel = false;
570  if ( fNt->dstop > 30.0 ) ismichel = false;
571  if ( fNt->tstop > 10.0 ) ismichel = false;
572  if (ismichel) this->DumpMichel(istop);
573  }
574  }
575 
576  //......................................................................
577 
578  void Leana::MuonCandidate(unsigned int iwin,
579  unsigned int sz,
580  unsigned int nx,
581  unsigned int ny)
582  {
583  if (sz<fMuonNhit || nx<fMuonNhitX || ny<fMuonNhitY) {
584  fMuon[fMuIndx1][iwin].fId = kNotMu;
585  return;
586  }
587 
589 
590  trk::RLFit trkfit(fDigi, fWin[iwin].first, fWin[iwin].second);
591  trkfit.Fit();
592  trkfit.AssumeDownGoing();
593 
594  //
595  // Does this look like an entering track?
596  //
597  bool isentering = false;
598  if (this->IsPointOnEdge(trkfit.fX1,trkfit.fY1,trkfit.fZ1)) {
599  isentering = true;
600  }
601 
602  if (isentering==false) {
603  fMuon[fMuIndx1][iwin].fId = kNotMu;
604  return;
605  }
606 
607  fMuon[fMuIndx1][iwin].fX1 = trkfit.fX1;
608  fMuon[fMuIndx1][iwin].fY1 = trkfit.fY1;
609  fMuon[fMuIndx1][iwin].fZ1 = trkfit.fZ1;
610 
611  fMuon[fMuIndx1][iwin].fX2 = trkfit.fX2;
612  fMuon[fMuIndx1][iwin].fY2 = trkfit.fY2;
613  fMuon[fMuIndx1][iwin].fZ2 = trkfit.fZ2;
614 
615  bool isstopping = false;
616  if (this->IsPointOnEdge(trkfit.fX2,trkfit.fY2,trkfit.fZ2)==false) {
617  isstopping = true;
618  }
619 
620  if (isstopping) fMuon[fMuIndx1][iwin].fId = kStopMu;
621  else fMuon[fMuIndx1][iwin].fId = kThruMu;
622 
623  unsigned int i;
624  double sumt = 0.0;
625  double sumw = 0.0;
626  for (i=0; i<trkfit.fXhit.size(); ++i) {
627  sumw += trkfit.fWx[i];
628  sumt += trkfit.fWx[i]*trkfit.fXhit[i]->TDC()*kUSEC_PER_TDC;
629  }
630  for (i=0; i<trkfit.fYhit.size(); ++i) {
631  sumw += trkfit.fWy[i];
632  sumt += trkfit.fWy[i]*trkfit.fYhit[i]->TDC()*kUSEC_PER_TDC;
633  }
634 
635  fMuon[fMuIndx1][iwin].fNOvAT = sumt/sumw;
636  fMuon[fMuIndx1][iwin].fT = fMuon[fMuIndx1][iwin].fNOvAT;
637  }
638 
639  //......................................................................
640 
641  void Leana::FindClosestMuon(unsigned int iwin,
642  double x, double y, double z, double t,
643  double* dend, double* tend, unsigned int* istop,
644  double* dtrk, double* ttrk, unsigned int* iperp,
645  bool assn)
646  {
647  //
648  // Decide if we're looking for real associations (assn=true) or are
649  // looking for random coinicidences
650  //
651  unsigned int imu = 0;
652  if (assn==true) imu = fMuIndx1;
653  else imu = fMuIndx2;
654 
655  unsigned int i;
656  double dstop;
657  double dperp;
658  *dend = 9E99;
659  *tend = 9E99;
660  *dtrk = 9E99;
661  *ttrk = 9E99;
662  *istop = 999999;
663  *iperp = 999999;
664  for (i=0; i<fMuon[imu].size(); ++i) {
665  if (fMuon[imu][i].fId==kNotMu) continue;
666  if (fMuon[imu][i].fT>=t) continue;
667 
668  if (fMuon[imu][i].fId==kStopMu) {
669  //
670  // Find distance to stop mu end point
671  //
672  dstop =
673  util::pythag((x-fMuon[imu][i].fX2),
674  (y-fMuon[imu][i].fY2),
675  (z-fMuon[imu][i].fZ2));
676  if (dstop<*dend) {
677  *dend = dstop;
678  *tend = t-fMuon[imu][i].fT;
679  *istop = i;
680  }
681  }
682  //
683  // Find perpendicular distance to muon track
684  //
685  double point[3] = {x,y,z};
686  double icept[3] = { fMuon[imu][i].fX1,
687  fMuon[imu][i].fY1,
688  fMuon[imu][i].fZ1
689  };
690  double slope[3] = {
691  fMuon[imu][i].fX2 - fMuon[imu][i].fX1,
692  fMuon[imu][i].fY2 - fMuon[imu][i].fY1,
693  fMuon[imu][i].fZ2 - fMuon[imu][i].fZ1,
694  };
695  double closest[3] = {-999999.,-999999.,-999999.};
696  dperp = geo::ClosestApproach(point, icept, slope, closest);
697  //
698  // If the point happens outside the limits of the track, compute
699  // the distance to the track end points
700  //
701  if (closest[1]<fMuon[imu][i].fY2) {
702  dperp = util::pythag((x-fMuon[imu][i].fX2),
703  (y-fMuon[imu][i].fY2),
704  (z-fMuon[imu][i].fZ2));
705  }
706  if (closest[1]>fMuon[imu][i].fY1) {
707  dperp = util::pythag((x-fMuon[imu][i].fX1),
708  (y-fMuon[imu][i].fY1),
709  (z-fMuon[imu][i].fZ1));
710  }
711  if (dperp<*dtrk) {
712  *dtrk = dperp;
713  *ttrk = t-fMuon[imu][i].fT;
714  *iperp = i;
715  }
716  }
717  }
718 
719  //......................................................................
720 
721  bool Leana::IsPointOnEdge(double x, double y, double z)
722  {
724 
725  if (y>geom->DetHalfHeight()-fEdgeBuffer) return true;
726  if (fabs(x)>geom->DetHalfWidth()-fEdgeBuffer) return true;
727  if (z<fEdgeBuffer) return true;
728  //
729  //! \todo Hardwired to the NDOS case...
730  //
731  if (z>810.0-fEdgeBuffer) return true;
732  // if (z>geom->DetLength()-fEdgeBuffer) return true;
733  return false;
734  }
735 } // end namespace comi
736 ////////////////////////////////////////////////////////////////////////
738 ////////////////////////////////////////////////////////////////////////
Leana(fhicl::ParameterSet const &pset)
void FindClosestMuon(unsigned int i, double x, double y, double z, double t, double *dstop, double *tstop, unsigned int *istop, double *dperp, double *tperp, unsigned int *iperp, bool assn=true)
SubRunNumber_t subRun() const
Definition: Event.h:72
unsigned int fNhitX
Number of hits req&#39;d in x view.
float fX1
Start x position (cm)
Definition: Leana_module.cc:38
TH1D * f3HitSamePlaneDcell
Max-min cell in same plane.
Definition: Leana_module.cc:94
TH1D * f1HitADC
ADC spectrum for 1 hit slices.
Definition: Leana_module.cc:87
float tstopna
time difference to muon stop point IDed above
Definition: LeanaNt.h:45
int dcellx
range of cells in x view
Definition: LeanaNt.h:20
unsigned int fMuonNhitY
required in y view
int dplaney
range of cells in y view
Definition: LeanaNt.h:23
double fY2
End y positon (cm)
Definition: RLFit.h:109
bool fDoFill
Should we be filling ntuple / histogram?
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
#define H1(NAME, LABEL, N, X1, X2)
TH2D * f1HitADCvsTDC
ADC vs TDC for 1 hit slices.
Definition: Leana_module.cc:89
Summary ntuple for low energy (<100 MeV) events.
int run
run number
Definition: LeanaNt.h:14
LeanaNt * fNt
Ntuple for hit clusters.
Definition: Leana_module.cc:97
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::vector< const rawdata::RawDigit * > fYhit
Y view hits.
Definition: RLFit.h:92
float fZ1
Start z position (cm)
Definition: Leana_module.cc:40
static const unsigned int MUON_BUFF_SZ
int ny
number of hits in y view
Definition: LeanaNt.h:19
const char * p
Definition: xmltok.h:285
float fX2
End x position (cm)
Definition: Leana_module.cc:41
std::vector< double > fWy
y view weights
Definition: RLFit.h:98
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< const rawdata::RawDigit * > fXhit
X view hits.
Definition: RLFit.h:91
static const char kThruMu
A through-going muon candidate.
#define H2(NAME, LABEL, NX, X1, X2, NY, Y1, Y2)
TH1D * fEventCount
Event counter.
Definition: Leana_module.cc:83
float dstop
minimum distance to a muon stop point
Definition: LeanaNt.h:35
float fY1
Start y position (cm)
Definition: Leana_module.cc:39
void abs(TH1 *hist)
int fDumpMichel
Dump table of Michel electrons?
TH1D * f1HitTDC
TDC spectrum for 1 hit slices.
Definition: Leana_module.cc:88
std::vector< double > fWx
x view weights
Definition: RLFit.h:95
float tperpna
time difference to muon track IDed above
Definition: LeanaNt.h:42
float compy
"completeness" in y view
Definition: LeanaNt.h:25
const PlaneGeo * Plane(unsigned int i) const
int nx
number of hits in x view
Definition: LeanaNt.h:18
DEFINE_ART_MODULE(TestTMapFile)
float fY2
End y position (cm)
Definition: Leana_module.cc:42
float tstop
time difference to muon stop point IDed above
Definition: LeanaNt.h:36
float tperp
time difference to muon track IDed above
Definition: LeanaNt.h:33
double fT0
Event NOvA trigger time (usec)
float z
location of collection (cm)
Definition: LeanaNt.h:28
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
unsigned int fMuIndx2
Old non-associated muons.
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
void CountXY(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int i1, unsigned int i2, unsigned int *nx, unsigned int *ny)
Count the number of digits in each detector view.
Definition: RawUtil.cxx:38
double fX1
Start x position (cm)
Definition: RLFit.h:105
float dperp
minimum perpendicular distance to a muon track
Definition: LeanaNt.h:32
Summary information about previous muons.
Definition: Leana_module.cc:36
c2
Definition: demo5.py:33
double fY1
Start y position (cm)
Definition: RLFit.h:106
float dstopna
minimum distance to a muon stop point
Definition: LeanaNt.h:44
Analysis module for low energy events and noise.
Definition: Leana_module.cc:49
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
void AssumeDownGoing()
Definition: RLFit.cxx:247
double fZ1
Start z position (cm)
Definition: RLFit.h:107
TH2D * fSliceSzVsEvent
Histogram of slice size vs. event.
Definition: Leana_module.cc:85
int dcelly
range of cells in y view
Definition: LeanaNt.h:21
void beginJob()
void reconfigure(const fhicl::ParameterSet &p)
Make a linear fit to a collection of raw hits.
bool IsPointOnEdge(double x, double y, double z)
void Fit()
Definition: RLFit.cxx:110
Commissioning files to look at the quality of our data.
Definition: Cana_module.cc:39
T get(std::string const &key) const
Definition: ParameterSet.h:231
float compx
"completeness" in x view
Definition: LeanaNt.h:24
void TimeSlice(const std::vector< art::Ptr< rawdata::RawDigit > > &d, unsigned int dt_tdc, unsigned int nhit, unsigned int nhitx, unsigned int nhity, std::vector< RawSlice > &slice)
Find windows in time that have significant activity in the detector.
Definition: RawUtil.cxx:111
TH1D * f2HitDplane
Distance from first to last plane.
Definition: Leana_module.cc:91
unsigned int fNhitY
Number of hits req&#39;d in y view.
Collect Geo headers and supply basic geometry functions.
unsigned int fNhit
Number of hits req&#39;d in window.
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
float y
location of collection (cm)
Definition: LeanaNt.h:27
void AnaNHit(unsigned int iwin)
double fDeltaTnamu
Time for non-association (s)
Functions for getting information about a collection of raw hits.
TH1D * fSliceSz
Histogram of Nhit per slice.
Definition: Leana_module.cc:84
double DetHalfHeight() const
static const char kStopMu
A stopping muon candidate.
unsigned char fId
Stop? Through-going?
Definition: Leana_module.cc:37
z
Definition: test.py:28
std::vector< util::RawSlice > fWin
The raw slices for this event.
float x
location of collection (cm)
Definition: LeanaNt.h:26
double fT
Time in spill (usec)
Definition: Leana_module.cc:44
void DumpMichel(unsigned int i)
OStream cout
Definition: OStream.cxx:6
int fDcellCut
Allowed gap between hit cells.
int dplanex
range of cells in x view
Definition: LeanaNt.h:22
EventNumber_t event() const
Definition: EventID.h:116
float fZ2
End z position (cm)
Definition: Leana_module.cc:43
int subrun
subrun number
Definition: LeanaNt.h:15
double fEdgeBuffer
[cm] Edge distance for stopping/entering decision
virtual ~Leana()
T * make(ARGS...args) const
float qy
charge of collection in y view (ADC)
Definition: LeanaNt.h:30
double fNOvAT
"Global" NOvA time (usec)
Definition: Leana_module.cc:45
double DetHalfWidth() const
void MuonCandidate(unsigned int iwin, unsigned int sz, unsigned int nx, unsigned int ny)
static const int kTDC_PER_USEC
Definition: Cana_module.cc:87
TH1D * f3HitDplane
Distance from first to last plane.
Definition: Leana_module.cc:92
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
double fZ2
End z position (cm)
Definition: RLFit.h:110
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void DumpMuon(unsigned int i)
void beginJob()
void Ana1Hit(unsigned int iwin)
void geom(int which=0)
Definition: geom.C:163
double fTimeWindow
[usec] Time window size
int Apply(art::PtrVector< rawdata::RawDigit > &rd)
static const double kUSEC_PER_TDC
TTree * fNtTree
Ntuple for hit clusters.
Definition: Leana_module.cc:98
std::vector< art::Ptr< rawdata::RawDigit > > fDigi
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
c1
Definition: demo5.py:24
unsigned int fMuonNhit
Required # hits on muon track.
int iwin
timing window number inside event
Definition: LeanaNt.h:17
void DumpMuonTable()
unsigned int fMuonNhitX
required in x view
static const char kNotMu
Not a good muon candidate.
float qx
charge of collection in x view (ADC)
Definition: LeanaNt.h:29
RunNumber_t run() const
Definition: Event.h:77
unsigned long long TDCT0() const
Return just the lower 32 bits of the timing marker. This is the reference "t0" for the RawDigit TDC c...
Definition: RawTrigger.cxx:44
void TimeSort(std::vector< art::Ptr< rawdata::RawDigit > > &d)
Arrange the list of raw hits in time order (early to late)
Definition: RawUtil.cxx:31
void analyze(art::Event const &evt)
Ntuple produced by the leana module.
Definition: LeanaNt.h:12
float t
time of collection (usec)
Definition: LeanaNt.h:31
unsigned int fMuIndx1
Current muon list.
Definition: fwd.h:28
std::vector< MuonCand > fMuon[MUON_BUFF_SZ]
Table of muon data.
int evt
event number
Definition: LeanaNt.h:16
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
double fX2
End x position (cm)
Definition: RLFit.h:108
float dperpna
minimum perpendicular distance to a muon track
Definition: LeanaNt.h:41
TH1D * f3HitDiffPlaneDcell
Max-min cell in different planes.
Definition: Leana_module.cc:95
int fDplaneCut
Allowed gap between hit planes.