NumiFilteringAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \file NumiFilteringAna_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 "TH1D.h"
14 #include "TH1I.h"
15 #include "TMath.h"
16 #include "TVector3.h"
17 
18 //Novasoft includes
19 #include "CosRej/CosRejFxs.h"
20 #include "Geometry/Geometry.h"
22 #include "NovaDAQConventions/DAQConventions.h"
23 #include "RecoBase/Cluster.h"
24 #include "RecoBase/FilterList.h"
25 #include "RecoBase/Prong.h"
26 #include "RecoBase/Track.h"
28 #include "Utilities/AssociationUtil.h"
29 
30 // Framework includes
41 #include "fhiclcpp/ParameterSet.h"
43 
44 
45 namespace comi
46 {
48  {
49  public:
50  explicit NumiFilteringAna(fhicl::ParameterSet const &pset);
51  virtual ~NumiFilteringAna();
52 
53  void beginRun(const art::Run& run);
54  void analyze(art::Event const& evt);
55  virtual void beginSubRun(const art::SubRun& sr);
56 
57  void beginJob();
58  void endJob();
59 
60  protected:
61 
65 
66  float fMinXFid;
67  float fMinYFid;
68  float fMinZFid;
69  float fMaxXFid;
70  float fMaxYFid;
71  float fMaxZFid;
72 
73 
75  private:
76 
77  //cosrej::CosRejFxs* fCosRejFxs;
78 
79  //event summary info
80  int run;
81  int subrun;
82  int event;
83 
85 
86  TH1D *fLivetime;
87  // Histograms
89 
90  TH1D *fSMinXHitUn;
91  TH1D *fSMinYHitUn;
92  TH1D *fSMinZHitUn;
93 
94  TH1D *fSMaxXHitUn;
95  TH1D *fSMaxYHitUn;
96  TH1D *fSMaxZHitUn;
97 
98  TH1D *fSMinX;
99  TH1D *fSMinY;
100  TH1D *fSMinZ;
101 
102  TH1D *fSMaxX;
103  TH1D *fSMaxY;
104  TH1D *fSMaxZ;
105 
106  TH1D *fFMinX;
107  TH1D *fFMinY;
108  TH1D *fFMinZ;
109 
110  TH1D *fFMaxX;
111  TH1D *fFMaxY;
112  TH1D *fFMaxZ;
113 
114  TH1D *fSNXplanes;
115  TH1D *fSNYplanes;
116  TH1D *fFNXplanes;
117  TH1D *fFNYplanes;
118 
119 
120 
121  };
122 
123 }//End namespace comi
124 
125 namespace comi
126 {
127  //----------------------------------------------------------------------
129  EDAnalyzer(pset),
130  fSlicerLabel (pset.get< std::string >("SlicerLabel") ),
131  fTrackType (pset.get< std::string >("TrackType") ),
132  fProngType (pset.get< std::string >("ProngType") ),
133  fPSetNDOS (pset.get< fhicl::ParameterSet >("ndos") ),
134  fPSetND (pset.get< fhicl::ParameterSet >("nd") ),
135  fPSetFD (pset.get< fhicl::ParameterSet >("fd") )
136  {
137  }
138 
139  //----------------------------------------------------------------------
141  {
142  }
143 
144  //----------------------------------------------------------------------
146  {
148 
149  double nhitmax = 100.0;
150  int nhitbins = 100;
151  double nplanemax = 500.0;
152  int nplanebins = 500;
153  double xmax = 1200.0;
154  double ymax = 1200.0;
155  double zmax = 6500.0;
156  int xbins = 480;
157  int ybins = 480;
158  int zbins = 1300;
159 
160  fLivetime = tfs->make<TH1D>("livetime", "TotalLiveTime;;Seconds", 1, 0, 1);
161  // # Uncontained hits in x/y/z min/max for Slicer Clusters
162  fSMaxAnyHitUn = tfs->make<TH1D>("smaxanyhitun",";# Slicer Cluster Uncontained Hits Max; # of Slices",nhitbins, 0., nhitmax);
163 
164  fSMinXHitUn = tfs->make<TH1D>("sminxhitun",";# Slicer Cluster Uncontained Hits minX; # of Slices",nhitbins, 0., nhitmax);
165  fSMinYHitUn = tfs->make<TH1D>("sminyhitun",";# Slicer Cluster Uncontained Hits minY; # of Slices",nhitbins, 0., nhitmax);
166  fSMinZHitUn = tfs->make<TH1D>("sminzhitun",";# Slicer Cluster Uncontained Hits minZ; # of Slices",nhitbins, 0., nhitmax);
167 
168  fSMaxXHitUn = tfs->make<TH1D>("smaxxhitun",";# Slicer Cluster Uncontained Hits maxX; # of Slices",nhitbins, 0., nhitmax);
169  fSMaxYHitUn = tfs->make<TH1D>("smaxyhitun",";# Slicer Cluster Uncontained Hits maxY; # of Slices",nhitbins, 0., nhitmax);
170  fSMaxZHitUn = tfs->make<TH1D>("smaxzhitun",";# Slicer Cluster Uncontained Hits maxZ; # of Slices",nhitbins, 0., nhitmax);
171 
172  // Min/Max X/Y/Z position of Slicer Cluster
173  fSMinX = tfs->make<TH1D>("sminx",";Min. X Slicer Cluster; # of Slices",xbins, -xmax, xmax);
174  fSMinY = tfs->make<TH1D>("sminy",";Min. Y Slicer Cluster; # of Slices",ybins, -ymax, ymax);
175  fSMinZ = tfs->make<TH1D>("sminz",";Min. Z Slicer Cluster; # of Slices",zbins, 0., zmax);
176 
177  // Min/Max X/Y/Z position of FuzzyK Cluster
178  fFMinX = tfs->make<TH1D>("fminx",";Min. X FuzzyK Cluster; # of Slices",xbins, -xmax, xmax);
179  fFMinY = tfs->make<TH1D>("fminy",";Min. Y FuzzyK Cluster; # of Slices",ybins, -ymax, ymax);
180  fFMinZ = tfs->make<TH1D>("fminz",";Min. Z FuzzyK Cluster; # of Slices",zbins, 0., zmax);
181 
182  // Min/Max X/Y/Z position of Slicer Cluster
183  fSMaxX = tfs->make<TH1D>("smaxx",";Max. X Slicer Cluster; # of Slices",xbins, -xmax, xmax);
184  fSMaxY = tfs->make<TH1D>("smaxy",";Max. Y Slicer Cluster; # of Slices",ybins, -ymax, ymax);
185  fSMaxZ = tfs->make<TH1D>("smaxz",";Max. Z Slicer Cluster; # of Slices",zbins, 0., zmax);
186 
187  // Min/Max X/Y/Z position of FuzzyK Cluster
188  fFMaxX = tfs->make<TH1D>("fmaxx",";Max. X FuzzyK Cluster; # of Slices",xbins, -xmax, xmax);
189  fFMaxY = tfs->make<TH1D>("fmaxy",";Max. Y FuzzyK Cluster; # of Slices",ybins, -ymax, ymax);
190  fFMaxZ = tfs->make<TH1D>("fmaxz",";Max. Z FuzzyK Cluster; # of Slices",zbins, 0., zmax);
191 
192  // # X/Y planes for Slicer Clusters
193  fSNXplanes = tfs->make<TH1D>("snxplanes",";# X Planes Slicer Cluster Extends; # of Slices",nplanebins, 0., nplanemax);
194  fSNYplanes = tfs->make<TH1D>("snyplanes",";# Y Planes Slicer Cluster Extends; # of Slices",nplanebins, 0., nplanemax);
195 
196  fFNXplanes = tfs->make<TH1D>("fnxplanes",";# X Planes FuzzyK Cluster Extends; # of Slices",nplanebins, 0., nplanemax);
197  fFNYplanes = tfs->make<TH1D>("fnyplanes",";# Y Planes FuzzyK Cluster Extends; # of Slices",nplanebins, 0., nplanemax);
198  }
199 
200  //----------------------------------------------------------------------
202  {
204 
205  fhicl::ParameterSet pset;
206 
207  switch(geom->DetId()){
208  case novadaq::cnv::kNDOS:
209  pset = fPSetNDOS;
210  break;
212  pset = fPSetND;
213  break;
215  pset = fPSetFD;
216  break;
217  default:
218  assert(0 && "Unknown detector");
219  }
220 
221  fMinXFid = pset.get< float >("MinXFid");
222  fMaxXFid = pset.get< float >("MaxXFid");
223  fMinYFid = pset.get< float >("MinYFid");
224  fMaxYFid = pset.get< float >("MaxYFid");
225  fMinZFid = pset.get< float >("MinZFid");
226  fMaxZFid = pset.get< float >("MaxZFid");
227  }
228 
229  //----------------------------------------------------------------------
231  {
232 
233  run = evt.run();
234  subrun = evt.subRun();
235  event = evt.id().event();
236 
237  //Grab event's Slicer slices
239  evt.getByLabel(fSlicerLabel, slices);
240 
241  // Object holding prongs associated with cluster
242  art::FindManyP<rb::Prong> fmpProng(slices, evt, fProngType);
243 
244  // Object holding tracks associated with cluster
245  art::FindManyP<rb::Track> fmpTrack(slices, evt, fTrackType);
246 
247  const int sliceMax = slices->size();
248 
249  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
250  if(rb::IsFiltered(evt, slices, sliceIdx)) continue;
251 
252  const rb::Cluster& slice = (*slices)[sliceIdx];
253 
254  // Skip noise slices
255  if(slice.IsNoise()) continue;
256 
257  std::vector< art::Ptr<rb::Prong> > prongs;
258 
259  double longestObj = 0;
260  unsigned int longObjIdx = 0;
261 
262  // If there are associations between Prong and Cluster use them
263  if(fmpProng.isValid()){
264  prongs = fmpProng.at(sliceIdx);
265  }
266  else continue; // otherwise skip this slice
267 
268  // 3D prongs container
269  std::vector< art::Ptr<rb::Prong> > prongs3D;
270  for(unsigned int itr = 0; itr < prongs.size(); ++itr){
271  if(prongs[itr]->Is3D()) prongs3D.push_back(prongs[itr]);
272  }
273 
274  //Find the longest 3D prong in the slice:
275  for(unsigned int n = 0; n < prongs3D.size(); ++n){
276  const double L = prongs3D[n]->TotalLength();
277  if(L > longestObj) {
278  longestObj = L;
279  longObjIdx = n;
280  }
281  }
282 
283  const double nCells = slice.NCell(); //Number of cells in slice
284 
285  unsigned int PXhi = 0, PXhiFK = 0;
286  unsigned int PXlo = UINT_MAX, PXloFK = UINT_MAX;
287  unsigned int PYhi = 0, PYhiFK = 0;
288  unsigned int PYlo = UINT_MAX, PYloFK = UINT_MAX;
289 
290  bool is3d = false;
291 
292 
293  if(fmpTrack.isValid()){
294  std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(sliceIdx);
295  if(tracks.size()>0) {
296  const art::Ptr<rb::Track> track = tracks[0];
297  is3d = track->Is3D();
298  }
299  }
300  else {
301  mf::LogWarning("NumiFilteringAna") << "No Tracks. Skip this slice!";
302  continue;
303  }
304 
305  //----------------Make any cut selections here----------------
306  // Do the NCell cut
307  if(nCells <= 20 || !is3d) continue;
308  if(prongs3D.size() < 1) continue;
309 
310  double sminx = -9999., sminy = -9999., sminz = -9999.;
311  double smaxx = 9999., smaxy = 9999., smaxz = 9999.;
312 
313  double fminx = -9999., fminy = -9999., fminz = -9999.;
314  double fmaxx = 9999., fmaxy = 9999., fmaxz = 9999.;
315 
316  sminx = slice.MinXYZ().X();
317  sminy = slice.MinXYZ().Y();
318  sminz = slice.MinXYZ().Z();
319 
320  smaxx = slice.MaxXYZ().X();
321  smaxy = slice.MaxXYZ().Y();
322  smaxz = slice.MaxXYZ().Z();
323 
324  // Skip 4 diblock slices
325  if(sminz > 1702 || smaxz > 1702) continue;
326 
327 
328  fminx = prongs3D[longObjIdx]->MinXYZ().X();
329  fminy = prongs3D[longObjIdx]->MinXYZ().Y();
330  fminz = prongs3D[longObjIdx]->MinXYZ().Z();
331 
332  fmaxx = prongs3D[longObjIdx]->MaxXYZ().X();
333  fmaxy = prongs3D[longObjIdx]->MaxXYZ().Y();
334  fmaxz = prongs3D[longObjIdx]->MaxXYZ().Z();
335 
336  unsigned int nXplanes = 0, nYplanes = 0;
337 
338  unsigned int stot[6] = {0,0,0,0,0,0};
339 
340  //loop over all hits
341  for(int view = geo::kX; view <= geo::kY; ++view) {
342  const geo::View_t geoview = geo::View_t(view);
343 
344  for(unsigned int j = 0; j < slice.NCell(geoview); ++j) {
345 
346  const rb::CellHit* h = slice.Cell(geoview,j).get();
347 
349 
350  double xyz[3];
351  geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
352  const double x = xyz[0];
353  const double y = xyz[1];
354  const double z = xyz[2];
355 
356  if(h->View() == geo::kX) {
357  if(h->Plane() > PXhi) PXhi = h->Plane();
358  if(h->Plane() < PXlo) PXlo = h->Plane();
359  }
360  if(h->View() == geo::kY) {
361  if(h->Plane() > PYhi) PYhi = h->Plane();
362  if(h->Plane() < PYlo) PYlo = h->Plane();
363  }
364  // Total up number of cells > Max (X/Y/Z) and < Min (X/Y/Z)
365  //TOP
366  if(y > fMaxYFid) {
367  stot[0]++;
368  continue;
369  }
370  //BOTTOM
371  if(y < fMinYFid) {
372  stot[1]++;
373  continue;
374  }
375  //LEFT
376  if(x > fMaxXFid) {
377  stot[2]++;
378  continue;
379  }
380  //RIGHT
381  if(x < fMinXFid) {
382  stot[3]++;
383  continue;
384  }
385  //FRONT
386  if(z < fMinZFid) {
387  stot[4]++;
388  continue;
389  }
390  //BACK
391  if(z > fMaxZFid) {
392  stot[5]++;
393  continue;
394  }
395 
396  }
397  }
398 
399 
400  nXplanes = PXhi - PXlo+1;
401  nYplanes = PYhi - PYlo+1;
402 
403  unsigned int nXplanesFK = 0, nYplanesFK = 0;
404 
405  //loop over all hits
406  for(int view = geo::kX; view <= geo::kY; ++view) {
407  const geo::View_t geoview = geo::View_t(view);
408  const unsigned int ncellview = prongs3D[longObjIdx]->NCell(geoview);
409 
410  for(unsigned int j = 0; j < ncellview; ++j) {
411 
413  const rb::CellHit* h = prongs3D[longObjIdx]->Cell(geoview,j).get();
414 
415  // Find max/min plane number
416  if(h->View() == geo::kX) {
417  if(h->Plane() > PXhiFK) PXhiFK = h->Plane();
418  if(h->Plane() < PXloFK) PXloFK = h->Plane();
419  }
420  if(h->View() == geo::kY) {
421  if(h->Plane() > PYhiFK) PYhiFK = h->Plane();
422  if(h->Plane() < PYloFK) PYloFK = h->Plane();
423  }
424  }
425  }
426 
427  nXplanesFK = PXhiFK - PXloFK+1;
428  nYplanesFK = PYhiFK - PYloFK+1;
429 
430 
431  // Fill Histograms
432  fSNXplanes->Fill(nXplanes);
433  fSNYplanes->Fill(nYplanes);
434 
435  fFNXplanes->Fill(nXplanesFK);
436  fFNYplanes->Fill(nYplanesFK);
437 
438  // Min/Max postions
439  fSMinX->Fill(sminx);
440  fSMinY->Fill(sminy);
441  fSMinZ->Fill(sminz);
442 
443  fFMinX->Fill(fminx);
444  fFMinY->Fill(fminy);
445  fFMinZ->Fill(fminz);
446 
447  fSMaxX->Fill(smaxx);
448  fSMaxY->Fill(smaxy);
449  fSMaxZ->Fill(smaxz);
450 
451  fFMaxX->Fill(fmaxx);
452  fFMaxY->Fill(fmaxy);
453  fFMaxZ->Fill(fmaxz);
454 
455  // Total # of hits above Max, below Min X/Y/Z
456  fSMaxYHitUn->Fill(stot[0]);
457  fSMinYHitUn->Fill(stot[1]);
458  fSMaxXHitUn->Fill(stot[2]);
459  fSMinXHitUn->Fill(stot[3]);
460  fSMinZHitUn->Fill(stot[4]);
461  fSMaxZHitUn->Fill(stot[5]);
462 
463  //std::cout << "The largest number is: " << *std::max_element(stot,stot+6) << "\n";
464 
465  const unsigned int stotmax = *std::max_element(stot,stot+6);
466  fSMaxAnyHitUn->Fill(stotmax);
467 
468  } // end for sliceIdx
469  }
470 
471 
473  {
475 
476  sr.getByLabel("exposure", expo);
477  if(!expo.failedToGet())
478  {
479  fTotalLivetime += expo->totlivetime;
480  mf::LogInfo("CAFMaker") << "Exposure in subrun " << sr.subRun()
481  << " = " << expo->totlivetime;
482  }
483  }
484 
486  {
487  fLivetime->Fill(.5, fTotalLivetime);
488  }
489 
491 }//end namespace comi
SubRunNumber_t subRun() const
Definition: Event.h:72
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::map< std::string, double > xmax
SubRunNumber_t subRun() const
Definition: SubRun.h:44
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
Definition: CellHit.h:39
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
NumiFilteringAna(fhicl::ParameterSet const &pset)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Definition: Run.h:31
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
Double_t ymax
Definition: plot.C:25
unsigned short Cell() const
Definition: CellHit.h:40
Far Detector at Ash River, MN.
const int xbins
Definition: MakePlots.C:82
static constexpr double L
Prototype Near Detector on the surface at FNAL.
Commissioning files to look at the quality of our data.
Definition: Cana_module.cc:39
void beginRun(const art::Run &run)
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
Near Detector in the NuMI cavern.
caf::StandardRecord * sr
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
z
Definition: test.py:28
Definition: run.py:1
const int ybins
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
void analyze(art::Event const &evt)
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
T const * get() const
Definition: Ptr.h:321
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
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
virtual void beginSubRun(const art::SubRun &sr)
RunNumber_t run() const
Definition: Event.h:77
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
bool failedToGet() const
Definition: Handle.h:196
const Binning zbins
enum BeamMode string