NumiFiltering_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \file NumiFiltering_module.cc
3 // \brief Module to filter events for neutrino search.
4 // Output files used for neutrino event scanning in Far Detector
5 // Select Contained/Exiting/Entering tracks after cuts.
6 // \author Gavin Davies - gsdavies@iastate.edu
7 ///////////////////////////////////////////////////////////////////////////
8 #include <algorithm>
9 #include <cassert>
10 #include <string>
11 
12 //ROOT includes
13 #include "TMath.h"
14 #include "TVector3.h"
15 
16 //Novasoft includes
17 #include "CosRej/CosRejFxs.h"
18 #include "Geometry/Geometry.h"
19 #include "Geometry/LiveGeometry.h"
21 #include "MCCheater/BackTracker.h"
22 #include "NovaDAQConventions/DAQConventions.h"
23 #include "RecoBase/Cluster.h"
24 #include "RecoBase/FilterList.h"
25 #include "RecoBase/Track.h"
26 
27 // Framework includes
32 #include "art_root_io/TFileService.h"
34 #include "fhiclcpp/ParameterSet.h"
35 
36 
37 namespace comi
38 {
40  {
41  public:
42  explicit NumiFiltering(fhicl::ParameterSet const &pset);
44  bool beginRun(art::Run& run);
45  bool filter(art::Event& evt);
46 
47  protected:
48 
51  bool fDoCleanup;
58  unsigned int fNCellLow;
59  unsigned int fMaxUncontainedHits;
60  double fMinTCut;
61  double fMaxTCut;
63 
68 
69  float fMinXFid;
70  float fMinYFid;
71  float fMinZ;
72  float fMaxXFid;
73  float fMaxYFid;
74  float fMaxZFid;
75  float fZCut;
76 
77 
79  private:
80 
82  };
83 
84 }//End namespace comi
85 
86 namespace comi
87 {
88  //----------------------------------------------------------------------
90  EDFilter(pset),
91  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
92  fTrackType (pset.get< std::string >("TrackType") ),
93  fDoCleanup (0),
94  fDoTrackAngleCut (0),
95  fDoXYContainment (0),
96  fDoTrackContain (0),
97  fUseStrongCut (0),
98  fDoInTimeCut (0),
100  fNCellLow (0),
102  fMinTCut (0),
103  fMaxTCut (0),
104  fUseWeakSel (0),
105  fSelectContained (0),
107  fSelectEntering (0),
108  fSelectExiting (0),
109  fPSetNDOS (pset.get< fhicl::ParameterSet >("ndos")),
110  fPSetND (pset.get< fhicl::ParameterSet >("nd") ),
111  fPSetFD (pset.get< fhicl::ParameterSet >("fd") )
112  {
113  produces< rb::FilterList<rb::Cluster> >();
114  }
115 
116  //----------------------------------------------------------------------
118  {
119  }
120 
121  //----------------------------------------------------------------------
123  {
125 
126  fhicl::ParameterSet pset;
127 
128  switch(geom->DetId()){
129  case novadaq::cnv::kNDOS:
130  pset = fPSetNDOS;
131  break;
133  pset = fPSetND;
134  break;
136  pset = fPSetFD;
137  break;
138  default:
139  assert(0 && "Unknown detector");
140  }
141 
142  fDoCleanup = pset.get< bool >("DoCleanup");
143  fDoTrackAngleCut = pset.get< bool >("DoTrackAngleCut");
144  fDoXYContainment = pset.get< bool >("DoXYContainment");
145  fDoTrackContain = pset.get< bool >("DoTrackContain");
146  fUseStrongCut = pset.get< bool >("UseStrongCut");
147  fDoInTimeCut = pset.get< bool >("DoInTimeCut");
148  fDoUnContainedHitCut = pset.get< bool >("DoUnContainedHitCut");
149  fNCellLow = pset.get< unsigned int >("NCellLow");
150  fMaxUncontainedHits = pset.get< unsigned int >("MaxUncontainedHits");
151  fMinTCut = pset.get< double >("MinTCut");
152  fMaxTCut = pset.get< double >("MaxTCut");
153  fUseWeakSel = pset.get< bool >("UseWeakSel");
154  fSelectContained = pset.get< bool >("SelectContained");
155  fSelectThroughGoing = pset.get< bool >("SelectThroughGoing");
156  fSelectEntering = pset.get< bool >("SelectEntering");
157  fSelectExiting = pset.get< bool >("SelectExiting");
158 
159  fMinXFid = pset.get< float >("MinXFid");
160  fMaxXFid = pset.get< float >("MaxXFid");
161  fMinYFid = pset.get< float >("MinYFid");
162  fMaxYFid = pset.get< float >("MaxYFid");
163  fMinZ = pset.get< float >("MinZ");
164  fZCut = pset.get< float >("ZCut");
165 
166  return true;
167  }
168 
169  //----------------------------------------------------------------------
171  {
172 
173  std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
174 
175  bool selected = false;
176  //Grab event's slices
178  evt.getByLabel(fSlicerLabel, slices);
179 
180  // Object holding tracks associated with slice
181  art::FindManyP<rb::Track> fmpTrack(slices, evt, fTrackType);
182 
183  const int sliceMax = slices->size();
184 
185  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
186  const rb::Cluster& slice = (*slices)[sliceIdx];
187 
188  if(slice.IsNoise())
189  {
190  filtcol->Add(slices, sliceIdx);
191  continue;
192  }
193 
194  double ratio = 1, diry = 0, angle = 1;
195  double minx = -9999, miny = -9999, minz = -9999;
196  double maxx = 9999, maxy = 9999, maxz = 9999;
197 
198  TVector3 start,stop;
199 
200  double meant = 0;
201 
202  unsigned int nXplanes = 0, nYplanes = 0;
203  bool is3d = false;
204 
205  unsigned int stot[6] = {0,0,0,0,0,0};
206 
207  const double nCells = slice.NCell(); //Number of cells in slice
208  unsigned int PXhi = 0;
209  unsigned int PXlo = UINT_MAX;
210  unsigned int PYhi = 0;
211  unsigned int PYlo = UINT_MAX;
212  //loop over all hits
213 
217 
218  double fMinZFid = 0.;
219  double zfront1 = 0.,zback1 = 0.,zfront2 = 0.,zback2 = 0.;
220 
221  if(geom->DetId() == novadaq::cnv::kFARDET){
222  int ndets = livegeom->InstrumentedDetEnds(zfront1,zback1,zfront2,zback2);
223  // take the most downstream zfront and zback
224  if(ndets > 0){
225  fMinZFid = zfront1 + fMinZ;
226  }
227  }
228  else fMinZFid = fMinZ;
229 
230 
231  fMaxZFid = zback1 - fZCut;
232 
233  // std::cout << zback1 << " " << zfront1 << "\n";
234  for(int view = geo::kX; view <= geo::kY; ++view) {
235  const geo::View_t geoview = geo::View_t(view);
236 
237  for(unsigned int j = 0; j < slice.NCell(geoview); ++j) {
238 
239  const rb::CellHit* h = slice.Cell(geoview,j).get();
240 
241  // Get Hit position
242  double xyz[3];
243  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
244  const double x = xyz[0];
245  const double y = xyz[1];
246  const double z = xyz[2];
247 
248  if(h->View() == geo::kX) {
249  if(h->Plane() > PXhi) PXhi = h->Plane();
250  if(h->Plane() < PXlo) PXlo = h->Plane();
251  }
252  if(h->View() == geo::kY) {
253  if(h->Plane() > PYhi) PYhi = h->Plane();
254  if(h->Plane() < PYlo) PYlo = h->Plane();
255  }
256 
257 
258  // Total up number of cells > Max (X/Y/Z) and < Min (X/Y/Z)
259  //TOP
260  if(y > fMaxYFid) {
261  stot[0]++;
262  // to avoid double-counting, stop counting for subsequent boundaries
263  continue;
264  }
265  //BOTTOM
266  if(y < fMinYFid) {
267  stot[1]++;
268  continue;
269  }
270  //LEFT
271  if(x > fMaxXFid) {
272  stot[2]++;
273  continue;
274  }
275  //RIGHT
276  if(x < fMinXFid) {
277  stot[3]++;
278  continue;
279  }
280  //FRONT
281  if(z < fMinZFid) {
282  stot[4]++;
283  continue;
284  }
285  //BACK
286  if(z > fMaxZFid) {
287  stot[5]++;
288  continue;
289  }
290 
291  }
292  }
293 
294 
295  nXplanes = PXhi - PXlo+1;
296  nYplanes = PYhi - PYlo+1;
297 
298  if(fmpTrack.isValid()){
299  std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(sliceIdx);
300  if(tracks.size()>0) {
301  const art::Ptr<rb::Track> track = tracks[0];
302 
303  meant = slice.MeanTNS();
304 
305 
306 
307  minx = slice.MinXYZ().X();
308  miny = slice.MinXYZ().Y();
309  minz = slice.MinXYZ().Z();
310 
311  maxx = slice.MaxXYZ().X();
312  maxy = slice.MaxXYZ().Y();
313  maxz = slice.MaxXYZ().Z();
314 
315  start = track->Start();
316  stop = track->Stop();
317 
318  diry = track->Dir().Y();
319  angle = fCosRejFxs->getAngle(track->Dir());
320  ratio = (track->NCell()*1.0) / (1.0*nCells);
321 
322 
323  is3d = track->Is3D();
324 
325 
326  }
327  }
328 
329  // Selected unless we fail any cuts below
330  bool sel = true;
331 
332  if(fDoInTimeCut){
333  if(meant < fMinTCut || meant > fMaxTCut){
334  sel = false;
335  }
336  }
337  if(sel == true) std::cout<< "True: Timing";
338  //----------------Actually start the cut selections here----------------
339  //Do the NCell cut
340  if(fDoCleanup){
341  if(nCells <= fNCellLow || !is3d) {
342  sel = false;
343  }
344  }
345  if(sel == true) std::cout<< "True: Cleanup";
346  //Do the XY slice.box containment cut + the track start/stop containment
347  //Change this to use the fwddist/bckdist at a later date
348  if(fDoXYContainment){
349  if(minx < fMinXFid || miny < fMinYFid ||
350  maxx > fMaxXFid || maxy > fMaxYFid)
351  sel = false;
352  }
353  if(sel == true) std::cout<< "True: XYContainment";
354 
356  const unsigned int smaxelement = *std::max_element(stot,stot+6);
357  if(smaxelement > fMaxUncontainedHits)
358  sel = false;
359  }
360  if(sel == true) std::cout<< "True: BBC";
361 
362  if(fDoTrackContain){
363  if(start.X() < fMinXFid || start.Y() < fMinYFid ||
364  stop.X() > fMaxXFid || stop.Y() > fMaxYFid ||
365  stop.X() < fMinXFid || stop.Y() < fMinYFid ||
366  start.X() > fMaxXFid || start.Y() > fMaxYFid) {
367  sel = false;
368  }
369  }
370  if(sel == true) std::cout<< "True: Contain";
371 
372  //Do the Track Angle cut:
373  if (fDoTrackAngleCut){
374  if(fmpTrack.isValid()) {
375 
376  double cosmicvar = (angle*angle) - (diry*diry);
377 
378  if(fUseStrongCut == false) { // weak cut
379  if(ratio > 0.88) {
380  if(nCells > 50 && nCells <= 135 &&
381  cosmicvar < ((1.1/85)*nCells-1.65)) {
382  sel = false;
383  }
384 
385  if(nCells > 135 && nCells <= 250 &&
386  cosmicvar < ((.6/115)*nCells-.6)) {
387  sel = false;
388  }
389 
390  if(nCells > 250 && nCells <= 500 &&
391  cosmicvar < ((.1/250)*nCells+0.6)) {
392  sel = false;
393  }
394 
395  if(nCells > 500 && cosmicvar > 0.8) {
396  sel = false;
397  }
398  }
399  }
400 
401  if(fUseStrongCut == true) { // strong cut
402  if(nCells > 30 && nCells <= 135 && ratio > 0.8 &&
403  cosmicvar < ((1/100)*nCells-1.0)) {
404  sel = false;
405  }
406 
407  if(nCells > 135 && nCells <= 250 && ratio > 0.8 &&
408  cosmicvar < ((.45/100)*nCells-.2575)) {
409  sel = false;
410  }
411 
412  if(nCells > 250 && nCells <= 500 && ratio > 0.83 &&
413  cosmicvar < ((.1/200)*nCells+0.6825)) {
414  sel = false;
415  }
416 
417  if(nCells > 500 && ratio > 0.8 &&
418  cosmicvar > 0.9) {
419  sel = false;
420  }
421  }
422  }
423  }
424  if(sel == true) std::cout<< "True: AngVar";
425 
426  bool iscont = false, isexit = false;
427  bool isenter = false, isthrough = false;
428 
429  if(fUseWeakSel == true){
430  if(fSelectContained){
431  if(minz > fMinZFid && maxz < fMaxZFid) iscont = true;
432  }
434  if(minz < fMinZFid && maxz > fMaxZFid) isthrough = true;
435  }
436  if(fSelectEntering){
437  if(minz < fMinZFid && maxz < fMaxZFid) isexit = true;
438  }
439  if(fSelectExiting){
440  if(minz > fMinZFid && maxz > fMaxZFid) isenter = true;
441  }
442 
443  if(!(iscont || isexit || isenter || isthrough)) {
444  sel = false;
445  }
446  }
447  /*
448  std::cout << "MinZ: " << fMinZFid
449  << " MaxZ: " << fMaxZFid
450  << " nXplanes: " << nXplanes
451  << " nYplanes: " << nYplanes
452  << " angle: " << angle
453  << "\n";
454  */
455  if(fUseWeakSel == false){
456  if(fSelectContained){
457  // std::cout << "Contained - start: " << start.Z()
458  // << " stop: " << stop.Z()
459  // << "\n";
460  if(((start.Z() > fMinZFid && stop.Z() < fMaxZFid) ||
461  (stop.Z() > fMinZFid && start.Z() < fMaxZFid)) &&
462  nXplanes > 4 && nYplanes > 4 &&
463  (angle > 0.5 || angle < -0.5)) iscont = true;
464  }
465 
467  if((start.Z() < fMinZFid && stop.Z() > fMaxZFid) ||
468  (stop.Z() < fMinZFid && start.Z() > fMaxZFid))
469  isthrough = true;
470  }
471 
472  if(fSelectEntering){
473  if(((start.Z() < fMinZFid && stop.Z() < fMaxZFid) ||
474  (stop.Z() < fMinZFid && start.Z() < fMaxZFid)) &&
475  nXplanes > 4 && nYplanes > 4 &&
476  (angle > 0.95 || angle < -0.95)) isenter = true;
477  }
478  if(fSelectExiting){
479  if(sel==true)std::cout<< "True: Exiting check";
480  if(((start.Z() > fMinZFid && stop.Z() > start.Z()) ||
481  (stop.Z() > fMinZFid && start.Z() > stop.Z())) &&
482  nXplanes > 4 && nYplanes > 4 &&
483  (angle > 0.75 || angle < -0.75)) {
484  std::cout << "Exiting - start: " << start.Z()
485  << " stop: " << stop.Z()
486  << " angle " << angle << "\n";
487  isexit = true;
488  }
489  }
490  std::cout << " IS EXIT " << isexit << " sel: " << sel << std::endl;
491  if(!isexit) sel = false;
492  /*
493  if(sel ==true){
494  std::cout << "isexit: " << isexit << "\n";
495  if(!(iscont || isexit || isenter || isthrough)) {
496  sel = false;
497  }
498  }
499  */
500  }
501 
502  if(sel == true){
503  std::cout << "Event: " << evt.id().event() << ", subrun: " << evt.subRun() << "\n";
504  selected = sel;
505  }
506 
507  if(sel == false){
508  filtcol->Add(slices, sliceIdx);
509  }
510  } // end for sliceIdx
511 
512  evt.put(std::move(filtcol));
513 
514  return selected;
515  }
516 
518 }//end namespace comi
cosrej::CosRejFxs * fCosRejFxs
A simple list of products that have been marked "filtered out".
Definition: FilterList.h:74
back track the reconstruction to the simulation
Double_t angle
Definition: plot.C:86
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
double maxy
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
const PlaneGeo * Plane(unsigned int i) const
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
fhicl::ParameterSet fPSetND
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:21
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
unsigned short Cell() const
Definition: CellHit.h:40
Far Detector at Ash River, MN.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
Prototype Near Detector on the surface at FNAL.
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 getAngle(TVector3 trackdir)
Definition: CosRejFxs.cxx:33
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
Near Detector in the NuMI cavern.
SubRunNumber_t subRun() const
const double j
Definition: BetheBloch.cxx:29
z
Definition: test.py:28
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
fhicl::ParameterSet fPSetNDOS
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
bool filter(art::Event &evt)
fhicl::ParameterSet fPSetFD
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.h:22
assert(nhit_max >=nhit_nbins)
bool beginRun(art::Run &run)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
T const * get() const
Definition: Ptr.h:149
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
NumiFiltering(fhicl::ParameterSet const &pset)
int InstrumentedDetEnds(double &frontDwnstrm, double &backDwnstrm, double &frontUpstrm, double &backUpstrm)
give all detector ends information; most general
unsigned int fMaxUncontainedHits
EventID id() const
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
enum BeamMode string