#include <ConversionFastHelix.h>

Public Member Functions
	ConversionFastHelix (const GlobalPoint &outerHit, const GlobalPoint &middleHit, const GlobalPoint &aVertex, const MagneticField *field)

FTS	helixStateAtVertex ()

bool	isValid ()

void	makeHelix ()

FTS	stateAtVertex ()

FTS	straightLineStateAtVertex ()

	~ConversionFastHelix ()

Private Types
typedef FreeTrajectoryState	FTS

Private Attributes
const MagneticField *	mField

FastCircle	theCircle

FTS	theHelix_

GlobalPoint	theMiddleHit

GlobalPoint	theOuterHit

GlobalPoint	theVertex

bool	validStateAtVertex

Detailed Description

Generation of track parameters at a vertex using two hits and a vertex.

Definition at line 16 of file ConversionFastHelix.h.

Member Typedef Documentation

typedef FreeTrajectoryState ConversionFastHelix::FTS

private

Definition at line 20 of file ConversionFastHelix.h.

Constructor & Destructor Documentation

ConversionFastHelix::ConversionFastHelix	(	const GlobalPoint &	outerHit,
		const GlobalPoint &	middleHit,
		const GlobalPoint &	aVertex,
		const MagneticField *	field
	)

Definition at line 10 of file ConversionFastHelix.cc.

References makeHelix(), and validStateAtVertex.

                                                   : 
   theOuterHit(outerHit),
   theMiddleHit(middleHit),
   theVertex(aVertex),
   theCircle(outerHit,
         middleHit,
         aVertex),
   mField(field) {
   
   validStateAtVertex=false;
 
 
   makeHelix();
  
   
 }

ConversionFastHelix::~ConversionFastHelix ( )

inline

Definition at line 30 of file ConversionFastHelix.h.

30 {}

Member Function Documentation

FreeTrajectoryState ConversionFastHelix::helixStateAtVertex ( )

Definition at line 51 of file ConversionFastHelix.cc.

References FastLine::c(), funct::C, MagneticField::inTesla(), edm::detail::isnan(), mField, lumiQueryAPI::q, rho, FastCircle::rho(), dbtoconf::root, mathSSE::sqrt(), theCircle, theMiddleHit, theOuterHit, theVertex, FreeTrajectoryState::transverseCurvature(), v, validStateAtVertex, PV3DBase< T, PVType, FrameType >::x(), FastCircle::x0(), PV3DBase< T, PVType, FrameType >::y(), FastCircle::y0(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by makeHelix().

                                                              {
   
   
   
   GlobalPoint pMid(theMiddleHit);
   GlobalPoint v(theVertex);
   FTS atVertex;
   
   double dydx = 0.;
   double pt = 0., px = 0., py = 0.;
   
   //remember (radius rho in cm):
   //rho = 
   //100. * pt * 
   //(10./(3.*MagneticField::inTesla(GlobalPoint(0., 0., 0.)).z()));
   
   double rho = theCircle.rho();
   pt = 0.01 * rho * (0.3*mField->inTesla(GlobalPoint(0,0,0)).z());
   
   // (py/px)|x=v.x() = (dy/dx)|x=v.x()
   //remember:
   //y(x) = +-sqrt(rho^2 - (x-x0)^2) + y0 
   //y(x) =  sqrt(rho^2 - (x-x0)^2) + y0  if y(x) >= y0 
   //y(x) = -sqrt(rho^2 - (x-x0)^2) + y0  if y(x) < y0
   //=> (dy/dx) = -(x-x0)/sqrt(Q)  if y(x) >= y0
   //   (dy/dx) =  (x-x0)/sqrt(Q)  if y(x) < y0
   //with Q = rho^2 - (x-x0)^2
   
   
   double arg=rho*rho - ( (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) );
   
   if ( arg >= 0 ) { 
     
     
     //  double root = sqrt(  rho*rho - ( (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) )  );
     double root = sqrt(  arg );
     
     if((v.y() - theCircle.y0()) > 0.)
       dydx = -(v.x() - theCircle.x0()) / root;
     else
       dydx = (v.x() - theCircle.x0()) / root;
     
     px = pt/sqrt(1. + dydx*dydx);
     py = px*dydx;
     // check sign with scalar product
     if(px*(pMid.x() - v.x()) + py*(pMid.y() - v.y()) < 0.) {
       px *= -1.;
       py *= -1.;
     } 
     
     //std::cout << " ConversionFastHelix:helixStateAtVertex  rho " << rho  << " pt " << pt  << " v " <<  v << " theCircle.x0() " <<theCircle.x0() << " theCircle.y0() "  <<  theCircle.y0() << " v.x()-theCircle.x0() "  << v.x()-theCircle.x0() << " rho^2 " << rho*rho << "  v.x()-theCircle.x0()^2 " <<   (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) <<  " root " << root << " arg " << arg <<  " dydx " << dydx << std::endl;
     //calculate z0, pz
     //(z, R*phi) linear relation in a helix
     //with R, phi defined as radius and angle w.r.t. centre of circle
     //in transverse plane
     //pz = pT*(dz/d(R*phi)))
     
     FastLine flfit(theOuterHit, theMiddleHit, theCircle.rho());
    
     
     
     double z_0 = 0; 
     
     //std::cout << " ConversionFastHelix:helixStateAtVertex  flfit.n2() " <<  flfit.n2() << " flfit.c() " << flfit.c() << " flfit.n2() " << flfit.n2() << std::endl;
     if ( flfit.n2() !=0 && !std::isnan( flfit.c()) && !std::isnan(flfit.n2())   ) {
       //  std::cout << " Accepted " << std::endl;
       z_0 = -flfit.c()/flfit.n2();
       double dzdrphi = -flfit.n1()/flfit.n2();
       double pz = pt*dzdrphi;
       
       //get sign of particle
       
       GlobalVector magvtx=mField->inTesla(v);
       TrackCharge q = 
     ((theCircle.x0()*py - theCircle.y0()*px) / 
      (magvtx.z()) < 0.) ? 
     -1 : 1;
       
 
       AlgebraicSymMatrix55 C = AlgebraicMatrixID();
       //MP
       
       atVertex = FTS(GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0),
                         GlobalVector(px, py, pz),
                         q,
                         mField),
              CurvilinearTrajectoryError(C));
       
       //std::cout << " ConversionFastHelix:helixStateAtVertex globalPoint " << GlobalPoint(v.x(), v.y(), z_0) << " GlobalVector " << GlobalVector(px, py, pz)  << " q " << q << " MField " << mField->inTesla(v) << std::endl;
       //std::cout << " ConversionFastHelix:helixStateAtVertex atVertex.transverseCurvature() " << atVertex.transverseCurvature() << std::endl;
       if( atVertex.transverseCurvature() !=0 ) {
     
     validStateAtVertex=true;    
     
     //std::cout << " ConversionFastHelix:helixStateAtVertex validHelixStateAtVertex status " << validStateAtVertex << std::endl;
     return atVertex;
       }else
     return atVertex;
     } else {
       //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted  validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
       return atVertex;
     }
     
     
     
   } else {
     
     //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted because arg <0 validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
     return atVertex;
   }
   
 
 
   
   
 }

