99 void endJob()
override;
118 int overlapCounts_[3];
129 unsigned short hitCounts_[2];
131 unsigned int overlapIds_[2];
132 float predictedPositions_[3][2];
133 float predictedLocalParameters_[5][2];
134 float predictedLocalErrors_[5][2];
139 float hitPositions_[2];
141 float hitPositionsY_[2];
142 float hitErrorsY_[2];
143 float simHitPositions_[2];
144 float simHitPositionsY_[2];
145 float clusterWidthX_[2];
146 float clusterWidthY_[2];
147 float clusterSize_[2];
165 float localxdotglobalphi_[2];
166 float localxdotglobalr_[2];
167 float localxdotglobalz_[2];
168 float localxdotglobalx_[2];
169 float localxdotglobaly_[2];
170 float localydotglobalphi_[2];
171 float localydotglobalr_[2];
172 float localydotglobalz_[2];
173 float localydotglobalx_[2];
174 float localydotglobaly_[2];
195 FileInPath_(
"CalibTracker/SiStripCommon/data/SiStripDetInfo.dat"),
196 compressionSettings_(iConfig.getUntrackedParameter<
int>(
"compressionSettings", -1)),
197 addExtraBranches_(
false),
199 chi2ProbCut_(0.001) {
227 rootTree_ = fs->make<TTree>(
"Overlaps",
"Overlaps");
328 for (
const auto& trajectory : *trajectoryCollection)
339 typedef std::pair<const TrajectoryMeasurement*, const TrajectoryMeasurement*>
Overlap;
340 typedef vector<Overlap> OverlapContainer;
343 OverlapContainer overlapHits;
356 vector<TrajectoryMeasurement> measurements(trajectory.
measurements());
357 for (vector<TrajectoryMeasurement>::const_iterator itm = measurements.begin(); itm != measurements.end(); ++itm) {
364 int subDet =
id.subdetId();
366 if (!
hit->isValid()) {
382 for (vector<TrajectoryMeasurement>::const_iterator itmCompare = itm - 1;
383 itmCompare >= measurements.begin() && itmCompare > itm - 4;
385 DetId compareId = itmCompare->recHit()->geographicalId();
389 if (!itmCompare->recHit()->isValid())
393 overlapHits.push_back(std::make_pair(&(*itmCompare), &(*itm)));
400 edm::LogInfo(
"Overlaps") <<
"BAD GLUED: Have glued layer with id = " <<
id.rawId()
404 edm::LogInfo(
"Overlaps") <<
"BAD GLUED: Have glued layer with id = " << compareId.
rawId()
421 for (
const auto& overlapHit : overlapHits) {
447 std::pair<TrajectoryStateOnSurface, double> tsosWithS =
propagator.propagateWithPath(comb1, fwdPred2.
surface());
466 for (
int i = 0;
i < 5; ++
i) {
478 for (
int i = 0;
i < 5; ++
i) {
499 double c00 = covComb1(3, 3);
502 for (
int i = 1;
i < 5; ++
i) {
503 c10 += jacobian(3,
i) * covComb1(
i, 3);
504 for (
int j = 1;
j < 5; ++
j)
505 c11 += jacobian(3,
i) * covComb1(
i,
j) * jacobian(3,
j);
508 double diff = c00 - 2 * fabs(c10) + c11;
514 double c00Y = covComb1(4, 4);
517 for (
int i = 1;
i < 5; ++
i) {
518 c10Y += jacobian(4,
i) * covComb1(
i, 4);
519 for (
int j = 1;
j < 5; ++
j)
520 c11Y += jacobian(4,
i) * covComb1(
i,
j) * jacobian(4,
j);
522 double diffY = c00Y - 2 * fabs(c10Y) + c11Y;
523 diffY = diffY > 0 ?
sqrt(diffY) : -
sqrt(-diffY);
528 overlapIds_[0] = overlapHit.first->recHit()->geographicalId().rawId();
529 overlapIds_[1] = overlapHit.second->recHit()->geographicalId().rawId();
532 moduleX_[0] = overlapHit.first->recHit()->det()->surface().position().x();
533 moduleX_[1] = overlapHit.second->recHit()->det()->surface().position().x();
534 moduleY_[0] = overlapHit.first->recHit()->det()->surface().position().y();
535 moduleY_[1] = overlapHit.second->recHit()->det()->surface().position().y();
536 moduleZ_[0] = overlapHit.first->recHit()->det()->surface().position().z();
537 moduleZ_[1] = overlapHit.second->recHit()->det()->surface().position().z();
538 subdetID = overlapHit.first->recHit()->geographicalId().subdetId();
540 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).phi();
542 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).phi();
545 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).perp();
547 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).perp();
549 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).z();
551 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).z();
553 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).x();
555 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).x();
557 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).y();
559 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).y();
561 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).perp();
563 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).perp();
565 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).z();
567 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).z();
569 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).x();
571 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).x();
573 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).y();
575 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).y();
577 overlapHit.first->recHit()->det()->surface().toGlobal(
zerozero).phi();
579 overlapHit.second->recHit()->det()->surface().toGlobal(
zerozero).phi();
582 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo);
584 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 4;
586 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 10;
588 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 13;
590 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 20;
592 layer_ =
layerFromId(overlapHit.first->recHit()->geographicalId().rawId(), tTopo) + 30;
598 <<
SiStripDetId(overlapHit.first->recHit()->geographicalId()).glued()
600 <<
SiStripDetId(overlapHit.first->recHit()->geographicalId()).stereo() << endl;
603 <<
SiStripDetId(overlapHit.second->recHit()->geographicalId()).glued()
605 <<
SiStripDetId(overlapHit.second->recHit()->geographicalId()).stereo() << endl;
621 DetId id1 = overlapHit.first->recHit()->geographicalId();
622 DetId id2 = overlapHit.second->recHit()->geographicalId();
624 int subDet1 =
id1.subdetId();
626 int subDet2 =
id2.subdetId();
628 edm::LogInfo(
"Overlaps") <<
"BAD: Bad hit position: Id = " <<
id1.rawId()
631 <<
" and layer = " << layer1 << endl;
633 edm::LogInfo(
"Overlaps") <<
"BAD: Bad hit position: Id = " <<
id2.rawId()
636 <<
" and layer = " << layer2 << endl;
644 const SiStripRecHit1D* hit1 = dynamic_cast<const SiStripRecHit1D*>((*thit1).hit());
653 uint16_t firstStrip = cluster1->firstStrip();
654 uint16_t lastStrip = firstStrip + (cluster1->amplitudes()).
size() - 1;
655 unsigned short Nstrips;
658 if (firstStrip == 0 || lastStrip == (Nstrips - 1))
666 const auto& stripCharges = cluster1->amplitudes();
668 for (
uint i = 0;
i < stripCharges.size();
i++) {
675 const SiStripRecHit1D* hit2 = dynamic_cast<const SiStripRecHit1D*>((*thit2).hit());
682 uint16_t firstStrip = cluster2->firstStrip();
683 uint16_t lastStrip = firstStrip + (cluster2->amplitudes()).
size() - 1;
684 unsigned short Nstrips;
687 if (firstStrip == 0 || lastStrip == (Nstrips - 1))
695 const auto& stripCharges = cluster2->amplitudes();
697 for (
uint i = 0;
i < stripCharges.size();
i++) {
708 const SiPixelRecHit* recHitPix1 = dynamic_cast<const SiPixelRecHit*>((*thit1).hit());
718 const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>((*trackerGeometry_).idToDet(
id1));
721 int minPixelRow = cluster1->minPixelRow();
722 int maxPixelRow = cluster1->maxPixelRow();
723 int minPixelCol = cluster1->minPixelCol();
724 int maxPixelCol = cluster1->maxPixelCol();
726 bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
727 bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
728 if (edgeHitX || edgeHitY)
741 const SiPixelRecHit* recHitPix2 = dynamic_cast<const SiPixelRecHit*>((*thit2).hit());
750 const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>((*trackerGeometry_).idToDet(
id2));
753 int minPixelRow = cluster2->minPixelRow();
754 int maxPixelRow = cluster2->maxPixelRow();
755 int minPixelCol = cluster2->minPixelCol();
756 int maxPixelCol = cluster2->maxPixelCol();
758 bool edgeHitX = (topol->isItEdgePixelInX(minPixelRow)) || (topol->isItEdgePixelInX(maxPixelRow));
759 bool edgeHitY = (topol->isItEdgePixelInY(minPixelCol)) || (topol->isItEdgePixelInY(maxPixelCol));
760 if (edgeHitX || edgeHitY)
777 std::vector<PSimHit> psimHits1;
778 std::vector<PSimHit> psimHits2;
780 DetId id = overlapHit.first->recHit()->geographicalId();
783 int subDet =
id.subdetId();
786 psimHits1 =
associator.associateHit(*(firstRecHit->hit()));
788 edm::LogVerbatim(
"OverlapValidation") <<
"length of psimHits1: " << psimHits1.size();
789 if (!psimHits1.empty()) {
790 float closest_dist = 99999.9;
791 std::vector<PSimHit>::const_iterator closest_simhit = psimHits1.begin();
792 for (std::vector<PSimHit>::const_iterator
m = psimHits1.begin();
m < psimHits1.end();
m++) {
794 float simX = (*m).localPosition().x();
795 float dist = fabs(simX - (overlapHit.first->recHit()->localPosition().x()));
797 <<
"simHit1 simX = " << simX <<
" hitX = " << overlapHit.first->recHit()->localPosition().x()
798 <<
" distX = " << dist <<
" layer = " <<
layer;
799 if (dist < closest_dist) {
810 (
const GluedGeomDet*)(*trackerGeometry_).idToDet((*firstRecHit).hit()->geographicalId());
814 LocalVector localdirection = (*closest_simhit).localDirection();
817 float scale = -lp.
z() / direction.
z();
834 psimHits2 =
associator.associateHit(*(secondRecHit->hit()));
835 if (!psimHits2.empty()) {
836 float closest_dist = 99999.9;
837 std::vector<PSimHit>::const_iterator closest_simhit = psimHits2.begin();
838 for (std::vector<PSimHit>::const_iterator
m = psimHits2.begin();
m < psimHits2.end();
m++) {
839 float simX = (*m).localPosition().x();
840 float dist = fabs(simX - (overlapHit.second->recHit()->localPosition().x()));
841 if (dist < closest_dist) {
850 (
const GluedGeomDet*)(*trackerGeometry_).idToDet((*secondRecHit).hit()->geographicalId());
854 LocalVector localdirection = (*closest_simhit).localDirection();
857 float scale = -lp.
z() / direction.
z();
877 int layer = subdetandlayer.second;