Conv4HitsReco2.cc

Go to the documentation of this file.

//
//  Simple photon conversion seeding class (src)
//
//  Author: E Song
//
//  Version: 1;     6 Aug 2012
//

#define Conv4HitsReco2_cxx

//#include "Conv4HitsReco2.h"
//#include "FWCore/MessegeLogger/interface/MessegeLogger.h"
#include "Conv4HitsReco2.h"
#include <ctime>

Conv4HitsReco2::Conv4HitsReco2(
    math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4) {
  Refresh(vPhotVertex, h1, h2, h3, h4);
}

Conv4HitsReco2::~Conv4HitsReco2() {}
//Conv4HitsReco2::~Conv4HitsReco() { }

void Conv4HitsReco2::Refresh(
    math::XYZVector &vPhotVertex, math::XYZVector &h1, math::XYZVector &h2, math::XYZVector &h3, math::XYZVector &h4) {
  // Fix 2D plane, make primary vertex the original point
  fPV = vPhotVertex;
  fPV.SetZ(0.);
  fHitv11 = h3;
  fHitv11.SetZ(0.);
  fHitv11 = fHitv11 - fPV;
  fHitv21 = h4;
  fHitv21.SetZ(0.);
  fHitv21 = fHitv21 - fPV;
  fHitv12 = h2;
  fHitv12.SetZ(0.);
  fHitv12 = fHitv12 - fPV;
  fHitv22 = h1;
  fHitv22.SetZ(0.);
  fHitv22 = fHitv22 - fPV;

  // DEFAULT setup
  fMaxNumberOfIterations = 40;
  fFixedNumberOfIterations = 0;
  fRadiusECut = 10.0;  //cm
  fPhiECut = 0.03;     //rad
  fRECut = 0.5;        //cm
  fBField = 3.8;       //T

  // TRIVIAL initialization
  fCutSatisfied = 0;
  fSignSatisfied = 0;
  fSolved = 0;

  fRecPhi = 0.;
  fRecR = 0.;
  fRecR1 = 0.;
  fRecR2 = 0.;
  fRadiusE = 0.;
  fRE = 0.;
  fPhiE = 0.;
}

void Conv4HitsReco2::LocalTransformation(const math::XYZVector &v11,
                                         const math::XYZVector &v12,
                                         const math::XYZVector &v21,
                                         const math::XYZVector &v22,
                                         math::XYZVector &V11,
                                         math::XYZVector &V12,
                                         math::XYZVector &V21,
                                         math::XYZVector &V22,
                                         double NextPhi) {
  double x11, x12, x21, x22, y11, y12, y21, y22;

  x11 = v11.X();
  y11 = v11.Y();
  x12 = v12.X();
  y12 = v12.Y();
  x21 = v21.X();
  y21 = v21.Y();
  x22 = v22.X();
  y22 = v22.Y();

  double SINP = std::sin(NextPhi);
  double COSP = std::cos(NextPhi);
  double SignCOSP = 1.;
  if (COSP < 0.)
    SignCOSP = -1.;
  double AbsCOSP = std::fabs(COSP);

  double X11 = -std::fabs(x11 * SINP * SignCOSP - y11 * AbsCOSP);
  double Y11 = std::fabs(y11 * SINP * SignCOSP + x11 * AbsCOSP);

  double X21 = -std::fabs(x21 * SINP * SignCOSP - y21 * AbsCOSP);
  double Y21 = std::fabs(y21 * SINP * SignCOSP + x21 * AbsCOSP);

  double X12 = std::fabs(x12 * SINP * SignCOSP - y12 * AbsCOSP);
  double Y12 = std::fabs(y12 * SINP * SignCOSP + x12 * AbsCOSP);

  double X22 = std::fabs(x22 * SINP * SignCOSP - y22 * AbsCOSP);
  double Y22 = std::fabs(y22 * SINP * SignCOSP + x22 * AbsCOSP);

  V11.SetXYZ(X11, Y11, 0.);
  V12.SetXYZ(X12, Y12, 0.);
  V21.SetXYZ(X21, Y21, 0.);
  V22.SetXYZ(X22, Y22, 0.);
}

int Conv4HitsReco2::ConversionCandidate(math::XYZVector &vtx, double &ptplus, double &ptminus) {
  Reconstruct();

  ptplus = fRecR2 * fBField * 0.01 * 0.3;   // 2 - positron
  ptminus = fRecR1 * fBField * 0.01 * 0.3;  // 1 - electron
  vtx = fRecV;
  //std::cout << ".";
  return fLoop;
}

