Conv4HitsReco.cc

Go to the documentation of this file.

#define Conv4HitsReco_cxx
#include "Conv4HitsReco.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <sstream>
// 1 -> 4
// 2 -> 3
// 3 -> 2
// 4 -> 1

double Conv4HitsReco::GetDm() {
  return (Tq * O0 - T0 * Oq) / (Tq * Om - Tm * Oq);  //dm
}

double Conv4HitsReco::GetDq() {
  return (Tm * O0 - T0 * Om) / (Tm * Oq - Tq * Om);  //dq
}

void Conv4HitsReco::SetLinSystCoeff(double m, double q) {
  // dq*Tq + dm*Tm = T0 // Tangent condition
  // dq*Oq + dm*Om = O0 // Orthogonality condition

  double sqrtEta2mm = sqrt(_eta2 + m * m);
  double sqrtPi2qq = sqrt(_pi2 + q * q);

  double signT = 1.;

  Tq = -2. * pPN + 2. * m * pn + signT * 2. * q * sqrtEta2mm / sqrtPi2qq;
  Tm = 2. * nPN + 2. * q * pn + signT * 2. * m * sqrtPi2qq / sqrtEta2mm;
  T0 = PN2 - _eta2 - _pi2 - 2. * q * m * pn + 2. * q * pPN - 2. * m * nPN - signT * 2. * sqrtEta2mm * sqrtPi2qq;

  TVector3 vQminusM = q * unitVp - m * unitVn + vPminusN;
  double QM = vQminusM.Mag();
  double pQM = unitVp * vQminusM;
  double nQM = unitVn * vQminusM;
  double NVQM = vNminusV * vQminusM;

  double signO = 1.;

  Oq = sqrtEta2mm * pQM / QM + m * pn + pNV;
  Om = m * QM / sqrtEta2mm - signO * sqrtEta2mm * nQM / QM + nQM - nNV - m;
  O0 = -signO * sqrtEta2mm * QM - m * nQM - NVQM;
}

std::string Conv4HitsReco::qFromM_print(double m) {
  std::stringstream ss;
  TVector3 vPminusM = vPminusN - m * unitVn;
  double m2 = m * m;
  double nPM = unitVn * vPminusM;
  double pPM = unitVp * vPminusM;
  double NVPM = vNminusV * vPminusM;

  double alpha = (m * pn + pNV) * (m * pn + pNV) - _eta2 - m2;
  double beta = m2 * pn * nPM + m * pn * NVPM + m * nPM * pNV - pPM * (_eta2 + m2) + pNV * NVPM;
  double gamma = m2 * nPM * nPM + NVPM * NVPM + 2. * m * nPM * NVPM - vPminusM.Mag2() * (_eta2 + m2);

  double disc = sqrt(beta * beta - alpha * gamma);

  double q01 = (-beta + disc) / alpha;
  double q02 = (-beta - disc) / alpha;

  ss << " m: " << m << " q01: " << std::setw(20) << q01 << " q02: " << std::setw(20) << q02 << "/n";
  return ss.str();
}

double Conv4HitsReco::qFromM(double m) {
  LogDebug("Conv4HitsReco") << qFromM_print(m);
  return 0.;
}

TVector3 Conv4HitsReco::GetPlusCenter(double &radius) {
  radius = plusRadius;
  return plusCenter;
}

TVector3 Conv4HitsReco::GetMinusCenter(double &radius) {
  radius = minusRadius;
  return minusCenter;
}

//
//Point of intersection between two lines (each identified by a vector and a point)
TVector3 Conv4HitsReco::GetIntersection(TVector3 &V1, TVector3 &p1, TVector3 &V2, TVector3 &p2) {
  TVector3 v1 = V1.Unit();
  TVector3 v2 = V2.Unit();

  double v1v2 = v1 * v2;
  return v2 * ((p1 - p2) * (v1v2 * v1 - v2) / (v1v2 * v1v2 - 1.)) + p2;
}

double Conv4HitsReco::GetPtFromParamAndHitPair(double &m, double &eta) {
  return 0.01 * 0.3 * 3.8 * sqrt(m * m + eta * eta);
}

double Conv4HitsReco::GetPtMinusFromParam(double &m) { return GetPtFromParamAndHitPair(m, _eta); }

double Conv4HitsReco::GetPtPlusFromParam(double &q) { return GetPtFromParamAndHitPair(q, _pi); }