bool ConversionFastHelix::isValid ( void )

inline

Definition at line 35 of file ConversionFastHelix.h.

References validStateAtVertex.

35 {return validStateAtVertex; }

ConversionFastHelix::validStateAtVertex

bool validStateAtVertex

Definition: ConversionFastHelix.h:46

void ConversionFastHelix::makeHelix ( )

Definition at line 31 of file ConversionFastHelix.cc.

References helixStateAtVertex(), FastCircle::isValid(), straightLineStateAtVertex(), theCircle, and theHelix_.

Referenced by ConversionFastHelix().

                                      {
 
 
   if(   theCircle.isValid()) {
      theHelix_=helixStateAtVertex();
   } else {
      theHelix_= straightLineStateAtVertex();
   }
   
 
 }

FreeTrajectoryState ConversionFastHelix::stateAtVertex ( )

Definition at line 44 of file ConversionFastHelix.cc.

References theHelix_.

Referenced by InOutConversionSeedFinder::findSeeds(), and InOutConversionSeedFinder::startSeed().

                                                         {
 
   return theHelix_;
 
 }

FreeTrajectoryState ConversionFastHelix::straightLineStateAtVertex ( )

Definition at line 168 of file ConversionFastHelix.cc.

References FastLine::c(), edm::detail::isnan(), mField, FastLine::n1(), FastCircle::n1(), FastLine::n2(), FastCircle::n2(), lumiQueryAPI::q, mathSSE::sqrt(), theCircle, theMiddleHit, theOuterHit, theVertex, FreeTrajectoryState::transverseCurvature(), v, validStateAtVertex, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

Referenced by makeHelix().

                                                                    {
 
   FTS atVertex;
 
   //calculate FTS assuming straight line...
 
   GlobalPoint pMid(theMiddleHit);
   GlobalPoint v(theVertex);
 
   double dydx = 0.;
   double pt = 0., px = 0., py = 0.;
   
   if(fabs(theCircle.n1()) > 0. || fabs(theCircle.n2()) > 0.)
     pt = 1.e+4;// 10 TeV //else no pt
   if(fabs(theCircle.n2()) > 0.) {
     dydx = -theCircle.n1()/theCircle.n2(); //else px = 0 
   }
 
   if ( pt==0 && dydx==0. ) {
     validStateAtVertex=false;
     return atVertex; 
   }
   px = pt/sqrt(1. + dydx*dydx);
   py = px*dydx;
 
   // check sign with scalar product
   if (px*(pMid.x() - v.x()) + py*(pMid.y() - v.y()) < 0.) {
     px *= -1.;
     py *= -1.;
   } 
 
   //calculate z_0 and pz at vertex using weighted mean
   //z = z(r) = z0 + (dz/dr)*r
   //tan(theta) = dr/dz = (dz/dr)^-1
   //theta = atan(1./dzdr)
   //p = pt/sin(theta)
   //pz = p*cos(theta) = pt/tan(theta) 
 
 
   FastLine flfit(theOuterHit, theMiddleHit);
 
   double z_0 = 0;
   if (flfit.n2() !=0  && !std::isnan( flfit.c()) && !std::isnan(flfit.n2())   ) {
     z_0 = -flfit.c()/flfit.n2();
 
     double dzdr = -flfit.n1()/flfit.n2();
     double pz = pt*dzdr; 
     
     TrackCharge q = 1;
  
     atVertex = FTS(GlobalPoint(v.x(), v.y(), z_0),
            GlobalVector(px, py, pz),
            q,
            mField
            );
 
     //    std::cout << "  ConversionFastHelix::straightLineStateAtVertex curvature " << atVertex.transverseCurvature() << "   signedInverseMomentum " << atVertex.signedInverseMomentum() << std::endl;
     if ( atVertex.transverseCurvature() == -0 ) {
       return atVertex;
     } else {
       validStateAtVertex=true;        
       return atVertex;
     }
   
   } else {
 
 
     return atVertex;
   
   }
 
 }