#include <PFGsfHelper.h>

Public Member Functions
GlobalVector	computeP (bool ComputeMode) const

double	fittedDP () const

bool	isValid () const

	PFGsfHelper (const TrajectoryMeasurement &)

double	sigmafittedDP () const

	~PFGsfHelper ()

Private Member Functions
void	computeQpMode (const TrajectoryStateOnSurface tsos, AlgebraicVector5 &parameters, AlgebraicSymMatrix55 &covariance) const

Private Attributes
double	dp

float	mode_Px

float	mode_Py

float	mode_Pz

double	sigmaDp

TrajectoryStateOnSurface	theBackwardState

TrajectoryStateOnSurface	theForwardState

TrajectoryStateOnSurface	theUpdateState

bool	Valid

Detailed Description

Definition at line 29 of file PFGsfHelper.h.

Constructor & Destructor Documentation

PFGsfHelper::PFGsfHelper ( const TrajectoryMeasurement & tm )

Definition at line 32 of file PFGsfHelper.cc.

References TrajectoryMeasurement::backwardPredictedState(), makeMuonMisalignmentScenario::components, CollinearFitAtTM::fit(), fitChi2(), TrajectoryMeasurement::forwardPredictedState(), SingleGaussianState1D::mean(), GaussianSumUtilities1D::mode(), TrajectoryMeasurement::recHit(), mathSSE::sqrt(), and TrajectoryMeasurement::updatedState().

                                                         {
   /* LogInfo("PFGsfHelper")<<" PFGsfHelper  built"; */
 
   // TrajectoryStateOnSurface theUpdateState = tm.forwardPredictedState();
   theUpdateState = tm.updatedState();
   theForwardState = tm.forwardPredictedState();
   theBackwardState = tm.backwardPredictedState();
 
   Valid = true;
   if (!theUpdateState.isValid() || !theForwardState.isValid() || !theBackwardState.isValid())
     Valid = false;
 
   if (Valid) {
     mode_Px = 0.;
     mode_Py = 0.;
     mode_Pz = 0.;
     GetComponents comps(theUpdateState);
     auto const& components = comps();
     auto numb = components.size();
     std::vector<SingleGaussianState1D> pxStates;
     pxStates.reserve(numb);
     std::vector<SingleGaussianState1D> pyStates;
     pyStates.reserve(numb);
     std::vector<SingleGaussianState1D> pzStates;
     pzStates.reserve(numb);
     for (auto ic = components.begin(); ic != components.end(); ++ic) {
       GlobalVector momentum(ic->globalMomentum());
       AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
       pxStates.push_back(SingleGaussianState1D(momentum.x(), cov(3, 3), ic->weight()));
       pyStates.push_back(SingleGaussianState1D(momentum.y(), cov(4, 4), ic->weight()));
       pzStates.push_back(SingleGaussianState1D(momentum.z(), cov(5, 5), ic->weight()));
       //    cout<<"COMP "<<momentum<<endl;
     }
     MultiGaussianState1D pxState(pxStates);
     MultiGaussianState1D pyState(pyStates);
     MultiGaussianState1D pzState(pzStates);
     GaussianSumUtilities1D pxUtils(pxState);
     GaussianSumUtilities1D pyUtils(pyState);
     GaussianSumUtilities1D pzUtils(pzState);
     mode_Px = pxUtils.mode().mean();
     mode_Py = pyUtils.mode().mean();
     mode_Pz = pzUtils.mode().mean();
 
     dp = 0.;
     sigmaDp = 0.;
 
     //
     // prepare input parameter vectors and covariance matrices
     //
 
     AlgebraicVector5 fwdPars = theForwardState.localParameters().vector();
     AlgebraicSymMatrix55 fwdCov = theForwardState.localError().matrix();
     computeQpMode(theForwardState, fwdPars, fwdCov);
     AlgebraicVector5 bwdPars = theBackwardState.localParameters().vector();
     AlgebraicSymMatrix55 bwdCov = theBackwardState.localError().matrix();
     computeQpMode(theBackwardState, bwdPars, bwdCov);
     LocalPoint hitPos(0., 0., 0.);
     LocalError hitErr(-1., -1., -1.);
     if (tm.recHit()->isValid()) {
       hitPos = tm.recHit()->localPosition();
       hitErr = tm.recHit()->localPositionError();
     }
 
     CollinearFitAtTM collinearFit;
     CollinearFitAtTM::ResultVector fitParameters;
     CollinearFitAtTM::ResultMatrix fitCovariance;
     double fitChi2;
     bool CollFit = true;
     if (!collinearFit.fit(fwdPars, fwdCov, bwdPars, bwdCov, hitPos, hitErr, fitParameters, fitCovariance, fitChi2))
       CollFit = false;
 
     if (CollFit) {
       double qpIn = fitParameters(0);
       double sig2In = fitCovariance(0, 0);
       double qpOut = fitParameters(1);
       double sig2Out = fitCovariance(1, 1);
       double corrInOut = fitCovariance(0, 1);
       double pIn = 1. / fabs(qpIn);
       double pOut = 1. / fabs(qpOut);
       double sig2DeltaP = pIn / qpIn * pIn / qpIn * sig2In - 2 * pIn / qpIn * pOut / qpOut * corrInOut +
                           pOut / qpOut * pOut / qpOut * sig2Out;
       //   std::cout << "fitted delta p = " << pOut-pIn << " sigma = "
       //        << sqrt(sig2DeltaP) << std::endl;
       dp = pOut - pIn;
       sigmaDp = sqrt(sig2DeltaP);
     }
   }
 }