TVector3 Conv4HitsReco::GetConvVertexFromParams(double &m, double &q) {
  TVector3 unitVQminusM = (plusCenter - minusCenter).Unit();
  TVector3 vtxViaPlus = vP + q * unitVp - plusRadius * unitVQminusM;
  TVector3 vtxViaMinus = vN + m * unitVn + minusRadius * unitVQminusM;

  //  return 0.5*(vN+m*unitVn+m*unitVQminusM+vP+q*unitVp-q*unitVQminusM);
  LogDebug("Conv4HitsReco") << ">>>>>>>> Conversion vertex computed via Plus pair\n"
                            << vtxViaPlus.x() << "," << vtxViaPlus.y() << "," << vtxViaPlus.z()
                            << ">>>>>>>> Conversion vertex computed via Minus pair\n"
                            << vtxViaMinus.x() << "," << vtxViaMinus.y() << "," << vtxViaMinus.z();

  return 0.5 * (vtxViaPlus + vtxViaMinus);
}

int Conv4HitsReco::ConversionCandidate(TVector3 &vtx, double &ptplus, double &ptminus) {
  double m;
  double q;
  int nits = 0;

  int isNotValidBefore = ComputeMaxLimits() + ComputeMinLimits();
  int isNotValidAfter = 0;
  if (!isNotValidBefore) {
    GuessStartingValues(m, q);
    nits = mqFindByIteration(m, q);

    if (q > qMaxLimit || q < qMinLimit) {
      isNotValidAfter = 1;
      LogDebug("Conv4HitsReco") << ">>>>>>>> quad result not valid for q: qMin= " << qMinLimit << " q= " << q
                                << " qMax=  " << qMaxLimit << "\n";
    }
    if (m > mMaxLimit || m < mMinLimit) {
      isNotValidAfter = 1;
      LogDebug("Conv4HitsReco") << ">>>>>>>> quad result not valid for m: mMin= " << mMinLimit << " m= " << m
                                << " mMax=  " << mMaxLimit << "\n";
    }

    ptminus = GetPtMinusFromParam(m);
    ptplus = GetPtPlusFromParam(q);
    minusCenter = vN + m * unitVn;
    minusRadius = (hit4 - minusCenter).Mag();
    plusCenter = vP + q * unitVp;
    plusRadius = (hit1 - plusCenter).Mag();
    convVtx = GetConvVertexFromParams(m, q);
    vtx = convVtx;
  }

  if (isNotValidBefore)
    return 0;
  if (IsNotValidForPtLimit(m, q, ptLegMinCut, ptLegMaxCut))
    return -1000;
  if (IsNotValidForVtxPosition(maxVtxDistance))
    return -2000;
  if (isNotValidAfter)
    return -1 * nits;
  return nits;
}

int Conv4HitsReco::mqFindByIteration(double &m, double &q) {
  int maxIte = maxNumberOfIterations;
  double err = iterationStopRelThreshold;
  double edm = 1.;
  double edq = 1.;
  int i = 0;
  while (((edq > err) || (edm > err)) && (i < maxIte)) {
    SetLinSystCoeff(m, q);
    double dm = GetDm();
    double dq = GetDq();
    /*    
      while( m+dm > mMaxLimit || m+dm < mMinLimit || q+dq > qMaxLimit || q+dq < qMinLimit ){

      LogDebug("Conv4HitsReco")<<">>>>>>>> Going outside limits, reducing increments \n";
      dm=dm/2.;
      dq=dq/2.;
      }
    */
    m += dm;
    q += dq;
    edm = fabs(dm / m);
    edq = fabs(dq / q);
    LogDebug("Conv4HitsReco") << ">>>>>>>> Iteration " << i << " m: " << m << " q: " << q << " dm: " << dm
                              << " dq: " << dq << " edm: " << edm << " edq: " << edq << "\n";
    i++;
  }

  return i;
}

