CosmicMuonUtilities.cc

Go to the documentation of this file.

#include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonUtilities.h"

#include "TrackingTools/GeomPropagators/interface/Propagator.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
#include "DataFormats/Math/interface/deltaPhi.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace edm;
using namespace std;

//
// constructor
//
CosmicMuonUtilities::CosmicMuonUtilities() {}

//
// destructor
//
CosmicMuonUtilities::~CosmicMuonUtilities() {}

void CosmicMuonUtilities::reverseDirection(TrajectoryStateOnSurface& tsos, const MagneticField* mgfield) const {
  GlobalTrajectoryParameters gtp(tsos.globalPosition(), -tsos.globalMomentum(), -tsos.charge(), mgfield);
  TrajectoryStateOnSurface newTsos(gtp, tsos.cartesianError(), tsos.surface());
  tsos = newTsos;

  return;
}

//
//
//
string CosmicMuonUtilities::print(const MuonTransientTrackingRecHit::ConstMuonRecHitContainer& hits) const {
  stringstream output;

  for (ConstMuonRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
    if (!(*ir)->isValid()) {
      output << "invalid RecHit" << endl;
      continue;
    }

    const GlobalPoint& pos = (*ir)->globalPosition();
    output << "pos" << pos << "radius " << pos.perp() << "  dim " << (*ir)->dimension() << "  det "
           << (*ir)->det()->geographicalId().det() << "  sub det " << (*ir)->det()->subDetector() << endl;
  }

  return output.str();
}

//
//
//
string CosmicMuonUtilities::print(const MuonTransientTrackingRecHit::MuonRecHitContainer& hits) const {
  stringstream output;

  for (MuonRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
    if (!(*ir)->isValid()) {
      output << "invalid RecHit" << endl;
      continue;
    }

    const GlobalPoint& pos = (*ir)->globalPosition();
    output << "pos" << pos << "radius " << pos.perp() << "  dim " << (*ir)->dimension() << "  det "
           << (*ir)->det()->geographicalId().det() << "  sub det " << (*ir)->det()->subDetector() << endl;
  }

  return output.str();
}

//
//
//
string CosmicMuonUtilities::print(const TransientTrackingRecHit::ConstRecHitContainer& hits) const {
  stringstream output;

  for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
    if (!(*ir)->isValid()) {
      output << "invalid RecHit" << endl;
      continue;
    }

    const GlobalPoint& pos = (*ir)->globalPosition();
    output << "pos" << pos << "radius " << pos.perp() << "  dim " << (*ir)->dimension() << "  det "
           << (*ir)->det()->geographicalId().det() << "  sub det " << (*ir)->det()->subDetector() << endl;
  }

  return output.str();
}

//
//
//
bool CosmicMuonUtilities::isTraversing(const Trajectory& t) const {
  TransientTrackingRecHit::ConstRecHitContainer hits = t.recHits();

  if (hits.empty())
    return false;

  ConstRecHitContainer::const_iterator frontHit = hits.begin();
  ConstRecHitContainer::const_iterator backHit = hits.end() - 1;

  // find first valid hit at both ends of the trajectory
  while (!(*frontHit)->isValid() && frontHit != backHit) {
    ++frontHit;
  }
  while (!(*backHit)->isValid() && backHit != frontHit) {
    --backHit;
  }

  if (frontHit == backHit)
    return false;

  GlobalPoint frontPos = (*frontHit)->globalPosition();
  GlobalPoint backPos = (*backHit)->globalPosition();

  // are there separate muon hits in 2 different hemispheres
  GlobalVector deltaPos(frontPos - backPos);
  if (deltaPos.mag() < 100.0)
    return false;
  if (fabs(deltaPos.z()) > 500.0)
    return true;
  if (deltaPos.perp() > 350.0)
    return true;
  GlobalPoint middle(
      (frontPos.x() + backPos.x()) / 2, (frontPos.y() + backPos.y()) / 2, (frontPos.z() + backPos.z()) / 2);

  return ((middle.perp() < frontPos.perp()) && (middle.perp() < backPos.perp()));
}

//
//
//
TrajectoryStateOnSurface CosmicMuonUtilities::stepPropagate(const TrajectoryStateOnSurface& tsos,
                                                            const ConstRecHitPointer& hit,
                                                            const Propagator& prop) const {
  const std::string metname = "Muon|RecoMuon|CosmicMuonUtilities";

  GlobalPoint start = tsos.globalPosition();
  GlobalPoint dest = hit->globalPosition();
  GlobalVector StepVector = dest - start;
  GlobalVector UnitStepVector = StepVector.unit();
  GlobalPoint GP = start;
  TrajectoryStateOnSurface currTsos(tsos);
  TrajectoryStateOnSurface predTsos;
  float totalDis = StepVector.mag();
  LogTrace(metname) << "stepPropagate: propagate from: " << start << " to " << dest;
  LogTrace(metname) << "stepPropagate: their distance: " << totalDis;

  int steps = 3;  // need to optimize

  float oneStep = totalDis / steps;
  Basic3DVector<float> Basic3DV(StepVector.x(), StepVector.y(), StepVector.z());
  for (int istep = 0; istep < steps - 1; istep++) {
    GP += oneStep * UnitStepVector;
    Surface::PositionType pos(GP.x(), GP.y(), GP.z());
    LogTrace(metname) << "stepPropagate: a middle plane: " << pos << endl;
    Surface::RotationType rot(Basic3DV, float(0));
    PlaneBuilder::ReturnType SteppingPlane = PlaneBuilder().plane(pos, rot);
    TrajectoryStateOnSurface predTsos = prop.propagate(currTsos, *SteppingPlane);
    if (predTsos.isValid()) {
      currTsos = predTsos;
      LogTrace(metname) << "stepPropagate: middle state " << currTsos.globalPosition() << endl;
    }
  }

  predTsos = prop.propagate(currTsos, hit->det()->surface());

  return predTsos;
}