CachedTrajectory Class Reference

#include <CachedTrajectory.h>

Public Member Functions
const std::vector< SteppingHelixStateInfo > &	getEcalTrajectory () const
const std::vector< SteppingHelixStateInfo > &	getHcalTrajectory () const
const std::vector< SteppingHelixStateInfo > &	getHOTrajectory () const
const std::vector< SteppingHelixStateInfo > &	getPreshowerTrajectory () const

Protected Member Functions
std::pair< float, float >	delta (const double &theta1, const double &theta2, const double &phi1, const double &phi2)
float	distance (const Plane *plane, int index)

Static Protected Member Functions
static int	sign (float number)

Protected Attributes
std::vector< SteppingHelixStateInfo >	ecalTrajectory_
std::deque< SteppingHelixStateInfo >	fullTrajectory_
bool	fullTrajectoryFilled_
std::vector< SteppingHelixStateInfo >	hcalTrajectory_
float	HOmaxRho_
float	HOmaxZ_
std::vector< SteppingHelixStateInfo >	hoTrajectory_
float	maxRho_
float	maxZ_
float	minRho_
float	minZ_
std::vector< SteppingHelixStateInfo >	preshowerTrajectory_
const Propagator *	propagator_
SteppingHelixStateInfo	stateAtIP_
float	step_
std::vector< GlobalPoint >	wideEcalTrajectory_
std::vector< GlobalPoint >	wideHcalTrajectory_
std::vector< GlobalPoint >	wideHOTrajectory_

Private Types
enum	TrajectorType { IpToEcal, IpToHcal, IpToHO, FullTrajectory }
enum	WideTrajectoryType { Ecal, Hcal, HO }

Private Member Functions
	CachedTrajectory ()
void	findEcalTrajectory (const FiducialVolume &)
void	findHcalTrajectory (const FiducialVolume &)
void	findHOTrajectory (const FiducialVolume &)
void	findPreshowerTrajectory (const FiducialVolume &)
SteppingHelixStateInfo	getInnerState ()
SteppingHelixStateInfo	getOuterState ()
float	getPropagationStep () const
SteppingHelixStateInfo	getStateAtEcal ()
SteppingHelixStateInfo	getStateAtHcal ()
SteppingHelixStateInfo	getStateAtHO ()
SteppingHelixStateInfo	getStateAtPreshower ()
void	getTrajectory (std::vector< SteppingHelixStateInfo > &, const FiducialVolume &, int steps=4)
std::vector< GlobalPoint > *	getWideTrajectory (const std::vector< SteppingHelixStateInfo > &, WideTrajectoryType)
void	propagate (SteppingHelixStateInfo &state, const Plane &plane)
void	propagate (SteppingHelixStateInfo &state, const Cylinder &cylinder)
TrajectoryStateOnSurface	propagate (const Plane *plane)
	get fast to a given DetId surface using cached trajectory More...
bool	propagateAll (const SteppingHelixStateInfo &initialState)
	propagate through the whole detector, returns true if successful More...
void	propagateForward (SteppingHelixStateInfo &state, float distance)
void	reset_trajectory ()
void	setMaxDetectorLength (float l=2200.)
void	setMaxDetectorRadius (float r=800.)
void	setMaxHOLength (float l=2200.)
void	setMaxHORadius (float r=800.)
void	setMinDetectorLength (float l=0.)
void	setMinDetectorRadius (float r=0.)
void	setPropagationStep (float s=20.)
void	setPropagator (const Propagator *ptr)
void	setStateAtIP (const SteppingHelixStateInfo &state)
std::pair< float, float >	trajectoryDelta (TrajectorType)

Friends
std::vector< SteppingHelixStateInfo >	propagateThoughFromIP (const SteppingHelixStateInfo &state, const Propagator *ptr, const FiducialVolume &volume, int nsteps, float step, float minR, float minZ, float maxR, float maxZ)
class	TrackDetectorAssociator

Detailed Description

Definition at line 44 of file CachedTrajectory.h.

Member Enumeration Documentation

TrajectorType

enum CachedTrajectory::TrajectorType

private

Enumerator
IpToEcal
IpToHcal
IpToHO
FullTrajectory

Definition at line 64 of file CachedTrajectory.h.

{ IpToEcal, IpToHcal, IpToHO, FullTrajectory };

WideTrajectoryType

enum CachedTrajectory::WideTrajectoryType

private

Enumerator
Ecal
Hcal
HO

Definition at line 65 of file CachedTrajectory.h.

{ Ecal, Hcal, HO };

Constructor & Destructor Documentation

CachedTrajectory()

CachedTrajectory::CachedTrajectory ( )

private

Definition at line 49 of file CachedTrajectory.cc.

References reset_trajectory(), setMaxDetectorLength(), setMaxDetectorRadius(), setMinDetectorLength(), setMinDetectorRadius(), and setPropagationStep().

                                   : propagator_(nullptr) {
  reset_trajectory();
  setMaxDetectorRadius();
  setMaxDetectorLength();
  setMinDetectorRadius();
  setMinDetectorLength();
  setPropagationStep();
}

Member Function Documentation

delta()

std::pair< float, float > CachedTrajectory::delta ( const double &  theta1, const double &  theta2, const double &  phi1, const double &  phi2  )

protected

Definition at line 319 of file CachedTrajectory.cc.

References M_PI, and mps_fire::result.

Referenced by trajectoryDelta().

                                                                    {
  std::pair<float, float> result(theta2 - theta1, phi2 - phi1);
  // this won't work for loopers, since deltaPhi cannot be larger than Pi.
  if (fabs(result.second) > 2 * M_PI - fabs(result.second)) {
    if (result.second > 0)
      result.second -= 2 * M_PI;
    else
      result.second += 2 * M_PI;
  }
  return result;
}

