Solver.h File Reference

#include <vector>
#include "QuadExpr.h"

Go to the source code of this file.

Classes
class	WireHit
class	InductionWireHit
class	Neighbour
class	SpaceCharge
class	CollectionWireHit

Functions
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 xmin, double xmax, 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().

{
  // Consider all pairs of crossings
  const unsigned int N = cwire->fCrossings.size();

  for(unsigned int i = 0; i+1 < N; ++i){
    SpaceCharge* sci = cwire->fCrossings[i];

    InductionWireHit* iwire1 = sci->fWire1;
    InductionWireHit* iwire2 = sci->fWire2;

    for(unsigned int j = i+1; j < N; ++j){
      SpaceCharge* scj = cwire->fCrossings[j];

      InductionWireHit* jwire1 = scj->fWire1;
      InductionWireHit* jwire2 = scj->fWire2;

      // There can't be any cross-overs of U and V views like this
      if(iwire1 == jwire2 || iwire2 == jwire1) abort();

      // Driving all the same wires, no move can have any effect
      if(iwire1 == jwire1 && iwire2 == jwire2) continue;

      const double x = SolvePair(cwire, sci, scj, alpha);

      if(x == 0) continue;

      // Actually make the update
      sci->AddCharge(+x);
      scj->AddCharge(-x);
    } // end for j
  } // end for i
}

void Iterate	(	const std::vector< CollectionWireHit * > &	cwires,
		double	alpha
	)

Definition at line 270 of file Solver.cxx.

References Iterate().

{
  for(CollectionWireHit* cwire: cwires) Iterate(cwire, alpha);
}

double Metric	(	const std::vector< SpaceCharge * > &	scs,
		double	alpha
	)

Definition at line 82 of file Solver.cxx.

References Metric(), and sqr().

{
  double ret = 0;

  std::set<InductionWireHit*> iwires;
  for(const SpaceCharge* sc: scs){
    if(sc->fWire1) iwires.insert(sc->fWire1);
    if(sc->fWire2) iwires.insert(sc->fWire2);

    if(alpha != 0){
      ret -= alpha*sqr(sc->fPred);
      // "Double-counting" of the two ends of the connection is
      // intentional. Otherwise we'd have a half in the line above.
      ret -= alpha * sc->fPred * sc->fNeiPotential;
    }
  }

  for(const InductionWireHit* iwire: iwires){
    ret += Metric(iwire->fCharge, iwire->fPred);
  }

  return ret;
}

double Metric	(	const std::vector< CollectionWireHit * > &	cwires,
		double	alpha
	)

Definition at line 107 of file Solver.cxx.

References Metric().

{
  std::vector<SpaceCharge*> scs;
  for(CollectionWireHit* cwire: cwires)
    scs.insert(scs.end(), cwire->fCrossings.begin(), cwire->fCrossings.end());
  return Metric(scs, alpha);
}

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().

{
  QuadExpr ret = 0;

  // How much charge moves from scj to sci
  QuadExpr x = QuadExpr::X();

  if(alpha != 0){
    const double scip = sci->fPred;
    const double scjp = scj->fPred;

    // Self energy. SpaceCharges are never the same object
    ret -= alpha*sqr(scip + x);
    ret -= alpha*sqr(scjp - x);

    // Interaction. We're only seeing one end of the double-ended connection
    // here, so multiply by two.
    ret -= 2 * alpha * (scip + x) * sci->fNeiPotential;
    ret -= 2 * alpha * (scjp - x) * scj->fNeiPotential;

    // This miscounts if i and j are neighbours of each other
    for(const Neighbour& n: sci->fNeighbours){
      if(n.fSC == scj){
        // If we detect that case, remove the erroneous terms
        ret += 2 * alpha * (scip + x) * scjp * n.fCoupling;
        ret += 2 * alpha * (scjp - x) * scip * n.fCoupling;

        // And replace with the correct interaction terms
        ret -= 2 * alpha * (scip + x) * (scjp - x) * n.fCoupling;
        break;
      }
    }
  }


  const InductionWireHit* iwire1 = sci->fWire1;
  const InductionWireHit* jwire1 = scj->fWire1;
  if(iwire1 && jwire1){
    const double qi1 = iwire1->fCharge;
    const double pi1 = iwire1->fPred;

    const double qj1 = jwire1->fCharge;
    const double pj1 = jwire1->fPred;

    if(iwire1 == jwire1){
      ret += Metric(qi1, pi1);
    }
    else{
      ret += Metric(qi1, pi1 + x);
      ret += Metric(qj1, pj1 - x);
    }
  }

  const InductionWireHit* iwire2 = sci->fWire2;
  const InductionWireHit* jwire2 = scj->fWire2;
  if(iwire2 && jwire2){
    const double qi2 = iwire2->fCharge;
    const double pi2 = iwire2->fPred;

    const double qj2 = jwire2->fCharge;
    const double pj2 = jwire2->fPred;

    if(iwire2 == jwire2){
      ret += Metric(qi2, pi2);
    }
    else{
      ret += Metric(qi2, pi2 + x);
      ret += Metric(qj2, pj2 - x);
    }
  }

  return ret;
}

double SolvePair	(	CollectionWireHit *	cwire,
		SpaceCharge *	sci,
		SpaceCharge *	scj,
		double	xmin,
		double	xmax,
		double	alpha
	)