CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoVertex/AdaptiveVertexFinder/src/AdaptiveVertexReconstructor.cc

Go to the documentation of this file.
00001 #include "RecoVertex/AdaptiveVertexFinder/interface/AdaptiveVertexReconstructor.h"
00002 #include "RecoVertex/VertexTools/interface/GeometricAnnealing.h"
00003 #include "FWCore/Utilities/interface/EDMException.h"
00004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00005 #include "RecoVertex/VertexTools/interface/DummyVertexSmoother.h"
00006 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
00007 #include <algorithm>
00008 
00009 using namespace std;
00010 
00011 TransientVertex AdaptiveVertexReconstructor::cleanUp ( const TransientVertex & old ) const
00012 {
00013   vector < reco::TransientTrack > trks = old.originalTracks();
00014   vector < reco::TransientTrack > newtrks;
00015   TransientVertex::TransientTrackToFloatMap mp;
00016   static const float minweight = 1.e-8; // discard all tracks with lower weight
00017   for ( vector< reco::TransientTrack >::const_iterator i=trks.begin();
00018         i!=trks.end() ; ++i )
00019   {
00020     if ( old.trackWeight ( *i ) > minweight )
00021     {
00022       newtrks.push_back ( *i );
00023       mp[*i]=old.trackWeight ( *i );
00024     }
00025   }
00026 
00027   TransientVertex ret;
00028 
00029   if ( old.hasPrior() )
00030   {
00031     VertexState priorstate ( old.priorPosition(), old.priorError() );
00032     ret=TransientVertex ( priorstate, old.vertexState(), newtrks,
00033         old.totalChiSquared(), old.degreesOfFreedom() );
00034   } else {
00035     ret=TransientVertex ( old.vertexState(), newtrks,
00036                           old.totalChiSquared(), old.degreesOfFreedom() );
00037   }
00038   ret.weightMap ( mp ); // set weight map
00039 
00040   if ( old.hasRefittedTracks() )
00041   {
00042     // we have refitted tracks -- copy them!
00043     vector < reco::TransientTrack > newrfs;
00044     vector < reco::TransientTrack > oldrfs=old.refittedTracks();
00045     vector < reco::TransientTrack >::const_iterator origtrkiter=trks.begin();
00046     for ( vector< reco::TransientTrack >::const_iterator i=oldrfs.begin(); i!=oldrfs.end() ; ++i )
00047     {
00048       if ( old.trackWeight ( *origtrkiter ) > minweight )
00049       {
00050         newrfs.push_back ( *i );
00051       }
00052       origtrkiter++;
00053     }
00054     if ( newrfs.size() ) ret.refittedTracks ( newrfs ); // copy refitted tracks
00055   }
00056 
00057   if ( ret.refittedTracks().size() > ret.originalTracks().size() )
00058   {
00059     edm::LogError("AdaptiveVertexReconstructor" )
00060       << "More refitted tracks than original tracks!";
00061   }
00062 
00063   return ret;
00064 }
00065 
00066 void AdaptiveVertexReconstructor::erase (
00067     const TransientVertex & newvtx,
00068     set < reco::TransientTrack > & remainingtrks,
00069     float w ) const
00070 {
00071   /*
00072    * Erase tracks that are in newvtx from remainingtrks 
00073    * But erase only if trackweight > w
00074    */
00075   const vector < reco::TransientTrack > & origtrks = newvtx.originalTracks();
00076   for ( vector< reco::TransientTrack >::const_iterator i=origtrks.begin();
00077         i!=origtrks.end(); ++i )
00078   {
00079     double weight = newvtx.trackWeight ( *i );
00080     if ( weight > w )
00081     {
00082       remainingtrks.erase ( *i );
00083     };
00084   };
00085 }
00086 
00087 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor(
00088     float primcut, float seccut, float min_weight, bool smoothing ) :
00089   thePrimaryFitter ( 0 ), theSecondaryFitter ( 0 ),
00090        theMinWeight( min_weight ), theWeightThreshold ( 0.001 )
00091 {
00092   setupFitters ( primcut, 256., 0.25, seccut, 256., 0.25, smoothing );
00093 }
00094 
00095 AdaptiveVertexReconstructor::~AdaptiveVertexReconstructor()
00096 {
00097   if ( thePrimaryFitter ) delete thePrimaryFitter;
00098   if ( theSecondaryFitter ) delete theSecondaryFitter;
00099 }
00100 
00101 void AdaptiveVertexReconstructor::setupFitters ( float primcut, 
00102     float primT, float primr, float seccut, float secT,
00103     float secr, bool smoothing )
00104 {
00105   VertexSmoother<5> * smoother ;
00106   if ( smoothing )
00107   {
00108     smoother = new KalmanVertexSmoother();
00109   } else {
00110     smoother = new DummyVertexSmoother<5>();
00111   }
00112 
00113   if ( thePrimaryFitter ) delete thePrimaryFitter;
00114   if ( theSecondaryFitter ) delete theSecondaryFitter;
00115 
00116   /*
00117   edm::LogError ("AdaptiveVertexReconstructor" )
00118     << "Tini and r are hardcoded now!";
00119     */
00120   thePrimaryFitter = new AdaptiveVertexFitter ( GeometricAnnealing ( primcut, primT, primr ), DefaultLinearizationPointFinder(),
00121       KalmanVertexUpdator<5>(), KalmanVertexTrackCompatibilityEstimator<5>(), *smoother );
00122   thePrimaryFitter->setWeightThreshold ( theWeightThreshold );
00123   // if the primary fails, sth is wrong, so here we set a threshold on the weight.
00124   theSecondaryFitter = new AdaptiveVertexFitter ( GeometricAnnealing ( seccut, secT, secr ), DefaultLinearizationPointFinder(),
00125       KalmanVertexUpdator<5>(), 
00126       KalmanVertexTrackCompatibilityEstimator<5>(), *smoother );
00127   theSecondaryFitter->setWeightThreshold ( 0. );
00128   // need to set it or else we have 
00129   // unwanted exceptions to deal with.
00130   // cleanup can come later!
00131   delete smoother;
00132 }
00133 
00134 AdaptiveVertexReconstructor::AdaptiveVertexReconstructor( const edm::ParameterSet & m )
00135   : thePrimaryFitter(0), theSecondaryFitter(0), theMinWeight(0.5), theWeightThreshold ( 0.001 )
00136 {
00137   float primcut = 2.0;
00138   float seccut = 6.0;
00139   bool smoothing=false;
00140   // float primT = 4096.;
00141   // float primr = 0.125;
00142   float primT = 256.;
00143   float primr = 0.25;
00144   float secT = 256.;
00145   float secr = 0.25;
00146   
00147   try {
00148     primcut =  m.getParameter<double>("primcut");
00149     primT =  m.getParameter<double>("primT");
00150     primr =  m.getParameter<double>("primr");
00151     seccut =  m.getParameter<double>("seccut");
00152     secT =  m.getParameter<double>("secT");
00153     secr =  m.getParameter<double>("secr");
00154     theMinWeight = m.getParameter<double>("minweight");
00155     theWeightThreshold = m.getParameter<double>("weightthreshold");
00156     smoothing =  m.getParameter<bool>("smoothing");
00157   } catch ( edm::Exception & e ) {
00158     edm::LogError ("AdaptiveVertexReconstructor") << e.what();
00159   }
00160 
00161   setupFitters ( primcut, primT, primr, seccut, secT, secr, smoothing );
00162 }
00163 
00164 vector<TransientVertex> 
00165     AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack> & t, 
00166         const reco::BeamSpot & s ) const
00167 {
00168   return vertices ( vector<reco::TransientTrack>(), t, s, false, true );
00169 }
00170 
00171 vector<TransientVertex> 
00172     AdaptiveVertexReconstructor::vertices(const vector<reco::TransientTrack> & primaries,
00173        const vector<reco::TransientTrack > & tracks, const reco::BeamSpot & s ) const
00174 {
00175   return vertices ( primaries, tracks, s, true, true );
00176 }
00177 
00178 
00179 vector<TransientVertex> AdaptiveVertexReconstructor::vertices (
00180     const vector<reco::TransientTrack> & tracks ) const
00181 {
00182   return vertices ( vector<reco::TransientTrack>(), tracks, reco::BeamSpot(), 
00183                     false, false );
00184 }
00185 
00186 vector<TransientVertex> AdaptiveVertexReconstructor::vertices (
00187     const vector<reco::TransientTrack> & primaries,
00188     const vector<reco::TransientTrack> & tracks,
00189     const reco::BeamSpot & s, bool has_primaries, bool usespot ) const
00190 {
00191   vector < TransientVertex > ret;
00192   set < reco::TransientTrack > remainingtrks;
00193 
00194   copy(tracks.begin(), tracks.end(), 
00195             inserter(remainingtrks, remainingtrks.begin()));
00196 
00197   int ctr=0;
00198   unsigned int n_tracks = remainingtrks.size();
00199 
00200   // cout << "[AdaptiveVertexReconstructor] DEBUG ::vertices!!" << endl;
00201   try {
00202     while ( remainingtrks.size() > 1 )
00203     {
00204       /*
00205       cout << "[AdaptiveVertexReconstructor] next round: "
00206            << remainingtrks.size() << endl;
00207            */
00208       ctr++;
00209       const AdaptiveVertexFitter * fitter = theSecondaryFitter;
00210       if ( ret.size() == 0 )
00211       {
00212         fitter = thePrimaryFitter;
00213       };
00214       vector < reco::TransientTrack > fittrks;
00215       fittrks.reserve ( remainingtrks.size() );
00216 
00217       copy(remainingtrks.begin(), remainingtrks.end(), back_inserter(fittrks));
00218 
00219       TransientVertex tmpvtx;
00220       if ( (ret.size() == 0) && has_primaries )
00221       {
00222         // add the primaries to the fitted tracks.
00223         copy ( primaries.begin(), primaries.end(), back_inserter(fittrks) );
00224       }
00225       if ( (ret.size() == 0) && usespot )
00226       {
00227         tmpvtx=fitter->vertex ( fittrks, s );
00228       } else {
00229         tmpvtx=fitter->vertex ( fittrks );
00230       }
00231       TransientVertex newvtx = cleanUp ( tmpvtx );
00232       ret.push_back ( newvtx );
00233       erase ( newvtx, remainingtrks, theMinWeight );
00234       if ( n_tracks == remainingtrks.size() )
00235       {
00236         if ( usespot )
00237         {
00238           // try once more without beamspot constraint!
00239           usespot=false;
00240           LogDebug("AdaptiveVertexReconstructor") 
00241             << "no tracks in vertex. trying again without beamspot constraint!";
00242           continue;
00243         }
00244         LogDebug("AdaptiveVertexReconstructor") << "all tracks (" << n_tracks
00245              << ") would be recycled for next fit. Trying with low threshold!";
00246         erase ( newvtx, remainingtrks, 1.e-5 );
00247         if ( n_tracks == remainingtrks.size() )
00248         {
00249           LogDebug("AdaptiveVertexReconstructor") << "low threshold didnt help! "
00250                                                   << "Discontinue procedure!";
00251           break;
00252         }
00253       };
00254 
00255       // cout << "[AdaptiveVertexReconstructor] erased" << endl;
00256       n_tracks = remainingtrks.size();
00257     };
00258   } catch ( VertexException & v ) {
00259     // Will catch all (not enough significant tracks exceptions.
00260     // in this case, the iteration can safely terminate.
00261   };
00262 
00263   return cleanUpVertices ( ret );
00264 }
00265 
00266 vector<TransientVertex> AdaptiveVertexReconstructor::cleanUpVertices (
00267     const vector < TransientVertex > & old ) const
00268 {
00269   vector < TransientVertex > ret;
00270   for ( vector< TransientVertex >::const_iterator i=old.begin(); i!=old.end() ; ++i )
00271   {
00272     if (!(i->hasTrackWeight()))
00273     { // if we dont have track weights, we take the vtx
00274       ret.push_back ( *i );
00275       continue;
00276     }
00277 
00278     // maybe this should be replaced with asking for the ndf ...
00279     // e.g. if ( ndf > - 1. )
00280     int nsig=0; // number of significant tracks. 
00281     TransientVertex::TransientTrackToFloatMap wm = i->weightMap();
00282     for ( TransientVertex::TransientTrackToFloatMap::const_iterator w=wm.begin(); w!=wm.end() ; ++w )
00283     {
00284       if (w->second > theWeightThreshold) nsig++;
00285     }
00286     if ( nsig > 1 ) ret.push_back ( *i );
00287   }
00288 
00289   return ret;
00290 }