/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoMuon/TrackerSeedGenerator/src/L1MuonPixelTrackFitter.cc

Go to the documentation of this file.

#include "RecoMuon/TrackerSeedGenerator/interface/L1MuonPixelTrackFitter.h"

#include "DataFormats/L1GlobalMuonTrigger/interface/L1MuGMTCand.h"

#include "RecoTracker/TkMSParametrization/interface/PixelRecoLineRZ.h"
#include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
#include "RecoPixelVertexing/PixelTrackFitting/interface/PixelTrackBuilder.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "RecoTracker/TkSeedingLayers/interface/SeedingHitSet.h"

#include <math.h>

template <class T> T sqr( T t) {return t*t;}

L1MuonPixelTrackFitter::L1MuonPixelTrackFitter(const edm::ParameterSet& cfg)
  : theConfig(cfg)
{ }

void L1MuonPixelTrackFitter::setL1Constraint(const L1MuGMTCand & muon)
{
  thePhiL1 = muon.phiValue()+0.021817;
  theEtaL1 = muon.etaValue();
  theChargeL1 = muon.charge();
}

void L1MuonPixelTrackFitter::setPxConstraint(const SeedingHitSet & hits)
{
  theHit1 = hits[0]->globalPosition();
  theHit1 = hits[1]->globalPosition();
}

reco::Track* L1MuonPixelTrackFitter::run( const edm::EventSetup& es, 
    const std::vector<const TrackingRecHit *>& hits, const TrackingRegion& region) const
{
  edm::ESHandle<MagneticField> fieldESH;
  es.get<IdealMagneticFieldRecord>().get(fieldESH);

  double phi_vtx_fromHits = (theHit2-theHit1).phi();

  double invPt    = valInversePt(phi_vtx_fromHits, thePhiL1, theEtaL1);
  double invPtErr = errInversePt(invPt, theEtaL1);

  int charge = (invPt > 0.) ? 1 : -1;
  double  curvature = PixelRecoUtilities::curvature(invPt, es);

  Circle circle(theHit1, theHit2, curvature);

  double valPhi = this->valPhi(circle, charge);
  double valTip = this->valTip(circle, curvature); 
  double valZip = this->valZip(curvature, theHit1,theHit2);
  double valCotTheta = this->valCotTheta(PixelRecoLineRZ(theHit1,theHit2)); 

  static double invPtErrorScale = theConfig.getParameter<double>("invPtErrorScale");
  static double phiErrorScale= theConfig.getParameter<double>("phiErrorScale");
  static double cotThetaErrorScale = theConfig.getParameter<double>("cotThetaErrorScale");
  static double tipErrorScale = theConfig.getParameter<double>("tipErrorScale");
  static double zipErrorScale = theConfig.getParameter<double>("zipErrorScale");

//  if ( (fabs(invPt)-0.1)/invPtErr > 3.) return 0;

  PixelTrackBuilder builder;
  double valPt =  (fabs(invPt) > 1.e-4) ? 1./fabs(invPt) : 1.e4;
  double errPt = invPtErrorScale*invPtErr * valPt*valPt;
  double eta = asinh(valCotTheta);
  double errPhi = this->errPhi(invPt, eta)           * phiErrorScale;
  double errCotTheta = this->errCotTheta(invPt, eta) * cotThetaErrorScale;
  double errTip = this->errTip(invPt, eta)           * tipErrorScale; 
  double errZip = this->errZip(invPt, eta)           * zipErrorScale;
    
  Measurement1D pt(valPt, errPt);
  Measurement1D phi(valPhi, errPhi);
  Measurement1D cotTheta(valCotTheta, errCotTheta);
  Measurement1D tip(valTip, errTip);
  Measurement1D zip(valZip, errZip);

  float chi2 = 0.;

/*
  std::ostringstream str;
    str <<"\t ipt: " << invPt <<"+/-"<<invPtErr
        <<"\t pt:  " << pt.value() <<"+/-"<<pt.error()
        <<"\t phi: " << phi.value() <<"+/-"<<phi.error()
        <<"\t cot: " << cotTheta.value() <<"+/-"<<cotTheta.error()
        <<"\t tip: " << tip.value() <<"+/-"<<tip.error()
        <<"\t zip: " << zip.value() <<"+/-"<<zip.error()
        <<"\t charge: " << charge;
  std::cout <<str.str()<< std::endl;
*/

 
  return  builder.build(pt, phi, cotTheta, tip, zip, chi2, charge, hits, fieldESH.product());
}


