CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch2/src/Alignment/ReferenceTrajectories/src/TwoBodyDecayTrajectory.cc

Go to the documentation of this file.
00001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00002 
00003 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h" 
00004 #include "DataFormats/GeometrySurface/interface/Surface.h" 
00005 #include "Alignment/ReferenceTrajectories/interface/TwoBodyDecayTrajectory.h"
00006 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h" 
00007 #include "DataFormats/Math/interface/Error.h" 
00008 #include "Alignment/TwoBodyDecay/interface/TwoBodyDecayParameters.h"
00009 
00010 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00011 
00012 
00013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00014 
00015 TwoBodyDecayTrajectory::TwoBodyDecayTrajectory( const TwoBodyDecayTrajectoryState& trajectoryState,
00016                                                 const ConstRecHitCollection & recHits,
00017                                                 const MagneticField* magField,
00018                                                 MaterialEffects materialEffects,
00019                                                 PropagationDirection propDir,
00020                                                 bool hitsAreReverse,
00021                                                 const reco::BeamSpot &beamSpot,
00022                                                 bool useRefittedState,
00023                                                 bool constructTsosWithErrors )
00024 
00025   : ReferenceTrajectoryBase( 
00026      TwoBodyDecayParameters::dimension, recHits.first.size() + recHits.second.size(),
00027   (materialEffects >= breakPoints) ? 2*(recHits.first.size() + recHits.second.size())-4 : 0,
00028   (materialEffects >= breakPoints) ? 2*(recHits.first.size() + recHits.second.size())-3 : 1 )
00029 {
00030   if ( hitsAreReverse )
00031   {
00032     TransientTrackingRecHit::ConstRecHitContainer::const_reverse_iterator itRecHits;
00033     ConstRecHitCollection fwdRecHits;
00034 
00035     fwdRecHits.first.reserve( recHits.first.size() );
00036     for ( itRecHits = recHits.first.rbegin(); itRecHits != recHits.first.rend(); ++itRecHits )
00037     {
00038       fwdRecHits.first.push_back( *itRecHits );
00039     }
00040 
00041     fwdRecHits.second.reserve( recHits.second.size() );
00042     for ( itRecHits = recHits.second.rbegin(); itRecHits != recHits.second.rend(); ++itRecHits )
00043     {
00044       fwdRecHits.second.push_back( *itRecHits );
00045     }
00046 
00047     theValidityFlag = this->construct( trajectoryState, fwdRecHits, magField, materialEffects, propDir,
00048                                        beamSpot, useRefittedState, constructTsosWithErrors );
00049   }
00050   else
00051   {
00052     theValidityFlag = this->construct( trajectoryState, recHits, magField, materialEffects, propDir,
00053                                        beamSpot, useRefittedState, constructTsosWithErrors );
00054   }
00055 }
00056 
00057 
00058 TwoBodyDecayTrajectory::TwoBodyDecayTrajectory( void )
00059   : ReferenceTrajectoryBase( 0, 0, 0, 0)
00060 {}
00061 
00062 
00063 bool TwoBodyDecayTrajectory::construct( const TwoBodyDecayTrajectoryState& state,
00064                                         const ConstRecHitCollection& recHits,
00065                                         const MagneticField* field,
00066                                         MaterialEffects materialEffects,
00067                                         PropagationDirection propDir,
00068                                         const reco::BeamSpot &beamSpot,
00069                                         bool useRefittedState,
00070                                         bool constructTsosWithErrors )
00071 {  
00072   const TwoBodyDecayTrajectoryState::TsosContainer& tsos = state.trajectoryStates( useRefittedState );
00073   const TwoBodyDecayTrajectoryState::Derivatives& deriv = state.derivatives();
00074   double mass = state.particleMass();
00075 
00076   //
00077   // first track
00078   //
00079 
00080   // construct a trajectory (hits should be already in correct order)
00081   ReferenceTrajectory trajectory1( tsos.first, recHits.first, false, field, materialEffects,
00082                                    propDir, mass, false, beamSpot);
00083 
00084   // check if construction of trajectory was successful
00085   if ( !trajectory1.isValid() ) return false;
00086   
00087   int nLocal = deriv.first.num_row();
00088   int nTbd   = deriv.first.num_col();
00089   unsigned int nHitMeas1     = trajectory1.numberOfHitMeas();
00090   unsigned int nVirtualMeas1 = trajectory1.numberOfVirtualMeas(); 
00091   unsigned int nPar1         = trajectory1.numberOfPar();
00092   unsigned int nVirtualPar1  = trajectory1.numberOfVirtualPar(); 
00093      
00094   // derivatives of the trajectory w.r.t. to the decay parameters
00095   AlgebraicMatrix fullDeriv1 = trajectory1.derivatives().sub(1,nHitMeas1+nVirtualMeas1,1,nLocal) * trajectory1.localToTrajectory() * deriv.first;
00096 
00097   //
00098   // second track
00099   //
00100 
00101   ReferenceTrajectory trajectory2( tsos.second, recHits.second, false, field, materialEffects,
00102                                    propDir, mass, false, beamSpot );
00103 
00104   if ( !trajectory2.isValid() ) return false;
00105   
00106   unsigned int nHitMeas2     = trajectory2.numberOfHitMeas();
00107   unsigned int nVirtualMeas2 = trajectory2.numberOfVirtualMeas();  
00108   unsigned int nPar2         = trajectory2.numberOfPar();
00109   unsigned int nVirtualPar2  = trajectory2.numberOfVirtualPar();
00110 
00111   AlgebraicMatrix fullDeriv2 = trajectory2.derivatives().sub(1,nHitMeas2+nVirtualMeas2,1,nLocal) * trajectory2.localToTrajectory() * deriv.second;
00112 
00113   //
00114   // combine both tracks
00115   //
00116   
00117   theNumberOfRecHits.first = recHits.first.size();
00118   theNumberOfRecHits.second = recHits.second.size();
00119 
00120   theNumberOfHits = trajectory1.numberOfHits() + trajectory2.numberOfHits(); 
00121   theNumberOfPars = nPar1 + nPar2;
00122   theNumberOfVirtualPars = nVirtualPar1 + nVirtualPar2;
00123   theNumberOfVirtualMeas = nVirtualMeas1 + nVirtualMeas2 + 1; // add virtual mass measurement
00124   
00125   // hit measurements from trajectory 1
00126   int rowOffset = 1;
00127   int colOffset = 1; 
00128   theDerivatives.sub( rowOffset, colOffset,                fullDeriv1.sub(            1, nHitMeas1,                     1, nTbd ) );
00129   colOffset += nTbd;
00130   theDerivatives.sub( rowOffset, colOffset, trajectory1.derivatives().sub(            1, nHitMeas1,            nLocal + 1, nPar1 + nVirtualPar1 ) );
00131   // hit measurements from trajectory 2
00132   rowOffset += nHitMeas1;
00133   colOffset = 1; 
00134   theDerivatives.sub( rowOffset, colOffset,                fullDeriv2.sub(            1, nHitMeas2,                     1, nTbd ) );
00135   colOffset += (nPar1 + nVirtualPar1 + nTbd - nLocal);
00136   theDerivatives.sub( rowOffset, colOffset, trajectory2.derivatives().sub(            1, nHitMeas2,            nLocal + 1, nPar2 + nVirtualPar2 ) );  
00137   // MS measurements from trajectory 1
00138   rowOffset += nHitMeas2;  
00139   colOffset = 1; 
00140   theDerivatives.sub( rowOffset, colOffset,                fullDeriv1.sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1,          1, nTbd ) );  
00141   colOffset += nTbd;
00142   theDerivatives.sub( rowOffset, colOffset, trajectory1.derivatives().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1, nLocal + 1, nPar1 + nVirtualPar1 ) ); 
00143   // MS measurements from trajectory 2
00144   rowOffset += nVirtualMeas1;  
00145   colOffset = 1; 
00146   theDerivatives.sub( rowOffset, colOffset,                fullDeriv2.sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2,          1, nTbd ) );
00147   colOffset += (nPar1 + nVirtualPar1 + nTbd - nLocal);  
00148   theDerivatives.sub( rowOffset, colOffset, trajectory2.derivatives().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2, nLocal + 1, nPar2 + nVirtualPar2 ) ); 
00149         
00150   theMeasurements.sub(                                         1, trajectory1.measurements().sub(            1, nHitMeas1 ) );
00151   theMeasurements.sub( nHitMeas1                             + 1, trajectory2.measurements().sub(            1, nHitMeas2 ) );
00152   theMeasurements.sub( nHitMeas1 + nHitMeas2                 + 1, trajectory1.measurements().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1 ) );
00153   theMeasurements.sub( nHitMeas1 + nHitMeas2 + nVirtualMeas1 + 1, trajectory2.measurements().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2 ) );
00154 
00155   theMeasurementsCov.sub(                                         1, trajectory1.measurementErrors().sub(            1, nHitMeas1 ) );
00156   theMeasurementsCov.sub( nHitMeas1                             + 1, trajectory2.measurementErrors().sub(            1, nHitMeas2 ) );
00157   theMeasurementsCov.sub( nHitMeas1 + nHitMeas2                 + 1, trajectory1.measurementErrors().sub(nHitMeas1 + 1, nHitMeas1 + nVirtualMeas1 ) );
00158   theMeasurementsCov.sub( nHitMeas1 + nHitMeas2 + nVirtualMeas1 + 1, trajectory2.measurementErrors().sub(nHitMeas2 + 1, nHitMeas2 + nVirtualMeas2 ) );
00159 
00160   theTrajectoryPositions.sub(             1, trajectory1.trajectoryPositions() );
00161   theTrajectoryPositions.sub( nHitMeas1 + 1, trajectory2.trajectoryPositions() );
00162 
00163   theTrajectoryPositionCov = state.decayParameters().covariance().similarity( theDerivatives.sub(1, nHitMeas1 + nHitMeas2, 1, 9) );
00164 
00165   theParameters = state.decayParameters().parameters();
00166 
00167   theRecHits.insert( theRecHits.end(), recHits.first.begin(), recHits.first.end() );
00168   theRecHits.insert( theRecHits.end(), recHits.second.begin(), recHits.second.end() );
00169 
00170   // add virtual mass measurement
00171   rowOffset += nVirtualMeas2;
00172   int indMass = rowOffset-1;
00173   theMeasurements[indMass] = state.primaryMass() - state.decayParameters()[TwoBodyDecayParameters::mass];
00174   theMeasurementsCov[indMass][indMass] = state.primaryWidth() * state.primaryWidth();
00175   theDerivatives[indMass][TwoBodyDecayParameters::mass] = 1.0;
00176 
00177   if ( constructTsosWithErrors )
00178   {
00179     constructTsosVecWithErrors( trajectory1, trajectory2, field );
00180   }
00181   else
00182   {
00183     theTsosVec.insert( theTsosVec.end(),
00184                        trajectory1.trajectoryStates().begin(),
00185                        trajectory1.trajectoryStates().end() );
00186 
00187     theTsosVec.insert( theTsosVec.end(),
00188                        trajectory2.trajectoryStates().begin(),
00189                        trajectory2.trajectoryStates().end() );
00190   }
00191 
00192   return true;
00193 }
00194 
00195 
00196 void TwoBodyDecayTrajectory::constructTsosVecWithErrors( const ReferenceTrajectory& traj1,
00197                                                          const ReferenceTrajectory& traj2,
00198                                                          const MagneticField* field )
00199 {
00200   int iTsos = 0;
00201 
00202   std::vector< TrajectoryStateOnSurface >::const_iterator itTsos;
00203 
00204   for ( itTsos = traj1.trajectoryStates().begin(); itTsos != traj1.trajectoryStates().end(); itTsos++ )
00205   {
00206     constructSingleTsosWithErrors( *itTsos, iTsos, field );
00207     iTsos++;
00208   }
00209 
00210   for ( itTsos = traj2.trajectoryStates().begin(); itTsos != traj2.trajectoryStates().end(); itTsos++ )
00211   {
00212     constructSingleTsosWithErrors( *itTsos, iTsos, field );
00213     iTsos++;
00214   }
00215 }
00216 
00217 
00218 void TwoBodyDecayTrajectory::constructSingleTsosWithErrors( const TrajectoryStateOnSurface& tsos,
00219                                                             int iTsos,
00220                                                             const MagneticField* field )
00221 {
00222   AlgebraicSymMatrix55 localErrors;
00223   const double coeff = 1e-2;
00224 
00225   double invP = tsos.localParameters().signedInverseMomentum();
00226   LocalVector p = tsos.localParameters().momentum();
00227 
00228   // rough estimate for the errors of q/p, dx/dz and dy/dz, assuming that
00229   // sigma(px) = sigma(py) = sigma(pz) = coeff*p.
00230   float dpinv = coeff*( fabs(p.x()) + fabs(p.y()) + fabs(p.z()) )*invP*invP;
00231   float dxdir = coeff*( fabs(p.x()) + fabs(p.z()) )/p.z()/p.z();
00232   float dydir = coeff*( fabs(p.y()) + fabs(p.z()) )/p.z()/p.z();
00233   localErrors[0][0] = dpinv*dpinv;
00234   localErrors[1][1] = dxdir*dxdir;
00235   localErrors[2][2] = dydir*dydir;
00236 
00237   // exact values for the errors on local x and y
00238   localErrors[3][3] = theTrajectoryPositionCov[nMeasPerHit*iTsos][nMeasPerHit*iTsos];
00239   localErrors[3][4] = theTrajectoryPositionCov[nMeasPerHit*iTsos][nMeasPerHit*iTsos+1];
00240   localErrors[4][4] = theTrajectoryPositionCov[nMeasPerHit*iTsos+1][nMeasPerHit*iTsos+1];
00241 
00242   // construct tsos with local errors
00243   theTsosVec[iTsos] =  TrajectoryStateOnSurface( tsos.localParameters(),
00244                                                  LocalTrajectoryError( localErrors ),
00245                                                  tsos.surface(),
00246                                                  field,
00247                                                  tsos.surfaceSide() );
00248 }