distance()

float CachedTrajectory::distance ( const Plane *  plane, int  index  )

inlineprotected

Definition at line 136 of file CachedTrajectory.h.

References fullTrajectory_, createfilelist::int, and position.

Referenced by propagate(), and propagateForward().

                                                             {
    if (index < 0 || fullTrajectory_.empty() || (unsigned int)index >= fullTrajectory_.size())
      return 0;
    return plane->localZ(fullTrajectory_[index].position());
  }

findEcalTrajectory()

void CachedTrajectory::findEcalTrajectory ( const FiducialVolume &  volume )

private

Definition at line 443 of file CachedTrajectory.cc.

References ecalTrajectory_, getTrajectory(), and LogTrace.

                                                                      {
  LogTrace("TrackAssociator") << "getting trajectory in ECAL";
  getTrajectory(ecalTrajectory_, volume, 4);
  LogTrace("TrackAssociator") << "# of points in ECAL trajectory:" << ecalTrajectory_.size();
}

findHcalTrajectory()

void CachedTrajectory::findHcalTrajectory ( const FiducialVolume &  volume )

private

Definition at line 461 of file CachedTrajectory.cc.

References getTrajectory(), hcalTrajectory_, and LogTrace.

                                                                      {
  LogTrace("TrackAssociator") << "getting trajectory in HCAL";
  getTrajectory(hcalTrajectory_, volume, 4);  // more steps to account for different depth
  LogTrace("TrackAssociator") << "# of points in HCAL trajectory:" << hcalTrajectory_.size();
}

findHOTrajectory()

void CachedTrajectory::findHOTrajectory ( const FiducialVolume &  volume )

private

Definition at line 469 of file CachedTrajectory.cc.

References getTrajectory(), hoTrajectory_, and LogTrace.

                                                                    {
  LogTrace("TrackAssociator") << "getting trajectory in HO";
  getTrajectory(hoTrajectory_, volume, 2);
  LogTrace("TrackAssociator") << "# of points in HO trajectory:" << hoTrajectory_.size();
}

findPreshowerTrajectory()

void CachedTrajectory::findPreshowerTrajectory ( const FiducialVolume &  volume )

private

Definition at line 449 of file CachedTrajectory.cc.

References getTrajectory(), LogTrace, and preshowerTrajectory_.

                                                                           {
  LogTrace("TrackAssociator") << "getting trajectory in Preshower";
  getTrajectory(preshowerTrajectory_, volume, 2);
  LogTrace("TrackAssociator") << "# of points in Preshower trajectory:" << preshowerTrajectory_.size();
}

getEcalTrajectory()

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getEcalTrajectory

(

)

const

Definition at line 455 of file CachedTrajectory.cc.

References ecalTrajectory_.

{ return ecalTrajectory_; }

getHcalTrajectory()

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getHcalTrajectory

(

)

const

Definition at line 467 of file CachedTrajectory.cc.

References hcalTrajectory_.

{ return hcalTrajectory_; }

getHOTrajectory()

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getHOTrajectory

(

)

const

Definition at line 475 of file CachedTrajectory.cc.

References hoTrajectory_.

{ return hoTrajectory_; }

getInnerState()

SteppingHelixStateInfo CachedTrajectory::getInnerState ( )

private

Definition at line 598 of file CachedTrajectory.cc.

