FindMatches.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file FindMatches.cxx
3 /// \brief Core of the LEM match-finding algorithm
4 /// \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include "LEM/func/FindMatches.h"
8 
9 #include "LEM/func/DistanceMap.h"
10 #include "LEM/func/Heads.h"
12 
13 #include <algorithm>
14 #include <cfloat>
15 #include <iostream>
16 #include <set>
17 
18 namespace lem
19 {
20  //......................................................................
21  float CalcEnergy(const Potential& Va, const MatchableEvent& a, const MatchableEvent& b)
22  {
23  float energy = a.selfEnergy + b.selfEnergy;
24 
25  const LiteHit* hb = b.hits;
26  const unsigned int Nb = b.nhits;
27 
28  for(unsigned int i = 0; i < Nb; ++i){
29  energy -= hb->pecorr*Va.V[hb->cellIdx];
30  ++hb;
31  }
32 
33  return energy;
34  }
35 
36  //......................................................................
38  bool flipEven, bool flipOdd)
39  {
40  // TODO: Duplicate hits definitely screw this up
41 
42  float qmatched = 0;
43 
44  // Quicker to spend this (not much) memory than the more obvious
45  // nested-loops implementation.
46  float acharge[256*256] = {0,};
47 
48  const unsigned int Na = a.nhits;
49  const unsigned int Nb = b.nhits;
50  for(unsigned int ia = 0; ia < Na; ++ia){
51  const LiteHit& ha = a.hits[ia];
52  acharge[ha.cellIdx] += ha.pecorr;
53  } // end for ia
54 
55  for(unsigned int ib = 0; ib < Nb; ++ib){
56  const LiteHit& hb = b.hits[ib];
57 
58  int hbc = hb.Cell();
59 
60  const int v = abs(hb.Plane())%2; // Matches definition in FillPotential
61  if((v == 0 && flipEven) ||
62  (v == 1 && flipOdd)) hbc = 2*lem::kVertexCell-hbc;
63 
64  const int bidx = 256*hb.Plane()+hbc;
65 
66  // This is the MINOS definition
67  // if(acharge[bidx]) qmatched += acharge[bidx]+hb.pecorr;
68 
69  qmatched += 2*std::min(acharge[bidx], hb.pecorr);
70  } // end for ib
71 
72  return qmatched/(a.totalPE+b.totalPE);
73  }
74 
75  //......................................................................
76  void FillPotential(const EventSummary& trial, Potential& V,
77  bool flipEven, bool flipOdd)
78  {
79  const DistanceMap& dm = DistanceMap::Instance();
80 
81  const unsigned int N = trial.nhits;
82 
83  for(unsigned int i = 0; i < N; ++i){
84  const LiteHit& h = trial.hits[i];
85 
86  // Doesn't matter? Has to match in qfrac definition
87  const int view = abs(h.Plane())%2;
88 
89  int hc = h.Cell();
90  if((flipEven && view == 0) || (flipOdd && view == 1))
91  hc = 2*lem::kVertexCell-hc;
92 
93  for(int plane = 0; plane < 256; ++plane){
94  // Only in the same view
95  if(abs(plane-h.Plane())%2) continue;
96  for(int cell = 0; cell < 256; ++cell){
97 
98  const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
99  cell, hc);
100 
101  V.V[256*plane+cell] += contrib;
102  }
103  }
104  }
105  }
106 
107  //......................................................................
109  {
110  // We could just call FillPotential 4 times with each flip variant like so
111  /*
112  for(bool flipEven: {false, true}){
113  for(bool flipOdd: {false, true}){
114  FillPotential(trial, Vs.V[flipEven][flipOdd], flipEven, flipOdd);
115  }
116  }
117  */
118 
119  // But this function actually takes noticeable runtime, so instead just
120  // fill it once, and then copy into the over three flips.
121  const DistanceMap& dm = DistanceMap::Instance();
122 
123  const unsigned int N = trial.nhits;
124 
125  // There's a complication because 256 is odd there isn't a perfect
126  // centreline to flip round. The way things are defined, we end up needing
127  // to know what cell 256 in the original map would have been, which we
128  // don't have. Make a place to store it here.
129  float overflow[256] = {0,};
130 
131  for(unsigned int i = 0; i < N; ++i){
132  const LiteHit& h = trial.hits[i];
133 
134  for(int plane = 0; plane < 256; ++plane){
135  // Only in the same view
136  if(abs(plane-h.Plane())%2) continue;
137  for(int cell = 0; cell < 256; ++cell){
138 
139  const double contrib = h.pecorr*dm.InvDist(plane, h.Plane(),
140  cell, h.Cell());
141 
142  Vs.V[false][false].V[256*plane+cell] += contrib;
143  }
144  // And one for cell 256, stored seperately
145  overflow[plane] += h.pecorr*dm.InvDist(plane, h.Plane(),
146  256, h.Cell());
147  }
148  }
149 
150  // Now copy over into the flipped variants
151 
152  for(bool flipEven: {false, true}){
153  for(bool flipOdd: {false, true}){
154  if(!flipEven && !flipOdd) continue;
155 
156  for(int plane = 0; plane < 256; ++plane){
157  const int view = plane%2;
158 
159  for(int cell = 0; cell < 256; ++cell){
160 
161  int fc = cell; // flipped cell
162  if((flipEven && view == 0) || (flipOdd && view == 1))
163  fc = 2*lem::kVertexCell-fc;
164 
165  if(fc >= 0 && fc < 256)
166  Vs.V[flipEven][flipOdd].V[plane*256+cell] = Vs.V[false][false].V[plane*256+fc];
167  if(fc == 256) Vs.V[flipEven][flipOdd].V[plane*256+cell] = overflow[plane];
168  }
169  }
170  }
171  }
172  }
173 
174  //......................................................................
175  std::vector<Match> FindMatches(const MatchableEvent& trial,
176  unsigned int libSize,
177  const MatchableEvent* lib,
178  const MatchableEvent* libDownsample,
179  int factor,
180  const std::multiset<float>& alreadyInput,
181  unsigned int numMatches,
182  int enrich,
183  const Heads* heads,
184  int headIdx,
185  bool flipHeadEven,
186  bool flipHeadOdd,
187  bool useDownsample)
188  {
189  // assert(libsize >= kMaxNumMatches);
190 
191  assert(!heads || headIdx >= 0);
192 
193  float target = FLT_MAX; // infinity
194 
195  // "De-normalize" values back into the range most of this function deals
196  // with. Need to copy to a new container, since set iterators are const
197  std::multiset<float> already;
198  for(auto a: alreadyInput) already.insert(a*trial.selfEnergy);
199 
200  // Truncate to the right length. Update target if possible.
201  while(already.size() > numMatches){
202  auto itLargest = already.end();
203  --itLargest;
204  already.erase(itLargest);
205  target = *already.rbegin();
206  }
207 
208  // Precalculate potential
209  // Try it with either/both of x&y flipped. Effectively gain factor of 4 library stats.
210  // TODO: This isn't entirely right because of attenuation. It /would/
211  // actually be possible to introduce this properly during library building
212  // when truth information is available.
213 
215  FillPotential(trial, Vs);
216 
217  // Make downsampled versions of the four potential cases
218  FlippedPotentials Vs_down;
219  if(useDownsample){
220  for(int flipEven = 0; flipEven <= 1; ++flipEven){
221  for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
222  Vs_down.V[flipEven][flipOdd] = Vs.V[flipEven][flipOdd].Downsampled(factor);
223  }
224  }
225  }
226 
227  // Candidates for return, remains sorted
228  std::multiset<Match> Elist;
229 
230  // This is the quickest ordering of the loops (locality into lib?)
231  for(unsigned int j = 0; j < libSize; ++j){
232  int i = j;
233  if(heads){
234  if(int(j) >= heads->HeadSeqsLen()) break;
235  i = heads->ChildIdxFor(headIdx, j);
236  }
237 
238  if(enrich != 2){
239  if(lib[i].enrich != bool(enrich)) continue;
240  }
241 
242  const MatchableEvent& downEvt = libDownsample[i];
243 
244  // There's beam nue in the library, but we don't want to match to it
245  if( useDownsample && downEvt.weight == 0) continue;
246  if(!useDownsample && lib[i].weight == 0) continue;
247 
248  for(int flipEven = 0; flipEven <= 1; ++flipEven){
249  if(heads && (heads->FlipEvenFor(headIdx, j) ^ flipHeadEven) != flipEven) continue;
250  for(int flipOdd = 0; flipOdd <= 1; ++flipOdd){
251  if(heads && (heads->FlipOddFor(headIdx, j) ^ flipHeadOdd) != flipOdd) continue;
252 
253  // Compare the downsampled event to the downsampled potential. The
254  // resulting energy is a lower-bound on what the real match would
255  // have done.
256  float E = useDownsample ? CalcEnergy(Vs_down.V[flipEven][flipOdd], trial, downEvt) : 0;
257 
258  if(E < target){
259  // Calculate the energy fully
260  E = CalcEnergy(Vs.V[flipEven][flipOdd], trial, lib[i]);
261 
262  // If it's still good, use it
263  if(E < target){
264  Elist.insert(Match(E, &lib[i], flipEven, flipOdd, j));
265  already.insert(E);
266 
267  if(Elist.size() > numMatches){
268  // Always drop the last one, to keep list correct size
269  auto itLargest = Elist.end();
270  --itLargest;
271  Elist.erase(itLargest);
272  // In case "already" isn't long enough
273  target = Elist.rbegin()->energy;
274  }
275  if(already.size() > numMatches){
276  // Always drop the last one, to keep list correct size
277  auto itLargest = already.end();
278  --itLargest;
279  already.erase(itLargest);
280  target = *already.rbegin();
281  }
282  }
283  }
284  } // end for flipOdd
285  } // end for flipEven
286  } // end for i
287 
288  // Convert back to a vector, and calculate the qfrac variable for surviving
289  // matches.
290  std::vector<Match> ret;
291  ret.reserve(numMatches);
292  for(const auto& m: Elist){
293  ret.push_back(m);
294  ret.back().qfrac = FracChargeMatched(trial, *m.evt, m.flipEven, m.flipOdd);
295  ret.back().energy /= trial.selfEnergy; // Normalize to approx 0-1
296  }
297  return ret;
298  }
299 
300 } // namespace
unsigned short nhits
Definition: EventSummary.h:53
Simple representation of event for LEM use.
Definition: EventSummary.h:26
float CalcEnergy(const Potential &Va, const MatchableEvent &a, const MatchableEvent &b)
Definition: FindMatches.cxx:21
Attach some information used in matching to an EventSummary.
const XML_Char * target
Definition: expat.h:268
Map of electrostatic potential at each cell.
Definition: Potential.h:13
const Var weight
const int Na
int ChildIdxFor(int head, int child) const
Definition: Heads.cxx:78
void abs(TH1 *hist)
float pecorr
Definition: LiteHit.h:24
int Cell() const
Definition: LiteHit.h:39
bool FlipEvenFor(int head, int child) const
Definition: Heads.cxx:85
Mapping from a subset of the library to their best matches.
Definition: __init__.py:1
PID
Definition: FillPIDs.h:14
Core of the LEM match-finding algorithm.
Potential V[2][2]
Definition: FindMatches.h:23
std::vector< Match > FindMatches(const MatchableEvent &trial, unsigned int libSize, const MatchableEvent *lib, const MatchableEvent *libDownsample, int factor, const std::multiset< float > &alreadyInput, unsigned int numMatches, int enrich, const Heads *heads, int headIdx, bool flipHeadEven, bool flipHeadOdd, bool useDownsample)
bool FlipOddFor(int head, int child) const
Definition: Heads.cxx:92
int Plane() const
Definition: LiteHit.h:38
unsigned short cellIdx
Definition: LiteHit.h:25
static const DistanceMap & Instance()
Singleton.
Definition: DistanceMap.cxx:26
Float_t E
Definition: plot.C:20
const double a
const int kVertexCell
Definition: EventSummary.h:23
double energy
Definition: plottest35.C:25
double FracChargeMatched(const EventSummary &a, const EventSummary &b, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:37
Mapping from a subset of the library to their best matches.
Definition: Heads.h:16
const double j
Definition: BetheBloch.cxx:29
Information about a LEM match.
Definition: Match.h:17
int HeadSeqsLen() const
Definition: Heads.h:28
Attach some information used in matching to an EventSummary.
void FillPotential(const EventSummary &trial, Potential &V, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:76
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
Collection of Potential objects with odd and/or even view flipped.
Definition: FindMatches.h:21
const int Nb
T min(const caf::Proxy< T > &a, T b)
Potential Downsampled(int factor) const
Definition: Potential.cxx:23
Compressed hit info, basic component of LEM events.
Definition: LiteHit.h:18
Calculate and cache electrostatic potential between cells.
double InvDist(int aplane, int bplane, int acell, int bcell) const
Definition: DistanceMap.cxx:99
float V[256 *256]
Definition: Potential.h:18
Calculate and cache electrostatic potential between cells.
Definition: DistanceMap.h:20