Filter_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Filter.cxx
3 // \brief Filter events based on simple expressions
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "NovaDAQConventions/DAQConventions.h"
8 #include "Geometry/Geometry.h"
10 #include "RecoBase/CellHit.h"
11 #include "Simulation/FLSHitList.h"
12 #include "RecoBase/Cluster.h"
13 #include "RecoBase/Track.h"
15 #include "Utilities/func/MathUtil.h"
16 #include "MCCheater/BackTracker.h"
17 // Framework includes
23 
24 #include "TFormula.h"
25 
26 namespace filter
27 {
28  enum kPars{
29  kNEvts, ///< Number of true neutrino interactions
30  kPdg, ///< PDG code of the neutrino
31  kCcNc, ///< Value of Nu().CCNC()
32  kIact, ///< Value of Nu().InteractionType()
33  kEnu, ///< True neutrino energy (GeV)
34  kX, ///< True vertex x position (cm)
35  kY, ///< True vertex y position (cm)
36  kZ, ///< True vertex z position (cm)
37  kEdge, ///< Distance to detector edge (cm)
38  kInsideFiducialVolume, ///< Is event inside fiducial volume? (fiducial volume is defined in Geometry.fcl)
39  kNhit, ///< Total number of hits in the event
40  kTotAdc, ///< Total ADC in the event
41  kNSlice, ///< Total number of slices in the event (including noise slice)
42  kSliceNXHit, ///< Number of hits in X view of the biggest physics slice
43  kSliceNYHit, ///< Number of hits in X view of the biggest physics slice
44  kVertexInFiducial, ///< Is track vertex (start position) inside the fiducial volume defined in filter
47  };
48 
49  // This list must be kept in sync with the enum above. This sucks, but
50  // apparently C++ doesn't have designated initialization (C does).
51  const TString kParNames[kNumPars] = {"nevts", "pdg", "ccnc", "iact", "enu",
52  "x", "y", "z", "edge", "is_inside_fiducial_volume",
53  "nhit", "totadc",
54  "nslice", "slicenxhit", "slicenyhit","vertex_in_fiducial","nFLShit"};
55 
56  /// Filter events based on simple expressions
57  class Filter: public art::EDFilter
58  {
59  public:
60  explicit Filter(const fhicl::ParameterSet& pset);
61  ~Filter();
62 
63  bool filter(art::Event& evt);
64 
65  void reconfigure(const fhicl::ParameterSet& pset);
66 
67  void beginJob();
68 
69  protected:
70  void FillTruthVariables(const art::Event& evt, double* vars) const;
71  void FillDigitVariables(const art::Event& evt, double* vars) const;
72  void FillG4GenVariables(const art::Event& evt, double* vars) const;
73  void FillSliceVariables(const art::Event& evt, double* vars) const;
74  void FillTrackVariables(const art::Event& evt, double* vars) const;
75 
76  bool isVertexInFiducial(const rb::Track& trk) const;
77 
78  TFormula* fFormula;
79 
80  std::string fGenLabel; ///< Default "generator"
81  std::string fG4GenLabel; ///< Default "geantgen"
82  std::string fCellHitLabel; ///< Default "daq"
83  std::string fSlicerLabel; ///< Default "slicer"
84  std::string fTrackLabel; ///< Default "kalmantrack"
85  };
86 
87  //......................................................................
89  : fFormula(0)
90  {
91  reconfigure(pset);
92  }
93 
94  //......................................................................
96  {
97  delete fFormula;
98  }
99 
100  bool compareBySecondLength(const std::pair<int, std::string>& a,
101  const std::pair<int, std::string>& b)
102  {
103  return a.second.size() < b.second.size();
104  }
105 
106  //......................................................................
108  {
109  fGenLabel = pset.get<std::string>("GeneratorModuleLabel", "generator");
110  fG4GenLabel = pset.get<std::string>("G4GenModuleLabel", "geantgen");
111  fCellHitLabel = pset.get<std::string>("CalHitModuleLabel", "calhit");
112  fSlicerLabel = pset.get<std::string>("SlicerModuleLabel", "slicer");
113  fTrackLabel = pset.get<std::string>("TrackModuleLabel", "kalmantrack");
114 
115  TString expr = pset.get<std::string>("Expression");
116 
117  // Empty string evaluates true
118  if(expr == "") expr = "1";
119 
120  // Need to sort the parameters from longest to shortest string, doing the
121  // replacements in any other order causes everything to go wrong...
122  std::vector<std::pair<int, std::string> > sortedPars;
123  for(int n = 0; n < kNumPars; ++n){
124  sortedPars.emplace_back(n, std::string(kParNames[n]));
125  }
126  std::sort(sortedPars.rbegin(), sortedPars.rend(), compareBySecondLength);
127 
128  // Rewrite the expression with numeric indices in place of variable names
129  for(int n = 0; n < kNumPars; ++n){
130  expr.ReplaceAll(sortedPars[n].second,
131  TString::Format("[%d]", sortedPars[n].first));
132  }
133 
134  delete fFormula;
135  fFormula = new TFormula("form", expr);
136  }
137 
138  //......................................................................
140  {
141  }
142 
143  //......................................................................
145  {
146  double vars[kNumPars] = {0};
147 
148  FillTruthVariables(evt, vars);
149  FillG4GenVariables(evt, vars);
150  FillDigitVariables(evt, vars);
151  FillSliceVariables(evt, vars);
152  FillTrackVariables(evt, vars);
153 
154  return fFormula->EvalPar(0, vars);
155  }
156 
157  //......................................................................
158  void Filter::FillTruthVariables(const art::Event& evt, double* vars) const
159  {
161  if( !bt->HaveTruthInfo() ) return;
162 
164  evt.getByLabel(fGenLabel, mcTruths);
165  if(mcTruths.failedToGet()) return;
166 
167  vars[kNEvts] = mcTruths->size();
168  // Abort if no truth info (data?)
169  if(mcTruths->empty()) return;
170  // If there are two events, take the first one...
171  simb::MCTruth mcTruth = (*mcTruths)[0];
172  // Abort if not a neutrino-induced event (cosmics?)
173  if(!mcTruth.NeutrinoSet()) return;
174 
175  const simb::MCNeutrino& nu = mcTruth.GetNeutrino();
176  vars[kCcNc] = nu.CCNC();
177  vars[kIact] = nu.InteractionType();
178 
179  const simb::MCParticle& nupart = nu.Nu();
180  vars[kPdg] = nupart.PdgCode();
181  vars[kEnu] = nupart.E();
182 
183  const simb::MCParticle& lept = nu.Lepton();
184  const double x = lept.Vx();
185  const double y = lept.Vy();
186  const double z = lept.Vz();
187  vars[kX] = x; vars[kY] = y; vars[kZ] = z;
188 
190  const double dx = geom->DetHalfWidth()-fabs(x); // If outside, then dx is negative
191  const double dy = geom->DetHalfHeight()-fabs(y); // If outside, then dy is negative
192  const double dz = 0.5*geom->DetLength() - fabs(z + 0.5*geom->DetLength());// If outside, then dz is negative
193  const double kEdge_mod = std::min(std::min(fabs(dx), fabs(dy)), fabs(dz));// minimum of the absolute values
194  vars[kEdge] = (dx < 0.0 || dy <0.0 || dz<0.0) ? -kEdge_mod : kEdge_mod; // If outside, then kEdge is negative
195 
196  /// Fill whether the point is inside fiducial volume
197  vars[kInsideFiducialVolume] = geom->isInsideFiducialVolume(x, y, z) ? 1.0 : 0.0;
198  }
199 
200  //......................................................................
201  void Filter::FillDigitVariables(const art::Event& evt, double* vars) const
202  {
203  double totadc = 0;
204  // get the RawDigit list
206  evt.getByLabel(fCellHitLabel, digitcol);
207  if(digitcol.failedToGet()) return;
208 
209  const unsigned int nhit = digitcol->size();
210  for(unsigned int i = 0; i < nhit; ++i){
211  art::Ptr<rb::CellHit> ch(digitcol, i);
212  totadc += ch->ADC();
213  }
214 
215  vars[kNhit] = nhit;
216  vars[kTotAdc] = totadc;
217  }
218 
219  //......................................................................
220 
221  void Filter::FillG4GenVariables(const art::Event& evt, double* vars) const
222  {
223  // get the FLSHit list
225  evt.getByLabel(fG4GenLabel, FLScol);
226  if (FLScol.failedToGet()) return;
227 
228  const unsigned int nhit = FLScol->size();
229  vars[kNFLShit] = nhit;
230  }
231 
232  //......................................................................
233  void Filter::FillSliceVariables(const art::Event& evt, double* vars) const
234  {
236  evt.getByLabel(fSlicerLabel, slices);
237  if(slices.failedToGet()) return;
238 
239  const unsigned int sliceMax = slices->size();
240  vars[kNSlice] = sliceMax;
241 
242  // Record the properties of the biggest non-noise slice
243  unsigned int most = 0;
244  for(unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
245  const rb::Cluster& slice = (*slices)[sliceIdx];
246  if(slice.IsNoise()) continue;
247  if(slice.NCell() > most){
248  most = slice.NCell();
249  vars[kSliceNXHit] = slice.NXCell();
250  vars[kSliceNYHit] = slice.NYCell();
251  }
252  } // end for sliceIdx
253  }
254 
255  //......................................................................
256  void Filter::FillTrackVariables(const art::Event& evt, double* vars) const
257  {
259  evt.getByLabel(fTrackLabel, trk);
260  if(trk.failedToGet()) return;
261 
262  unsigned int trkMax = trk->size();
263  // Record the properties of tracks inside a defined fiducial volume
264  for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
265  const rb::Track& track = (*trk)[trkIdx];
266  //if(track.ExtentPlane() > 10) return;
267  const TVector3 dir(track.Dir());
268  const double coscosmic = TMath::Abs(-dir.Y()/util::pythag(dir.X(),dir.Y(),dir.Z()));
269  if(coscosmic < 0.7) return;
270  //if(track.TotalLength() > 70) return;
271  const bool check = isVertexInFiducial(track);
272  vars[kVertexInFiducial] = check;
273  } // end for trkIdx
274  }
275 
277  {
278  // Start and end don't mean much, pick whichever is highest in y.
279 
280  const TVector3 v1 = trk.Start();
281  const TVector3 v2 = trk.Stop();
282  // assume down-going in order to reject cosmics
283  const TVector3 start = (v1.Y() > v2.Y()) ? v1 : v2;
284  //const TVector3 start = v1;
285 
286  const double x = start.X();
287  const double y = start.Y();
288  const double z = start.Z();
289 
291 
292 
293  switch(geom->DetId()){
294  case novadaq::cnv::kNDOS:
296  if(fabs(x) > 115) return false; // Fiducial X
297  if(fabs(y) > 150) return false; // Fiducial Y
298  if(z < 50 || z > 770) return false; // Fiducial Z
299 
300  if(v1.X()==v2.X()||v1.Y()==v2.Y()) return false; // dx or dy = 0 means bad fit
301  break;
303  if(fabs(x) > 730 || fabs(y) > 730) return false;
304  if(z < 50 || z > 6300) return false;
305 
306  if(v1.X()==v2.X()||v1.Y()==v2.Y()) return false;
307  break;
308  default:
309  assert(0 && "Unknown detector");
310  }
311 
312  return true;
313  }
314 
316 
317 } // namespace
318 ////////////////////////////////////////////////////////////////////////
void FillTrackVariables(const art::Event &evt, double *vars) const
double E(const int i=0) const
Definition: MCParticle.h:232
Filter events based on simple expressions.
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
Is event inside fiducial volume? (fiducial volume is defined in Geometry.fcl)
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void reconfigure(const fhicl::ParameterSet &pset)
Number of true neutrino interactions.
Total ADC in the event.
Value of Nu().CCNC()
True neutrino energy (GeV)
Total number of hits in the event.
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Filter(const fhicl::ParameterSet &pset)
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
Total number of slices in the event (including noise slice)
nhit
Definition: demo1.py:25
std::string fSlicerLabel
Default "slicer".
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
TFormula * fFormula
Distance to detector edge (cm)
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
bool filter(art::Event &evt)
void FillSliceVariables(const art::Event &evt, double *vars) const
Track finder for cosmic rays.
std::string fTrackLabel
Default "kalmantrack".
int InteractionType() const
Definition: MCNeutrino.h:150
std::string fGenLabel
Default "generator".
Far Detector at Ash River, MN.
double dy[NP][NC]
double dx[NP][NC]
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
void FillDigitVariables(const art::Event &evt, double *vars) const
double dz[NP][NC]
bool compareBySecondLength(const std::pair< int, std::string > &a, const std::pair< int, std::string > &b)
std::string fG4GenLabel
Default "geantgen".
const double a
Prototype Near Detector on the surface at FNAL.
std::string fCellHitLabel
Default "daq".
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
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
Number of hits in X view of the biggest physics slice.
Near Detector in the NuMI cavern.
const TString kParNames[kNumPars]
True vertex y position (cm)
const std::map< std::pair< std::string, std::string >, Variable > vars
double DetHalfHeight() const
void FillG4GenVariables(const art::Event &evt, double *vars) const
Value of Nu().InteractionType()
z
Definition: test.py:28
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
True vertex x position (cm)
double DetHalfWidth() const
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
TDirectory * dir
Definition: macro.C:5
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
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
const hit & b
Definition: hits.cxx:21
double Vz(const int i=0) const
Definition: MCParticle.h:222
assert(nhit_max >=nhit_nbins)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void FillTruthVariables(const art::Event &evt, double *vars) const
bool NeutrinoSet() const
Definition: MCTruth.h:77
Event generator information.
Definition: MCTruth.h:32
T min(const caf::Proxy< T > &a, T b)
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
True vertex z position (cm)
bool isVertexInFiducial(const rb::Track &trk) const
Event generator information.
Definition: MCNeutrino.h:18
double Vy(const int i=0) const
Definition: MCParticle.h:221
PDG code of the neutrino.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?
bool failedToGet() const
Definition: Handle.h:196
Is track vertex (start position) inside the fiducial volume defined in filter.
Number of hits in X view of the biggest physics slice.