CMSSW_4_4_3_patch1/src/RecoEgamma/EgammaPhotonAlgos/src/ConversionFastHelix.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionFastHelix.h"
#include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
#include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
#include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
#include "TrackingTools/TrajectoryParametrization/interface/CartesianTrajectoryError.h"


#include <cfloat>

ConversionFastHelix::ConversionFastHelix(const GlobalPoint& outerHit,
                     const GlobalPoint& middleHit,
                                         const GlobalPoint& aVertex,
                                         const MagneticField* field ) : 
  theOuterHit(outerHit),
  theMiddleHit(middleHit),
  theVertex(aVertex),
  theCircle(outerHit,
            middleHit,
            aVertex),
  mField(field) {
  
  validStateAtVertex=false;


  makeHelix();
 
  
}


void ConversionFastHelix::makeHelix()  {


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

}


FreeTrajectoryState ConversionFastHelix::stateAtVertex()  {

  return theHelix_;

}


FreeTrajectoryState ConversionFastHelix::helixStateAtVertex()  {
  
  
  
  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;
  }
  


  
  
}

FreeTrajectoryState ConversionFastHelix::straightLineStateAtVertex() {

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

}