MultiTrajectoryStateCombiner Class Reference

#include <MultiTrajectoryStateCombiner.h>

Public Member Functions
TrajectoryStateOnSurface	combine (const std::vector< TrajectoryStateOnSurface > &tsos) const
	MultiTrajectoryStateCombiner ()
	~MultiTrajectoryStateCombiner ()

Detailed Description

Class which combines a set of components of a Gaussian mixture into a single component. Given all the components of a mixture, it calculates the mean and covariance matrix of the entire mixture. This combiner class can also be used in the process of transforming a Gaussian mixture into another Gaussian mixture with a smaller number of components. The relevant formulas can be found in R. Fruhwirth, Computer Physics Communications 100 (1997), 1.

Definition at line 15 of file MultiTrajectoryStateCombiner.h.

Constructor & Destructor Documentation

MultiTrajectoryStateCombiner::MultiTrajectoryStateCombiner ( )

inline

Definition at line 19 of file MultiTrajectoryStateCombiner.h.

{}

MultiTrajectoryStateCombiner::~MultiTrajectoryStateCombiner ( )

inline

Definition at line 20 of file MultiTrajectoryStateCombiner.h.

{}

Member Function Documentation

TrajectoryStateOnSurface MultiTrajectoryStateCombiner::combine

(

const std::vector< TrajectoryStateOnSurface > &

tsos

)

const

Definition at line 6 of file MultiTrajectoryStateCombiner.cc.

References diffTreeTool::diff, plotscripts::mean(), asciidump::s, and CommonMethods::weight().

Referenced by BasicMultiTrajectoryState::checkCombinedState().

                                                                                           {
  
  if (tsos.empty()) {
    edm::LogError("MultiTrajectoryStateCombiner") 
      << "Trying to collapse empty set of trajectory states!";
    return TrajectoryStateOnSurface();
  }

  double pzSign = tsos.front().localParameters().pzSign();
  for (std::vector<TrajectoryStateOnSurface>::const_iterator it = tsos.begin(); 
       it != tsos.end(); it++) {
    if (it->localParameters().pzSign() != pzSign) {
      edm::LogError("MultiTrajectoryStateCombiner") 
    << "Trying to collapse trajectory states with different signs on p_z!";
      return TrajectoryStateOnSurface();
    }
  }
  
  if (tsos.size() == 1) {
    return TrajectoryStateOnSurface(tsos.front());
  }
  
  double sumw = 0.;
  //int dim = tsos.front().localParameters().vector().num_row();
  AlgebraicVector5 mean;
  AlgebraicSymMatrix55 covarPart1, covarPart2;
  for (std::vector<TrajectoryStateOnSurface>::const_iterator it1 = tsos.begin(); 
       it1 != tsos.end(); it1++) {
    double weight = it1->weight();
    AlgebraicVector5 param = it1->localParameters().vector();
    sumw += weight;
    mean += weight * param;
    covarPart1 += weight * it1->localError().matrix();
    for (std::vector<TrajectoryStateOnSurface>::const_iterator it2 = it1 + 1; 
     it2 != tsos.end(); it2++) {
      AlgebraicVector5 diff = param - it2->localParameters().vector();
      AlgebraicSymMatrix11 s = AlgebraicMatrixID(); //stupid trick to make CLHEP work decently
      covarPart2 += weight * it2->weight() * 
                    ROOT::Math::Similarity(AlgebraicMatrix51(diff.Array(), 5), s);
                        //FIXME: we can surely write this thing in a better way
    }   
  }
  double sumwI = 1.0/sumw;
  mean *= sumwI;
  covarPart1 *= sumwI; covarPart2 *= (sumwI*sumwI);
  AlgebraicSymMatrix55 covar = covarPart1 + covarPart2;

  return TrajectoryStateOnSurface(LocalTrajectoryParameters(mean, pzSign), 
                  LocalTrajectoryError(covar), 
                  tsos.front().surface(), 
                  &(tsos.front().globalParameters().magneticField()),
                  tsos.front().surfaceSide(), 
                  sumw);
}