LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Solver.cxx File Reference
#include "larreco/SpacePointSolver/Solver.h"
#include <algorithm>
#include <cstdlib>
#include <iostream>
#include <set>
#include <string>

Go to the source code of this file.

Functions

template<class T >
sqr (T x)
 
double Metric (double q, double p)
 
QuadExpr Metric (double q, QuadExpr p)
 
double Metric (const std::vector< SpaceCharge * > &scs, double alpha)
 
double Metric (const std::vector< CollectionWireHit * > &cwires, double alpha)
 
QuadExpr Metric (const SpaceCharge *sci, const SpaceCharge *scj, double alpha)
 
QuadExpr Metric (const SpaceCharge *sc, double alpha)
 
double SolvePair (SpaceCharge *sci, SpaceCharge *scj, double alpha)
 
void Iterate (CollectionWireHit *cwire, double alpha)
 
void Iterate (SpaceCharge *sc, double alpha)
 
void Iterate (const std::vector< CollectionWireHit * > &cwires, const std::vector< SpaceCharge * > &orphanSCs, double alpha)
 

Function Documentation

void Iterate ( CollectionWireHit cwire,
double  alpha 
)

Definition at line 260 of file Solver.cxx.

References SpaceCharge::AddCharge(), CollectionWireHit::fCrossings, SolvePair(), and x.

Referenced by Iterate(), and reco3d::SpacePointSolver::Minimize().

261 {
262  // Consider all pairs of crossings
263  const unsigned int N = cwire->fCrossings.size();
264 
265  for (unsigned int i = 0; i + 1 < N; ++i) {
266  SpaceCharge* sci = cwire->fCrossings[i];
267 
268  for (unsigned int j = i + 1; j < N; ++j) {
269  SpaceCharge* scj = cwire->fCrossings[j];
270 
271  const double x = SolvePair(sci, scj, alpha);
272 
273  if (x == 0) continue;
274 
275  // Actually make the update
276  sci->AddCharge(+x);
277  scj->AddCharge(-x);
278  } // end for j
279  } // end for i
280 }
Float_t x
Definition: compare.C:6
void AddCharge(double dq)
Definition: Solver.cxx:32
double SolvePair(SpaceCharge *sci, SpaceCharge *scj, double alpha)
Definition: Solver.cxx:215
std::vector< SpaceCharge * > fCrossings
Definition: Solver.h:71
void Iterate ( SpaceCharge sc,
double  alpha 
)

Definition at line 283 of file Solver.cxx.

References SpaceCharge::AddCharge(), QuadExpr::Eval(), SpaceCharge::fPred, QuadExpr::Linear(), Metric(), QuadExpr::Quadratic(), and x.

284 {
285  const QuadExpr chisq = Metric(sc, alpha);
286 
287  // Find the minimum of a quadratic expression
288  double x = -chisq.Linear() / (2 * chisq.Quadratic());
289 
290  // Don't allow the SpaceCharge to go negative
291  const double xmin = -sc->fPred;
292 
293  // Clamp to allowed range
294  x = std::max(xmin, x);
295 
296  const double chisq_new = chisq.Eval(x);
297 
298  // Should try here too, because the function might be convex not concave, so
299  // d/dx=0 gives the max not the min, and the true min is at one extreme of
300  // the range.
301  const double chisq_n = chisq.Eval(xmin);
302 
303  if (chisq_n < chisq_new)
304  sc->AddCharge(xmin);
305  else
306  sc->AddCharge(x);
307 }
Float_t x
Definition: compare.C:6
void AddCharge(double dq)
Definition: Solver.cxx:32
double Linear() const
Definition: QuadExpr.h:15
double Metric(double q, double p)
Definition: Solver.cxx:70
double Eval(double x) const
Definition: QuadExpr.cxx:67
double fPred
Definition: Solver.h:57
double Quadratic() const
Definition: QuadExpr.h:14
void Iterate ( const std::vector< CollectionWireHit * > &  cwires,
const std::vector< SpaceCharge * > &  orphanSCs,
double  alpha 
)

Definition at line 310 of file Solver.cxx.

References Iterate(), and sc.

313 {
314  // Visiting in a "random" order helps prevent local artefacts that are slow
315  // to break up.
316  unsigned int cwireIdx = 0;
317  if (!cwires.empty()) {
318  do {
319  Iterate(cwires[cwireIdx], alpha);
320 
321  const unsigned int prime = 1299827;
322  cwireIdx = (cwireIdx + prime) % cwires.size();
323  } while (cwireIdx != 0);
324  }
325 
326  for (SpaceCharge* sc : orphanSCs)
327  Iterate(sc, alpha);
328 }
void Iterate(CollectionWireHit *cwire, double alpha)
Definition: Solver.cxx:260
Float_t sc
Definition: plot.C:23
double Metric ( double  q,
double  p 
)

Definition at line 70 of file Solver.cxx.

References sqr().

Referenced by Iterate(), Metric(), reco3d::SpacePointSolver::Minimize(), and SolvePair().

71 {
72  return sqr(q - p);
73 }
T sqr(T x)
Definition: Solver.cxx:10
QuadExpr Metric ( double  q,
QuadExpr  p 
)

