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 39 of file CachedTrajectory.h.

Member Enumeration Documentation

enum CachedTrajectory::TrajectorType

private

Enumerator
IpToEcal
IpToHcal
IpToHO
FullTrajectory

Definition at line 55 of file CachedTrajectory.h.

{ IpToEcal, IpToHcal, IpToHO, FullTrajectory };

enum CachedTrajectory::WideTrajectoryType

private

Enumerator
Ecal
Hcal
HO

Definition at line 56 of file CachedTrajectory.h.

{ Ecal, Hcal, HO };

Constructor & Destructor Documentation

CachedTrajectory::CachedTrajectory ( )

private

Definition at line 47 of file CachedTrajectory.cc.

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

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

Member Function Documentation

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

protected

Definition at line 332 of file CachedTrajectory.cc.

References M_PI, and query::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;
}

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

inlineprotected

Definition at line 131 of file CachedTrajectory.h.

References fullTrajectory_, cmsHarvester::index, 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());
  }

void CachedTrajectory::findEcalTrajectory ( const FiducialVolume &  volume )

private

Definition at line 456 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();
}

void CachedTrajectory::findHcalTrajectory ( const FiducialVolume &  volume )

private

Definition at line 476 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();
}

void CachedTrajectory::findHOTrajectory ( const FiducialVolume &  volume )

private

Definition at line 486 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();
}

void CachedTrajectory::findPreshowerTrajectory ( const FiducialVolume &  volume )

private

Definition at line 462 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();
}

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

(

)

const

Definition at line 468 of file CachedTrajectory.cc.

References ecalTrajectory_.

Referenced by MuonCaloDistanceProducer::fillDistanceMap().

                                                                                  {
   return ecalTrajectory_;
}

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

(

)

const

Definition at line 482 of file CachedTrajectory.cc.

References hcalTrajectory_.

Referenced by MuonCaloDistanceProducer::fillDistanceMap().

                                                                                  {
   return hcalTrajectory_;
}

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

(

)

const

Definition at line 492 of file CachedTrajectory.cc.

References hoTrajectory_.

Referenced by MuonCaloDistanceProducer::fillDistanceMap().

                                                                                 {
   return hoTrajectory_;
}

SteppingHelixStateInfo CachedTrajectory::getInnerState ( )

private

Definition at line 623 of file CachedTrajectory.cc.

References fullTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getOuterState ( )

private

Definition at line 631 of file CachedTrajectory.cc.

References fullTrajectory_.

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

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

(

)

const

Definition at line 472 of file CachedTrajectory.cc.

References preshowerTrajectory_.

Referenced by MuonCaloDistanceProducer::fillDistanceMap().

                                                                                       {
   return preshowerTrajectory_;
}

float CachedTrajectory::getPropagationStep ( ) const

inlineprivate

Definition at line 114 of file CachedTrajectory.h.

References step_.

{ return step_;}

SteppingHelixStateInfo CachedTrajectory::getStateAtEcal ( )

private

Definition at line 590 of file CachedTrajectory.cc.

References ecalTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtHcal ( )

private

Definition at line 606 of file CachedTrajectory.cc.

References hcalTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtHO ( )

private

Definition at line 614 of file CachedTrajectory.cc.

References hoTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtPreshower ( )

private

Definition at line 598 of file CachedTrajectory.cc.

References preshowerTrajectory_.

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

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 348 of file CachedTrajectory.cc.

References Reference_intrackfit_cff::barrel, Plane::build(), Cylinder::build(), PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, Reference_intrackfit_cff::endcap, Exception, fullTrajectory_, fullTrajectoryFilled_, i, FiducialVolume::isValid(), SteppingHelixStateInfo::isValid(), LogTrace, bookConverter::max, FiducialVolume::maxR(), FiducialVolume::maxZ(), FiducialVolume::minR(), FiducialVolume::minZ(), SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::perp(), SteppingHelixStateInfo::position(), position, propagate(), propagateForward(), relval_parameters_module::step, 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);
     }
}

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

private

Definition at line 497 of file CachedTrajectory.cc.

References funct::cos(), Ecal, Hcal, HO, cmsHarvester::index, LogTrace, M_PI, M_PI_2, phi, funct::pow(), dttmaxenums::r32, LocalError::rotate(), idealTransformation::rotation, LocalError::scale(), funct::sin(), mathSSE::sqrt(), filterCSVwithJSON::target, Surface::toGlobal(), wideEcalTrajectory_, wideHcalTrajectory_, wideHOTrajectory_, LocalError::xx(), LocalError::xy(), create_public_lumi_plots::xy, LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

                                                                           {
   std::vector<GlobalPoint>* wideTrajectory = 0;
   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 0;

   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;
}

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

private

Definition at line 78 of file CachedTrajectory.cc.

References cms::Exception::category(), cms::Exception::explainSelf(), TrajectoryStateOnSurface::freeState(), SteppingHelixStateInfo::getFreeState(), SteppingHelixStateInfo::momentum(), Propagator::propagate(), propagator_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

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()) );
     }
}

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

private

Definition at line 102 of file CachedTrajectory.cc.

References cms::Exception::category(), cms::Exception::explainSelf(), TrajectoryStateOnSurface::freeState(), SteppingHelixStateInfo::getFreeState(), SteppingHelixStateInfo::momentum(), Propagator::propagate(), propagator_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

{
   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()) );
     }
}

TrajectoryStateOnSurface CachedTrajectory::propagate ( const Plane *  plane )

private

get fast to a given DetId surface using cached trajectory

Definition at line 195 of file CachedTrajectory.cc.

References cms::Exception::category(), distance(), PV3DBase< T, PVType, FrameType >::eta(), cms::Exception::explainSelf(), fullTrajectory_, SteppingHelixStateInfo::getStateOnSurface(), LogTrace, SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), GloballyPositioned< T >::position(), Propagator::propagate(), propagator_, sign(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), 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);
     }
}

