GEMCSCSegment/plugins/GEMCSCSegFit.h

0001 #ifndef GEMCSCSegment_GEMCSCSegFit_h
0002 #define GEMCSCSegment_GEMCSCSegFit_h
0003
0004 // GEMCSCSegFit.h extension to fit both CSCHits and GEMHits - Piet Verwilligen
0005 // Created: 21.04.2015
0006
0007 /* Work with general Tracking Rechits, such that GEM and CSC rechits
0008  * can be treated at the same level. In case needed the Tracking Rechit can be
0009  * cast to the CSCRecHit2D or the GEMRecHit for access to more specific info.
0010  */
0011
0012 // CSCSegFit.h - Segment fitting factored out of CSC segment builders - Tim Cox
0013 // Last mod: 03.02.2015
0014
0015 /* This as an object which is initialized by a set of rechits (2 to 6) in a
0016  * specific CSC and has the functionality to make a least squares fit to a
0017  * straight line in 2-dim for those rechits.
0018  * The covariance matrix and chi2 of the fit are calculated.
0019  * The original code made use of CLHEP matrices but this version uses
0020  * ROOT SMatrices because they are  multithreading compatible.
0021  * Because of this, the no. of rechits that can be handled is limited to
0022  * a maximum of 6, one per layer of a CSC. This means maximum dimensions
0023  * can be specified at compile time and hence satisfies SMatrix constraints.
0024  * For 2 hits of course there is no fit - just draw a straight line between them.
0025  * Details of the algorithm are in the .cc file
0026  *
0027  */
0028
0029 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0030 #include <DataFormats/CSCRecHit/interface/CSCRecHit2D.h>
0031 #include <DataFormats/GEMRecHit/interface/GEMRecHit.h>
0032 #include <DataFormats/MuonDetId/interface/CSCDetId.h>
0033
0034 #include <Geometry/CSCGeometry/interface/CSCChamber.h>
0035 #include <Geometry/CSCGeometry/interface/CSCLayer.h>
0036 #include <Geometry/GEMGeometry/interface/GEMChamber.h>
0037
0038 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0039 #include "FWCore/Utilities/interface/Exception.h"
0040
0041 #include <Math/Functions.h>
0042 #include <Math/SVector.h>
0043 #include <Math/SMatrix.h>
0044
0045 #include <vector>
0046
0047 class GEMCSCSegFit {
0048 public:
0049   // TYPES
0050
0051   // 16 x 16 Symmetric
0052   typedef ROOT::Math::SMatrix<double, 16, 16, ROOT::Math::MatRepSym<double, 16> > SMatrixSym16;
0053
0054   // 16 x 4
0055   typedef ROOT::Math::SMatrix<double, 16, 4> SMatrix16by4;
0056
0057   // 4 x 4 General + Symmetric
0058   typedef ROOT::Math::SMatrix<double, 4> SMatrix4;
0059   typedef ROOT::Math::SMatrix<double, 4, 4, ROOT::Math::MatRepSym<double, 4> > SMatrixSym4;
0060
0061   // 2 x 2 Symmetric
0062   typedef ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > SMatrixSym2;
0063
0064   // 4-dim vector
0065   typedef ROOT::Math::SVector<double, 4> SVector4;
0066
0067   // PUBLIC FUNCTIONS
0068
0069   //@@ WANT OBJECT TO CACHE THE SET OF HITS SO CANNOT PASS BY REF
0070   GEMCSCSegFit(std::map<uint32_t, const CSCLayer*> csclayermap,
0071                std::map<uint32_t, const GEMEtaPartition*> gemrollmap,
0072                const std::vector<const TrackingRecHit*> hits)
0073       : csclayermap_(csclayermap),
0074         gemetapartmap_(gemrollmap),
0075         hits_(hits),
0076         scaleXError_(1.0),
0077         refid_(csclayermap_.begin()->first),
0078         fitdone_(false) {
0079     // --- LogDebug info about reading of CSC Layer map and GEM Eta Partition map -----------------------
0080     edm::LogVerbatim("GEMCSCSegFit") << "[GEMCSCSegFit::ctor] cached the csclayermap and the gemrollmap";
0081
0082     // --- LogDebug for CSC Layer map -------------------------------------------------------------------
0083     std::stringstream csclayermapss;
0084     csclayermapss << "[GEMCSCSegFit::ctor] :: csclayermap :: elements [" << std::endl;
0085     for (std::map<uint32_t, const CSCLayer*>::const_iterator mapIt = csclayermap_.begin(); mapIt != csclayermap_.end();
0086          ++mapIt) {
0087       csclayermapss << "[CSC DetId " << mapIt->first << " =" << CSCDetId(mapIt->first) << ", CSC Layer "
0088                     << mapIt->second << " =" << (mapIt->second)->id() << "]," << std::endl;
0089     }
0090     csclayermapss << "]" << std::endl;
0091     std::string csclayermapstr = csclayermapss.str();
0092     edm::LogVerbatim("GEMCSCSegFit") << csclayermapstr;
0093     // --- End LogDebug -----------------------------------------------------------------------------------
0094
0095     // --- LogDebug for GEM Eta Partition map ------------------------------------------------------------
0096     std::stringstream gemetapartmapss;
0097     gemetapartmapss << "[GEMCSCSegFit::ctor] :: gemetapartmap :: elements [" << std::endl;
0098     for (std::map<uint32_t, const GEMEtaPartition*>::const_iterator mapIt = gemetapartmap_.begin();
0099          mapIt != gemetapartmap_.end();
0100          ++mapIt) {
0101       gemetapartmapss << "[GEM DetId " << mapIt->first << " =" << GEMDetId(mapIt->first) << ", GEM EtaPart "
0102                       << mapIt->second << "]," << std::endl;
0103     }
0104     gemetapartmapss << "]" << std::endl;
0105     std::string gemetapartmapstr = gemetapartmapss.str();
0106     edm::LogVerbatim("GEMCSCSegFit") << gemetapartmapstr;
0107     // --- End LogDebug -----------------------------------------------------------------------------------
0108   }
0109
0110   virtual ~GEMCSCSegFit() {}
0111
0112   // Least-squares fit
0113   void fit(void);  // fill uslope_, vslope_, intercept_  @@ FKA fitSlopes()
0114   // Calculate covariance matrix of fitted parameters
0115   AlgebraicSymMatrix covarianceMatrix(void);
0116
0117   // Change scale factor of rechit x error
0118   // - expert use only!
0119   void setScaleXError(double factor) { scaleXError_ = factor; }
0120
0121   // Fit values
0122   float xfit(float z) const;
0123   float yfit(float z) const;
0124
0125   // Deviations from fit for given input (local w.r.t. chamber)
0126   float xdev(float x, float z) const;
0127   float ydev(float y, float z) const;
0128   float Rdev(float x, float y, float z) const;
0129
0130   // Other public functions are accessors
0131   std::vector<const TrackingRecHit*> hits(void) const { return hits_; }
0132   double scaleXError(void) const { return scaleXError_; }
0133   size_t nhits(void) const { return hits_.size(); }
0134   double chi2(void) const { return chi2_; }
0135   int ndof(void) const { return ndof_; }
0136   LocalPoint intercept() const { return intercept_; }
0137   LocalVector localdir() const { return localdir_; }
0138
0139   const CSCChamber* cscchamber(uint32_t id) const {
0140     // Test whether id is found
0141     if (csclayermap_.find(id) == csclayermap_.end()) {  // id is not found
0142       throw cms::Exception("InvalidDetId") << "[GEMCSCSegFit] Failed to find CSCChamber in CSCLayerMap" << std::endl;
0143     }       // chamber is not found and exception is thrown
0144     else {  // id is found
0145       return (csclayermap_.find(id)->second)->chamber();
0146     }  // chamber found and returned
0147   }
0148   const CSCLayer* csclayer(uint32_t id) const {
0149     if (csclayermap_.find(id) == csclayermap_.end()) {
0150       throw cms::Exception("InvalidDetId") << "[GEMCSCSegFit] Failed to find CSCLayer in CSCLayerMap" << std::endl;
0151     } else {
0152       return csclayermap_.find(id)->second;
0153     }
0154   }
0155   const GEMEtaPartition* gemetapartition(uint32_t id) const {
0156     if (gemetapartmap_.find(id) == gemetapartmap_.end()) {
0157       throw cms::Exception("InvalidDetId")
0158           << "[GEMCSCSegFit] Failed to find GEMEtaPartition in GEMEtaPartMap" << std::endl;
0159     } else {
0160       return gemetapartmap_.find(id)->second;
0161     }
0162   }
0163   const CSCChamber* refcscchamber() const {
0164     if (csclayermap_.find(refid_) == csclayermap_.end()) {
0165       throw cms::Exception("InvalidDetId")
0166           << "[GEMCSCSegFit] Failed to find Reference CSCChamber in CSCLayerMap" << std::endl;
0167     } else {
0168       return (csclayermap_.find(refid_)->second)->chamber();
0169     }
0170   }
0171   bool fitdone() const { return fitdone_; }
0172
0173 private:
0174   // PRIVATE FUNCTIONS
0175
0176   void fit2(void);     // fit for 2 hits
0177   void fitlsq(void);   // least-squares fit for 3-6 hits
0178   void setChi2(void);  // fill chi2_ & ndof_ @@ FKA fillChiSquared()
0179
0180 protected:
0181   // PROTECTED FUNCTIONS - derived class needs access
0182
0183   // Set segment direction 'out' from IP
0184   void setOutFromIP(void);  // fill localdir_  @@ FKA fillLocalDirection()
0185
0186   SMatrix16by4 derivativeMatrix(void);
0187   SMatrixSym16 weightMatrix(void);
0188   AlgebraicSymMatrix flipErrors(const SMatrixSym4&);
0189
0190   // PROTECTED MEMBER VARIABLES - derived class needs access
0191
0192   std::map<uint32_t, const CSCLayer*> csclayermap_;
0193   std::map<uint32_t, const GEMEtaPartition*> gemetapartmap_;
0194   const CSCChamber* refcscchamber_;
0195
0196   std::vector<const TrackingRecHit*> hits_;  //@@ FKA protoSegment
0197   float uslope_;                             //@@ FKA protoSlope_u
0198   float vslope_;                             //@@ FKA protoSlope_v
0199   LocalPoint intercept_;                     //@@ FKA protoIntercept
0200   LocalVector localdir_;                     //@@ FKA protoDirection
0201   double chi2_;                              //@@ FKA protoChi2
0202   int ndof_;                                 //@@ FKA protoNDF, which was double!!
0203   double scaleXError_;
0204   uint32_t refid_;
0205   bool fitdone_;
0206   // order needs to be the same here when filled by the constructor
0207 };
0208
0209 #endif