FindMatchesTranspose_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file FindMatchesTranspose.cxx
3 // \brief TODO
4 // \author Christopher Backhouse - bckhouse@caltech.edu
5 ////////////////////////////////////////////////////////////////////////
6 
7 // Framework includes
14 #include "fhiclcpp/ParameterSet.h"
16 
17 #include "LEM/LEMInput.h"
18 #include "LEM/Util.h"
19 
20 #include "LEM/func/DistanceMap.h"
21 #include "LEM/func/EventSummary.h"
22 #include "LEM/func/FindMatches.h"
23 #include "LEM/func/Heads.h"
24 #include "LEM/func/Library.h"
25 #include "LEM/func/MatchList.h"
26 #include "LEM/func/OnDemand.h"
27 
30 #include "Utilities/func/MathUtil.h"
31 
32 #include <iostream>
33 
34 namespace lem
35 {
37  {
38  public:
39  explicit FindMatchesTranspose(const fhicl::ParameterSet& pset);
41 
42  virtual void beginJob();
43 
44  virtual void reconfigure(const fhicl::ParameterSet& pset);
45 
46  virtual void produce(art::Event& evt);
47 
48  virtual void endRun(art::Run& run);
49 
50  protected:
53  bool fOnDemand;
54 
56 
57  std::vector<MatchableEvent> fTrials;
58  std::vector<lem::MatchableEvent> fTrialsDownsample;
59  std::vector<Potential> fVs;
60  std::vector<Potential> fVsDownsample;
61  };
62 
63  //.......................................................................
65  {
66  reconfigure(pset);
67 
68  produces<std::vector<geo::OfflineChan>>();
69 
70  produces<lem::LibrarySummary, art::InRun>();
71 
72  produces<std::vector<lem::MatchList>, art::InRun>();
73  }
74 
75  //......................................................................
77  {
78  }
79 
80  //......................................................................
82  {
83  std::string libpath = pset.get<std::string>("LibraryPath");
84  fLibPath = util::EnvExpansion(libpath);
85 
86  fDownsampleFactor = pset.get<int>("DownsampleFactor");
87 
88  fInputLabel = pset.get<std::string>("InputLabel");
89 
90  // This needs to be before anyone asks the DistanceMap anything, should be
91  // safe enough here.
92  DistanceMap::SetPlaneScale(pset.get<double>("DistPlaneScale"));
93  DistanceMap::SetCellScale(pset.get<double>("DistCellScale"));
94  DistanceMap::SetDecayPower(pset.get<double>("DistDecayPower"));
95  DistanceMap::SetExpMode(pset.get<bool>("DistExpMode"));
96 
97  EventSummary::fgDecayRate = pset.get<double>("HitDecayRate");
98  EventSummary::fgChargePower = pset.get<double>("HitChargePower");
99 
100  SetOnDemandMemoryLimit(pset.get<int>("LibraryMemoryLimit"));
101  fOnDemand = (pset.get<int>("LibraryMemoryLimit") > 0);
102  }
103 
104  //......................................................................
106  {
107  }
108 
109  //......................................................................
111  {
112  std::unique_ptr<std::vector<geo::OfflineChan> > chancol(new std::vector<geo::OfflineChan>);
113 
115  evt.getByLabel(fInputLabel, inputs);
116 
117  const int inputMax = inputs->size();
118  for(int inputIdx = 0; inputIdx < inputMax; ++inputIdx){
119  const LEMInput& input = (*inputs)[inputIdx];
120 
121  EventSummary sum(input.hits, -1, std::vector<int>());
122 
123  sum.run = evt.run();
124  sum.subrun = evt.subRun();
125  sum.event = evt.event();
126 
127  sum.origPdg = 0;
128 
129  sum.pdg = 0;
130  sum.ccnc = 0;
131  sum.y = 0;
132  sum.trueEVis = 0;
133 
134 
135  // chancol->push_back(geo::OfflineChan(vtxPlane, vtxCell));
136  // Should really be nextPlane, but that's lost inside of Util.cxx. This
137  // is almost always right, and only for display purposes anyway.
138  // chancol->push_back(geo::OfflineChan(vtxPlane+1, vtxCellOther));
139 
140  const MatchableEvent evt(sum);
141  fTrials.push_back(evt);
143 
144  Potential V;
145  FillPotential(sum, V, false, false);
146  fVs.push_back(V);
147 
148  Potential V_down;
149  V_down = V.Downsampled(fDownsampleFactor);
150  fVsDownsample.push_back(V_down);
151 
152  } // end for sliceIdx
153 
154  evt.put(std::move(chancol));
155  }
156 
157  // TODO: don't copy and paste with LoadLibrary
158  //......................................................................
160  {
161  if(evt.IsSig()){
162  // Scale up to keep as much signal as possible in sample, weighted
163  // correctly internally. The missing prefactor is (ss2th13*ssth23)
164  return util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
165  }
166  if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 14){
167  // CC mu
168  return 1-util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
169  }
170  if(evt.ccnc == 0 && abs(evt.pdg) == 12 && abs(evt.origPdg) == 12){
171  // Beam nue
172  // return 1-.1*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
173  // Want library matching to ignore nues
174  return 0;
175  }
176  if(evt.ccnc == 1 && evt.pdg == evt.origPdg){
177  // NC
178  return 1;
179  }
180  if(evt.ccnc == 1 && evt.pdg != evt.origPdg){
181  // Swapped NC. We play clever reweighting games so this is OK
182  return 1;
183  }
184  if(evt.ccnc == 0 && abs(evt.pdg) == 14 && abs(evt.origPdg) == 12){
185  // CC mu appearance, take ss2th13=0.1, maximal th23
186  return .1*.5*util::sqr(sin(1.267*2.35e-3*810/evt.trueEVis));
187  }
188  // Don't know what this is
189  return 0;
190  }
191 
192  //......................................................................
194  {
195  // One thing that can happen is that there's an exception thrown before
196  // anything's even happened. But endRun gets called first. Don't want to
197  // stall confusingly in that case.
198  if(fTrials.empty()) return;
199 
200  const unsigned int Ntrials = fTrials.size();
201  std::cout << "Have " << Ntrials << " trial events" << std::endl;
202 
204  std::cout << "Loaded " << heads << std::endl;
206  std::cout << "Loaded " << lib->Summary() << std::endl;
207 
208  std::cout << "Matching heads..." << std::endl;
209 
210  std::vector<int> bestHeads(Ntrials, -1);
211  std::vector<float> bestHeadEs(Ntrials, FLT_MAX);
212 
213  const int nHeads = heads->NHeads();
214  for(int i = 0; i < nHeads; ++i){
215  if(i%100 == 0) std::cout << i/100 << "/" << nHeads/100 << std::endl;
216 
217  const int headIdx = heads->HeadIdx(i);
218 
219  for(unsigned int trialIdx = 0; trialIdx < Ntrials; ++trialIdx){
220  const MatchableEvent& trial = fTrials[trialIdx];
221  const MatchableEvent& trialDown = fTrialsDownsample[trialIdx];
222 
223  for(bool flipEven: {false, true}){
224  for(bool flipOdd: {false, true}){
225 
226  const MatchableEvent headevt_down = lib->DownsampledEvent(headIdx).Flipped(flipEven, flipOdd);
227 
228  // Compare the downsampled event to the downsampled potential. The
229  // resulting energy is a lower-bound on what the real match would
230  // have done.
231  float E = CalcEnergy(fVsDownsample[trialIdx],
232  trialDown, headevt_down);
233 
234  if(E < bestHeadEs[trialIdx]){
235  const MatchableEvent headevt = lib->Event(headIdx).Flipped(flipEven, flipOdd);
236 
237  // Calculate the energy fully
238  E = CalcEnergy(fVs[trialIdx], trial, headevt);
239 
240  // If it's still good, use it
241  if(E < bestHeadEs[trialIdx]){
242  bestHeads[trialIdx] = i;
243  bestHeadEs[trialIdx] = E;
244  }
245  }
246  } // end for flipOdd
247  } // end for flipEven
248  } // end for trialIdx
249  } // end for i
250 
251  for(unsigned int trialIdx = 0; trialIdx < Ntrials; ++trialIdx){
252  std::cout << "Best head " << bestHeads[trialIdx] << " with " << bestHeadEs[trialIdx]/fTrials[trialIdx].selfEnergy << std::endl;
253  }
254 
255  std::cout << "Main PID..." << std::endl;
256 
257  const unsigned int kNumMatches = 1001;
258 
259  // Candidates for storing, remains sorted
260  std::vector<std::multiset<MatchSummary>> Elists(Ntrials);
261 
262  int* pHeadSeqs = heads->pHeadSeqs();
263 
264  const unsigned int libSize = lib->NEvents();
265  for(unsigned int libIdx = 0; libIdx < libSize; ++libIdx){
266  if(libIdx%10000 == 0) std::cout << libIdx/10000 << "/" << libSize/10000 << std::endl;
267 
268  const int numHeadsForLibEvt = *pHeadSeqs;
269  ++pHeadSeqs;
270  std::vector<bool> isHeadForLibEvt(nHeads, false);
271  for(int i = 0; i < numHeadsForLibEvt; ++i){
272  isHeadForLibEvt[*pHeadSeqs] = true;
273  ++pHeadSeqs;
274  }
275 
276  for(unsigned int trialIdx = 0; trialIdx < Ntrials; ++trialIdx){
277  // Don't try to match this library event and this trial if none of this
278  // library events heads are the one the trial event wants.
279  if(!isHeadForLibEvt[bestHeads[trialIdx]]) continue;
280 
281  const MatchableEvent& trial = fTrials[trialIdx];
282  const MatchableEvent& trialDown = fTrialsDownsample[trialIdx];
283  std::multiset<MatchSummary>& Elist = Elists[trialIdx];
284 
285  float target = (Elist.size() == kNumMatches) ? Elist.rbegin()->potential : FLT_MAX;
286 
287  for(bool flipEven: {false, true}){
288  for(bool flipOdd: {false, true}){
289  const lem::MatchableEvent libevt_down = lib->DownsampledEvent(libIdx).Flipped(flipEven, flipOdd);
290 
291  // Compare the downsampled event to the downsampled potential. The
292  // resulting energy is a lower-bound on what the real match would
293  // have done.
294  float E = CalcEnergy(fVsDownsample[trialIdx],
295  trialDown, libevt_down);
296 
297  if(E < target){
298  const lem::MatchableEvent libevt = lib->Event(libIdx).Flipped(flipEven, flipOdd);
299  // Calculate the energy fully
300  E = CalcEnergy(fVs[trialIdx], trial, libevt);
301 
302  // If it's still good, use it
303  if(E < target){
304  Match m(E, &libevt, flipEven, flipOdd);
305  // This is extra work for matches we're not even sure we're
306  // going to keep. But it's better than having to do a whole
307  // extra pass through the library at the end.
308  m.qfrac = FracChargeMatched(trial, libevt,
309  flipEven, flipOdd);
310  MatchSummary ms(m);
311  Elist.insert(ms);
312 
313  if(Elist.size() > kNumMatches){
314  // Always drop the last one, to keep list correct size
315  auto itLargest = Elist.end();
316  --itLargest;
317  Elist.erase(itLargest);
318  target = Elist.rbegin()->potential;
319  }
320  }
321  }
322  } // end for flipOdd
323  } // end for flipEven
324 
325  } // end for trialIdx
326  } // end for n
327 
328  std::cout << "Writing results out..." << std::endl;
329 
330  std::unique_ptr<std::vector<lem::MatchList>> ms(new std::vector<lem::MatchList>);
331  ms->resize(Ntrials);
332 
333  for(unsigned int trialIdx = 0; trialIdx < Ntrials; ++trialIdx){
334  std::cout << trialIdx << "/" << Ntrials << std::endl;
335 
336  (*ms)[trialIdx].trial = MatchSummary(Match(0, &fTrials[trialIdx]));
337  for(MatchSummary m: Elists[trialIdx]){
338  m.potential /= fTrials[trialIdx].selfEnergy;
339  (*ms)[trialIdx].matches.push_back(m);
340  }
341  }
342 
343  run.put(std::move(ms));
344 
345  std::unique_ptr<lem::LibrarySummary> ls(new lem::LibrarySummary);
346  *ls = lib->Summary();
347  run.put(std::move(ls));
348 
349  delete lib;
350  delete heads;
351  }
352 
354 
355 } // end namespace lem
356 ////////////////////////////////////////////////////////////////////////
static constexpr Double_t ms
Definition: Munits.h:192
SubRunNumber_t subRun() const
Definition: Event.h:72
Simple representation of event for LEM use.
Definition: EventSummary.h:26
static double fgChargePower
Definition: EventSummary.h:76
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
const MatchableEvent & DownsampledEvent(int i) const
Definition: Library.h:33
Map of electrostatic potential at each cell.
Definition: Potential.h:13
Collection of events for matching.
Collection of events for matching.
Definition: Library.h:18
MatchableEvent Flipped(bool flipEven, bool flipOdd) const
std::vector< Potential > fVsDownsample
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
const LibrarySummary & Summary() const
Definition: Library.h:30
void abs(TH1 *hist)
art::ProductID put(std::unique_ptr< PROD > &&)
Definition: Run.h:149
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Mapping from a subset of the library to their best matches.
Definition: Run.h:31
Collection of MatchSummary objects.
static HeadsTranspose * FromMMapOnDemand(const std::string &libPath)
Definition: Heads.cxx:116
double SimpleSurvivalProb(const EventSummary &evt)
Definition: __init__.py:1
PID
Definition: FillPIDs.h:14
Core of the LEM match-finding algorithm.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Simplified Match information, suitable for serialization.
Definition: MatchSummary.h:16
virtual void endRun(art::Run &run)
static void SetPlaneScale(double ps)
Must call before first call to Instance()
Definition: DistanceMap.h:32
bool IsSig() const
Definition: EventSummary.h:36
string lib
Definition: cvnie.py:9
Float_t E
Definition: plot.C:20
int HeadIdx(int i) const
Definition: Heads.h:57
virtual void reconfigure(const fhicl::ParameterSet &pset)
double qfrac
Definition: Match.h:25
unsigned int NEvents() const
Definition: Library.h:31
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
double FracChargeMatched(const EventSummary &a, const EventSummary &b, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:37
static void SetDecayPower(double dp)
Definition: DistanceMap.h:34
int NHeads() const
Definition: Heads.h:56
MatchableEvent Downsampled(int factor) const
int * pHeadSeqs() const
Definition: Heads.h:59
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.
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
Utilities for mapping in library files on-demand.
std::vector< MatchableEvent > fTrials
void FillPotential(const EventSummary &trial, Potential &V, bool flipEven, bool flipOdd)
Definition: FindMatches.cxx:76
void SetOnDemandMemoryLimit(long limit_mb)
Definition: OnDemand.cxx:117
T sin(T number)
Definition: d0nt_math.hpp:132
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
static Library * FromMMapOnDemand(const std::string &libPath)
Definition: Library.cxx:82
std::vector< LiteHit > hits
Definition: LEMInput.h:14
std::vector< lem::MatchableEvent > fTrialsDownsample
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
virtual void produce(art::Event &evt)
static HeadsTranspose * FromMMap(const std::string &libPath)
Definition: Heads.cxx:99
def ls(target="")
Definition: g4zmq.py:69
static Library * FromMMap(const std::string &libPath, bool touchAll)
Definition: Library.cxx:62
Simple representation of event for LEM use.
Potential Downsampled(int factor) const
Definition: Potential.cxx:23
Double_t sum
Definition: plot.C:31
Calculate and cache electrostatic potential between cells.
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
FindMatchesTranspose(const fhicl::ParameterSet &pset)
Simple object representing a (plane, cell) pair.
enum BeamMode string