TrackingTools/src/MuonSeedFromRecHits.cc

0001 /**
0002  *  See header file for a description of this class.
0003  *
0004  *
0005  *  \author A. Vitelli - INFN Torino, V.Palichik
0006  *  \author porting  R. Bellan
0007  *
0008  */
0009 #include "RecoMuon/TrackingTools/interface/MuonSeedFromRecHits.h"
0010
0011 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0012
0013 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0014
0015 #include "MagneticField/Engine/interface/MagneticField.h"
0016 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0017
0018 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0019 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0020
0021 #include "DataFormats/TrajectoryState/interface/PTrajectoryStateOnDet.h"
0022 #include "DataFormats/Common/interface/OwnVector.h"
0023 #include "FWCore/Framework/interface/EventSetup.h"
0024 #include "FWCore/Framework/interface/ESHandle.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026
0027 #include "gsl/gsl_statistics.h"
0028
0029 using namespace std;
0030
0031 template <class T>
0032 T sqr(const T& t) {
0033   return t * t;
0034 }
0035
0036 MuonSeedFromRecHits::MuonSeedFromRecHits() : theField(nullptr) {}
0037
0038 TrajectorySeed MuonSeedFromRecHits::createSeed(float ptmean, float sptmean, ConstMuonRecHitPointer last) const {
0039   const std::string metname = "Muon|RecoMuon|MuonSeedFromRecHits";
0040
0041   MuonPatternRecoDumper debug;
0042
0043   // FIXME: put it into a parameter set!
0044   double theMinMomentum = 3.0;
0045   int charge = std::copysign(1, ptmean);
0046
0047   // Minimal pt
0048   if (fabs(ptmean) < theMinMomentum)
0049     ptmean = theMinMomentum * charge;
0050
0051   AlgebraicVector t(4);
0052   AlgebraicSymMatrix mat(5, 0);
0053
0054   // Fill the LocalTrajectoryParameters
0055   LocalPoint segPos = last->localPosition();
0056   GlobalVector mom = last->globalPosition() - GlobalPoint();
0057   GlobalVector polar(GlobalVector::Spherical(mom.theta(), last->globalDirection().phi(), 1.));
0058   polar *= fabs(ptmean) / polar.perp();
0059   LocalVector segDirFromPos = last->det()->toLocal(polar);
0060
0061   LocalTrajectoryParameters param(segPos, segDirFromPos, charge);
0062
0063   // this perform H.T() * parErr * H, which is the projection of the
0064   // the measurement error (rechit rf) to the state error (TSOS rf)
0065   // Legenda:
0066   // H => is the 4x5 projection matrix
0067   // parError the 4x4 parameter error matrix of the RecHit
0068
0069   // LogTrace(metname) << "Projection matrix:\n" << last->projectionMatrix();
0070   // LogTrace(metname) << "Error matrix:\n" << last->parametersError();
0071
0072   mat = last->parametersError().similarityT(last->projectionMatrix());
0073
0074   float p_err = sqr(sptmean / (ptmean * ptmean));
0075   mat[0][0] = p_err;
0076
0077   LocalTrajectoryError error(asSMatrix<5>(mat));
0078
0079   // Create the TrajectoryStateOnSurface
0080   TrajectoryStateOnSurface tsos(param, error, last->det()->surface(), theField);
0081
0082   // The following LogTraces must be moved somewhere else (StandAloneTrajectoryBuilder)
0083   // Here the TSOS does not have the magnetic field set, so dumpTSOS causes a crash
0084   // (when LogTrace/LogDebug is activated)
0085   //LogTrace(metname) << "Trajectory State on Surface before the extrapolation"<<endl;
0086   //LogTrace(metname) << debug.dumpTSOS(tsos);
0087
0088   // Take the DetLayer on which relies the rechit
0089   DetId id = last->geographicalId();
0090   // Segment layer
0091   //LogTrace(metname) << "The RecSegment relies on: "<<endl;
0092   //LogTrace(metname) << debug.dumpMuonId(id);
0093   //LogTrace(metname) << debug.dumpTSOS(tsos);
0094
0095   // Transform it in a TrajectoryStateOnSurface
0096
0097   PTrajectoryStateOnDet const& seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0098
0099   edm::OwnVector<TrackingRecHit> container;
0100   for (unsigned l = 0; l < theRhits.size(); l++) {
0101     container.push_back(theRhits[l]->hit()->clone());
0102   }
0103
0104   TrajectorySeed theSeed(seedTSOS, container, alongMomentum);
0105
0106   return theSeed;
0107 }