Go to the documentation of this file.00001
00002 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayLinearizationPointFinder.h"
00003 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
00004 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"
00005
00006 const TwoBodyDecayParameters
00007 TwoBodyDecayLinearizationPointFinder::getLinearizationPoint( const std::vector< RefCountedLinearizedTrackState > & tracks,
00008 const double primaryMass,
00009 const double secondaryMass ) const
00010 {
00011 GlobalPoint linPoint = tracks[0]->linearizationPoint();
00012 PerigeeLinearizedTrackState* linTrack1 = dynamic_cast<PerigeeLinearizedTrackState*>( tracks[0].get() );
00013 PerigeeLinearizedTrackState* linTrack2 = dynamic_cast<PerigeeLinearizedTrackState*>( tracks[1].get() );
00014 if (!linTrack1 || !linTrack2) return TwoBodyDecayParameters();
00015
00016 GlobalVector firstMomentum = linTrack1->predictedState().momentum();
00017 GlobalVector secondMomentum = linTrack2->predictedState().momentum();
00018
00019 AlgebraicVector secondaryMomentum1( 3 );
00020 secondaryMomentum1[0] = firstMomentum.x();
00021 secondaryMomentum1[1] = firstMomentum.y();
00022 secondaryMomentum1[2] = firstMomentum.z();
00023
00024 AlgebraicVector secondaryMomentum2( 3 );
00025 secondaryMomentum2[0] = secondMomentum.x();
00026 secondaryMomentum2[1] = secondMomentum.y();
00027 secondaryMomentum2[2] = secondMomentum.z();
00028
00029 AlgebraicVector primaryMomentum = secondaryMomentum1 + secondaryMomentum2;
00030
00031 TwoBodyDecayModel decayModel( primaryMass, secondaryMass );
00032 AlgebraicMatrix rotMat = decayModel.rotationMatrix( primaryMomentum[0], primaryMomentum[1], primaryMomentum[2] );
00033 AlgebraicMatrix invRotMat = rotMat.T();
00034
00035 double p = primaryMomentum.norm();
00036 double pSquared = p*p;
00037 double gamma = sqrt( pSquared + primaryMass*primaryMass )/primaryMass;
00038 double betaGamma = p/primaryMass;
00039 AlgebraicSymMatrix lorentzTransformation( 4, 1 );
00040 lorentzTransformation[0][0] = gamma;
00041 lorentzTransformation[3][3] = gamma;
00042 lorentzTransformation[0][3] = -betaGamma;
00043
00044 double p1 = secondaryMomentum1.norm();
00045 AlgebraicVector boostedLorentzMomentum1( 4 );
00046 boostedLorentzMomentum1[0] = sqrt( p1*p1 + secondaryMass*secondaryMass );
00047 boostedLorentzMomentum1.sub( 2, invRotMat*secondaryMomentum1 );
00048
00049 AlgebraicVector restFrameLorentzMomentum1 = lorentzTransformation*boostedLorentzMomentum1;
00050 AlgebraicVector restFrameMomentum1 = restFrameLorentzMomentum1.sub( 2, 4 );
00051 double perp1 = sqrt( restFrameMomentum1[0]*restFrameMomentum1[0] + restFrameMomentum1[1]*restFrameMomentum1[1] );
00052 double theta1 = atan2( perp1, restFrameMomentum1[2] );
00053 double phi1 = atan2( restFrameMomentum1[1], restFrameMomentum1[0] );
00054
00055 double p2 = secondaryMomentum2.norm();
00056 AlgebraicVector boostedLorentzMomentum2( 4 );
00057 boostedLorentzMomentum2[0] = sqrt( p2*p2 + secondaryMass*secondaryMass );
00058 boostedLorentzMomentum2.sub( 2, invRotMat*secondaryMomentum2 );
00059
00060 AlgebraicVector restFrameLorentzMomentum2 = lorentzTransformation*boostedLorentzMomentum2;
00061 AlgebraicVector restFrameMomentum2 = restFrameLorentzMomentum2.sub( 2, 4 );
00062 double perp2 = sqrt( restFrameMomentum2[0]*restFrameMomentum2[0] + restFrameMomentum2[1]*restFrameMomentum2[1] );
00063 double theta2 = atan2( perp2, restFrameMomentum2[2] );
00064 double phi2 = atan2( restFrameMomentum2[1], restFrameMomentum2[0] );
00065
00066 double pi = 3.141592654;
00067 double relSign = -1.;
00068
00069 if ( phi1 < 0 ) phi1 += 2*pi;
00070 if ( phi2 < 0 ) phi2 += 2*pi;
00071 if ( phi1 > phi2 ) relSign = 1.;
00072
00073 double momentumSquared1 = secondaryMomentum1.normsq();
00074 double energy1 = sqrt( secondaryMass*secondaryMass + momentumSquared1 );
00075 double momentumSquared2 = secondaryMomentum2.normsq();
00076 double energy2 = sqrt( secondaryMass*secondaryMass + momentumSquared2 );
00077 double sumMomentaSquared = ( secondaryMomentum1 + secondaryMomentum2 ).normsq();
00078 double sumEnergiesSquared = ( energy1 + energy2 )*( energy1 + energy2 );
00079 double estimatedPrimaryMass = sqrt( sumEnergiesSquared - sumMomentaSquared );
00080
00081 AlgebraicVector linParam( TwoBodyDecayParameters::dimension, 0 );
00082 linParam[TwoBodyDecayParameters::x] = linPoint.x();
00083 linParam[TwoBodyDecayParameters::y] = linPoint.y();
00084 linParam[TwoBodyDecayParameters::z] = linPoint.z();
00085 linParam[TwoBodyDecayParameters::px] = primaryMomentum[0];
00086 linParam[TwoBodyDecayParameters::py] = primaryMomentum[1];
00087 linParam[TwoBodyDecayParameters::pz] = primaryMomentum[2];
00088 linParam[TwoBodyDecayParameters::theta] = 0.5*( theta1 - theta2 + pi ) ;
00089 linParam[TwoBodyDecayParameters::phi] = 0.5*( phi1 + phi2 + relSign*pi );
00090 linParam[TwoBodyDecayParameters::mass] = estimatedPrimaryMass;
00091
00092 return TwoBodyDecayParameters( linParam );
00093 }
1.7.3