bool CachedTrajectory::propagateAll ( const SteppingHelixStateInfo &  initialState )

private

propagate through the whole detector, returns true if successful

Definition at line 126 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_==0) 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();
}

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

private

Definition at line 56 of file CachedTrajectory.cc.

References Plane::build(), distance(), SteppingHelixStateInfo::momentum(), SteppingHelixStateInfo::position(), funct::pow(), propagate(), dttmaxenums::r32, idealTransformation::rotation, mathSSE::sqrt(), filterCSVwithJSON::target, Vector3DBase< T, FrameTag >::unit(), and PV3DBase< T, PVType, FrameType >::z().

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);
}

void CachedTrajectory::reset_trajectory ( )

private

Definition at line 444 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;
}

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

inlineprivate

Definition at line 107 of file CachedTrajectory.h.

References cmsLHEtoEOSManager::l, and maxZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ maxZ_ = l/2.;}

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

inlineprivate

Definition at line 106 of file CachedTrajectory.h.

References maxRho_, and alignCSCRings::r.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ maxRho_ = r;}

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

inlineprivate

Definition at line 109 of file CachedTrajectory.h.

References HOmaxZ_, and cmsLHEtoEOSManager::l.

{ HOmaxZ_ = l/2.;}

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

inlineprivate

Definition at line 108 of file CachedTrajectory.h.

References HOmaxRho_, and alignCSCRings::r.

{ HOmaxRho_ = r;}

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

inlineprivate

Definition at line 111 of file CachedTrajectory.h.

References cmsLHEtoEOSManager::l, and minZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ minZ_ = l/2.;}

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

inlineprivate

Definition at line 110 of file CachedTrajectory.h.

References minRho_, and alignCSCRings::r.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ minRho_ = r;}

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

inlineprivate

Definition at line 113 of file CachedTrajectory.h.

References alignCSCRings::s, and step_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ step_ = s;}

void CachedTrajectory::setPropagator ( const Propagator *  ptr )

inlineprivate

Definition at line 74 of file CachedTrajectory.h.

References propagator_.

Referenced by propagateThoughFromIP().

{ propagator_ = ptr; }

void CachedTrajectory::setStateAtIP ( const SteppingHelixStateInfo &  state )

inlineprivate

Definition at line 75 of file CachedTrajectory.h.

References stateAtIP_.

Referenced by propagateThoughFromIP().

{ stateAtIP_ = state; }

static int CachedTrajectory::sign ( float  number )

inlinestaticprotected

Definition at line 118 of file CachedTrajectory.h.

References contentValuesFiles::number.

Referenced by propagate().

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

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

private

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

Definition at line 281 of file CachedTrajectory.cc.

References delta(), ecalTrajectory_, FullTrajectory, fullTrajectory_, hcalTrajectory_, hoTrajectory_, IpToEcal, IpToHcal, IpToHO, SteppingHelixStateInfo::isValid(), SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::phi(), query::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

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 19 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;

}

friend class TrackDetectorAssociator

friend

Definition at line 48 of file CachedTrajectory.h.

Member Data Documentation

std::vector<SteppingHelixStateInfo> CachedTrajectory::ecalTrajectory_

protected

Definition at line 137 of file CachedTrajectory.h.

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

std::deque<SteppingHelixStateInfo> CachedTrajectory::fullTrajectory_

protected

Definition at line 136 of file CachedTrajectory.h.

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

bool CachedTrajectory::fullTrajectoryFilled_

protected

Definition at line 146 of file CachedTrajectory.h.

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

std::vector<SteppingHelixStateInfo> CachedTrajectory::hcalTrajectory_

protected

Definition at line 138 of file CachedTrajectory.h.

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

float CachedTrajectory::HOmaxRho_

protected

Definition at line 152 of file CachedTrajectory.h.

Referenced by setMaxHORadius().

float CachedTrajectory::HOmaxZ_

protected

Definition at line 153 of file CachedTrajectory.h.

Referenced by setMaxHOLength().

std::vector<SteppingHelixStateInfo> CachedTrajectory::hoTrajectory_

protected

Definition at line 139 of file CachedTrajectory.h.

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

float CachedTrajectory::maxRho_

protected

Definition at line 150 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMaxDetectorRadius().

float CachedTrajectory::maxZ_

protected

Definition at line 151 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMaxDetectorLength().

float CachedTrajectory::minRho_

protected

Definition at line 154 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMinDetectorRadius().

float CachedTrajectory::minZ_

protected

Definition at line 155 of file CachedTrajectory.h.

Referenced by propagateAll(), and setMinDetectorLength().

std::vector<SteppingHelixStateInfo> CachedTrajectory::preshowerTrajectory_

protected

Definition at line 140 of file CachedTrajectory.h.

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

const Propagator* CachedTrajectory::propagator_

protected

Definition at line 148 of file CachedTrajectory.h.

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

SteppingHelixStateInfo CachedTrajectory::stateAtIP_

protected

Definition at line 144 of file CachedTrajectory.h.

Referenced by setStateAtIP(), and trajectoryDelta().

float CachedTrajectory::step_

protected

Definition at line 156 of file CachedTrajectory.h.

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

std::vector<GlobalPoint> CachedTrajectory::wideEcalTrajectory_

protected

Definition at line 141 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().

std::vector<GlobalPoint> CachedTrajectory::wideHcalTrajectory_

protected

Definition at line 142 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().

std::vector<GlobalPoint> CachedTrajectory::wideHOTrajectory_

protected

Definition at line 143 of file CachedTrajectory.h.

Referenced by getWideTrajectory(), and reset_trajectory().