SiStripHitEfficiency/src/TrajectoryAtInvalidHit.cc

0001 #include "CalibTracker/SiStripHitEfficiency/interface/TrajectoryAtInvalidHit.h"
0002 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0003 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementVector.h"
0004 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0005 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0008 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0009 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0010 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0011 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
0012 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "TrackingTools/TransientTrackingRecHit/interface/TrackingRecHitProjector.h"
0016 // #include "RecoTracker/MeasurementDet/interface/RecHitPropagator.h"
0017
0018 using namespace std;
0019 TrajectoryAtInvalidHit::TrajectoryAtInvalidHit(const TrajectoryMeasurement& tm,
0020                                                const TrackerTopology* tTopo,
0021                                                const TrackerGeometry* tracker,
0022                                                const Propagator& propagator,
0023                                                const unsigned int mono) {
0024   if (tm.backwardPredictedState().isValid())
0025     theCombinedPredictedState =
0026         TrajectoryStateCombiner().combine(tm.forwardPredictedState(), tm.backwardPredictedState());
0027   else
0028     theCombinedPredictedState = tm.forwardPredictedState();
0029
0030   if (!theCombinedPredictedState.isValid()) {
0031     return;
0032   }
0033   theHit = tm.recHit();
0034   iidd = theHit->geographicalId().rawId();
0035   StripSubdetector strip = StripSubdetector(iidd);
0036   unsigned int subid = strip.subdetId();
0037   // xB and yB are for absolute borders on the trajectories included in the study, sigmaX sigmaY are
0038   // significance cuts on the distance from the detector surface
0039   float xB = 0.;
0040   float sigmaX = 5.0;
0041   float yB = 0.;
0042   float sigmaY = 5.0;
0043   float sigmaYBond = 0.;
0044   //set bounds for point to be within to be counted in the study
0045   if (subid == StripSubdetector::TOB) {
0046     sigmaYBond = 5.0;
0047   }
0048   const GeomDetUnit* monodet;
0049
0050   // if module is from a double sided layer, write out info for either the
0051   // rphi surface (mono = 1) or the stereo surface (mono = 2)--not the matched hit surface
0052   if ((mono > 0) && isDoubleSided(iidd, tTopo)) {
0053     // find matched det id, that is the matched hit surface between the two sensors
0054     unsigned int matched_iidd = iidd - (iidd & 0x3);
0055     DetId matched_id(matched_iidd);
0056
0057     const GluedGeomDet* gdet = static_cast<const GluedGeomDet*>(tracker->idToDet(matched_id));
0058
0059     // get the sensor det indicated by mono
0060     if (mono == 1)
0061       monodet = gdet->stereoDet();
0062     else
0063       monodet = gdet->monoDet();  // this should only be mono == 2
0064
0065     // set theCombinedPredictedState to be on the sensor surface, not the matched surface
0066     DetId mono_id = monodet->geographicalId();
0067     const Surface& surface = tracker->idToDet(mono_id)->surface();
0068     theCombinedPredictedState = propagator.propagate(theCombinedPredictedState, surface);
0069
0070     if (!theCombinedPredictedState.isValid()) {
0071       edm::LogWarning("TrajectoryAtInvalidHit") << "found invalid combinedpredictedstate after propagation" << endl;
0072       return;
0073     }
0074
0075     //check if hit on trajectory is valid for this sensor
0076     if ((iidd & 0x3) == 0) {
0077       // passed id is on matched surface (ie neither rphi, nor stereo surface)
0078       // so the rechits on both sensors should be invalid
0079       hasValidHit = theHit->isValid();
0080       //cout << "should always be invalid = " << theHit->isValid() << endl;
0081     } else {
0082       if (iidd == mono_id) {
0083         //mono sensor for this TrajAtInvalidHit is same as originally passed in from trajectory
0084         //normally should only happen for valid hits, but sometimes is invalid
0085         hasValidHit = theHit->isValid();
0086         //cout << "check always true = " << theHit->isValid() << endl;
0087       } else {
0088         // mono sensor for this TrajAtInvalidHit is not as originally passed in from trajectory
0089         // so this is a partner missing from the trajectory and is thus invalid
0090         hasValidHit = false;
0091         //cout << "setting invalid by fiat" << endl;
0092       }
0093     }
0094
0095     //set module id to be mono det
0096     iidd = monodet->geographicalId().rawId();
0097   } else {
0098     monodet = theHit->det();
0099     hasValidHit = theHit->isValid();
0100   }
0101
0102   locX = theCombinedPredictedState.localPosition().x();
0103   locY = theCombinedPredictedState.localPosition().y();
0104   locZ = theCombinedPredictedState.localPosition().z();
0105   locXError = sqrt(theCombinedPredictedState.localError().positionError().xx());
0106   locYError = sqrt(theCombinedPredictedState.localError().positionError().yy());
0107   locDxDz = theCombinedPredictedState.localParameters().vector()[1];
0108   locDyDz = theCombinedPredictedState.localParameters().vector()[2];
0109   globX = theCombinedPredictedState.globalPosition().x();
0110   globY = theCombinedPredictedState.globalPosition().y();
0111   globZ = theCombinedPredictedState.globalPosition().z();
0112
0113   // this should never be a glued det, only rphi or stero
0114   //cout << "From TrajAtValidHit module " << iidd << "   matched/stereo/rphi = " << ((iidd & 0x3)==0) << "/" << ((iidd & 0x3)==1) << "/" << ((iidd & 0x3)==2) << endl;
0115
0116   // Restrict the bound regions for better understanding of the modul assignment.
0117
0118   LocalPoint BoundedPoint;
0119   float xx, yy, zz;
0120
0121   // Insert the bounded values
0122   if (locX < 0.)
0123     xx = min(locX - xB, locX - sigmaX * locXError);
0124   else
0125     xx = max(locX + xB, locX + sigmaX * locXError);
0126
0127   if (locY < 0.)
0128     yy = min(locY - yB, locY - sigmaY * locYError);
0129   else
0130     yy = max(locY + yB, locY + sigmaY * locYError);
0131
0132   zz = theCombinedPredictedState.localPosition().z();
0133
0134   BoundedPoint = LocalPoint(xx, yy, zz);
0135
0136   if (monodet->surface().bounds().inside(BoundedPoint) && abs(locY) > sigmaYBond * locYError) {
0137     acceptance = true;
0138   } else {
0139     // hit is within xB, yB from the edge of the detector, so throw it out
0140     acceptance = false;
0141   }
0142 }
0143
0144 double TrajectoryAtInvalidHit::localX() const { return locX; }
0145 double TrajectoryAtInvalidHit::localY() const { return locY; }
0146 double TrajectoryAtInvalidHit::localZ() const { return locZ; }
0147 double TrajectoryAtInvalidHit::localErrorX() const { return locXError; }
0148 double TrajectoryAtInvalidHit::localErrorY() const { return locYError; }
0149 double TrajectoryAtInvalidHit::localDxDz() const { return locDxDz; }
0150 double TrajectoryAtInvalidHit::localDyDz() const { return locDyDz; }
0151 double TrajectoryAtInvalidHit::globalX() const { return globX; }
0152 double TrajectoryAtInvalidHit::globalY() const { return globY; }
0153 double TrajectoryAtInvalidHit::globalZ() const { return globZ; }
0154
0155 unsigned int TrajectoryAtInvalidHit::monodet_id() const { return iidd; }
0156
0157 bool TrajectoryAtInvalidHit::withinAcceptance() const { return acceptance; }
0158
0159 bool TrajectoryAtInvalidHit::validHit() const { return hasValidHit; }
0160
0161 bool TrajectoryAtInvalidHit::isDoubleSided(unsigned int iidd, const TrackerTopology* tTopo) const {
0162   StripSubdetector strip = StripSubdetector(iidd);
0163   unsigned int subid = strip.subdetId();
0164   unsigned int layer = 0;
0165   if (subid == StripSubdetector::TIB) {
0166     layer = tTopo->tibLayer(iidd);
0167     if (layer == 1 || layer == 2)
0168       return true;
0169     else
0170       return false;
0171   } else if (subid == StripSubdetector::TOB) {
0172     layer = tTopo->tobLayer(iidd) + 4;
0173     if (layer == 5 || layer == 6)
0174       return true;
0175     else
0176       return false;
0177   } else if (subid == StripSubdetector::TID) {
0178     layer = tTopo->tidRing(iidd) + 10;
0179     if (layer == 11 || layer == 12)
0180       return true;
0181     else
0182       return false;
0183   } else if (subid == StripSubdetector::TEC) {
0184     layer = tTopo->tecRing(iidd) + 13;
0185     if (layer == 14 || layer == 15 || layer == 18)
0186       return true;
0187     else
0188       return false;
0189   } else
0190     return false;
0191 }
0192
0193 TrajectoryStateOnSurface TrajectoryAtInvalidHit::tsos() const { return theCombinedPredictedState; }