PixelVertexFinding/interface/Cluster1DMerger.h

0001 #ifndef RecoVertex_PixelVertexFinding_Cluster1DMerger_h
0002 #define RecoVertex_PixelVertexFinding_Cluster1DMerger_h
0003
0004 #include "CommonTools/Clustering1D/interface/Cluster1D.h"
0005 #include "CommonTools/Clustering1D/interface/WeightEstimator.h"
0006
0007 #include "DataFormats/Math/interface/Point3D.h"
0008
0009 /**
0010  *  The class that should always be used to merge
0011  *  two Cluster1D into a single Cluster1D.
0012  */
0013
0014 namespace pixeltemp {
0015
0016   template <class T>
0017   class Cluster1DMerger {
0018   public:
0019     Cluster1DMerger(const WeightEstimator<T> &, const math::XYZPoint &bs = math::XYZPoint(0., 0., 0.));
0020     void setBeamSpot(const math::XYZPoint &bs) { theBS = bs; }
0021     ~Cluster1DMerger();
0022     Cluster1DMerger(const Cluster1DMerger &);
0023     Cluster1D<T> operator()(const Cluster1D<T> &first, const Cluster1D<T> &second) const;
0024
0025   private:
0026     WeightEstimator<T> *theEstimator;
0027     math::XYZPoint theBS;
0028   };
0029
0030   /*
0031  *                                implementation
0032  */
0033
0034   template <class T>
0035   Cluster1DMerger<T>::Cluster1DMerger(const WeightEstimator<T> &est, const math::XYZPoint &bs)
0036       : theEstimator(est.clone()), theBS(bs) {}
0037
0038   template <class T>
0039   Cluster1DMerger<T>::~Cluster1DMerger() {
0040     delete theEstimator;
0041   }
0042
0043   template <class T>
0044   Cluster1DMerger<T>::Cluster1DMerger(const Cluster1DMerger &other)
0045       : theEstimator(other.theEstimator->clone()), theBS(other.theBS) {}
0046
0047   template <class T>
0048   Cluster1D<T> Cluster1DMerger<T>::operator()(const Cluster1D<T> &first, const Cluster1D<T> &second) const {
0049     std::vector<const T *> tracks = first.tracks();
0050     std::vector<const T *> sectracks = second.tracks();
0051     for (typename std::vector<const T *>::const_iterator i = sectracks.begin(); i != sectracks.end(); ++i) {
0052       tracks.push_back(*i);
0053     };
0054     float newpos =
0055         (first.position().value() * first.weight() / first.position().error() / first.position().error() +
0056          second.position().value() * second.weight() / second.position().error() / second.position().error()) /
0057         (first.weight() / first.position().error() / first.position().error() +
0058          second.weight() / second.position().error() / second.position().error());
0059
0060     float newerr = sqrt(first.position().error() * first.position().error() +
0061                         second.position().error() * second.position().error());
0062     float newWeight = theEstimator->weight(tracks);
0063
0064     // Now try a different method
0065     float sumUp = 0;
0066     float sumDown = 0;
0067     float err = 0;
0068     for (unsigned int i = 0; i < tracks.size(); ++i) {
0069       //      float err2 = tracks[i]->covariance( reco::TrackBase::i_dz, reco::TrackBase::i_dz );
0070       float err2 = tracks[i]->dzError();
0071       err2 *= err2;
0072
0073       if (err2 != 0) {
0074         sumUp += tracks[i]->dz(theBS) * 1 / err2;  // error-weighted average of Z at IP
0075         sumDown += 1 / err2;
0076       }
0077       err += sqrt(err2);
0078     }
0079
0080     newpos = sumUp / sumDown;
0081     newerr = err / tracks.size() / sqrt(1.0 * tracks.size());
0082
0083     Measurement1D newmeas(newpos, newerr);
0084     return Cluster1D<T>(newmeas, tracks, newWeight);
0085   }
0086 }  // namespace pixeltemp
0087 #endif