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/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h"
00024 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexUpdator.h"
00025 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexTrackCompatibilityEstimator.h"
00026 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexSmoother.h"
00027 #include "RecoVertex/MultiVertexFit/interface/MultiVertexFitter.h"
00028 
00029 class InclusiveVertexFinder : public edm::EDProducer {
00030     public:
00031         InclusiveVertexFinder(const edm::ParameterSet &params);
00032 
00033         virtual void produce(edm::Event &event, const edm::EventSetup &es);
00034 
00035     private:
00036         bool trackFilter(const reco::TrackRef &track) const;
00037         std::pair<std::vector<reco::TransientTrack>,GlobalPoint> nearTracks(const reco::TransientTrack &seed, const std::vector<reco::TransientTrack> & tracks, const reco::Vertex & primaryVertex) const;
00038 
00039         edm::InputTag                           beamSpotCollection;
00040         edm::InputTag                           primaryVertexCollection;
00041         edm::InputTag                           trackCollection;
00042         unsigned int                            minHits;
00043         double                                  minPt;
00044         double                                  min3DIPSignificance;
00045         double                                  min3DIPValue;
00046         double                                  clusterMaxDistance;
00047         double                                  clusterMaxSignificance;
00048         double                                  clusterScale;
00049         double                                  clusterMinAngleCosine;
00050         double                                  vertexMinAngleCosine;
00051         double                                  vertexMinDLen2DSig;
00052         double                                  vertexMinDLenSig;
00053 
00054         std::auto_ptr<VertexReconstructor>      vtxReco;
00055 
00056 };
00057 
00058 InclusiveVertexFinder::InclusiveVertexFinder(const edm::ParameterSet &params) :
00059         beamSpotCollection(params.getParameter<edm::InputTag>("beamSpot")),
00060         primaryVertexCollection(params.getParameter<edm::InputTag>("primaryVertices")),
00061         trackCollection(params.getParameter<edm::InputTag>("tracks")),
00062         minHits(params.getParameter<unsigned int>("minHits")),
00063         minPt(params.getParameter<double>("minPt")), //0.8
00064         min3DIPSignificance(params.getParameter<double>("seedMin3DIPSignificance")),
00065         min3DIPValue(params.getParameter<double>("seedMin3DIPValue")),
00066         clusterMaxDistance(params.getParameter<double>("clusterMaxDistance")),
00067         clusterMaxSignificance(params.getParameter<double>("clusterMaxSignificance")), //3
00068         clusterScale(params.getParameter<double>("clusterScale")),//10.
00069         clusterMinAngleCosine(params.getParameter<double>("clusterMinAngleCosine")), //0.0
00070         vertexMinAngleCosine(params.getParameter<double>("vertexMinAngleCosine")), //0.98
00071         vertexMinDLen2DSig(params.getParameter<double>("vertexMinDLen2DSig")), //2.5
00072         vertexMinDLenSig(params.getParameter<double>("vertexMinDLenSig")), //0.5
00073         vtxReco(new ConfigurableVertexReconstructor(params.getParameter<edm::ParameterSet>("vertexReco")))
00074 
00075 {
00076         produces<reco::VertexCollection>();
00077         //produces<reco::VertexCollection>("multi");
00078 }
00079 
00080 bool InclusiveVertexFinder::trackFilter(const reco::TrackRef &track) const
00081 {
00082         if (track->hitPattern().trackerLayersWithMeasurement() < (int)minHits)
00083                 return false;
00084         if (track->pt() < minPt )
00085                 return false;
00086         
00087         return true;
00088 }
00089 
00090 std::pair<std::vector<reco::TransientTrack>,GlobalPoint> InclusiveVertexFinder::nearTracks(const reco::TransientTrack &seed, const std::vector<reco::TransientTrack> & tracks, const  reco::Vertex & primaryVertex) const
00091 {
00092       VertexDistance3D distanceComputer;
00093       GlobalPoint pv(primaryVertex.position().x(),primaryVertex.position().y(),primaryVertex.position().z());
00094       std::vector<reco::TransientTrack> result;
00095       TwoTrackMinimumDistance dist;
00096       GlobalPoint seedingPoint;
00097       float sumWeights=0;
00098       std::pair<bool,Measurement1D> ipSeed = IPTools::absoluteImpactParameter3D(seed,primaryVertex);
00099       float pvDistance = ipSeed.second.value();
00100 
00101       for(std::vector<reco::TransientTrack>::const_iterator tt = tracks.begin();tt!=tracks.end(); ++tt )   {
00102 
00103        if(*tt==seed) continue;
00104 
00105        std::pair<bool,Measurement1D> ip = IPTools::absoluteImpactParameter3D(*tt,primaryVertex);
00106        if(dist.calculate(tt->impactPointState(),seed.impactPointState()))
00107             {
00108                  GlobalPoint ttPoint          = dist.points().first;
00109                  GlobalError ttPointErr       = tt->impactPointState().cartesianError().position();
00110                  GlobalPoint seedPosition     = dist.points().second;
00111                  GlobalError seedPositionErr  = seed.impactPointState().cartesianError().position();
00112                  Measurement1D m = distanceComputer.distance(VertexState(seedPosition,seedPositionErr), VertexState(ttPoint, ttPointErr));
00113                  GlobalPoint cp(dist.crossingPoint()); 
00114 
00115 
00116                  float distanceFromPV =  (dist.points().second-pv).mag();
00117                  float distance = dist.distance();
00118                  GlobalVector trackDir2D(tt->impactPointState().globalDirection().x(),tt->impactPointState().globalDirection().y(),0.); 
00119                  GlobalVector seedDir2D(seed.impactPointState().globalDirection().x(),seed.impactPointState().globalDirection().y(),0.); 
00120                  float dotprodTrackSeed2D = trackDir2D.unit().dot(seedDir2D.unit());
00121 
00122                  float dotprodTrack = (dist.points().first-pv).unit().dot(tt->impactPointState().globalDirection().unit());
00123                  float dotprodSeed = (dist.points().second-pv).unit().dot(seed.impactPointState().globalDirection().unit());
00124 
00125                  float w = distanceFromPV*distanceFromPV/(pvDistance*distance);
00126 
00127                  if(m.significance() < clusterMaxSignificance && 
00128                     dotprodSeed > clusterMinAngleCosine && //Angles between PV-PCAonSeed vectors and seed directions
00129                     dotprodTrack > clusterMinAngleCosine && //Angles between PV-PCAonTrack vectors and track directions
00130                     dotprodTrackSeed2D > clusterMinAngleCosine && //Angle between track and seed
00131                     distance*clusterScale < distanceFromPV*distanceFromPV/pvDistance && // cut scaling with track density
00132                     distance < clusterMaxDistance)  // absolute distance cut
00133                  {
00134                      result.push_back(*tt);
00135                      seedingPoint = GlobalPoint(cp.x()*w+seedingPoint.x(),cp.y()*w+seedingPoint.y(),cp.z()*w+seedingPoint.z());  
00136                      sumWeights+=w; 
00137                  }
00138             }
00139        }
00140 
00141 seedingPoint = GlobalPoint(seedingPoint.x()/sumWeights,seedingPoint.y()/sumWeights,seedingPoint.z()/sumWeights);
00142 return std::pair<std::vector<reco::TransientTrack>,GlobalPoint>(result,seedingPoint);
00143 
00144 }
00145 
00146 
00147 void InclusiveVertexFinder::produce(edm::Event &event, const edm::EventSetup &es)
00148 {
00149         using namespace reco;
00150 
00151   double sigmacut = 3.0;
00152   double Tini = 256.;
00153   double ratio = 0.25;
00154   VertexDistance3D vdist;
00155   VertexDistanceXY vdist2d;
00156   MultiVertexFitter theMultiVertexFitter;
00157   AdaptiveVertexFitter theAdaptiveFitter(
00158                                             GeometricAnnealing(sigmacut, Tini, ratio),
00159                                             DefaultLinearizationPointFinder(),
00160                                             KalmanVertexUpdator<5>(),
00161                                             KalmanVertexTrackCompatibilityEstimator<5>(),
00162                                             KalmanVertexSmoother() );
00163 
00164 
00165         edm::Handle<BeamSpot> beamSpot;
00166         event.getByLabel(beamSpotCollection, beamSpot);
00167 
00168         edm::Handle<VertexCollection> primaryVertices;
00169         event.getByLabel(primaryVertexCollection, primaryVertices);
00170 
00171         edm::Handle<TrackCollection> tracks;
00172         event.getByLabel(trackCollection, tracks);
00173 
00174         edm::ESHandle<TransientTrackBuilder> trackBuilder;
00175         es.get<TransientTrackRecord>().get("TransientTrackBuilder",
00176                                            trackBuilder);
00177 
00178 
00179         const reco::Vertex &pv = (*primaryVertices)[0];
00180         
00181         std::vector<TransientTrack> tts;
00182         std::vector<TransientTrack> seeds;
00183         //Fill transient track vector and find seeds
00184         for(TrackCollection::const_iterator track = tracks->begin();
00185             track != tracks->end(); ++track) {
00186                 TrackRef ref(tracks, track - tracks->begin());
00187                 if (!trackFilter(ref))
00188                         continue;
00189 
00190                 TransientTrack tt = trackBuilder->build(ref);
00191                 tt.setBeamSpot(*beamSpot);
00192                 tts.push_back(tt);
00193                 std::pair<bool,Measurement1D> ip = IPTools::absoluteImpactParameter3D(tt,pv);
00194 //                std::cout << "track: " << ip.second.value() << " " << ip.second.significance() << " " << track->hitPattern().trackerLayersWithMeasurement() << " " << track->pt() << " " << track->eta() << std::endl;
00195                 if(ip.first && ip.second.value() >= min3DIPValue && ip.second.significance() >= min3DIPSignificance)
00196                   { 
00197   //                  std::cout << "new seed " << ip.second.value() << " " << ip.second.significance() << " " << track->hitPattern().trackerLayersWithMeasurement() << " " << track->pt() << " " << track->eta() << std::endl;
00198                     seeds.push_back(tt);  
00199                   }
00200  
00201         }
00202 
00203 
00204         //Create BS object from PV to feed in the AVR
00205         BeamSpot::CovarianceMatrix cov;
00206         for(unsigned int i = 0; i < 7; i++) {
00207                 for(unsigned int j = 0; j < 7; j++) {
00208                         if (i < 3 && j < 3)
00209                                 cov(i, j) = pv.covariance(i, j);
00210                         else
00211                                 cov(i, j) = 0.0;
00212                 }
00213         }
00214         BeamSpot bs(pv.position(), 0.0, 0.0, 0.0, 0.0, cov, BeamSpot::Unknown);
00215         std::auto_ptr<VertexCollection> recoVertices(new VertexCollection);
00216         int i=0;
00217         std::vector< std::vector<TransientTrack> > clusters;
00218         for(std::vector<TransientTrack>::const_iterator s = seeds.begin();
00219             s != seeds.end(); ++s,++i)
00220         {
00221                 
00222 //              std::cout << "Match pvd = "<< pvd[i] <<   std::endl;
00223                 std::pair<std::vector<reco::TransientTrack>,GlobalPoint>  ntracks = nearTracks(*s,tts,pv);
00224                 if(ntracks.first.size() == 0 ) continue;
00225                 ntracks.first.push_back(*s);
00226                 clusters.push_back(ntracks.first); 
00227                 std::vector<TransientVertex> vertices;
00228 //              try {
00229                         vertices = vtxReco->vertices(ntracks.first, bs);
00230                         TransientVertex singleFitVertex;
00231                         singleFitVertex = theAdaptiveFitter.vertex(ntracks.first,ntracks.second); //edPoint);
00232                         if(singleFitVertex.isValid())
00233                           vertices.push_back(singleFitVertex);
00234 //              } catch(...) {
00235 //                      vertices.clear();
00236 //              }
00237 
00238 //              std::cout << "for this seed I found  " << vertices.size() << " vertices"<< std::endl;
00239                 for(std::vector<TransientVertex>::const_iterator v = vertices.begin();
00240                     v != vertices.end(); ++v) {
00241 //                      if(v->degreesOfFreedom() > 0.2)
00242                         {
00243                          Measurement1D dlen= vdist.distance(pv,*v);
00244                          Measurement1D dlen2= vdist2d.distance(pv,*v);
00245                          reco::Vertex vv(*v);
00246   //                       std::cout << "V chi2/n: " << v->normalisedChiSquared() << " ndof: " <<v->degreesOfFreedom() ;
00247     //                     std::cout << " dlen: " << dlen.value() << " error: " << dlen.error() << " signif: " << dlen.significance();
00248       //                   std::cout << " dlen2: " << dlen2.value() << " error2: " << dlen2.error() << " signif2: " << dlen2.significance();
00249         //                 std::cout << " pos: " << vv.position() << " error: " <<vv.xError() << " " << vv.yError() << " " << vv.zError() << std::endl;
00250                          GlobalVector dir;  
00251                          std::vector<reco::TransientTrack> ts = v->originalTracks();
00252                         for(std::vector<reco::TransientTrack>::const_iterator i = ts.begin();
00253                             i != ts.end(); ++i) {
00254                                 reco::TrackRef t = i->trackBaseRef().castTo<reco::TrackRef>();
00255                                 float w = v->trackWeight(*i);
00256                                 if (w > 0.5) dir+=i->impactPointState().globalDirection();
00257           //                      std::cout << "\t[" << (*t).pt() << ": "
00258             //                              << (*t).eta() << ", "
00259               //                            << (*t).phi() << "], "
00260                 //                          << w << std::endl;
00261                         }
00262                        GlobalPoint ppv(pv.position().x(),pv.position().y(),pv.position().z());
00263                        GlobalPoint sv((*v).position().x(),(*v).position().y(),(*v).position().z());
00264                        float vscal = dir.unit().dot((sv-ppv).unit()) ;
00265 //                        std::cout << "Vscal: " <<  vscal << std::endl;
00266                         if(dlen.significance() > vertexMinDLenSig  && vscal > vertexMinAngleCosine &&  v->normalisedChiSquared() < 10 && dlen2.significance() > vertexMinDLen2DSig)
00267                          recoVertices->push_back(*v);
00268 
00269 
00270                         }
00271                    }
00272         }
00273         event.put(recoVertices);
00274 
00275 }
00276 
00277 DEFINE_FWK_MODULE(InclusiveVertexFinder);
/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/RecoVertex/AdaptiveVertexFinder/plugins/InclusiveVertexFinder.cc
