File indexing completed on 2021-10-01 22:41:00
0001 #include "HitQuadrupletGeneratorFromLayerPairForPhotonConversion.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0005 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0006 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0007
0008 #include "RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0010 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
0011 #include "RecoTracker/TkHitPairs/interface/OrderedHitPairs.h"
0012
0013 using namespace GeomDetEnumerators;
0014 using namespace std;
0015 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0016
0017
0018
0019 typedef PixelRecoRange<float> Range;
0020 template <class T>
0021 inline T sqr(T t) {
0022 return t * t;
0023 }
0024
0025 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0026 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030 #include "DataFormats/Math/interface/deltaPhi.h"
0031
0032 HitQuadrupletGeneratorFromLayerPairForPhotonConversion::HitQuadrupletGeneratorFromLayerPairForPhotonConversion(
0033 unsigned int inner, unsigned int outer, LayerCacheType *layerCache, unsigned int max)
0034 : theLayerCache(*layerCache), theOuterLayer(outer), theInnerLayer(inner), theMaxElement(max) {
0035 ss = new std::stringstream;
0036 }
0037
0038 void HitQuadrupletGeneratorFromLayerPairForPhotonConversion::hitPairs(const TrackingRegion ®ion,
0039 OrderedHitPairs &result,
0040 const Layers &layers) {
0041 ss->str("");
0042
0043 typedef OrderedHitPair::InnerRecHit InnerHit;
0044 typedef OrderedHitPair::OuterRecHit OuterHit;
0045 typedef RecHitsSortedInPhi::Hit Hit;
0046
0047 Layer innerLayerObj = layers[theInnerLayer];
0048 Layer outerLayerObj = layers[theOuterLayer];
0049
0050 size_t totCountP2 = 0, totCountP1 = 0, totCountM2 = 0, totCountM1 = 0, selCount = 0;
0051 #ifdef mydebug_QSeed
0052 (*ss) << "In " << innerLayerObj.name() << " Out " << theOuterLayer.name() << std::endl;
0053 #endif
0054
0055
0056 const RecHitsSortedInPhi &innerHitsMap = theLayerCache(innerLayerObj, region);
0057 if (innerHitsMap.empty())
0058 return;
0059
0060 const RecHitsSortedInPhi &outerHitsMap = theLayerCache(outerLayerObj, region);
0061 if (outerHitsMap.empty())
0062 return;
0063
0064
0065
0066
0067
0068 vector<RecHitsSortedInPhi::Hit> innerHits;
0069
0070 float innerPhimin, innerPhimax;
0071 float maxDeltaPhi = 1.;
0072
0073
0074 RecHitsSortedInPhi::Range outerHits = outerHitsMap.all();
0075
0076 RecHitsSortedInPhi::HitIter nextoh;
0077 for (RecHitsSortedInPhi::HitIter oh = outerHits.first; oh != outerHits.second; ++oh) {
0078 RecHitsSortedInPhi::Hit ohit = (*oh).hit();
0079 GlobalPoint oPos = ohit->globalPosition();
0080
0081 totCountP2++;
0082 std::unique_ptr<const HitRZCompatibility> checkRZ = region.checkRZ(outerLayerObj.detLayer(), ohit);
0083 for (nextoh = oh + 1; nextoh != outerHits.second; ++nextoh) {
0084 RecHitsSortedInPhi::Hit nohit = (*nextoh).hit();
0085 GlobalPoint noPos = nohit->globalPosition();
0086
0087 #ifdef mydebug_QSeed
0088 (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0089 << oPos.z() / oPos.perp() << std::endl;
0090 (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0091 << noPos.z() / noPos.perp() << std::endl;
0092 (*ss) << "\tDeltaPhi " << reco::deltaPhi(noPos.barePhi(), oPos.barePhi()) << std::endl;
0093 #endif
0094
0095 if (fabs(reco::deltaPhi(noPos.barePhi(), oPos.barePhi())) > maxDeltaPhi)
0096 break;
0097
0098 totCountM2++;
0099
0100
0101 if (failCheckRZCompatibility(nohit, *outerLayerObj.detLayer(), checkRZ.get(), region))
0102 continue;
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113 innerPhimin = ohit->globalPosition().phi();
0114 innerPhimax = nohit->globalPosition().phi();
0115
0116
0117 innerHits.clear();
0118 innerHitsMap.hits(innerPhimin, innerPhimax, innerHits);
0119
0120 #ifdef mydebug_QSeed
0121 (*ss) << "\tiphimin, iphimax " << innerPhimin << " " << innerPhimax << std::endl;
0122 #endif
0123
0124 std::unique_ptr<const HitRZCompatibility> checkRZb = region.checkRZ(innerLayerObj.detLayer(), ohit);
0125 std::unique_ptr<const HitRZCompatibility> checkRZc = region.checkRZ(innerLayerObj.detLayer(), nohit);
0126
0127
0128 vector<RecHitsSortedInPhi::Hit>::const_iterator ieh = innerHits.end();
0129 for (vector<RecHitsSortedInPhi::Hit>::const_iterator ih = innerHits.begin(); ih < ieh; ++ih) {
0130 RecHitsSortedInPhi::Hit ihit = *ih;
0131
0132 #ifdef mydebug_QSeed
0133 GlobalPoint innPos = (*ih)->globalPosition();
0134 (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0135 << oPos.z() / oPos.perp() << std::endl;
0136 (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0137 << noPos.z() / noPos.perp() << std::endl;
0138 (*ss) << "\tinnPos " << innPos << " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta "
0139 << innPos.z() / innPos.perp() << std::endl;
0140 #endif
0141
0142 totCountP1++;
0143
0144
0145 if (failCheckRZCompatibility(ihit, *innerLayerObj.detLayer(), checkRZb.get(), region) ||
0146 failCheckRZCompatibility(ihit, *innerLayerObj.detLayer(), checkRZc.get(), region))
0147 continue;
0148
0149 for (vector<RecHitsSortedInPhi::Hit>::const_iterator nextih = ih + 1; nextih != ieh; ++nextih) {
0150 RecHitsSortedInPhi::Hit nihit = *nextih;
0151
0152 #ifdef mydebug_QSeed
0153 GlobalPoint ninnPos = (*nextih)->globalPosition();
0154 (*ss) << "\toPos " << oPos << " r " << oPos.perp() << " phi " << oPos.phi() << " cotTheta "
0155 << oPos.z() / oPos.perp() << std::endl;
0156 (*ss) << "\tnoPos " << noPos << " r " << noPos.perp() << " phi " << noPos.phi() << " cotTheta "
0157 << noPos.z() / noPos.perp() << std::endl;
0158 (*ss) << "\tinnPos " << innPos << " r " << innPos.perp() << " phi " << innPos.phi() << " cotTheta "
0159 << innPos.z() / innPos.perp() << std::endl;
0160 (*ss) << "\tninnPos " << ninnPos << " r " << ninnPos.perp() << " phi " << ninnPos.phi() << " cotTheta "
0161 << ninnPos.z() / ninnPos.perp() << std::endl;
0162 #endif
0163
0164 totCountM1++;
0165
0166
0167 if (failCheckRZCompatibility(nihit, *innerLayerObj.detLayer(), checkRZb.get(), region) ||
0168 failCheckRZCompatibility(nihit, *innerLayerObj.detLayer(), checkRZc.get(), region))
0169 continue;
0170
0171
0172 if (failCheckSlopeTest(ohit, nohit, ihit, nihit, region))
0173 continue;
0174
0175 if (theMaxElement != 0 && result.size() >= theMaxElement) {
0176 result.clear();
0177 edm::LogError("TooManyQuads") << "number of Quad combinations exceed maximum, no quads produced";
0178 #ifdef mydebug_QSeed
0179 (*ss) << "In " << innerLayerObj.name() << " Out " << outerLayerObj.name() << "\tP2 " << totCountP2
0180 << "\tM2 " << totCountM2 << "\tP1 " << totCountP1 << "\tM1 " << totCountM1 << "\tsel " << selCount
0181 << std::endl;
0182 std::cout << (*ss).str();
0183 #endif
0184 return;
0185 }
0186 selCount++;
0187 result.push_back(OrderedHitPair(ihit, ohit));
0188 result.push_back(OrderedHitPair(nihit, nohit));
0189
0190
0191
0192 }
0193 }
0194 }
0195 }
0196 #ifdef mydebug_QSeed
0197 (*ss) << "In " << innerLayerObj.name() << " Out " << outerLayerObj.name() << "\tP2 " << totCountP2 << "\tM2 "
0198 << totCountM2 << "\tP1 " << totCountP1 << "\tM1 " << totCountM1 << "\tsel " << selCount << std::endl;
0199 std::cout << (*ss).str();
0200 #endif
0201 }
0202
0203 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckRZCompatibility(const RecHitsSortedInPhi::Hit &hit,
0204 const DetLayer &layer,
0205 const HitRZCompatibility *checkRZ,
0206 const TrackingRegion ®ion) {
0207 if (!checkRZ) {
0208 #ifdef mydebug_QSeed
0209 (*ss) << "*******\nNo valid checkRZ\n*******" << std::endl;
0210 #endif
0211 return true;
0212 }
0213
0214 static const float nSigmaRZ = std::sqrt(12.f);
0215 float r_reduced = std::sqrt(sqr(hit->globalPosition().x() - region.origin().x()) +
0216 sqr(hit->globalPosition().y() - region.origin().y()));
0217 Range allowed;
0218 Range hitRZ;
0219 if (layer.location() == barrel) {
0220 allowed = checkRZ->range(r_reduced);
0221 float zErr = nSigmaRZ * hit->errorGlobalZ();
0222 hitRZ = Range(hit->globalPosition().z() - zErr, hit->globalPosition().z() + zErr);
0223 } else {
0224 allowed = checkRZ->range(hit->globalPosition().z());
0225 float rErr = nSigmaRZ * hit->errorGlobalR();
0226 hitRZ = Range(r_reduced - rErr, r_reduced + rErr);
0227 }
0228 Range crossRange = allowed.intersection(hitRZ);
0229
0230 #ifdef mydebug_QSeed
0231 (*ss) << "\n\t\t allowed Range " << allowed.min() << " \t, " << allowed.max() << "\n\t\t hitRz Range "
0232 << hitRZ.min() << " \t, " << hitRZ.max() << "\n\t\t Cross Range " << crossRange.min() << " \t, "
0233 << crossRange.max() << std::endl;
0234 if (!crossRange.empty())
0235 (*ss) << "\n\t\t !!!!ACCEPTED!!! \n\n";
0236 #endif
0237
0238 return crossRange.empty();
0239 }
0240
0241 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckSlopeTest(const RecHitsSortedInPhi::Hit &ohit,
0242 const RecHitsSortedInPhi::Hit &nohit,
0243 const RecHitsSortedInPhi::Hit &ihit,
0244 const RecHitsSortedInPhi::Hit &nihit,
0245 const TrackingRegion ®ion) {
0246 double r[5], z[5], ez[5];
0247
0248
0249
0250
0251 r[0] = ohit->globalPosition().perp();
0252 z[0] = ohit->globalPosition().z();
0253 ez[0] = getEffectiveErrorOnZ(ohit, region);
0254
0255 r[1] = nohit->globalPosition().perp();
0256 z[1] = nohit->globalPosition().z();
0257 ez[1] = getEffectiveErrorOnZ(nohit, region);
0258
0259 r[2] = nihit->globalPosition().perp();
0260 z[2] = nihit->globalPosition().z();
0261 ez[2] = getEffectiveErrorOnZ(nihit, region);
0262
0263 r[3] = ihit->globalPosition().perp();
0264 z[3] = ihit->globalPosition().z();
0265 ez[3] = getEffectiveErrorOnZ(ihit, region);
0266
0267
0268 bubbleSortVsR(4, r, z, ez);
0269
0270
0271 r[4] = region.origin().perp();
0272 z[4] = region.origin().z();
0273 double vError = region.originZBound();
0274 if (vError > 15.)
0275 vError = 1.;
0276 ez[4] = 3. * vError;
0277
0278
0279
0280
0281
0282 double rInn = r[4];
0283 double zInnMin = z[4] - ez[4];
0284 double zInnMax = z[4] + ez[4];
0285
0286
0287 double rOut = r[3];
0288 double zOutMin = z[3] - ez[3];
0289 double zOutMax = z[3] + ez[3];
0290 double rInt = r[2];
0291 double zIntMin = z[2] - ez[2];
0292 double zIntMax = z[2] + ez[2];
0293 if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0294 return true;
0295
0296
0297 rOut = rInt;
0298 zOutMin = zIntMin;
0299 zOutMax = zIntMax;
0300 rInt = r[1];
0301 zIntMin = z[1] - ez[1];
0302 zIntMax = z[1] + ez[1];
0303 if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0304 return true;
0305
0306
0307 rOut = rInt;
0308 zOutMin = zIntMin;
0309 zOutMax = zIntMax;
0310 rInt = r[0];
0311 zIntMin = z[0] - ez[0];
0312 zIntMax = z[0] + ez[0];
0313 if (failCheckSegmentZCompatibility(rInn, zInnMin, zInnMax, rInt, zIntMin, zIntMax, rOut, zOutMin, zOutMax))
0314 return true;
0315
0316
0317
0318 return false;
0319 }
0320
0321 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::verySimpleFit(
0322 int size, double *ax, double *ay, double *e2y, double &p0, double &e2p0, double &p1) {
0323
0324 return 0;
0325 }
0326
0327 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getSqrEffectiveErrorOnZ(
0328 const RecHitsSortedInPhi::Hit &hit, const TrackingRegion ®ion) {
0329
0330
0331
0332
0333 double sqrProjFactor =
0334 sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
0335 return (hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
0336 }
0337
0338 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getEffectiveErrorOnZ(const RecHitsSortedInPhi::Hit &hit,
0339 const TrackingRegion ®ion) {
0340
0341
0342
0343 double sqrProjFactor =
0344 sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
0345 double effErr =
0346 sqrt(hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
0347 if (effErr > 2.) {
0348 effErr *= 1.8;
0349
0350
0351
0352 } else {
0353 effErr *= 2.;
0354
0355 }
0356 return effErr;
0357 }
0358
0359 void HitQuadrupletGeneratorFromLayerPairForPhotonConversion::bubbleSortVsR(int n, double *ar, double *az, double *aez) {
0360 bool swapped = true;
0361 int j = 0;
0362 double tmpr, tmpz, tmpez;
0363 while (swapped) {
0364 swapped = false;
0365 j++;
0366 for (int i = 0; i < n - j; i++) {
0367 if (ar[i] > ar[i + 1]) {
0368 tmpr = ar[i];
0369 ar[i] = ar[i + 1];
0370 ar[i + 1] = tmpr;
0371 tmpz = az[i];
0372 az[i] = az[i + 1];
0373 az[i + 1] = tmpz;
0374 tmpez = aez[i];
0375 aez[i] = aez[i + 1];
0376 aez[i + 1] = tmpez;
0377 swapped = true;
0378 }
0379 }
0380 }
0381 }
0382
0383 bool HitQuadrupletGeneratorFromLayerPairForPhotonConversion::failCheckSegmentZCompatibility(double &rInn,
0384 double &zInnMin,
0385 double &zInnMax,
0386 double &rInt,
0387 double &zIntMin,
0388 double &zIntMax,
0389 double &rOut,
0390 double &zOutMin,
0391 double &zOutMax) {
0392
0393
0394
0395
0396
0397 double zLeft = getZAtR(rInn, zInnMin, rInt, rOut, zOutMin);
0398 if (zIntMax < zLeft)
0399 return true;
0400
0401 double zRight = getZAtR(rInn, zInnMax, rInt, rOut, zOutMax);
0402 if (zIntMin > zRight)
0403 return true;
0404 if (zIntMin < zLeft && zIntMax < zRight) {
0405 zIntMax = zLeft;
0406 return false;
0407 }
0408 if (zIntMin > zLeft && zIntMax > zRight) {
0409 zIntMax = zRight;
0410 return false;
0411 }
0412
0413
0414 return false;
0415 }
0416
0417 double HitQuadrupletGeneratorFromLayerPairForPhotonConversion::getZAtR(
0418 double &rInn, double &zInn, double &r, double &rOut, double &zOut) {
0419
0420
0421
0422
0423 return zInn + (zOut - zInn) * (r - rInn) / (rOut - rInn);
0424 }