| 1 |
|
| 2 |
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayLinearizationPointFinder.h"
|
| 3 |
#include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
|
| 4 |
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
|
| 5 |
|
| 6 |
const TwoBodyDecayParameters
|
| 7 |
TwoBodyDecayLinearizationPointFinder::getLinearizationPoint( const std::vector< RefCountedLinearizedTrackState > & tracks,
|
| 8 |
const double primaryMass,
|
| 9 |
const double secondaryMass ) const
|
| 10 |
{
|
| 11 |
GlobalPoint linPoint = tracks[0]->linearizationPoint();
|
| 12 |
PerigeeLinearizedTrackState* linTrack1 = dynamic_cast<PerigeeLinearizedTrackState*>( tracks[0].get() );
|
| 13 |
GlobalVector firstMomentum = linTrack1->predictedState().momentum();
|
| 14 |
PerigeeLinearizedTrackState* linTrack2 = dynamic_cast<PerigeeLinearizedTrackState*>( tracks[1].get() );
|
| 15 |
GlobalVector secondMomentum = linTrack2->predictedState().momentum();
|
| 16 |
|
| 17 |
AlgebraicVector secondaryMomentum1( 3 );
|
| 18 |
secondaryMomentum1[0] = firstMomentum.x();
|
| 19 |
secondaryMomentum1[1] = firstMomentum.y();
|
| 20 |
secondaryMomentum1[2] = firstMomentum.z();
|
| 21 |
|
| 22 |
AlgebraicVector secondaryMomentum2( 3 );
|
| 23 |
secondaryMomentum2[0] = secondMomentum.x();
|
| 24 |
secondaryMomentum2[1] = secondMomentum.y();
|
| 25 |
secondaryMomentum2[2] = secondMomentum.z();
|
| 26 |
|
| 27 |
AlgebraicVector primaryMomentum = secondaryMomentum1 + secondaryMomentum2;
|
| 28 |
|
| 29 |
TwoBodyDecayModel decayModel( primaryMass, secondaryMass );
|
| 30 |
AlgebraicMatrix rotMat = decayModel.rotationMatrix( primaryMomentum[0], primaryMomentum[1], primaryMomentum[2] );
|
| 31 |
AlgebraicMatrix invRotMat = rotMat.T();
|
| 32 |
|
| 33 |
double p = primaryMomentum.norm();
|
| 34 |
double pSquared = p*p;
|
| 35 |
double gamma = sqrt( pSquared + primaryMass*primaryMass )/primaryMass;
|
| 36 |
double betaGamma = p/primaryMass;
|
| 37 |
AlgebraicSymMatrix lorentzTransformation( 4, 1 );
|
| 38 |
lorentzTransformation[0][0] = gamma;
|
| 39 |
lorentzTransformation[3][3] = gamma;
|
| 40 |
lorentzTransformation[0][3] = -betaGamma;
|
| 41 |
|
| 42 |
double p1 = secondaryMomentum1.norm();
|
| 43 |
AlgebraicVector boostedLorentzMomentum1( 4 );
|
| 44 |
boostedLorentzMomentum1[0] = sqrt( p1*p1 + secondaryMass*secondaryMass );
|
| 45 |
boostedLorentzMomentum1.sub( 2, invRotMat*secondaryMomentum1 );
|
| 46 |
|
| 47 |
AlgebraicVector restFrameLorentzMomentum1 = lorentzTransformation*boostedLorentzMomentum1;
|
| 48 |
AlgebraicVector restFrameMomentum1 = restFrameLorentzMomentum1.sub( 2, 4 );
|
| 49 |
double perp1 = sqrt( restFrameMomentum1[0]*restFrameMomentum1[0] + restFrameMomentum1[1]*restFrameMomentum1[1] );
|
| 50 |
double theta1 = atan2( perp1, restFrameMomentum1[2] );
|
| 51 |
double phi1 = atan2( restFrameMomentum1[1], restFrameMomentum1[0] );
|
| 52 |
|
| 53 |
double p2 = secondaryMomentum2.norm();
|
| 54 |
AlgebraicVector boostedLorentzMomentum2( 4 );
|
| 55 |
boostedLorentzMomentum2[0] = sqrt( p2*p2 + secondaryMass*secondaryMass );
|
| 56 |
boostedLorentzMomentum2.sub( 2, invRotMat*secondaryMomentum2 );
|
| 57 |
|
| 58 |
AlgebraicVector restFrameLorentzMomentum2 = lorentzTransformation*boostedLorentzMomentum2;
|
| 59 |
AlgebraicVector restFrameMomentum2 = restFrameLorentzMomentum2.sub( 2, 4 );
|
| 60 |
double perp2 = sqrt( restFrameMomentum2[0]*restFrameMomentum2[0] + restFrameMomentum2[1]*restFrameMomentum2[1] );
|
| 61 |
double theta2 = atan2( perp2, restFrameMomentum2[2] );
|
| 62 |
double phi2 = atan2( restFrameMomentum2[1], restFrameMomentum2[0] );
|
| 63 |
|
| 64 |
double pi = 3.141592654;
|
| 65 |
double relSign = -1.;
|
| 66 |
|
| 67 |
if ( phi1 < 0 ) phi1 += 2*pi;
|
| 68 |
if ( phi2 < 0 ) phi2 += 2*pi;
|
| 69 |
if ( phi1 > phi2 ) relSign = 1.;
|
| 70 |
|
| 71 |
double momentumSquared1 = secondaryMomentum1.normsq();
|
| 72 |
double energy1 = sqrt( secondaryMass*secondaryMass + momentumSquared1 );
|
| 73 |
double momentumSquared2 = secondaryMomentum2.normsq();
|
| 74 |
double energy2 = sqrt( secondaryMass*secondaryMass + momentumSquared2 );
|
| 75 |
double sumMomentaSquared = ( secondaryMomentum1 + secondaryMomentum2 ).normsq();
|
| 76 |
double sumEnergiesSquared = ( energy1 + energy2 )*( energy1 + energy2 );
|
| 77 |
double estimatedPrimaryMass = sqrt( sumEnergiesSquared - sumMomentaSquared );
|
| 78 |
|
| 79 |
AlgebraicVector linParam( TwoBodyDecayParameters::dimension, 0 );
|
| 80 |
linParam[TwoBodyDecayParameters::x] = linPoint.x();
|
| 81 |
linParam[TwoBodyDecayParameters::y] = linPoint.y();
|
| 82 |
linParam[TwoBodyDecayParameters::z] = linPoint.z();
|
| 83 |
linParam[TwoBodyDecayParameters::px] = primaryMomentum[0];
|
| 84 |
linParam[TwoBodyDecayParameters::py] = primaryMomentum[1];
|
| 85 |
linParam[TwoBodyDecayParameters::pz] = primaryMomentum[2];
|
| 86 |
linParam[TwoBodyDecayParameters::theta] = 0.5*( theta1 - theta2 + pi ) ;
|
| 87 |
linParam[TwoBodyDecayParameters::phi] = 0.5*( phi1 + phi2 + relSign*pi );
|
| 88 |
linParam[TwoBodyDecayParameters::mass] = estimatedPrimaryMass;
|
| 89 |
|
| 90 |
return TwoBodyDecayParameters( linParam );
|
| 91 |
}
|
