TwoBodyDecay/src/TwoBodyDecayLinearizationPointFinder.cc

0001
0002 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayLinearizationPointFinder.h"
0003 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
0004 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0005
0006 const TwoBodyDecayParameters TwoBodyDecayLinearizationPointFinder::getLinearizationPoint(
0007     const std::vector<RefCountedLinearizedTrackState> &tracks,
0008     const double primaryMass,
0009     const double secondaryMass) const {
0010   GlobalPoint linPoint = tracks[0]->linearizationPoint();
0011   PerigeeLinearizedTrackState *linTrack1 = dynamic_cast<PerigeeLinearizedTrackState *>(tracks[0].get());
0012   PerigeeLinearizedTrackState *linTrack2 = dynamic_cast<PerigeeLinearizedTrackState *>(tracks[1].get());
0013   if (!linTrack1 || !linTrack2)
0014     return TwoBodyDecayParameters();
0015
0016   GlobalVector firstMomentum = linTrack1->predictedState().momentum();
0017   GlobalVector secondMomentum = linTrack2->predictedState().momentum();
0018
0019   AlgebraicVector secondaryMomentum1(3);
0020   secondaryMomentum1[0] = firstMomentum.x();
0021   secondaryMomentum1[1] = firstMomentum.y();
0022   secondaryMomentum1[2] = firstMomentum.z();
0023
0024   AlgebraicVector secondaryMomentum2(3);
0025   secondaryMomentum2[0] = secondMomentum.x();
0026   secondaryMomentum2[1] = secondMomentum.y();
0027   secondaryMomentum2[2] = secondMomentum.z();
0028
0029   AlgebraicVector primaryMomentum = secondaryMomentum1 + secondaryMomentum2;
0030
0031   TwoBodyDecayModel decayModel(primaryMass, secondaryMass);
0032   AlgebraicMatrix rotMat = decayModel.rotationMatrix(primaryMomentum[0], primaryMomentum[1], primaryMomentum[2]);
0033   AlgebraicMatrix invRotMat = rotMat.T();
0034
0035   double p = primaryMomentum.norm();
0036   double pSquared = p * p;
0037   double gamma = sqrt(pSquared + primaryMass * primaryMass) / primaryMass;
0038   double betaGamma = p / primaryMass;
0039   AlgebraicSymMatrix lorentzTransformation(4, 1);
0040   lorentzTransformation[0][0] = gamma;
0041   lorentzTransformation[3][3] = gamma;
0042   lorentzTransformation[0][3] = -betaGamma;
0043
0044   double p1 = secondaryMomentum1.norm();
0045   AlgebraicVector boostedLorentzMomentum1(4);
0046   boostedLorentzMomentum1[0] = sqrt(p1 * p1 + secondaryMass * secondaryMass);
0047   boostedLorentzMomentum1.sub(2, invRotMat * secondaryMomentum1);
0048
0049   AlgebraicVector restFrameLorentzMomentum1 = lorentzTransformation * boostedLorentzMomentum1;
0050   AlgebraicVector restFrameMomentum1 = restFrameLorentzMomentum1.sub(2, 4);
0051   double perp1 = sqrt(restFrameMomentum1[0] * restFrameMomentum1[0] + restFrameMomentum1[1] * restFrameMomentum1[1]);
0052   double theta1 = atan2(perp1, restFrameMomentum1[2]);
0053   double phi1 = atan2(restFrameMomentum1[1], restFrameMomentum1[0]);
0054
0055   double p2 = secondaryMomentum2.norm();
0056   AlgebraicVector boostedLorentzMomentum2(4);
0057   boostedLorentzMomentum2[0] = sqrt(p2 * p2 + secondaryMass * secondaryMass);
0058   boostedLorentzMomentum2.sub(2, invRotMat * secondaryMomentum2);
0059
0060   AlgebraicVector restFrameLorentzMomentum2 = lorentzTransformation * boostedLorentzMomentum2;
0061   AlgebraicVector restFrameMomentum2 = restFrameLorentzMomentum2.sub(2, 4);
0062   double perp2 = sqrt(restFrameMomentum2[0] * restFrameMomentum2[0] + restFrameMomentum2[1] * restFrameMomentum2[1]);
0063   double theta2 = atan2(perp2, restFrameMomentum2[2]);
0064   double phi2 = atan2(restFrameMomentum2[1], restFrameMomentum2[0]);
0065
0066   double pi = 3.141592654;
0067   double relSign = -1.;
0068
0069   if (phi1 < 0)
0070     phi1 += 2 * pi;
0071   if (phi2 < 0)
0072     phi2 += 2 * pi;
0073   if (phi1 > phi2)
0074     relSign = 1.;
0075
0076   double momentumSquared1 = secondaryMomentum1.normsq();
0077   double energy1 = sqrt(secondaryMass * secondaryMass + momentumSquared1);
0078   double momentumSquared2 = secondaryMomentum2.normsq();
0079   double energy2 = sqrt(secondaryMass * secondaryMass + momentumSquared2);
0080   double sumMomentaSquared = (secondaryMomentum1 + secondaryMomentum2).normsq();
0081   double sumEnergiesSquared = (energy1 + energy2) * (energy1 + energy2);
0082   double estimatedPrimaryMass = sqrt(sumEnergiesSquared - sumMomentaSquared);
0083
0084   AlgebraicVector linParam(TwoBodyDecayParameters::dimension, 0);
0085   linParam[TwoBodyDecayParameters::x] = linPoint.x();
0086   linParam[TwoBodyDecayParameters::y] = linPoint.y();
0087   linParam[TwoBodyDecayParameters::z] = linPoint.z();
0088   linParam[TwoBodyDecayParameters::px] = primaryMomentum[0];
0089   linParam[TwoBodyDecayParameters::py] = primaryMomentum[1];
0090   linParam[TwoBodyDecayParameters::pz] = primaryMomentum[2];
0091   linParam[TwoBodyDecayParameters::theta] = 0.5 * (theta1 - theta2 + pi);
0092   linParam[TwoBodyDecayParameters::phi] = 0.5 * (phi1 + phi2 + relSign * pi);
0093   linParam[TwoBodyDecayParameters::mass] = estimatedPrimaryMass;
0094
0095   return TwoBodyDecayParameters(linParam);
0096 }