FindMatchesAlg.cxx
Go to the documentation of this file.
2 
3 #include "LEM/func/FindMatches.h" // For FillPotentials
4 
5 #include <cassert>
6 #include <cfloat>
7 
8 namespace lem
9 {
10  // TODO investigate other data structures
11  class SortedMatchList: public std::multiset<Match>
12  {
13  public:
14  SortedMatchList(unsigned int numMatches)
15  : fNumMatches(numMatches), fTarget(FLT_MAX) {}
16 
17  template<class... T> void emplace(double E, T... args)
18  {
19  if(E >= fTarget) return;
20 
21  std::multiset<Match>::emplace(E, args...);
22 
23  if(size() > fNumMatches){
24  // Always drop the last one, to keep list correct size
25  auto itLargest = end();
26  --itLargest;
27  erase(itLargest);
28  fTarget = rbegin()->energy;
29  }
30  }
31 
32  protected:
33  void insert(Match m); // not implemented
34 
35  unsigned int fNumMatches;
36  double fTarget;
37  };
38 
39  //......................................................................
41  bool preload,
42  bool useHeads)
43  {
44  fLib = Library::FromMMap(libDir, preload);
45  assert(fLib);
46 
47  if(useHeads){
48  fHeads = Heads::FromMMap(libDir, preload);
49  assert(fHeads);
50  }
51  else{
52  fHeads = 0;
53  }
54  }
55 
56  //......................................................................
58  {
59  delete fLib;
60  delete fHeads;
61  }
62 
63  //......................................................................
64  std::vector<Match> FindMatchesAlg::
65  PackageMatches(const std::multiset<Match>& matches,
66  const MatchableEvent& trial) const
67  {
68  // Convert back to a vector, and calculate the qfrac variable for surviving
69  // matches.
70  std::vector<Match> ret;
71  ret.reserve(matches.size());
72  for(const Match& m: matches){
73  ret.push_back(m);
74  ret.back().qfrac = FracChargeMatched(trial, *m.evt, m.flipEven, m.flipOdd);
75  ret.back().energy /= trial.selfEnergy; // Normalize to approx 0-1
76  }
77 
78  return ret;
79  }
80 
81  //......................................................................
82  std::vector<Match> FindMatchesAlg::
84  unsigned int numMatches,
85  int enrich) const
86  {
87  if(fHeads)
88  return FindMatchesHeads (trial, numMatches, enrich);
89  else
90  return FindMatchesNoHeads(trial, numMatches, enrich);
91  }
92 
93  //......................................................................
94  std::vector<Match> FindMatchesAlg::
96  unsigned int numMatches,
97  int enrich) const
98  {
99  FlippedPotentials trialVs;
100  FillPotential(trial, trialVs);
101 
102  // Candidates for return, remains sorted
103  SortedMatchList Elist(numMatches);
104 
105  // This is the quickest ordering of the loops (locality into lib?)
106  const unsigned int I = fLib->NEvents();
107  for(unsigned int i = 0; i < I; ++i){
108  const MatchableEvent& libEvt = fLib->Event(i);
109 
110  if(enrich != 2){
111  if(libEvt.enrich != bool(enrich)) continue;
112  }
113 
114  // There's beam nue in the library, but we don't want to match to it
115  if(libEvt.weight == 0) continue;
116 
117  for(int flipEven: {0, 1}){
118  for(int flipOdd: {0, 1}){
119 
120  const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial, libEvt);
121  Elist.emplace(E, &libEvt, flipEven, flipOdd, i);
122  } // end for flipOdd
123  } // end for flipEven
124  } // end for i
125 
126  return PackageMatches(Elist, trial);
127  }
128 
129  //......................................................................
131  const MatchableEvent& trial,
132  bool& bestFlipEven, bool& bestFlipOdd,
133  double& bestEnergy) const
134  {
135  int bestHead = -1;
136  bestEnergy = FLT_MAX;
137 
138  const int nHeads = fHeads->NHeads();
139  for(int iHead = 0; iHead < nHeads; ++iHead){
140  const MatchableEvent& headEvt = fLib->Event(fHeads->HeadIdx(iHead));
141 
142  for(int flipEven: {0, 1}){
143  for(int flipOdd: {0, 1}){
144  const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial,
145  headEvt);
146  if(E < bestEnergy){
147  bestEnergy = E;
148  bestHead = iHead;
149  bestFlipEven = flipEven;
150  bestFlipOdd = flipOdd;
151  }
152  }
153  }
154  }
155  assert(bestHead >= 0);
156 
157  return bestHead;
158  }
159 
160  //......................................................................
161  std::vector<Match> FindMatchesAlg::
163  unsigned int numMatches,
164  int enrich) const
165  {
166  FlippedPotentials trialVs;
167  FillPotential(trial, trialVs);
168 
169  bool headFlipEven, headFlipOdd;
170  double headE;
171  const int headIdx = BestHeadFor(trialVs, trial,
172  headFlipEven, headFlipOdd, headE);
173  // const MatchableEvent& head = fLib->Event(fHeads->HeadIdx(headIdx));
174 
175  // Candidates for return, remains sorted
176  SortedMatchList Elist(numMatches);
177 
178  // This is the quickest ordering of the loops (locality into lib?)
179  const unsigned int J = fHeads->HeadSeqsLen();
180  for(unsigned int j = 0; j < J; ++j){
181  const int i = fHeads->ChildIdxFor(headIdx, j);
182 
183  const MatchableEvent& libEvt = fLib->Event(i);
184 
185  if(enrich != 2){
186  if(libEvt.enrich != bool(enrich)) continue;
187  }
188 
189  // There's beam nue in the library, but we don't want to match to it
190  if(libEvt.weight == 0) continue;
191 
192  const bool flipEven = (fHeads->FlipEvenFor(headIdx, j) ^ headFlipEven);
193  const bool flipOdd = (fHeads->FlipOddFor (headIdx, j) ^ headFlipOdd );
194 
195  const float E = CalcEnergy(trialVs.V[flipEven][flipOdd], trial, libEvt);
196  Elist.emplace(E, &libEvt, flipEven, flipOdd, j);
197  } // end for i/j
198 
199  return PackageMatches(Elist, trial);
200  }
201 } // namespace
bool enrich
Was this event generated in a filtered-for-pi0 job?
Definition: EventSummary.h:73
unsigned int fNumMatches
float CalcEnergy(const Potential &Va, const MatchableEvent &a, const MatchableEvent &b)
Definition: FindMatches.cxx:21
Attach some information used in matching to an EventSummary.
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
std::vector< Match > FindMatches(const MatchableEvent &trial, unsigned int numMatches, int enrich) const
FindMatchesAlg(const std::string &libDir, bool preload, bool useHeads)
std::vector< Match > PackageMatches(const std::multiset< Match > &matches, const MatchableEvent &trial) const
PID
Definition: FillPIDs.h:14
Core of the LEM match-finding algorithm.
Potential V[2][2]
Definition: FindMatches.h:23
void emplace(double E, T...args)
Float_t E
Definition: plot.C:20
double FracChargeMatched(const EventSummary &a, const EventSummary &b, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:37
SortedMatchList(unsigned int numMatches)
const double j
Definition: BetheBloch.cxx:29
Information about a LEM match.
Definition: Match.h:17
void FillPotential(const EventSummary &trial, Potential &V, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:76
int BestHeadFor(const FlippedPotentials &trialVs, const MatchableEvent &trial, bool &bestFlipEven, bool &bestFlipOdd, double &bestEnergy) const
assert(nhit_max >=nhit_nbins)
Collection of Potential objects with odd and/or even view flipped.
Definition: FindMatches.h:21
static Heads * FromMMap(const std::string &libPath, bool touchAll)
Definition: Heads.cxx:23
void insert(Match m)
std::vector< Match > FindMatchesHeads(const MatchableEvent &trial, unsigned int numMatches, int enrich) const
double T
Definition: Xdiff_gwt.C:5
static Library * FromMMap(const std::string &libPath, bool touchAll)
Definition: Library.cxx:62
std::vector< Match > FindMatchesNoHeads(const MatchableEvent &trial, unsigned int numMatches, int enrich) const
enum BeamMode string