estimate_field.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     Tracks
// Class  :     estimate_field
//

#include "Math/Vector3D.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "Fireworks/Tracks/interface/estimate_field.h"

double
fw::estimate_field( const reco::Track& track, bool highQuality )
{
   if ( !track.extra().isAvailable() ) return -1;
   // We have 3 states so we can form 3 pairs to check for curvature
   // In practice we want to avoid estimates over lost of material and
   // in regions with small field. Therefore we should look at the pair
   // based on POCA and closest state to it. Than we will estimate the
   // filed inside the solenoid. 
   double estimate = -1;
   if ( pow(track.outerPosition().x()-track.vx(),2)+pow(track.outerPosition().y()-track.vy(),2) <
    pow(track.innerPosition().x()-track.vx(),2)+pow(track.innerPosition().y()-track.vy(),2) )
     {
       estimate = estimate_field( track.outerPosition().x(), 
                  track.outerPosition().y(),
                  track.vx(),
                  track.vy(),
                  track.px(),
                  track.py() );
     }
   else
     {
       estimate = estimate_field( track.innerPosition().x(), 
                  track.innerPosition().y(),
                  track.vx(),
                  track.vy(),
                  track.px(),
                  track.py() );
     }
   
   if ( highQuality ){
     if ( fabs( track.outerMomentum().rho()/track.pt()-1 )<0.1 && 
      track.innerPosition().rho()<129 &&
      track.outerPosition().rho()<129 ){
       /*printf("outer-poca: \tPOCA Pt: %0.2f, \tinner Rho: %0.1f, \touter Rho: %0.1f, \tchange in pt: %0.1f%%, \testimate: %0.2f\n",
          track.pt(), track.innerPosition().rho(), track.outerPosition().rho(), 
          fabs( track.outerMomentum().rho()/track.pt()-1 )*100, estimate);*/
       return estimate;
     } else {
       return -1;
     }
   } else
     return estimate;
}

double
fw::estimate_field( double vx1, double vy1, double vx2, double vy2, double px, double py )
{
  math::XYZVector displacement(vx2-vx1, vy2-vy1, 0);
  math::XYZVector transverseMomentum(px, py, 0);
  double cosAlpha = transverseMomentum.Dot(displacement)/transverseMomentum.r()/displacement.r();
  return 200*sqrt(1-cosAlpha*cosAlpha)/0.2998*transverseMomentum.r()/displacement.r();
}