CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_1/src/RecoVertex/GaussianSumVertexFit/src/GsfVertexSmoother.cc

Go to the documentation of this file.
00001 #include "RecoVertex/GaussianSumVertexFit/interface/GsfVertexSmoother.h"
00002 #include "RecoVertex/GaussianSumVertexFit/interface/BasicMultiVertexState.h"
00003 #include "RecoVertex/GaussianSumVertexFit/interface/MultiRefittedTS.h"
00004 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00005 
00006 GsfVertexSmoother::GsfVertexSmoother(bool limit, const GsfVertexMerger * merger) :
00007   limitComponents (limit)
00008 {
00009   if (limitComponents) theMerger = merger->clone();
00010 }
00011 
00012 CachingVertex<5>
00013 GsfVertexSmoother::smooth(const CachingVertex<5> & vertex) const
00014 {
00015 
00016   std::vector<RefCountedVertexTrack> tracks = vertex.tracks();
00017   int numberOfTracks = tracks.size();
00018   if (numberOfTracks<1) return vertex;
00019 
00020   // Initial vertex for ascending fit
00021   GlobalPoint priorVertexPosition = tracks[0]->linearizedTrack()->linearizationPoint();
00022   AlgebraicSymMatrix33 we = ROOT::Math::SMatrixIdentity();
00023   GlobalError priorVertexError(we*10000);
00024 
00025   std::vector<RefCountedVertexTrack> initialTracks;
00026   CachingVertex<5> fitVertex(priorVertexPosition,priorVertexError,initialTracks,0);
00027   //In case prior vertex was used.
00028   if (vertex.hasPrior()) {
00029     VertexState priorVertexState = vertex.priorVertexState();
00030     fitVertex = CachingVertex<5>(priorVertexState, priorVertexState,
00031                 initialTracks,0);
00032   }
00033 
00034   // vertices from ascending fit
00035   std::vector<CachingVertex<5> > ascendingFittedVertices;
00036   ascendingFittedVertices.reserve(numberOfTracks);
00037   ascendingFittedVertices.push_back(fitVertex);
00038 
00039   // ascending fit
00040   for (std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin();
00041           i != (tracks.end()-1); ++i) {
00042     fitVertex = theUpdator.add(fitVertex,*i);
00043     if (limitComponents) fitVertex = theMerger->merge(fitVertex);
00044     ascendingFittedVertices.push_back(fitVertex);
00045   }
00046 
00047   // Initial vertex for descending fit
00048   priorVertexPosition = tracks[0]->linearizedTrack()->linearizationPoint();
00049   priorVertexError = GlobalError(we*10000);
00050   fitVertex = CachingVertex<5>(priorVertexPosition,priorVertexError,initialTracks,0);
00051 
00052   // vertices from descending fit
00053   std::vector<CachingVertex<5> > descendingFittedVertices;
00054   descendingFittedVertices.reserve(numberOfTracks);
00055   descendingFittedVertices.push_back(fitVertex);
00056 
00057   // descending fit
00058   for (std::vector<RefCountedVertexTrack>::const_iterator i = (tracks.end()-1);
00059           i != (tracks.begin()); --i) {
00060     fitVertex = theUpdator.add(fitVertex,*i);
00061     if (limitComponents) fitVertex = theMerger->merge(fitVertex);
00062     descendingFittedVertices.insert(descendingFittedVertices.begin(), fitVertex);
00063   }
00064 
00065   std::vector<RefCountedVertexTrack> newTracks;
00066   double smoothedChi2 = 0.;  // Smoothed chi2
00067 
00068   // Track refit
00069   for(std::vector<RefCountedVertexTrack>::const_iterator i = tracks.begin();
00070         i != tracks.end();i++)
00071   {
00072     int indexNumber = i-tracks.begin();
00073     //First, combine the vertices:
00074     VertexState meanedVertex = 
00075          meanVertex(ascendingFittedVertices[indexNumber].vertexState(), 
00076                     descendingFittedVertices[indexNumber].vertexState());
00077     if (limitComponents) meanedVertex = theMerger->merge(meanedVertex);
00078     // Add the vertex and smooth the track:
00079     TrackChi2Pair thePair = vertexAndTrackUpdate(meanedVertex, *i, vertex.position());
00080     smoothedChi2 += thePair.second.second;
00081     newTracks.push_back( theVTFactory.vertexTrack((**i).linearizedTrack(),
00082         vertex.vertexState(), thePair.first, thePair.second.second,
00083         AlgebraicSymMatrixOO(), (**i).weight()) );
00084   }
00085 
00086   if  (vertex.hasPrior()) {
00087     smoothedChi2 += priorVertexChi2(vertex.priorVertexState(), vertex.vertexState());
00088     return CachingVertex<5>(vertex.priorVertexState(), vertex.vertexState(),
00089                 newTracks, smoothedChi2);
00090   } else {
00091     return CachingVertex<5>(vertex.vertexState(), newTracks, smoothedChi2);
00092   }
00093 }
00094 
00095 GsfVertexSmoother::TrackChi2Pair 
00096 GsfVertexSmoother::vertexAndTrackUpdate(const VertexState & oldVertex,
00097         const RefCountedVertexTrack track, const GlobalPoint & referencePosition) const
00098 {
00099 
00100   VSC prevVtxComponents = oldVertex.components();
00101 
00102   if (prevVtxComponents.empty()) {
00103   throw VertexException
00104     ("GsfVertexSmoother::(Previous) Vertex to update has no components");
00105   }
00106 
00107   LTC ltComponents = track->linearizedTrack()->components();
00108   if (ltComponents.empty()) {
00109   throw VertexException
00110     ("GsfVertexSmoother::Track to add to vertex has no components");
00111   }
00112   float trackWeight = track->weight();
00113 
00114   std::vector<RefittedTrackComponent> newTrackComponents;
00115   newTrackComponents.reserve(prevVtxComponents.size()*ltComponents.size());
00116 
00117   for (VSC::iterator vertexCompIter = prevVtxComponents.begin();
00118         vertexCompIter != prevVtxComponents.end(); vertexCompIter++ ) {
00119 
00120     for (LTC::iterator trackCompIter = ltComponents.begin();
00121         trackCompIter != ltComponents.end(); trackCompIter++ ) {
00122       newTrackComponents.push_back
00123         (createNewComponent(*vertexCompIter, *trackCompIter, trackWeight));
00124     }
00125   }
00126 
00127   return assembleTrackComponents(newTrackComponents, referencePosition);
00128 }
00129 
00135 GsfVertexSmoother::TrackChi2Pair GsfVertexSmoother::assembleTrackComponents(
00136         const std::vector<GsfVertexSmoother::RefittedTrackComponent> & trackComponents,
00137         const GlobalPoint & referencePosition)
00138         const
00139 {
00140 
00141   //renormalize weights
00142 
00143   double totalWeight = 0.;
00144   double totalVtxChi2 = 0., totalTrkChi2 = 0.;
00145 
00146   for (std::vector<RefittedTrackComponent>::const_iterator iter = trackComponents.begin();
00147     iter != trackComponents.end(); ++iter ) {
00148     totalWeight += iter->first.second;
00149     totalVtxChi2 += iter->second.first  * iter->first.second ;
00150     totalTrkChi2 += iter->second.second * iter->first.second ;
00151   }
00152 
00153   totalVtxChi2 /= totalWeight ;
00154   totalTrkChi2 /= totalWeight ;
00155 
00156   std::vector<RefCountedRefittedTrackState> reWeightedRTSC;
00157   reWeightedRTSC.reserve(trackComponents.size());
00158   
00159 
00160   for (std::vector<RefittedTrackComponent>::const_iterator iter = trackComponents.begin();
00161     iter != trackComponents.end(); ++iter ) {
00162     if (iter->second.first!=0) {
00163       reWeightedRTSC.push_back(iter->first.first->stateWithNewWeight(iter->second.first/totalWeight));
00164     }
00165   }
00166 
00167   RefCountedRefittedTrackState finalRTS = 
00168     RefCountedRefittedTrackState(new MultiRefittedTS(reWeightedRTSC, referencePosition));
00169   return TrackChi2Pair(finalRTS, VtxTrkChi2Pair(totalVtxChi2, totalTrkChi2));
00170 }
00171 
00172 
00178 GsfVertexSmoother::RefittedTrackComponent 
00179 GsfVertexSmoother::createNewComponent(const VertexState & oldVertex,
00180          const RefCountedLinearizedTrackState linTrack, float trackWeight) const
00181 {
00182 
00183   int sign =+1;
00184 
00185   // Weight of the component in the mixture (non-normalized)
00186   double weightInMixture = theWeightCalculator.calculate(oldVertex, linTrack, 1000000000.);
00187 
00188   // position estimate of the component
00189   VertexState newVertex = kalmanVertexUpdator.positionUpdate(oldVertex,
00190         linTrack, trackWeight, sign);
00191 
00192   KalmanVertexTrackUpdator<5>::trackMatrixPair thePair = 
00193         theVertexTrackUpdator.trackRefit(newVertex, linTrack, trackWeight);
00194 
00195   //Chi**2 contribution of the track component
00196   double vtxChi2 = helper.vertexChi2(oldVertex, newVertex);
00197   std::pair<bool, double> trkCi2 = helper.trackParameterChi2(linTrack, thePair.first);
00198 
00199   return RefittedTrackComponent(TrackWeightPair(thePair.first, weightInMixture), 
00200                         VtxTrkChi2Pair(vtxChi2, trkCi2.second));
00201 }
00202 
00203 
00204 VertexState
00205 GsfVertexSmoother::meanVertex(const VertexState & vertexA,
00206                               const VertexState & vertexB) const
00207 {
00208   std::vector<VertexState> vsCompA = vertexA.components();
00209   std::vector<VertexState> vsCompB = vertexB.components();
00210   std::vector<VertexState> finalVS;
00211   finalVS.reserve(vsCompA.size()*vsCompB.size());
00212   for (std::vector<VertexState>::iterator iA = vsCompA.begin(); iA!= vsCompA.end(); ++iA)
00213   {
00214     for (std::vector<VertexState>::iterator iB = vsCompB.begin(); iB!= vsCompB.end(); ++iB)
00215     {
00216       AlgebraicSymMatrix33 newWeight = iA->weight().matrix_new() +
00217                                      iB->weight().matrix_new();
00218       AlgebraicVector3 newWtP = iA->weightTimesPosition() +
00219                                iB->weightTimesPosition();
00220       double newWeightInMixture = iA->weightInMixture() *
00221                                   iB->weightInMixture();
00222       finalVS.push_back( VertexState(newWtP, newWeight, newWeightInMixture) );
00223     }
00224   }
00225   #ifndef CMS_NO_COMPLEX_RETURNS
00226     return VertexState(new BasicMultiVertexState(finalVS));
00227   #else
00228     VertexState theFinalVM(new BasicMultiVertexState(finalVS));
00229     return theFinalVM;
00230   #endif
00231 }
00232 
00233 
00234 double GsfVertexSmoother::priorVertexChi2(
00235         const VertexState priorVertex, const VertexState fittedVertex) const
00236 {
00237   std::vector<VertexState> priorVertexComp  = priorVertex.components();
00238   std::vector<VertexState> fittedVertexComp = fittedVertex.components();
00239   double vetexChi2 = 0.;
00240   for (std::vector<VertexState>::iterator pvI = priorVertexComp.begin(); 
00241         pvI!= priorVertexComp.end(); ++pvI)
00242   {
00243     for (std::vector<VertexState>::iterator fvI = fittedVertexComp.begin(); 
00244         fvI!= fittedVertexComp.end(); ++fvI)
00245     {
00246       vetexChi2 += (pvI->weightInMixture())*(fvI->weightInMixture())*
00247                         helper.vertexChi2(*pvI, *fvI);
00248     }
00249   }
00250   return vetexChi2;
00251 }