ThreeHitFilters.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <TVector3.h>
#include "tracking/vxdCaTracking/TwoHitFilters.h"
#include "tracking/vxdCaTracking/FilterExceptions.h"
#include <framework/logging/Logger.h>
 
namespace Belle2 {
  class ThreeHitFilters {
  public:
 
    ThreeHitFilters():
      m_circleCenterCalculated(false),
      m_radiusCalculated(false),
      m_radius(0.),
      m_x2(0.),
      m_y2(0.),
      m_z2(0.)
    {
      m_hitA.SetXYZ(0., 0., 0.);
      m_hitB.SetXYZ(0., 0., 0.);
      m_hitC.SetXYZ(0., 0., 0.);
      m_vecAB.SetXYZ(0., 0., 0.);
      m_vecBC.SetXYZ(0., 0., 0.);
      resetMagneticField(1.5);
    }
 
    ThreeHitFilters(const TVector3& outerHit, const TVector3& centerHit, const TVector3& innerHit,
                    const double magneticFieldStrength = 1.5):
      m_circleCenterCalculated(false),
      m_radiusCalculated(false),
      m_radius(0.),
      m_hitA(outerHit),
      m_hitB(centerHit),
      m_hitC(innerHit),
      m_vecAB(outerHit - centerHit),
      m_vecBC(centerHit - innerHit)
    {
      m_x2 = m_vecAB[0] * m_vecBC[0]; 
      m_y2 = m_vecAB[1] * m_vecBC[1]; 
      m_z2 = m_vecAB[2] * m_vecBC[2]; 
      resetMagneticField(magneticFieldStrength);
    }
 
 
    ~ThreeHitFilters() {}
 
    void resetValues(const TVector3& outerHit, const TVector3& centerHit, const TVector3& innerHit)
    {
      m_radiusCalculated = false;
      m_circleCenterCalculated = false;
      m_radius = 0.;
      m_hitA = outerHit;
      m_hitB = centerHit;
      m_hitC = innerHit;
      m_vecAB = outerHit - centerHit;
      m_vecBC = centerHit - innerHit;
 
      m_x2 = m_vecAB[0] * m_vecBC[0]; 
      m_y2 = m_vecAB[1] * m_vecBC[1]; 
      m_z2 = m_vecAB[2] * m_vecBC[2]; 
    }
 
 
 
    void resetMagneticField(const double magneticFieldStrength = 1.5) { m_magneticFieldFactor = magneticFieldStrength * 0.00299710; }
 
 
 
    double getMagneticField() { return m_magneticFieldFactor / 0.00299710; }
 
 
 
    double filterNan(double value) { return m_twoHitFilter.filterNan(value); }
 
 
 
    double calcAngle3D()
    {
      double angle = ((m_x2 + m_y2 + m_z2) / (m_vecAB.Mag2() *
                                              m_vecBC.Mag2())); // fullCalc would be acos(m_vecAB.Dot(m_vecBC) / m_vecAB.Mag()*m_vecBC.Mag())
      return filterNan(angle);
    } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngle3D instead)
 
 
 
    double fullAngle3D()
    {
      double angle = acos(m_vecAB.Dot(m_vecBC) / (m_vecAB.Mag() * m_vecBC.Mag())); // 0-pi
      angle = (angle * (180. / M_PI));
      return filterNan(angle);
    } // return unit: ° (0 - 180°)
 
 
 