PFGsfHelper::~PFGsfHelper ( )

Definition at line 121 of file PFGsfHelper.cc.

121 {}

Member Function Documentation

GlobalVector PFGsfHelper::computeP ( bool ComputeMode ) const

Definition at line 122 of file PFGsfHelper.cc.

Referenced by PFTrackTransformer::addPointsAndBrems().

                                                          {
   GlobalVector gsfp;
   if (ComputeMode)
     gsfp = GlobalVector(mode_Px, mode_Py, mode_Pz);
   else
     gsfp = theUpdateState.globalMomentum();
   return gsfp;
 }

void PFGsfHelper::computeQpMode	(	const TrajectoryStateOnSurface	tsos,
		AlgebraicVector5 &	parameters,
		AlgebraicSymMatrix55 &	covariance
	)		const

private

Definition at line 134 of file PFGsfHelper.cc.

References mps_fire::i, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), LocalTrajectoryError::matrix(), SingleGaussianState1D::mean(), GaussianSumUtilities1D::mode(), GaussianSumUtilities1D::modeIsValid(), MultiGaussianStateTransform::multiState1D(), mathSSE::sqrt(), SingleGaussianState1D::variance(), and LocalTrajectoryParameters::vector().

                                                                         {
   //
   // parameters and errors from combined state
   //
   parameters = tsos.localParameters().vector();
   covariance = tsos.localError().matrix();
   //
   // mode for parameter 0 (q/p)
   //
   MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(tsos, 0));
   GaussianSumUtilities1D qpGS(qpState);
   if (!qpGS.modeIsValid())
     return;
   double qp = qpGS.mode().mean();
   double varQp = qpGS.mode().variance();
   //
   // replace q/p value and variance, rescale correlation terms
   //   (heuristic procedure - alternative would be mode in 5D ...)
   //
   double VarQpRatio = sqrt(varQp / covariance(0, 0));
   parameters(0) = qp;
   covariance(0, 0) = varQp;
   for (int i = 1; i < 5; ++i)
     covariance(i, 0) *= VarQpRatio;
   //   std::cout << "covariance " << VarQpRatio << " "
   //        << covariance(1,0) << " " << covariance(0,1) << std::endl;
 }