Tracks/src/estimate_field.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:     Tracks
0004 // Class  :     estimate_field
0005 //
0006
0007 #include "Math/Vector3D.h"
0008 #include "DataFormats/TrackReco/interface/Track.h"
0009 #include "Fireworks/Tracks/interface/estimate_field.h"
0010
0011 double fw::estimate_field(const reco::Track& track, bool highQuality) {
0012   if (!track.extra().isAvailable())
0013     return -1;
0014   // We have 3 states so we can form 3 pairs to check for curvature
0015   // In practice we want to avoid estimates over lost of material and
0016   // in regions with small field. Therefore we should look at the pair
0017   // based on POCA and closest state to it. Than we will estimate the
0018   // filed inside the solenoid.
0019   double estimate = -1;
0020   if (pow(track.outerPosition().x() - track.vx(), 2) + pow(track.outerPosition().y() - track.vy(), 2) <
0021       pow(track.innerPosition().x() - track.vx(), 2) + pow(track.innerPosition().y() - track.vy(), 2)) {
0022     estimate = estimate_field(
0023         track.outerPosition().x(), track.outerPosition().y(), track.vx(), track.vy(), track.px(), track.py());
0024   } else {
0025     estimate = estimate_field(
0026         track.innerPosition().x(), track.innerPosition().y(), track.vx(), track.vy(), track.px(), track.py());
0027   }
0028
0029   if (highQuality) {
0030     if (fabs(track.outerMomentum().rho() / track.pt() - 1) < 0.1 && track.innerPosition().rho() < 129 &&
0031         track.outerPosition().rho() < 129) {
0032       /*printf("outer-poca: \tPOCA Pt: %0.2f, \tinner Rho: %0.1f, \touter Rho: %0.1f, \tchange in pt: %0.1f%%, \testimate: %0.2f\n",
0033           track.pt(), track.innerPosition().rho(), track.outerPosition().rho(),
0034           fabs( track.outerMomentum().rho()/track.pt()-1 )*100, estimate);*/
0035       return estimate;
0036     } else {
0037       return -1;
0038     }
0039   } else
0040     return estimate;
0041 }
0042
0043 double fw::estimate_field(double vx1, double vy1, double vx2, double vy2, double px, double py) {
0044   math::XYZVector displacement(vx2 - vx1, vy2 - vy1, 0);
0045   math::XYZVector transverseMomentum(px, py, 0);
0046   double cosAlpha = transverseMomentum.Dot(displacement) / transverseMomentum.r() / displacement.r();
0047   return 200 * sqrt(1 - cosAlpha * cosAlpha) / 0.2998 * transverseMomentum.r() / displacement.r();
0048 }