References fullTrajectory_.

                                                       {
  if (fullTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return fullTrajectory_.front();
}

getOuterState()

SteppingHelixStateInfo CachedTrajectory::getOuterState ( )

private

Definition at line 605 of file CachedTrajectory.cc.

References fullTrajectory_.

                                                       {
  if (fullTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return fullTrajectory_.back();
}

getPreshowerTrajectory()

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getPreshowerTrajectory

(

)

const

Definition at line 457 of file CachedTrajectory.cc.

References preshowerTrajectory_.

                                                                                        {
  return preshowerTrajectory_;
}

getPropagationStep()

float CachedTrajectory::getPropagationStep ( ) const

inlineprivate

Definition at line 121 of file CachedTrajectory.h.

References step_.

{ return step_; }

getStateAtEcal()

SteppingHelixStateInfo CachedTrajectory::getStateAtEcal ( )

private

Definition at line 570 of file CachedTrajectory.cc.

References ecalTrajectory_.

                                                        {
  if (ecalTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return ecalTrajectory_.front();
}

getStateAtHcal()

SteppingHelixStateInfo CachedTrajectory::getStateAtHcal ( )

private

Definition at line 584 of file CachedTrajectory.cc.

References hcalTrajectory_.

                                                        {
  if (hcalTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return hcalTrajectory_.front();
}

getStateAtHO()

SteppingHelixStateInfo CachedTrajectory::getStateAtHO ( )

private

Definition at line 591 of file CachedTrajectory.cc.

References hoTrajectory_.

                                                      {
  if (hoTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return hoTrajectory_.front();
}

getStateAtPreshower()

SteppingHelixStateInfo CachedTrajectory::getStateAtPreshower ( )

private

Definition at line 577 of file CachedTrajectory.cc.

References preshowerTrajectory_.

                                                             {
  if (preshowerTrajectory_.empty())
    return SteppingHelixStateInfo();
  else
    return preshowerTrajectory_.front();
}

getTrajectory()

void CachedTrajectory::getTrajectory ( std::vector< SteppingHelixStateInfo > &  trajectory, const FiducialVolume &  volume, int  steps = 4  )

private

get a set of points representing the trajectory between two cylinders of radius R1 and R2 and length L1 and L2. Parameter steps defines maximal number of steps in the detector.

Definition at line 334 of file CachedTrajectory.cc.

References Reference_intrackfit_cff::barrel, Plane::build(), Cylinder::build(), HGC3DClusterGenMatchSelector_cfi::dR, l1ctLayer1_cff::dZ, makeMuonMisalignmentScenario::endcap, Exception, fullTrajectory_, fullTrajectoryFilled_, mps_fire::i, FiducialVolume::isValid(), SteppingHelixStateInfo::isValid(), LogTrace, SiStripPI::max, FiducialVolume::maxR(), FiducialVolume::maxZ(), FiducialVolume::minR(), FiducialVolume::minZ(), SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::perp(), SteppingHelixStateInfo::position(), position, propagate(), propagateForward(), customisers::steps, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), z, and PV3DBase< T, PVType, FrameType >::z().

Referenced by findEcalTrajectory(), findHcalTrajectory(), findHOTrajectory(), findPreshowerTrajectory(), and propagateThoughFromIP().

                                                {
  if (!fullTrajectoryFilled_)
    throw cms::Exception("FatalError") << "trajectory is not defined yet. Please use propagateAll first.";
  if (fullTrajectory_.empty()) {
    LogTrace("TrackAssociator") << "Trajectory is empty. Move on";
    return;
  }
  if (!volume.isValid()) {
    LogTrace("TrackAssociator") << "no trajectory is expected to be found since the fiducial volume is not valid";
    return;
  }
  double step = std::max(volume.maxR() - volume.minR(), volume.maxZ() - volume.minZ()) / steps;

  int closestPointOnLeft = -1;

  // check whether the trajectory crossed the region
  if (!((fullTrajectory_.front().position().perp() < volume.maxR() &&
         fabs(fullTrajectory_.front().position().z()) < volume.maxZ()) &&
        (fullTrajectory_.back().position().perp() > volume.minR() ||
         fabs(fullTrajectory_.back().position().z()) > volume.minZ()))) {
    LogTrace("TrackAssociator") << "Track didn't cross the region (R1,R2,L1,L2): " << volume.minR() << ", "
                                << volume.maxR() << ", " << volume.minZ() << ", " << volume.maxZ();
    return;
  }

  // get distance along momentum to the surface.

  // the following code can be made faster, but it'll hardly be a significant improvement
  // simplifications:
  //   1) direct loop over stored trajectory points instead of some sort
  //      of fast root search (Newton method)
  //   2) propagate from the closest point outside the region with the
  //      requested step ignoring stored trajectory points.
  double dZ(-1.);
  double dR(-1.);
  int firstPointInside(-1);
  for (unsigned int i = 0; i < fullTrajectory_.size(); i++) {
    // LogTrace("TrackAssociator") << "Trajectory info (i,perp,r1,r2,z,z1,z2): " << i << ", " << fullTrajectory_[i].position().perp() <<
    //  ", " << volume.minR() << ", " << volume.maxR() << ", " << fullTrajectory_[i].position().z() << ", " << volume.minZ() << ", " <<
    //  volume.maxZ() << ", " << closestPointOnLeft;
    dR = fullTrajectory_[i].position().perp() - volume.minR();
    dZ = fabs(fullTrajectory_[i].position().z()) - volume.minZ();
    if (dR > 0 || dZ > 0) {
      if (i > 0) {
        firstPointInside = i;
        closestPointOnLeft = i - 1;
      } else {
        firstPointInside = 0;
        closestPointOnLeft = 0;
      }
      break;
    }
  }
  if (closestPointOnLeft == -1)
    throw cms::Exception("FatalError") << "This shouls never happen - internal logic error";

  SteppingHelixStateInfo currentState(fullTrajectory_[closestPointOnLeft]);
  if (currentState.position().x() * currentState.momentum().x() +
          currentState.position().y() * currentState.momentum().y() +
          currentState.position().z() * currentState.momentum().z() <
      0)
    step = -step;

  // propagate to the inner surface of the active volume

  if (firstPointInside != closestPointOnLeft) {
    if (dR > 0) {
      Cylinder::CylinderPointer barrel =
          Cylinder::build(volume.minR(), Cylinder::PositionType(0, 0, 0), Cylinder::RotationType());
      propagate(currentState, *barrel);
    } else {
      Plane::PlanePointer endcap =
          Plane::build(Plane::PositionType(0, 0, currentState.position().z() > 0 ? volume.minZ() : -volume.minZ()),
                       Plane::RotationType());
      propagate(currentState, *endcap);
    }
    if (currentState.isValid())
      trajectory.push_back(currentState);
  } else
    LogTrace("TrackAssociator") << "Weird message\n";

  while (currentState.isValid() && currentState.position().perp() < volume.maxR() &&
         fabs(currentState.position().z()) < volume.maxZ()) {
    propagateForward(currentState, step);
    if (!currentState.isValid()) {
      LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
      break;
    }
    // LogTrace("TrackAssociator") << "New state (perp, z): " << currentState.position().perp() << ", " << currentState.position().z();
    //if ( ( currentState.position().perp() < volume.maxR() && fabs(currentState.position().z()) < volume.maxZ() ) &&
    //     ( currentState.position().perp()-volume.minR() > 0  || fabs(currentState.position().z()) - volume.minZ() >0 ) )
    trajectory.push_back(currentState);
  }
}

getWideTrajectory()

std::vector< GlobalPoint > * CachedTrajectory::getWideTrajectory ( const std::vector< SteppingHelixStateInfo > &  states, WideTrajectoryType  wideTrajectoryType  )

private

Definition at line 477 of file CachedTrajectory.cc.

References funct::cos(), makePileupJSON::denom, Ecal, Hcal, HO, LogTrace, M_PI, M_PI_2, phi, funct::pow(), dttmaxenums::r32, LocalError::rotate(), idealTransformation::rotation, LocalError::scale(), funct::sin(), mathSSE::sqrt(), filterCSVwithJSON::target, trackerHitRTTI::vector, wideEcalTrajectory_, wideHcalTrajectory_, wideHOTrajectory_, geometryCSVtoXML::xx, LocalError::xx(), geometryCSVtoXML::xy, LocalError::xy(), geometryCSVtoXML::yy, and LocalError::yy().

                                                                                                     {
  std::vector<GlobalPoint>* wideTrajectory = nullptr;
  switch (wideTrajectoryType) {
    case Ecal:
      LogTrace("TrackAssociator") << "Filling ellipses in Ecal trajectory";
      wideTrajectory = &wideEcalTrajectory_;
      break;
    case Hcal:
      LogTrace("TrackAssociator") << "Filling ellipses in Hcal trajectory";
      wideTrajectory = &wideHcalTrajectory_;
      break;
    case HO:
      LogTrace("TrackAssociator") << "Filling ellipses in HO trajectory";
      wideTrajectory = &wideHOTrajectory_;
      break;
  }
  if (!wideTrajectory)
    return nullptr;

  for (std::vector<SteppingHelixStateInfo>::const_iterator state = states.begin(); state != states.end(); state++) {
    // defined a normal plane wrt the particle trajectory direction
    // let's hope that I computed the rotation matrix correctly.
    GlobalVector vector(state->momentum().unit());
    float r21 = 0;
    float r22 = vector.z() / sqrt(1 - pow(vector.x(), 2));
    float r23 = -vector.y() / sqrt(1 - pow(vector.x(), 2));
    float r31 = vector.x();
    float r32 = vector.y();
    float r33 = vector.z();
    float r11 = r22 * r33 - r23 * r32;
    float r12 = r23 * r31;
    float r13 = -r22 * r31;

    Plane::RotationType rotation(r11, r12, r13, r21, r22, r23, r31, r32, r33);
    Plane* target = new Plane(state->position(), rotation);

    TrajectoryStateOnSurface tsos = state->getStateOnSurface(*target);

    if (!tsos.isValid()) {
      LogTrace("TrackAssociator") << "[getWideTrajectory] TSOS not valid";
      continue;
    }
    if (!tsos.hasError()) {
      LogTrace("TrackAssociator") << "[getWideTrajectory] TSOS does not have Errors";
      continue;
    }
    LocalError localErr = tsos.localError().positionError();
    localErr.scale(2);  // get the 2 sigma ellipse
    float xx = localErr.xx();
    float xy = localErr.xy();
    float yy = localErr.yy();

    float denom = yy - xx;
    float phi = 0.;
    if (xy == 0 && denom == 0)
      phi = M_PI_4;
    else
      phi = 0.5 * atan2(2. * xy, denom);  // angle of MAJOR axis
    // Unrotate the error ellipse to get the semimajor and minor axes. Then place points on
    // the endpoints of semiminor an seminajor axes on original(rotated) error ellipse.
    LocalError rotErr = localErr.rotate(-phi);  // xy covariance of rotErr should be zero
    float semi1 = sqrt(rotErr.xx());
    float semi2 = sqrt(rotErr.yy());

    // Just use one point if the ellipse is small
    // if(semi1 < 0.1 && semi2 < 0.1) {
    //   LogTrace("TrackAssociator") << "[getWideTrajectory] Error ellipse is small, using one trajectory point";
    //   wideTrajectory->push_back(state->position());
    //   continue;
    // }

    Local2DPoint bounds[4];
    bounds[0] = Local2DPoint(semi1 * cos(phi), semi1 * sin(phi));
    bounds[1] = Local2DPoint(semi1 * cos(phi + M_PI), semi1 * sin(phi + M_PI));
    phi += M_PI_2;  // add pi/2 for the semi2 axis
    bounds[2] = Local2DPoint(semi2 * cos(phi), semi2 * sin(phi));
    bounds[3] = Local2DPoint(semi2 * cos(phi + M_PI), semi2 * sin(phi + M_PI));

    // LogTrace("TrackAssociator") << "Axes " << semi1 <<","<< semi2 <<"   phi "<< phi;
    // LogTrace("TrackAssociator") << "Local error ellipse: " << bounds[0] << bounds[1] << bounds[2] << bounds[3];

    wideTrajectory->push_back(state->position());
    for (int index = 0; index < 4; ++index)
      wideTrajectory->push_back(target->toGlobal(bounds[index]));

    //LogTrace("TrackAssociator") <<"Global error ellipse: (" << target->toGlobal(bounds[0]) <<","<< target->toGlobal(bounds[1])
    //         <<","<< target->toGlobal(bounds[2]) <<","<< target->toGlobal(bounds[3]) <<","<<state->position() <<")";
  }

  return wideTrajectory;
}

propagate() [1/3]

void CachedTrajectory::propagate ( SteppingHelixStateInfo &  state, const Plane &  plane  )

private

Definition at line 77 of file CachedTrajectory.cc.

References cms::Exception::category(), cms::Exception::explainSelf(), TrajectoryStateOnSurface::freeState(), Propagator::propagate(), and propagator_.

Referenced by getTrajectory(), and propagateForward().

                                                                                  {
  if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
    try {
      shp->propagate(state, plane, state);
    } catch (cms::Exception& ex) {
      edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
      edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
                                         << state.momentum().x() << ", " << state.momentum().y() << ", "
                                         << state.momentum().z();
      state = SteppingHelixStateInfo();
    }
  } else {
    FreeTrajectoryState fts;
    state.getFreeState(fts);
    TrajectoryStateOnSurface stateOnSurface = propagator_->propagate(fts, plane);
    state = SteppingHelixStateInfo(*(stateOnSurface.freeState()));
  }
}

propagate() [2/3]

void CachedTrajectory::propagate ( SteppingHelixStateInfo &  state, const Cylinder &  cylinder  )

private

Definition at line 96 of file CachedTrajectory.cc.

References cms::Exception::category(), cms::Exception::explainSelf(), TrajectoryStateOnSurface::freeState(), Propagator::propagate(), and propagator_.

                                                                                        {
  if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
    try {
      shp->propagate(state, cylinder, state);
    } catch (cms::Exception& ex) {
      edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
      edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
                                         << state.momentum().x() << ", " << state.momentum().y() << ", "
                                         << state.momentum().z();
      state = SteppingHelixStateInfo();
    }
  } else {
    FreeTrajectoryState fts;
    state.getFreeState(fts);
    TrajectoryStateOnSurface stateOnSurface = propagator_->propagate(fts, cylinder);
    state = SteppingHelixStateInfo(*(stateOnSurface.freeState()));
  }
}

propagate() [3/3]

TrajectoryStateOnSurface CachedTrajectory::propagate ( const Plane *  plane )

private

get fast to a given DetId surface using cached trajectory

Definition at line 179 of file CachedTrajectory.cc.

References cms::Exception::category(), distance(), PV3DBase< T, PVType, FrameType >::eta(), cms::Exception::explainSelf(), fullTrajectory_, createfilelist::int, LogTrace, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), GloballyPositioned< T >::position(), Propagator::propagate(), propagator_, sign(), and PV3DBase< T, PVType, FrameType >::z().

                                                                       {
  // TimerStack timers(TimerStack::Disableable);
  // timers.benchmark("CachedTrajectory::propagate::benchmark");
  // timers.push("CachedTrajectory::propagate",TimerStack::FastMonitoring);
  // timers.push("CachedTrajectory::propagate::findClosestPoint",TimerStack::FastMonitoring);

  // Assume that all points along the trajectory are equally spread out.
  // For simplication assume that the trajectory is just a straight
  // line and find a point closest to the target plane. Propagate to
  // the plane from the point.

  const float matchingDistance = 1;
  // find the closest point to the plane
  int leftIndex = 0;
  int rightIndex = fullTrajectory_.size() - 1;
  int closestPointOnLeft = 0;

  // check whether the trajectory crossed the plane (signs should be different)
  if (sign(distance(plane, leftIndex)) * sign(distance(plane, rightIndex)) != -1) {
    LogTrace("TrackAssociator") << "Track didn't cross the plane:\n\tleft distance: " << distance(plane, leftIndex)
                                << "\n\tright distance: " << distance(plane, rightIndex);
    return TrajectoryStateOnSurface();
  }

  while (leftIndex + 1 < rightIndex) {
    closestPointOnLeft = int((leftIndex + rightIndex) / 2);
    float dist = distance(plane, closestPointOnLeft);
    // LogTrace("TrackAssociator") << "Closest point on left: " << closestPointOnLeft << "\n"
    //    << "Distance to the plane: " << dist;
    if (fabs(dist) < matchingDistance) {
      // found close match, verify that we are on the left side
      if (closestPointOnLeft > 0 && sign(distance(plane, closestPointOnLeft - 1)) * dist == -1)
        closestPointOnLeft--;
      break;
    }

    // check where is the plane
    if (sign(distance(plane, leftIndex) * dist) == -1)
      rightIndex = closestPointOnLeft;
    else
      leftIndex = closestPointOnLeft;

    // LogTrace("TrackAssociator") << "Distance on left: " << distance(plane, leftIndex) << "\n"
    //  << "Distance to closest point: " <<  distance(plane, closestPointOnLeft) << "\n"
    //  << "Left index: " << leftIndex << "\n"
    //  << "Right index: " << rightIndex;
  }
  LogTrace("TrackAssociator") << "closestPointOnLeft: " << closestPointOnLeft << "\n\ttrajectory point (z,R,eta,phi): "
                              << fullTrajectory_[closestPointOnLeft].position().z() << ", "
                              << fullTrajectory_[closestPointOnLeft].position().perp() << " , "
                              << fullTrajectory_[closestPointOnLeft].position().eta() << " , "
                              << fullTrajectory_[closestPointOnLeft].position().phi()
                              << "\n\tplane center (z,R,eta,phi): " << plane->position().z() << ", "
                              << plane->position().perp() << " , " << plane->position().eta() << " , "
                              << plane->position().phi();

  // propagate to the plane
  // timers.pop_and_push("CachedTrajectory::propagate::localPropagation",TimerStack::FastMonitoring);
  if (const SteppingHelixPropagator* shp = dynamic_cast<const SteppingHelixPropagator*>(propagator_)) {
    SteppingHelixStateInfo state;
    try {
      shp->propagate(fullTrajectory_[closestPointOnLeft], *plane, state);
    } catch (cms::Exception& ex) {
      edm::LogWarning("TrackAssociator") << "Caught exception " << ex.category() << ": " << ex.explainSelf();
      edm::LogWarning("TrackAssociator") << "An exception is caught during the track propagation\n"
                                         << state.momentum().x() << ", " << state.momentum().y() << ", "
                                         << state.momentum().z();
      return TrajectoryStateOnSurface();
    }
    return state.getStateOnSurface(*plane);
  } else {
    FreeTrajectoryState fts;
    fullTrajectory_[closestPointOnLeft].getFreeState(fts);
    return propagator_->propagate(fts, *plane);
  }
}

propagateAll()

bool CachedTrajectory::propagateAll ( const SteppingHelixStateInfo &  initialState )

private

propagate through the whole detector, returns true if successful

Definition at line 115 of file CachedTrajectory.cc.

References Exception, fullTrajectory_, fullTrajectoryFilled_, SteppingHelixStateInfo::isValid(), LogTrace, PV3DBase< T, PVType, FrameType >::mag(), maxRho_, maxZ_, minRho_, minZ_, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), SteppingHelixStateInfo::position(), propagateForward(), propagator_, reset_trajectory(), step_, and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagateThoughFromIP().

                                                                              {
  if (fullTrajectoryFilled_) {
    edm::LogWarning("TrackAssociator")
        << "Reseting all trajectories. Please call reset_trajectory() explicitely to avoid this message";
    reset_trajectory();
  }

  //   TimerStack timers(TimerStack::Disableable);

  reset_trajectory();
  if (propagator_ == nullptr)
    throw cms::Exception("FatalError") << "Track propagator is not defined\n";
  SteppingHelixStateInfo currentState(initialState);
  fullTrajectory_.push_back(currentState);

  while (currentState.position().perp() < maxRho_ && fabs(currentState.position().z()) < maxZ_) {
    LogTrace("TrackAssociator") << "[propagateAll] Propagate outward from (rho, r, z, phi) ("
                                << currentState.position().perp() << ", " << currentState.position().mag() << ", "
                                << currentState.position().z() << ", " << currentState.position().phi() << ")";
    propagateForward(currentState, step_);
    if (!currentState.isValid()) {
      LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
      break;
    }
    LogTrace("TrackAssociator") << "\treached (rho, r, z, phi) (" << currentState.position().perp() << ", "
                                << currentState.position().mag() << ", " << currentState.position().z() << ", "
                                << currentState.position().phi() << ")";
    fullTrajectory_.push_back(currentState);
  }

  SteppingHelixStateInfo currentState2(initialState);
  SteppingHelixStateInfo previousState;
  while (currentState2.position().perp() > minRho_ || fabs(currentState2.position().z()) > minZ_) {
    previousState = currentState2;
    propagateForward(currentState2, -step_);
    if (!currentState2.isValid()) {
      LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n";
      break;
    }
    if (previousState.position().perp() - currentState2.position().perp() < 0) {
      LogTrace("TrackAssociator")
          << "Error: TrackAssociator has propogated the particle past the point of closest approach to IP" << std::endl;
      break;
    }
    LogTrace("TrackAssociator") << "[propagateAll] Propagated inward from (rho, r, z, phi) ("
                                << previousState.position().perp() << ", " << previousState.position().mag() << ", "
                                << previousState.position().z() << "," << previousState.position().phi() << ") to ("
                                << currentState2.position().perp() << ", " << currentState2.position().mag() << ", "
                                << currentState2.position().z() << ", " << currentState2.position().phi() << ")";
    fullTrajectory_.push_front(currentState2);
  }

  // LogTrace("TrackAssociator") << "fullTrajectory_ has " << fullTrajectory_.size() << " states with (R, z):\n";
  // for(unsigned int i=0; i<fullTrajectory_.size(); i++) {
  //  LogTrace("TrackAssociator") << "state " << i << ": (" << fullTrajectory_[i].position().perp() << ", "
  //    << fullTrajectory_[i].position().z() << ")\n";
  // }

  LogTrace("TrackAssociator") << "Done with the track propagation in the detector. Number of steps: "
                              << fullTrajectory_.size();
  fullTrajectoryFilled_ = true;
  return !fullTrajectory_.empty();
}

