#include <CachedTrajectory.h>

Protected Member Functions
std::pair< float, float >	delta (const double &theta1, const double &theta2, const double &phi1, const double &phi2) dso_internal
float	distance (const Plane *plane, int index) dso_internal
Static Protected Member Functions
static int	sign (float number) dso_internal
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 &) dso_internal
void	findHcalTrajectory (const FiducialVolume &) dso_internal
void	findHOTrajectory (const FiducialVolume &) dso_internal
void	findPreshowerTrajectory (const FiducialVolume &) dso_internal
const std::vector < SteppingHelixStateInfo > &	getEcalTrajectory () dso_internal
const std::vector < SteppingHelixStateInfo > &	getHcalTrajectory () dso_internal
const std::vector < SteppingHelixStateInfo > &	getHOTrajectory () dso_internal
SteppingHelixStateInfo	getInnerState () dso_internal
SteppingHelixStateInfo	getOuterState () dso_internal
const std::vector < SteppingHelixStateInfo > &	getPreshowerTrajectory () dso_internal
float	getPropagationStep () const dso_internal
SteppingHelixStateInfo	getStateAtEcal () dso_internal
SteppingHelixStateInfo	getStateAtHcal () dso_internal
SteppingHelixStateInfo	getStateAtHO () dso_internal
SteppingHelixStateInfo	getStateAtPreshower () dso_internal
void	getTrajectory (std::vector< SteppingHelixStateInfo > &, const FiducialVolume &, int steps=4) dso_internal
std::vector< GlobalPoint > *	getWideTrajectory (const std::vector< SteppingHelixStateInfo > &, WideTrajectoryType) dso_internal
void	propagate (SteppingHelixStateInfo &state, const Plane &plane) dso_internal
void	propagate (SteppingHelixStateInfo &state, const Cylinder &cylinder) dso_internal
TrajectoryStateOnSurface	propagate (const Plane *plane) dso_internal
	get fast to a given DetId surface using cached trajectory
bool	propagateAll (const SteppingHelixStateInfo &initialState) dso_internal
	propagate through the whole detector, returns true if successful
void	propagateForward (SteppingHelixStateInfo &state, float distance) dso_internal
void	reset_trajectory () dso_internal
void	setMaxDetectorLength (float l=2200.) dso_internal
void	setMaxDetectorRadius (float r=800.) dso_internal
void	setMaxHOLength (float l=2200.) dso_internal
void	setMaxHORadius (float r=800.) dso_internal
void	setMinDetectorLength (float l=0.) dso_internal
void	setMinDetectorRadius (float r=0.) dso_internal
void	setPropagationStep (float s=20.)
void	setPropagator (const Propagator *ptr) dso_internal
void	setStateAtIP (const SteppingHelixStateInfo &state) dso_internal
std::pair< float, float >	trajectoryDelta (TrajectorType) dso_internal
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 41 of file CachedTrajectory.h.

Member Enumeration Documentation

enum CachedTrajectory::TrajectorType [private]

Enumerator:

IpToEcal
IpToHcal
IpToHO
FullTrajectory

Definition at line 50 of file CachedTrajectory.h.

{ IpToEcal, IpToHcal, IpToHO, FullTrajectory };

enum CachedTrajectory::WideTrajectoryType [private]

Enumerator:

Ecal
Hcal
HO

Definition at line 51 of file CachedTrajectory.h.

{ Ecal, Hcal, HO };

Constructor & Destructor Documentation

CachedTrajectory::CachedTrajectory ( ) [private]

Definition at line 48 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 333 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
	)		`[inline, protected]`

Definition at line 131 of file CachedTrajectory.h.

References fullTrajectory_, getHLTprescales::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 457 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 477 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 487 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 463 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 ( ) [private]

Definition at line 469 of file CachedTrajectory.cc.

References ecalTrajectory_.

                                                                             {
   return ecalTrajectory_;
}

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getHcalTrajectory ( ) [private]

Definition at line 483 of file CachedTrajectory.cc.

References hcalTrajectory_.

                                                                             {
   return hcalTrajectory_;
}

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getHOTrajectory ( ) [private]

Definition at line 493 of file CachedTrajectory.cc.

References hoTrajectory_.

                                                                           {
   return hoTrajectory_;
}

SteppingHelixStateInfo CachedTrajectory::getInnerState ( ) [private]

Definition at line 624 of file CachedTrajectory.cc.

References fullTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getOuterState ( ) [private]

Definition at line 632 of file CachedTrajectory.cc.

References fullTrajectory_.

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

const std::vector< SteppingHelixStateInfo > & CachedTrajectory::getPreshowerTrajectory ( ) [private]

Definition at line 473 of file CachedTrajectory.cc.

References preshowerTrajectory_.

                                                                                  {
   return preshowerTrajectory_;
}

float CachedTrajectory::getPropagationStep ( ) const [inline, private]

Definition at line 114 of file CachedTrajectory.h.

References step_.

{ return step_;}

SteppingHelixStateInfo CachedTrajectory::getStateAtEcal ( ) [private]

Definition at line 591 of file CachedTrajectory.cc.

References ecalTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtHcal ( ) [private]

Definition at line 607 of file CachedTrajectory.cc.

References hcalTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtHO ( ) [private]

Definition at line 615 of file CachedTrajectory.cc.

References hoTrajectory_.

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

SteppingHelixStateInfo CachedTrajectory::getStateAtPreshower ( ) [private]

Definition at line 599 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 349 of file CachedTrajectory.cc.

References Reference_intrackfit_cff::barrel, newFWLiteAna::build, Reference_intrackfit_cff::endcap, Exception, fullTrajectory_, fullTrajectoryFilled_, i, FiducialVolume::isValid(), SteppingHelixStateInfo::isValid(), LogTrace, max(), FiducialVolume::maxR(), FiducialVolume::maxZ(), FiducialVolume::minR(), FiducialVolume::minZ(), SteppingHelixStateInfo::momentum(), PV3DBase< T, PVType, FrameType >::perp(), SteppingHelixStateInfo::position(), position, propagate(), propagateForward(), launcher::step, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and 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( Cylinder::PositionType (0, 0, 0), Cylinder::RotationType (), volume.minR());
         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 498 of file CachedTrajectory.cc.

References funct::cos(), Ecal, Hcal, HO, getHLTprescales::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 79 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 {
           state = shp->propagate(state, plane);
        }
        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 103 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 {
           state = shp->propagate(state, cylinder);
        }
        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 196 of file CachedTrajectory.cc.

References cms::Exception::category(), distance(), cms::Exception::explainSelf(), fullTrajectory_, SteppingHelixStateInfo::getStateOnSurface(), LogTrace, SteppingHelixStateInfo::momentum(), GloballyPositioned< T >::position(), Propagator::propagate(), propagator_, sign(), evf::utils::state, 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 { 
           state = shp->propagate(fullTrajectory_[closestPointOnLeft], *plane);
        }
        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 127 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 57 of file CachedTrajectory.cc.

References newFWLiteAna::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 445 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. ) [inline, private]

Definition at line 107 of file CachedTrajectory.h.

References prof2calltree::l, and maxZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ maxZ_ = l/2.;}

void CachedTrajectory::setMaxDetectorRadius ( float r = 800. ) [inline, private]

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. ) [inline, private]

Definition at line 109 of file CachedTrajectory.h.

References HOmaxZ_, and prof2calltree::l.

{ HOmaxZ_ = l/2.;}

void CachedTrajectory::setMaxHORadius ( float r = 800. ) [inline, private]

Definition at line 108 of file CachedTrajectory.h.

References HOmaxRho_, and alignCSCRings::r.

{ HOmaxRho_ = r;}

void CachedTrajectory::setMinDetectorLength ( float l = 0. ) [inline, private]

Definition at line 111 of file CachedTrajectory.h.

References prof2calltree::l, and minZ_.

Referenced by CachedTrajectory(), and propagateThoughFromIP().

{ minZ_ = l/2.;}

void CachedTrajectory::setMinDetectorRadius ( float r = 0. ) [inline, private]

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. ) [inline, private]

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 ) [inline, private]

Definition at line 69 of file CachedTrajectory.h.

References propagator_.

Referenced by propagateThoughFromIP().

{ propagator_ = ptr; }

void CachedTrajectory::setStateAtIP ( const SteppingHelixStateInfo & state ) [inline, private]

Definition at line 70 of file CachedTrajectory.h.

References evf::utils::state, and stateAtIP_.

Referenced by propagateThoughFromIP().

{ stateAtIP_ = state; }

static int CachedTrajectory::sign ( float number ) [inline, static, protected]

Definition at line 118 of file CachedTrajectory.h.

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 282 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 20 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;

}