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
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  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
91  fTrackType (pset.get< std::string >("TrackType") ),
92  fDoCleanup (0),
93  fDoTrackAngleCut (0),
94  fDoXYContainment (0),
95  fDoTrackContain (0),
96  fUseStrongCut (0),
97  fDoInTimeCut (0),
99  fNCellLow (0),
101  fMinTCut (0),
102  fMaxTCut (0),
103  fUseWeakSel (0),
104  fSelectContained (0),
106  fSelectEntering (0),
107  fSelectExiting (0),
108  fPSetNDOS (pset.get< fhicl::ParameterSet >("ndos")),
109  fPSetND (pset.get< fhicl::ParameterSet >("nd") ),
110  fPSetFD (pset.get< fhicl::ParameterSet >("fd") )
111  {
112  produces< rb::FilterList<rb::Cluster> >();
113  }
114 
115  //----------------------------------------------------------------------
117  {
118  }
119 
120  //----------------------------------------------------------------------
122  {
124 
125  fhicl::ParameterSet pset;
126 
127  switch(geom->DetId()){
128  case novadaq::cnv::kNDOS:
129  pset = fPSetNDOS;
130  break;
132  pset = fPSetND;
133  break;
135  pset = fPSetFD;
136  break;
137  default:
138  assert(0 && "Unknown detector");
139  }
140 
141  fDoCleanup = pset.get< bool >("DoCleanup");
142  fDoTrackAngleCut = pset.get< bool >("DoTrackAngleCut");
143  fDoXYContainment = pset.get< bool >("DoXYContainment");
144  fDoTrackContain = pset.get< bool >("DoTrackContain");
145  fUseStrongCut = pset.get< bool >("UseStrongCut");
146  fDoInTimeCut = pset.get< bool >("DoInTimeCut");
147  fDoUnContainedHitCut = pset.get< bool >("DoUnContainedHitCut");
148  fNCellLow = pset.get< unsigned int >("NCellLow");
149  fMaxUncontainedHits = pset.get< unsigned int >("MaxUncontainedHits");
150  fMinTCut = pset.get< double >("MinTCut");
151  fMaxTCut = pset.get< double >("MaxTCut");
152  fUseWeakSel = pset.get< bool >("UseWeakSel");
153  fSelectContained = pset.get< bool >("SelectContained");
154  fSelectThroughGoing = pset.get< bool >("SelectThroughGoing");
155  fSelectEntering = pset.get< bool >("SelectEntering");
156  fSelectExiting = pset.get< bool >("SelectExiting");
157 
158  fMinXFid = pset.get< float >("MinXFid");
159  fMaxXFid = pset.get< float >("MaxXFid");
160  fMinYFid = pset.get< float >("MinYFid");
161  fMaxYFid = pset.get< float >("MaxYFid");
162  fMinZ = pset.get< float >("MinZ");
163  fZCut = pset.get< float >("ZCut");
164 
165  return true;
166  }
167 
168  //----------------------------------------------------------------------
170  {
171 
172  std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
173 
174  bool selected = false;
175  //Grab event's slices
177  evt.getByLabel(fSlicerLabel, slices);
178 
179  // Object holding tracks associated with slice
180  art::FindManyP<rb::Track> fmpTrack(slices, evt, fTrackType);
181 
182  const int sliceMax = slices->size();
183 
184  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
185  const rb::Cluster& slice = (*slices)[sliceIdx];
186 
187  if(slice.IsNoise())
188  {
189  filtcol->Add(slices, sliceIdx);
190  continue;
191  }
192 
193  double ratio = 1, diry = 0, angle = 1;
194  double minx = -9999, miny = -9999, minz = -9999;
195  double maxx = 9999, maxy = 9999, maxz = 9999;
196 
197  TVector3 start,stop;
198 
199  double meant = 0;
200 
201  unsigned int nXplanes = 0, nYplanes = 0;
202  bool is3d = false;
203 
204  unsigned int stot[6] = {0,0,0,0,0,0};
205 
206  const double nCells = slice.NCell(); //Number of cells in slice
207  unsigned int PXhi = 0;
208  unsigned int PXlo = UINT_MAX;
209  unsigned int PYhi = 0;
210  unsigned int PYlo = UINT_MAX;
211  //loop over all hits
212 
216 
217  double fMinZFid = 0.;
218  double zfront1 = 0.,zback1 = 0.,zfront2 = 0.,zback2 = 0.;
219 
220  if(geom->DetId() == novadaq::cnv::kFARDET){
221  int ndets = livegeom->InstrumentedDetEnds(zfront1,zback1,zfront2,zback2);
222  // take the most downstream zfront and zback
223  if(ndets > 0){
224  fMinZFid = zfront1 + fMinZ;
225  }
226  }
227  else fMinZFid = fMinZ;
228 
229 
230  fMaxZFid = zback1 - fZCut;
231 
232  // std::cout << zback1 << " " << zfront1 << "\n";
233  for(int view = geo::kX; view <= geo::kY; ++view) {
234  const geo::View_t geoview = geo::View_t(view);
235 
236  for(unsigned int j = 0; j < slice.NCell(geoview); ++j) {
237 
238  const rb::CellHit* h = slice.Cell(geoview,j).get();
239 
240  // Get Hit position
241  double xyz[3];
242  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
243  const double x = xyz[0];
244  const double y = xyz[1];
245  const double z = xyz[2];
246 
247  if(h->View() == geo::kX) {
248  if(h->Plane() > PXhi) PXhi = h->Plane();
249  if(h->Plane() < PXlo) PXlo = h->Plane();
250  }
251  if(h->View() == geo::kY) {
252  if(h->Plane() > PYhi) PYhi = h->Plane();
253  if(h->Plane() < PYlo) PYlo = h->Plane();
254  }
255 
256 
257  // Total up number of cells > Max (X/Y/Z) and < Min (X/Y/Z)
258  //TOP
259  if(y > fMaxYFid) {
260  stot[0]++;
261  // to avoid double-counting, stop counting for subsequent boundaries
262  continue;
263  }
264  //BOTTOM
265  if(y < fMinYFid) {
266  stot[1]++;
267  continue;
268  }
269  //LEFT
270  if(x > fMaxXFid) {
271  stot[2]++;
272  continue;
273  }
274  //RIGHT
275  if(x < fMinXFid) {
276  stot[3]++;
277  continue;
278  }
279  //FRONT
280  if(z < fMinZFid) {
281  stot[4]++;
282  continue;
283  }
284  //BACK
285  if(z > fMaxZFid) {
286  stot[5]++;
287  continue;
288  }
289 
290  }
291  }
292 
293 
294  nXplanes = PXhi - PXlo+1;
295  nYplanes = PYhi - PYlo+1;
296 
297  if(fmpTrack.isValid()){
298  std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(sliceIdx);
299  if(tracks.size()>0) {
300  const art::Ptr<rb::Track> track = tracks[0];
301 
302  meant = slice.MeanTNS();
303 
304 
305 
306  minx = slice.MinXYZ().X();
307  miny = slice.MinXYZ().Y();
308  minz = slice.MinXYZ().Z();
309 
310  maxx = slice.MaxXYZ().X();
311  maxy = slice.MaxXYZ().Y();
312  maxz = slice.MaxXYZ().Z();
313 
314  start = track->Start();
315  stop = track->Stop();
316 
317  diry = track->Dir().Y();
318  angle = fCosRejFxs->getAngle(track->Dir());
319  ratio = (track->NCell()*1.0) / (1.0*nCells);
320 
321 
322  is3d = track->Is3D();
323 
324 
325  }
326  }
327 
328  // Selected unless we fail any cuts below
329  bool sel = true;
330 
331  if(fDoInTimeCut){
332  if(meant < fMinTCut || meant > fMaxTCut){
333  sel = false;
334  }
335  }
336  if(sel == true) std::cout<< "True: Timing";
337  //----------------Actually start the cut selections here----------------
338  //Do the NCell cut
339  if(fDoCleanup){
340  if(nCells <= fNCellLow || !is3d) {
341  sel = false;
342  }
343  }
344  if(sel == true) std::cout<< "True: Cleanup";
345  //Do the XY slice.box containment cut + the track start/stop containment
346  //Change this to use the fwddist/bckdist at a later date
347  if(fDoXYContainment){
348  if(minx < fMinXFid || miny < fMinYFid ||
349  maxx > fMaxXFid || maxy > fMaxYFid)
350  sel = false;
351  }
352  if(sel == true) std::cout<< "True: XYContainment";
353 
355  const unsigned int smaxelement = *std::max_element(stot,stot+6);
356  if(smaxelement > fMaxUncontainedHits)
357  sel = false;
358  }
359  if(sel == true) std::cout<< "True: BBC";
360 
361  if(fDoTrackContain){
362  if(start.X() < fMinXFid || start.Y() < fMinYFid ||
363  stop.X() > fMaxXFid || stop.Y() > fMaxYFid ||
364  stop.X() < fMinXFid || stop.Y() < fMinYFid ||
365  start.X() > fMaxXFid || start.Y() > fMaxYFid) {
366  sel = false;
367  }
368  }
369  if(sel == true) std::cout<< "True: Contain";
370 
371  //Do the Track Angle cut:
372  if (fDoTrackAngleCut){
373  if(fmpTrack.isValid()) {
374 
375  double cosmicvar = (angle*angle) - (diry*diry);
376 
377  if(fUseStrongCut == false) { // weak cut
378  if(ratio > 0.88) {
379  if(nCells > 50 && nCells <= 135 &&
380  cosmicvar < ((1.1/85)*nCells-1.65)) {
381  sel = false;
382  }
383 
384  if(nCells > 135 && nCells <= 250 &&
385  cosmicvar < ((.6/115)*nCells-.6)) {
386  sel = false;
387  }
388 
389  if(nCells > 250 && nCells <= 500 &&
390  cosmicvar < ((.1/250)*nCells+0.6)) {
391  sel = false;
392  }
393 
394  if(nCells > 500 && cosmicvar > 0.8) {
395  sel = false;
396  }
397  }
398  }
399 
400  if(fUseStrongCut == true) { // strong cut
401  if(nCells > 30 && nCells <= 135 && ratio > 0.8 &&
402  cosmicvar < ((1/100)*nCells-1.0)) {
403  sel = false;
404  }
405 
406  if(nCells > 135 && nCells <= 250 && ratio > 0.8 &&
407  cosmicvar < ((.45/100)*nCells-.2575)) {
408  sel = false;
409  }
410 
411  if(nCells > 250 && nCells <= 500 && ratio > 0.83 &&
412  cosmicvar < ((.1/200)*nCells+0.6825)) {
413  sel = false;
414  }
415 
416  if(nCells > 500 && ratio > 0.8 &&
417  cosmicvar > 0.9) {
418  sel = false;
419  }
420  }
421  }
422  }
423  if(sel == true) std::cout<< "True: AngVar";
424 
425  bool iscont = false, isexit = false;
426  bool isenter = false, isthrough = false;
427 
428  if(fUseWeakSel == true){
429  if(fSelectContained){
430  if(minz > fMinZFid && maxz < fMaxZFid) iscont = true;
431  }
433  if(minz < fMinZFid && maxz > fMaxZFid) isthrough = true;
434  }
435  if(fSelectEntering){
436  if(minz < fMinZFid && maxz < fMaxZFid) isexit = true;
437  }
438  if(fSelectExiting){
439  if(minz > fMinZFid && maxz > fMaxZFid) isenter = true;
440  }
441 
442  if(!(iscont || isexit || isenter || isthrough)) {
443  sel = false;
444  }
445  }
446  /*
447  std::cout << "MinZ: " << fMinZFid
448  << " MaxZ: " << fMaxZFid
449  << " nXplanes: " << nXplanes
450  << " nYplanes: " << nYplanes
451  << " angle: " << angle
452  << "\n";
453  */
454  if(fUseWeakSel == false){
455  if(fSelectContained){
456  // std::cout << "Contained - start: " << start.Z()
457  // << " stop: " << stop.Z()
458  // << "\n";
459  if(((start.Z() > fMinZFid && stop.Z() < fMaxZFid) ||
460  (stop.Z() > fMinZFid && start.Z() < fMaxZFid)) &&
461  nXplanes > 4 && nYplanes > 4 &&
462  (angle > 0.5 || angle < -0.5)) iscont = true;
463  }
464 
466  if((start.Z() < fMinZFid && stop.Z() > fMaxZFid) ||
467  (stop.Z() < fMinZFid && start.Z() > fMaxZFid))
468  isthrough = true;
469  }
470 
471  if(fSelectEntering){
472  if(((start.Z() < fMinZFid && stop.Z() < fMaxZFid) ||
473  (stop.Z() < fMinZFid && start.Z() < fMaxZFid)) &&
474  nXplanes > 4 && nYplanes > 4 &&
475  (angle > 0.95 || angle < -0.95)) isenter = true;
476  }
477  if(fSelectExiting){
478  if(sel==true)std::cout<< "True: Exiting check";
479  if(((start.Z() > fMinZFid && stop.Z() > start.Z()) ||
480  (stop.Z() > fMinZFid && start.Z() > stop.Z())) &&
481  nXplanes > 4 && nYplanes > 4 &&
482  (angle > 0.75 || angle < -0.75)) {
483  std::cout << "Exiting - start: " << start.Z()
484  << " stop: " << stop.Z()
485  << " angle " << angle << "\n";
486  isexit = true;
487  }
488  }
489  std::cout << " IS EXIT " << isexit << " sel: " << sel << std::endl;
490  if(!isexit) sel = false;
491  /*
492  if(sel ==true){
493  std::cout << "isexit: " << isexit << "\n";
494  if(!(iscont || isexit || isenter || isthrough)) {
495  sel = false;
496  }
497  }
498  */
499  }
500 
501  if(sel == true){
502  std::cout << "Event: " << evt.id().event() << ", subrun: " << evt.subRun() << "\n";
503  selected = sel;
504  }
505 
506  if(sel == false){
507  filtcol->Add(slices, sliceIdx);
508  }
509  } // end for sliceIdx
510 
511  evt.put(std::move(filtcol));
512 
513  return selected;
514  }
515 
517 }//end namespace comi
cosrej::CosRejFxs * fCosRejFxs
A simple list of products that have been marked "filtered out".
Definition: FilterList.h:74
SubRunNumber_t subRun() const
Definition: Event.h:72
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
TH1 * ratio(TH1 *h1, TH1 *h2)
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
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:31
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
Far Detector at Ash River, MN.
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.
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T const * get() const
Definition: Ptr.h:321
fhicl::ParameterSet fPSetNDOS
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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
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
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
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71