propagateForward()

void CachedTrajectory::propagateForward ( SteppingHelixStateInfo &  state, float  distance  )

private

Definition at line 58 of file CachedTrajectory.cc.

References Plane::build(), distance(), funct::pow(), propagate(), dttmaxenums::r32, idealTransformation::rotation, mathSSE::sqrt(), filterCSVwithJSON::target, and trackerHitRTTI::vector.

Referenced by getTrajectory(), and propagateAll().

                                                                                     {
  // defined a normal plane wrt the particle trajectory direction
  // let's hope that I computed the rotation matrix correctly.
  GlobalVector vector(state.momentum().unit());
  float r21 = 0;
  float r22 = vector.z() / sqrt(1 - pow(vector.x(), 2));
  float r23 = -vector.y() / sqrt(1 - pow(vector.x(), 2));
  float r31 = vector.x();
  float r32 = vector.y();
  float r33 = vector.z();
  float r11 = r22 * r33 - r23 * r32;
  float r12 = r23 * r31;
  float r13 = -r22 * r31;

  Surface::RotationType rotation(r11, r12, r13, r21, r22, r23, r31, r32, r33);
  Plane::PlanePointer target = Plane::build(state.position() + vector * distance, rotation);
  propagate(state, *target);
}

reset_trajectory()

