GEMSegment/plugins/MuonSegFit.h

0001 #ifndef GEMSegment_MuonSegFit_h
0002 #define GEMSegment_MuonSegFit_h
0003
0004 // MuonSegFit.h - Segment fitting factored out of MuonRecHit segment builder based on
0005 // CSCSegFit.h  - Segment fitting factored out of CSC segment builders - Tim Cox
0006 // Last mod: 15.04.2016 Jason Lee
0007
0008 /* This as an object which is initialized by a set of rechits (2 to 6) in a
0009  * specific MuonRecHit and has the functionality to make a least squares fit
0010  * to a straight line in 2-dim for those rechits.
0011  * The covariance matrix and chi2 of the fit are calculated.
0012  * The original code made use of CLHEP matrices but this version uses
0013  * ROOT SMatrices because they are  multithreading compatible.
0014  * Because of this, the no. of rechits that can be handled is limited to
0015  * a maximum of 6, one per layer of a Muon chamber. This means maximum dimensions
0016  * can be specified at compile time and hence satisfies SMatrix constraints.
0017  * This means that if at a later stage we would change the geometry to have
0018  * for instance 10 detection layers, we will have to modify this code too.
0019  * For 2 hits of course there is no fit - just draw a straight line between them.
0020  * Details of the algorithm are in the .cc file
0021  *
0022  */
0023
0024 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0025
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027
0028 #include <Math/Functions.h>
0029 #include <Math/SVector.h>
0030 #include <Math/SMatrix.h>
0031
0032 #include <vector>
0033
0034 class MuonSegFit {
0035 public:
0036   // TYPES
0037   //typedef std::vector<const TrackingRecHit*> MuonRecHitContainer;
0038   typedef std::shared_ptr<TrackingRecHit> MuonRecHitPtr;
0039   typedef std::vector<MuonRecHitPtr> MuonRecHitContainer;
0040
0041   static const int MaxHits2 = 22;  // maxHit2 = 2*max hits
0042   // 12 x12 Symmetric
0043   typedef ROOT::Math::SMatrix<double, MaxHits2, MaxHits2, ROOT::Math::MatRepSym<double, MaxHits2> > SMatrixSym12;
0044
0045   // 12 x 4
0046   typedef ROOT::Math::SMatrix<double, MaxHits2, 4> SMatrix12by4;
0047
0048   // 4 x 4 General + Symmetric
0049   typedef ROOT::Math::SMatrix<double, 4> SMatrix4;
0050   typedef ROOT::Math::SMatrix<double, 4, 4, ROOT::Math::MatRepSym<double, 4> > SMatrixSym4;
0051
0052   // 2 x 2 Symmetric
0053   typedef ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > SMatrixSym2;
0054
0055   // 4-dim vector
0056   typedef ROOT::Math::SVector<double, 4> SVector4;
0057
0058   // PUBLIC FUNCTIONS
0059
0060   //@@ WANT OBJECT TO CACHE THE SET OF HITS SO CANNOT PASS BY REF
0061   // NOTE - We need local position of a RecHit w.r.t. the CHAMBER
0062   MuonSegFit(MuonRecHitContainer hits)
0063       : hits_(hits), uslope_(.0), vslope_(.0), chi2_(.0), ndof_(0), scaleXError_(1.0), fitdone_(false) {}
0064
0065   virtual ~MuonSegFit() {}
0066
0067   // Least-squares fit
0068   bool fit(void);  // fill uslope_, vslope_, intercept_  @@ FKA fitSlopes()
0069   // Calculate covariance matrix of fitted parameters
0070   AlgebraicSymMatrix covarianceMatrix(void);
0071
0072   // Change scale factor of rechit x error
0073   // - expert use only!
0074   void setScaleXError(double factor) { scaleXError_ = factor; }
0075
0076   // Fit values
0077   float xfit(float z) const;
0078   float yfit(float z) const;
0079
0080   // Deviations from fit for given input (local w.r.t. chamber)
0081   float xdev(float x, float z) const;
0082   float ydev(float y, float z) const;
0083   float Rdev(float x, float y, float z) const;
0084
0085   // Other public functions are accessors
0086   MuonRecHitContainer hits(void) const { return hits_; }
0087   double scaleXError(void) const { return scaleXError_; }
0088   size_t nhits(void) const { return hits_.size(); }
0089   double chi2(void) const { return chi2_; }
0090   int ndof(void) const { return ndof_; }
0091   LocalPoint intercept() const { return intercept_; }
0092   LocalVector localdir() const { return localdir_; }
0093   bool fitdone() const { return fitdone_; }
0094
0095 private:
0096   // PRIVATE FUNCTIONS
0097
0098   void fit2(void);     // fit for 2 hits
0099   void fitlsq(void);   // least-squares fit for 3-6 hits
0100   void setChi2(void);  // fill chi2_ & ndof_ @@ FKA fillChiSquared()
0101
0102 protected:
0103   // PROTECTED FUNCTIONS - derived class needs access
0104
0105   // Set segment direction 'out' from IP
0106   void setOutFromIP(void);  // fill localdir_  @@ FKA fillLocalDirection()
0107
0108   SMatrix12by4 derivativeMatrix(void);
0109   SMatrixSym12 weightMatrix(void);
0110   AlgebraicSymMatrix flipErrors(const SMatrixSym4&);
0111
0112   // PROTECTED MEMBER VARIABLES - derived class needs access
0113
0114   MuonRecHitContainer hits_;  //@@ FKA protoSegment
0115   float uslope_;              //@@ FKA protoSlope_u
0116   float vslope_;              //@@ FKA protoSlope_v
0117   LocalPoint intercept_;      //@@ FKA protoIntercept
0118   LocalVector localdir_;      //@@ FKA protoDirection
0119   double chi2_;               //@@ FKA protoChi2
0120   int ndof_;                  //@@ FKA protoNDF, which was double!!
0121   double scaleXError_;
0122   bool fitdone_;
0123 };
0124
0125 #endif