PFTracking/src/PFGsfHelper.cc

0001 //
0002 // -*- C++ -*-
0003 // Package:    PFTracking
0004 // Class:      PFGsfHelper
0005 //
0006 // Original Author:  Daniele Benedetti
0007
0008 #include "RecoParticleFlow/PFTracking/interface/PFGsfHelper.h"
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012
0013 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0014 #include "DataFormats/ParticleFlowReco/interface/PFRecTrackFwd.h"
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016
0017 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0018 //#include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0019
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0021 // Add by Daniele
0022 #include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
0023 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
0024 #include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"
0025 #include "RecoParticleFlow/PFTracking/interface/CollinearFitAtTM.h"
0026 #include "TrackingTools/GsfTools/interface/GetComponents.h"
0027
0028 using namespace std;
0029 using namespace reco;
0030 using namespace edm;
0031
0032 PFGsfHelper::PFGsfHelper(const TrajectoryMeasurement& tm) {
0033   /* LogInfo("PFGsfHelper")<<" PFGsfHelper  built"; */
0034
0035   // TrajectoryStateOnSurface theUpdateState = tm.forwardPredictedState();
0036   theUpdateState = tm.updatedState();
0037   theForwardState = tm.forwardPredictedState();
0038   theBackwardState = tm.backwardPredictedState();
0039
0040   Valid = true;
0041   if (!theUpdateState.isValid() || !theForwardState.isValid() || !theBackwardState.isValid())
0042     Valid = false;
0043
0044   if (Valid) {
0045     mode_Px = 0.;
0046     mode_Py = 0.;
0047     mode_Pz = 0.;
0048     GetComponents comps(theUpdateState);
0049     auto const& components = comps();
0050     auto numb = components.size();
0051     std::vector<SingleGaussianState1D> pxStates;
0052     pxStates.reserve(numb);
0053     std::vector<SingleGaussianState1D> pyStates;
0054     pyStates.reserve(numb);
0055     std::vector<SingleGaussianState1D> pzStates;
0056     pzStates.reserve(numb);
0057     for (auto ic = components.begin(); ic != components.end(); ++ic) {
0058       GlobalVector momentum(ic->globalMomentum());
0059       AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
0060       pxStates.push_back(SingleGaussianState1D(momentum.x(), cov(3, 3), ic->weight()));
0061       pyStates.push_back(SingleGaussianState1D(momentum.y(), cov(4, 4), ic->weight()));
0062       pzStates.push_back(SingleGaussianState1D(momentum.z(), cov(5, 5), ic->weight()));
0063       //    cout<<"COMP "<<momentum<<endl;
0064     }
0065     MultiGaussianState1D pxState(pxStates);
0066     MultiGaussianState1D pyState(pyStates);
0067     MultiGaussianState1D pzState(pzStates);
0068     GaussianSumUtilities1D pxUtils(pxState);
0069     GaussianSumUtilities1D pyUtils(pyState);
0070     GaussianSumUtilities1D pzUtils(pzState);
0071     mode_Px = pxUtils.mode().mean();
0072     mode_Py = pyUtils.mode().mean();
0073     mode_Pz = pzUtils.mode().mean();
0074
0075     dp = 0.;
0076     sigmaDp = 0.;
0077
0078     //
0079     // prepare input parameter vectors and covariance matrices
0080     //
0081
0082     AlgebraicVector5 fwdPars = theForwardState.localParameters().vector();
0083     AlgebraicSymMatrix55 fwdCov = theForwardState.localError().matrix();
0084     computeQpMode(theForwardState, fwdPars, fwdCov);
0085     AlgebraicVector5 bwdPars = theBackwardState.localParameters().vector();
0086     AlgebraicSymMatrix55 bwdCov = theBackwardState.localError().matrix();
0087     computeQpMode(theBackwardState, bwdPars, bwdCov);
0088     LocalPoint hitPos(0., 0., 0.);
0089     LocalError hitErr(-1., -1., -1.);
0090     if (tm.recHit()->isValid()) {
0091       hitPos = tm.recHit()->localPosition();
0092       hitErr = tm.recHit()->localPositionError();
0093     }
0094
0095     CollinearFitAtTM collinearFit;
0096     CollinearFitAtTM::ResultVector fitParameters;
0097     CollinearFitAtTM::ResultMatrix fitCovariance;
0098     double fitChi2;
0099     bool CollFit = true;
0100     if (!collinearFit.fit(fwdPars, fwdCov, bwdPars, bwdCov, hitPos, hitErr, fitParameters, fitCovariance, fitChi2))
0101       CollFit = false;
0102
0103     if (CollFit) {
0104       double qpIn = fitParameters(0);
0105       double sig2In = fitCovariance(0, 0);
0106       double qpOut = fitParameters(1);
0107       double sig2Out = fitCovariance(1, 1);
0108       double corrInOut = fitCovariance(0, 1);
0109       double pIn = 1. / fabs(qpIn);
0110       double pOut = 1. / fabs(qpOut);
0111       double sig2DeltaP = pIn / qpIn * pIn / qpIn * sig2In - 2 * pIn / qpIn * pOut / qpOut * corrInOut +
0112                           pOut / qpOut * pOut / qpOut * sig2Out;
0113       //   std::cout << "fitted delta p = " << pOut-pIn << " sigma = "
0114       //        << sqrt(sig2DeltaP) << std::endl;
0115       dp = pOut - pIn;
0116       sigmaDp = sqrt(sig2DeltaP);
0117     }
0118   }
0119 }
0120
0121 PFGsfHelper::~PFGsfHelper() {}
0122 GlobalVector PFGsfHelper::computeP(bool ComputeMode) const {
0123   GlobalVector gsfp;
0124   if (ComputeMode)
0125     gsfp = GlobalVector(mode_Px, mode_Py, mode_Pz);
0126   else
0127     gsfp = theUpdateState.globalMomentum();
0128   return gsfp;
0129 }
0130 double PFGsfHelper::fittedDP() const { return dp; }
0131 double PFGsfHelper::sigmafittedDP() const { return sigmaDp; }
0132 bool PFGsfHelper::isValid() const { return Valid; }
0133
0134 void PFGsfHelper::computeQpMode(const TrajectoryStateOnSurface tsos,
0135                                 AlgebraicVector5& parameters,
0136                                 AlgebraicSymMatrix55& covariance) const {
0137   //
0138   // parameters and errors from combined state
0139   //
0140   parameters = tsos.localParameters().vector();
0141   covariance = tsos.localError().matrix();
0142   //
0143   // mode for parameter 0 (q/p)
0144   //
0145   MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(tsos, 0));
0146   GaussianSumUtilities1D qpGS(qpState);
0147   if (!qpGS.modeIsValid())
0148     return;
0149   double qp = qpGS.mode().mean();
0150   double varQp = qpGS.mode().variance();
0151   //
0152   // replace q/p value and variance, rescale correlation terms
0153   //   (heuristic procedure - alternative would be mode in 5D ...)
0154   //
0155   double VarQpRatio = sqrt(varQp / covariance(0, 0));
0156   parameters(0) = qp;
0157   covariance(0, 0) = varQp;
0158   for (int i = 1; i < 5; ++i)
0159     covariance(i, 0) *= VarQpRatio;
0160   //   std::cout << "covariance " << VarQpRatio << " "
0161   //        << covariance(1,0) << " " << covariance(0,1) << std::endl;
0162 }