void CachedTrajectory::reset_trajectory ( )

private

Definition at line 431 of file CachedTrajectory.cc.

References ecalTrajectory_, fullTrajectory_, fullTrajectoryFilled_, hcalTrajectory_, hoTrajectory_, preshowerTrajectory_, wideEcalTrajectory_, wideHcalTrajectory_, and wideHOTrajectory_.

Referenced by CachedTrajectory(), propagateAll(), and propagateThoughFromIP().

                                        {
  fullTrajectory_.clear();
  ecalTrajectory_.clear();
  hcalTrajectory_.clear();
  hoTrajectory_.clear();
  preshowerTrajectory_.clear();
  wideEcalTrajectory_.clear();
  wideHcalTrajectory_.clear();
  wideHOTrajectory_.clear();
  fullTrajectoryFilled_ = false;
}

setMaxDetectorLength()

void CachedTrajectory::setMaxDetectorLength ( float  l = 2200. )

inlineprivate

Definition at line 114 of file CachedTrajectory.h.

References MainPageGenerator::l, and maxZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ maxZ_ = l / 2.; }

setMaxDetectorRadius()

void CachedTrajectory::setMaxDetectorRadius ( float  r = 800. )

inlineprivate

Definition at line 113 of file CachedTrajectory.h.

References maxRho_, and alignCSCRings::r.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ maxRho_ = r; }

