#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.

{}

Member Function Documentation

FreeTrajectoryState ConversionFastHelix::helixStateAtVertex ( )

Definition at line 51 of file ConversionFastHelix.cc.

References cmsCodeRulesChecker::arg, funct::C, FastLine::c(), MagneticField::inTesla(), edm::detail::isnan(), mField, lumiQueryAPI::q, FastCircle::rho(), 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.

{return validStateAtVertex;  }

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 funct::C, FastLine::c(), edm::detail::isnan(), mField, FastCircle::n1(), FastLine::n1(), FastCircle::n2(), FastLine::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;
 
    AlgebraicSymMatrix66 C = AlgebraicMatrixID();
    //MP
    atVertex = FTS(GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0),
                                                  GlobalVector(px, py, pz),
                                                  q,
                                                  mField),
                       CartesianTrajectoryError(C));

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

}