int Conv4HitsReco::ComputeMaxLimits() {
  // Q max limit
  TVector3 vVminusHit2Orthogonal = (vV - hit2).Orthogonal();
  TVector3 medianPointHit2V = 0.5 * (vV + hit2);
  vQMaxLimit = GetIntersection(vVminusHit2Orthogonal, medianPointHit2V, unitVp, vP);
  qMaxLimit = (vQMaxLimit - vP).Mag();

  // M max limit
  TVector3 vVminusHit3Orthogonal = (vV - hit3).Orthogonal();
  TVector3 medianPointHit3V = 0.5 * (vV + hit3);
  vMMaxLimit = GetIntersection(vVminusHit3Orthogonal, medianPointHit3V, unitVn, vN);
  mMaxLimit = (vMMaxLimit - vN).Mag();

  LogDebug("Conv4HitsReco") << " >>>>>> Limits: qMax= " << qMaxLimit << " ==>pt "
                            << GetPtFromParamAndHitPair(qMaxLimit, _pi) << " mMax= " << mMaxLimit << " ==>pt "
                            << GetPtFromParamAndHitPair(mMaxLimit, _eta) << "\n";

  return IsNotValidForPtLimit(mMaxLimit, qMaxLimit, ptLegMinCut, 100000.);  //Max limit not applied here
}

int Conv4HitsReco::IsNotValidForPtLimit(double m, double q, double ptmin, double ptmax) {
  if (GetPtFromParamAndHitPair(q, _pi) < ptmin)
    return 1;
  if (GetPtFromParamAndHitPair(m, _eta) < ptmin)
    return 1;
  if (GetPtFromParamAndHitPair(q, _pi) > ptmax)
    return 1;
  if (GetPtFromParamAndHitPair(m, _eta) > ptmax)
    return 1;
  return 0;
}

int Conv4HitsReco::IsNotValidForVtxPosition(double &maxDist) {
  TVector3 hitAve = 0.25 * (hit1 + hit2 + hit3 + hit4);
  if ((convVtx - hitAve).Mag() > maxDist)
    return 1;
  return 0;
}

int Conv4HitsReco::ComputeMinLimits() {
  //Evaluate if quad is valid and compute min limits
  if (((vV - vQMaxLimit).Cross(vMMaxLimit - vQMaxLimit)).Z() > 0.) {
    //
    //Quad is invalid
    LogDebug("Conv4HitsReco") << " >>>>>> Quad is invalid\n";
    return 1;
  } else {
    //
    // Compute q and m Min limits
    TVector3 vQMinLimit = GetIntersection(v1minus2, vMMaxLimit, unitVp, vP);
    qMinLimit = (vQMinLimit - vP) * unitVp;
    TVector3 vMMinLimit = GetIntersection(v3minus4, vQMaxLimit, unitVn, vN);
    mMinLimit = (vMMinLimit - vN) * unitVn;
    if (mMinLimit > mMaxLimit || qMinLimit > qMaxLimit) {
      LogDebug("Conv4HitsReco") << " >>>>>> Quad is invalid. qMin= " << qMinLimit << " mMin= " << mMinLimit << "\n";
      return 2;
    }
    if (IsNotValidForPtLimit(mMinLimit, qMinLimit, -1000., ptLegMaxCut)) {  //Min limit not applied here
      return 2;
    }

    LogDebug("Conv4HitsReco") << " >>>>>> Quad is valid. qMin= " << qMinLimit << " mMin= " << mMinLimit << "\n";
    return 0;
  }
}

int Conv4HitsReco::GuessStartingValues(double &m, double &q) {
  /*
  m = 0.5*(mMinLimit+mMaxLimit);
  q = 0.5*(qMinLimit+qMaxLimit);
  */

  m = mMaxLimit;
  q = qMaxLimit;

  LogDebug("Conv4HitsReco") << " >>>>>> Starting values: q= " << q << " m= " << m << "\n";

  return 0;
}

void Conv4HitsReco::Dump() {
  LogDebug("Conv4HitsReco") << " ======================================= "
                            << "\n"
                            << " Photon Vertex: " << vV.x() << "," << vV.y() << "," << vV.z() << " Hit1: " << hit1.x()
                            << "," << hit1.y() << "," << hit1.z() << " Hit2: " << hit2.x() << "," << hit2.y() << ","
                            << hit2.z() << " Hit3: " << hit3.x() << "," << hit3.y() << "," << hit3.z()
                            << " Hit4: " << hit4.x() << "," << hit4.y() << "," << hit4.z() << " N: " << vN.x() << ","
                            << vN.y() << "," << vN.z() << " P: " << vP.x() << "," << vP.y() << "," << vP.z()
                            << " P-N: " << vPminusN.x() << "," << vP.y() << "," << vP.z() << " n: " << unitVn.x() << ","
                            << unitVn.y() << "," << unitVn.z() << " p: " << unitVp.x() << "," << unitVp.y() << ","
                            << unitVp.z() << " eta: " << _eta << " pi: " << _pi << "\n";
}