Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 11:57:05

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