#include <RhoFindOmittedParticle.h>

Inheritance diagram for RhoFindOmittedParticle:

Public Member Functions
	RhoFindOmittedParticle (const TLorentzVector &ip4Reson, const double imChild, const double imSought, const Bool_t izBoostApprox=kTRUE)

virtual	~RhoFindOmittedParticle ()

TLorentzVector	FitToSeen (const TLorentzVector &p4Obs, const TVector3 &p3Seen)

TLorentzVector	SecondVector () const

Protected Member Functions
void	MakeCone (const TLorentzVector &p4Obs)

TLorentzVector	ClosestFit (const TVector3 &p3Seen)

Protected Attributes
float	mSought2

Bool_t	zBoostApprox
	squared mass of Sought More...

float	cmEChild
	Flag: do we approximate boost in. More...

float	cmpChild2

TLorentzVector	p4ObsCache
	energy, squared 3-momentum of Child. More...

TLorentzVector	secondChoice

TVector3	beta
	The other vector, stored. More...

TVector3	cmAxis
	From lab to CM (-beta CM to lab) More...

double	openCos
	Axis of CM cone. More...

double	cmESought
	cos of opening angle. More...

double	cmpSought2
	CM energy of sought particle. More...

Detailed Description

Definition at line 64 of file RhoFindOmittedParticle.h.

Constructor & Destructor Documentation

RhoFindOmittedParticle::RhoFindOmittedParticle	(	const TLorentzVector &	ip4Reson,
		const double	imChild,
		const double	imSought,
		const Bool_t	izBoostApprox = `kTRUE`
	)

Definition at line 87 of file RhoFindOmittedParticle.cxx.

References beta, cmEChild, cmpChild2, fabs(), and mSought2.

   : mSought2 ( imSought* imSought )
   , zBoostApprox ( izBoostApprox )
   , p4ObsCache ( 0.0,0.0,0.0,-1.0 )   // Non-physical initial value
 {
   if ( imChild == 0.0 ) {
     cerr << "RhoFindOmittedParticle constructed with mChild=0" << endl;
   }
 
   if ( mSought2 == 0.0 ) {
     cerr << "RhoFindOmittedParticle constructed with mSought=0" << endl;
   }
 
   // Figure out the boost vectors to and from the center of mass
 
   //p4Reson = ip4Reson;
   beta = -ip4Reson.BoostVector();
 
   // Figure out the decay energies, gamma factors, beta factors.
   // Assume that resonance decays to two "child" particles (e.g. B + Bbar)
 
   cmEChild = ip4Reson.Mag() / 2;
 
   if ( cmEChild < fabs ( imChild ) ) {
     cerr << "RhoFindOmittedParticle: mChild too large for decay" << endl;
   }
 
   // The square of the three-momentum of the child (B0).
   cmpChild2 = cmEChild*cmEChild - imChild*imChild;
 
 }

RhoFindOmittedParticle::~RhoFindOmittedParticle ( )

virtual

Definition at line 130 of file RhoFindOmittedParticle.cxx.

131 {

132 }

Member Function Documentation

TLorentzVector RhoFindOmittedParticle::ClosestFit ( const TVector3 & p3Seen )

protected

Definition at line 247 of file RhoFindOmittedParticle.cxx.

References beta, cmAxis, cmESought, cmpSought2, fabs(), openCos, secondChoice, sqrt(), X, Y, Z, and zBoostApprox.