double L1MuonPixelTrackFitter::valPhi(const Circle &circle, int charge) const
{
  Circle::Point zero(0.,0.,0.);
  Circle::Vector center = circle.center()-zero;
  long double radius = center.perp();
  center /= radius;
  Circle::Vector dir = center.cross(-charge*Circle::Vector(0.,0.,1.));
  return dir.phi();
}

double L1MuonPixelTrackFitter::errPhi(double invPt, double eta) const {
  //
  // sigma = p0+p1/|pt|;
  //
  double p0,p1;
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p0 = 0.000597506; p1 = 0.00221057; break; }
  case 1:  { p0 = 0.000591867; p1 = 0.00278744; break; }
  case 2:  { p0 = 0.000635666; p1 = 0.00207433; break; }
  case 3:  { p0 = 0.000619086; p1 = 0.00255121; break; }
  case 4:  { p0 = 0.000572067; p1 = 0.00310618; break; }
  case 5:  { p0 = 0.000596239; p1 = 0.00288442; break; }
  case 6:  { p0 = 0.000607608; p1 = 0.00282996; break; }
  case 7:  { p0 = 0.000606446; p1 = 0.00281118; break; }
  case 8:  { p0 = 0.000594076; p1 = 0.00280546; break; }
  case 9:  { p0 = 0.000579615; p1 = 0.00335534; break; }
  case 10:  { p0 = 0.000659546; p1 = 0.00340443; break; }
  case 11:  { p0 = 0.000659031; p1 = 0.00343151; break; }
  case 12:  { p0 = 0.000738391; p1 = 0.00337297; break; }
  case 13:  { p0 = 0.000798966; p1 = 0.00330008; break; }
  case 14:  { p0 = 0.000702997; p1 = 0.00562643; break; }
  case 15:  { p0 = 0.000973417; p1 = 0.00312666; break; }
  case 16:  { p0 = 0.000995213; p1 = 0.00564278; break; }
  case 17:  { p0 = 0.00121436; p1 = 0.00572704; break; }
  case 18:  { p0 = 0.00119216; p1 = 0.00760204; break; }
  case 19:  { p0 = 0.00141204; p1 = 0.0093777; break; }
  default:  { p0 = 0.00153161; p1 = 0.00940265; break; }
  }
  return p0+p1*fabs(invPt);
}


double L1MuonPixelTrackFitter::valCotTheta(const PixelRecoLineRZ& line) const
{
  return line.cotLine();
}

double L1MuonPixelTrackFitter::errCotTheta(double invPt, double eta) const 
{
  //
  // sigma = p0+p1/|pt|;
  //
  double p0,p1;
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p0 = 0.00166115; p1 = 5.75533e-05; break; }
  case 1:  { p0 = 0.00157525; p1 = 0.000707437; break; }
  case 2:  { p0 = 0.00122246; p1 = 0.000325456; break; }
  case 3:  { p0 = 0.000852422; p1 = 0.000429216; break; }
  case 4:  { p0 = 0.000637561; p1 = 0.00122298; break; }
  case 5:  { p0 = 0.000555766; p1 = 0.00158096; break; }
  case 6:  { p0 = 0.000641202; p1 = 0.00143339; break; }
  case 7:  { p0 = 0.000803207; p1 = 0.000648816; break; }
  case 8:  { p0 = 0.000741394; p1 = 0.0015289; break; }
  case 9:  { p0 = 0.000652019; p1 = 0.00168873; break; }
  case 10:  { p0 = 0.000716902; p1 = 0.00257556; break; }
  case 11:  { p0 = 0.000800409; p1 = 0.00190563; break; }
  case 12:  { p0 = 0.000808778; p1 = 0.00264139; break; }
  case 13:  { p0 = 0.000775757; p1 = 0.00318478; break; }
  case 14:  { p0 = 0.000705781; p1 = 0.00460576; break; }
  case 15:  { p0 = 0.000580679; p1 = 0.00748248; break; }
  case 16:  { p0 = 0.000561667; p1 = 0.00767487; break; }
  case 17:  { p0 = 0.000521626; p1 = 0.0100178; break; }
  case 18:  { p0 = 0.00064253; p1 = 0.0106062; break; }
  case 19:  { p0 = 0.000636868; p1 = 0.0140047; break; }
  default:  { p0 = 0.000682478; p1 = 0.0163569; break; }
  }
  return p0+p1*fabs(invPt);
}