void Conv4HitsReco2::Reconstruct() {
  double x11, x12, x21, x22, y11, y12, y21, y22;
  //double X11,X12,X21,X22,Y11,Y12,Y21,Y22;
  x11 = fHitv11.X();
  y11 = fHitv11.Y();
  x12 = fHitv12.X();
  y12 = fHitv12.Y();
  x21 = fHitv21.X();
  y21 = fHitv21.Y();
  x22 = fHitv22.X();
  y22 = fHitv22.Y();

  if (fFixedNumberOfIterations == 0)
    fLoop = fMaxNumberOfIterations;
  else
    fLoop = fFixedNumberOfIterations;

  // Setting Phi1, Phi2 initial guess range, and first guess
  double tempr1 = std::sqrt(y11 * y11 + x11 * x11);
  double tempr2 = std::sqrt(y12 * y12 + x12 * x12);

  double Phi1 = 2.0 * std::atan(y11 / (x11 + tempr1));
  double Phi2 = 2.0 * std::atan(y12 / (x12 + tempr2));

  if (Phi1 < Phi2)
    Phi1 += 2.0 * 3.141592653;  // stupid Atan correction

  fPhiE = std::fabs((Phi1 - Phi2)) / std::pow(2.0, fLoop + 1);

  double NextPhi = (Phi1 + Phi2) / 2.0;  // first guess
  double D1, D2 = 0.0;
  double prevR1 = 0;
  double prevR2 = 0;
  double R1 = 0;
  double R2 = 0;

  // Iterations
  for (int i = 0; i < fLoop; i++) {
    // LOCAL TRANFORMATION & EXTRACTION
    double SINP = std::sin(NextPhi);
    double COSP = std::cos(NextPhi);
    double SignCOSP = 1.;
    if (COSP < 0.)
      SignCOSP = -1.;
    double AbsCOSP = std::fabs(COSP);

    double X11 = -std::fabs(x11 * SINP * SignCOSP - y11 * AbsCOSP);
    double Y11 = std::fabs(y11 * SINP * SignCOSP + x11 * AbsCOSP);

    double X21 = -std::fabs(x21 * SINP * SignCOSP - y21 * AbsCOSP);
    double Y21 = std::fabs(y21 * SINP * SignCOSP + x21 * AbsCOSP);

    double X12 = std::fabs(x12 * SINP * SignCOSP - y12 * AbsCOSP);
    double Y12 = std::fabs(y12 * SINP * SignCOSP + x12 * AbsCOSP);

    double X22 = std::fabs(x22 * SINP * SignCOSP - y22 * AbsCOSP);
    double Y22 = std::fabs(y22 * SINP * SignCOSP + x22 * AbsCOSP);
    // I'm not using LocalTransform() function because this direct way turns out to be faster

    // SOLVING EQUATIONS
    double d1 = Y21 - Y11;
    double d2 = Y22 - Y12;

    if (((X11 * X11 * d1 * d1 / (X21 - X11) / (X21 - X11) + X11 * X21 + X11 * d1 * d1 / (X21 - X11)) < 0) ||
        ((X12 * X12 * d2 * d2 / (X22 - X12) / (X22 - X12) + X12 * X22 + X12 * d2 * d2 / (X22 - X12)) < 0)) {
      fSolved = -1;
      fLoop = i;
      return;
    }  // No real root.  Break out.

    else {
      fSolved = 1;
      D1 = X11 * d1 / (X21 - X11);
      D1 = D1 + std::sqrt(X11 * X11 * d1 * d1 / (X21 - X11) / (X21 - X11) + X11 * X21 + X11 * d1 * d1 / (X21 - X11));
      D2 = X12 * d2 / (X22 - X12);
      D2 = D2 + std::sqrt(X12 * X12 * d2 * d2 / (X22 - X12) / (X22 - X12) + X12 * X22 + X12 * d2 * d2 / (X22 - X12));

      R1 = std::fabs((X11 + X21) / 2.0 + (D1 + d1 / 2.0) * d1 / (X21 - X11));
      R2 = std::fabs((X12 + X22) / 2.0 + (D2 + d2 / 2.0) * d2 / (X22 - X12));

      if ((Y11 - D1) >= (Y12 - D2)) {  // Moving RIGHT
        Phi1 = NextPhi;
        // Phi2 remains the same
        NextPhi = (Phi1 + Phi2) / 2.0;
      } else if ((Y11 - D1) < (Y12 - D2)) {  // Moving LEFT
        // Phi1 remains the same
        Phi2 = NextPhi;
        NextPhi = (Phi1 + Phi2) / 2.0;
      }

      // CHECK STOP CONDITION
      double tmpPhiE = std::fabs(Phi1 - Phi2);
      double tmpRE = std::fabs((Y11 - D1) - (Y12 - D2));
      double tmpRadiusE = (std::fabs(R1 - prevR1) + std::fabs(R2 - prevR2)) / 2.;

      // A. Cut threshold satisfied - STOP - record
      if ((tmpPhiE <= fPhiECut) && (tmpRE <= fRECut) && (tmpRadiusE <= fRadiusECut) &&
          (fFixedNumberOfIterations == 0)) {
        fSolved = 1;
        fCutSatisfied = 1;
        fLoop = i + 1;
        fPhiE = tmpPhiE;
        fRE = tmpRE;
        fRadiusE = tmpRadiusE;
        fRecR1 = R1;
        fRecR2 = R2;
        fRecR = ((Y11 - D1) + (Y12 - D2)) / 2.0;
        fRecPhi = NextPhi;
        fRecV.SetX(fRecR * cos(fRecPhi));
        fRecV.SetY(fRecR * sin(fRecPhi));
        fRecC1.SetXYZ(fRecV.X() - fRecR1 * sin(fRecPhi), fRecV.Y() + fRecR1 * cos(fRecPhi), 0.);
        fRecC2.SetXYZ(fRecV.X() + fRecR2 * sin(fRecPhi), fRecV.Y() - fRecR2 * cos(fRecPhi), 0.);
        fRecV = fRecV + fPV;
        fRecC1 = fRecC1 + fPV;
        fRecC2 = fRecC2 + fPV;
        fCutSatisfied = 1;
        if ((R1 > 0) && (R2 > 0) && (D1 > 0) && (D2 > 0) && ((Y11 - D1) > 0) && ((Y12 - D2) > 0))
          fSignSatisfied = 1;
        else
          fSignSatisfied = 0;
      } else if (i == fLoop - 1) {
        fSolved = 1;
        fCutSatisfied = 1;
        fLoop = i + 1;
        fPhiE = tmpPhiE;
        fRE = tmpRE;
        fRadiusE = tmpRadiusE;
        fRecR1 = R1;
        fRecR2 = R2;
        fRecR = ((Y11 - D1) + (Y12 - D2)) / 2.0;
        fRecPhi = NextPhi;
        fRecV.SetX(fRecR * cos(fRecPhi));
        fRecV.SetY(fRecR * sin(fRecPhi));
        fRecC1.SetXYZ(fRecV.X() - fRecR1 * sin(fRecPhi), fRecV.Y() + fRecR1 * cos(fRecPhi), 0.);
        fRecC2.SetXYZ(fRecV.X() + fRecR2 * sin(fRecPhi), fRecV.Y() - fRecR2 * cos(fRecPhi), 0.);
        fRecV = fRecV + fPV;
        fRecC1 = fRecC1 + fPV;
        fRecC2 = fRecC2 + fPV;
        fCutSatisfied = 0;
        if ((R1 > 0) && (R2 > 0) && (D1 > 0) && (D2 > 0) && ((Y11 - D1) > 0) && ((Y12 - D2) > 0))
          fSignSatisfied = 1;
        else
          fSignSatisfied = 0;
      }
      // B. Cut threshold NOT satisfied - prepare for next loop
      prevR1 = R1;
      prevR2 = R2;
    }
  }
}

void Conv4HitsReco2::Dump() {
  std::cout << std::endl << "================================================" << std::endl;
  std::cout << "    Nothing happend here.";
  if (fSolved == 1)
    std::cout << "Solved.";
  if (fCutSatisfied == 1)
    std::cout << "Cut good.";
  if (fSignSatisfied == 1)
    std::cout << "Sign good.";
}

math::XYZVector Conv4HitsReco2::GetPlusCenter(double &plusR) {
  plusR = fRecR1;
  return fRecC1;
}
math::XYZVector Conv4HitsReco2::GetMinusCenter(double &minusR) {
  minusR = fRecR2;
  return fRecC2;
}