setMaxHOLength()

void CachedTrajectory::setMaxHOLength ( float  l = 2200. )

inlineprivate

Definition at line 116 of file CachedTrajectory.h.

References HOmaxZ_, and MainPageGenerator::l.

{ HOmaxZ_ = l / 2.; }

setMaxHORadius()

void CachedTrajectory::setMaxHORadius ( float  r = 800. )

inlineprivate

Definition at line 115 of file CachedTrajectory.h.

References HOmaxRho_, and alignCSCRings::r.

{ HOmaxRho_ = r; }

setMinDetectorLength()

void CachedTrajectory::setMinDetectorLength ( float  l = 0. )

inlineprivate

Definition at line 118 of file CachedTrajectory.h.

References MainPageGenerator::l, and minZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ minZ_ = l / 2.; }

setMinDetectorRadius()

void CachedTrajectory::setMinDetectorRadius ( float  r = 0. )

inlineprivate

Definition at line 117 of file CachedTrajectory.h.

References minRho_, and alignCSCRings::r.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ minRho_ = r; }

setPropagationStep()

void CachedTrajectory::setPropagationStep ( float  s = 20. )

inlineprivate

Definition at line 120 of file CachedTrajectory.h.

References alignCSCRings::s, and step_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ step_ = s; }

setPropagator()

