OfflineValidation/plugins/OverlapValidation.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/OfflineValidation
0004 // Class:      OverlapValidation
0005 //
0006 /**\ Alignment/OfflineValidation/src/OverlapValidation.cc
0007
0008  Description: <one line class summary>
0009
0010 */
0011 //
0012 // Original Authors:  Wolfgang Adam, Keith Ulmer
0013 //         Created:  Thu Oct 11 14:53:32 CEST 2007
0014 // $Id: OverlapValidation.cc,v 1.16 2009/11/04 19:40:27 kaulmer Exp $
0015 //
0016 //
0017
0018 // system include files
0019 #include <memory>
0020
0021 // user include files
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0024
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/Utilities/interface/EDGetToken.h"
0030 #include "FWCore/Utilities/interface/ESGetToken.h"
0031
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0034 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0035 #include "MagneticField/Engine/interface/MagneticField.h"
0036 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0039
0040 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0041 #include "FWCore/ParameterSet/interface/FileInPath.h"
0042 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0043
0044 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0045 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0046 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0048 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0049
0050 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0053 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0054 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0055 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0056 #include "DataFormats/DetId/interface/DetId.h"
0057 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0058 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0059 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0060 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0061 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0062 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0063 #include "RecoTracker/TransientTrackingRecHit/interface/TSiStripMatchedRecHit.h"
0064 #include "RecoTracker/TransientTrackingRecHit/interface/ProjectedRecHit2D.h"
0065 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0066 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0067
0068 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0069 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0070 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0071 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0072 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0073 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0074 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0075 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0076 #include "Alignment/TrackerAlignment/interface/TrackerAlignableId.h"
0077 #include "TFile.h"
0078 #include "TTree.h"
0079
0080 #include <vector>
0081 #include <utility>
0082
0083 using namespace std;
0084 //
0085 // class decleration
0086 //
0087
0088 class OverlapValidation : public edm::one::EDAnalyzer<> {
0089 public:
0090   explicit OverlapValidation(const edm::ParameterSet&);
0091   ~OverlapValidation() override;
0092
0093 private:
0094   typedef vector<Trajectory> TrajectoryCollection;
0095
0096   edm::EDGetTokenT<TrajectoryCollection> trajectoryToken_;
0097   typedef TransientTrackingRecHit::ConstRecHitPointer ConstRecHitPointer;
0098   void analyze(const edm::Event&, const edm::EventSetup&) override;
0099   void endJob() override;
0100
0101   virtual void analyzeTrajectory(const Trajectory&,
0102                                  const Propagator&,
0103                                  TrackerHitAssociator&,
0104                                  const TrackerTopology* const tTopo);
0105   int layerFromId(const DetId&, const TrackerTopology* const tTopo) const;
0106
0107   // ----------member data ---------------------------
0108
0109   const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> geomToken_;
0110   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0111   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0112
0113   edm::ParameterSet config_;
0114   edm::InputTag trajectoryTag_;
0115   SiStripDetInfo detInfo_;
0116   bool doSimHit_;
0117   const TrackerGeometry* trackerGeometry_;
0118   const MagneticField* magField_;
0119
0120   TrajectoryStateCombiner combiner_;
0121   int overlapCounts_[3];
0122
0123   TTree* rootTree_;
0124   edm::FileInPath FileInPath_;
0125   const int compressionSettings_;
0126   const bool addExtraBranches_;
0127   const int minHitsCut_;
0128   const float chi2ProbCut_;
0129
0130   float overlapPath_;
0131   uint layer_;
0132   unsigned short hitCounts_[2];
0133   float chi2_[2];
0134   unsigned int overlapIds_[2];
0135   float predictedPositions_[3][2];
0136   float predictedLocalParameters_[5][2];
0137   float predictedLocalErrors_[5][2];
0138   float predictedDeltaXError_;
0139   float predictedDeltaYError_;
0140   char relativeXSign_;
0141   char relativeYSign_;
0142   float hitPositions_[2];
0143   float hitErrors_[2];
0144   float hitPositionsY_[2];
0145   float hitErrorsY_[2];
0146   float simHitPositions_[2];
0147   float simHitPositionsY_[2];
0148   float clusterWidthX_[2];
0149   float clusterWidthY_[2];
0150   float clusterSize_[2];
0151   uint clusterCharge_[2];
0152   int edge_[2];
0153
0154   vector<bool> acceptLayer;
0155   float momentum_;
0156   uint run_, event_;
0157   bool barrelOnly_;
0158
0159   //added by Heshy and Jared
0160   float moduleX_[2];
0161   float moduleY_[2];
0162   float moduleZ_[2];
0163   int subdetID;
0164   const Point2DBase<float, LocalTag> zerozero = Point2DBase<float, LocalTag>(0, 0);
0165   const Point2DBase<float, LocalTag> onezero = Point2DBase<float, LocalTag>(1, 0);
0166   const Point2DBase<float, LocalTag> zeroone = Point2DBase<float, LocalTag>(0, 1);
0167   //added by Jason
0168   float localxdotglobalphi_[2];
0169   float localxdotglobalr_[2];
0170   float localxdotglobalz_[2];
0171   float localxdotglobalx_[2];
0172   float localxdotglobaly_[2];
0173   float localydotglobalphi_[2];
0174   float localydotglobalr_[2];
0175   float localydotglobalz_[2];
0176   float localydotglobalx_[2];
0177   float localydotglobaly_[2];
0178 };
0179
0180 //
0181 // constants, enums and typedefs
0182 //
0183
0184 //
0185 // static data member definitions
0186 //
0187
0188 using std::vector;
0189 //
0190 // constructors and destructor
0191 //
0192 OverlapValidation::OverlapValidation(const edm::ParameterSet& iConfig)
0193     : geomToken_(esConsumes()),
0194       magFieldToken_(esConsumes()),
0195       topoToken_(esConsumes()),
0196       config_(iConfig),
0197       rootTree_(nullptr),
0198       FileInPath_(SiStripDetInfoFileReader::kDefaultFile),
0199       compressionSettings_(iConfig.getUntrackedParameter<int>("compressionSettings", -1)),
0200       addExtraBranches_(false),
0201       minHitsCut_(6),
0202       chi2ProbCut_(0.001) {
0203   edm::ConsumesCollector&& iC = consumesCollector();
0204   //now do what ever initialization is needed
0205   trajectoryTag_ = iConfig.getParameter<edm::InputTag>("trajectories");
0206   trajectoryToken_ = iC.consumes<TrajectoryCollection>(trajectoryTag_);
0207   doSimHit_ = iConfig.getParameter<bool>("associateStrip");
0208   detInfo_ = SiStripDetInfoFileReader::read(FileInPath_.fullPath());
0209
0210   overlapCounts_[0] = 0;  // #trajectories
0211   overlapCounts_[1] = 0;  // #hits
0212   overlapCounts_[2] = 0;  // #overlap hits
0213   acceptLayer.resize(7, false);
0214   acceptLayer[PixelSubdetector::PixelBarrel] = iConfig.getParameter<bool>("usePXB");
0215   acceptLayer[PixelSubdetector::PixelEndcap] = iConfig.getParameter<bool>("usePXF");
0216   acceptLayer[StripSubdetector::TIB] = iConfig.getParameter<bool>("useTIB");
0217   acceptLayer[StripSubdetector::TOB] = iConfig.getParameter<bool>("useTOB");
0218   acceptLayer[StripSubdetector::TID] = iConfig.getParameter<bool>("useTID");
0219   acceptLayer[StripSubdetector::TEC] = iConfig.getParameter<bool>("useTEC");
0220   barrelOnly_ = iConfig.getParameter<bool>("barrelOnly");
0221
0222   edm::Service<TFileService> fs;
0223   //
0224   // root output
0225   //
0226   if (compressionSettings_ > 0) {
0227     fs->file().SetCompressionSettings(compressionSettings_);
0228   }
0229
0230   rootTree_ = fs->make<TTree>("Overlaps", "Overlaps");
0231   if (addExtraBranches_) {
0232     rootTree_->Branch("hitCounts", hitCounts_, "found/s:lost/s");
0233     rootTree_->Branch("chi2", chi2_, "chi2/F:ndf/F");
0234     rootTree_->Branch("path", &overlapPath_, "path/F");
0235   }
0236   rootTree_->Branch("layer", &layer_, "layer/i");
0237   rootTree_->Branch("detids", overlapIds_, "id[2]/i");
0238   rootTree_->Branch("predPos", predictedPositions_, "gX[2]/F:gY[2]/F:gZ[2]/F");
0239   rootTree_->Branch("predPar", predictedLocalParameters_, "predQP[2]/F:predDX[2]/F:predDY[2]/F:predX[2]/F:predY[2]/F");
0240   rootTree_->Branch("predErr", predictedLocalErrors_, "predEQP[2]/F:predEDX[2]/F:predEDY[2]/F:predEX[2]/F:predEY[2]/F");
0241   rootTree_->Branch("predEDeltaX", &predictedDeltaXError_, "sigDeltaX/F");
0242   rootTree_->Branch("predEDeltaY", &predictedDeltaYError_, "sigDeltaY/F");
0243   rootTree_->Branch("relSignX", &relativeXSign_, "relSignX/B");
0244   rootTree_->Branch("relSignY", &relativeYSign_, "relSignY/B");
0245   rootTree_->Branch("hitX", hitPositions_, "hitX[2]/F");
0246   rootTree_->Branch("hitEX", hitErrors_, "hitEX[2]/F");
0247   rootTree_->Branch("hitY", hitPositionsY_, "hitY[2]/F");
0248   rootTree_->Branch("hitEY", hitErrorsY_, "hitEY[2]/F");
0249   if (addExtraBranches_) {
0250     rootTree_->Branch("simX", simHitPositions_, "simX[2]/F");
0251     rootTree_->Branch("simY", simHitPositionsY_, "simY[2]/F");
0252     rootTree_->Branch("clusterSize", clusterSize_, "clusterSize[2]/F");
0253     rootTree_->Branch("clusterWidthX", clusterWidthX_, "clusterWidthX[2]/F");
0254     rootTree_->Branch("clusterWidthY", clusterWidthY_, "clusterWidthY[2]/F");
0255     rootTree_->Branch("clusterCharge", clusterCharge_, "clusterCharge[2]/i");
0256     rootTree_->Branch("edge", edge_, "edge[2]/I");
0257   }
0258   rootTree_->Branch("momentum", &momentum_, "momentum/F");
0259   rootTree_->Branch("run", &run_, "run/i");
0260   rootTree_->Branch("event", &event_, "event/i");
0261   rootTree_->Branch("subdetID", &subdetID, "subdetID/I");
0262   rootTree_->Branch("moduleX", moduleX_, "moduleX[2]/F");
0263   rootTree_->Branch("moduleY", moduleY_, "moduleY[2]/F");
0264   rootTree_->Branch("moduleZ", moduleZ_, "moduleZ[2]/F");
0265   rootTree_->Branch("localxdotglobalphi", localxdotglobalphi_, "localxdotglobalphi[2]/F");
0266   rootTree_->Branch("localxdotglobalr", localxdotglobalr_, "localxdotglobalr[2]/F");
0267   rootTree_->Branch("localxdotglobalz", localxdotglobalz_, "localxdotglobalz[2]/F");
0268   rootTree_->Branch("localxdotglobalx", localxdotglobalx_, "localxdotglobalx[2]/F");
0269   rootTree_->Branch("localxdotglobaly", localxdotglobaly_, "localxdotglobaly[2]/F");
0270   rootTree_->Branch("localydotglobalphi", localydotglobalphi_, "localydotglobalphi[2]/F");
0271   rootTree_->Branch("localydotglobalr", localydotglobalr_, "localydotglobalr[2]/F");
0272   rootTree_->Branch("localydotglobalz", localydotglobalz_, "localydotglobalz[2]/F");
0273   rootTree_->Branch("localydotglobalx", localydotglobalx_, "localydotglobalx[2]/F");
0274   rootTree_->Branch("localydotglobaly", localydotglobaly_, "localydotglobaly[2]/F");
0275 }
0276
0277 OverlapValidation::~OverlapValidation() {
0278   edm::LogWarning w("Overlaps");
0279   // do anything here that needs to be done at desctruction time
0280   // (e.g. close files, deallocate resources etc.)
0281
0282   w << "Counters =";
0283   w << " Number of tracks: " << overlapCounts_[0];
0284   w << " Number of valid hits: " << overlapCounts_[1];
0285   w << " Number of overlaps: " << overlapCounts_[2];
0286 }
0287
0288 //
0289 // member functions
0290 //
0291
0292 // ------------ method called to for each event  ------------
0293 void OverlapValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0294   using namespace edm;
0295   //
0296   // mag field & search tracker
0297   //
0298   magField_ = &iSetup.getData(magFieldToken_);
0299   //
0300   // propagator
0301   //
0302   AnalyticalPropagator propagator(magField_, anyDirection);
0303   //
0304   // geometry
0305   //
0306   trackerGeometry_ = &iSetup.getData(geomToken_);
0307   //
0308   // make associator for SimHits
0309   //
0310   TrackerHitAssociator* associator;
0311   if (doSimHit_) {
0312     TrackerHitAssociator::Config hitassociatorconfig(config_, consumesCollector());
0313     associator = new TrackerHitAssociator(iEvent, hitassociatorconfig);
0314   } else {
0315     associator = nullptr;
0316   }
0317
0318   //if(doSimHit_) associator = new TrackerHitAssociator(iEvent, config_); else associator = 0;
0319
0320   //
0321   // trajectories (from refit)
0322   //
0323   //typedef vector<Trajectory> TrajectoryCollection;
0324   edm::Handle<TrajectoryCollection> trajectoryCollectionHandle;
0325   iEvent.getByToken(trajectoryToken_, trajectoryCollectionHandle);
0326   const TrajectoryCollection* const trajectoryCollection = trajectoryCollectionHandle.product();
0327
0328   //
0329   // loop over trajectories from refit
0330   const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0331   for (const auto& trajectory : *trajectoryCollection)
0332     analyzeTrajectory(trajectory, propagator, *associator, tTopo);
0333
0334   run_ = iEvent.id().run();
0335   event_ = iEvent.id().event();
0336 }
0337
0338 void OverlapValidation::analyzeTrajectory(const Trajectory& trajectory,
0339                                           const Propagator& propagator,
0340                                           TrackerHitAssociator& associator,
0341                                           const TrackerTopology* const tTopo) {
0342   typedef std::pair<const TrajectoryMeasurement*, const TrajectoryMeasurement*> Overlap;
0343   typedef vector<Overlap> OverlapContainer;
0344   ++overlapCounts_[0];
0345
0346   OverlapContainer overlapHits;
0347
0348   // quality cuts on trajectory
0349   // min. # hits / matched hits
0350
0351   if (trajectory.foundHits() < minHitsCut_)
0352     return;
0353   if (ChiSquaredProbability((double)(trajectory.chiSquared()), (double)(trajectory.ndof(false))) < chi2ProbCut_)
0354     return;
0355   //
0356   // loop over measurements in the trajectory and calculate residuals
0357   //
0358
0359   vector<TrajectoryMeasurement> measurements(trajectory.measurements());
0360   for (vector<TrajectoryMeasurement>::const_iterator itm = measurements.begin(); itm != measurements.end(); ++itm) {
0361     //
0362     // skip "invalid" (i.e. missing) hits
0363     //
0364     ConstRecHitPointer hit = itm->recHit();
0365     DetId id = hit->geographicalId();
0366     int layer(layerFromId(id, tTopo));
0367     int subDet = id.subdetId();
0368
0369     if (!hit->isValid()) {
0370       edm::LogVerbatim("OverlapValidation") << "Invalid";
0371       continue;
0372     }
0373     if (barrelOnly_ && (subDet == StripSubdetector::TID || subDet == StripSubdetector::TEC))
0374       return;
0375
0376     //edm::LogVerbatim("OverlapValidation")  << "Check " <<subDet << ", layer = " << layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(id).stereo()):2);
0377     //cout << "Check SubID " <<subDet << ", layer = " << layer<<" stereo: "<< ((subDet > 2)?(SiStripDetId(id).stereo()):2) << endl;
0378
0379     //
0380     // check for overlap: same subdet-id && layer number for
0381     // two consecutive hits
0382     //
0383     ++overlapCounts_[1];
0384     if ((layer != -1) && (acceptLayer[subDet])) {
0385       for (vector<TrajectoryMeasurement>::const_iterator itmCompare = itm - 1;
0386            itmCompare >= measurements.begin() && itmCompare > itm - 4;
0387            --itmCompare) {
0388         DetId compareId = itmCompare->recHit()->geographicalId();
0389
0390         if (subDet != compareId.subdetId() || layer != layerFromId(compareId, tTopo))
0391           break;
0392         if (!itmCompare->recHit()->isValid())
0393           continue;
0394         if ((subDet == PixelSubdetector::PixelBarrel || subDet == PixelSubdetector::PixelEndcap) ||
0395             (SiStripDetId(id).stereo() == SiStripDetId(compareId).stereo())) {
0396           overlapHits.push_back(std::make_pair(&(*itmCompare), &(*itm)));
0397           //edm::LogVerbatim("OverlapValidation") << "adding pair "<< ((subDet > 2)?(SiStripDetId(id).stereo()) : 2)
0398           //                 << " from layer = " << layer;
0399           //cout << "adding pair "<< ((subDet > 2)?(SiStripDetId(id).stereo()) : 2) << " from subDet = " << subDet << " and layer = " << layer;
0400           //cout << " \t"<<run_<< "\t"<<event_<<"\t";
0401           //cout << min(id.rawId(),compareId.rawId())<<"\t"<<max(id.rawId(),compareId.rawId())<<endl;
0402           if (SiStripDetId(id).glued() == id.rawId())
0403             edm::LogInfo("Overlaps") << "BAD GLUED: Have glued layer with id = " << id.rawId()
0404                                      << " and glued id = " << SiStripDetId(id).glued()
0405                                      << "  and stereo = " << SiStripDetId(id).stereo() << endl;
0406           if (SiStripDetId(compareId).glued() == compareId.rawId())
0407             edm::LogInfo("Overlaps") << "BAD GLUED: Have glued layer with id = " << compareId.rawId()
0408                                      << " and glued id = " << SiStripDetId(compareId).glued()
0409                                      << "  and stereo = " << SiStripDetId(compareId).stereo() << endl;
0410           break;
0411         }
0412       }
0413     }
0414   }
0415
0416   //
0417   // Loop over all overlap pairs.
0418   //
0419   hitCounts_[0] = trajectory.foundHits();
0420   hitCounts_[1] = trajectory.lostHits();
0421   chi2_[0] = trajectory.chiSquared();
0422   chi2_[1] = trajectory.ndof(false);
0423
0424   for (const auto& overlapHit : overlapHits) {
0425     //
0426     // create reference state @ module 1 (no info from overlap hits)
0427     //
0428     ++overlapCounts_[2];
0429     // backward predicted state at module 1
0430     TrajectoryStateOnSurface bwdPred1 = overlapHit.first->backwardPredictedState();
0431     if (!bwdPred1.isValid())
0432       continue;
0433     //cout << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl;
0434     // forward predicted state at module 2
0435     TrajectoryStateOnSurface fwdPred2 = overlapHit.second->forwardPredictedState();
0436     //cout << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl;
0437     if (!fwdPred2.isValid())
0438       continue;
0439     // extrapolate fwdPred2 to module 1
0440     TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2, bwdPred1.surface());
0441     if (!fwdPred2At1.isValid())
0442       continue;
0443     // combine fwdPred2At1 with bwdPred1 (ref. state, best estimate without hits 1 and 2)
0444     TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1, fwdPred2At1);
0445     if (!comb1.isValid())
0446       continue;
0447     //
0448     // propagation of reference parameters to module 2
0449     //
0450     std::pair<TrajectoryStateOnSurface, double> tsosWithS = propagator.propagateWithPath(comb1, fwdPred2.surface());
0451     TrajectoryStateOnSurface comb1At2 = tsosWithS.first;
0452     if (!comb1At2.isValid())
0453       continue;
0454     //distance of propagation from one surface to the next==could cut here
0455     overlapPath_ = tsosWithS.second;
0456     if (abs(overlapPath_) > 15)
0457       continue;  //cut to remove hit pairs > 15 cm apart
0458     // global position on module 1
0459     GlobalPoint position = comb1.globalPosition();
0460     predictedPositions_[0][0] = position.x();
0461     predictedPositions_[1][0] = position.y();
0462     predictedPositions_[2][0] = position.z();
0463     momentum_ = comb1.globalMomentum().mag();
0464     //cout << "momentum from combination = " << momentum_ << endl;
0465     //cout << "magnetic field from TSOS = " << comb1.magneticField()->inTesla(position).mag() << endl;
0466     // local parameters and errors on module 1
0467     AlgebraicVector5 pars = comb1.localParameters().vector();
0468     AlgebraicSymMatrix55 errs = comb1.localError().matrix();
0469     for (int i = 0; i < 5; ++i) {
0470       predictedLocalParameters_[i][0] = pars[i];
0471       predictedLocalErrors_[i][0] = sqrt(errs(i, i));
0472     }
0473     // global position on module 2
0474     position = comb1At2.globalPosition();
0475     predictedPositions_[0][1] = position.x();
0476     predictedPositions_[1][1] = position.y();
0477     predictedPositions_[2][1] = position.z();
0478     // local parameters and errors on module 2
0479     pars = comb1At2.localParameters().vector();
0480     errs = comb1At2.localError().matrix();
0481     for (int i = 0; i < 5; ++i) {
0482       predictedLocalParameters_[i][1] = pars[i];
0483       predictedLocalErrors_[i][1] = sqrt(errs(i, i));
0484     }
0485
0486     //print out local errors in X to check
0487     //cout << "Predicted local error in X at 1 = " << predictedLocalErrors_[3][0] << "   and predicted local error in X at 2 is = " <<  predictedLocalErrors_[3][1] << endl;
0488     //cout << "Predicted local error in Y at 1 = " << predictedLocalErrors_[4][0] << "   and predicted local error in Y at 2 is = " <<  predictedLocalErrors_[4][1] << endl;
0489
0490     //
0491     // jacobians (local-to-global@1,global 1-2,global-to-local@2)
0492     //
0493     JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(), comb1.localParameters(), *magField_);
0494     AnalyticalCurvilinearJacobian jacCurvToCurv(
0495         comb1.globalParameters(), comb1At2.globalPosition(), comb1At2.globalMomentum(), tsosWithS.second);
0496     JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(), comb1At2.localParameters(), *magField_);
0497     // combined jacobian local-1-to-local-2
0498     AlgebraicMatrix55 jacobian = jacLocToCurv.jacobian() * jacCurvToCurv.jacobian() * jacCurvToLoc.jacobian();
0499     // covariance on module 1
0500     AlgebraicSymMatrix55 covComb1 = comb1.localError().matrix();
0501     // variance and correlations for predicted local_x on modules 1 and 2
0502     double c00 = covComb1(3, 3);
0503     double c10(0.);
0504     double c11(0.);
0505     for (int i = 1; i < 5; ++i) {
0506       c10 += jacobian(3, i) * covComb1(i, 3);
0507       for (int j = 1; j < 5; ++j)
0508         c11 += jacobian(3, i) * covComb1(i, j) * jacobian(3, j);
0509     }
0510     // choose relative sign in order to minimize error on difference
0511     double diff = c00 - 2 * fabs(c10) + c11;
0512     diff = diff > 0 ? sqrt(diff) : -sqrt(-diff);
0513     predictedDeltaXError_ = diff;
0514     relativeXSign_ = c10 > 0 ? -1 : 1;
0515     //
0516     // now find variance and correlations for predicted local_y
0517     double c00Y = covComb1(4, 4);
0518     double c10Y(0.);
0519     double c11Y(0.);
0520     for (int i = 1; i < 5; ++i) {
0521       c10Y += jacobian(4, i) * covComb1(i, 4);
0522       for (int j = 1; j < 5; ++j)
0523         c11Y += jacobian(4, i) * covComb1(i, j) * jacobian(4, j);
0524     }
0525     double diffY = c00Y - 2 * fabs(c10Y) + c11Y;
0526     diffY = diffY > 0 ? sqrt(diffY) : -sqrt(-diffY);
0527     predictedDeltaYError_ = diffY;
0528     relativeYSign_ = c10Y > 0 ? -1 : 1;
0529
0530     // information on modules and hits
0531     overlapIds_[0] = overlapHit.first->recHit()->geographicalId().rawId();
0532     overlapIds_[1] = overlapHit.second->recHit()->geographicalId().rawId();
0533
0534     //added by Heshy and Jared
0535     moduleX_[0] = overlapHit.first->recHit()->det()->surface().position().x();
0536     moduleX_[1] = overlapHit.second->recHit()->det()->surface().position().x();
0537     moduleY_[0] = overlapHit.first->recHit()->det()->surface().position().y();
0538     moduleY_[1] = overlapHit.second->recHit()->det()->surface().position().y();
0539     moduleZ_[0] = overlapHit.first->recHit()->det()->surface().position().z();
0540     moduleZ_[1] = overlapHit.second->recHit()->det()->surface().position().z();
0541     subdetID = overlapHit.first->recHit()->geographicalId().subdetId();
0542     localxdotglobalphi_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(onezero).phi() -
0543                              overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).phi();
0544     localxdotglobalphi_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(onezero).phi() -
0545                              overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).phi();
0546     //added by Jason
0547     localxdotglobalr_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(onezero).perp() -
0548                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).perp();
0549     localxdotglobalr_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(onezero).perp() -
0550                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).perp();
0551     localxdotglobalz_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(onezero).z() -
0552                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).z();
0553     localxdotglobalz_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(onezero).z() -
0554                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).z();
0555     localxdotglobalx_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(onezero).x() -
0556                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).x();
0557     localxdotglobalx_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(onezero).x() -
0558                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).x();
0559     localxdotglobaly_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(onezero).y() -
0560                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).y();
0561     localxdotglobaly_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(onezero).y() -
0562                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).y();
0563     localydotglobalr_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(zeroone).perp() -
0564                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).perp();
0565     localydotglobalr_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(zeroone).perp() -
0566                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).perp();
0567     localydotglobalz_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(zeroone).z() -
0568                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).z();
0569     localydotglobalz_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(zeroone).z() -
0570                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).z();
0571     localydotglobalx_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(zeroone).x() -
0572                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).x();
0573     localydotglobalx_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(zeroone).x() -
0574                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).x();
0575     localydotglobaly_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(zeroone).y() -
0576                            overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).y();
0577     localydotglobaly_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(zeroone).y() -
0578                            overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).y();
0579     localydotglobalphi_[0] = overlapHit.first->recHit()->det()->surface().toGlobal(zeroone).phi() -
0580                              overlapHit.first->recHit()->det()->surface().toGlobal(zerozero).phi();
0581     localydotglobalphi_[1] = overlapHit.second->recHit()->det()->surface().toGlobal(zeroone).phi() -
0582                              overlapHit.second->recHit()->det()->surface().toGlobal(zerozero).phi();
0583
0584     if (overlapHit.first->recHit()->geographicalId().subdetId() == StripSubdetector::TIB)
0585       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo);
0586     else if (overlapHit.first->recHit()->geographicalId().subdetId() == StripSubdetector::TOB)
0587       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 4;
0588     else if (overlapHit.first->recHit()->geographicalId().subdetId() == StripSubdetector::TID)
0589       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 10;
0590     else if (overlapHit.first->recHit()->geographicalId().subdetId() == StripSubdetector::TEC)
0591       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 13;
0592     else if (overlapHit.first->recHit()->geographicalId().subdetId() == PixelSubdetector::PixelBarrel)
0593       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 20;
0594     else if (overlapHit.first->recHit()->geographicalId().subdetId() == PixelSubdetector::PixelEndcap)
0595       layer_ = layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 30;
0596     else
0597       layer_ = 99;
0598
0599     if (overlapIds_[0] == SiStripDetId(overlapHit.first->recHit()->geographicalId()).glued())
0600       edm::LogWarning("Overlaps") << "BAD GLUED: First Id = " << overlapIds_[0] << " has glued = "
0601                                   << SiStripDetId(overlapHit.first->recHit()->geographicalId()).glued()
0602                                   << "  and stereo = "
0603                                   << SiStripDetId(overlapHit.first->recHit()->geographicalId()).stereo() << endl;
0604     if (overlapIds_[1] == SiStripDetId(overlapHit.second->recHit()->geographicalId()).glued())
0605       edm::LogWarning("Overlaps") << "BAD GLUED: Second Id = " << overlapIds_[1] << " has glued = "
0606                                   << SiStripDetId(overlapHit.second->recHit()->geographicalId()).glued()
0607                                   << "  and stereo = "
0608                                   << SiStripDetId(overlapHit.second->recHit()->geographicalId()).stereo() << endl;
0609
0610     const TransientTrackingRecHit::ConstRecHitPointer firstRecHit = overlapHit.first->recHit();
0611     const TransientTrackingRecHit::ConstRecHitPointer secondRecHit = overlapHit.second->recHit();
0612     hitPositions_[0] = firstRecHit->localPosition().x();
0613     hitErrors_[0] = sqrt(firstRecHit->localPositionError().xx());
0614     hitPositions_[1] = secondRecHit->localPosition().x();
0615     hitErrors_[1] = sqrt(secondRecHit->localPositionError().xx());
0616
0617     hitPositionsY_[0] = firstRecHit->localPosition().y();
0618     hitErrorsY_[0] = sqrt(firstRecHit->localPositionError().yy());
0619     hitPositionsY_[1] = secondRecHit->localPosition().y();
0620     hitErrorsY_[1] = sqrt(secondRecHit->localPositionError().yy());
0621
0622     //cout << "printing local X hit position and error for the overlap hits. Hit 1 = " << hitPositions_[0] << "+-" << hitErrors_[0] << "  and hit 2 is "  << hitPositions_[1] << "+-" << hitErrors_[1] << endl;
0623
0624     DetId id1 = overlapHit.first->recHit()->geographicalId();
0625     DetId id2 = overlapHit.second->recHit()->geographicalId();
0626     int layer1 = layerFromId(id1, tTopo);
0627     int subDet1 = id1.subdetId();
0628     int layer2 = layerFromId(id2, tTopo);
0629     int subDet2 = id2.subdetId();
0630     if (abs(hitPositions_[0]) > 5)
0631       edm::LogInfo("Overlaps") << "BAD: Bad hit position: Id = " << id1.rawId()
0632                                << " stereo = " << SiStripDetId(id1).stereo()
0633                                << "  glued = " << SiStripDetId(id1).glued() << " from subdet = " << subDet1
0634                                << " and layer = " << layer1 << endl;
0635     if (abs(hitPositions_[1]) > 5)
0636       edm::LogInfo("Overlaps") << "BAD: Bad hit position: Id = " << id2.rawId()
0637                                << " stereo = " << SiStripDetId(id2).stereo()
0638                                << "  glued = " << SiStripDetId(id2).glued() << " from subdet = " << subDet2
0639                                << " and layer = " << layer2 << endl;
0640
0641     // get track momentum
0642     momentum_ = comb1.globalMomentum().mag();
0643
0644     // get cluster size
0645     if (!(subDet1 == PixelSubdetector::PixelBarrel || subDet1 == PixelSubdetector::PixelEndcap)) {  //strip
0646       const TransientTrackingRecHit::ConstRecHitPointer thit1 = overlapHit.first->recHit();
0647       const SiStripRecHit1D* hit1 = dynamic_cast<const SiStripRecHit1D*>((*thit1).hit());
0648       if (hit1) {
0649         //check cluster width
0650         const SiStripRecHit1D::ClusterRef& cluster1 = hit1->cluster();
0651         clusterSize_[0] = (cluster1->amplitudes()).size();
0652         clusterWidthX_[0] = (cluster1->amplitudes()).size();
0653         clusterWidthY_[0] = -1;
0654
0655         //check if cluster at edge of sensor
0656         uint16_t firstStrip = cluster1->firstStrip();
0657         uint16_t lastStrip = firstStrip + (cluster1->amplitudes()).size() - 1;
0658         unsigned short Nstrips;
0659         Nstrips = detInfo_.getNumberOfApvsAndStripLength(id1).first * 128;
0660         bool atEdge = false;
0661         if (firstStrip == 0 || lastStrip == (Nstrips - 1))
0662           atEdge = true;
0663         if (atEdge)
0664           edge_[0] = 1;
0665         else
0666           edge_[0] = -1;
0667
0668         // get cluster total charge
0669         const auto& stripCharges = cluster1->amplitudes();
0670         uint16_t charge = 0;
0671         for (uint i = 0; i < stripCharges.size(); i++) {
0672           charge += stripCharges[i];
0673         }
0674         clusterCharge_[0] = charge;
0675       }
0676
0677       const TransientTrackingRecHit::ConstRecHitPointer thit2 = overlapHit.second->recHit();
0678       const SiStripRecHit1D* hit2 = dynamic_cast<const SiStripRecHit1D*>((*thit2).hit());
0679       if (hit2) {
0680         const SiStripRecHit1D::ClusterRef& cluster2 = hit2->cluster();
0681         clusterSize_[1] = (cluster2->amplitudes()).size();
0682         clusterWidthX_[1] = (cluster2->amplitudes()).size();
0683         clusterWidthY_[1] = -1;
0684
0685         uint16_t firstStrip = cluster2->firstStrip();
0686         uint16_t lastStrip = firstStrip + (cluster2->amplitudes()).size() - 1;
0687         unsigned short Nstrips;
0688         Nstrips = detInfo_.getNumberOfApvsAndStripLength(id2).first * 128;
0689         bool atEdge = false;
0690         if (firstStrip == 0 || lastStrip == (Nstrips - 1))
0691           atEdge = true;
0692         if (atEdge)
0693           edge_[1] = 1;
0694         else
0695           edge_[1] = -1;
0696
0697         // get cluster total charge
0698         const auto& stripCharges = cluster2->amplitudes();
0699         uint16_t charge = 0;
0700         for (uint i = 0; i < stripCharges.size(); i++) {
0701           charge += stripCharges[i];
0702         }
0703         clusterCharge_[1] = charge;
0704       }
0705       //cout << "strip cluster size2 = " << clusterWidthX_[0] << "  and size 2 = " << clusterWidthX_[1] << endl;
0706     }
0707
0708     if (subDet2 == PixelSubdetector::PixelBarrel || subDet2 == PixelSubdetector::PixelEndcap) {  //pixel
0709
0710       const TransientTrackingRecHit::ConstRecHitPointer thit1 = overlapHit.first->recHit();
0711       const SiPixelRecHit* recHitPix1 = dynamic_cast<const SiPixelRecHit*>((*thit1).hit());
0712       if (recHitPix1) {
0713         // check for cluster size and width
0714         SiPixelRecHit::ClusterRef const& cluster1 = recHitPix1->cluster();
0715
0716         clusterSize_[0] = cluster1->size();
0717         clusterWidthX_[0] = cluster1->sizeX();
0718         clusterWidthY_[0] = cluster1->sizeY();
0719
0720         // check for cluster at edge
0721         const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>((*trackerGeometry_).idToDet(id1));
0722         const PixelTopology* topol = (&(theGeomDet->specificTopology()));
0723
0724         int minPixelRow = cluster1->minPixelRow();  //x
0725         int maxPixelRow = cluster1->maxPixelRow();
0726         int minPixelCol = cluster1->minPixelCol();  //y
0727         int maxPixelCol = cluster1->maxPixelCol();
0728
0729         bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
0730         bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
0731         if (edgeHitX || edgeHitY)
0732           edge_[0] = 1;
0733         else
0734           edge_[0] = -1;
0735
0736         clusterCharge_[0] = (uint)cluster1->charge();
0737
0738       } else {
0739         edm::LogWarning("Overlaps") << "didn't find pixel cluster" << endl;
0740         continue;
0741       }
0742
0743       const TransientTrackingRecHit::ConstRecHitPointer thit2 = overlapHit.second->recHit();
0744       const SiPixelRecHit* recHitPix2 = dynamic_cast<const SiPixelRecHit*>((*thit2).hit());
0745       if (recHitPix2) {
0746         SiPixelRecHit::ClusterRef const& cluster2 = recHitPix2->cluster();
0747
0748         clusterSize_[1] = cluster2->size();
0749         clusterWidthX_[1] = cluster2->sizeX();
0750         clusterWidthY_[1] = cluster2->sizeY();
0751         //cout << "second pixel cluster is " << clusterSize_[1] << " pixels with x width = " << clusterWidthX_[1] << " and y width = " << clusterWidthY_[1] << endl;
0752
0753         const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>((*trackerGeometry_).idToDet(id2));
0754         const PixelTopology* topol = (&(theGeomDet->specificTopology()));
0755
0756         int minPixelRow = cluster2->minPixelRow();  //x
0757         int maxPixelRow = cluster2->maxPixelRow();
0758         int minPixelCol = cluster2->minPixelCol();  //y
0759         int maxPixelCol = cluster2->maxPixelCol();
0760
0761         bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
0762         bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
0763         if (edgeHitX || edgeHitY)
0764           edge_[1] = 1;
0765         else
0766           edge_[1] = -1;
0767
0768         clusterCharge_[1] = (uint)cluster2->charge();
0769
0770       } else {
0771         edm::LogWarning("Overlaps") << "didn't find pixel cluster" << endl;
0772         continue;
0773       }
0774     }
0775
0776     //also check for edge pixels
0777
0778     //try writing out the SimHit info (for MC only)
0779     if (doSimHit_) {
0780       std::vector<PSimHit> psimHits1;
0781       std::vector<PSimHit> psimHits2;
0782       //calculate layer
0783       DetId id = overlapHit.first->recHit()->geographicalId();
0784       int layer(-1);
0785       layer = layerFromId(id, tTopo);
0786       int subDet = id.subdetId();
0787       edm::LogVerbatim("OverlapValidation") << "Subdet = " << subDet << " ; layer = " << layer;
0788
0789       psimHits1 = associator.associateHit(*(firstRecHit->hit()));
0790       edm::LogVerbatim("OverlapValidation") << "single hit ";
0791       edm::LogVerbatim("OverlapValidation") << "length of psimHits1: " << psimHits1.size();
0792       if (!psimHits1.empty()) {
0793         float closest_dist = 99999.9;
0794         std::vector<PSimHit>::const_iterator closest_simhit = psimHits1.begin();
0795         for (std::vector<PSimHit>::const_iterator m = psimHits1.begin(); m < psimHits1.end(); m++) {
0796           //find closest simHit to the recHit
0797           float simX = (*m).localPosition().x();
0798           float dist = fabs(simX - (overlapHit.first->recHit()->localPosition().x()));
0799           edm::LogVerbatim("OverlapValidation")
0800               << "simHit1 simX = " << simX << "   hitX = " << overlapHit.first->recHit()->localPosition().x()
0801               << "   distX = " << dist << "   layer = " << layer;
0802           if (dist < closest_dist) {
0803             //cout << "found newest closest dist for simhit1" << endl;
0804             closest_dist = dist;
0805             closest_simhit = m;
0806           }
0807         }
0808         //if glued layer, convert sim hit position to matchedhit surface
0809         //layer index from 1-4 for TIB, 1-6 for TOB
0810         // Are the sim hits on the glued layers or are they split???
0811         if (subDet > 2 && !SiStripDetId(id).glued()) {
0812           const GluedGeomDet* gluedDet =
0813               (const GluedGeomDet*)(*trackerGeometry_).idToDet((*firstRecHit).hit()->geographicalId());
0814           const StripGeomDetUnit* stripDet = (StripGeomDetUnit*)gluedDet->monoDet();
0815           GlobalPoint gp = stripDet->surface().toGlobal((*closest_simhit).localPosition());
0816           LocalPoint lp = gluedDet->surface().toLocal(gp);
0817           LocalVector localdirection = (*closest_simhit).localDirection();
0818           GlobalVector globaldirection = stripDet->surface().toGlobal(localdirection);
0819           LocalVector direction = gluedDet->surface().toLocal(globaldirection);
0820           float scale = -lp.z() / direction.z();
0821           LocalPoint projectedPos = lp + scale * direction;
0822           simHitPositions_[0] = projectedPos.x();
0823           edm::LogVerbatim("OverlapValidation") << "simhit position from matched layer = " << simHitPositions_[0];
0824           simHitPositionsY_[0] = projectedPos.y();
0825         } else {
0826           simHitPositions_[0] = (*closest_simhit).localPosition().x();
0827           simHitPositionsY_[0] = (*closest_simhit).localPosition().y();
0828           edm::LogVerbatim("OverlapValidation") << "simhit position from non-matched layer = " << simHitPositions_[0];
0829         }
0830         edm::LogVerbatim("OverlapValidation") << "hit position = " << hitPositions_[0];
0831       } else {
0832         simHitPositions_[0] = -99.;
0833         simHitPositionsY_[0] = -99.;
0834         //cout << " filling simHitX: " << -99 << endl;
0835       }
0836
0837       psimHits2 = associator.associateHit(*(secondRecHit->hit()));
0838       if (!psimHits2.empty()) {
0839         float closest_dist = 99999.9;
0840         std::vector<PSimHit>::const_iterator closest_simhit = psimHits2.begin();
0841         for (std::vector<PSimHit>::const_iterator m = psimHits2.begin(); m < psimHits2.end(); m++) {
0842           float simX = (*m).localPosition().x();
0843           float dist = fabs(simX - (overlapHit.second->recHit()->localPosition().x()));
0844           if (dist < closest_dist) {
0845             closest_dist = dist;
0846             closest_simhit = m;
0847           }
0848         }
0849         //if glued layer, convert sim hit position to matchedhit surface
0850         // if no sim hits on matched layers then this section can be removed
0851         if (subDet > 2 && !SiStripDetId(id).glued()) {
0852           const GluedGeomDet* gluedDet =
0853               (const GluedGeomDet*)(*trackerGeometry_).idToDet((*secondRecHit).hit()->geographicalId());
0854           const StripGeomDetUnit* stripDet = (StripGeomDetUnit*)gluedDet->monoDet();
0855           GlobalPoint gp = stripDet->surface().toGlobal((*closest_simhit).localPosition());
0856           LocalPoint lp = gluedDet->surface().toLocal(gp);
0857           LocalVector localdirection = (*closest_simhit).localDirection();
0858           GlobalVector globaldirection = stripDet->surface().toGlobal(localdirection);
0859           LocalVector direction = gluedDet->surface().toLocal(globaldirection);
0860           float scale = -lp.z() / direction.z();
0861           LocalPoint projectedPos = lp + scale * direction;
0862           simHitPositions_[1] = projectedPos.x();
0863           simHitPositionsY_[1] = projectedPos.y();
0864         } else {
0865           simHitPositions_[1] = (*closest_simhit).localPosition().x();
0866           simHitPositionsY_[1] = (*closest_simhit).localPosition().y();
0867         }
0868       } else {
0869         simHitPositions_[1] = -99.;
0870         simHitPositionsY_[1] = -99.;
0871       }
0872     }
0873     rootTree_->Fill();
0874   }
0875 }
0876
0877 int OverlapValidation::layerFromId(const DetId& id, const TrackerTopology* const tTopo) const {
0878   TrackerAlignableId aliid;
0879   std::pair<int, int> subdetandlayer = aliid.typeAndLayerFromDetId(id, tTopo);
0880   int layer = subdetandlayer.second;
0881
0882   return layer;
0883 }
0884
0885 void OverlapValidation::endJob() {
0886   if (rootTree_) {
0887     rootTree_->GetDirectory()->cd();
0888     rootTree_->Write();
0889     delete rootTree_;
0890   }
0891 }
0892
0893 //define this as a plug-in
0894 DEFINE_FWK_MODULE(OverlapValidation);