Referenced by FitToSeen().

 {
   // Take care of these now to save clock cycles.
   if ( fabs ( openCos ) > 1.0 ) {
     return TLorentzVector ( 0.0,0.0,0.0,0.0 );  // Error condition.
   }
 
   TLorentzVector cmp4Sought[2];  // Two intersection points...
 
   // Set up local coordinate system:
   //  z = normal z.
   //  x-z plane contains p3Seen.  Since we only keep phi information,
   //  we needn't bother boosting.
   // I'm doing more stuff by hand than in previous versions in order
   // to speed things up.
   double cmCosPlaneAxis2;   // Squared cosine between the plane and axis.
   {
     // Original, most abstract method:
     //const TVector3 cmZ(0.0, 0.0, 1.0);
     //const TVector3 cmX = (p3Seen - p3Seen.dot(cmZ) * cmZ).unit();
     //const TVector3 cmY = cmZ.cross(cmX);
 
     // Slightly faster method:
     //const TVector3 cmX = TVector3(p3Seen.x(), p3Seen.y(), 0.0).unit();
     //const TVector3 cmY = TVector3(-cmX.y(), cmX.x(), 0.0);
 
     TVector3 cmY;
     // Fastest method:
     if ( zBoostApprox ) {
       cmY = TVector3 ( -p3Seen.Y(), p3Seen.X(), 0.0 ).Unit();
     } else {
       cmY = beta.Cross ( p3Seen ).Unit();
     }
 
     const TVector3 unitAxis = cmAxis.Unit();
 
     // Projection of unit vector of axis into x-z plane.
 
     // Project unitAxis into the xz plane.
     const TVector3 unitProj = unitAxis - unitAxis.Dot ( cmY ) * cmY;
     // While we know cmY.z() == 0, doing it by hand would make the
     // code too baroque for not a great speed advantage.
 
     // Square of cosine of angle between plane and axis:
     cmCosPlaneAxis2 = unitProj.Mag2();
 
 
     if ( cmCosPlaneAxis2 < openCos*openCos ) {
       return TLorentzVector ( 0.0,0.0,0.0,0.0 );  // Error condition.
     }
 
     // Projection of full axis into x-z plane.
     // Actually, it's longer; touches the base of the cone.
     // cmAxis is projection of unitProj, hence 1/cosine.
 
     const TVector3 axisProj =
       ( cmAxis.Mag() / ( cmCosPlaneAxis2 ) ) * unitProj;
 
     // vector difference between axisProj and the true
     // vector, which could lie to either side.
 
     //const TVector3  difference = cmAxis.mag() *
     //sqrt(1/(openCos*openCos) - 1/(cmCosPlaneAxis2)) *
     //unitProj.cross(cmY).unit();
 
     const TVector3 difference =
       sqrt ( cmpSought2 - axisProj.Mag2() ) *
       unitProj.Cross ( cmY ).Unit();
 
     cmp4Sought[0] = TLorentzVector ( axisProj + difference, cmESought );
     cmp4Sought[1] = TLorentzVector ( axisProj - difference, cmESought );
   }
 
 
   // Boost the results back into the lab frame.
   // Return whichever one is closer in angle to the seen vector.
 
   cmp4Sought[0].Boost ( -beta );
   cmp4Sought[1].Boost ( -beta );
 
   if ( TVector3 ( cmp4Sought[0].X(),cmp4Sought[0].Y(),cmp4Sought[0].Z() ).Unit().Dot ( p3Seen ) >
        TVector3 ( cmp4Sought[1].X(),cmp4Sought[1].Y(),cmp4Sought[1].Z() ).Unit().Dot ( p3Seen ) ) {
     secondChoice = cmp4Sought[1];
     return cmp4Sought[0];
   } else {
     secondChoice = cmp4Sought[0];
     return cmp4Sought[1];
   }
 }

TLorentzVector RhoFindOmittedParticle::FitToSeen	(	const TLorentzVector &	p4Obs,
		const TVector3 &	p3Seen
	)

Definition at line 148 of file RhoFindOmittedParticle.cxx.

References ClosestFit(), MakeCone(), and p4ObsCache.

 {
   if ( p4Obs != p4ObsCache ) {
     this->MakeCone ( p4Obs );
     p4ObsCache = p4Obs;
   }
   return this->ClosestFit ( p3Seen );
 }

void RhoFindOmittedParticle::MakeCone ( const TLorentzVector & p4Obs )

protected

Definition at line 166 of file RhoFindOmittedParticle.cxx.

References beta, cmAxis, cmEChild, cmESought, cmpChild2, cmpSought2, mSought2, openCos, and sqrt().

Referenced by FitToSeen().

 {
   // Get the observed mass, figure out the decay momentum of the
   // child -> sought + observed.
 
   // Square of mass of observed particle.
   //const double mObs2 = p4Obs.Mag2();
 
   /* This was an old, labor-intensive way of calculating the
      decay momentum.  Using conservation of energy is much simpler to
      understand and faster to compute.
   // The following is a solution of sqrt(pDecay*pDecay + mSought*mSought)
   // + sqrt(pDecay*pDecay + mObs*mObs) == mChild
   double pDecay;
   {
   double temp;
 
   temp = (mObs2 - mSought*mSought) / mChild;
   temp *= temp;
   temp += mChild*mChild - 2 * (mObs2 + mSought*mSought);
 
   pDecay = 0.5 * sqrt(temp);
   }
 
   //const double EObsDecay = sqrt(mObs2 + pDecay*pDecay);
   //const double ESoughtDecay = sqrt(mSought*mSought + pDecay*pDecay);
   */
 
 
   // Boost vectors into the CM frame; extract vital statistics.
 
   TLorentzVector cmp4Obs ( p4Obs );
   cmp4Obs.Boost ( beta );
 
   //const double cmEObs = cmp4Obs.t();
   const double cmpObs2 = TVector3 ( cmp4Obs.X(),cmp4Obs.Y(),cmp4Obs.Z() ).Mag2();
 
   cmESought = cmEChild - cmp4Obs.T();
   double cmESought2 = cmESought*cmESought;
   if ( cmESought2 > mSought2 ) {
     cmpSought2 = cmESought2 - mSought2;
   } else {
     cmpSought2 = 0.0;
   }
 
 
 
   /*
    * Because we know the lengths of the 3 momentum 3-vectors, we can
    * get the angles using the law of cosines.
    * abs(openCos) > 1 signifies an error, causing closestFit to return
    * a zero vector, signifying an error.
    */
 
   const double cmpSought = sqrt ( cmpSought2 );
 
   if ( cmpSought != 0 )
     openCos =
       ( cmpObs2 + cmpSought2 - cmpChild2 ) /
       ( 2 * sqrt ( cmpObs2 ) * cmpSought );
   else {
     openCos = 2;  // Error condition.
   }
 
   cmAxis = -TVector3 ( cmp4Obs.X(),cmp4Obs.Y(),cmp4Obs.Z() ).Unit() * openCos * cmpSought;
 
 }