void CachedTrajectory::setPropagator ( const Propagator *  ptr )

inlineprivate

Definition at line 83 of file CachedTrajectory.h.

References propagator_.

Referenced by propagateThoughFromIP().

{ propagator_ = ptr; }

setStateAtIP()

void CachedTrajectory::setStateAtIP ( const SteppingHelixStateInfo &  state )

inlineprivate

Definition at line 84 of file CachedTrajectory.h.

References stateAtIP_.

Referenced by propagateThoughFromIP().

{ stateAtIP_ = state; }

sign()

static int CachedTrajectory::sign ( float  number )

inlinestaticprotected

Definition at line 124 of file CachedTrajectory.h.

References contentValuesFiles::number.

Referenced by propagate().

                                             {
    if (number == 0)
      return 0;
    if (number > 0)
      return 1;
    else
      return -1;
  }

trajectoryDelta()

std::pair< float, float > CachedTrajectory::trajectoryDelta ( TrajectorType  trajectoryType )

private

calculate trajectory change (Theta,Phi) delta = final - original

Definition at line 256 of file CachedTrajectory.cc.

References delta(), ecalTrajectory_, FullTrajectory, fullTrajectory_, hcalTrajectory_, hoTrajectory_, IpToEcal, IpToHcal, IpToHO, SteppingHelixStateInfo::isValid(), SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::phi(), mps_fire::result, stateAtIP_, and PV3DBase< T, PVType, FrameType >::theta().

                                                                                    {
  // MEaning of trajectory change depends on its usage. In most cases we measure
  // change in a trajectory as difference between final track position and initial
  // direction. In some cases such as change of trajectory in the muon detector we
  // might want to compare theta-phi of two points or even find local maximum and
  // mimimum. In general it's not essential what defenition of the trajectory change
  // is used since we use these numbers only as a rough estimate on how much wider
  // we should make the preselection region.
  std::pair<float, float> result(0, 0);
  if (!stateAtIP_.isValid()) {
    edm::LogWarning("TrackAssociator") << "State at IP is not known set. Cannot estimate trajectory change. "
                                       << "Trajectory change is not taken into account in matching";
    return result;
  }
  switch (trajectoryType) {
    case IpToEcal:
      if (ecalTrajectory_.empty())
        edm::LogWarning("TrackAssociator") << "ECAL trajector is empty. Cannot estimate trajectory change. "
                                           << "Trajectory change is not taken into account in matching";
      else
        return delta(stateAtIP_.momentum().theta(),
                     ecalTrajectory_.front().position().theta(),
                     stateAtIP_.momentum().phi(),
                     ecalTrajectory_.front().position().phi());
      break;
    case IpToHcal:
      if (hcalTrajectory_.empty())
        edm::LogWarning("TrackAssociator") << "HCAL trajector is empty. Cannot estimate trajectory change. "
                                           << "Trajectory change is not taken into account in matching";
      else
        return delta(stateAtIP_.momentum().theta(),
                     hcalTrajectory_.front().position().theta(),
                     stateAtIP_.momentum().phi(),
                     hcalTrajectory_.front().position().phi());
      break;
    case IpToHO:
      if (hoTrajectory_.empty())
        edm::LogWarning("TrackAssociator") << "HO trajector is empty. Cannot estimate trajectory change. "
                                           << "Trajectory change is not taken into account in matching";
      else
        return delta(stateAtIP_.momentum().theta(),
                     hoTrajectory_.front().position().theta(),
                     stateAtIP_.momentum().phi(),
                     hoTrajectory_.front().position().phi());
      break;
    case FullTrajectory:
      if (fullTrajectory_.empty())
        edm::LogWarning("TrackAssociator") << "Full trajector is empty. Cannot estimate trajectory change. "
                                           << "Trajectory change is not taken into account in matching";
      else
        return delta(stateAtIP_.momentum().theta(),
                     fullTrajectory_.back().position().theta(),
                     stateAtIP_.momentum().phi(),
                     fullTrajectory_.back().position().phi());
      break;
    default:
      edm::LogWarning("TrackAssociator")
          << "Unkown or not supported trajector type. Cannot estimate trajectory change. "
          << "Trajectory change is not taken into account in matching";
  }
  return result;
}

