Line Code
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
#include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
#include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/CommonDetUnit/interface/GeomDet.h"
#include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
#include "Geometry/CommonTopologies/interface/PixelTopology.h"
#include "Geometry/CommonTopologies/interface/StripTopology.h"
#include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
#include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
// #include "RecoTracker/MeasurementDet/interface/RecHitPropagator.h"

using namespace std;
TrajectoryAtInvalidHit::TrajectoryAtInvalidHit(const TrajectoryMeasurement& tm,
                                               const TrackerTopology* tTopo,
                                               const TrackerGeometry* tracker,
                                               const Propagator& propagator,
                                               const unsigned int mono) {
  if (tm.backwardPredictedState().isValid())
    theCombinedPredictedState =
        TrajectoryStateCombiner().combine(tm.forwardPredictedState(), tm.backwardPredictedState());
  else
    theCombinedPredictedState = tm.forwardPredictedState();

  if (!theCombinedPredictedState.isValid()) {
    return;
  }
  theHit = tm.recHit();
  iidd = theHit->geographicalId().rawId();
  StripSubdetector strip = StripSubdetector(iidd);
  unsigned int subid = strip.subdetId();
  // xB and yB are for absolute borders on the trajectories included in the study, sigmaX sigmaY are
  // significance cuts on the distance from the detector surface
  float xB = 0.;
  float sigmaX = 5.0;
  float yB = 0.;
  float sigmaY = 5.0;
  float sigmaYBond = 0.;
  //set bounds for point to be within to be counted in the study
  if (subid == StripSubdetector::TOB) {
    sigmaYBond = 5.0;
  }
  const GeomDetUnit* monodet;

  // if module is from a double sided layer, write out info for either the
  // rphi surface (mono = 1) or the stereo surface (mono = 2)--not the matched hit surface
  if ((mono > 0) && isDoubleSided(iidd, tTopo)) {
    // find matched det id, that is the matched hit surface between the two sensors
    unsigned int matched_iidd = iidd - (iidd & 0x3);
    DetId matched_id(matched_iidd);

    const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(tracker->idToDet(matched_id));

    // get the sensor det indicated by mono
    if (mono == 1)
      monodet = gdet->stereoDet();
    else
      monodet = gdet->monoDet();  // this should only be mono == 2

    // set theCombinedPredictedState to be on the sensor surface, not the matched surface
    DetId mono_id = monodet->geographicalId();
    const Surface& surface = tracker->idToDet(mono_id)->surface();
    theCombinedPredictedState = propagator.propagate(theCombinedPredictedState, surface);

    if (!theCombinedPredictedState.isValid()) {
      edm::LogWarning("TrajectoryAtInvalidHit") << "found invalid combinedpredictedstate after propagation" << endl;
      return;
    }

    //check if hit on trajectory is valid for this sensor
    if ((iidd & 0x3) == 0) {
      // passed id is on matched surface (ie neither rphi, nor stereo surface)
      // so the rechits on both sensors should be invalid
      hasValidHit = theHit->isValid();
      //cout << "should always be invalid = " << theHit->isValid() << endl;
    } else {
      if (iidd == mono_id) {
        //mono sensor for this TrajAtInvalidHit is same as originally passed in from trajectory
        //normally should only happen for valid hits, but sometimes is invalid
        hasValidHit = theHit->isValid();
        //cout << "check always true = " << theHit->isValid() << endl;
      } else {
        // mono sensor for this TrajAtInvalidHit is not as originally passed in from trajectory
        // so this is a partner missing from the trajectory and is thus invalid
        hasValidHit = false;
        //cout << "setting invalid by fiat" << endl;
      }
    }

    //set module id to be mono det
    iidd = monodet->geographicalId().rawId();
  } else {
    monodet = theHit->det();
    hasValidHit = theHit->isValid();
  }

  locX = theCombinedPredictedState.localPosition().x();
  locY = theCombinedPredictedState.localPosition().y();
  locZ = theCombinedPredictedState.localPosition().z();
  locXError = sqrt(theCombinedPredictedState.localError().positionError().xx());
  locYError = sqrt(theCombinedPredictedState.localError().positionError().yy());
  locDxDz = theCombinedPredictedState.localParameters().vector()[1];
  locDyDz = theCombinedPredictedState.localParameters().vector()[2];
  globX = theCombinedPredictedState.globalPosition().x();
  globY = theCombinedPredictedState.globalPosition().y();
  globZ = theCombinedPredictedState.globalPosition().z();

  // this should never be a glued det, only rphi or stero
  //cout << "From TrajAtValidHit module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;

  // Restrict the bound regions for better understanding of the modul assignment.

  LocalPoint BoundedPoint;
  float xx, yy, zz;

  // Insert the bounded values
  if (locX < 0.)
    xx = min(locX - xB, locX - sigmaX * locXError);
  else
    xx = max(locX + xB, locX + sigmaX * locXError);

  if (locY < 0.)
    yy = min(locY - yB, locY - sigmaY * locYError);
  else
    yy = max(locY + yB, locY + sigmaY * locYError);

  zz = theCombinedPredictedState.localPosition().z();

  BoundedPoint = LocalPoint(xx, yy, zz);

  if (monodet->surface().bounds().inside(BoundedPoint) && abs(locY) > sigmaYBond * locYError) {
    acceptance = true;
  } else {
    // hit is within xB, yB from the edge of the detector, so throw it out
    acceptance = false;
  }
}

double TrajectoryAtInvalidHit::localX() const { return locX; }
double TrajectoryAtInvalidHit::localY() const { return locY; }
double TrajectoryAtInvalidHit::localZ() const { return locZ; }
double TrajectoryAtInvalidHit::localErrorX() const { return locXError; }
double TrajectoryAtInvalidHit::localErrorY() const { return locYError; }
double TrajectoryAtInvalidHit::localDxDz() const { return locDxDz; }
double TrajectoryAtInvalidHit::localDyDz() const { return locDyDz; }
double TrajectoryAtInvalidHit::globalX() const { return globX; }
double TrajectoryAtInvalidHit::globalY() const { return globY; }
double TrajectoryAtInvalidHit::globalZ() const { return globZ; }

unsigned int TrajectoryAtInvalidHit::monodet_id() const { return iidd; }

bool TrajectoryAtInvalidHit::withinAcceptance() const { return acceptance; }

bool TrajectoryAtInvalidHit::validHit() const { return hasValidHit; }

bool TrajectoryAtInvalidHit::isDoubleSided(unsigned int iidd, const TrackerTopology* tTopo) const {
  StripSubdetector strip = StripSubdetector(iidd);
  unsigned int subid = strip.subdetId();
  unsigned int layer = 0;
  if (subid == StripSubdetector::TIB) {
    layer = tTopo->tibLayer(iidd);
    if (layer == 1 || layer == 2)
      return true;
    else
      return false;
  } else if (subid == StripSubdetector::TOB) {
    layer = tTopo->tobLayer(iidd) + 4;
    if (layer == 5 || layer == 6)
      return true;
    else
      return false;
  } else if (subid == StripSubdetector::TID) {
    layer = tTopo->tidRing(iidd) + 10;
    if (layer == 11 || layer == 12)
      return true;
    else
      return false;
  } else if (subid == StripSubdetector::TEC) {
    layer = tTopo->tecRing(iidd) + 13;
    if (layer == 14 || layer == 15 || layer == 18)
      return true;
    else
      return false;
  } else
    return false;
}

TrajectoryStateOnSurface TrajectoryAtInvalidHit::tsos() const { return theCombinedPredictedState; }