#include <Conv4HitsReco.h>

Public Member Functions
int	ComputeMaxLimits ()
int	ComputeMinLimits ()
	Conv4HitsReco (TVector3 &, TVector3 &, TVector3 &, TVector3 &, TVector3 &)
int	ConversionCandidate (TVector3 &, double &, double &)
void	Dump ()
TVector3	GetConvVertexFromParams (double &, double &)
double	GetDm ()
double	GetDq ()
TVector3	GetIntersection (TVector3 &, TVector3 &, TVector3 &, TVector3 &)
TVector3	GetMinusCenter (double &)
TVector3	GetPlusCenter (double &)
double	GetPtFromParamAndHitPair (double &, double &)
double	GetPtMinusFromParam (double &)
double	GetPtPlusFromParam (double &)
int	GuessStartingValues (double &, double &)
int	IsNotValidForPtLimit (double, double, double, double)
int	IsNotValidForVtxPosition (double &)
int	mqFindByIteration (double &, double &)
double	qFromM (double)
std::string	qFromM_print (double m)
void	Set (double)
void	SetIterationStopRelThreshold (double val)
void	SetLinSystCoeff (double, double)
void	SetMaxNumberOfIterations (int val)
void	SetMaxVtxDistance (int val)
void	SetPtLegMaxCut (double val)
void	SetPtLegMinCut (double val)
void	SetPtPhotMaxCut (double val)
	~Conv4HitsReco ()
Public Attributes
double	_eta
double	_eta2
double	_pi
double	_pi2
TVector3	convVtx
TVector3	hit1
TVector3	hit2
TVector3	hit3
TVector3	hit4
double	iterationStopRelThreshold
int	maxNumberOfIterations
double	maxVtxDistance
TVector3	minusCenter
double	minusRadius
double	mMaxLimit
double	mMinLimit
double	nNV
double	nPN
double	O0
double	Om
double	Oq
TVector3	plusCenter
double	plusRadius
double	PN
double	pn
double	PN2
double	pNV
double	pPN
double	ptLegMaxCut
double	ptLegMinCut
double	ptPhotMaxCut
double	qMaxLimit
double	qMinLimit
double	T0
double	Tm
double	Tq
TVector3	unitVn
TVector3	unitVp
TVector3	v1minus2
TVector3	v3minus4
TVector3	vMMaxLimit
TVector3	vN
TVector3	vNminusV
TVector3	vP
TVector3	vPminusN
TVector3	vPminusV
TVector3	vQMaxLimit
TVector3	vV

Detailed Description

Definition at line 10 of file Conv4HitsReco.h.

Constructor & Destructor Documentation

Conv4HitsReco::Conv4HitsReco	(	TVector3 &	,
		TVector3 &	,
		TVector3 &	,
		TVector3 &	,
		TVector3 &
	)

Conv4HitsReco::~Conv4HitsReco ( )

Member Function Documentation

int Conv4HitsReco::ComputeMaxLimits ( )

Definition at line 208 of file Conv4HitsReco.cc.

References _eta, _pi, GetIntersection(), GetPtFromParamAndHitPair(), hit2, hit3, IsNotValidForPtLimit(), LogDebug, mMaxLimit, ptLegMinCut, qMaxLimit, unitVn, unitVp, vMMaxLimit, vN, vP, vQMaxLimit, and vV.

