src/smear/PlanarTracker.cxx

0001 /**
0002  \file
0003  Implementation of class Smear::PlanarTracker.
0004
0005  \author    Will Foreman
0006  \date      2011-08-19
0007  \copyright 2011 Brookhaven National Lab
0008  */
0009
0010 #include "eicsmear/smear/PlanarTracker.h"
0011
0012 #include <algorithm>
0013 #include <cmath>
0014 #include <limits>
0015 #include <list>
0016
0017 #include <TMath.h>
0018
0019 namespace Smear {
0020
0021 PlanarTracker::PlanarTracker()
0022 : Tracker(2., 0.03, 0.001)
0023 , mNPlanes(25.)
0024 , mInnerRadius(0.1)
0025 , mOuterRadius(1.)
0026 , mZMin(-1.5)
0027 , mZMax(1.5) {
0028 }
0029
0030 PlanarTracker::PlanarTracker(double inner, double outer,
0031                              double zmin, double zmax,
0032                              double Bb, double nrl,
0033                              double sigmaRPhi, double N)
0034 : Tracker(Bb, nrl, sigmaRPhi)
0035 , mNPlanes(N)
0036 , mInnerRadius(inner)
0037 , mOuterRadius(outer)
0038 // Require zmin < zmax
0039 , mZMin(std::min(zmin, zmax))
0040 , mZMax(std::max(zmin, zmax)) {
0041 }
0042
0043 PlanarTracker::~PlanarTracker() {
0044 }
0045
0046 void PlanarTracker::Print(Option_t* /* option */) const {
0047   std::cout << ClassName() << " with:" << std::endl <<
0048   "\tinner radius " << mInnerRadius << " m\n" <<
0049   "\touter radius " << mOuterRadius << " m\n" <<
0050   "\tlength " << mZMin << " to " << mZMax <<
0051   " (= " << mZMax - mZMin << ") m\n" <<
0052   "\tmagnetic field " << mMagField << " Tesla\n" <<
0053   "\t" << mNRadLengths << " radiation lengths\n" <<
0054   "\tpoint resolution " << mSigmaRPhi * 1.e6 << " microns\n" <<
0055   "\t" << mNPlanes << " planes" << std::endl;
0056 }
0057
0058 double PlanarTracker::GetThetaMin() const {
0059   if (mZMax > 0.) {
0060     return atan2(mInnerRadius, mZMax);
0061   } else {
0062     return atan2(mOuterRadius, mZMax);
0063   }  // if
0064 }
0065
0066 double PlanarTracker::GetThetaMax() const {
0067   if (mZMin > 0.) {
0068     return atan2(mOuterRadius, mZMin);
0069   } else {
0070     return atan2(mInnerRadius, mZMin);
0071   }  // if
0072 }
0073
0074 TVector3 PlanarTracker::ComputeIntersectionWithRadius(
0075       const erhic::VirtualParticle& p, double radius) const {
0076   // Compute the longitudinal position at which the intersection occurs.
0077   // Adjust for any z offset in the vertex of the particle.
0078   const double z = radius / tan(p.GetTheta()) + p.GetVertex().z();
0079   TVector3 intersection(0., 0., std::numeric_limits<double>::quiet_NaN());
0080   if (z > mZMin && z < mZMax) {
0081     // Don't care about phi angle, so set x and y arbitrarily.
0082     intersection.SetXYZ(radius, 0., z);
0083   }  // if
0084   return intersection;
0085 }
0086
0087 TVector3 PlanarTracker::ComputeIntersectionWithPlane(
0088       const erhic::VirtualParticle& p, double z) const {
0089   // Compute the radius of intersection, adusting for any z ooffset
0090   // in the particle vertex.
0091   const double r = (z - p.GetVertex().z()) * tan(p.GetTheta());
0092   TVector3 intersection(0., 0., std::numeric_limits<double>::quiet_NaN());
0093   if (r > mInnerRadius && r < mOuterRadius) {
0094     // Don't care about phi angle, so set x and y arbitrarily.
0095     intersection.SetXYZ(r, 0., z);
0096   }  // if
0097   return intersection;
0098 }
0099
0100 double PlanarTracker::L(const erhic::VirtualParticle& p) const {
0101   float r1 = 0., r2 = 0., pi = 3.1415926535;
0102   double zPlane1 = 0., zPlane2 = 0., Length = 0.;
0103   // Find radial/z points of first and last plane intersection.
0104   // (While there are certainly eaiser ways to carry out this
0105   // calculation, I used a manual check at each z-plane since
0106   // it is easier to generalize to arbitrarily spaced planes.)
0107   for (int i = 0; i < mNPlanes; i++) {
0108     zPlane1 = mZMin + i * (mZMax - mZMin) / (mNPlanes - 1);
0109     r1 = fabs(tan(p.GetTheta()) * zPlane1);
0110     if ((r1 > mInnerRadius) && (r1 < mOuterRadius)) {
0111       if (((p.GetTheta() < pi / 2.) && (zPlane1 > 0.)) ||
0112          ((p.GetTheta() > pi / 2.) && (zPlane1 < 0.))) {
0113         zPlane2 = zPlane1 + (NPoints(p) - 1) * (mZMax - mZMin) / (mNPlanes - 1);
0114         r2 = fabs(tan(p.GetTheta()) * zPlane2);
0115         break;
0116       }  // if
0117     }  // if
0118   }  // for
0119   Length = sqrt((zPlane2 - zPlane1) * (zPlane2 - zPlane1) +
0120                 (r2 - r1) * (r2 - r1));
0121   return Length;
0122 }
0123
0124 double PlanarTracker::LPrime(const erhic::VirtualParticle& p) const {
0125   float r1 = 0., r2 = 0., pi = 3.1415926535;
0126   double zPlane1 = 0., zPlane2 = 0.;
0127   // Find radial points of first/ and last plane intersection.
0128   for (int i = 0; i < mNPlanes; i++) {
0129     zPlane1 = mZMin + i * (mZMax - mZMin) / (mNPlanes-1);
0130     r1 = fabs(tan(p.GetTheta()) * zPlane1);
0131     if ((r1 > mInnerRadius) && (r1 < mOuterRadius)) {
0132       if (((p.GetTheta() < pi / 2.) && (zPlane1 > 0.)) ||
0133          ((p.GetTheta() > pi / 2.) && (zPlane1 < 0.))) {
0134         zPlane2 = zPlane1 + (NPoints(p) - 1) * (mZMax - mZMin) / (mNPlanes - 1);
0135         r2 = fabs(tan(p.GetTheta()) * zPlane2);
0136         break;
0137       }  // if
0138     }  // if
0139   }  // for
0140   return fabs(r2 - r1);
0141 }
0142
0143 int PlanarTracker::NPoints(const erhic::VirtualParticle& p) const {
0144   double zPlane = 0., pi = 3.1415926535;
0145   float r = 0.;
0146   int i, n = 0;
0147   for (i = 0; i < mNPlanes; i++) {
0148     zPlane = mZMin + i * (mZMax - mZMin) / (mNPlanes - 1);
0149     // radial intersect of particle at plane position zPlane
0150     r = fabs(tan(p.GetTheta()) * zPlane);
0151     // At each plane, check if particle passes through
0152     if ((r > mInnerRadius) && (r < mOuterRadius)) {
0153       if ((p.GetTheta() < pi / 2.) && (zPlane > 0)) n++;
0154       if ((p.GetTheta() > pi / 2.) && (zPlane < 0)) n++;
0155     }  // if
0156   }  // for
0157   return n;
0158 }
0159
0160 bool PlanarTracker::Accepts(const erhic::VirtualParticle& p) const {
0161   if (NPoints(p) >= 3) {
0162     return true;
0163   } else {
0164     return false;
0165   }  // if
0166 }
0167
0168 }  // namespace Smear