LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
Solver.cxx File Reference
#include "Solver.h"
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <set>
#include <iostream>

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)
 
double SolvePair (CollectionWireHit *cwire, SpaceCharge *sci, SpaceCharge *scj, double alpha)
 
void Iterate (CollectionWireHit *cwire, double alpha)
 
void Iterate (const std::vector< CollectionWireHit * > &cwires, double alpha)
 

Function Documentation

void Iterate ( CollectionWireHit cwire,
double  alpha 
)

Definition at line 235 of file Solver.cxx.

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

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

236 {
237  // Consider all pairs of crossings
238  const unsigned int N = cwire->fCrossings.size();
239 
240  for(unsigned int i = 0; i+1 < N; ++i){
241  SpaceCharge* sci = cwire->fCrossings[i];
242 
243  InductionWireHit* iwire1 = sci->fWire1;
244  InductionWireHit* iwire2 = sci->fWire2;
245 
246  for(unsigned int j = i+1; j < N; ++j){
247  SpaceCharge* scj = cwire->fCrossings[j];
248 
249  InductionWireHit* jwire1 = scj->fWire1;
250  InductionWireHit* jwire2 = scj->fWire2;
251 
252  // There can't be any cross-overs of U and V views like this
253  if(iwire1 == jwire2 || iwire2 == jwire1) abort();
254 
255  // Driving all the same wires, no move can have any effect
256  if(iwire1 == jwire1 && iwire2 == jwire2) continue;
257 
258  const double x = SolvePair(cwire, sci, scj, alpha);
259 
260  if(x == 0) continue;
261 
262  // Actually make the update
263  sci->AddCharge(+x);
264  scj->AddCharge(-x);
265  } // end for j
266  } // end for i
267 }
Float_t x
Definition: compare.C:6
void AddCharge(double dq)
Definition: Solver.cxx:35
std::vector< SpaceCharge * > fCrossings
Definition: Solver.h:73
double SolvePair(CollectionWireHit *cwire, SpaceCharge *sci, SpaceCharge *scj, double alpha)
Definition: Solver.cxx:191
InductionWireHit * fWire2
Definition: Solver.h:54
InductionWireHit * fWire1
Definition: Solver.h:54
void Iterate ( const std::vector< CollectionWireHit * > &  cwires,
double  alpha 
)

Definition at line 270 of file Solver.cxx.

References Iterate().

272 {
273  for(CollectionWireHit* cwire: cwires) Iterate(cwire, alpha);
274 }
void Iterate(CollectionWireHit *cwire, double alpha)
Definition: Solver.cxx:235
double Metric ( double  q,
double  p 
)

Definition at line 70 of file Solver.cxx.

References sqr().

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

71 {
72  return sqr(q-p);
73 }
T sqr(T x)
Definition: Solver.cxx:9
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:9
double Metric ( const std::vector< SpaceCharge * > &  scs,
double  alpha 
)

Definition at line 82 of file Solver.cxx.

References Metric(), 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:9
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 
151  const InductionWireHit* iwire1 = sci->fWire1;
152  const InductionWireHit* jwire1 = scj->fWire1;
153  if(iwire1 && jwire1){
154  const double qi1 = iwire1->fCharge;
155  const double pi1 = iwire1->fPred;
156 
157  const double qj1 = jwire1->fCharge;
158  const double pj1 = jwire1->fPred;
159 
160  if(iwire1 == jwire1){
161  ret += Metric(qi1, pi1);
162  }
163  else{
164  ret += Metric(qi1, pi1 + x);
165  ret += Metric(qj1, pj1 - x);
166  }
167  }
168 
169  const InductionWireHit* iwire2 = sci->fWire2;
170  const InductionWireHit* jwire2 = scj->fWire2;
171  if(iwire2 && jwire2){
172  const double qi2 = iwire2->fCharge;
173  const double pi2 = iwire2->fPred;
174 
175  const double qj2 = jwire2->fCharge;
176  const double pj2 = jwire2->fPred;
177 
178  if(iwire2 == jwire2){
179  ret += Metric(qi2, pi2);
180  }
181  else{
182  ret += Metric(qi2, pi2 + x);
183  ret += Metric(qj2, pj2 - x);
184  }
185  }
186 
187  return ret;
188 }
Float_t x
Definition: compare.C:6
double fCharge
Definition: Solver.h:25
double fPred
Definition: Solver.h:27
double fCoupling
Definition: Solver.h:37
std::vector< Neighbour > fNeighbours
Definition: Solver.h:56
double fNeiPotential
Neighbour-induced potential.
Definition: Solver.h:59
InductionWireHit * fWire2
Definition: Solver.h:54
static QuadExpr X()
Definition: QuadExpr.cxx:6
InductionWireHit * fWire1
Definition: Solver.h:54
double Metric(double q, double p)
Definition: Solver.cxx:70
T sqr(T x)
Definition: Solver.cxx:9
Char_t n[5]
double fPred
Definition: Solver.h:58
SpaceCharge * fSC
Definition: Solver.h:36
double SolvePair ( CollectionWireHit cwire,
SpaceCharge sci,
SpaceCharge scj,
double  alpha 
)

Definition at line 191 of file Solver.cxx.

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

Referenced by Iterate().

194 {
195  const QuadExpr chisq = Metric(sci, scj, alpha);
196  const double chisq0 = chisq.Eval(0);
197 
198  // Find the minimum of a quadratic expression
199  double x = -chisq.Linear()/(2*chisq.Quadratic());
200 
201  // Don't allow either SpaceCharge to go negative
202  const double xmin = -sci->fPred;
203  const double xmax = scj->fPred;
204 
205  // Clamp to allowed range
206  x = std::min(xmax, x);
207  x = std::max(xmin, x);
208 
209  const double chisq_new = chisq.Eval(x);
210 
211  // Should try these too, because the function might be convex not concave, so
212  // d/dx=0 gives the max not the min, and the true min is at one extreme of
213  // the range.
214  const double chisq_p = chisq.Eval(xmax);
215  const double chisq_n = chisq.Eval(xmin);
216 
217  if(std::min(std::min(chisq_p, chisq_n), chisq_new) > chisq0+1){
218  for(double x = xmin; x < xmax; x += .01*(xmax-xmin)){
219  std::cout << x << " " << chisq.Eval(x) << std::endl;
220  }
221 
222  std::cout << chisq_new << " " << chisq0 << " " << chisq_p << " " << chisq_n << std::endl;
223  abort();
224  }
225 
226  if(std::min(chisq_n, chisq_p) < chisq_new){
227  if(chisq_n < chisq_p) return xmin;
228  return xmax;
229  }
230 
231  return x;
232 }
Float_t x
Definition: compare.C:6
double Linear() const
Definition: QuadExpr.h:16
Int_t max
Definition: plot.C:27
double Metric(double q, double p)
Definition: Solver.cxx:70
Int_t min
Definition: plot.C:26
double Eval(double x) const
Definition: QuadExpr.cxx:67
double fPred
Definition: Solver.h:58
double Quadratic() const
Definition: QuadExpr.h:15
template<class T >
T sqr ( x)

Definition at line 9 of file Solver.cxx.

References x.

Referenced by Metric().

9 {return x*x;}
Float_t x
Definition: compare.C:6