Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-22 04:02:33 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!