ConformalMappingFit.h

Go to the documentation of this file.

#ifndef ConformalMappingFit_H
#define ConformalMappingFit_H

#include "DataFormats/GeometryVector/interface/Basic2DVector.h"
#include "DataFormats/GeometryVector/interface/Basic3DVector.h"
#include "DataFormats/GeometrySurface/interface/TkRotation.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"

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

namespace edm {
  class ParameterSet;
}

class ConformalMappingFit {
public:
  typedef TkRotation<double> Rotation;
  typedef Basic2DVector<double> PointXY;

  ConformalMappingFit(const std::vector<PointXY> &hits,
                      const std::vector<float> &errRPhi2,
                      const Rotation *rotation = nullptr);

  ~ConformalMappingFit();

  Measurement1D curvature() const;
  Measurement1D directionPhi() const;
  Measurement1D impactParameter() const;

  int charge() const;
  double chi2() const { return theFit.chi2(); }

  const Rotation *rotation() const { return theRotation; }

  void fixImpactParmaeter(double ip) { theFit.fixParC(ip); }
  void skipErrorCalculation() { theFit.skipErrorCalculationByDefault(); }

private:
  double phiRot() const;
  void findRot(const PointXY &);

private:
  const Rotation *theRotation;
  bool myRotation;
  ParabolaFit theFit;

  template <class T>
  class MappedPoint {
  public:
    typedef Basic2DVector<T> PointXY;
    MappedPoint() : theU(0), theV(0), theW(0), pRot(0) {}
    MappedPoint(const T &aU, const T &aV, const T &aWeight, const TkRotation<T> *aRot)
        : theU(aU), theV(aV), theW(aWeight), pRot(aRot) {}
    MappedPoint(const PointXY &point, const T &weight, const TkRotation<T> *aRot) : pRot(aRot) {
      T radius2 = point.mag2();
      Basic3DVector<T> rotated = (*pRot) * point;
      theU = rotated.x() / radius2;
      theV = rotated.y() / radius2;
      theW = weight * radius2 * radius2;
    }
    T u() const { return theU; }
    T v() const { return theV; }
    T weight() const { return theW; }
    PointXY unmap() const {
      T invRadius2 = theU * theU + theV * theV;
      Basic3DVector<T> tmp = (*pRot).multiplyInverse(Basic2DVector<T>(theU, theV));
      return PointXY(tmp.x() / invRadius2, tmp.y() / invRadius2);
    }
    T unmappedWeight() const {
      T invRadius2 = theU * theU + theV * theV;
      return theW * invRadius2 * invRadius2;
    }

  private:
    T theU, theV, theW;
    const TkRotation<T> *pRot;
  };
};

#endif