EgammaPhotonAlgos/src/ConversionFastHelix.cc

0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionFastHelix.h"
0002 #include "RecoTracker/TkSeedGenerator/interface/FastLine.h"
0003 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0004 #include "TrackingTools/TrajectoryParametrization/interface/CurvilinearTrajectoryError.h"
0005 #include "TrackingTools/TrajectoryParametrization/interface/CartesianTrajectoryError.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007
0008 #include <cfloat>
0009
0010 ConversionFastHelix::ConversionFastHelix(const GlobalPoint& outerHit,
0011                                          const GlobalPoint& middleHit,
0012                                          const GlobalPoint& aVertex,
0013                                          const MagneticField* field)
0014     : theOuterHit(outerHit),
0015       theMiddleHit(middleHit),
0016       theVertex(aVertex),
0017       theCircle(outerHit, middleHit, aVertex),
0018       mField(field) {
0019   validStateAtVertex = false;
0020
0021   makeHelix();
0022 }
0023
0024 void ConversionFastHelix::makeHelix() {
0025   if (theCircle.isValid()) {
0026     theHelix_ = helixStateAtVertex();
0027   } else {
0028     theHelix_ = straightLineStateAtVertex();
0029   }
0030 }
0031
0032 FreeTrajectoryState ConversionFastHelix::stateAtVertex() { return theHelix_; }
0033
0034 FreeTrajectoryState ConversionFastHelix::helixStateAtVertex() {
0035   GlobalPoint pMid(theMiddleHit);
0036   GlobalPoint v(theVertex);
0037   FTS atVertex;
0038
0039   double dydx = 0.;
0040   double pt = 0., px = 0., py = 0.;
0041
0042   //remember (radius rho in cm):
0043   //rho =
0044   //100. * pt *
0045   //(10./(3.*MagneticField::inTesla(GlobalPoint(0., 0., 0.)).z()));
0046
0047   double rho = theCircle.rho();
0048   pt = 0.01 * rho * (0.3 * mField->inTesla(GlobalPoint(0, 0, 0)).z());
0049
0050   // (py/px)|x=v.x() = (dy/dx)|x=v.x()
0051   //remember:
0052   //y(x) = +-sqrt(rho^2 - (x-x0)^2) + y0
0053   //y(x) =  sqrt(rho^2 - (x-x0)^2) + y0  if y(x) >= y0
0054   //y(x) = -sqrt(rho^2 - (x-x0)^2) + y0  if y(x) < y0
0055   //=> (dy/dx) = -(x-x0)/sqrt(Q)  if y(x) >= y0
0056   //   (dy/dx) =  (x-x0)/sqrt(Q)  if y(x) < y0
0057   //with Q = rho^2 - (x-x0)^2
0058
0059   double arg = rho * rho - ((v.x() - theCircle.x0()) * (v.x() - theCircle.x0()));
0060
0061   if (arg > 0) {
0062     //  double root = sqrt(  rho*rho - ( (v.x()-theCircle.x0())*(v.x()-theCircle.x0()) )  );
0063     double root = sqrt(arg);
0064
0065     if ((v.y() - theCircle.y0()) > 0.)
0066       dydx = -(v.x() - theCircle.x0()) / root;
0067     else
0068       dydx = (v.x() - theCircle.x0()) / root;
0069
0070     px = pt / sqrt(1. + dydx * dydx);
0071     py = px * dydx;
0072     // check sign with scalar product
0073     if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
0074       px *= -1.;
0075       py *= -1.;
0076     }
0077
0078     //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;
0079     //calculate z0, pz
0080     //(z, R*phi) linear relation in a helix
0081     //with R, phi defined as radius and angle w.r.t. centre of circle
0082     //in transverse plane
0083     //pz = pT*(dz/d(R*phi)))
0084
0085     FastLine flfit(theOuterHit, theMiddleHit, theCircle.rho());
0086
0087     double z_0 = 0;
0088
0089     //std::cout << " ConversionFastHelix:helixStateAtVertex  flfit.n2() " <<  flfit.n2() << " flfit.c() " << flfit.c() << " flfit.n2() " << flfit.n2() << std::endl;
0090     if (flfit.n2() != 0 && !edm::isNotFinite(flfit.c()) && !edm::isNotFinite(flfit.n2())) {
0091       //  std::cout << " Accepted " << std::endl;
0092       z_0 = -flfit.c() / flfit.n2();
0093       double dzdrphi = -flfit.n1() / flfit.n2();
0094       double pz = pt * dzdrphi;
0095
0096       //get sign of particle
0097
0098       GlobalVector magvtx = mField->inTesla(v);
0099       TrackCharge q = ((theCircle.x0() * py - theCircle.y0() * px) / (magvtx.z()) < 0.) ? -1 : 1;
0100
0101       AlgebraicSymMatrix55 C = AlgebraicMatrixID();
0102       //MP
0103
0104       atVertex = FTS(GlobalTrajectoryParameters(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, mField),
0105                      CurvilinearTrajectoryError(C));
0106
0107       //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;
0108       //std::cout << " ConversionFastHelix:helixStateAtVertex atVertex.transverseCurvature() " << atVertex.transverseCurvature() << std::endl;
0109       if (atVertex.transverseCurvature() != 0) {
0110         validStateAtVertex = true;
0111
0112         //std::cout << " ConversionFastHelix:helixStateAtVertex validHelixStateAtVertex status " << validStateAtVertex << std::endl;
0113         return atVertex;
0114       } else
0115         return atVertex;
0116     } else {
0117       //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted  validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
0118       return atVertex;
0119     }
0120
0121   } else {
0122     //std::cout << " ConversionFastHelix:helixStateAtVertex not accepted because arg <0 validHelixStateAtVertex status  " << validStateAtVertex << std::endl;
0123     return atVertex;
0124   }
0125 }
0126
0127 FreeTrajectoryState ConversionFastHelix::straightLineStateAtVertex() {
0128   FTS atVertex;
0129
0130   //calculate FTS assuming straight line...
0131
0132   GlobalPoint pMid(theMiddleHit);
0133   GlobalPoint v(theVertex);
0134
0135   double dydx = 0.;
0136   double pt = 0., px = 0., py = 0.;
0137
0138   if (fabs(theCircle.n1()) > 0. || fabs(theCircle.n2()) > 0.)
0139     pt = 1.e+4;  // 10 TeV //else no pt
0140   if (fabs(theCircle.n2()) > 0.) {
0141     dydx = -theCircle.n1() / theCircle.n2();  //else px = 0
0142   }
0143
0144   if (pt == 0 && dydx == 0.) {
0145     validStateAtVertex = false;
0146     return atVertex;
0147   }
0148   px = pt / sqrt(1. + dydx * dydx);
0149   py = px * dydx;
0150
0151   // check sign with scalar product
0152   if (px * (pMid.x() - v.x()) + py * (pMid.y() - v.y()) < 0.) {
0153     px *= -1.;
0154     py *= -1.;
0155   }
0156
0157   //calculate z_0 and pz at vertex using weighted mean
0158   //z = z(r) = z0 + (dz/dr)*r
0159   //tan(theta) = dr/dz = (dz/dr)^-1
0160   //theta = atan(1./dzdr)
0161   //p = pt/sin(theta)
0162   //pz = p*cos(theta) = pt/tan(theta)
0163
0164   FastLine flfit(theOuterHit, theMiddleHit);
0165
0166   double z_0 = 0;
0167   if (flfit.n2() != 0 && !edm::isNotFinite(flfit.c()) && !edm::isNotFinite(flfit.n2())) {
0168     z_0 = -flfit.c() / flfit.n2();
0169
0170     double dzdr = -flfit.n1() / flfit.n2();
0171     double pz = pt * dzdr;
0172
0173     TrackCharge q = 1;
0174
0175     atVertex = FTS(GlobalPoint(v.x(), v.y(), z_0), GlobalVector(px, py, pz), q, mField);
0176
0177     //    std::cout << "  ConversionFastHelix::straightLineStateAtVertex curvature " << atVertex.transverseCurvature() << "   signedInverseMomentum " << atVertex.signedInverseMomentum() << std::endl;
0178     if (atVertex.transverseCurvature() == -0) {
0179       return atVertex;
0180     } else {
0181       validStateAtVertex = true;
0182       return atVertex;
0183     }
0184
0185   } else {
0186     return atVertex;
0187   }
0188 }