File indexing completed on 2024-04-06 12:01:02
0001
0002 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0003 #include <array>
0004
0005 BaseParticlePropagator::BaseParticlePropagator() : rCyl(0.), rCyl2(0.), zCyl(0.), bField(0), properDecayTime(1E99) { ; }
0006
0007 BaseParticlePropagator::BaseParticlePropagator(const RawParticle& myPart, double R, double Z, double B, double t)
0008 : particle_(myPart), rCyl(R), rCyl2(R * R), zCyl(Z), bField(B), properDecayTime(t) {
0009 init();
0010 }
0011
0012 BaseParticlePropagator::BaseParticlePropagator(const RawParticle& myPart, double R, double Z, double B)
0013 : particle_(myPart), rCyl(R), rCyl2(R * R), zCyl(Z), bField(B), properDecayTime(1E99) {
0014 init();
0015 }
0016
0017 void BaseParticlePropagator::init() {
0018
0019 success = 0;
0020
0021 firstLoop = true;
0022
0023 fiducial = true;
0024
0025 decayed = false;
0026
0027 properTime = 0.;
0028
0029 propDir = 1;
0030
0031 debug = false;
0032 }
0033
0034 bool BaseParticlePropagator::propagate() {
0035
0036
0037
0038 double rPos2 = particle_.R2();
0039
0040 if (onBarrel(rPos2)) {
0041 success = 1;
0042 return true;
0043 }
0044
0045 if (onEndcap(rPos2)) {
0046 success = 2;
0047 return true;
0048 }
0049
0050
0051
0052 if (fabs(particle_.charge()) < 1E-12 || bField < 1E-12) {
0053
0054
0055
0056 double pT2 = particle_.Perp2();
0057
0058
0059 double b2 = rCyl2 * pT2;
0060 double ac = fabs(particle_.X() * particle_.Py() - particle_.Y() * particle_.Px());
0061 double ac2 = ac * ac;
0062
0063
0064
0065
0066 if (ac2 > b2)
0067 return false;
0068
0069
0070 double delta = std::sqrt(b2 - ac2);
0071 double tplus = -(particle_.X() * particle_.Px() + particle_.Y() * particle_.Py()) + delta;
0072 double tminus = -(particle_.X() * particle_.Px() + particle_.Y() * particle_.Py()) - delta;
0073
0074
0075
0076
0077 double solution = tminus < 0 ? tplus / pT2 : tminus / pT2;
0078 if (solution < 0.)
0079 return false;
0080
0081
0082
0083 double zProp = particle_.Z() + particle_.Pz() * solution;
0084 if (fabs(zProp) > zCyl) {
0085 tplus = (zCyl - particle_.Z()) / particle_.Pz();
0086 tminus = (-zCyl - particle_.Z()) / particle_.Pz();
0087 solution = tminus < 0 ? tplus : tminus;
0088 if (solution < 0.)
0089 return false;
0090 success = 2;
0091 } else {
0092 success = 1;
0093 }
0094
0095
0096
0097
0098 double delTime = propDir * particle_.mass() * solution;
0099 double factor = 1.;
0100 properTime += delTime;
0101 if (properTime > properDecayTime) {
0102 factor = 1. - (properTime - properDecayTime) / delTime;
0103 properTime = properDecayTime;
0104 decayed = true;
0105 }
0106
0107
0108
0109
0110 double xProp = particle_.X() + particle_.Px() * solution * factor;
0111 double yProp = particle_.Y() + particle_.Py() * solution * factor;
0112 zProp = particle_.Z() + particle_.Pz() * solution * factor;
0113 double tProp = particle_.T() + particle_.E() * solution * factor;
0114
0115
0116
0117
0118
0119 if (success == 2 && xProp * xProp + yProp * yProp > rCyl2) {
0120 success = -1;
0121 return true;
0122 }
0123
0124
0125
0126
0127 particle_.setVertex(XYZTLorentzVector(xProp, yProp, zProp, tProp));
0128
0129 return true;
0130 }
0131
0132
0133
0134 else {
0135
0136
0137
0138 double pT = particle_.Pt();
0139 double radius = helixRadius(pT);
0140 double phi0 = helixStartPhi();
0141 double xC = helixCentreX(radius, phi0);
0142 double yC = helixCentreY(radius, phi0);
0143 double dist = helixCentreDistToAxis(xC, yC);
0144
0145
0146
0147 if (fabs(fabs(radius) - dist) > rCyl)
0148 return false;
0149
0150
0151
0152
0153 if (particle_.Z() * particle_.Pz() > zCyl * fabs(particle_.Pz())) {
0154 success = -2;
0155 return true;
0156 }
0157
0158 double pZ = particle_.Pz();
0159 double phiProp, zProp;
0160
0161
0162
0163 double rProp = std::min(fabs(radius) + dist - 0.000001, rCyl);
0164
0165
0166 double sinPhiProp = (rProp * rProp - radius * radius - dist * dist) / (2. * dist * radius);
0167
0168
0169 double deltaPhi = 1E99;
0170
0171
0172 if (1. - fabs(sinPhiProp) < 1E-9) {
0173 double cphi0 = std::cos(phi0);
0174 double sphi0 = std::sin(phi0);
0175 double r0 = (particle_.X() * cphi0 + particle_.Y() * sphi0) / radius;
0176 double q0 = (particle_.X() * sphi0 - particle_.Y() * cphi0) / radius;
0177 double rcyl2 = (rCyl2 - particle_.X() * particle_.X() - particle_.Y() * particle_.Y()) / (radius * radius);
0178 double delta = r0 * r0 + rcyl2 * (1. - q0);
0179
0180
0181
0182 deltaPhi = radius > 0 ? (-r0 + std::sqrt(delta)) / (1. - q0) : (-r0 - std::sqrt(delta)) / (1. - q0);
0183 }
0184
0185
0186 if (fabs(deltaPhi) > 1E-3) {
0187
0188
0189 phiProp = std::asin(sinPhiProp);
0190
0191
0192
0193
0194
0195 phiProp = helixCentrePhi(xC, yC) + (inside(rPos2) || onSurface(rPos2) ? phiProp : M_PI - phiProp);
0196
0197
0198
0199
0200 if (fabs(phiProp - phi0) > 2. * M_PI)
0201 phiProp += (phiProp > phi0 ? -2. * M_PI : +2. * M_PI);
0202
0203
0204
0205
0206
0207 if ((phiProp - phi0) * radius < 0.0)
0208 radius > 0.0 ? phiProp += 2 * M_PI : phiProp -= 2 * M_PI;
0209
0210 } else {
0211
0212 phiProp = phi0 + deltaPhi;
0213 }
0214
0215
0216
0217
0218 zProp = particle_.Z() + (phiProp - phi0) * pZ * radius / pT;
0219 if (fabs(zProp) > zCyl) {
0220 zProp = zCyl * fabs(pZ) / pZ;
0221 phiProp = phi0 + (zProp - particle_.Z()) / radius * pT / pZ;
0222 success = 2;
0223
0224 } else {
0225
0226
0227
0228
0229 if (rProp < rCyl) {
0230 if (firstLoop || fabs(pZ) / pT < 1E-10) {
0231 success = 0;
0232
0233 } else {
0234 zProp = zCyl * fabs(pZ) / pZ;
0235 phiProp = phi0 + (zProp - particle_.Z()) / radius * pT / pZ;
0236 success = 2;
0237 }
0238
0239
0240
0241 } else {
0242 success = 1;
0243 }
0244 }
0245
0246
0247
0248
0249
0250
0251 double delTime = propDir * (phiProp - phi0) * radius * particle_.mass() / pT;
0252 double factor = 1.;
0253 properTime += delTime;
0254 if (properTime > properDecayTime) {
0255 factor = 1. - (properTime - properDecayTime) / delTime;
0256 properTime = properDecayTime;
0257 decayed = true;
0258 }
0259
0260 zProp = particle_.Z() + (zProp - particle_.Z()) * factor;
0261
0262 phiProp = phi0 + (phiProp - phi0) * factor;
0263
0264 double sProp = std::sin(phiProp);
0265 double cProp = std::cos(phiProp);
0266 double xProp = xC + radius * sProp;
0267 double yProp = yC - radius * cProp;
0268 double tProp = particle_.T() + (phiProp - phi0) * radius * particle_.E() / pT;
0269 double pxProp = pT * cProp;
0270 double pyProp = pT * sProp;
0271
0272
0273
0274
0275
0276 if (success == 2 && xProp * xProp + yProp * yProp > rCyl2) {
0277 success = -1;
0278 return true;
0279 }
0280
0281
0282
0283 particle_.setVertex(XYZTLorentzVector(xProp, yProp, zProp, tProp));
0284 particle_.setMomentum(pxProp, pyProp, pZ, particle_.E());
0285
0286
0287 return true;
0288 }
0289 }
0290
0291 bool BaseParticlePropagator::backPropagate() {
0292
0293 particle_.setMomentum(-particle_.Px(), -particle_.Py(), -particle_.Pz(), particle_.E());
0294 particle_.setCharge(-particle_.charge());
0295 propDir = -1;
0296 bool done = propagate();
0297 particle_.setMomentum(-particle_.Px(), -particle_.Py(), -particle_.Pz(), particle_.E());
0298 particle_.setCharge(-particle_.charge());
0299 propDir = +1;
0300
0301 return done;
0302 }
0303
0304 BaseParticlePropagator BaseParticlePropagator::propagated() const {
0305
0306 BaseParticlePropagator myPropPart(*this);
0307
0308 myPropPart.firstLoop = false;
0309
0310 myPropPart.propagate();
0311
0312 if (myPropPart.success < 0)
0313 myPropPart.backPropagate();
0314
0315 return myPropPart;
0316 }
0317
0318 void BaseParticlePropagator::setPropagationConditions(double R, double Z, bool f) {
0319 rCyl = R;
0320 rCyl2 = R * R;
0321 zCyl = Z;
0322 firstLoop = f;
0323 }
0324
0325 bool BaseParticlePropagator::propagateToClosestApproach(double x0, double y0, bool first) {
0326
0327 double pT = particle_.Pt();
0328 double radius = helixRadius(pT);
0329 double phi0 = helixStartPhi();
0330
0331
0332
0333 double xC = helixCentreX(radius, phi0);
0334 double yC = helixCentreY(radius, phi0);
0335 double distz = helixCentreDistToAxis(xC - x0, yC - y0);
0336 double dist0 = helixCentreDistToAxis(xC, yC);
0337
0338
0339
0340
0341 double rz, r0, z;
0342 if (particle_.charge() != 0.0 && bField != 0.0) {
0343 rz = fabs(fabs(radius) - distz) + std::sqrt(x0 * x0 + y0 * y0) + 0.0000001;
0344 r0 = fabs(fabs(radius) - dist0) + 0.0000001;
0345 } else {
0346 rz = fabs(particle_.Px() * (particle_.Y() - y0) - particle_.Py() * (particle_.X() - x0)) / particle_.Pt();
0347 r0 = fabs(particle_.Px() * particle_.Y() - particle_.Py() * particle_.X()) / particle_.Pt();
0348 }
0349
0350 z = 999999.;
0351
0352
0353
0354 setPropagationConditions(rz, z, first);
0355 bool done = backPropagate();
0356
0357
0358 if (!done)
0359 return done;
0360
0361
0362 if (fabs(rz - r0) < 1E-10)
0363 return done;
0364 double dist1 = (particle_.X() - x0) * (particle_.X() - x0) + (particle_.Y() - y0) * (particle_.Y() - y0);
0365
0366
0367 if (dist1 < 1E-10)
0368 return done;
0369
0370
0371 XYZTLorentzVector vertex1 = particle_.vertex();
0372 XYZTLorentzVector momentum1 = particle_.momentum();
0373
0374
0375 setPropagationConditions(r0, z, first);
0376 done = backPropagate();
0377 if (!done)
0378 return done;
0379
0380
0381
0382 setPropagationConditions(rz, z, first);
0383 done = backPropagate();
0384 if (!done)
0385 return done;
0386 double dist2 = (particle_.X() - x0) * (particle_.X() - x0) + (particle_.Y() - y0) * (particle_.Y() - y0);
0387
0388
0389 if (dist2 > dist1) {
0390 particle_.setVertex(vertex1);
0391 particle_.setMomentum(momentum1.X(), momentum1.Y(), momentum1.Z(), momentum1.E());
0392 }
0393
0394
0395 return done;
0396 }
0397
0398 bool BaseParticlePropagator::propagateToEcal(bool first) {
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408 setPropagationConditions(129.0, 303.353, first);
0409 return propagate();
0410 }
0411
0412 bool BaseParticlePropagator::propagateToPreshowerLayer1(bool first) {
0413
0414
0415
0416
0417
0418
0419
0420
0421 setPropagationConditions(129.0, 303.353, first);
0422 bool done = propagate();
0423
0424
0425 if (done && (particle_.R2() > 125.0 * 125.0 || particle_.R2() < 45.0 * 45.0))
0426 success = 0;
0427
0428 return done;
0429 }
0430
0431 bool BaseParticlePropagator::propagateToPreshowerLayer2(bool first) {
0432
0433
0434
0435
0436
0437
0438
0439
0440 setPropagationConditions(129.0, 307.838, first);
0441 bool done = propagate();
0442
0443
0444 if (done && (particle_.R2() > 125.0 * 125.0 || particle_.R2() < 45.0 * 45.0))
0445 success = 0;
0446
0447 return done;
0448 }
0449
0450 bool BaseParticlePropagator::propagateToEcalEntrance(bool first) {
0451
0452
0453
0454
0455
0456
0457 setPropagationConditions(129.0, 320.9, first);
0458 bool done = propagate();
0459
0460
0461
0462
0463 if (done && particle_.cos2ThetaV() > 0.81230 && success == 1) {
0464 setPropagationConditions(152.6, 320.9, first);
0465 done = propagate();
0466 }
0467
0468
0469
0470 if (particle_.cos2ThetaV() > 0.99014)
0471 success = 0;
0472
0473 return done;
0474 }
0475
0476 bool BaseParticlePropagator::propagateToHcalEntrance(bool first) {
0477
0478
0479
0480
0481
0482
0483
0484 setPropagationConditions(177.5, 335.0, first);
0485 propDir = 0;
0486 bool done = propagate();
0487 propDir = 1;
0488
0489
0490 if (done && success == 2) {
0491 setPropagationConditions(300.0, 400.458, first);
0492 propDir = 0;
0493 done = propagate();
0494 propDir = 1;
0495 }
0496
0497
0498
0499 if (done && particle_.cos2ThetaV() > 0.99014)
0500 success = 0;
0501
0502 return done;
0503 }
0504
0505 bool BaseParticlePropagator::propagateToVFcalEntrance(bool first) {
0506
0507
0508
0509
0510 setPropagationConditions(400.0, 1114.95, first);
0511 propDir = 0;
0512 bool done = propagate();
0513 propDir = 1;
0514
0515 if (!done)
0516 success = 0;
0517
0518
0519
0520
0521 double c2teta = particle_.cos2ThetaV();
0522 if (done && (c2teta < 0.99014 || c2teta > 0.9998755))
0523 success = 0;
0524
0525 return done;
0526 }
0527
0528 bool BaseParticlePropagator::propagateToHcalExit(bool first) {
0529
0530
0531
0532
0533
0534
0535
0536 setPropagationConditions(263.9, 554.1, first);
0537
0538 propDir = 0;
0539 bool done = propagate();
0540 propDir = 1;
0541
0542 return done;
0543 }
0544
0545 bool BaseParticlePropagator::propagateToHOLayer(bool first) {
0546
0547
0548
0549
0550
0551
0552
0553 setPropagationConditions(387.6, 1000.0, first);
0554
0555
0556 propDir = 0;
0557 bool done = propagate();
0558 propDir = 1;
0559
0560 if (done && std::abs(particle_.Z()) > 700.25)
0561 success = 0;
0562
0563 return done;
0564 }
0565
0566 bool BaseParticlePropagator::propagateToNominalVertex(const XYZTLorentzVector& v) {
0567
0568 if (particle_.charge() == 0. || bField == 0.)
0569 return false;
0570
0571
0572 double dx = particle_.X() - v.X();
0573 double dy = particle_.Y() - v.Y();
0574 double phi = std::atan2(dy, dx) + std::asin(std::sqrt(dx * dx + dy * dy) / (2. * helixRadius()));
0575
0576
0577 double pT = particle_.pt();
0578
0579
0580 particle_.SetPx(pT * std::cos(phi));
0581 particle_.SetPy(pT * std::sin(phi));
0582
0583 return propagateToClosestApproach(v.X(), v.Y());
0584 }
0585
0586 bool BaseParticlePropagator::propagateToBeamCylinder(const XYZTLorentzVector& v, double radius) {
0587
0588 if (particle_.charge() == 0. || bField == 0.)
0589 return fabs(particle_.Px() * particle_.Y() - particle_.Py() * particle_.X()) / particle_.Pt() < radius;
0590
0591
0592
0593
0594
0595 double r = radius;
0596 double r2 = r * r;
0597 double r4 = r2 * r2;
0598
0599 double dx = particle_.X() - v.X();
0600 double dy = particle_.Y() - v.Y();
0601 double dz = particle_.Z() - v.Z();
0602 double Sxy = particle_.X() * v.X() + particle_.Y() * v.Y();
0603 double Dxy = particle_.Y() * v.X() - particle_.X() * v.Y();
0604 double Dxy2 = Dxy * Dxy;
0605 double Sxy2 = Sxy * Sxy;
0606 double SDxy = dx * dx + dy * dy;
0607 double SSDxy = std::sqrt(SDxy);
0608 double ATxy = std::atan2(dy, dx);
0609
0610 double a = r2 - Dxy2 / SDxy;
0611 double b = r * (r2 - Sxy);
0612 double c = r4 - 2. * Sxy * r2 + Sxy2 + Dxy2;
0613
0614
0615
0616 std::array<double, 4> helixRs = {{0.}};
0617 helixRs[0] = (b - std::sqrt(b * b - a * c)) / (2. * a);
0618 helixRs[1] = (b + std::sqrt(b * b - a * c)) / (2. * a);
0619 helixRs[2] = -helixRs[0];
0620 helixRs[3] = -helixRs[1];
0621
0622 std::array<double, 4> helixPhis = {{0.}};
0623 helixPhis[0] = std::asin(SSDxy / (2. * helixRs[0]));
0624 helixPhis[1] = std::asin(SSDxy / (2. * helixRs[1]));
0625 helixPhis[2] = -helixPhis[0];
0626 helixPhis[3] = -helixPhis[1];
0627
0628 double solution1 = 0.;
0629 double solution2 = 0.;
0630 double phi1 = 0.;
0631 double phi2 = 0.;
0632
0633 for (unsigned int i = 0; i < 4; ++i) {
0634 helixPhis[i] = ATxy + helixPhis[i];
0635
0636 double xC = helixCentreX(helixRs[i], helixPhis[i]);
0637 double yC = helixCentreY(helixRs[i], helixPhis[i]);
0638
0639 double dist = helixCentreDistToAxis(xC, yC);
0640
0641
0642
0643
0644
0645
0646
0647 if (fabs(fabs(fabs(helixRs[i]) - dist) - fabs(radius)) < 1e-6) {
0648
0649 if (solution1 == 0.) {
0650 solution1 = helixRs[i];
0651 phi1 = helixPhis[i];
0652
0653 } else if (solution2 == 0.) {
0654 if (helixRs[i] < solution1) {
0655 solution2 = solution1;
0656 solution1 = helixRs[i];
0657 phi2 = phi1;
0658 phi1 = helixPhis[i];
0659 } else {
0660 solution2 = helixRs[i];
0661 phi2 = helixPhis[i];
0662 }
0663
0664 } else {
0665 std::cout << "warning !!! More than two solutions for this track !!! " << std::endl;
0666 }
0667 }
0668 }
0669
0670
0671 double pT = 0.;
0672 double helixPhi = 0.;
0673 double helixR = 0.;
0674
0675 if (solution1 * solution2 == 0.) {
0676 std::cout << "warning !!! Less than two solution for this track! " << solution1 << " " << solution2 << std::endl;
0677 return false;
0678
0679 } else if (solution1 * solution2 < 0.) {
0680 pT = 1000.;
0681 double norm = pT / SSDxy;
0682 particle_.setCharge(+1.);
0683 particle_.setMomentum(dx * norm, dy * norm, dz * norm, 0.);
0684 particle_.SetE(std::sqrt(particle_.momentum().Vect().Mag2()));
0685
0686 } else {
0687 if (solution1 < 0.) {
0688 helixR = solution1;
0689 helixPhi = phi1;
0690 particle_.setCharge(+1.);
0691 } else {
0692 helixR = solution2;
0693 helixPhi = phi2;
0694 particle_.setCharge(-1.);
0695 }
0696 pT = fabs(helixR) * 1e-5 * c_light() * bField;
0697 double norm = pT / SSDxy;
0698 particle_.setMomentum(pT * std::cos(helixPhi), pT * std::sin(helixPhi), dz * norm, 0.);
0699 particle_.SetE(std::sqrt(particle_.momentum().Vect().Mag2()));
0700 }
0701
0702
0703 return propagateToClosestApproach();
0704 }
0705
0706 double BaseParticlePropagator::xyImpactParameter(double x0, double y0) const {
0707 double ip = 0.;
0708 double pT = particle_.Pt();
0709
0710 if (particle_.charge() != 0.0 && bField != 0.0) {
0711 double radius = helixRadius(pT);
0712 double phi0 = helixStartPhi();
0713
0714
0715 double xC = helixCentreX(radius, phi0);
0716 double yC = helixCentreY(radius, phi0);
0717 double distz = helixCentreDistToAxis(xC - x0, yC - y0);
0718 ip = distz - fabs(radius);
0719 } else {
0720 ip = fabs(particle_.Px() * (particle_.Y() - y0) - particle_.Py() * (particle_.X() - x0)) / pT;
0721 }
0722
0723 return ip;
0724 }