Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:49:54

0001 //FAMOS headers
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   // The status of the propagation
0019   success = 0;
0020   // Propagate only for the first half-loop
0021   firstLoop = true;
0022   //
0023   fiducial = true;
0024   // The particle has not yet decayed
0025   decayed = false;
0026   // The proper time is set to zero
0027   properTime = 0.;
0028   // The propagation direction is along the momentum
0029   propDir = 1;
0030   //
0031   debug = false;
0032 }
0033 
0034 bool BaseParticlePropagator::propagate() {
0035   //
0036   // Check that the particle is not already on the cylinder surface
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   // Treat neutral particles (or no magnetic field) first
0051   //
0052   if (fabs(particle_.charge()) < 1E-12 || bField < 1E-12) {
0053     //
0054     // First check that the particle crosses the cylinder
0055     //
0056     double pT2 = particle_.Perp2();
0057     //    double pT2 = pT*pT;
0058     //    double b2 = rCyl * pT;
0059     double b2 = rCyl2 * pT2;
0060     double ac = fabs(particle_.X() * particle_.Py() - particle_.Y() * particle_.Px());
0061     double ac2 = ac * ac;
0062     //
0063     // The particle never crosses (or never crossed) the cylinder
0064     //
0065     //    if ( ac > b2 ) return false;
0066     if (ac2 > b2)
0067       return false;
0068 
0069     //    double delta  = std::sqrt(b2*b2 - ac*ac);
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     // Find the first (time-wise) intersection point with the cylinder
0075     //
0076 
0077     double solution = tminus < 0 ? tplus / pT2 : tminus / pT2;
0078     if (solution < 0.)
0079       return false;
0080     //
0081     // Check that the particle does (did) not exit from the endcaps first.
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     // Check the decay time
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     // Compute the coordinates of the RawParticle after propagation
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     // Last check : Although propagated to the endcaps, R could still be
0117     // larger than rCyl because the cylinder was too short...
0118     //
0119     if (success == 2 && xProp * xProp + yProp * yProp > rCyl2) {
0120       success = -1;
0121       return true;
0122     }
0123 
0124     //
0125     // ... and update the particle with its propagated coordinates
0126     //
0127     particle_.setVertex(XYZTLorentzVector(xProp, yProp, zProp, tProp));
0128 
0129     return true;
0130   }
0131   //
0132   // Treat charged particles with magnetic field next.
0133   //
0134   else {
0135     //
0136     // First check that the particle can cross the cylinder
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     // The helix either is outside or includes the cylinder -> no intersection
0146     //
0147     if (fabs(fabs(radius) - dist) > rCyl)
0148       return false;
0149     //
0150     // The particle is already away from the endcaps, and will never come back
0151     // Could be back-propagated, so warn the user that it is so.
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     // Does the helix cross the cylinder barrel ? If not, do the first half-loop
0162     //
0163     double rProp = std::min(fabs(radius) + dist - 0.000001, rCyl);
0164 
0165     // Check for rounding errors in the ArcSin.
0166     double sinPhiProp = (rProp * rProp - radius * radius - dist * dist) / (2. * dist * radius);
0167 
0168     //
0169     double deltaPhi = 1E99;
0170 
0171     // Taylor development up to third order for large momenta
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       // This is a solution of a second order equation, assuming phi = phi0 + epsilon
0181       // and epsilon is small.
0182       deltaPhi = radius > 0 ? (-r0 + std::sqrt(delta)) / (1. - q0) : (-r0 - std::sqrt(delta)) / (1. - q0);
0183     }
0184 
0185     // Use complete calculation otherwise, or if the delta phi is large anyway.
0186     if (fabs(deltaPhi) > 1E-3) {
0187       //    phiProp =  fabs(sinPhiProp) < 0.99999999  ?
0188       //      std::asin(sinPhiProp) : M_PI/2. * sinPhiProp;
0189       phiProp = std::asin(sinPhiProp);
0190 
0191       //
0192       // Compute phi after propagation (two solutions, but the asin definition
0193       // (between -pi/2 and +pi/2) takes care of the wrong one.
0194       //
0195       phiProp = helixCentrePhi(xC, yC) + (inside(rPos2) || onSurface(rPos2) ? phiProp : M_PI - phiProp);
0196 
0197       //
0198       // Solve the obvious two-pi ambiguities: more than one turn!
0199       //
0200       if (fabs(phiProp - phi0) > 2. * M_PI)
0201         phiProp += (phiProp > phi0 ? -2. * M_PI : +2. * M_PI);
0202 
0203       //
0204       // Check that the phi difference has the right sign, according to Q*B
0205       // (Another two pi ambuiguity, slightly more subtle)
0206       //
0207       if ((phiProp - phi0) * radius < 0.0)
0208         radius > 0.0 ? phiProp += 2 * M_PI : phiProp -= 2 * M_PI;
0209 
0210     } else {
0211       // Use Taylor
0212       phiProp = phi0 + deltaPhi;
0213     }
0214 
0215     //
0216     // Check that the particle does not exit from the endcaps first.
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       // If the particle does not cross the barrel, either process the
0227       // first-half loop only or propagate all the way down to the endcaps
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         // The particle crossed the barrel
0240         //
0241       } else {
0242         success = 1;
0243       }
0244     }
0245     //
0246     // Compute the coordinates of the RawParticle after propagation
0247     //
0248     //
0249     // Check the decay time
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     // Last check : Although propagated to the endcaps, R could still be
0274     // larger than rCyl because the cylinder was too short...
0275     //
0276     if (success == 2 && xProp * xProp + yProp * yProp > rCyl2) {
0277       success = -1;
0278       return true;
0279     }
0280     //
0281     // ... and update the particle vertex and momentum
0282     //
0283     particle_.setVertex(XYZTLorentzVector(xProp, yProp, zProp, tProp));
0284     particle_.setMomentum(pxProp, pyProp, pZ, particle_.E());
0285     //    particle_.SetPx(pxProp);
0286     //    particle_.SetPy(pyProp);
0287     return true;
0288   }
0289 }
0290 
0291 bool BaseParticlePropagator::backPropagate() {
0292   // Backpropagate
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   // Copy the input particle in a new instance
0306   BaseParticlePropagator myPropPart(*this);
0307   // Allow for many loops
0308   myPropPart.firstLoop = false;
0309   // Propagate the particle !
0310   myPropPart.propagate();
0311   // Back propagate if the forward propagation was not a success
0312   if (myPropPart.success < 0)
0313     myPropPart.backPropagate();
0314   // Return the new instance
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   // Pre-computed quantities
0327   double pT = particle_.Pt();
0328   double radius = helixRadius(pT);
0329   double phi0 = helixStartPhi();
0330 
0331   // The Helix centre coordinates are computed wrt to the z axis
0332   // Actually the z axis might be centred on 0,0: it's the beam spot position (x0,y0)!
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   // Propagation to point of clostest approach to z axis
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   // Propagate to the first interesection point
0353   // with cylinder of radius sqrt(x0**2+y0**2)
0354   setPropagationConditions(rz, z, first);
0355   bool done = backPropagate();
0356 
0357   // Unsuccessful propagation - should not happen!
0358   if (!done)
0359     return done;
0360 
0361   // The z axis is (0,0) - no need to go further
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   // We are already at closest approach - no need to go further
0367   if (dist1 < 1E-10)
0368     return done;
0369 
0370   // Keep for later if it happens to be the right solution
0371   XYZTLorentzVector vertex1 = particle_.vertex();
0372   XYZTLorentzVector momentum1 = particle_.momentum();
0373 
0374   // Propagate to the distance of closest approach to (0,0)
0375   setPropagationConditions(r0, z, first);
0376   done = backPropagate();
0377   if (!done)
0378     return done;
0379 
0380   // Propagate to the first interesection point
0381   // with cylinder of radius sqrt(x0**2+y0**2)
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   // Keep the good solution.
0389   if (dist2 > dist1) {
0390     particle_.setVertex(vertex1);
0391     particle_.setMomentum(momentum1.X(), momentum1.Y(), momentum1.Z(), momentum1.E());
0392   }
0393 
0394   // Done
0395   return done;
0396 }
0397 
0398 bool BaseParticlePropagator::propagateToEcal(bool first) {
0399   //
0400   // Propagation to Ecal (including preshower) after the
0401   // last Tracker layer
0402   // TODO: include proper geometry
0403   // Geometry taken from CMS ECAL TDR
0404   //
0405   // First propagate to global barrel / endcap cylinder
0406   //  setPropagationConditions(1290. , 3045 , first);
0407   //  New position of the preshower as of CMSSW_1_7_X
0408   setPropagationConditions(129.0, 303.353, first);
0409   return propagate();
0410 }
0411 
0412 bool BaseParticlePropagator::propagateToPreshowerLayer1(bool first) {
0413   //
0414   // Propagation to Preshower Layer 1
0415   // TODO: include proper geometry
0416   // Geometry taken from CMS ECAL TDR
0417   //
0418   // First propagate to global barrel / endcap cylinder
0419   //  setPropagationConditions(1290., 3045 , first);
0420   //  New position of the preshower as of CMSSW_1_7_X
0421   setPropagationConditions(129.0, 303.353, first);
0422   bool done = propagate();
0423 
0424   // Check that were are on the Layer 1
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   // Propagation to Preshower Layer 2
0434   // TODO: include proper geometry
0435   // Geometry taken from CMS ECAL TDR
0436   //
0437   // First propagate to global barrel / endcap cylinder
0438   //  setPropagationConditions(1290. , 3090 , first);
0439   //  New position of the preshower as of CMSSW_1_7_X
0440   setPropagationConditions(129.0, 307.838, first);
0441   bool done = propagate();
0442 
0443   // Check that we are on Layer 2
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   // Propagation to ECAL entrance
0453   // TODO: include proper geometry
0454   // Geometry taken from CMS ECAL TDR
0455   //
0456   // First propagate to global barrel / endcap cylinder
0457   setPropagationConditions(129.0, 320.9, first);
0458   bool done = propagate();
0459 
0460   // Go to endcap cylinder in the "barrel cut corner"
0461   // eta = 1.479 -> cos^2(theta) = 0.81230
0462   //  if ( done && eta > 1.479 && success == 1 ) {
0463   if (done && particle_.cos2ThetaV() > 0.81230 && success == 1) {
0464     setPropagationConditions(152.6, 320.9, first);
0465     done = propagate();
0466   }
0467 
0468   // We are not in the ECAL acceptance
0469   // eta = 3.0 -> cos^2(theta) = 0.99013
0470   if (particle_.cos2ThetaV() > 0.99014)
0471     success = 0;
0472 
0473   return done;
0474 }
0475 
0476 bool BaseParticlePropagator::propagateToHcalEntrance(bool first) {
0477   //
0478   // Propagation to HCAL entrance
0479   // TODO: include proper geometry
0480   // Geometry taken from CMSSW_3_1_X xml (Sunanda)
0481   //
0482 
0483   // First propagate to global barrel / endcap cylinder
0484   setPropagationConditions(177.5, 335.0, first);
0485   propDir = 0;
0486   bool done = propagate();
0487   propDir = 1;
0488 
0489   // If went through the bottom of HB cylinder -> re-propagate to HE surface
0490   if (done && success == 2) {
0491     setPropagationConditions(300.0, 400.458, first);
0492     propDir = 0;
0493     done = propagate();
0494     propDir = 1;
0495   }
0496 
0497   // out of the HB/HE acceptance
0498   // eta = 3.0 -> cos^2(theta) = 0.99014
0499   if (done && particle_.cos2ThetaV() > 0.99014)
0500     success = 0;
0501 
0502   return done;
0503 }
0504 
0505 bool BaseParticlePropagator::propagateToVFcalEntrance(bool first) {
0506   //
0507   // Propagation to VFCAL (HF) entrance
0508   // Geometry taken from xml: Geometry/CMSCommonData/data/cms.xml
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   // We are not in the VFCAL acceptance
0519   // eta = 3.0  -> cos^2(theta) = 0.99014
0520   // eta = 5.2  -> cos^2(theta) = 0.9998755
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   // Propagation to HCAL exit
0531   // TODO: include proper geometry
0532   // Geometry taken from CMSSW_3_1_X xml (Sunanda)
0533   //
0534 
0535   // Approximate it to a single cylinder as it is not that crucial.
0536   setPropagationConditions(263.9, 554.1, first);
0537   //  this->rawPart().particle_.setCharge(0.0); ?? Shower Propagation ??
0538   propDir = 0;
0539   bool done = propagate();
0540   propDir = 1;
0541 
0542   return done;
0543 }
0544 
0545 bool BaseParticlePropagator::propagateToHOLayer(bool first) {
0546   //
0547   // Propagation to Layer0 of HO (Ring-0)
0548   // TODO: include proper geometry
0549   // Geometry taken from CMSSW_3_1_X xml (Sunanda)
0550   //
0551 
0552   // Approximate it to a single cylinder as it is not that crucial.
0553   setPropagationConditions(387.6, 1000.0, first);  //Dia is the middle of HO \sqrt{384.8**2+46.2**2}
0554   //  setPropagationConditions(410.2, 1000.0, first); //sqrt{406.6**2+54.2**2} //for outer layer
0555   //  this->rawPart().setCharge(0.0); ?? Shower Propagation ??
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   // Not implemented for neutrals (used for electrons only)
0568   if (particle_.charge() == 0. || bField == 0.)
0569     return false;
0570 
0571   // Define the proper pT direction to meet the point (vx,vy) in (x,y)
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   // The absolute pT (kept to the original value)
0577   double pT = particle_.pt();
0578 
0579   // Set the new pT
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   // For neutral (or BField = 0, simply check that the track passes through the cylinder
0588   if (particle_.charge() == 0. || bField == 0.)
0589     return fabs(particle_.Px() * particle_.Y() - particle_.Py() * particle_.X()) / particle_.Pt() < radius;
0590 
0591   // Now go to the charged particles
0592 
0593   // A few needed intermediate variables (to make the code understandable)
0594   // The inner cylinder
0595   double r = radius;
0596   double r2 = r * r;
0597   double r4 = r2 * r2;
0598   // The two hits
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   // The coefficients of the second order equation for the trajectory radius
0610   double a = r2 - Dxy2 / SDxy;
0611   double b = r * (r2 - Sxy);
0612   double c = r4 - 2. * Sxy * r2 + Sxy2 + Dxy2;
0613 
0614   // Here are the four possible solutions for
0615   // 1) The trajectory radius
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   // 2) The azimuthal direction at the second point
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   // Only two solutions are valid though
0628   double solution1 = 0.;
0629   double solution2 = 0.;
0630   double phi1 = 0.;
0631   double phi2 = 0.;
0632   // Loop over the four possibilities
0633   for (unsigned int i = 0; i < 4; ++i) {
0634     helixPhis[i] = ATxy + helixPhis[i];
0635     // The centre of the helix
0636     double xC = helixCentreX(helixRs[i], helixPhis[i]);
0637     double yC = helixCentreY(helixRs[i], helixPhis[i]);
0638     // The distance of that centre to the origin
0639     double dist = helixCentreDistToAxis(xC, yC);
0640     /*
0641     std::cout << "Solution "<< i << " = " << helixRs[i] << " accepted ? " 
0642           << fabs(fabs(helixRs[i]) - dist ) << " - " << radius << " = " 
0643           << fabs(fabs(helixRs[i]) - dist ) - fabs(radius) << " " 
0644           << ( fabs( fabs(fabs(helixRs[i]) - dist) -fabs(radius) ) < 1e-6)  << std::endl;
0645     */
0646     // Check that the propagation will lead to a trajectroy tangent to the inner cylinder
0647     if (fabs(fabs(fabs(helixRs[i]) - dist) - fabs(radius)) < 1e-6) {
0648       // Fill first solution
0649       if (solution1 == 0.) {
0650         solution1 = helixRs[i];
0651         phi1 = helixPhis[i];
0652         // Fill second solution (ordered)
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         // Must not happen!
0664       } else {
0665         std::cout << "warning !!! More than two solutions for this track !!! " << std::endl;
0666       }
0667     }
0668   }
0669 
0670   // Find the largest possible pT compatible with coming from within the inne cylinder
0671   double pT = 0.;
0672   double helixPhi = 0.;
0673   double helixR = 0.;
0674   // This should not happen
0675   if (solution1 * solution2 == 0.) {
0676     std::cout << "warning !!! Less than two solution for this track! " << solution1 << " " << solution2 << std::endl;
0677     return false;
0678     // If the two solutions have opposite sign, it means that pT = infinity is a possibility.
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     // Otherwise take the solution that gives the largest transverse momentum
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   // Propagate to closest approach to get the Z value (a bit of an overkill)
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     // The Helix centre coordinates are computed wrt to the z axis
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 }