double L1MuonPixelTrackFitter::valTip(const Circle &circle, double curvature) const
{
  Circle::Point zero(0.,0.,0.);
  Circle::Vector center = circle.center()-zero;
  long double radius = center.perp();
  return radius-1./fabs(curvature);
}

double L1MuonPixelTrackFitter::errTip(double invPt, double eta) const {
  //
  // sigma = p0+p1/|pt|;
  //
  double p0,p1;
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p0 = 0.00392416; p1 = 0.00551809; break; }
  case 1:  { p0 = 0.00390391; p1 = 0.00543244; break; }
  case 2:  { p0 = 0.0040651; p1 = 0.00406496; break; }
  case 3:  { p0 = 0.00387782; p1 = 0.00797637; break; }
  case 4:  { p0 = 0.00376798; p1 = 0.00866894; break; }
  case 5:  { p0 = 0.0042131; p1 = 0.00462184; break; }
  case 6:  { p0 = 0.00392579; p1 = 0.00784685; break; }
  case 7:  { p0 = 0.00370472; p1 = 0.00790174; break; }
  case 8:  { p0 = 0.00364433; p1 = 0.00928368; break; }
  case 9:  { p0 = 0.00387578; p1 = 0.00640431; break; }
  case 10:  { p0 = 0.00382464; p1 = 0.00960763; break; }
  case 11:  { p0 = 0.0038907; p1 = 0.0104562; break; }
  case 12:  { p0 = 0.00392525; p1 = 0.0106442; break; }
  case 13:  { p0 = 0.00400634; p1 = 0.011218; break; }
  case 14:  { p0 = 0.0036229; p1 = 0.0156403; break; }
  case 15:  { p0 = 0.00444317; p1 = 0.00832987; break; }
  case 16:  { p0 = 0.00465492; p1 = 0.0179908; break; }
  case 17:  { p0 = 0.0049652; p1 = 0.0216647; break; }
  case 18:  { p0 = 0.0051395; p1 = 0.0233692; break; }
  case 19:  { p0 = 0.0062917; p1 = 0.0262175; break; }
  default:  { p0 = 0.00714444; p1 = 0.0253856; break; }
  }
  return p0+p1*fabs(invPt);
}


double L1MuonPixelTrackFitter::valZip(double curv,
    const GlobalPoint& pinner, const GlobalPoint& pouter) const
{
  //
  // phi = asin(r*rho/2) with asin(x) ~= x+x**3/(2*3)
  //
  double rho3 = curv*curv*curv;
  double r1 = pinner.perp();
  double phi1 = r1*curv/2 + pinner.perp2()*r1*rho3/48.;
  double r2 = pouter.perp();
  double phi2 = r2*curv/2 + pouter.perp2()*r2*rho3/48.;
  double z1 = pinner.z();
  double z2 = pouter.z();

  return z1 - phi1/(phi1-phi2)*(z1-z2);
}


double L1MuonPixelTrackFitter::errZip(double invPt, double eta) const {
  //
  // sigma = p0+p1/pt;
  //
  double p0,p1;
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p0 = 0.0120743; p1 = 0; break; }
  case 1:  { p0 = 0.0110343; p1 = 0.0051199; break; }
  case 2:  { p0 = 0.00846487; p1 = 0.00570084; break; }
  case 3:  { p0 = 0.00668726; p1 = 0.00331165; break; }
  case 4:  { p0 = 0.00467126; p1 = 0.00578239; break; }
  case 5:  { p0 = 0.0043042; p1 = 0.00598517; break; }
  case 6:  { p0 = 0.00515392; p1 = 0.00495422; break; }
  case 7:  { p0 = 0.0060843; p1 = 0.00320512; break; }
  case 8:  { p0 = 0.00564942; p1 = 0.00478876; break; }
  case 9:  { p0 = 0.00532111; p1 = 0.0073239; break; }
  case 10:  { p0 = 0.00579429; p1 = 0.00952782; break; }
  case 11:  { p0 = 0.00614229; p1 = 0.00977795; break; }
  case 12:  { p0 = 0.00714661; p1 = 0.00550482; break; }
  case 13:  { p0 = 0.0066593; p1 = 0.00999362; break; }
  case 14:  { p0 = 0.00634922; p1 = 0.0148156; break; }
  case 15:  { p0 = 0.00600586; p1 = 0.0318022; break; }
  case 16:  { p0 = 0.00676919; p1 = 0.027456; break; }
  case 17:  { p0 = 0.00670066; p1 = 0.0317005; break; }
  case 18:  { p0 = 0.00752392; p1 = 0.0347714; break; }
  case 19:  { p0 = 0.00791425; p1 = 0.0566665; break; }
  default:  { p0 = 0.00882372; p1 = 0.0596858; break; }
  }
  return  p0+p1*fabs(invPt);
}
  