Referenced by ConversionCandidate().

                                   {

  // 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::ComputeMinLimits ( )

Definition at line 246 of file Conv4HitsReco.cc.

References GetIntersection(), IsNotValidForPtLimit(), LogDebug, mMaxLimit, mMinLimit, ptLegMaxCut, qMaxLimit, qMinLimit, unitVn, unitVp, v1minus2, v3minus4, vMMaxLimit, vN, vP, vQMaxLimit, vV, and Gflash::Z.

Referenced by ConversionCandidate().

                                   {

  //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::ConversionCandidate	(	TVector3 &	vtx,
		double &	ptplus,
		double &	ptminus
	)

Definition at line 135 of file Conv4HitsReco.cc.

References ComputeMaxLimits(), ComputeMinLimits(), convVtx, GetConvVertexFromParams(), GetPtMinusFromParam(), GetPtPlusFromParam(), GuessStartingValues(), hit1, hit4, IsNotValidForPtLimit(), IsNotValidForVtxPosition(), LogDebug, m, maxVtxDistance, minusCenter, minusRadius, mMaxLimit, mMinLimit, mqFindByIteration(), plusCenter, plusRadius, ptLegMaxCut, ptLegMinCut, lumiQueryAPI::q, qMaxLimit, qMinLimit, unitVn, unitVp, vN, and vP.

Referenced by SeedForPhotonConversionFromQuadruplets::trajectorySeed().

                                                                                    {

  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;

}

void Conv4HitsReco::Dump ( void )

Definition at line 292 of file Conv4HitsReco.cc.

References _eta, _pi, hit1, hit2, hit3, hit4, LogDebug, unitVn, unitVp, vN, vP, vPminusN, and vV.

Referenced by SeedForPhotonConversionFromQuadruplets::trajectorySeed().

                        {

  LogDebug("Conv4HitsReco")
    << " ======================================= " << "\n"
    << " Photon Vertex: "
    << vV
    << " Hit1: "
    << hit1
    << " Hit2: "
    << hit2
    << " Hit3: "
    << hit3
    << " Hit4: "
    << hit4
    << " N: "
    << vN
    << " P: "
    << vP
    << " P-N: "
    << vPminusN
    << " n: "
    << unitVn
    << " p: "
    << unitVp
    << " eta: " << _eta << " pi: " << _pi << "\n";

}

TVector3 Conv4HitsReco::GetConvVertexFromParams	(	double &	m,
		double &	q
	)

Definition at line 119 of file Conv4HitsReco.cc.

References LogDebug, minusCenter, minusRadius, plusCenter, plusRadius, unitVn, unitVp, vN, and vP.

Referenced by ConversionCandidate().

                                                                   {

  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
                             << ">>>>>>>> Conversion vertex computed via Minus pair\n"
                             << vtxViaMinus;

  return 0.5*(vtxViaPlus+vtxViaMinus);

}

double Conv4HitsReco::GetDm ( )

Definition at line 10 of file Conv4HitsReco.cc.

References O0, Om, Oq, T0, and Tq.

Referenced by mqFindByIteration().

                           {

  return (Tq*O0-T0*Oq)/(Tq*Om-Tm*Oq); //dm
  
}

double Conv4HitsReco::GetDq ( )

Definition at line 16 of file Conv4HitsReco.cc.

References O0, Om, Oq, T0, and Tq.

Referenced by mqFindByIteration().

                           {

  return (Tm*O0-T0*Om)/(Tm*Oq-Tq*Om); //dq
  
}

TVector3 Conv4HitsReco::GetIntersection	(	TVector3 &	V1,
		TVector3 &	p1,
		TVector3 &	V2,
		TVector3 &	p2
	)

Definition at line 91 of file Conv4HitsReco.cc.

References p2.

Referenced by ComputeMaxLimits(), and ComputeMinLimits().

                                                                                             {

  TVector3 v1 = V1.Unit();
  TVector3 v2 = V2.Unit();
  
  double v1v2=v1*v2;
  return v2*((p1 - p2)*(v1v2*v1-v2)/(v1v2*v1v2-1.))+p2;

}

TVector3 Conv4HitsReco::GetMinusCenter ( double & radius )

Definition at line 84 of file Conv4HitsReco.cc.

References minusCenter, and minusRadius.

Referenced by SeedForPhotonConversionFromQuadruplets::trajectorySeed().

                                                    {
  radius = minusRadius;
  return minusCenter;
}

TVector3 Conv4HitsReco::GetPlusCenter ( double & radius )

Definition at line 79 of file Conv4HitsReco.cc.

References plusCenter, and plusRadius.

Referenced by SeedForPhotonConversionFromQuadruplets::trajectorySeed().

                                                   {
  radius = plusRadius;
  return plusCenter;
}

double Conv4HitsReco::GetPtFromParamAndHitPair	(	double &	m,
		double &	eta
	)

Definition at line 101 of file Conv4HitsReco.cc.

References mathSSE::sqrt().

Referenced by ComputeMaxLimits(), GetPtMinusFromParam(), GetPtPlusFromParam(), and IsNotValidForPtLimit().

                                                                    {

  return 0.01*0.3*3.8*sqrt(m*m+eta*eta);

}

double Conv4HitsReco::GetPtMinusFromParam ( double & m )

Definition at line 107 of file Conv4HitsReco.cc.

References _eta, and GetPtFromParamAndHitPair().

Referenced by ConversionCandidate().

                                                  {

  return GetPtFromParamAndHitPair(m, _eta);

}

double Conv4HitsReco::GetPtPlusFromParam ( double & q )

Definition at line 113 of file Conv4HitsReco.cc.

References _pi, and GetPtFromParamAndHitPair().

Referenced by ConversionCandidate().

                                                 {

  return GetPtFromParamAndHitPair(q, _pi);

}

int Conv4HitsReco::GuessStartingValues	(	double &	m,
		double &	q
	)

Definition at line 275 of file Conv4HitsReco.cc.

References LogDebug, mMaxLimit, and qMaxLimit.

Referenced by ConversionCandidate().

                                                          {

  /*
  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;

}

int Conv4HitsReco::IsNotValidForPtLimit	(	double	m,
		double	q,
		double	ptmin,
		double	ptmax
	)

Definition at line 228 of file Conv4HitsReco.cc.

References _eta, _pi, and GetPtFromParamAndHitPair().

Referenced by ComputeMaxLimits(), ComputeMinLimits(), and ConversionCandidate().

                                                                                     {

  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 )

Definition at line 238 of file Conv4HitsReco.cc.

References convVtx, hit1, hit2, hit3, hit4, and insertMaterial::maxDist.

Referenced by ConversionCandidate().

                                                          {

  TVector3 hitAve = 0.25*(hit1+hit2+hit3+hit4); 
  if ( (convVtx-hitAve).Mag() > maxDist ) return 1;
  return 0;

}

int Conv4HitsReco::mqFindByIteration	(	double &	m,
		double &	q
	)

Definition at line 177 of file Conv4HitsReco.cc.

References GetDm(), GetDq(), i, iterationStopRelThreshold, LogDebug, maxNumberOfIterations, and SetLinSystCoeff().

Referenced by ConversionCandidate().

                                                        {

  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;

}

double Conv4HitsReco::qFromM ( double m )

Definition at line 72 of file Conv4HitsReco.cc.

References LogDebug, and qFromM_print().

                                    {

  LogDebug("Conv4HitsReco")<<qFromM_print(m);
  return 0.;

}

std::string Conv4HitsReco::qFromM_print ( double m )

Definition at line 50 of file Conv4HitsReco.cc.

References _eta2, alpha, beta, m, pn, pNV, mathSSE::sqrt(), unitVn, unitVp, vNminusV, and vPminusN.

Referenced by qFromM().

                                             {
  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();
}

void Conv4HitsReco::Set ( double )

void Conv4HitsReco::SetIterationStopRelThreshold ( double val ) [inline]

Definition at line 78 of file Conv4HitsReco.h.

References iterationStopRelThreshold.

{iterationStopRelThreshold = val;};

void Conv4HitsReco::SetLinSystCoeff	(	double	m,
		double	q
	)

Definition at line 22 of file Conv4HitsReco.cc.

References _eta2, _pi2, m, nNV, nPN, O0, Om, Oq, pn, PN2, pNV, pPN, mathSSE::sqrt(), T0, Tq, unitVn, unitVp, vNminusV, and vPminusN.

Referenced by mqFindByIteration().

                                                     {
  
  // 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;

}