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