Preselection_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \file Preselection_module.cc
3 // \brief Preselection to skip events that fail cuts for optimisation
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ///////////////////////////////////////////////////////////////////////////
6 
8 
9 #include "Geometry/Geometry.h"
12 
13 #include "RecoBase/CellHit.h"
14 #include "RecoBase/RecoHit.h"
15 #include "RecoBase/Track.h"
16 #include "RecoBase/Cluster.h"
17 #include "RecoBase/FilterList.h"
18 #include "RecoBase/PID.h"
19 
22 
23 #include "Utilities/AssociationUtil.h"
24 
28 
33 
34 namespace lem
35 {
37  {
38  public:
39  explicit Preselection(const fhicl::ParameterSet& pset);
40  virtual ~Preselection();
41 
42  virtual void produce(art::Event& evt);
43 
44  protected:
46 
47  bool fTruthMode;
54  double fRecoELow;
55  double fRecoEHigh;
56  bool fInFidCut;
57  };
58 
59  //----------------------------------------------------------------------
61  fSlicerLabel(pset.get<std::string>("SlicerLabel")),
62  fTruthMode(pset.get<bool>("TruthMode")),
63  fSkipTrackCut(pset.get<bool>("SkipTrackCut")),
64  fSkipEnergyCut(pset.get<bool>("SkipEnergyCut")),
65  fSliceNCellLow(pset.get<double>("SliceNCellLow")),
66  fSliceNCellHigh(pset.get<double>("SliceNCellHigh")),
67  fSliceExtentPlane(pset.get<double>("SliceExtentPlane")),
68  fTrackLengthCut(pset.get<double>("TrackLengthCut")),
69  fRecoELow(pset.get<double>("RecoELow")),
70  fRecoEHigh(pset.get<double>("RecoEHigh")),
71  fInFidCut(pset.get<bool>("InFidCut"))
72  {
73 
74  produces<rb::FilterList<rb::Cluster> >();
75  produces<std::vector<rb::PID> >();
76  produces<art::Assns<rb::PID, rb::Cluster> >();
77 
78  }
79 
80  //----------------------------------------------------------------------
82  {
83  }
84 
85  //----------------------------------------------------------------------
87  {
88  std::unique_ptr<rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
89 
90  std::unique_ptr<std::vector<rb::PID> > pidcol(new std::vector<rb::PID>);
91  std::unique_ptr<art::Assns<rb::PID, rb::Cluster> > assns(new art::Assns<rb::PID, rb::Cluster>);
92 
95 
96  // TODO: break truth selection out into its own module
97  /*
98  if(fTruthMode){
99  art::Handle<std::vector<simb::MCTruth> > truths;
100  evt.getByLabel("generator", truths);
101  if(truths->empty()) return false;
102  if(truths->size() > 1) return true;
103 
104  const simb::MCTruth& truth = (*truths)[0];
105 
106  const simb::MCNeutrino& nu = truth.GetNeutrino();
107  const double trueE = nu.Nu().E();
108  // These cuts are empirically observed to remove almost nothing that
109  // would survive the preselection below. Approximately 0.1% lost to the
110  // low cuts and 0.5% to the high cuts.
111  if(nu.CCNC() == simb::kCC){
112  return (trueE > 1 && trueE < 3.5);
113  }
114  else{
115  const double hadY = nu.Y();
116  return (trueE*hadY > .75 && trueE*hadY < 4.5);
117  }
118 
119  assert(0 && "Not reached");
120  }
121  */
122 
124  evt.getByLabel(fSlicerLabel, slices);
125 
126  const int sliceMax = slices->size();
127  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
128  if(rb::IsFiltered(evt, slices, sliceIdx)) continue;
129 
130  const rb::Cluster& slice = (*slices)[sliceIdx];
131 
132  if(slice.IsNoise()) continue;
133 
134  double longestTrack = 0;
135  if(!fSkipTrackCut){
136  // This is maybe inefficient, but for now FindManyP throws an exception
137  // on finding no assns. We don't want that.
138  art::FindManyP<rb::Track> fmt(slices, evt, std::string("discretetrack"));
139  // Get all the tracks associated with this slice
140  std::vector<art::Ptr<rb::Track> > tracks = fmt.at(sliceIdx);
141 
142  for(unsigned int n = 0; n < tracks.size(); ++n){
143  const double L = tracks[n]->TotalLength();
144  if(L > longestTrack) longestTrack = L;
145  }
146  }
147 
148  const double nCells = slice.NCell();
149  const int lengthPlanes = slice.ExtentPlane();
150 
151  const TVector3 mean = slice.MeanXYZ();
152  const bool inFid = geom->isInsideFiducialVolume(mean.X(), mean.Y(), mean.Z());
153 
154  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
155  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
156 
157  double totalRecoGeV = 0;
158  for(unsigned int hitIdx = 0; hitIdx < slice.NCell(); ++hitIdx){
159  // This is slow because we have to keep recalculating the mean
160  // positions. The manual version below is much better.
161  // const rb::RecoHit rhit = slice.RecoHit(hitIdx);
162 
163  const art::Ptr<rb::CellHit>& chit = slice.Cell(hitIdx);
164  const rb::RecoHit rhit = cal->MakeRecoHit(*chit,
165  chit->View() == geo::kX ? wx : wy);
166  if(!rhit.IsCalibrated()){
167  std::cout << "Not calibrated?!" << std::endl;
168  continue;
169  }
170  totalRecoGeV += rhit.GeV();
171  }
172 
173  bool sel = true;
174  if(nCells < fSliceNCellLow || nCells > fSliceNCellHigh) sel = false;
175  if(longestTrack > fTrackLengthCut) sel = false;
176  if(lengthPlanes > fSliceExtentPlane) sel = false;
177  if(!fSkipEnergyCut) {
178  if(2*totalRecoGeV < fRecoELow || 2*totalRecoGeV > fRecoEHigh) sel = false;
179  }
180  if(!inFid && fInFidCut == true) sel = false;
181 
182  if(!sel) filtcol->Add(slices, sliceIdx);
183 
184  // We're a nue preselection
185  pidcol->push_back(rb::PID(12, sel));
186 
187  util::CreateAssn(*this, evt, *pidcol, art::Ptr<rb::Cluster>(slices, sliceIdx), *assns);
188  } // end for sliceIdx
189 
190  evt.put(std::move(pidcol));
191  evt.put(std::move(assns));
192  evt.put(std::move(filtcol));
193  }
194 
196 
197 } // end namespace
A simple list of products that have been marked "filtered out".
Definition: FilterList.h:74
A pid value and corresponding pdg code.
Definition: PID.h:13
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
pdg code and pid value
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
DEFINE_ART_MODULE(TestTMapFile)
virtual void produce(art::Event &evt)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
PID
Definition: FillPIDs.h:13
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
static constexpr double L
Encapsulate the geometry of one entire detector (near, far, ndos)
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
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
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
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
float GeV() const
Definition: RecoHit.cxx:69
Preselection(const fhicl::ParameterSet &pset)
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
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
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?