CosmicMuonProducer/src/CosmicMuonUtilities.cc

0001 /** \file CosmicMuonUtilities
0002  *
0003  *
0004  *  \author Chang Liu  -  Purdue University
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 // constructor
0021 //
0022 CosmicMuonUtilities::CosmicMuonUtilities() {}
0023
0024 //
0025 // destructor
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   // find first valid hit at both ends of the trajectory
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   // are there separate muon hits in 2 different hemispheres
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;  // need to optimize
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 }