    double calcAngleXY()
    {
      double angle = ((m_x2 + m_y2) / (m_vecAB.Perp2() * m_vecBC.Perp2())); // fullAngle:
      return filterNan(angle);
    } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngleXY instead)
 
 
 
    double fullAngleXY()
    {
      double angle = fullAngle2D(m_vecAB, m_vecBC); // 0-pi
      angle = (angle * (180. / M_PI));
      return filterNan(angle);
    } // return unit: ° (0 - 180°)
 
 
 
    double calcAngleRZ()
    {
      TVector3 rzVecAB(m_vecAB.Perp(), m_vecAB[2], 0.);
      TVector3 rzVecBC(m_vecBC.Perp(), m_vecBC[2], 0.);
      return calcAngle2D(rzVecAB, rzVecBC);
    } // return unit: none (calculation for degrees is incomplete, if you want readable numbers, use fullAngleRZ instead)
 
 
 
    double fullAngleRZ()
    {
      TVector3 rzVecAB(m_vecAB.Perp(), m_vecAB[2], 0.);
      TVector3 rzVecBC(m_vecBC.Perp(), m_vecBC[2], 0.);
      double angle = fullAngle2D(rzVecAB, rzVecBC); // 0-pi
      angle = (angle * (180. / M_PI));
      return filterNan(angle);
    } // return unit: ° (0 - 180°)
 
 
 
    double calcCircleDist2IP()
    {
      checkCalcRadius();
 
      return (fabs(m_centerABC.Perp() - m_radius)); // distance of closest approach of circle to the IP
    } // return unit: cm
 
 
 
    double calcPt()
    {
      checkCalcRadius();
 
      return calcPt(m_radius);
    } // return unit: GeV/c
 
 
 
    double calcPt(double radius)
    {
      sanityCheckRadius(radius);
 
      return (m_magneticFieldFactor * radius);
    } // return unit: GeV/c
 
 
 
    double calcDeltaSlopeRZ()
    {
      m_twoHitFilter.resetValues(m_hitA, m_hitB);
      double slopeAB = m_twoHitFilter.calcSlopeRZ();
      m_twoHitFilter.resetValues(m_hitB, m_hitC);
      double slopeBC = m_twoHitFilter.calcSlopeRZ();
 
      return filterNan(slopeBC - slopeAB);
    } // return unit: none
 
 
 
    double calcDeltaSOverZ()
    {
      checkCalcCircleCenter();
 
      TVector3 points2hitA = m_hitA - m_centerABC;
      TVector3 points2hitB = m_hitB - m_centerABC;
      TVector3 points2hitC = m_hitC - m_centerABC;
      double alfaAB = fullAngle2D(points2hitA, points2hitB);
      double alfaBC = fullAngle2D(points2hitB, points2hitC);
      //return filterNan( (alfaAB * m_vecBC[2]) - (alfaBC *m_vecAB[2]) );
      return (alfaAB * m_vecBC[2]) - (alfaBC *
                                      m_vecAB[2]); // equals to alfaAB/dZAB and alfaBC/dZBC, but this solution here can not produce a division by zero
    } // return unit: radians*cm
 
 
 
    double calcDeltaSlopeZOverS()
    {
      checkCalcRadius();
 
      TVector3 points2hitA = m_hitA - m_centerABC;
      TVector3 points2hitB = m_hitB - m_centerABC;
      TVector3 points2hitC = m_hitC - m_centerABC;
      double alfaABr = fullAngle2D(points2hitA, points2hitB) * m_radius;
      double alfaBCr = fullAngle2D(points2hitB, points2hitC) * m_radius;
 
      return filterNan((asin(m_vecAB[2] / alfaABr)) - asin(m_vecBC[2] / alfaBCr)); // Beware of z>r!
    }
 
 
 
    double calcHelixParameterFit()
    {
      checkCalcCircleCenter();
      TVector3 points2hitA = m_hitA - m_centerABC;
      TVector3 points2hitB = m_hitB - m_centerABC;
      TVector3 points2hitC = m_hitC - m_centerABC;
      double alfaAB = calcAngle2D(points2hitA, points2hitB);
      double alfaBC = calcAngle2D(points2hitB, points2hitC);
      // real calculation: ratio is (m_vecij[2] = deltaZ): alfaAB/deltaZab : alfaBC/deltaZbc, the following equation saves two times '/'
      return filterNan((alfaAB * m_vecBC[2]) / (alfaBC * m_vecAB[2]));
    } // return unit: none
 
 
 
    double calcHelixFit() { return calcHelixParameterFit(); }
 
 
 
    double calcAngle2D(const TVector3& vecA, const TVector3& vecB)
    {
      double angle = ((vecA[0] * vecB[0] + vecA[1] * vecB[1]) / sqrt(vecA.Perp2() * vecB.Perp2()));
      return filterNan(angle);
    }
 
    double fullAngle2D(const TVector3& vecA, const TVector3& vecB)
    {
      return acos(calcAngle2D(vecA, vecB));
      //return filterNan(angle);
    }
 
 
 
    double calcRadius(const TVector3& a, const TVector3& b, const TVector3& c, const TVector3& circleCenter)
    {
      return ((circleCenter - a).Perp() + (circleCenter - b).Perp() + (circleCenter - c).Perp()) /
             3.;   // = radius in [cm], sign here not needed. normally: signKappaAB/normAB1
    } // used by calcPt() and calcCircleDist2IP()
 
 
 
    void calcCircleCenter(const TVector3& a, const TVector3& b, const TVector3& c, TVector3& circleCenter)
    {
      // calculates the intersection point using Cramer's rule.
      // x_1+s*n_1==x_2+t*n_2 --> n_1 *s - n_2 *t == x_2 - x_1 --> http://en.wikipedia.org/wiki/Cramer%27s_rule
      double inX = b[0] - c[0]; // x value of the normal vector of the inner segment (b-c)
      double inY = b[1] - c[1]; // y value of the normal vector of the inner segment (b-c)
      double outX = a[0] - b[0]; // x value of the normal vector of the outer segment (a-b)
      double outY = a[1] - b[1]; // y value of the normal vector of the outer segment (a-b)
 
      //searching solution for Ax = b, aij are the matrix elements of A, bi are elements of b
      double a11 = inY;
      double a12 = -inX;
      double a21 = -outY;
      double a22 = outX;
      double b1 = b[0] + outX * 0.5 - (c[0] + inX * 0.5);
      double b2 = b[1] + outY * 0.5 - (c[1] + inY * 0.5);
 
      if (a11 * a22 == a12 * a21) { throw FilterExceptions::Straight_Line(); }
 
      double s = (b1 * a22 - b2 * a21) / (a11 * a22 - a12 * a21); //the determinant is zero iff the three hits are on a line in (x,y).
 
      circleCenter.SetXYZ(c[0] + inX * 0.5 + s * inY, c[1] + inY * 0.5 - s * inX, 0.);
    }
 
 
 
    int calcSign(const TVector3& a, const TVector3& b, const TVector3& c);
 
 
 
    int calcSign(const TVector3& a, const TVector3& b, const TVector3& c, const TVector3& sigma_a, const TVector3& sigma_b,
                 const TVector3& sigma_c)
    {
      TVector3 c2b = b - c;   c2b.SetZ(0.);
      TVector3 b2a = a - b;   b2a.SetZ(0.);
      double angle = atan2(b2a[0], b2a[1]) - atan2(c2b[0], c2b[1]);
      double sigmaan = (sigma_a.Mag() + sigma_b.Mag() + sigma_c.Mag()) / (3.*(c2b.Mag() +
                       b2a.Mag())); //TODO 1/3...mean of the sigmas. Possible improvement: Use a parameter instead, and determine with simulated events.
      if (angle < (-sigmaan)) { return -1; }
      else if (angle > sigmaan)  {return 1; }
      else  { return 0; }
    }
 
 
 
  protected:
 
 
    void checkCalcCircleCenter()
    {
      if (m_circleCenterCalculated == false) {
        calcCircleCenter(m_hitA, m_hitB, m_hitC, m_centerABC);
        m_circleCenterCalculated = true;
      }
    }
 
 
    void checkCalcRadius()
    {
      checkCalcCircleCenter();
      if (m_radiusCalculated == false) {
        m_radius = calcRadius(m_hitA, m_hitB, m_hitC, m_centerABC);
        sanityCheckRadius(m_radius);
        m_radiusCalculated = true;
      }
    }
 
 
    void sanityCheckRadius(double radius)
    {
      if (fabs(radius) <
          0.0000001) { // WARNING hardcoded value, is there a quasi-global value for such cases (this case, minimal accepted radius)
        m_radiusCalculated = false;
        m_circleCenterCalculated = false;
        m_radius = 0.;
        B2ERROR("sanityCheckRadius: given radius is " << radius << ", which is below threshold of " << 0.0000001 << ", throwing exception");
        throw FilterExceptions::Circle_too_small();
      }
    }
 
 
 
    TwoHitFilters m_twoHitFilter; 
    bool m_circleCenterCalculated; 
    bool m_radiusCalculated; 
    double m_radius; 
    double m_x2; 
    double m_y2; 
    double m_z2; 
    double m_magneticFieldFactor; 
    TVector3 m_centerABC;  
    TVector3 m_hitA; 
    TVector3 m_hitB; 
    TVector3 m_hitC; 
    TVector3 m_vecAB; 
    TVector3 m_vecBC; 
  }; //end class ThreeHitFilters
} //end namespace Belle2