File indexing completed on 2024-04-06 12:26:52
0001
0002
0003
0004
0005
0006
0007 #include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonUtilities.h"
0008
0009 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0010 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0011 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 using namespace edm;
0017 using namespace std;
0018
0019
0020
0021
0022 CosmicMuonUtilities::CosmicMuonUtilities() {}
0023
0024
0025
0026
0027 CosmicMuonUtilities::~CosmicMuonUtilities() {}
0028
0029 void CosmicMuonUtilities::reverseDirection(TrajectoryStateOnSurface& tsos, const MagneticField* mgfield) const {
0030 GlobalTrajectoryParameters gtp(tsos.globalPosition(), -tsos.globalMomentum(), -tsos.charge(), mgfield);
0031 TrajectoryStateOnSurface newTsos(gtp, tsos.cartesianError(), tsos.surface());
0032 tsos = newTsos;
0033
0034 return;
0035 }
0036
0037
0038
0039
0040 string CosmicMuonUtilities::print(const MuonTransientTrackingRecHit::ConstMuonRecHitContainer& hits) const {
0041 stringstream output;
0042
0043 for (ConstMuonRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0044 if (!(*ir)->isValid()) {
0045 output << "invalid RecHit" << endl;
0046 continue;
0047 }
0048
0049 const GlobalPoint& pos = (*ir)->globalPosition();
0050 output << "pos" << pos << "radius " << pos.perp() << " dim " << (*ir)->dimension() << " det "
0051 << (*ir)->det()->geographicalId().det() << " sub det " << (*ir)->det()->subDetector() << endl;
0052 }
0053
0054 return output.str();
0055 }
0056
0057
0058
0059
0060 string CosmicMuonUtilities::print(const MuonTransientTrackingRecHit::MuonRecHitContainer& hits) const {
0061 stringstream output;
0062
0063 for (MuonRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0064 if (!(*ir)->isValid()) {
0065 output << "invalid RecHit" << endl;
0066 continue;
0067 }
0068
0069 const GlobalPoint& pos = (*ir)->globalPosition();
0070 output << "pos" << pos << "radius " << pos.perp() << " dim " << (*ir)->dimension() << " det "
0071 << (*ir)->det()->geographicalId().det() << " sub det " << (*ir)->det()->subDetector() << endl;
0072 }
0073
0074 return output.str();
0075 }
0076
0077
0078
0079
0080 string CosmicMuonUtilities::print(const TransientTrackingRecHit::ConstRecHitContainer& hits) const {
0081 stringstream output;
0082
0083 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0084 if (!(*ir)->isValid()) {
0085 output << "invalid RecHit" << endl;
0086 continue;
0087 }
0088
0089 const GlobalPoint& pos = (*ir)->globalPosition();
0090 output << "pos" << pos << "radius " << pos.perp() << " dim " << (*ir)->dimension() << " det "
0091 << (*ir)->det()->geographicalId().det() << " sub det " << (*ir)->det()->subDetector() << endl;
0092 }
0093
0094 return output.str();
0095 }
0096
0097
0098
0099
0100 bool CosmicMuonUtilities::isTraversing(const Trajectory& t) const {
0101 TransientTrackingRecHit::ConstRecHitContainer hits = t.recHits();
0102
0103 if (hits.empty())
0104 return false;
0105
0106 ConstRecHitContainer::const_iterator frontHit = hits.begin();
0107 ConstRecHitContainer::const_iterator backHit = hits.end() - 1;
0108
0109
0110 while (!(*frontHit)->isValid() && frontHit != backHit) {
0111 ++frontHit;
0112 }
0113 while (!(*backHit)->isValid() && backHit != frontHit) {
0114 --backHit;
0115 }
0116
0117 if (frontHit == backHit)
0118 return false;
0119
0120 GlobalPoint frontPos = (*frontHit)->globalPosition();
0121 GlobalPoint backPos = (*backHit)->globalPosition();
0122
0123
0124 GlobalVector deltaPos(frontPos - backPos);
0125 if (deltaPos.mag() < 100.0)
0126 return false;
0127 if (fabs(deltaPos.z()) > 500.0)
0128 return true;
0129 if (deltaPos.perp() > 350.0)
0130 return true;
0131 GlobalPoint middle(
0132 (frontPos.x() + backPos.x()) / 2, (frontPos.y() + backPos.y()) / 2, (frontPos.z() + backPos.z()) / 2);
0133
0134 return ((middle.perp() < frontPos.perp()) && (middle.perp() < backPos.perp()));
0135 }
0136
0137
0138
0139
0140 TrajectoryStateOnSurface CosmicMuonUtilities::stepPropagate(const TrajectoryStateOnSurface& tsos,
0141 const ConstRecHitPointer& hit,
0142 const Propagator& prop) const {
0143 const std::string metname = "Muon|RecoMuon|CosmicMuonUtilities";
0144
0145 GlobalPoint start = tsos.globalPosition();
0146 GlobalPoint dest = hit->globalPosition();
0147 GlobalVector StepVector = dest - start;
0148 GlobalVector UnitStepVector = StepVector.unit();
0149 GlobalPoint GP = start;
0150 TrajectoryStateOnSurface currTsos(tsos);
0151 TrajectoryStateOnSurface predTsos;
0152 float totalDis = StepVector.mag();
0153 LogTrace(metname) << "stepPropagate: propagate from: " << start << " to " << dest;
0154 LogTrace(metname) << "stepPropagate: their distance: " << totalDis;
0155
0156 int steps = 3;
0157
0158 float oneStep = totalDis / steps;
0159 Basic3DVector<float> Basic3DV(StepVector.x(), StepVector.y(), StepVector.z());
0160 for (int istep = 0; istep < steps - 1; istep++) {
0161 GP += oneStep * UnitStepVector;
0162 Surface::PositionType pos(GP.x(), GP.y(), GP.z());
0163 LogTrace(metname) << "stepPropagate: a middle plane: " << pos << endl;
0164 Surface::RotationType rot(Basic3DV, float(0));
0165 PlaneBuilder::ReturnType SteppingPlane = PlaneBuilder().plane(pos, rot);
0166 TrajectoryStateOnSurface predTsos = prop.propagate(currTsos, *SteppingPlane);
0167 if (predTsos.isValid()) {
0168 currTsos = predTsos;
0169 LogTrace(metname) << "stepPropagate: middle state " << currTsos.globalPosition() << endl;
0170 }
0171 }
0172
0173 predTsos = prop.propagate(currTsos, hit->det()->surface());
0174
0175 return predTsos;
0176 }