ConversionFastHelix.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionFastHelix.h"
#include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
#include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
#include "TrackingTools/TrajectoryParametrization/interface/CartesianTrajectoryError.h"
#include "FWCore/Utilities/interface/isFinite.h"

#include <cfloat>

ConversionFastHelix::ConversionFastHelix(const GlobalPoint& outerHit,
                                         const GlobalPoint& middleHit,
                                         const GlobalPoint& aVertex,
                                         const MagneticField* field)
    : theOuterHit(outerHit),
      theMiddleHit(middleHit),
      theVertex(aVertex),
      theCircle(outerHit, middleHit, aVertex),
      mField(field) {
  validStateAtVertex = false;

  makeHelix();
}

void ConversionFastHelix::makeHelix() {
  if (theCircle.isValid()) {
    theHelix_ = helixStateAtVertex();
  } else {
    theHelix_ = straightLineStateAtVertex();
  }
}

FreeTrajectoryState ConversionFastHelix::stateAtVertex() { return theHelix_; }

FreeTrajectoryState ConversionFastHelix::helixStateAtVertex() {
  GlobalPoint pMid(theMiddleHit);
  GlobalPoint v(theVertex);
  FTS atVertex;

  double dydx = 0.;
  double pt = 0., px = 0., py = 0.;

  //remember (radius rho in cm):
  //rho =
  //100. * pt *
  //(10./(3.*MagneticField::inTesla(GlobalPoint(0., 0., 0.)).z()));

  double rho = theCircle.rho();
  pt = 0.01 * rho * (0.3 * mField->inTesla(GlobalPoint(0, 0, 0)).z());

  // (py/px)|x=v.x() = (dy/dx)|x=v.x()
  //remember:
  //y(x) = +-sqrt(rho^2 - (x-x0)^2) + y0
  //y(x) =  sqrt(rho^2 - (x-x0)^2) + y0  if y(x) >= y0
  //y(x) = -sqrt(rho^2 - (x-x0)^2) + y0  if y(x) < y0
  //=> (dy/dx) = -(x-x0)/sqrt(Q)  if y(x) >= y0
  //   (dy/dx) =  (x-x0)/sqrt(Q)  if y(x) < y0
  //with Q = rho^2 - (x-x0)^2

  double arg = rho * rho - ((v.x() - theCircle.x0()) * (v.x() - theCircle.x0()));

  if (arg > 0) {
    //  double root = sqrt(  rho*rho - ( (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) )  );
    double root = sqrt(arg);

    if ((v.y() - theCircle.y0()) > 0.)
      dydx = -(v.x() - theCircle.x0()) / root;
    else
      dydx = (v.x() - theCircle.x0()) / root;

    px = pt / sqrt(1. + dydx * dydx);
    py = px * dydx;
    // check sign with scalar product
    if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
      px *= -1.;
      py *= -1.;
    }

    //std::cout << " ConversionFastHelix:helixStateAtVertex  rho " << rho  << " pt " << pt  << " v " <<  v << " theCircle.x0() " <<theCircle.x0() << " theCircle.y0() "  <<  theCircle.y0() << " v.x()-theCircle.x0() "  << v.x()-theCircle.x0() << " rho^2 " << rho*rho << "  v.x()-theCircle.x0()^2 " <<   (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) <<  " root " << root << " arg " << arg <<  " dydx " << dydx << std::endl;
    //calculate z0, pz
    //(z, R*phi) linear relation in a helix
    //with R, phi defined as radius and angle w.r.t. centre of circle
    //in transverse plane
    //pz = pT*(dz/d(R*phi)))

    FastLine flfit(theOuterHit, theMiddleHit, theCircle.rho());

    double z_0 = 0;

    //std::cout << " ConversionFastHelix:helixStateAtVertex  flfit.n2() " <<  flfit.n2() << " flfit.c() " << flfit.c() << " flfit.n2() " << flfit.n2() << std::endl;
    if (flfit.n2() != 0 && !edm::isNotFinite(flfit.c()) && !edm::isNotFinite(flfit.n2())) {
      //  std::cout << " Accepted " << std::endl;
      z_0 = -flfit.c() / flfit.n2();
      double dzdrphi = -flfit.n1() / flfit.n2();
      double pz = pt * dzdrphi;

      //get sign of particle

      GlobalVector magvtx = mField->inTesla(v);
      TrackCharge q = ((theCircle.x0() * py - theCircle.y0() * px) / (magvtx.z()) < 0.) ? -1 : 1;

      AlgebraicSymMatrix55 C = AlgebraicMatrixID();
      //MP

      atVertex = FTS(GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, mField),
                     CurvilinearTrajectoryError(C));

      //std::cout << " ConversionFastHelix:helixStateAtVertex globalPoint " << GlobalPoint(v.x(), v.y(), z_0) << " GlobalVector " << GlobalVector(px, py, pz)  << " q " << q << " MField " << mField->inTesla(v) << std::endl;
      //std::cout << " ConversionFastHelix:helixStateAtVertex atVertex.transverseCurvature() " << atVertex.transverseCurvature() << std::endl;
      if (atVertex.transverseCurvature() != 0) {
        validStateAtVertex = true;

        //std::cout << " ConversionFastHelix:helixStateAtVertex validHelixStateAtVertex status " << validStateAtVertex << std::endl;
        return atVertex;
      } else
        return atVertex;
    } else {
      //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted  validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
      return atVertex;
    }

  } else {
    //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted because arg <0 validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
    return atVertex;
  }
}

FreeTrajectoryState ConversionFastHelix::straightLineStateAtVertex() {
  FTS atVertex;

  //calculate FTS assuming straight line...

  GlobalPoint pMid(theMiddleHit);
  GlobalPoint v(theVertex);

  double dydx = 0.;
  double pt = 0., px = 0., py = 0.;

  if (fabs(theCircle.n1()) > 0. || fabs(theCircle.n2()) > 0.)
    pt = 1.e+4;  // 10 TeV //else no pt
  if (fabs(theCircle.n2()) > 0.) {
    dydx = -theCircle.n1() / theCircle.n2();  //else px = 0
  }

  if (pt == 0 && dydx == 0.) {
    validStateAtVertex = false;
    return atVertex;
  }
  px = pt / sqrt(1. + dydx * dydx);
  py = px * dydx;

  // check sign with scalar product
  if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
    px *= -1.;
    py *= -1.;
  }

  //calculate z_0 and pz at vertex using weighted mean
  //z = z(r) = z0 + (dz/dr)*r
  //tan(theta) = dr/dz = (dz/dr)^-1
  //theta = atan(1./dzdr)
  //p = pt/sin(theta)
  //pz = p*cos(theta) = pt/tan(theta)

  FastLine flfit(theOuterHit, theMiddleHit);

  double z_0 = 0;
  if (flfit.n2() != 0 && !edm::isNotFinite(flfit.c()) && !edm::isNotFinite(flfit.n2())) {
    z_0 = -flfit.c() / flfit.n2();

    double dzdr = -flfit.n1() / flfit.n2();
    double pz = pt * dzdr;

    TrackCharge q = 1;

    atVertex = FTS(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, mField);

    //    std::cout << "  ConversionFastHelix::straightLineStateAtVertex curvature " << atVertex.transverseCurvature() << "   signedInverseMomentum " << atVertex.signedInverseMomentum() << std::endl;
    if (atVertex.transverseCurvature() == -0) {
      return atVertex;
    } else {
      validStateAtVertex = true;
      return atVertex;
    }

  } else {
    return atVertex;
  }
}