FindLEMMatches_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file FindMatches_module.cc
3 // \brief Module to find best LEM matches
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "LEM/func/DistanceMap.h"
9 #include "LEM/func/Library.h"
10 #include "LEM/func/MatchIndices.h"
11 #include "LEM/func/MatchList.h"
12 
13 #include "LEM/LEMInput.h"
14 #include "LEM/Util.h"
15 
16 // ROOT includes
17 #include "TH2.h"
18 
19 // NOvASoft includes
21 #include "Geometry/Geometry.h"
22 #include "RecoBase/CellHit.h"
23 #include "RecoBase/Cluster.h"
24 #include "RecoBase/FilterList.h"
25 #include "RecoBase/Vertex.h"
26 #include "Utilities/AssociationUtil.h"
28 #include "Utilities/func/MathUtil.h"
29 
30 // Framework includes
38 #include "fhiclcpp/ParameterSet.h"
39 
40 /// \f$ \nu_e \f$ %PID
41 namespace lem
42 {
43  /// Module to find best %LEM matches
44  class FindLEMMatches: public art::EDProducer // TODO rename once FindMatches is gone
45  {
46  public:
47  explicit FindLEMMatches(const fhicl::ParameterSet& pset);
49 
50  virtual void reconfigure(const fhicl::ParameterSet& pset);
51  virtual void beginRun(art::Run& run);
52  virtual void produce(art::Event& evt);
53  protected:
54  void AddHists(std::vector<TH2F>* th2col, const Match& bestMatch,
55  const EventSummary& trial,
56  const EventSummary* head,
57  bool headFlipEven, bool headFlipOdd) const;
58 
60  int vtxPlane, int vtxCell, int vtxCellOther) const;
61 
63 
65 
67  };
68 
69  //......................................................................
71  {
72  reconfigure(pset);
73 
74  produces<std::vector<lem::MatchList> >();
75  // MatchLists are associated with the slice they are against
76  produces<art::Assns<lem::MatchList, rb::Cluster> >();
77 
78  produces<std::vector<lem::MatchIndices> >();
79  // MatchLists are associated with the slice they are against
80  produces<art::Assns<lem::MatchIndices, rb::Cluster> >();
81 
82  // True vertex of best match. Associate to slice
83  produces<std::vector<rb::Vertex> >();
84  produces<art::Assns<rb::Vertex, rb::Cluster> >();
85 
86  // Really should be file level, not run level, but don't know how to do
87  // that
88  produces<lem::LibrarySummary, art::InRun>();
89 
90  // Just to mark where the vertex is so you can check it's right
91  produces<std::vector<geo::OfflineChan> >();
92 
93  // Images of best match
95  produces<std::vector<TH2F> >();
96  }
97 
98  //......................................................................
100  {
101  delete fAlg;
102  }
103 
104  //......................................................................
106  {
107  fAlg = new FindMatchesAlg(util::EnvExpansion(pset.get<std::string>("LibraryPath")),
108  pset.get<bool>("PreloadLib"),
109  pset.get<bool>("UseHeads"));
110 
111  fInputLabel = pset.get<std::string>("InputLabel");
112 
113  // This needs to be before anyone asks the DistanceMap anything, should be
114  // safe enough here.
115  DistanceMap::SetPlaneScale(pset.get<double>("DistPlaneScale"));
116  DistanceMap::SetCellScale(pset.get<double>("DistCellScale"));
117  DistanceMap::SetDecayPower(pset.get<double>("DistDecayPower"));
118  DistanceMap::SetExpMode(pset.get<bool>("DistExpMode"));
119 
120  EventSummary::fgDecayRate = pset.get<double>("HitDecayRate");
121  EventSummary::fgChargePower = pset.get<double>("HitChargePower");
122 
123  fSaveEventImages = pset.get<bool>("SaveEventImages", true);
124  fSaveEventPotentials = pset.get<bool>("SaveEventPotentials", false);
125  }
126 
127  //......................................................................
129  {
130  std::unique_ptr<lem::LibrarySummary> ls(new lem::LibrarySummary);
131  *ls = fAlg->GetLibrary()->Summary();
132  run.put(std::move(ls));
133  }
134 
135  //......................................................................
137  {
138  std::unique_ptr<std::vector<lem::MatchList> > matchcol(new std::vector<lem::MatchList>);
139  std::unique_ptr<std::vector<lem::MatchIndices>> matchidxs(new std::vector<lem::MatchIndices>);
140 
141  // Associate matches with the slice they came from
142  std::unique_ptr<art::Assns<lem::MatchList, rb::Cluster> > assns(new art::Assns<lem::MatchList, rb::Cluster>);
143  std::unique_ptr<art::Assns<lem::MatchIndices, rb::Cluster> > idxassns(new art::Assns<lem::MatchIndices, rb::Cluster>);
144 
145  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
146  std::unique_ptr<art::Assns<rb::Vertex, rb::Cluster> > vtxassns(new art::Assns<rb::Vertex, rb::Cluster>);
147 
148  std::unique_ptr<std::vector<geo::OfflineChan> > chancol(new std::vector<geo::OfflineChan>);
149 
150  std::unique_ptr<std::vector<TH2F>> th2col(new std::vector<TH2F>);
151 
152 
154  evt.getByLabel(fInputLabel, inputs);
156 
157  const int inputMax = inputs->size();
158  for(int inputIdx = 0; inputIdx < inputMax; ++inputIdx){
159  const LEMInput& input = (*inputs)[inputIdx];
160 
161  EventSummary sum(input.hits, -1, std::vector<int>());
162 
163  sum.run = evt.run();
164  sum.subrun = evt.subRun();
165  sum.event = evt.event();
166 
167  sum.origPdg = 0;
168 
169  sum.pdg = 0;
170  sum.ccnc = 0;
171  sum.y = 0;
172  sum.trueEVis = 0;
173 
174 
175  // chancol->push_back(geo::OfflineChan(vtxPlane, vtxCell));
176  // Should really be nextPlane, but that's lost inside of Util.cxx. This
177  // is almost always right, and only for display purposes anyway.
178  // chancol->push_back(geo::OfflineChan(vtxPlane+1, vtxCellOther));
179 
180  MatchableEvent trial(sum);
181 
182  const std::vector<Match> matches = fAlg->FindMatches(trial, kMaxNumMatches, false);
183  const std::vector<Match> matchesEnrich = fAlg->FindMatches(trial, kMaxNumMatches, true);
184 
185 /*
186  if(!matches.empty()){
187  vtxcol->push_back(MatchToVertex(matches[0], vtxPlane, vtxCell, vtxCellOther));
188  vtxcol->back().SetT(slice.MinTNS());
189  util::CreateAssn(*this, evt, *vtxcol,
190  art::Ptr<rb::Cluster>(slices, sliceIdx), *vtxassns);
191  }
192 */
193 
194  // TODO TODO TODO
195 /*
196  if(!matches.empty()){
197  AddHists(th2col.get(), matches[0], sum,
198  (bestHead >= 0) ? &fLib->Event(fHeads->HeadIdx(bestHead)) : 0,
199  headFlipEven, headFlipOdd);
200  }
201 */
202 
203  matchcol->push_back(MatchList());
204  MatchList& ml = matchcol->back();
205 
206  matchidxs->push_back(MatchIndices());
207  MatchIndices& mi = matchidxs->back();
208  // TODO TODO TODO
209  // mi.headIdx = (bestHead >= 0) ? fHeads->HeadIdx(bestHead) : -1;
210 
211  for(unsigned int i = 0; i < matches.size(); ++i){
212  ml.matches.push_back(MatchSummary(matches[i]));
213  mi.matches.push_back(matches[i].evt - &fAlg->GetLibrary()->Event(0));
214  } // end for i
215 
216  for(unsigned int i = 0; i < matchesEnrich.size(); ++i){
217  ml.matches.push_back(MatchSummary(matchesEnrich[i]));
218  mi.matchesEnrich.push_back(matchesEnrich[i].evt - &fAlg->GetLibrary()->Event(0));
219  } // end for i
220 
221  ml.trial = MatchSummary(Match(0, &trial, false, false));
222 
223  // TODO TODO TODO
224  // ml.headIdx = (bestHead >= 0) ? fHeads->HeadIdx(bestHead) : -1;
225 
226  // If the original input hits were associated to a cluster, then so
227  // should our matches be.
228  try{
229  std::vector<art::Ptr<rb::Cluster>> slices = fmt.at(inputIdx);
230  if(slices.size() == 1){
231  util::CreateAssn(*this, evt, *matchcol, slices[0], *assns);
232  util::CreateAssn(*this, evt, *matchidxs, slices[0], *idxassns);
233  }
234  }
235  catch(...){
236  // ugh
237  }
238  } // end for sliceIdx
239 
240  evt.put(std::move(matchcol));
241  evt.put(std::move(assns));
242 
243  evt.put(std::move(vtxcol));
244  evt.put(std::move(vtxassns));
245 
246  evt.put(std::move(chancol));
247  if(fSaveEventImages || fSaveEventPotentials) evt.put(std::move(th2col));
248  }
249 
250  //......................................................................
251  void FindLEMMatches::AddHists(std::vector<TH2F>* th2col,
252  const Match& bestMatch,
253  const EventSummary& trial,
254  const EventSummary* head,
255  bool headFlipEven, bool headFlipOdd) const
256  {
257  // Best
258  const EventSummary* bestEvt = bestMatch.evt;
259 
260  const std::string bestTitle = bestEvt->Description()+";Plane;Cell";
261 
262  if(fSaveEventImages){
263  TH2F evenBest("bestEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192);
264  TH2F oddBest ("bestOdd", bestTitle.c_str(), 64, 0, 128, 128, 64, 192);
265 
266  bestEvt->FillHists(&evenBest, &oddBest, bestMatch.flipEven, bestMatch.flipOdd);
267 
268  th2col->push_back(evenBest);
269  th2col->push_back(oddBest);
270  }
271 
273  TH2F hVsBest[2] = {TH2F("bestVEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192),
274  TH2F("bestVOdd", bestTitle.c_str(), 64, 0, 128, 128, 64, 192)};
275  Potential VBest;
276  FillPotential(*bestEvt, VBest, bestMatch.flipEven, bestMatch.flipOdd);
277  for(int plane = 0; plane < 256; ++plane){
278  for(int cell = 0; cell < 256; ++cell){
279  const double p = VBest.V[256*plane+cell];
280  hVsBest[plane%2].Fill(plane, cell, p);
281  }
282  }
283  th2col->push_back(hVsBest[0]);
284  th2col->push_back(hVsBest[1]);
285  }
286 
287 
288  // Head
289  if(head){
290  const char* headTitle = (head->Description()+";Plane;Cell").c_str();
291 
292  if(fSaveEventImages){
293  TH2F evenHead("headEven", headTitle, 64, 0, 128, 128, 64, 192);
294  TH2F oddHead ("headOdd", headTitle, 64, 0, 128, 128, 64, 192);
295 
296  head->FillHists(&evenHead, &oddHead, headFlipEven, headFlipOdd);
297 
298  th2col->push_back(evenHead);
299  th2col->push_back(oddHead);
300  }
301 
303  TH2F hVsHead[2] = {TH2F("headVEven", bestTitle.c_str(), 64, 0, 128, 128, 64, 192),
304  TH2F("headVOdd", bestTitle.c_str(), 64, 0, 128, 128, 64, 192)};
305  Potential VHead;
306  FillPotential(*head, VHead, headFlipEven, headFlipOdd);
307  for(int plane = 0; plane < 256; ++plane){
308  for(int cell = 0; cell < 256; ++cell){
309  const double p = VHead.V[256*plane+cell];
310  hVsHead[plane%2].Fill(plane, cell, p);
311  }
312  }
313  th2col->push_back(hVsHead[0]);
314  th2col->push_back(hVsHead[1]);
315  }
316  }
317 
318 
319  // Trial
320  const char* trialTitle = "Trial event;Plane;Cell";
321 
322  if(fSaveEventImages){
323  TH2F evenTrial("trialEven", trialTitle, 64, 0, 128, 128, 64, 192);
324  TH2F oddTrial ("trialOdd", trialTitle, 64, 0, 128, 128, 64, 192);
325 
326  trial.FillHists(&evenTrial, &oddTrial);
327 
328  th2col->push_back(evenTrial);
329  th2col->push_back(oddTrial);
330  }
331 
333  TH2F hVsTrial[2] = {TH2F("trialVEven", trialTitle, 64, 0, 128, 128, 64, 192),
334  TH2F("trialVOdd", trialTitle, 64, 0, 128, 128, 64, 192)};
335  Potential VTrial;
336  FillPotential(trial, VTrial, false, false);
337  for(int plane = 0; plane < 256; ++plane){
338  for(int cell = 0; cell < 256; ++cell){
339  const double p = VTrial.V[256*plane+cell];
340  hVsTrial[plane%2].Fill(plane, cell, p);
341  }
342  }
343  th2col->push_back(hVsTrial[0]);
344  th2col->push_back(hVsTrial[1]);
345  }
346  }
347 
348  //......................................................................
350  int vtxPlane, int vtxCell,
351  int vtxCellOther) const
352  {
354 
355  int plane = vtxPlane+match.evt->trueVtxPlane-kVertexPlane;
356 
357  if(plane < 0) plane = 0;
358  if(plane > int(geom->NPlanes())) plane = geom->NPlanes()-1;
359 
360  int cell, cellOther;
361  if(abs(match.evt->trueVtxPlane-kVertexPlane)%2 == 0){
362  cell = vtxCell+match.evt->trueVtxCell-kVertexCell;
363  cellOther = vtxCellOther+match.evt->trueVtxCellOther-kVertexCell;
364  }
365  else{
366  cell = vtxCellOther+match.evt->trueVtxCell-kVertexCell;
367  cellOther = vtxCell+match.evt->trueVtxCellOther-kVertexCell;
368  }
369 
370  int planeOther = 0;
371  if(geom->Plane(plane)){
372  geom->NextPlaneOtherView(plane, +1);
373  if(planeOther == geo::kPLANE_NOT_FOUND) planeOther = geom->NextPlaneOtherView(plane, -1);
374  }
375 
376  if(cell < 0) cell = 0;
377  if(cellOther < 0) cellOther = 0;
378  if(cell >= int(geom->Plane(plane)->Ncells())) cell = geom->Plane(plane)->Ncells()-1;
379  if(cellOther >= int(geom->Plane(planeOther)->Ncells())) cellOther = geom->Plane(planeOther)->Ncells()-1;
380 
381  double xyz[3];
382  geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
383 
384  double xyz2[3];
385  geom->Plane(planeOther)->Cell(cellOther)->GetCenter(xyz2);
386 
387  const geo::View_t v = geom->Plane(plane)->View();
388 
389  return rb::Vertex((v == geo::kX) ? xyz[0] : xyz2[0],
390  (v == geo::kY) ? xyz[1] : xyz2[1],
391  xyz[2],
392  -1); // Don't bother to fill the time
393  }
394 
396 
397 } //namespace lem
398 ////////////////////////////////////////////////////////////////////////
virtual void reconfigure(const fhicl::ParameterSet &pset)
std::string Description() const
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
SubRunNumber_t subRun() const
Definition: Event.h:72
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.
std::vector< MatchSummary > matches
Definition: MatchList.h:24
virtual void beginRun(art::Run &run)
Simple representation of event for LEM use.
Definition: EventSummary.h:26
static double fgChargePower
Definition: EventSummary.h:76
std::vector< int > matches
Definition: MatchIndices.h:19
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
Attach some information used in matching to an EventSummary.
const unsigned int kMaxNumMatches
Definition: FindMatches.h:29
Map of electrostatic potential at each cell.
Definition: Potential.h:13
Collection of events for matching.
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
rb::Vertex MatchToVertex(const Match &match, int vtxPlane, int vtxCell, int vtxCellOther) const
const MatchableEvent * evt
Definition: Match.h:26
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
bool flipOdd
Definition: Match.h:27
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
const LibrarySummary & Summary() const
Definition: Library.h:30
const int kVertexPlane
Definition: EventSummary.h:22
std::vector< int > matchesEnrich
Definition: MatchIndices.h:20
void AddHists(std::vector< TH2F > *th2col, const Match &bestMatch, const EventSummary &trial, const EventSummary *head, bool headFlipEven, bool headFlipOdd) const
void abs(TH1 *hist)
std::vector< Match > FindMatches(const MatchableEvent &trial, unsigned int numMatches, int enrich) const
const PlaneGeo * Plane(unsigned int i) const
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
Eigen::Matrix< T, Eigen::Dynamic, 1 > head(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &v, size_t n)
Definition: head.hpp:24
Definition: Run.h:31
Collection of MatchSummary objects.
PID
Definition: FillPIDs.h:14
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Simplified Match information, suitable for serialization.
Definition: MatchSummary.h:16
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
static void SetPlaneScale(double ps)
Must call before first call to Instance()
Definition: DistanceMap.h:32
Collection of MatchSummary objects.
Definition: MatchList.h:17
const Library * GetLibrary() const
const int kVertexCell
Definition: EventSummary.h:23
T get(std::string const &key) const
Definition: ParameterSet.h:231
static void SetDecayPower(double dp)
Definition: DistanceMap.h:34
Information about a LEM match.
Definition: Match.h:17
EventNumber_t event() const
Definition: Event.h:67
Details of the library LEM matches were made against.
virtual void produce(art::Event &evt)
Vertex location in position and time.
Definition: run.py:1
void FillPotential(const EventSummary &trial, Potential &V, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:76
Module to find best LEM matches.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
FindLEMMatches(const fhicl::ParameterSet &pset)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
static void SetExpMode(bool em=true)
Definition: DistanceMap.h:35
std::vector< LiteHit > hits
Definition: LEMInput.h:14
void geom(int which=0)
Definition: geom.C:163
const MatchableEvent & Event(int i) const
Definition: Library.h:32
static double fgDecayRate
Definition: EventSummary.h:75
static void SetCellScale(double cs)
Definition: DistanceMap.h:33
def ls(target="")
Definition: g4zmq.py:69
MatchSummary trial
Definition: MatchList.h:25
unsigned int NPlanes() const
Double_t sum
Definition: plot.C:31
Calculate and cache electrostatic potential between cells.
RunNumber_t run() const
Definition: Event.h:77
bool flipEven
Definition: Match.h:27
List of indices of matched events.
float V[256 *256]
Definition: Potential.h:18
void FillHists(TH2 *h1, TH2 *h2, bool flipEven=false, bool flipOdd=false, bool idxs=false) const
Encapsulate the geometry of one entire detector (near, far, ndos)
List of indices of matched events.
Definition: MatchIndices.h:15
const unsigned int NextPlaneOtherView(unsigned int p, int d=+1) const