Friends And Related Function Documentation

propagateThoughFromIP

std::vector<SteppingHelixStateInfo> propagateThoughFromIP ( const SteppingHelixStateInfo &  state, const Propagator *  ptr, const FiducialVolume &  volume, int  nsteps, float  step, float  minR, float  minZ, float  maxR, float  maxZ  )

friend

Definition at line 18 of file CachedTrajectory.cc.

                                                                      {
  CachedTrajectory neckLace;
  neckLace.setStateAtIP(state);
  neckLace.reset_trajectory();
  neckLace.setPropagator(prop);
  neckLace.setPropagationStep(0.1);
  neckLace.setMinDetectorRadius(minR);
  neckLace.setMinDetectorLength(minZ * 2.);
  neckLace.setMaxDetectorRadius(maxR);
  neckLace.setMaxDetectorLength(maxZ * 2.);

  // Propagate track
  bool isPropagationSuccessful = neckLace.propagateAll(state);

  if (!isPropagationSuccessful)
    return std::vector<SteppingHelixStateInfo>();

  std::vector<SteppingHelixStateInfo> complicatePoints;
  neckLace.getTrajectory(complicatePoints, volume, nsteps);

  return complicatePoints;
}

TrackDetectorAssociator

friend class TrackDetectorAssociator

friend

Definition at line 52 of file CachedTrajectory.h.

Member Data Documentation

ecalTrajectory_

std::vector<SteppingHelixStateInfo> CachedTrajectory::ecalTrajectory_

protected

Definition at line 143 of file CachedTrajectory.h.

Referenced by findEcalTrajectory(), getEcalTrajectory(), getStateAtEcal(), reset_trajectory(), and trajectoryDelta().

fullTrajectory_

std::deque<SteppingHelixStateInfo> CachedTrajectory::fullTrajectory_

protected

Definition at line 142 of file CachedTrajectory.h.

Referenced by distance(), getInnerState(), getOuterState(), getTrajectory(), propagate(), propagateAll(), reset_trajectory(), and trajectoryDelta().

fullTrajectoryFilled_

bool CachedTrajectory::fullTrajectoryFilled_

protected

Definition at line 152 of file CachedTrajectory.h.

Referenced by getTrajectory(), propagateAll(), and reset_trajectory().

hcalTrajectory_

std::vector<SteppingHelixStateInfo> CachedTrajectory::hcalTrajectory_

protected

Definition at line 144 of file CachedTrajectory.h.

Referenced by findHcalTrajectory(), getHcalTrajectory(), getStateAtHcal(), reset_trajectory(), and trajectoryDelta().

HOmaxRho_

float CachedTrajectory::HOmaxRho_

protected

Definition at line 158 of file CachedTrajectory.h.

Referenced by setMaxHORadius().

HOmaxZ_

float CachedTrajectory::HOmaxZ_

protected

Definition at line 159 of file CachedTrajectory.h.

Referenced by setMaxHOLength().

hoTrajectory_

std::vector<SteppingHelixStateInfo> CachedTrajectory::hoTrajectory_

protected

Definition at line 145 of file CachedTrajectory.h.

Referenced by findHOTrajectory(), getHOTrajectory(), getStateAtHO(), reset_trajectory(), and trajectoryDelta().

maxRho_

float CachedTrajectory::maxRho_

protected

Definition at line 156 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMaxDetectorRadius().

maxZ_

float CachedTrajectory::maxZ_

protected

Definition at line 157 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMaxDetectorLength().

minRho_

float CachedTrajectory::minRho_

protected

Definition at line 160 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMinDetectorRadius().

minZ_

float CachedTrajectory::minZ_

protected

Definition at line 161 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMinDetectorLength().

preshowerTrajectory_

std::vector<SteppingHelixStateInfo> CachedTrajectory::preshowerTrajectory_

protected

Definition at line 146 of file CachedTrajectory.h.

Referenced by findPreshowerTrajectory(), getPreshowerTrajectory(), getStateAtPreshower(), and reset_trajectory().

propagator_

const Propagator* CachedTrajectory::propagator_

protected

Definition at line 154 of file CachedTrajectory.h.

Referenced by propagate(), propagateAll(), and setPropagator().

stateAtIP_

SteppingHelixStateInfo CachedTrajectory::stateAtIP_

protected

Definition at line 150 of file CachedTrajectory.h.

Referenced by setStateAtIP(), and trajectoryDelta().

step_

float CachedTrajectory::step_

protected

Definition at line 162 of file CachedTrajectory.h.

Referenced by getPropagationStep(), propagateAll(), and setPropagationStep().

wideEcalTrajectory_

std::vector<GlobalPoint> CachedTrajectory::wideEcalTrajectory_

protected

Definition at line 147 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().

wideHcalTrajectory_

std::vector<GlobalPoint> CachedTrajectory::wideHcalTrajectory_

protected

Definition at line 148 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().

wideHOTrajectory_

std::vector<GlobalPoint> CachedTrajectory::wideHOTrajectory_

protected

Definition at line 149 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().