double L1MuonPixelTrackFitter::valInversePt( double phi0, double phiL1, double eta) const
{
   double result = 0.;   
   
  // solve equtaion p3*result^3 + p1*result + dphi = 0;
  // where:  result = 1/pt 
  //         dphi = phiL1 - phi0
  // Cardan way is used
  
  double p1,p2,p3;        //parameters p1,p3 are negative by parameter fix in fit!
  param(eta, p1, p2, p3);

  double dphi = deltaPhi(phiL1,phi0);    // phi0-phiL1
  if ( fabs(dphi) < 0.01) {
    result = -dphi/p1; 
  } else {
    double q = dphi/2./p3;
    double p = p1/3./p3;                 // positive
    double D = q*q + p*p*p;              // positive 
    double u = pow(-q+sqrt(D), 1./3.);
    double v = -pow(q+sqrt(D), 1./3.);
    result = u+v;
  }
  return result;
}


double L1MuonPixelTrackFitter::errInversePt(double invPt, double eta) const {
  //
  // pt*sigma(1/pt) = p0+p1*pt;
  //
  double p0,p1;
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p0 = 0.0196835; p1 = 0.00517533; break; }
  case 1:  { p0 = 0.0266583; p1 = 0.00478101; break; }
  case 2:  { p0 = 0.0217164; p1 = 0.00545425; break; }
  case 3:  { p0 = 0.0197547; p1 = 0.00552263; break; }
  case 4:  { p0 = 0.0208778; p1 = 0.00536009; break; }
  case 5:  { p0 = 0.024192; p1 = 0.00521709; break; }
  case 6:  { p0 = 0.0265315; p1 = 0.0051897; break; }
  case 7:  { p0 = 0.0198071; p1 = 0.00566822; break; }
  case 8:  { p0 = 0.0361955; p1 = 0.00486352; break; }
  case 9:  { p0 = 0.037864; p1 = 0.00509094; break; }
  case 10:  { p0 = 0.0382968; p1 = 0.00612354; break; }
  case 11:  { p0 = 0.0308326; p1 = 0.0074234; break; }
  case 12:  { p0 = 0.0248577; p1 = 0.00883049; break; }
  case 13:  { p0 = 0.0279965; p1 = 0.00888293; break; }
  case 14:  { p0 = 0.0372582; p1 = 0.00950252; break; }
  case 15:  { p0 = 0.0281366; p1 = 0.0111501; break; }
  case 16:  { p0 = 0.0421483; p1 = 0.0109413; break; }
  case 17:  { p0 = 0.0461798; p1 = 0.0125824; break; }
  case 18:  { p0 = 0.0530603; p1 = 0.0132638; break; }
  case 19:  { p0 = 0.0601148; p1 = 0.0147911; break; }
  default:  { p0 = 0.0552377; p1 = 0.0155574; break; }
  }
  return p1+p0*fabs(invPt);
}

double L1MuonPixelTrackFitter::findPt(double phi0, double phiL1, double eta, int charge) const {

  double dphi_min = fabs(deltaPhi(phi0,phiL1));
  double pt_best = 1.;
  double pt_cur = 1;
  while ( pt_cur < 10000.) {
    double phi_exp = phi0+getBending(1./pt_cur, eta, charge);
    double dphi = fabs(deltaPhi(phi_exp,phiL1));
    if ( dphi < dphi_min) {
      pt_best = pt_cur;
      dphi_min = dphi;
    }
    if (pt_cur < 10.) pt_cur += 0.01; 
    else if (pt_cur < 20.) pt_cur+= 0.025; 
    else if( pt_cur < 100.) pt_cur+= 0.1;
    else pt_cur +=1;
  };
  return pt_best;
}


double L1MuonPixelTrackFitter::getBending(double invPt, double eta, int charge) 
{
  double p1, p2, p3;
  param(eta,p1,p2,p3);
  return charge*p1*invPt + charge*p2*invPt*invPt + charge*p3*invPt*invPt*invPt;
}