Definition at line 76 of file Solver.cxx.

References sqr().

77 {
78  return sqr(q - p);
79 }
T sqr(T x)
Definition: Solver.cxx:10
double Metric ( const std::vector< SpaceCharge * > &  scs,
double  alpha 
)

Definition at line 82 of file Solver.cxx.

References Metric(), sc, and sqr().

83 {
84  double ret = 0;
85 
86  std::set<InductionWireHit*> iwires;
87  for (const SpaceCharge* sc : scs) {
88  if (sc->fWire1) iwires.insert(sc->fWire1);
89  if (sc->fWire2) iwires.insert(sc->fWire2);
90 
91  if (alpha != 0) {
92  ret -= alpha * sqr(sc->fPred);
93  // "Double-counting" of the two ends of the connection is
94  // intentional. Otherwise we'd have a half in the line above.
95  ret -= alpha * sc->fPred * sc->fNeiPotential;
96  }
97  }
98 
99  for (const InductionWireHit* iwire : iwires) {
100  ret += Metric(iwire->fCharge, iwire->fPred);
101  }
102 
103  return ret;
104 }
double Metric(double q, double p)
Definition: Solver.cxx:70
T sqr(T x)
Definition: Solver.cxx:10
Float_t sc
Definition: plot.C:23
double Metric ( const std::vector< CollectionWireHit * > &  cwires,
double  alpha 
)

Definition at line 107 of file Solver.cxx.

References Metric().

108 {
109  std::vector<SpaceCharge*> scs;
110  for (CollectionWireHit* cwire : cwires)
111  scs.insert(scs.end(), cwire->fCrossings.begin(), cwire->fCrossings.end());
112  return Metric(scs, alpha);
113 }
double Metric(double q, double p)
Definition: Solver.cxx:70
QuadExpr Metric ( const SpaceCharge sci,
const SpaceCharge scj,
double  alpha 
)

Definition at line 116 of file Solver.cxx.

References InductionWireHit::fCharge, Neighbour::fCoupling, SpaceCharge::fNeighbours, SpaceCharge::fNeiPotential, InductionWireHit::fPred, SpaceCharge::fPred, Neighbour::fSC, SpaceCharge::fWire1, SpaceCharge::fWire2, Metric(), n, sqr(), x, and QuadExpr::X().

117 {
118  QuadExpr ret = 0;
119 
120  // How much charge moves from scj to sci
121  QuadExpr x = QuadExpr::X();
122 
123  if (alpha != 0) {
124  const double scip = sci->fPred;
125  const double scjp = scj->fPred;
126 
127  // Self energy. SpaceCharges are never the same object
128  ret -= alpha * sqr(scip + x);
129  ret -= alpha * sqr(scjp - x);
130 
131  // Interaction. We're only seeing one end of the double-ended connection
132  // here, so multiply by two.
133  ret -= 2 * alpha * (scip + x) * sci->fNeiPotential;
134  ret -= 2 * alpha * (scjp - x) * scj->fNeiPotential;
135 
136  // This miscounts if i and j are neighbours of each other
137  for (const Neighbour& n : sci->fNeighbours) {
138  if (n.fSC == scj) {
139  // If we detect that case, remove the erroneous terms
140  ret += 2 * alpha * (scip + x) * scjp * n.fCoupling;
141  ret += 2 * alpha * (scjp - x) * scip * n.fCoupling;
142 
143  // And replace with the correct interaction terms
144  ret -= 2 * alpha * (scip + x) * (scjp - x) * n.fCoupling;
145  break;
146  }
147  }
148  }
149 
150  const InductionWireHit* iwire1 = sci->fWire1;
151  const InductionWireHit* jwire1 = scj->fWire1;
152 
153  const double qi1 = iwire1 ? iwire1->fCharge : 0;
154  const double pi1 = iwire1 ? iwire1->fPred : 0;
155 
156  const double qj1 = jwire1 ? jwire1->fCharge : 0;
157  const double pj1 = jwire1 ? jwire1->fPred : 0;
158 
159  if (iwire1 == jwire1) {
160  // Same wire means movement of charge cancels itself out
161  if (iwire1) ret += Metric(qi1, pi1);
162  }
163  else {
164  if (iwire1) ret += Metric(qi1, pi1 + x);
165  if (jwire1) ret += Metric(qj1, pj1 - x);
166  }
167 
168  const InductionWireHit* iwire2 = sci->fWire2;
169  const InductionWireHit* jwire2 = scj->fWire2;
170 
171  const double qi2 = iwire2 ? iwire2->fCharge : 0;
172  const double pi2 = iwire2 ? iwire2->fPred : 0;
173 
174  const double qj2 = jwire2 ? jwire2->fCharge : 0;
175  const double pj2 = jwire2 ? jwire2->fPred : 0;
176 
177  if (iwire2 == jwire2) {
178  if (iwire2) ret += Metric(qi2, pi2);
179  }
180  else {
181  if (iwire2) ret += Metric(qi2, pi2 + x);
182  if (jwire2) ret += Metric(qj2, pj2 - x);
183  }
184 
185  return ret;
186 }
Float_t x
Definition: compare.C:6
double fCharge
Definition: Solver.h:23
double fPred
Definition: Solver.h:25
double fCoupling
Definition: Solver.h:34
std::vector< Neighbour > fNeighbours
Definition: Solver.h:55
double fNeiPotential
Neighbour-induced potential.
Definition: Solver.h:58
InductionWireHit * fWire2
Definition: Solver.h:53
static QuadExpr X()
Definition: QuadExpr.cxx:8
InductionWireHit * fWire1
Definition: Solver.h:53
double Metric(double q, double p)
Definition: Solver.cxx:70
T sqr(T x)
Definition: Solver.cxx:10
Char_t n[5]
double fPred
Definition: Solver.h:57
SpaceCharge * fSC
Definition: Solver.h:33
QuadExpr Metric ( const SpaceCharge sc,
double  alpha 
)

