Reference Manual

Go to the documentation of this file.
00001 #include <memory>
00002 
00003 #include "FWCore/Framework/interface/EDProducer.h"
00004 #include "FWCore/Framework/interface/Event.h"
00005 #include "FWCore/Framework/interface/MakerMacros.h"
00006 #include "FWCore/Utilities/interface/InputTag.h"
00007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00008 
00009 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00010 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00011 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00012 
00013 #include "DataFormats/Common/interface/Handle.h"
00014 #include "DataFormats/TrackReco/interface/Track.h"
00015 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00016 #include "DataFormats/VertexReco/interface/Vertex.h"
00017 #include "DataFormats/VertexReco/interface/VertexFwd.h"
00018 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
00019 
00020 #include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h"
00021 #include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h"
00022 #include "TrackingTools/IPTools/interface/IPTools.h"
00023 #include "RecoVertex/AdaptiveVertexFinder/interface/TracksClusteringFromDisplacedSeed.h"
00024 
00025 #include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
00026 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
00027 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h"
00028 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
00029 #include "RecoVertex/MultiVertexFit/interface/MultiVertexFitter.h"
00030 
00031 //#define VTXDEBUG 1
00032 
00033 class InclusiveVertexFinder : public edm::EDProducer {
00034     public:
00035         InclusiveVertexFinder(const edm::ParameterSet &params);
00036 
00037         virtual void produce(edm::Event &event, const edm::EventSetup &es);
00038 
00039     private:
00040         bool trackFilter(const reco::TrackRef &track) const;
00041         std::pair<std::vector<reco::TransientTrack>,GlobalPoint> nearTracks(const reco::TransientTrack &seed, const std::vector<reco::TransientTrack> & tracks, const reco::Vertex & primaryVertex) const;
00042 
00043         edm::InputTag                           beamSpotCollection;
00044         edm::InputTag                           primaryVertexCollection;
00045         edm::InputTag                           trackCollection;
00046         unsigned int                            minHits;
00047         unsigned int                            maxNTracks;
00048         double                                  maxLIP;
00049         double                                  minPt;
00050         double                                  vertexMinAngleCosine;
00051         double                                  vertexMinDLen2DSig;
00052         double                                  vertexMinDLenSig;
00053 
00054         std::auto_ptr<VertexReconstructor>      vtxReco;
00055         std::auto_ptr<TracksClusteringFromDisplacedSeed>        clusterizer;
00056 
00057 };
00058 
00059 InclusiveVertexFinder::InclusiveVertexFinder(const edm::ParameterSet &params) :
00060         beamSpotCollection(params.getParameter<edm::InputTag>("beamSpot")),
00061         primaryVertexCollection(params.getParameter<edm::InputTag>("primaryVertices")),
00062         trackCollection(params.getParameter<edm::InputTag>("tracks")),
00063         minHits(params.getParameter<unsigned int>("minHits")),
00064         maxNTracks(params.getParameter<unsigned int>("maxNTracks")),
00065         maxLIP(params.getParameter<double>("maximumLongitudinalImpactParameter")),
00066         minPt(params.getParameter<double>("minPt")), //0.8
00067         vertexMinAngleCosine(params.getParameter<double>("vertexMinAngleCosine")), //0.98
00068         vertexMinDLen2DSig(params.getParameter<double>("vertexMinDLen2DSig")), //2.5
00069         vertexMinDLenSig(params.getParameter<double>("vertexMinDLenSig")), //0.5
00070         vtxReco(new ConfigurableVertexReconstructor(params.getParameter<edm::ParameterSet>("vertexReco"))),
00071         clusterizer(new TracksClusteringFromDisplacedSeed(params.getParameter<edm::ParameterSet>("clusterizer")))
00072 
00073 {
00074         produces<reco::VertexCollection>();
00075         //produces<reco::VertexCollection>("multi");
00076 }
00077 
00078 bool InclusiveVertexFinder::trackFilter(const reco::TrackRef &track) const
00079 {
00080         if (track->hitPattern().numberOfValidHits() < (int)minHits)
00081 //      if (track->hitPattern().trackerLayersWithMeasurement() < (int)minHits)
00082                 return false;
00083         if (track->pt() < minPt )
00084                 return false;
00085  
00086         return true;
00087 }
00088 
00089 void InclusiveVertexFinder::produce(edm::Event &event, const edm::EventSetup &es)
00090 {
00091         using namespace reco;
00092 
00093   double sigmacut = 3.0;
00094   double Tini = 256.;
00095   double ratio = 0.25;
00096   VertexDistance3D vdist;
00097   VertexDistanceXY vdist2d;
00098   MultiVertexFitter theMultiVertexFitter;
00099   AdaptiveVertexFitter theAdaptiveFitter(
00100                                             GeometricAnnealing(sigmacut, Tini, ratio),
00101                                             DefaultLinearizationPointFinder(),
00102                                             KalmanVertexUpdator<5>(),
00103                                             KalmanVertexTrackCompatibilityEstimator<5>(),
00104                                             KalmanVertexSmoother() );
00105 
00106 
00107         edm::Handle<BeamSpot> beamSpot;
00108         event.getByLabel(beamSpotCollection, beamSpot);
00109 
00110         edm::Handle<VertexCollection> primaryVertices;
00111         event.getByLabel(primaryVertexCollection, primaryVertices);
00112 
00113         edm::Handle<TrackCollection> tracks;
00114         event.getByLabel(trackCollection, tracks);
00115 
00116         edm::ESHandle<TransientTrackBuilder> trackBuilder;
00117         es.get<TransientTrackRecord>().get("TransientTrackBuilder",
00118                                            trackBuilder);
00119 
00120 
00121         const reco::Vertex &pv = (*primaryVertices)[0];
00122         
00123         std::vector<TransientTrack> tts;
00124         //Fill transient track vector 
00125         for(TrackCollection::const_iterator track = tracks->begin();
00126             track != tracks->end(); ++track) {
00127                 TrackRef ref(tracks, track - tracks->begin());
00128                 if (!trackFilter(ref))
00129                         continue;
00130                 if( std::abs(ref->dz(pv.position())) > maxLIP)
00131                         continue;
00132                 TransientTrack tt = trackBuilder->build(ref);
00133                 tt.setBeamSpot(*beamSpot);
00134                 tts.push_back(tt);
00135         }
00136         std::vector<TracksClusteringFromDisplacedSeed::Cluster> clusters = clusterizer->clusters(pv,tts);
00137 
00138         //Create BS object from PV to feed in the AVR
00139         BeamSpot::CovarianceMatrix cov;
00140         for(unsigned int i = 0; i < 7; i++) {
00141                 for(unsigned int j = 0; j < 7; j++) {
00142                         if (i < 3 && j < 3)
00143                                 cov(i, j) = pv.covariance(i, j);
00144                         else
00145                                 cov(i, j) = 0.0;
00146                 }
00147         }
00148         BeamSpot bs(pv.position(), 0.0, 0.0, 0.0, 0.0, cov, BeamSpot::Unknown);
00149 
00150 
00151         std::auto_ptr<VertexCollection> recoVertices(new VertexCollection);
00152         int i=0;
00153 #ifdef VTXDEBUG
00154 
00155         std::cout <<  "CLUSTERS " << clusters.size() << std::endl; 
00156 #endif
00157 
00158         for(std::vector<TracksClusteringFromDisplacedSeed::Cluster>::iterator cluster = clusters.begin();
00159             cluster != clusters.end(); ++cluster,++i)
00160         {
00161                 if(cluster->tracks.size() == 0 || cluster->tracks.size() > maxNTracks ) 
00162                      continue;
00163         
00164                 cluster->tracks.push_back(cluster->seedingTrack); //add the seed to the list of tracks to fit
00165                 std::vector<TransientVertex> vertices;
00166                 vertices = vtxReco->vertices(cluster->tracks, bs);  // attempt with config given reconstructor
00167                 TransientVertex singleFitVertex;
00168                 singleFitVertex = theAdaptiveFitter.vertex(cluster->tracks,cluster->seedPoint); //attempt with direct fitting
00169                 if(singleFitVertex.isValid())
00170                           vertices.push_back(singleFitVertex);
00171                 for(std::vector<TransientVertex>::const_iterator v = vertices.begin();
00172                     v != vertices.end(); ++v) {
00173 //                      if(v->degreesOfFreedom() > 0.2)
00174                         {
00175                          Measurement1D dlen= vdist.distance(pv,*v);
00176                          Measurement1D dlen2= vdist2d.distance(pv,*v);
00177                          reco::Vertex vv(*v);
00178 #ifdef VTXDEBUG
00179                        std::cout << "V chi2/n: " << v->normalisedChiSquared() << " ndof: " <<v->degreesOfFreedom() ;
00180                          std::cout << " dlen: " << dlen.value() << " error: " << dlen.error() << " signif: " << dlen.significance();
00181                          std::cout << " dlen2: " << dlen2.value() << " error2: " << dlen2.error() << " signif2: " << dlen2.significance();
00182                          std::cout << " pos: " << vv.position() << " error: " <<vv.xError() << " " << vv.yError() << " " << vv.zError() << std::endl;
00183 #endif
00184                          GlobalVector dir;  
00185                          std::vector<reco::TransientTrack> ts = v->originalTracks();
00186                         for(std::vector<reco::TransientTrack>::const_iterator i = ts.begin();
00187                             i != ts.end(); ++i) {
00188                                 reco::TrackRef t = i->trackBaseRef().castTo<reco::TrackRef>();
00189                                 float w = v->trackWeight(*i);
00190                                 if (w > 0.5) dir+=i->impactPointState().globalDirection();
00191 #ifdef VTXDEBUG
00192                                 std::cout << "\t[" << (*t).pt() << ": "
00193                                           << (*t).eta() << ", "
00194                                           << (*t).phi() << "], "
00195                                           << w << std::endl;
00196 #endif
00197                         }
00198                        GlobalPoint ppv(pv.position().x(),pv.position().y(),pv.position().z());
00199                        GlobalPoint sv((*v).position().x(),(*v).position().y(),(*v).position().z());
00200                        float vscal = dir.unit().dot((sv-ppv).unit()) ;
00201 //                        std::cout << "Vscal: " <<  vscal << std::endl;
00202                        if(dlen.significance() > vertexMinDLenSig  && vscal > vertexMinAngleCosine &&  v->normalisedChiSquared() < 10 && dlen2.significance() > vertexMinDLen2DSig)
00203                           {      
00204                                 recoVertices->push_back(*v);
00205 #ifdef VTXDEBUG
00206 
00207                                 std::cout << "ADDED" << std::endl;
00208 #endif
00209 
00210                          }
00211                       }
00212                    }
00213         }
00214 #ifdef VTXDEBUG
00215 
00216         std::cout <<  "Final put  " << recoVertices->size() << std::endl;
00217 #endif  
00218 
00219 
00220         event.put(recoVertices);
00221 
00222 }
00223 
00224 DEFINE_FWK_MODULE(InclusiveVertexFinder);
/data/refman/pasoursint/CMSSW_5_3_3/src/RecoVertex/AdaptiveVertexFinder/plugins/InclusiveVertexFinder.cc
