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         std::auto_ptr<VertexCollection> recoVertices(new VertexCollection);
00122         if(primaryVertices->size()!=0) {
00123      
00124         const reco::Vertex &pv = (*primaryVertices)[0];
00125         
00126         std::vector<TransientTrack> tts;
00127         //Fill transient track vector 
00128         for(TrackCollection::const_iterator track = tracks->begin();
00129             track != tracks->end(); ++track) {
00130                 TrackRef ref(tracks, track - tracks->begin());
00131                 if (!trackFilter(ref))
00132                         continue;
00133                 if( std::abs(ref->dz(pv.position())) > maxLIP)
00134                         continue;
00135                 TransientTrack tt = trackBuilder->build(ref);
00136                 tt.setBeamSpot(*beamSpot);
00137                 tts.push_back(tt);
00138         }
00139         std::vector<TracksClusteringFromDisplacedSeed::Cluster> clusters = clusterizer->clusters(pv,tts);
00140 
00141         //Create BS object from PV to feed in the AVR
00142         BeamSpot::CovarianceMatrix cov;
00143         for(unsigned int i = 0; i < 7; i++) {
00144                 for(unsigned int j = 0; j < 7; j++) {
00145                         if (i < 3 && j < 3)
00146                                 cov(i, j) = pv.covariance(i, j);
00147                         else
00148                                 cov(i, j) = 0.0;
00149                 }
00150         }
00151         BeamSpot bs(pv.position(), 0.0, 0.0, 0.0, 0.0, cov, BeamSpot::Unknown);
00152 
00153 
00154         int i=0;
00155 #ifdef VTXDEBUG
00156 
00157         std::cout <<  "CLUSTERS " << clusters.size() << std::endl; 
00158 #endif
00159 
00160         for(std::vector<TracksClusteringFromDisplacedSeed::Cluster>::iterator cluster = clusters.begin();
00161             cluster != clusters.end(); ++cluster,++i)
00162         {
00163                 if(cluster->tracks.size() == 0 || cluster->tracks.size() > maxNTracks ) 
00164                      continue;
00165         
00166                 cluster->tracks.push_back(cluster->seedingTrack); //add the seed to the list of tracks to fit
00167                 std::vector<TransientVertex> vertices;
00168                 vertices = vtxReco->vertices(cluster->tracks, bs);  // attempt with config given reconstructor
00169                 TransientVertex singleFitVertex;
00170                 singleFitVertex = theAdaptiveFitter.vertex(cluster->tracks,cluster->seedPoint); //attempt with direct fitting
00171                 if(singleFitVertex.isValid())
00172                           vertices.push_back(singleFitVertex);
00173                 for(std::vector<TransientVertex>::const_iterator v = vertices.begin();
00174                     v != vertices.end(); ++v) {
00175 //                      if(v->degreesOfFreedom() > 0.2)
00176                         {
00177                          Measurement1D dlen= vdist.distance(pv,*v);
00178                          Measurement1D dlen2= vdist2d.distance(pv,*v);
00179                          reco::Vertex vv(*v);
00180 #ifdef VTXDEBUG
00181                        std::cout << "V chi2/n: " << v->normalisedChiSquared() << " ndof: " <<v->degreesOfFreedom() ;
00182                          std::cout << " dlen: " << dlen.value() << " error: " << dlen.error() << " signif: " << dlen.significance();
00183                          std::cout << " dlen2: " << dlen2.value() << " error2: " << dlen2.error() << " signif2: " << dlen2.significance();
00184                          std::cout << " pos: " << vv.position() << " error: " <<vv.xError() << " " << vv.yError() << " " << vv.zError() << std::endl;
00185 #endif
00186                          GlobalVector dir;  
00187                          std::vector<reco::TransientTrack> ts = v->originalTracks();
00188                         for(std::vector<reco::TransientTrack>::const_iterator i = ts.begin();
00189                             i != ts.end(); ++i) {
00190                                 reco::TrackRef t = i->trackBaseRef().castTo<reco::TrackRef>();
00191                                 float w = v->trackWeight(*i);
00192                                 if (w > 0.5) dir+=i->impactPointState().globalDirection();
00193 #ifdef VTXDEBUG
00194                                 std::cout << "\t[" << (*t).pt() << ": "
00195                                           << (*t).eta() << ", "
00196                                           << (*t).phi() << "], "
00197                                           << w << std::endl;
00198 #endif
00199                         }
00200                        GlobalPoint ppv(pv.position().x(),pv.position().y(),pv.position().z());
00201                        GlobalPoint sv((*v).position().x(),(*v).position().y(),(*v).position().z());
00202                        float vscal = dir.unit().dot((sv-ppv).unit()) ;
00203 //                        std::cout << "Vscal: " <<  vscal << std::endl;
00204                        if(dlen.significance() > vertexMinDLenSig  && vscal > vertexMinAngleCosine &&  v->normalisedChiSquared() < 10 && dlen2.significance() > vertexMinDLen2DSig)
00205                           {      
00206                                 recoVertices->push_back(*v);
00207 #ifdef VTXDEBUG
00208 
00209                                 std::cout << "ADDED" << std::endl;
00210 #endif
00211 
00212                          }
00213                       }
00214                    }
00215         }
00216 #ifdef VTXDEBUG
00217 
00218         std::cout <<  "Final put  " << recoVertices->size() << std::endl;
00219 #endif  
00220         }
00221  
00222         event.put(recoVertices);
00223 
00224 }
00225 
00226 DEFINE_FWK_MODULE(InclusiveVertexFinder);
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/RecoVertex/AdaptiveVertexFinder/plugins/InclusiveVertexFinder.cc
