Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
#ifndef ConeDefinition_h
#define ConeDefinition_h

#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "DataFormats/GeometryVector/interface/GlobalVector.h"

#include "CommonTools/UtilAlgos/interface/DeltaR.h"

inline double getDistInPlaneSimple(const GlobalPoint caloPoint, const GlobalPoint rechitPoint) {
  // Simplified version of getDistInPlane
  // Assume track direction is origin -> point of hcal intersection

  const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());

  const GlobalVector caloIntersectUnitVector = caloIntersectVector.unit();

  const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());

  const GlobalVector rechitUnitVector = rechitVector.unit();
  double dotprod = caloIntersectUnitVector.dot(rechitUnitVector);
  double rechitdist = caloIntersectVector.mag() / dotprod;

  const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
  const GlobalPoint effectiveRechitPoint(
      effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());

  GlobalVector distance_vector = effectiveRechitPoint - caloPoint;

  if (dotprod > 0.) {
    return distance_vector.mag();
  } else {
    return 999999.;
  }
}

inline double getDistInPlaneTrackDir(const GlobalPoint caloPoint,
                                     const GlobalVector caloVector,
                                     const GlobalPoint rechitPoint) {
  // Simplified version of getDistInPlane : no cone "within" Hcal, but
  // don't assume track direction is origin -> point of hcal
  // intersection.
  const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(),
                                         caloPoint.z());  //p

  const GlobalVector caloUnitVector = caloVector.unit();
  const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
  const GlobalVector rechitUnitVector = rechitVector.unit();
  double dotprod_denominator = caloUnitVector.dot(rechitUnitVector);
  double dotprod_numerator = caloUnitVector.dot(caloIntersectVector);
  double rechitdist = dotprod_numerator / dotprod_denominator;
  //  double rechitdist=caloIntersectVector.dot(rechitUnitVector);
  const GlobalVector effectiveRechitVector = rechitdist * rechitUnitVector;
  const GlobalPoint effectiveRechitPoint(
      effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());
  GlobalVector distance_vector = effectiveRechitPoint - caloPoint;
  if (dotprod_denominator > 0. && dotprod_numerator > 0.) {
    return distance_vector.mag();
  } else {
    return 999999.;
  }
}

inline double getDistInPlane(const GlobalVector trackDirection,
                             const GlobalPoint caloPoint,
                             const GlobalPoint rechitPoint,
                             double coneHeight) {
  // The iso track candidate hits the Calo (Ecal or Hcal) at "caloPoint"
  // with direction "trackDirection".

  // "rechitPoint" is the position of the rechit.  We only care about
  // the direction of the rechit.

  // Consider the rechitPoint as characterized by angles theta and phi
  // wrt the origin which points at the calo cell of the rechit.  In
  // some sense the distance along the line theta/phi is arbitrary. A
  // simplified choice might be to put the rechit at the surface of
  // the Hcal.  Our approach is to see whether/where this line
  // intersects a cone of height "coneHeight" with vertex at caloPoint
  // aligned with trackDirection.
  // To that end, this function returns the distance between the
  // center of the base of the cone and the intersection of the rechit
  // line and the plane that contains the base of the cone.  This
  // distance is compared with the radius of the cone outside this
  // function.

  // Make vector of length cone height along track direction
  const GlobalVector heightVector = trackDirection * coneHeight;

  // Make vector from origin to point where iso track intersects the
  // calorimeter.
  const GlobalVector caloIntersectVector(caloPoint.x(), caloPoint.y(), caloPoint.z());

  // Make vector from origin to point in center of base of cone
  const GlobalVector coneBaseVector = heightVector + caloIntersectVector;

  // Make point in center of base of cone
  const GlobalPoint coneBasePoint(coneBaseVector.x(), coneBaseVector.y(), coneBaseVector.z());

  // Make unit vector pointing at rechit.
  const GlobalVector rechitVector(rechitPoint.x(), rechitPoint.y(), rechitPoint.z());
  const GlobalVector rechitDirection = rechitVector.unit();

  // Find distance "r" along "rechit line" (with angles theta2 and
  // phi2) where line intersects plane defined by base of cone.

  // Definition plane of that contains base of cone:
  // trackDirection.x() (x - coneBaseVector.x()) +
  // trackDirection.y() (y - coneBaseVector.y()) +
  // trackDirection.z() (z - coneBaseVector.z()) = 0

  // Point P_{rh} where rechit line intersects plane:
  // (rechitdist sin(theta2) cos(phi2),
  //  rechitdist sin(theta2) cos(phi2),
  //  rechitdist cos(theta2))

  // Substitute P_{rh} into equation for plane and solve for rechitdist.
  // rechitDist turns out to be the ratio of dot products:

  double rechitdist = trackDirection.dot(coneBaseVector) / trackDirection.dot(rechitDirection);

  // Now find distance between point at center of cone base and point
  // where the "rechit line" intersects the plane defined by the base
  // of the cone; i.e. the effectiveRecHitPoint.
  const GlobalVector effectiveRechitVector = rechitdist * rechitDirection;
  const GlobalPoint effectiveRechitPoint(
      effectiveRechitVector.x(), effectiveRechitVector.y(), effectiveRechitVector.z());

  GlobalVector distance_vector = effectiveRechitPoint - coneBasePoint;
  return distance_vector.mag();
}

#endif