HitInHyperBox.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 <tracking/trackFindingCDC/hough/boxes/Box.h>
#include <tracking/trackFindingCDC/hough/axes/StandardAxes.h>
#include <tracking/trackingUtilities/numerics/Weight.h>
#include <tracking/trackingUtilities/eventdata/hits/CDCRecoHit3D.h>
#include <tracking/trackFindingCDC/hough/baseelements/SameSignChecker.h>
 
#include <framework/gearbox/Unit.h>
 
#include <cmath>
 
namespace Belle2 {
  namespace TrackFindingCDC {

    class HitInHyperBox {
 
    public:
      using HoughBox = Box<DiscreteQ, DiscreteP, DiscreteZ0>;

      TrackingUtilities::Weight operator()(const TrackingUtilities::CDCRecoHit3D& recoHit,
                                           const HoughBox* hyperBox)
      {
        const float lowerQ = *(hyperBox->getLowerBound<DiscreteQ>()); //DiscreteValue is based on std::vector<T>::const_iterator
        const float upperQ = *(hyperBox->getUpperBound<DiscreteQ>());
        const float deltaQ = 0.5 * (upperQ - lowerQ);
 
        const float centerP = *(hyperBox->getCenter<1>()); //TODO getCenter(class T) is not implemented
 
        const float centerZ0 = *(hyperBox->getCenter<2>());
 
        const double perpS = recoHit.getArcLength2D();
        const double recoZ = recoHit.getRecoZ();
 
        const bool sameSign = SameSignChecker::sameSign(catZ(lowerQ, centerP, perpS) + centerZ0 - 100.0 * Unit::cm * deltaQ - recoZ,
                                                        catZ(upperQ, centerP, perpS) + centerZ0 + 100.0 * Unit::cm * deltaQ - recoZ); //100 is a reference - size of CDC in cm
        if (not sameSign) {
          return 1.0;
        } else {
          return NAN;
        }
      }

      static bool compareDistances(const HoughBox& hyperBox, const TrackingUtilities::CDCRecoHit3D& lhsRecoHit,
                                   const TrackingUtilities::CDCRecoHit3D& rhsRecoHit)
      {
        const float centerQ = *(hyperBox.getCenter<0>());
        const float centerP = *(hyperBox.getCenter<1>());
        const float centerZ0 = *(hyperBox.getCenter<2>());
 
        const double lhsZ = lhsRecoHit.getRecoZ();
        const double rhsZ = rhsRecoHit.getRecoZ();
 
        const double lhsS = lhsRecoHit.getArcLength2D();
        const double rhsS = rhsRecoHit.getArcLength2D();
 
        const double lhsZDistance = catZ(centerQ, centerP, lhsS) + centerZ0 - lhsZ;
        const double rhsZDistance = catZ(centerQ, centerP, rhsS) + centerZ0 - rhsZ;
 
        return lhsZDistance < rhsZDistance;
      }
 
    private:
      static double catZ(const double q, const double p, const double R)
      {
        //100 here is a reference  - size of CDC in cm
        return 100.0 * Unit::cm * q * (std::sqrt(1 - p * p) * std::cosh(R / (100.0 * Unit::cm) + std::asinh(p / std::sqrt(1 - p * p))) - 1);
      }
 
    public:
      static float centerX(const HoughBox& hyperBox)
      {
        return *(hyperBox.getCenter<0>());
      }

      static float centerY(const HoughBox& hyperBox)
      {
        return *(hyperBox.getCenter<1>());
      }

      static float centerZ(const HoughBox& hyperBox)
      {
        return *(hyperBox.getCenter<2>());
      }

      static float deltaX(const HoughBox& hyperBox)
      {
        const float lowerQ = *(hyperBox.getLowerBound<DiscreteQ>()); //DiscreteValue is based on std::vector<T>::const_iterator
        const float upperQ = *(hyperBox.getUpperBound<DiscreteQ>());
        return 0.5 * (upperQ - lowerQ);
      }

      static const char* debugLine() { return "100.0 * [0] * (TMath::Sqrt(1 - [1] * [1]) * TMath::CosH(x / 100.0 + TMath::ASinH([1] / TMath::Sqrt(1 - [1] * [1]))) - 1) + [2]";}
    };
  }
}