double L1MuonPixelTrackFitter::getBendingError(double invPt, double eta) 
{
  int ieta = int (10*fabs(eta));
  double p0,p1;
  switch (ieta) {
  case 0:  { p0 = 0.0196835; p1 = 0.00517533; break; }
  case 1:  { p0 = 0.0266583; p1 = 0.00478101; break; }
  case 2:  { p0 = 0.0217164; p1 = 0.00545425; break; }
  case 3:  { p0 = 0.0197547; p1 = 0.00552263; break; }
  case 4:  { p0 = 0.0208778; p1 = 0.00536009; break; }
  case 5:  { p0 = 0.024192; p1 = 0.00521709; break; }
  case 6:  { p0 = 0.0265315; p1 = 0.0051897; break; }
  case 7:  { p0 = 0.0198071; p1 = 0.00566822; break; }
  case 8:  { p0 = 0.0361955; p1 = 0.00486352; break; }
  case 9:  { p0 = 0.037864; p1 = 0.00509094; break; }
  case 10:  { p0 = 0.0382968; p1 = 0.00612354; break; }
  case 11:  { p0 = 0.0308326; p1 = 0.0074234; break; }
  case 12:  { p0 = 0.0248577; p1 = 0.00883049; break; }
  case 13:  { p0 = 0.0279965; p1 = 0.00888293; break; }
  case 14:  { p0 = 0.0372582; p1 = 0.00950252; break; }
  case 15:  { p0 = 0.0281366; p1 = 0.0111501; break; }
  case 16:  { p0 = 0.0421483; p1 = 0.0109413; break; }
  case 17:  { p0 = 0.0461798; p1 = 0.0125824; break; }
  case 18:  { p0 = 0.0530603; p1 = 0.0132638; break; }
  case 19:  { p0 = 0.0601148; p1 = 0.0147911; break; }
  default:  { p0 = 0.0552377; p1 = 0.0155574; break; }
  }
  return p0+p1*sqr(invPt);
}


void L1MuonPixelTrackFitter::param(double eta, double &p1, double& p2, double& p3) 
{
  int ieta = int (10*fabs(eta));
  switch (ieta) {
  case 0:  { p1 = -2.68016; p2 = 0; p3 = -12.9653; break; }
  case 1:  { p1 = -2.67864; p2 = 0; p3 = -12.0036; break; }
  case 2:  { p1 = -2.72997; p2 = 0; p3 = -10.3468; break; }
  case 3:  { p1 = -2.68836; p2 = 0; p3 = -12.3369; break; }
  case 4:  { p1 = -2.66885; p2 = 0; p3 = -11.589; break; }
  case 5:  { p1 = -2.64932; p2 = 0; p3 = -12.7176; break; }
  case 6:  { p1 = -2.64513; p2 = 0; p3 = -11.3912; break; }
  case 7:  { p1 = -2.64487; p2 = 0; p3 = -9.83239; break; }
  case 8:  { p1 = -2.58875; p2 = 0; p3 = -10.0541; break; }
  case 9:  { p1 = -2.5551; p2 = 0; p3 = -9.75757; break; }
  case 10:  { p1 = -2.55294; p2 = 0; p3 = -7.25238; break; }
  case 11:  { p1 = -2.45333; p2 = 0; p3 = -7.966; break; }
  case 12:  { p1 = -2.29283; p2 = 0; p3 = -7.03231; break; }
  case 13:  { p1 = -2.13167; p2 = 0; p3 = -4.29182; break; }
  case 14:  { p1 = -1.9102; p2 = 0; p3 = -3.21295; break; }
  case 15:  { p1 = -1.74552; p2 = 0; p3 = -0.531827; break; }
  case 16:  { p1 = -1.53642; p2 = 0; p3 = -1.74057; break; }
  case 17:  { p1 = -1.39446; p2 = 0; p3 = -1.56819; break; }
  case 18:  { p1 = -1.26176; p2 = 0; p3 = -0.598631; break; }
  case 19:  { p1 = -1.14133; p2 = 0; p3 = -0.0941055; break; }
  default:  { p1 = -1.02086; p2 = 0; p3 = -0.100491; break; }
  }
}

double L1MuonPixelTrackFitter::deltaPhi(double phi1, double phi2) const
{
  while ( phi1 >= 2*M_PI) phi1 -= 2*M_PI;
  while ( phi2 >= 2*M_PI) phi2 -= 2*M_PI;
  while ( phi1 < 0) phi1 += 2*M_PI;
  while ( phi2 < 0) phi2 += 2*M_PI;
  double dPhi = phi2-phi1;

  if ( dPhi > M_PI ) dPhi -= 2*M_PI;
  if ( dPhi < -M_PI ) dPhi += 2*M_PI;

  return dPhi;
}