Definition at line 189 of file Solver.cxx.

References InductionWireHit::fCharge, SpaceCharge::fNeiPotential, InductionWireHit::fPred, SpaceCharge::fPred, SpaceCharge::fWire1, SpaceCharge::fWire2, Metric(), sqr(), x, and QuadExpr::X().

190 {
191  QuadExpr ret = 0;
192 
193  // How much charge is added to sc
194  QuadExpr x = QuadExpr::X();
195 
196  if (alpha != 0) {
197  const double scp = sc->fPred;
198 
199  // Self energy
200  ret -= alpha * sqr(scp + x);
201 
202  // Interaction. We're only seeing one end of the double-ended connection
203  // here, so multiply by two.
204  ret -= 2 * alpha * (scp + x) * sc->fNeiPotential;
205  }
206 
207  // Prediction of the induction wires
208  ret += Metric(sc->fWire1->fCharge, sc->fWire1->fPred + x);
209  ret += Metric(sc->fWire2->fCharge, sc->fWire2->fPred + x);
210 
211  return ret;
212 }
Float_t x
Definition: compare.C:6
double fCharge
Definition: Solver.h:23
double fPred
Definition: Solver.h:25
double fNeiPotential
Neighbour-induced potential.
Definition: Solver.h:58
InductionWireHit * fWire2
Definition: Solver.h:53
static QuadExpr X()
Definition: QuadExpr.cxx:8
InductionWireHit * fWire1
Definition: Solver.h:53
double Metric(double q, double p)
Definition: Solver.cxx:70
T sqr(T x)
Definition: Solver.cxx:10
double fPred
Definition: Solver.h:57
double SolvePair ( SpaceCharge sci,
SpaceCharge scj,
double  alpha 
)

Definition at line 215 of file Solver.cxx.

References QuadExpr::Eval(), SpaceCharge::fPred, QuadExpr::Linear(), Metric(), QuadExpr::Quadratic(), and x.

Referenced by Iterate().

216 {
217  const QuadExpr chisq = Metric(sci, scj, alpha);
218  const double chisq0 = chisq.Eval(0);
219 
220  // Find the minimum of a quadratic expression
221  double x = -chisq.Linear() / (2 * chisq.Quadratic());
222 
223  // Don't allow either SpaceCharge to go negative
224  const double xmin = -sci->fPred;
225  const double xmax = scj->fPred;
226 
227  // Clamp to allowed range
228  x = std::min(xmax, x);
229  x = std::max(xmin, x);
230 
231  const double chisq_new = chisq.Eval(x);
232 
233  // Should try these too, because the function might be convex not concave, so
234  // d/dx=0 gives the max not the min, and the true min is at one extreme of
235  // the range.
236  const double chisq_p = chisq.Eval(xmax);
237  const double chisq_n = chisq.Eval(xmin);
238 
239  if (std::min(std::min(chisq_p, chisq_n), chisq_new) > chisq0 + 1) {
240  std::cout << "Solution at " << x << " is worse than current state! Scan from " << xmin << " to "
241  << xmax << std::endl;
242  for (double x = xmin; x < xmax; x += .01 * (xmax - xmin)) {
243  std::cout << x << " " << chisq.Eval(x) << std::endl;
244  }
245 
246  std::cout << "Soln, original, up edge, low edge:" << std::endl;
247  std::cout << chisq_new << " " << chisq0 << " " << chisq_p << " " << chisq_n << std::endl;
248  abort();
249  }
250 
251  if (std::min(chisq_n, chisq_p) < chisq_new) {
252  if (chisq_n < chisq_p) return xmin;
253  return xmax;
254  }
255 
256  return x;
257 }
Float_t x
Definition: compare.C:6
double Linear() const
Definition: QuadExpr.h:15
double Metric(double q, double p)
Definition: Solver.cxx:70
double Eval(double x) const
Definition: QuadExpr.cxx:67
double fPred
Definition: Solver.h:57
double Quadratic() const
Definition: QuadExpr.h:14
template<class T >
T sqr ( x)

Definition at line 10 of file Solver.cxx.

References x.

Referenced by Metric().

11 {
12  return x * x;
13 }
Float_t x
Definition: compare.C:6