#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 31 of file PFGsfHelper.h.

Constructor & Destructor Documentation

PFGsfHelper::PFGsfHelper ( const TrajectoryMeasurement & tm )

Definition at line 31 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.;
    std::vector<TrajectoryStateOnSurface> components(theUpdateState.components());
    unsigned int 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 ( std::vector<TrajectoryStateOnSurface>::const_iterator 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 128 of file PFGsfHelper.cc.

{
}

Member Function Documentation

GlobalVector PFGsfHelper::computeP ( bool ComputeMode ) const

Definition at line 130 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 149 of file PFGsfHelper.cc.

References i, TrajectoryStateOnSurface::localError(), TrajectoryStateOnSurface::localParameters(), LocalTrajectoryError::matrix(), SingleGaussianState1D::mean(), GaussianSumUtilities1D::mode(), GaussianSumUtilities1D::modeIsValid(), MultiGaussianStateTransform::multiState1D(), Parameters::parameters, 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;
}