File indexing completed on 2024-04-06 12:28:02
0001 #include "SeedForPhotonConversionFromQuadruplets.h"
0002
0003 #include "Conv4HitsReco2.h"
0004
0005 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0006 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0012 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0013 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0014 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0015 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0016 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0017 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0018 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0019
0020 #include "RecoTracker/TkTrackingRegions/interface/OrderedHitsGenerator.h"
0021 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
0022 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0023 #include "RecoTracker/TkSeedGenerator/interface/SeedCreator.h"
0024 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
0025 #include "RecoTracker/TkSeedGenerator/interface/SeedCreatorFactory.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToBeamLine.h"
0027 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0028 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0029
0030
0031
0032
0033 template <class T>
0034 T sqr(T t) {
0035 return t * t;
0036 }
0037
0038 SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets(
0039 const edm::ParameterSet& cfg, edm::ConsumesCollector& iC, const edm::ParameterSet& SeedComparitorPSet)
0040 : SeedForPhotonConversionFromQuadruplets(
0041 iC,
0042 SeedComparitorPSet,
0043 cfg.getParameter<std::string>("propagator"),
0044 cfg.existsAs<double>("SeedMomentumForBOFF") ? cfg.getParameter<double>("SeedMomentumForBOFF") : 5.0) {}
0045
0046 SeedForPhotonConversionFromQuadruplets::SeedForPhotonConversionFromQuadruplets(
0047 edm::ConsumesCollector& iC,
0048 const edm::ParameterSet& SeedComparitorPSet,
0049 const std::string& propagator,
0050 double seedMomentumForBOFF)
0051 : theBfieldToken(iC.esConsumes()),
0052 theTTRHBuilderToken(iC.esConsumes(edm::ESInputTag("", "WithTrackAngle"))),
0053 theTrackerToken(iC.esConsumes()),
0054 thePropagatorToken(iC.esConsumes(edm::ESInputTag("", propagator))),
0055 theBOFFMomentum(seedMomentumForBOFF) {
0056 std::string comparitorName = SeedComparitorPSet.getParameter<std::string>("ComponentName");
0057 if (comparitorName != "none") {
0058 theComparitor = SeedComparitorFactory::get()->create(comparitorName, SeedComparitorPSet, iC);
0059 }
0060 }
0061
0062 SeedForPhotonConversionFromQuadruplets::~SeedForPhotonConversionFromQuadruplets() {}
0063
0064 const TrajectorySeed* SeedForPhotonConversionFromQuadruplets::trajectorySeed(TrajectorySeedCollection& seedCollection,
0065 const SeedingHitSet& phits,
0066 const SeedingHitSet& mhits,
0067 const TrackingRegion& region,
0068 const edm::Event& ev,
0069 const edm::EventSetup& es,
0070 std::stringstream& ss,
0071 std::vector<Quad>& quadV,
0072 edm::ParameterSet& QuadCutPSet) {
0073
0074
0075 bool rejectAllQuads = QuadCutPSet.getParameter<bool>("rejectAllQuads");
0076 if (rejectAllQuads)
0077 return nullptr;
0078
0079 bool applyDeltaPhiCuts = QuadCutPSet.getParameter<bool>("apply_DeltaPhiCuts");
0080 bool ClusterShapeFiltering = QuadCutPSet.getParameter<bool>("apply_ClusterShapeFilter");
0081 bool applyArbitration = QuadCutPSet.getParameter<bool>("apply_Arbitration");
0082 bool applydzCAcut = QuadCutPSet.getParameter<bool>("apply_zCACut");
0083 double CleaningmaxRadialDistance = QuadCutPSet.getParameter<double>("Cut_DeltaRho");
0084 double BeamPipeRadiusCut = QuadCutPSet.getParameter<double>("Cut_BeamPipeRadius");
0085 double CleaningMinLegPt = QuadCutPSet.getParameter<double>("Cut_minLegPt");
0086 double maxLegPt = QuadCutPSet.getParameter<double>("Cut_maxLegPt");
0087 double dzcut = QuadCutPSet.getParameter<double>("Cut_zCA");
0088
0089 double toleranceFactorOnDeltaPhiCuts = 0.1;
0090
0091
0092
0093 pss = &ss;
0094
0095 if (phits.size() < 2)
0096 return nullptr;
0097 if (mhits.size() < 2)
0098 return nullptr;
0099
0100
0101
0102
0103
0104
0105
0106 #ifdef mydebug_sguazz
0107 std::cout << " --------------------------------------------------------------------------"
0108 << "\n";
0109 std::cout << " Starting a hit quad fast reco "
0110 << "\n";
0111 std::cout << " --------------------------------------------------------------------------"
0112 << "\n";
0113 #endif
0114
0115
0116
0117 SeedingHitSet::ConstRecHitPointer ptth1 = phits[0];
0118 SeedingHitSet::ConstRecHitPointer ptth2 = phits[1];
0119 SeedingHitSet::ConstRecHitPointer mtth1 = mhits[0];
0120 SeedingHitSet::ConstRecHitPointer mtth2 = mhits[1];
0121
0122 GlobalPoint vHit[4];
0123 vHit[0] = ptth2->globalPosition();
0124 vHit[1] = ptth1->globalPosition();
0125 vHit[2] = mtth1->globalPosition();
0126 vHit[3] = mtth2->globalPosition();
0127
0128
0129 const GlobalPoint& vgPhotVertex = region.origin();
0130 math::XYZVector vPhotVertex(vgPhotVertex.x(), vgPhotVertex.y(), vgPhotVertex.z());
0131
0132 math::XYZVector h1(vHit[0].x(), vHit[0].y(), vHit[0].z());
0133 math::XYZVector h2(vHit[1].x(), vHit[1].y(), vHit[1].z());
0134 math::XYZVector h3(vHit[2].x(), vHit[2].y(), vHit[2].z());
0135 math::XYZVector h4(vHit[3].x(), vHit[3].y(), vHit[3].z());
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155
0156 math::XYZVector P1;
0157 math::XYZVector P2;
0158 math::XYZVector M1;
0159 math::XYZVector M2;
0160
0161 if (h1.x() * h1.x() + h1.y() * h1.y() < h2.x() * h2.x() + h2.y() * h2.y()) {
0162 P1 = h1;
0163 P2 = h2;
0164 } else {
0165 P1 = h2;
0166 P2 = h1;
0167 }
0168
0169 if (h3.x() * h3.x() + h3.y() * h3.y() < h4.x() * h4.x() + h4.y() * h4.y()) {
0170 M1 = h3;
0171 M2 = h4;
0172 } else {
0173 M1 = h4;
0174 M2 = h3;
0175 }
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187 math::XYZVector IP(0, 0, 0);
0188
0189
0190 double kP = (P1.y() - P2.y()) / (P1.x() - P2.x());
0191 double dP = P1.y() - kP * P1.x();
0192
0193 double kM = (M1.y() - M2.y()) / (M1.x() - M2.x());
0194 double dM = M1.y() - kM * M1.x();
0195
0196 double IPx = (dM - dP) / (kP - kM);
0197 double IPy = kP * IPx + dP;
0198
0199 IP.SetXYZ(IPx, IPy, 0);
0200
0201 double IPrho = std::sqrt(IP.x() * IP.x() + IP.y() * IP.y());
0202 double P1rho2 = P1.x() * P1.x() + P1.y() * P1.y();
0203 double M1rho2 = M1.x() * M1.x() + M1.y() * M1.y();
0204 double maxIPrho2 = IPrho + CleaningmaxRadialDistance;
0205 maxIPrho2 *= maxIPrho2;
0206
0207 if (IPrho < BeamPipeRadiusCut || P1rho2 > maxIPrho2 || M1rho2 > maxIPrho2) {
0208 return nullptr;
0209 }
0210
0211 if (applyDeltaPhiCuts) {
0212 kPI_ = std::atan(1.0) * 4;
0213
0214 const auto& bfield = es.getData(theBfieldToken);
0215 math::XYZVector QuadMean(0, 0, 0);
0216 QuadMean.SetXYZ((M1.x() + M2.x() + P1.x() + P2.x()) / 4.,
0217 (M1.y() + M2.y() + P1.y() + P2.y()) / 4.,
0218 (M1.z() + M2.z() + P1.z() + P2.z()) / 4.);
0219
0220 double fBField = bfield.inTesla(GlobalPoint(QuadMean.x(), QuadMean.y(), QuadMean.z())).z();
0221
0222 double rMax = CleaningMinLegPt / (0.01 * 0.3 * fBField);
0223 double rMax_squared = rMax * rMax;
0224 double Mx = M1.x();
0225 double My = M1.y();
0226
0227 math::XYZVector B(0, 0, 0);
0228 math::XYZVector C(0, 0, 0);
0229
0230 if (rMax_squared * 4. > Mx * Mx + My * My) {
0231
0232
0233
0234
0235
0236
0237
0238
0239 double k = -Mx / My;
0240 double d = My / 2. - k * Mx / 2.;
0241
0242 #ifdef mydebug_knuenz
0243 std::cout << "k" << k << std::endl;
0244 std::cout << "d" << d << std::endl;
0245 #endif
0246
0247
0248 double CsolutionPart1 = -2 * k * d;
0249 double CsolutionPart2 = std::sqrt(4 * k * k * d * d - 4 * (1 + k * k) * (d * d - rMax_squared));
0250 double CsolutionPart3 = 2 * (1 + k * k);
0251 double Cx1 = (CsolutionPart1 + CsolutionPart2) / CsolutionPart3;
0252 double Cx2 = (CsolutionPart1 - CsolutionPart2) / CsolutionPart3;
0253 double Cy1 = k * Cx1 + d;
0254 double Cy2 = k * Cx2 + d;
0255
0256
0257 double Cx, Cy;
0258 math::XYZVector C1(Cx1, Cy1, 0);
0259 if (C1.x() * M1.y() - C1.y() * M1.x() < 0) {
0260 Cx = Cx1;
0261 Cy = Cy1;
0262 } else {
0263 Cx = Cx2;
0264 Cy = Cy2;
0265 }
0266 C.SetXYZ(Cx, Cy, 0);
0267
0268 #ifdef mydebug_knuenz
0269 std::cout << "Cx1" << Cx1 << std::endl;
0270 std::cout << "Cx2" << Cx2 << std::endl;
0271 std::cout << "Cy1" << Cy1 << std::endl;
0272 std::cout << "Cy2" << Cy2 << std::endl;
0273 std::cout << "Cx" << Cx << std::endl;
0274 std::cout << "Cy" << Cy << std::endl;
0275 #endif
0276
0277
0278 k = -Cx / Cy;
0279 d = 0;
0280 double Bx1 = std::sqrt(Mx * Mx + My * My / (1 + k * k));
0281 double Bx2 = -Bx1;
0282 double By1 = k * Bx1 + d;
0283 double By2 = k * Bx2 + d;
0284
0285 #ifdef mydebug_knuenz
0286 std::cout << "k" << k << std::endl;
0287 std::cout << "d" << d << std::endl;
0288 #endif
0289
0290
0291 double Bx, By;
0292 math::XYZVector B1(Bx1, By1, 0);
0293 if (M1.x() * B1.y() - M1.y() * B1.x() < 0) {
0294 Bx = Bx1;
0295 By = By1;
0296 } else {
0297 Bx = Bx2;
0298 By = By2;
0299 }
0300 B.SetXYZ(Bx, By, 0);
0301
0302 #ifdef mydebug_knuenz
0303 std::cout << "Bx1" << Bx1 << std::endl;
0304 std::cout << "Bx2" << Bx2 << std::endl;
0305 std::cout << "By1" << By1 << std::endl;
0306 std::cout << "By2" << By2 << std::endl;
0307 std::cout << "Bx" << Bx << std::endl;
0308 std::cout << "By" << By << std::endl;
0309 #endif
0310
0311 double DeltaPhiMaxM1P1 = DeltaPhiManual(M1, B) * 2;
0312
0313 #ifdef mydebug_knuenz
0314 std::cout << "DeltaPhiMaxM1P1 " << DeltaPhiMaxM1P1 << std::endl;
0315 std::cout << "M1.DeltaPhi(P1) " << DeltaPhiManual(M1, P1) << std::endl;
0316 std::cout << "rho P1: " << std::sqrt(P1.x() * P1.x() + P1.y() * P1.y()) << "phi P1: " << P1.Phi() << std::endl;
0317 std::cout << "rho M1: " << std::sqrt(M1.x() * M1.x() + M1.y() * M1.y()) << "phi M1: " << M1.Phi() << std::endl;
0318 #endif
0319
0320
0321
0322 double tol_DeltaPhiMaxM1P1 = DeltaPhiMaxM1P1 * toleranceFactorOnDeltaPhiCuts;
0323 double DeltaPhiManualM1P1 = DeltaPhiManual(M1, P1);
0324
0325 if (DeltaPhiManualM1P1 > DeltaPhiMaxM1P1 + tol_DeltaPhiMaxM1P1 || DeltaPhiManualM1P1 < 0 - tol_DeltaPhiMaxM1P1) {
0326 return nullptr;
0327 }
0328
0329 }
0330
0331
0332
0333
0334
0335
0336
0337 double rM2_squared = M2.x() * M2.x() + M2.y() * M2.y();
0338 if (rMax_squared * 4. > rM2_squared) {
0339
0340
0341 double k = -C.x() / C.y();
0342 double d = (rM2_squared - rMax_squared + C.x() * C.x() + C.y() * C.y()) / (2 * C.y());
0343
0344 double M2solutionPart1 = -2 * k * d;
0345 double M2solutionPart2 = std::sqrt(4 * k * k * d * d - 4 * (1 + k * k) * (d * d - rM2_squared));
0346 double M2solutionPart3 = 2 + 2 * k * k;
0347 double M2xMax1 = (M2solutionPart1 + M2solutionPart2) / M2solutionPart3;
0348 double M2xMax2 = (M2solutionPart1 - M2solutionPart2) / M2solutionPart3;
0349 double M2yMax1 = k * M2xMax1 + d;
0350 double M2yMax2 = k * M2xMax2 + d;
0351
0352
0353 math::XYZVector M2MaxVec1(M2xMax1, M2yMax1, 0);
0354 math::XYZVector M2MaxVec2(M2xMax2, M2yMax2, 0);
0355 math::XYZVector M2MaxVec(0, 0, 0);
0356 if (M2MaxVec1.x() * M2MaxVec2.y() - M2MaxVec1.y() * M2MaxVec2.x() < 0) {
0357 M2MaxVec.SetXYZ(M2xMax2, M2yMax2, 0);
0358 } else {
0359 M2MaxVec.SetXYZ(M2xMax1, M2yMax1, 0);
0360 }
0361
0362 double DeltaPhiMaxM2 = DeltaPhiManual(M2MaxVec, M1);
0363
0364 #ifdef mydebug_knuenz
0365 std::cout << "C.x() " << C.x() << std::endl;
0366 std::cout << "C.y() " << C.y() << std::endl;
0367 std::cout << "M1.x() " << M1.x() << std::endl;
0368 std::cout << "M1.y() " << M1.y() << std::endl;
0369 std::cout << "M2.x() " << M2.x() << std::endl;
0370 std::cout << "M2.y() " << M2.y() << std::endl;
0371 std::cout << "k " << k << std::endl;
0372 std::cout << "d " << d << std::endl;
0373 std::cout << "M2xMax1 " << M2xMax1 << std::endl;
0374 std::cout << "M2xMax2 " << M2xMax2 << std::endl;
0375 std::cout << "M2yMax1 " << M2yMax1 << std::endl;
0376 std::cout << "M2yMax2 " << M2yMax2 << std::endl;
0377 std::cout << "M2xMax " << M2MaxVec.x() << std::endl;
0378 std::cout << "M2yMax " << M2MaxVec.y() << std::endl;
0379 std::cout << "rM2_squared " << rM2_squared << std::endl;
0380 std::cout << "rMax " << rMax << std::endl;
0381 std::cout << "DeltaPhiMaxM2 " << DeltaPhiMaxM2 << std::endl;
0382 #endif
0383
0384 double tol_DeltaPhiMaxM2 = DeltaPhiMaxM2 * toleranceFactorOnDeltaPhiCuts;
0385 double DeltaPhiManualM2M1 = DeltaPhiManual(M2, M1);
0386
0387 if (DeltaPhiManualM2M1 > DeltaPhiMaxM2 + tol_DeltaPhiMaxM2 || DeltaPhiManualM2M1 < 0 - tol_DeltaPhiMaxM2) {
0388 return nullptr;
0389 }
0390
0391
0392 double DeltaPhiManualP1P2 = DeltaPhiManual(P1, P2);
0393 if (DeltaPhiManualP1P2 > DeltaPhiMaxM2 + tol_DeltaPhiMaxM2 || DeltaPhiManualP1P2 < 0 - tol_DeltaPhiMaxM2) {
0394 return nullptr;
0395 }
0396 }
0397 }
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407 Conv4HitsReco2 quad(vPhotVertex, h1, h2, h3, h4);
0408 quad.SetMaxNumberOfIterations(100);
0409 #ifdef mydebug_sguazz
0410 quad.Dump();
0411 #endif
0412 math::XYZVector candVtx;
0413 double candPtPlus, candPtMinus;
0414
0415 int nite = quad.ConversionCandidate(candVtx, candPtPlus, candPtMinus);
0416
0417 if (!(nite && abs(nite) < 25 && nite != -1000 && nite != -2000))
0418 return nullptr;
0419
0420
0421
0422
0423 #ifdef mydebug_sguazz
0424 std::cout << " >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
0425 << "\n";
0426 std::cout << " >>>>>>>>>>> Conv Cand: "
0427 << " Vertex X: " << candVtx.X() << " [cm] Y: " << candVtx.Y() << " [cm] pt+: " << candPtPlus
0428 << " [GeV] pt-: " << candPtMinus << " [GeV]; #its: " << nite << "\n";
0429 #endif
0430
0431
0432 double minLegPt = CleaningMinLegPt;
0433 double maxRadialDistance = CleaningmaxRadialDistance;
0434
0435
0436
0437 if (candPtPlus < minLegPt)
0438 return nullptr;
0439 if (candPtMinus < minLegPt)
0440 return nullptr;
0441
0442 if (candPtPlus > maxLegPt)
0443 return nullptr;
0444 if (candPtMinus > maxLegPt)
0445 return nullptr;
0446
0447
0448 double cr = std::sqrt(candVtx.Perp2());
0449 double maxr2 = (maxRadialDistance + cr);
0450 maxr2 *= maxr2;
0451 if (h2.Perp2() > maxr2)
0452 return nullptr;
0453 if (h3.Perp2() > maxr2)
0454 return nullptr;
0455
0456
0457
0458
0459 const auto& bfield = es.getData(theBfieldToken);
0460 float nomField = bfield.nominalValue();
0461
0462 if (ClusterShapeFiltering) {
0463 if (theComparitor)
0464 theComparitor->init(ev, es);
0465
0466 GlobalTrajectoryParameters pkine;
0467 GlobalTrajectoryParameters mkine;
0468
0469
0470 SeedingHitSet::ConstRecHitPointer ptth2 = phits[1];
0471 SeedingHitSet::ConstRecHitPointer mtth1 = mhits[0];
0472 SeedingHitSet::ConstRecHitPointer mtth2 = mhits[1];
0473
0474 GlobalPoint vertexPos(candVtx.x(), candVtx.y(), candVtx.z());
0475
0476 float ptMinReg = 0.1;
0477 GlobalTrackingRegion region(ptMinReg, vertexPos, 0, 0, true);
0478
0479 FastHelix phelix(ptth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0480 pkine = phelix.stateAtVertex();
0481 FastHelix mhelix(mtth2->globalPosition(), mtth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0482 mkine = mhelix.stateAtVertex();
0483
0484 if (theComparitor && !theComparitor->compatible(phits, pkine, phelix)) {
0485 return nullptr;
0486 }
0487 if (theComparitor && !theComparitor->compatible(mhits, mkine, mhelix)) {
0488 return nullptr;
0489 }
0490 }
0491
0492
0493 double quadPhotCotTheta = 0.;
0494 double quadZ0 = 0.;
0495 simpleGetSlope(ptth2, ptth1, mtth1, mtth2, region, quadPhotCotTheta, quadZ0);
0496
0497 double quadPhotPhi = (candVtx - vPhotVertex).Phi();
0498
0499 math::XYZVector fittedPrimaryVertex(vgPhotVertex.x(), vgPhotVertex.y(), quadZ0);
0500
0501 candVtx.SetZ(std::sqrt(candVtx.Perp2()) * quadPhotCotTheta + quadZ0);
0502 GlobalPoint convVtxGlobalPoint(candVtx.X(), candVtx.Y(), candVtx.Z());
0503
0504
0505
0506
0507
0508
0509 Quad thisQuad;
0510 thisQuad.x = candVtx.X();
0511 thisQuad.y = candVtx.Y();
0512 thisQuad.z = candVtx.Z();
0513 thisQuad.ptPlus = candPtPlus;
0514 thisQuad.ptMinus = candPtMinus;
0515 thisQuad.cot = quadPhotCotTheta;
0516 if (similarQuadExist(thisQuad, quadV) && applyArbitration)
0517 return nullptr;
0518
0519
0520
0521
0522 FastHelix helixPlus(
0523 ptth2->globalPosition(), ptth1->globalPosition(), convVtxGlobalPoint, nomField, &bfield, convVtxGlobalPoint);
0524 GlobalTrajectoryParameters kinePlus = helixPlus.stateAtVertex();
0525 kinePlus = GlobalTrajectoryParameters(
0526 convVtxGlobalPoint,
0527 GlobalVector(candPtPlus * cos(quadPhotPhi), candPtPlus * sin(quadPhotPhi), candPtPlus * quadPhotCotTheta),
0528 1,
0529 &kinePlus.magneticField());
0530
0531
0532
0533 FastHelix helixMinus(
0534 mtth2->globalPosition(), mtth1->globalPosition(), convVtxGlobalPoint, nomField, &bfield, convVtxGlobalPoint);
0535 GlobalTrajectoryParameters kineMinus = helixMinus.stateAtVertex();
0536 kineMinus = GlobalTrajectoryParameters(
0537 convVtxGlobalPoint,
0538 GlobalVector(candPtMinus * cos(quadPhotPhi), candPtMinus * sin(quadPhotPhi), candPtMinus * quadPhotCotTheta),
0539 -1,
0540 &kineMinus.magneticField());
0541
0542 float sinThetaPlus = sin(kinePlus.momentum().theta());
0543 float sinThetaMinus = sin(kineMinus.momentum().theta());
0544 float ptmin = region.ptMin();
0545
0546 GlobalVector vertexBounds(region.originRBound(), region.originRBound(), region.originZBound());
0547
0548 CurvilinearTrajectoryError errorPlus = initialError(vertexBounds, ptmin, sinThetaPlus);
0549 CurvilinearTrajectoryError errorMinus = initialError(vertexBounds, ptmin, sinThetaMinus);
0550 FreeTrajectoryState ftsPlus(kinePlus, errorPlus);
0551 FreeTrajectoryState ftsMinus(kineMinus, errorMinus);
0552
0553
0554
0555
0556
0557 #ifdef quadDispLine
0558 double vError = region.originZBound();
0559 if (vError > 15.)
0560 vError = 1.;
0561 std::cout << "QuadDispLine " << vgPhotVertex.x() << " " << vgPhotVertex.y() << " " << vgPhotVertex.z() << " "
0562 << vError << " " << vHit[0].x() << " " << vHit[0].y() << " " << vHit[0].z() << " "
0563 << sqrt(getSqrEffectiveErrorOnZ(ptth2, region)) << " " << vHit[1].x() << " " << vHit[1].y() << " "
0564 << vHit[1].z() << " " << sqrt(getSqrEffectiveErrorOnZ(ptth1, region)) << " " << vHit[2].x() << " "
0565 << vHit[2].y() << " " << vHit[2].z() << " " << sqrt(getSqrEffectiveErrorOnZ(mtth1, region)) << " "
0566 << vHit[3].x() << " " << vHit[3].y() << " " << vHit[3].z() << " "
0567 << sqrt(getSqrEffectiveErrorOnZ(mtth2, region)) << " " << candVtx.X() << " " << candVtx.Y() << " "
0568 << candVtx.Z() << " " << fittedPrimaryVertex.X() << " " << fittedPrimaryVertex.Y() << " "
0569 << fittedPrimaryVertex.Z()
0570 << " "
0571
0572
0573 << nite << " " << chi2 << "\n";
0574 #endif
0575 #ifdef mydebug_sguazz
0576 std::cout << " >>>>> Hit quad fast reco done >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
0577 << "\n";
0578 uint32_t detid;
0579 std::cout << "[SeedForPhotonConversionFromQuadruplets]\n ptth1 ";
0580 detid = ptth1->geographicalId().rawId();
0581
0582 std::cout << " \t " << detid << " " << ptth1->localPosition() << " " << ptth1->globalPosition();
0583 detid = ptth2->geographicalId().rawId();
0584 std::cout << " \n\t ptth2 ";
0585
0586 std::cout << " \t " << detid << " " << ptth2->localPosition() << " " << ptth2->globalPosition() << "\nhelix momentum "
0587 << kinePlus.momentum() << " pt " << kinePlus.momentum().perp() << " radius "
0588 << 1 / kinePlus.transverseCurvature() << " q " << kinePlus.charge();
0589 std::cout << " \n\t mtth1 ";
0590 detid = mtth1->geographicalId().rawId();
0591 std::cout << " \t " << detid << " " << mtth1->localPosition() << " " << mtth1->globalPosition();
0592 std::cout << " \n\t mtth2 ";
0593 detid = mtth2->geographicalId().rawId();
0594
0595 std::cout << " \t " << detid << " " << mtth2->localPosition() << " " << mtth2->globalPosition() << "\nhelix momentum "
0596 << kineMinus.momentum() << " pt " << kineMinus.momentum().perp() << " radius "
0597 << 1 / kineMinus.transverseCurvature() << " q " << kineMinus.charge();
0598 std::cout << "\n <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"
0599 << "\n";
0600 #endif
0601
0602
0603 auto builder = static_cast<TkTransientTrackingRecHitBuilder const*>(&es.getData(theTTRHBuilderToken));
0604 cloner = (*builder).cloner();
0605
0606 bool buildSeedBoolPos = buildSeedBool(seedCollection, phits, ftsPlus, es, applydzCAcut, region, dzcut);
0607 bool buildSeedBoolNeg = buildSeedBool(seedCollection, mhits, ftsMinus, es, applydzCAcut, region, dzcut);
0608
0609 if (buildSeedBoolPos && buildSeedBoolNeg) {
0610 buildSeed(seedCollection, phits, ftsPlus, es, false, region);
0611 buildSeed(seedCollection, mhits, ftsMinus, es, false, region);
0612 }
0613
0614 return nullptr;
0615 }
0616
0617 GlobalTrajectoryParameters SeedForPhotonConversionFromQuadruplets::initialKinematic(const SeedingHitSet& hits,
0618 const GlobalPoint& vertexPos,
0619 const edm::EventSetup& es,
0620 const float cotTheta) const {
0621 GlobalTrajectoryParameters kine;
0622
0623 SeedingHitSet::ConstRecHitPointer tth1 = hits[0];
0624 SeedingHitSet::ConstRecHitPointer tth2 = hits[1];
0625
0626 const auto& bfield = es.getData(theBfieldToken);
0627 float nomField = bfield.nominalValue();
0628 bool isBOFF = (0 == nomField);
0629
0630 FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0631 kine = helix.stateAtVertex();
0632
0633
0634 if (fabs(cotTheta) < cotTheta_Max)
0635 kine = GlobalTrajectoryParameters(
0636 kine.position(),
0637 GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta),
0638 kine.charge(),
0639 &kine.magneticField());
0640 else
0641 kine = GlobalTrajectoryParameters(
0642 kine.position(),
0643 GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta_Max),
0644 kine.charge(),
0645 &kine.magneticField());
0646
0647 #ifdef mydebug_seed
0648 uint32_t detid;
0649 (*pss) << "[SeedForPhotonConversionFromQuadruplets] initialKinematic tth1 ";
0650 detid = tth1->geographicalId().rawId();
0651 po.print(*pss, detid);
0652 (*pss) << " \t " << detid << " " << tth1->localPosition() << " " << tth1->globalPosition();
0653 detid = tth2->geographicalId().rawId();
0654 (*pss) << " \n\t tth2 ";
0655 po.print(*pss, detid);
0656 (*pss) << " \t " << detid << " " << tth2->localPosition() << " " << tth2->globalPosition() << "\nhelix momentum "
0657 << kine.momentum() << " pt " << kine.momentum().perp() << " radius " << 1 / kine.transverseCurvature();
0658 #endif
0659
0660 if (isBOFF && (theBOFFMomentum > 0)) {
0661 kine =
0662 GlobalTrajectoryParameters(kine.position(), kine.momentum().unit() * theBOFFMomentum, kine.charge(), &bfield);
0663 }
0664 return kine;
0665 }
0666
0667 CurvilinearTrajectoryError SeedForPhotonConversionFromQuadruplets::initialError(const GlobalVector& vertexBounds,
0668 float ptMin,
0669 float sinTheta) const {
0670
0671
0672 GlobalError vertexErr(sqr(vertexBounds.x()), 0, sqr(vertexBounds.y()), 0, 0, sqr(vertexBounds.z()));
0673
0674 AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
0675
0676
0677
0678
0679 float sin2th = sqr(sinTheta);
0680 float minC00 = 1.0;
0681 C[0][0] = std::max(sin2th / sqr(ptMin), minC00);
0682 float zErr = vertexErr.czz();
0683 float transverseErr = vertexErr.cxx();
0684 C[3][3] = transverseErr;
0685 C[4][4] = zErr * sin2th + transverseErr * (1 - sin2th);
0686
0687 return CurvilinearTrajectoryError(C);
0688 }
0689
0690 const TrajectorySeed* SeedForPhotonConversionFromQuadruplets::buildSeed(TrajectorySeedCollection& seedCollection,
0691 const SeedingHitSet& hits,
0692 const FreeTrajectoryState& fts,
0693 const edm::EventSetup& es,
0694 bool apply_dzCut,
0695 const TrackingRegion& region) const {
0696
0697 const auto& tracker = es.getData(theTrackerToken);
0698
0699
0700 const Propagator* propagator = &es.getData(thePropagatorToken);
0701
0702
0703 KFUpdator updator;
0704
0705
0706
0707 TrajectoryStateOnSurface updatedState;
0708 edm::OwnVector<TrackingRecHit> seedHits;
0709
0710 const TrackingRecHit* hit = nullptr;
0711 for (unsigned int iHit = 0; iHit < hits.size() && iHit < 2; iHit++) {
0712 hit = hits[iHit];
0713 TrajectoryStateOnSurface state =
0714 (iHit == 0) ? propagator->propagate(fts, tracker.idToDet(hit->geographicalId())->surface())
0715 : propagator->propagate(updatedState, tracker.idToDet(hit->geographicalId())->surface());
0716
0717 SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0718
0719 SeedingHitSet::RecHitPointer newtth = refitHit(tth, state);
0720
0721 updatedState = updator.update(state, *newtth);
0722
0723 seedHits.push_back(newtth);
0724 #ifdef mydebug_seed
0725 uint32_t detid = hit->geographicalId().rawId();
0726 (*pss) << "\n[SeedForPhotonConversionFromQuadruplets] hit " << iHit;
0727 po.print(*pss, detid);
0728 (*pss) << " " << detid << "\t lp " << hit->localPosition() << " tth " << tth->localPosition() << " newtth "
0729 << newtth->localPosition() << " state " << state.globalMomentum().perp();
0730 #endif
0731 }
0732
0733 if (!hit)
0734 return nullptr;
0735
0736 PTrajectoryStateOnDet const& PTraj =
0737 trajectoryStateTransform::persistentState(updatedState, hit->geographicalId().rawId());
0738
0739 seedCollection.push_back(TrajectorySeed(PTraj, seedHits, alongMomentum));
0740 return &seedCollection.back();
0741 }
0742
0743 bool SeedForPhotonConversionFromQuadruplets::buildSeedBool(TrajectorySeedCollection& seedCollection,
0744 const SeedingHitSet& hits,
0745 const FreeTrajectoryState& fts,
0746 const edm::EventSetup& es,
0747 bool apply_dzCut,
0748 const TrackingRegion& region,
0749 double dzcut) const {
0750
0751 const auto& tracker = es.getData(theTrackerToken);
0752
0753
0754 const Propagator* propagator = &es.getData(thePropagatorToken);
0755
0756
0757 KFUpdator updator;
0758
0759
0760
0761 TrajectoryStateOnSurface updatedState;
0762 edm::OwnVector<TrackingRecHit> seedHits;
0763
0764 const TrackingRecHit* hit = nullptr;
0765 for (unsigned int iHit = 0; iHit < hits.size() && iHit < 2; iHit++) {
0766 hit = hits[iHit]->hit();
0767 TrajectoryStateOnSurface state =
0768 (iHit == 0) ? propagator->propagate(fts, tracker.idToDet(hit->geographicalId())->surface())
0769 : propagator->propagate(updatedState, tracker.idToDet(hit->geographicalId())->surface());
0770 if (!state.isValid()) {
0771 return false;
0772 }
0773
0774 SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0775
0776 SeedingHitSet::RecHitPointer newtth = refitHit(tth, state);
0777
0778 if (!checkHit(state, newtth, es)) {
0779 return false;
0780 }
0781
0782 updatedState = updator.update(state, *newtth);
0783 if (!updatedState.isValid()) {
0784 return false;
0785 }
0786
0787 seedHits.push_back(newtth->hit()->clone());
0788 }
0789
0790 if (apply_dzCut) {
0791
0792
0793 double zCA;
0794
0795 math::XYZVector EstMomGam(
0796 updatedState.globalMomentum().x(), updatedState.globalMomentum().y(), updatedState.globalMomentum().z());
0797 math::XYZVector EstPosGam(
0798 updatedState.globalPosition().x(), updatedState.globalPosition().y(), updatedState.globalPosition().z());
0799
0800 double EstMomGamLength =
0801 std::sqrt(EstMomGam.x() * EstMomGam.x() + EstMomGam.y() * EstMomGam.y() + EstMomGam.z() * EstMomGam.z());
0802 math::XYZVector EstMomGamNorm(
0803 EstMomGam.x() / EstMomGamLength, EstMomGam.y() / EstMomGamLength, EstMomGam.z() / EstMomGamLength);
0804
0805
0806
0807 const GlobalPoint EstPosGamGlobalPoint(
0808 updatedState.globalPosition().x(), updatedState.globalPosition().y(), updatedState.globalPosition().z());
0809 const GlobalVector EstMomGamGlobalVector(
0810 updatedState.globalMomentum().x(), updatedState.globalMomentum().y(), updatedState.globalMomentum().z());
0811
0812 const MagneticField* magField = &es.getData(theBfieldToken);
0813 TrackCharge qCharge = 0;
0814
0815 const GlobalTrajectoryParameters myGlobalTrajectoryParameter(
0816 EstPosGamGlobalPoint, EstMomGamGlobalVector, qCharge, magField);
0817
0818 AlgebraicSymMatrix66 aCovarianceMatrix;
0819
0820 for (int i = 0; i < 6; ++i)
0821 for (int j = 0; j < 6; ++j)
0822 aCovarianceMatrix(i, j) = 1e-4;
0823
0824 CartesianTrajectoryError myCartesianError(aCovarianceMatrix);
0825
0826 const FreeTrajectoryState stateForProjectionToBeamLine(myGlobalTrajectoryParameter, myCartesianError);
0827
0828 const GlobalPoint BeamSpotGlobalPoint(0, 0, 0);
0829
0830 const reco::BeamSpot::Point BeamSpotPoint(region.origin().x(), region.origin().y(), region.origin().z());
0831
0832 TSCBLBuilderNoMaterial tscblBuilder;
0833
0834 double CovMatEntry = 0.;
0835 reco::BeamSpot::CovarianceMatrix cov;
0836 for (int i = 0; i < 3; ++i) {
0837 cov(i, i) = CovMatEntry;
0838 }
0839 reco::BeamSpot::BeamType BeamType_ = reco::BeamSpot::Unknown;
0840
0841 reco::BeamSpot myBeamSpot(BeamSpotPoint, 0., 0., 0., 0., cov, BeamType_);
0842
0843 TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine, myBeamSpot);
0844 if (tscbl.isValid() == false) {
0845 zCA = 0;
0846 } else {
0847 GlobalPoint v = tscbl.trackStateAtPCA().position();
0848 zCA = v.z();
0849 }
0850
0851
0852
0853
0854
0855
0856
0857
0858
0859
0860
0861
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873
0874
0875
0876
0877
0878
0879
0880
0881
0882
0883 #ifdef mydebug_knuenz
0884 std::cout << "zCA: " << zCA << std::endl;
0885 #endif
0886
0887 if (std::fabs(zCA) > dzcut)
0888 return false;
0889 }
0890
0891 return true;
0892 }
0893
0894 SeedingHitSet::RecHitPointer SeedForPhotonConversionFromQuadruplets::refitHit(
0895 SeedingHitSet::ConstRecHitPointer hit, const TrajectoryStateOnSurface& state) const {
0896
0897
0898
0899
0900
0901
0902
0903 return (SeedingHitSet::RecHitPointer)(cloner(*hit, state));
0904 }
0905
0906
0907
0908
0909
0910 void SeedForPhotonConversionFromQuadruplets::stupidPrint(std::string s, float* d) {
0911 (*pss) << "\n" << s << "\t";
0912 for (size_t i = 0; i < 2; ++i)
0913 (*pss) << std::setw(60) << d[i] << std::setw(1) << " | ";
0914 }
0915
0916 void SeedForPhotonConversionFromQuadruplets::stupidPrint(std::string s, double* d) {
0917 (*pss) << "\n" << s << "\t";
0918 for (size_t i = 0; i < 2; ++i)
0919 (*pss) << std::setw(60) << d[i] << std::setw(1) << " | ";
0920 }
0921
0922 void SeedForPhotonConversionFromQuadruplets::stupidPrint(const char* s, GlobalPoint* d) {
0923 (*pss) << "\n" << s << "\t";
0924 for (size_t i = 0; i < 2; ++i)
0925 (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << " | ";
0926 }
0927
0928 void SeedForPhotonConversionFromQuadruplets::stupidPrint(const char* s, GlobalPoint* d, int n) {
0929 (*pss) << "\n" << s << "\n";
0930 for (int i = 0; i < n; ++i)
0931 (*pss) << std::setw(20) << d[i] << " r " << d[i].perp() << " phi " << d[i].phi() << "\n";
0932 }
0933
0934 #include "DataFormats/Math/interface/deltaPhi.h"
0935
0936 void SeedForPhotonConversionFromQuadruplets::bubbleSortVsPhi(GlobalPoint arr[], int n, GlobalPoint vtx) {
0937 bool swapped = true;
0938 int j = 0;
0939 GlobalPoint tmp;
0940 while (swapped) {
0941 swapped = false;
0942 j++;
0943 for (int i = 0; i < n - j; i++) {
0944 if (reco::deltaPhi((arr[i] - vtx).barePhi(), (arr[i + 1] - vtx).barePhi()) > 0.) {
0945 tmp = arr[i];
0946 arr[i] = arr[i + 1];
0947 arr[i + 1] = tmp;
0948 swapped = true;
0949 }
0950 }
0951 }
0952 }
0953
0954 void SeedForPhotonConversionFromQuadruplets::bubbleReverseSortVsPhi(GlobalPoint arr[], int n, GlobalPoint vtx) {
0955 bool swapped = true;
0956 int j = 0;
0957 GlobalPoint tmp;
0958 while (swapped) {
0959 swapped = false;
0960 j++;
0961 for (int i = 0; i < n - j; i++) {
0962 if (reco::deltaPhi((arr[i] - vtx).barePhi(), (arr[i + 1] - vtx).barePhi()) < 0.) {
0963 tmp = arr[i];
0964 arr[i] = arr[i + 1];
0965 arr[i + 1] = tmp;
0966 swapped = true;
0967 }
0968 }
0969 }
0970 }
0971
0972 double SeedForPhotonConversionFromQuadruplets::simpleGetSlope(const SeedingHitSet::ConstRecHitPointer& ohit,
0973 const SeedingHitSet::ConstRecHitPointer& nohit,
0974 const SeedingHitSet::ConstRecHitPointer& ihit,
0975 const SeedingHitSet::ConstRecHitPointer& nihit,
0976 const TrackingRegion& region,
0977 double& cotTheta,
0978 double& z0) {
0979 double x[5], y[5], e2y[5];
0980
0981
0982
0983
0984 x[0] = ohit->globalPosition().perp();
0985 y[0] = ohit->globalPosition().z();
0986 e2y[0] = getSqrEffectiveErrorOnZ(ohit, region);
0987
0988 x[1] = nohit->globalPosition().perp();
0989 y[1] = nohit->globalPosition().z();
0990 e2y[1] = getSqrEffectiveErrorOnZ(nohit, region);
0991
0992 x[2] = nihit->globalPosition().perp();
0993 y[2] = nihit->globalPosition().z();
0994 e2y[2] = getSqrEffectiveErrorOnZ(nihit, region);
0995
0996 x[3] = ihit->globalPosition().perp();
0997 y[3] = ihit->globalPosition().z();
0998 e2y[3] = getSqrEffectiveErrorOnZ(ihit, region);
0999
1000
1001 x[4] = region.origin().perp();
1002 y[4] = region.origin().z();
1003 double vError = region.originZBound();
1004 if (vError > 15.)
1005 vError = 1.;
1006 e2y[4] = sqr(vError);
1007
1008 double e2z0;
1009 double chi2 = verySimpleFit(5, x, y, e2y, z0, e2z0, cotTheta);
1010
1011 return chi2;
1012 }
1013
1014 double SeedForPhotonConversionFromQuadruplets::verySimpleFit(
1015 int size, double* ax, double* ay, double* e2y, double& p0, double& e2p0, double& p1) {
1016
1017 return 0;
1018 }
1019
1020 double SeedForPhotonConversionFromQuadruplets::getSqrEffectiveErrorOnZ(const SeedingHitSet::ConstRecHitPointer& hit,
1021 const TrackingRegion& region) {
1022
1023
1024
1025
1026 double sqrProjFactor =
1027 sqr((hit->globalPosition().z() - region.origin().z()) / (hit->globalPosition().perp() - region.origin().perp()));
1028 return (hit->globalPositionError().czz() + sqrProjFactor * hit->globalPositionError().rerr(hit->globalPosition()));
1029 }
1030
1031 bool SeedForPhotonConversionFromQuadruplets::similarQuadExist(Quad& thisQuad, std::vector<Quad>& quadV) {
1032 for (auto const& quad : quadV) {
1033 double dx = thisQuad.x - quad.x;
1034 double dy = thisQuad.y - quad.y;
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052 if (sqrt(dx * dx + dy * dy) < 1. && fabs(thisQuad.ptPlus - quad.ptPlus) < 0.5 * quad.ptPlus &&
1053 fabs(thisQuad.ptMinus - quad.ptMinus) < 0.5 * quad.ptMinus) {
1054
1055 return true;
1056 }
1057 }
1058 quadV.push_back(thisQuad);
1059 return false;
1060 }
1061
1062 double SeedForPhotonConversionFromQuadruplets::DeltaPhiManual(const math::XYZVector& v1, const math::XYZVector& v2) {
1063 double kTWOPI = 2. * kPI_;
1064 double DeltaPhiMan = v1.Phi() - v2.Phi();
1065 while (DeltaPhiMan >= kPI_)
1066 DeltaPhiMan -= kTWOPI;
1067 while (DeltaPhiMan < -kPI_)
1068 DeltaPhiMan += kTWOPI;
1069
1070 return DeltaPhiMan;
1071 }