Cana_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 /// \file Cana_module.cc
3 /// \brief Commisionning analysis module - looks for interesting events
4 /// \author messier@indiana.edu
5 ////////////////////////////////////////////////////////////////////////
6 #include <string>
7 #include <vector>
8 
9 #include "TH1F.h"
10 #include "TH2F.h"
11 #include "TTree.h"
12 
13 // ART includes
19 
20 // NOvA includes
24 #include "Commissioning/CanaNt.h"
25 #include "Geometry/Geometry.h"
27 #include "RawData/RawDigit.h"
28 #include "RawData/RawTrigger.h"
29 #include "RecoBase/Prong.h"
32 #include "SummaryData/POTSum.h"
33 #include "SummaryData/SpillData.h"
34 #include "TrackFit/RLFit.h"
35 #include "Utilities/func/MathUtil.h"
36 #include "Utilities/RawUtil.h"
37 
38 
39 namespace comi {
40 
41  /// Analysis module to hunt for neutrino candidates
42  class Cana : public art::EDProducer {
43  public:
44  explicit Cana(fhicl::ParameterSet const& pset);
45  virtual ~Cana();
46 
47  void beginJob();
48  void endSubRun(art::SubRun &sr);
49 
50  void produce(art::Event& evt);
51  void reconfigure(const fhicl::ParameterSet& p);
52 
53  private:
54  void WipeNt();
55  void LoadMCInfo(art::Event& evt);
56 
57  private:
58  TTree* fNtTree;
59  TTree* fPotTree;
62  TH1F* fNhitTotal;
63  //TH1F* fNSlices;
64  TH1F* fNhitInTime;
65  TH1F* fNhitOutTime;
66  TH1F* fNtimeWin;
67  TH1F* fNhitTimeWin;
68  TH1F* fHitTimesAll;
72 
73  private:
74  int fApplyBadChan; ///< Apply bad channel masks?
75  std::string fRawDataLabel; ///< Module that produced the raw data
76  std::string fMCLabel; ///< Module that produced the mc truth
77  std::string fSpillLabel; ///< Module that produced the spill info
78  std::string fPotLabel; ///< Module that produced the POTSum object
79  float fTimeWindow; ///< Size of sliding time window
80  unsigned int fNhit; ///< Total number of intime hits req'd on track
81  unsigned int fNhitX; ///< Required number of hits in window on track
82  unsigned int fNhitY; ///< Required number of hits in window on track
83  };
84 
85  // NOvA uses a 64 Mhz clock
86  //static const double kUSEC_PER_TDC = (1.0E6/64.0E6);
87  static const int kTDC_PER_USEC = 64;
88 
89  //......................................................................
90 
92  {
93  this->reconfigure(p);
94  this->produces< std::vector<CanaNt> >();
95  this->produces< std::vector<rb::Prong> >();
96  }
97 
98  //......................................................................
99 
101  {
103 
104  // Some event-level histograms
105  fNhitTotal = f->make<TH1F>("fNhitTotal",
106  ";Nhit;Events",
107  100, 0.0, 50000.0);
108 
109  fNhitInTime = f->make<TH1F>("fNhitInTime",
110  ";Nhit;Events",
111  100, 0.0, 50000.0);
112 
113  fNhitOutTime = f->make<TH1F>("fNhitOutTime",
114  ";Nhit;Events",
115  100,0.0,50000.0);
116 
117  fNtimeWin = f->make<TH1F>("fNtimeWin",
118  ";N time slices;Events",
119  100, 0.0, 1000.0);
120 
121  fNhitTimeWin = f->make<TH1F>("fNhitTimeWin",
122  ";Nhit;Time slices",
123  200, 0.0, 2000.0);
124 
125  fHitTimesAll = f->make<TH1F>("fHitTimesAll",
126  ";time (usec);hits / 100 usec",
127  610, -1000.0, 60000.0);
128 
129  fHitTimesGood = f->make<TH1F>("fHitTimesGood",
130  ";time (uses);hits / 100 usec",
131  610, -1000.0, 60000.0);
132 
133  fNhitTotalVsEvent = f->make<TH2F>("fNhitTotalVsEvent",
134  ";Event Number;Nhit",
135  500,0.0,500000.0,
136  100,0.0,50000.0);
137 
138  fNtimeWinVsEvent = f->make<TH2F>("fNtimeWinVsEvent",
139  ";Event Number; Nwin",
140  500,0.0,50000.0,
141  50,0.0,500.0);
142  fNt = new CanaNt;
143  fNtTree = f->make<TTree>("nt","Cana Analysis Tree");
144  fNtTree->Branch("fCanaNt", "comi::CanaNt", &fNt, 1600, 99);
145 
146  fPOTSum = new sumdata::POTSum();
147  fPotTree = f->make<TTree>("pottree","Total POT Tree");
148  fPotTree->Branch("fPOTSum","sumdata::POTSum", &fPOTSum,1600,99);
149  }
150 
151  //......................................................................
152 
154 
155  //......................................................................
156 
158  {
160  try {
161  sr.getByLabel(fPotLabel,p);
162  fPOTSum->totpot = p->totpot;
166  fPotTree->Fill();
167  }
168  catch (...) {
169  }
170  }
171 
172  //......................................................................
173 
175  {
176  unsigned int i;
177  std::unique_ptr< std::vector<CanaNt> > ntvec(new std::vector<CanaNt>);
178  std::unique_ptr< std::vector<rb::Prong> > prong(new std::vector<rb::Prong>);
179 
181 
182  this->WipeNt();
183 
184  //
185  // Load event ID information into ntuple
186  //
187  fNt->run = evt.run();
188  fNt->subrun = evt.subRun();
189  fNt->evt = evt.id().event();
190 
191  //
192  // ART time stamps are 32-bits of UNIX-seconds (high part) and
193  // 32-bits of nanoseconds (low part)
194  //
195  fNt->utc.SetSec( (evt.time().value()>>32)&0xFFFFFFFFUL );
196  fNt->utc.SetNanoSec(evt.time().value() &0xFFFFFFFFUL );
197 
198  //
199  // If there is information from the simulation, load it
200  //
201  this->LoadMCInfo(evt);
202 
203  //
204  // Get trigger information to extract trigger type
205  //
207  evt.getByLabel(fRawDataLabel, trigv);
208  const rawdata::RawTrigger& trig = (*trigv)[0];
209 
212 
213  //
214  // Get spill information
215  //
217  try {
218  evt.getByLabel(fSpillLabel, spill);
219  fNt->pot = spill->spillpot;
220  fNt->dt = spill->deltaspilltimensec/1000.;
221  fNt->gbeam=spill->goodbeam;
222  }
223  catch (...) { }
224 
225  //
226  // Pull out the raw digits from the event
227  //
229  evt.getByLabel(fRawDataLabel, digidummy);
230 
231  //
232  // Shuffle the data into a vector we can manipulate
233  //
234  std::vector< art::Ptr<rawdata::RawDigit> > digi;
235  for (i=0; i<digidummy->size(); ++i) {
236  digi.push_back(art::Ptr<rawdata::RawDigit>(digidummy, i));
237  }
238  if (digi.size()==0) {
239  //
240  // Make sure each event gets at least one entry in the ntuple even
241  // if we failed to find anything interesting at all
242  //
243  if (ntvec->size()==0) ntvec->push_back(*fNt);
244  evt.put(std::move(ntvec));
245  evt.put(std::move(prong));
246  fNtTree->Fill();
247  return;
248  }
249 
250  //
251  // Fill the histogram showing the times of all hits
252  //
253  for (i=0; i<digi.size(); ++i) {
254  bool goodtime;
255  fHitTimesAll->Fill(cal->GetTNS(digi[i].get(), goodtime)/1000);
256  }
257 
258  unsigned int tdcwin = (unsigned int)rint(fTimeWindow*kTDC_PER_USEC);
259  double t;
260 
261  //
262  // Filter out the bad channels for this run
263  //
264  if (fApplyBadChan!=0) {
266  badc->Apply(digi);
267  }
268 
269  if (digi.size()==0) {
270  //
271  // Make sure each event gets at least one entry in the ntuple even
272  // if we failed to find anything interesting at all
273  //
274  if (ntvec->size()==0) ntvec->push_back(*fNt);
275  evt.put(std::move(ntvec));
276  evt.put(std::move(prong));
277  fNtTree->Fill();
278  return;
279  }
280  for (i=0; i<digi.size(); ++i) {
281  bool goodtime;
282  fHitTimesGood->Fill(cal->GetTNS(digi[i].get(), goodtime)/1000);
283  }
284 
285  //
286  // Time sort the hits and look for in-time activity
287  //
288  std::vector<util::RawSlice> win;
289  util::TimeSort(digi);
290  util::TimeSlice(digi,tdcwin,fNhit,fNhitX,fNhitY,win);
291 
292  // Log some information about the numbers of hits and the number of
293  // windows
294  fNtimeWin->Fill(win.size());
295 
296  unsigned int nhittot = digi.size();
297  unsigned int nhitwin = 0;
298  unsigned int nhitinwin = 0;
299  fNhitTotal->Fill(nhittot);
300  fNhitTotalVsEvent->Fill(fNt->evt, nhittot);
301  fNtimeWinVsEvent->Fill (fNt->evt, win.size());
302  for (i=0; i<win.size(); ++i) {
303  nhitwin = (win[i].second+1)-win[i].first;
304  fNhitTimeWin->Fill(nhitwin);
305  nhitinwin += nhitwin;
306  }
307  fNhitInTime->Fill(nhitinwin);
308  fNhitOutTime->Fill(nhittot-nhitinwin);
309 
310  //
311  // Evaluate the activity inside the time window
312  //
313  double tlast = -999999.;
314  for (i=0; i<win.size(); ++i) {
315  double sumt = 0.0;
316  double sumt2 = 0.0;
317  double sumw = 0.0;
318  double wcut = 0.01;
319  fNt->iwin = i;
320  fNt->nhit = 0;
321  fNt->nhitx = 0;
322  fNt->nhity = 0;
323  fNt->nunused = 0;
324  fNt->t1 = 999999.;
325  fNt->t2 =-999999.;
326  fNt->cell1x = 999999;
327  fNt->cell1y = 999999;
328  fNt->plane1x = 999999;
329  fNt->plane1y = 999999;
330  fNt->cell2x =-999999;
331  fNt->cell2y =-999999;
332  fNt->plane2x =-999999;
333  fNt->plane2y =-999999;
334 
335  //
336  // Skip the line fitting if we don't have enough hits in each view
337  //
338  unsigned int nx, ny;
339  util::CountXY(digi, win[i].first, win[i].second, &nx, &ny);
340  if (nx<2 || ny<2) continue;
341 
342  //
343  // Make a straight line fit to the hits in the window
344  //
345  trk::RLFit line(digi, win[i].first, win[i].second);
346  line.Fit();
347 
348  //
349  // Compute some information about the track fit
350  //
352  for (unsigned int j=0; j<line.fXhit.size(); ++j) {
353  //
354  // Check weight to see if hit contributed to the track
355  //
356  if (line.fWx[j]>wcut) {
357  ++fNt->nhit;
358  ++fNt->nhitx;
359 
360  int p = cmap->GetPlane(line.fXhit[j]);
361  int c = cmap->GetCell (line.fXhit[j]);
362  if (p<fNt->plane1x) fNt->plane1x = p;
363  if (p>fNt->plane2x) fNt->plane2x = p;
364  if (c<fNt->cell1x) fNt->cell1x = c;
365  if (c>fNt->cell2x) fNt->cell2x = c;
366 
367  //
368  // Record the hit time in units of usec
369  //
370  bool goodtime;
371  t = cal->GetTNS(line.fXhit[j], goodtime)/1000;
372 
373  if (t<fNt->t1) fNt->t1 = t;
374  if (t>fNt->t2) fNt->t2 = t;
375 
376  sumt += line.fWx[j]*t;
377  sumt2 += line.fWx[j]*t*t;
378  sumw += line.fWx[j];
379  }
380  else {
381  ++fNt->nunused;
382  }
383  }
384  for (unsigned int j=0; j<line.fYhit.size(); ++j) {
385  //
386  // Check weight to see if hit contributed to the track
387  //
388  if (line.fWy[j]>wcut) {
389  ++fNt->nhit;
390  ++fNt->nhity;
391 
392  int p = cmap->GetPlane(line.fYhit[j]);
393  int c = cmap->GetCell (line.fYhit[j]);
394  if (p<fNt->plane1y) fNt->plane1y = p;
395  if (p>fNt->plane2y) fNt->plane2y = p;
396  if (c<fNt->cell1y) fNt->cell1y = c;
397  if (c>fNt->cell2y) fNt->cell2y = c;
398 
399  //
400  // Record the hit time in units of usec
401  //
402  bool goodtime;
403  t = cal->GetTNS(line.fYhit[j], goodtime)/1000;
404 
405  if (t<fNt->t1) fNt->t1 = t;
406  if (t>fNt->t2) fNt->t2 = t;
407 
408  sumt += line.fWy[j]*t;
409  sumt2 += line.fWy[j]*t*t;
410  sumw += line.fWy[j];
411  }
412  else {
413  ++fNt->nunused;
414  }
415  }
416  // Leave at -999999 if there is somehow no weight
417  if (sumw>0.0) {
418  fNt->tave = sumt / sumw;
419  fNt->trms = sqrt(sumt2-sumt*sumt/sumw)/(sumw-1.0);
420  fNt->tprev = fNt->tave - tlast;
421  tlast = fNt->tave;
422  }
423 
424  fNt->cosdet = line.CosThetaBeam();
425  fNt->phidet = (180.0/M_PI)*line.PhiBeam();
426 
427  fNt->coscosmic = line.CosThetaCosmic();
428  fNt->phicosmic = (180.0/M_PI)*line.PhiCosmic();
429 
431  //
432  // Approximate directions of the BNB beam in the NDOS
433  // frame. NuMI direction from GeometryBase::NuMIBeamDirection(). BNB direction from
434  // Zelimir Djurcic <zdjurcic@hep.anl.gov> 2/26/2011, private
435  // communication.
436  //
437  static const double bnbdir[] = { -0.367, 0.012, 0.930};
438  static const double numidir[] = {geom->NuMIBeamDirection().X(), geom->NuMIBeamDirection().Y(), geom->NuMIBeamDirection().Z()};
439 
440  //
441  // Both the Booster and NuMI beamlines send neutrinos generally
442  // north so analyze the track with a north-going assumption
443  //
444  line.AssumeNorthGoing();
445 
446  //
447  // Compute angles wrt the NuMI and Booster beams
448  //
449  static const double magbnb = util::pythag(bnbdir[0],bnbdir[1],bnbdir[2]);
450  static const double magnumi= util::pythag(numidir[0],numidir[1],numidir[2]);
451 
452  double magtrk = util::pythag((line.fX2-line.fX1),
453  (line.fY2-line.fY1),
454  (line.fZ2-line.fZ1));
455 
456  fNt->cosbnb =
457  ((line.fX2-line.fX1)*bnbdir[0] +
458  (line.fY2-line.fY1)*bnbdir[1] +
459  (line.fZ2-line.fZ1)*bnbdir[2]) / magbnb / magtrk;
460 
461  fNt->cosnumi =
462  ((line.fX2-line.fX1)*numidir[0] +
463  (line.fY2-line.fY1)*numidir[1] +
464  (line.fZ2-line.fZ1)*numidir[2]) / magnumi / magtrk;
465 
466  //
467  // Record projected distances to detector edges
468  //
469  bool indetector = true;
470  if (line.fX1 < -geom->DetHalfWidth()) indetector = false;
471  if (line.fX1 > +geom->DetHalfWidth()) indetector = false;
472  if (line.fY1 < -geom->DetHalfHeight()) indetector = false;
473  if (line.fY1 > +geom->DetHalfHeight()) indetector = false;
474  if (line.fZ1 < 0) indetector = false;
475  if (line.fZ1 > geom->DetLength()) indetector = false;
476  if (indetector) {
477  double x[3] = {line.fX1, line.fY1, line.fZ1};
478  double dxf[3] = {line.fX2-line.fX1,
479  line.fY2-line.fY1,
480  line.fZ2-line.fZ1};
481  double dxb[3] = {-dxf[0], -dxf[1], -dxf[2]};
482  fNt->dedgef = geom->DEdge(x, dxf);
483  fNt->dedgeb = geom->DEdge(x, dxb);
484  }
485 
486  //
487  // Now analyze the track assuming its a down-going cosmic ray
488  //
489  line.AssumeDownGoing();
490  fNt->x1 = line.fX1;
491  fNt->y1 = line.fY1;
492  fNt->z1 = line.fZ1;
493  fNt->x2 = line.fX2;
494  fNt->y2 = line.fY2;
495  fNt->z2 = line.fZ2;
496 
497  fNtTree->Fill();
498  ntvec->push_back(*fNt);
499 
500  //
501  // Store a prong to show the cana fit results
502  // TODO CJB - prong doesn't store the endpoint which Cana knows
503  // but then Cana doesn't store the hits anyway, so we're equally
504  // at fault.
505  //
506  TVector3 xyz1(fNt->x1, fNt->y1, fNt->z1);
507  TVector3 xyz2(fNt->x2, fNt->y2, fNt->z2);
508  rb::Prong p(art::PtrVector<rb::CellHit>(), xyz1, xyz2-xyz1);
509  prong->push_back(p);
510  } // end loop i=0...win.size()
511 
512  //
513  // Make sure that every event gets at least one entry in the ntuple
514  //
515  if (win.size()==0) fNtTree->Fill();
516 
517  //
518  // Store the ntuple and the prongs into the event
519  //
520  evt.put(std::move(ntvec));
521  evt.put(std::move(prong));
522  }
523 
524  //......................................................................
525 
527  {
528  fApplyBadChan = p.get<int>("ApplyBadChan");
529  fTimeWindow = p.get<float>("TimeWindow");
530  fNhit = p.get<int>("Nhit");
531  fNhitX = p.get<int>("NhitX");
532  fNhitY = p.get<int>("NhitY");
533  fRawDataLabel = p.get<std::string>("RawDataLabel", "daq");
534  fPotLabel = p.get<std::string>("PotLabel");
535  fMCLabel = p.get<std::string>("MCLabel", "generator");
536  fSpillLabel = p.get<std::string>("SpillLabel");
537  }
538 
539  //......................................................................
540 
541  ///
542  /// Fill the ntuple with clearly invalid data
543  ///
544  void Cana::WipeNt()
545  {
546  fNt->run = -99;
547  fNt->subrun = -99;
548  fNt->evt = -99;
549  fNt->iwin = -99;
550  fNt->trgsrc = -99;
551  fNt->nhit = -99;
552  fNt->nhitx = -99;
553  fNt->nhity = -99;
554  fNt->plane1x = 999999;
555  fNt->plane2x = -999999;
556  fNt->plane1y = 999999;
557  fNt->plane2y = -999999;
558  fNt->cell1x = 999999;
559  fNt->cell2x = -999999;
560  fNt->cell1y = 999999;
561  fNt->cell2y = -999999;
562  fNt->t1 = 999999;
563  fNt->t2 = -999999;
564  fNt->tave = -999999;
565  fNt->trms = -999999;
566  fNt->cosnumi = -999999;
567  fNt->cosbnb = -999999;
568  fNt->cosdet = -999999;
569  fNt->phidet = -999999;
570  fNt->coscosmic = -999999;
571  fNt->phicosmic = -999999;
572  fNt->pot = -999999;
573  fNt->dt = -999999;
574  fNt->gbeam = -999999;
575  fNt->x1 = -999999;
576  fNt->y1 = -999999;
577  fNt->z1 = -999999;
578  fNt->x2 = -999999;
579  fNt->y2 = -999999;
580  fNt->z2 = -999999;
581  fNt->dedgef = -999999;
582  fNt->dedgeb = -999999;
583  fNt->nunused = -99;
584  fNt->nux = -999999;
585  fNt->nuy = -999999;
586  fNt->nuz = -999999;
587  fNt->nupx = -999999;
588  fNt->nupy = -999999;
589  fNt->nupz = -999999;
590  fNt->leppx = -999999;
591  fNt->leppy = -999999;
592  fNt->leppz = -999999;
593  fNt->nupdg = -99;
594  fNt->iscc = -99;
595  fNt->imode = -99;
596  fNt->tstart = 0;
597  fNt->utc.Set(0,0,0,0,0,0,0,1,0);
598  }
599 
600  //......................................................................
601 
603  {
604  //
605  // Pull the objects holding MC truth out of the event
606  //
608  try {
609  evt.getByLabel(fMCLabel, mc);
610  //
611  // The code only knows how to deal with 1 and exactly 1 MC truth
612  // particle list
613  //
614  if (mc->size()!=1) return;
615  }
616  catch (...) {
617  //
618  // Read failed. Skip setting MC information.
619  //
620  return;
621  }
622 
623  //
624  // Check if the neutrino kinematics are set. If yes, use them
625  //
626  if ((*mc)[0].NeutrinoSet()) {
627  const simb::MCNeutrino& gen = (*mc)[0].GetNeutrino();
628  const simb::MCParticle& nu = gen.Nu();
629  fNt->nux = nu.Vx();
630  fNt->nuy = nu.Vy();
631  fNt->nuz = nu.Vz();
632  fNt->nupx = nu.Px();
633  fNt->nupy = nu.Py();
634  fNt->nupz = nu.Pz();
635  fNt->nupdg = nu.PdgCode();
636  fNt->iscc = (gen.CCNC()==simb::kCC);
637  fNt->imode = gen.Mode();
638  fNt->leppx = gen.Lepton().Px();
639  fNt->leppy = gen.Lepton().Py();
640  fNt->leppz = gen.Lepton().Pz();
641  return;
642  }
643 
644  //
645  // Other cases fall through to here. Assume the first particle in
646  // the list is the one we're interested in...
647  //
648  const simb::MCParticle& p = (*mc)[0].GetParticle(0);
649  fNt->nux = p.Vx();
650  fNt->nuy = p.Vy();
651  fNt->nuz = p.Vz();
652  fNt->nupx = -999999;
653  fNt->nupy = -999999;
654  fNt->nupz = -999999;
655  fNt->nupdg = p.PdgCode();
656  fNt->iscc = -1;
657  fNt->imode = -1;
658  fNt->leppx = p.Px();
659  fNt->leppy = p.Py();
660  fNt->leppz = p.Pz();
661  }
662 } // end namespace comi
663 ////////////////////////////////////////////////////////////////////////
665 ////////////////////////////////////////////////////////////////////////
void produce(art::Event &evt)
Definition: Cana_module.cc:174
double PhiBeam()
Definition: RLFit.cxx:280
float nupy
MC neutrino y momentum (GeV)
Definition: CanaNt.h:59
Ntuple produced by the Cana module.
TTimeStamp utc
UTC calendar date of event.
Definition: CanaNt.h:68
int evt
event number
Definition: CanaNt.h:18
TH2F * fNtimeWinVsEvent
Definition: Cana_module.cc:71
SubRunNumber_t subRun() const
Definition: Event.h:72
float tave
time where peak occurs (ave over hits) (usec)
Definition: CanaNt.h:34
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
float x2
End position of track (cm)
Definition: CanaNt.h:49
double Py(const int i=0) const
Definition: MCParticle.h:230
float z1
Start position of track (cm)
Definition: CanaNt.h:48
unsigned int fNhitY
Required number of hits in window on track.
Definition: Cana_module.cc:82
double fY2
End y positon (cm)
Definition: RLFit.h:109
float y1
Start position of track (cm)
Definition: CanaNt.h:47
float pot
Proton on target.
Definition: CanaNt.h:43
float phicosmic
Azimuthal angle (point back wrt. vertical)
Definition: CanaNt.h:42
std::vector< const rawdata::RawDigit * > fYhit
Y view hits.
Definition: RLFit.h:92
int cell1x
1st cell in x view
Definition: CanaNt.h:28
int subrun
subrun number
Definition: CanaNt.h:17
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< double > fWy
y view weights
Definition: RLFit.h:98
std::vector< const rawdata::RawDigit * > fXhit
X view hits.
Definition: RLFit.h:91
double Px(const int i=0) const
Definition: MCParticle.h:229
float nupz
MC neutrino z momentum (GeV)
Definition: CanaNt.h:60
int nhitx
Number of x view hits on track in peak time window.
Definition: CanaNt.h:22
double DetLength() const
TH1F * fHitTimesAll
Definition: Cana_module.cc:68
TH1F * fNhitOutTime
Definition: Cana_module.cc:65
float trms
rms width of hits inside of window (usec)
Definition: CanaNt.h:35
void endSubRun(art::SubRun &sr)
Definition: Cana_module.cc:157
std::vector< double > fWx
x view weights
Definition: RLFit.h:95
float GetTNS(const rawdata::RawDigit *rawdigit, bool &goodTime, double *maxadc=0)
float nuz
MC neutrino interaction z location (cm)
Definition: CanaNt.h:57
double PhiCosmic()
Definition: RLFit.cxx:303
int trgsrc
trigger source
Definition: CanaNt.h:20
int nupdg
MC neutrino PDG code.
Definition: CanaNt.h:64
DEFINE_ART_MODULE(TestTMapFile)
unsigned long long tstart
Time start from raw trigger.
Definition: CanaNt.h:67
int iscc
1 = CC interaction, 0 = NC interaction
Definition: CanaNt.h:65
unsigned int fNhit
Total number of intime hits req&#39;d on track.
Definition: Cana_module.cc:80
int iwin
window number
Definition: CanaNt.h:19
int imode
interaction mode
Definition: CanaNt.h:66
#define M_PI
Definition: SbMath.h:34
constexpr TimeValue_t value() const
Definition: Timestamp.h:24
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
TH1F * fNhitTotal
Definition: Cana_module.cc:62
float t2
end time window (usec)
Definition: CanaNt.h:33
double fX1
Start x position (cm)
Definition: RLFit.h:105
uint8_t fTriggerMask_TriggerType
Definition: RawTrigger.h:43
void WipeNt()
Definition: Cana_module.cc:544
int gbeam
good beam
Definition: CanaNt.h:45
TH1F * fNtimeWin
Definition: Cana_module.cc:66
void LoadMCInfo(art::Event &evt)
Definition: Cana_module.cc:602
int cell2x
last cell in x view
Definition: CanaNt.h:29
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
double fY1
Start y position (cm)
Definition: RLFit.h:106
TH1F * fNhitInTime
Definition: Cana_module.cc:64
double DEdge(const double *x0, const double *dx, double *exitp=0) const
void AssumeDownGoing()
Definition: RLFit.cxx:247
double fZ1
Start z position (cm)
Definition: RLFit.h:107
TTree * fPotTree
Definition: Cana_module.cc:59
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
std::string fRawDataLabel
Module that produced the raw data.
Definition: Cana_module.cc:75
Analysis module to hunt for neutrino candidates.
Definition: Cana_module.cc:42
double CosThetaCosmic()
Definition: RLFit.cxx:291
signed long long int deltaspilltimensec
Definition: SpillData.h:24
void AssumeNorthGoing()
Definition: RLFit.cxx:236
std::string fSpillLabel
Module that produced the spill info.
Definition: Cana_module.cc:77
Make a linear fit to a collection of raw hits.
float leppy
MC out going lepton momentum (GeV)
Definition: CanaNt.h:62
float z2
End position of track (cm)
Definition: CanaNt.h:51
void Fit()
Definition: RLFit.cxx:110
unsigned int fNhitX
Required number of hits in window on track.
Definition: Cana_module.cc:81
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
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
float fTimeWindow
Size of sliding time window.
Definition: Cana_module.cc:79
float phidet
Azimuthal angle in detector x-y plane.
Definition: CanaNt.h:40
float dt
usec time difference to spill
Definition: CanaNt.h:44
Ntuple produced by the cana module.
Definition: CanaNt.h:14
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
int plane1y
1st plane in y view
Definition: CanaNt.h:26
float cosdet
Cosine of theta wrt to detector z-axis.
Definition: CanaNt.h:39
const double j
Definition: BetheBloch.cxx:29
float nupx
MC neutrino x momentum (GeV)
Definition: CanaNt.h:58
Functions for getting information about a collection of raw hits.
double DetHalfHeight() const
Cana(fhicl::ParameterSet const &pset)
Definition: Cana_module.cc:91
float y2
End position of track (cm)
Definition: CanaNt.h:50
int nhit
Number of hits on track in peak time window.
Definition: CanaNt.h:21
float x1
Start position of track (cm)
Definition: CanaNt.h:46
double CosThetaBeam()
Definition: RLFit.cxx:268
float t1
start time window (usec)
Definition: CanaNt.h:32
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
int nunused
Number of hits in time window but not on the track.
Definition: CanaNt.h:54
EventNumber_t event() const
Definition: EventID.h:116
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
TH1F * fHitTimesGood
Definition: Cana_module.cc:69
CanaNt * fNt
Definition: Cana_module.cc:60
A Cluster with defined start position and direction.
Definition: Prong.h:19
float dedgef
distance to detector edge going north along track (cm)
Definition: CanaNt.h:52
double DetHalfWidth() const
sumdata::POTSum * fPOTSum
Definition: Cana_module.cc:61
static const int kTDC_PER_USEC
Definition: Cana_module.cc:87
double fZ2
End z position (cm)
Definition: RLFit.h:110
float leppz
MC out going lepton momentum (GeV)
Definition: CanaNt.h:63
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
unsigned long long fTriggerTimingMarker_TimeStart
Definition: RawTrigger.h:38
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
float cosnumi
Cosine of angle of track wrt to NuMI beam axis.
Definition: CanaNt.h:37
int goodspills
Definition: POTSum.h:31
void reconfigure(const fhicl::ParameterSet &p)
Definition: Cana_module.cc:526
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
int Apply(art::PtrVector< rawdata::RawDigit > &rd)
float tprev
time to previous time window (usec)
Definition: CanaNt.h:36
float dedgeb
distance to detector edge going south along track (cm)
Definition: CanaNt.h:53
float nuy
MC neutrino interaction y location (cm)
Definition: CanaNt.h:56
float nux
MC neutrino interaction x location (cm)
Definition: CanaNt.h:55
int run
run number
Definition: CanaNt.h:16
TH2F * fNhitTotalVsEvent
Definition: Cana_module.cc:70
double Vz(const int i=0) const
Definition: MCParticle.h:222
Timestamp time() const
Definition: Event.h:61
TTree * fNtTree
Definition: Cana_module.cc:58
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
int cell2y
last cell in y view
Definition: CanaNt.h:31
void beginJob()
Definition: Cana_module.cc:100
float leppx
MC out going lepton momentum (GeV)
Definition: CanaNt.h:61
int cell1y
1st cell in y view
Definition: CanaNt.h:30
virtual ~Cana()
Definition: Cana_module.cc:153
int totspills
Definition: POTSum.h:30
int fApplyBadChan
Apply bad channel masks?
Definition: Cana_module.cc:74
int plane2x
last plane in x view
Definition: CanaNt.h:25
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
RunNumber_t run() const
Definition: Event.h:77
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
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
double totpot
normalized by 10^12 POT
Definition: POTSum.h:27
TH1F * fNhitTimeWin
Definition: Cana_module.cc:67
Definition: fwd.h:28
std::string fMCLabel
Module that produced the mc truth.
Definition: Cana_module.cc:76
int plane1x
1st plane in x view
Definition: CanaNt.h:24
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
EventID id() const
Definition: Event.h:56
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
int plane2y
last plane in y view
Definition: CanaNt.h:27
double fX2
End x position (cm)
Definition: RLFit.h:108
float cosbnb
Cosine of angle of track wrt to BNB beam axis.
Definition: CanaNt.h:38
static constexpr Double_t sr
Definition: Munits.h:164
std::string fPotLabel
Module that produced the POTSum object.
Definition: Cana_module.cc:78
float coscosmic
Cosine of zenith angle (point back wrt. vertical)
Definition: CanaNt.h:41
int nhity
Number of y view hits on track in peak time window.
Definition: CanaNt.h:23