Reference Manual

Go to the documentation of this file.
00001 #include <functional>
00002 #include <algorithm>
00003 #include <iterator>
00004 #include <cstddef>
00005 #include <string>
00006 #include <vector>
00007 #include <map>
00008 #include <set>
00009 
00010 #include <boost/iterator/transform_iterator.hpp>
00011 
00012 #include "FWCore/Framework/interface/EDProducer.h"
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/Framework/interface/EventSetup.h"
00015 #include "FWCore/Framework/interface/MakerMacros.h"
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 #include "FWCore/Utilities/interface/InputTag.h"
00018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00019 #include "FWCore/Utilities/interface/Exception.h"
00020 
00021 #include "DataFormats/TrackReco/interface/Track.h"
00022 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00023 #include "DataFormats/VertexReco/interface/Vertex.h"
00024 
00025 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
00026 #include "DataFormats/BTauReco/interface/TrackIPTagInfo.h"
00027 #include "DataFormats/BTauReco/interface/SecondaryVertexTagInfo.h"
00028 
00029 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
00030 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
00031 #include "RecoVertex/ConfigurableVertexReco/interface/ConfigurableVertexReconstructor.h"
00032 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackVertexFinder.h"
00033 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackPrediction.h"
00034 #include "RecoVertex/GhostTrackFitter/interface/GhostTrackState.h"
00035 #include "RecoVertex/GhostTrackFitter/interface/GhostTrack.h"
00036 
00037 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
00038 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00039 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00040 
00041 #include "RecoBTag/SecondaryVertex/interface/TrackSelector.h"
00042 #include "RecoBTag/SecondaryVertex/interface/TrackSorting.h"
00043 #include "RecoBTag/SecondaryVertex/interface/SecondaryVertex.h"
00044 #include "RecoBTag/SecondaryVertex/interface/VertexFilter.h"
00045 #include "RecoBTag/SecondaryVertex/interface/VertexSorting.h"
00046 
00047 using namespace reco;
00048 
00049 namespace {
00050         template<typename T>
00051         struct RefToBaseLess : public std::binary_function<edm::RefToBase<T>,
00052                                                            edm::RefToBase<T>,
00053                                                            bool> {
00054                 inline bool operator()(const edm::RefToBase<T> &r1,
00055                                        const edm::RefToBase<T> &r2) const
00056                 {
00057                         return r1.id() < r2.id() ||
00058                                (r1.id() == r2.id() && r1.key() < r2.key());
00059                 }
00060         };
00061 }
00062 
00063 class SecondaryVertexProducer : public edm::EDProducer {
00064     public:
00065         explicit SecondaryVertexProducer(const edm::ParameterSet &params);
00066         ~SecondaryVertexProducer();
00067 
00068         virtual void produce(edm::Event &event, const edm::EventSetup &es);
00069 
00070     private:
00071         enum ConstraintType {
00072                 CONSTRAINT_NONE = 0,
00073                 CONSTRAINT_BEAMSPOT,
00074                 CONSTRAINT_PV_BEAMSPOT_SIZE,
00075                 CONSTRAINT_PV_BS_Z_ERRORS_SCALED,
00076                 CONSTRAINT_PV_ERROR_SCALED,
00077                 CONSTRAINT_PV_PRIMARIES_IN_FIT
00078         };
00079 
00080         static ConstraintType getConstraintType(const std::string &name);
00081 
00082         const edm::InputTag             trackIPTagInfoLabel;
00083         edm::InputTag beamSpotTag;
00084         TrackIPTagInfo::SortCriteria    sortCriterium;
00085         TrackSelector                   trackSelector;
00086         ConstraintType                  constraint;
00087         double                          constraintScaling;
00088         edm::ParameterSet               vtxRecoPSet;
00089         bool                            useGhostTrack;
00090         bool                            withPVError;
00091         double                          minTrackWeight;
00092         VertexFilter                    vertexFilter;
00093         VertexSorting                   vertexSorting;
00094 };
00095 
00096 SecondaryVertexProducer::ConstraintType
00097 SecondaryVertexProducer::getConstraintType(const std::string &name)
00098 {
00099         if (name == "None")
00100                 return CONSTRAINT_NONE;
00101         else if (name == "BeamSpot")
00102                 return CONSTRAINT_BEAMSPOT;
00103         else if (name == "BeamSpot+PVPosition")
00104                 return CONSTRAINT_PV_BEAMSPOT_SIZE;
00105         else if (name == "BeamSpotZ+PVErrorScaledXY")
00106                 return CONSTRAINT_PV_BS_Z_ERRORS_SCALED;
00107         else if (name == "PVErrorScaled")
00108                 return CONSTRAINT_PV_ERROR_SCALED;
00109         else if (name == "BeamSpot+PVTracksInFit")
00110                 return CONSTRAINT_PV_PRIMARIES_IN_FIT;
00111         else
00112                 throw cms::Exception("InvalidArgument")
00113                         << "SecondaryVertexProducer: ``constraint'' parameter "
00114                            "value \"" << name << "\" not understood."
00115                         << std::endl;
00116 }
00117 
00118 static GhostTrackVertexFinder::FitType
00119 getGhostTrackFitType(const std::string &name)
00120 {
00121         if (name == "AlwaysWithGhostTrack")
00122                 return GhostTrackVertexFinder::kAlwaysWithGhostTrack;
00123         else if (name == "SingleTracksWithGhostTrack")
00124                 return GhostTrackVertexFinder::kSingleTracksWithGhostTrack;
00125         else if (name == "RefitGhostTrackWithVertices")
00126                 return GhostTrackVertexFinder::kRefitGhostTrackWithVertices;
00127         else
00128                 throw cms::Exception("InvalidArgument")
00129                         << "SecondaryVertexProducer: ``fitType'' "
00130                            "parameter value \"" << name << "\" for "
00131                            "GhostTrackVertexFinder settings not "
00132                            "understood." << std::endl;
00133 }
00134 
00135 SecondaryVertexProducer::SecondaryVertexProducer(
00136                                         const edm::ParameterSet &params) :
00137         trackIPTagInfoLabel(params.getParameter<edm::InputTag>("trackIPTagInfos")),
00138         sortCriterium(TrackSorting::getCriterium(params.getParameter<std::string>("trackSort"))),
00139         trackSelector(params.getParameter<edm::ParameterSet>("trackSelection")),
00140         constraint(getConstraintType(params.getParameter<std::string>("constraint"))),
00141         constraintScaling(1.0),
00142         vtxRecoPSet(params.getParameter<edm::ParameterSet>("vertexReco")),
00143         useGhostTrack(vtxRecoPSet.getParameter<std::string>("finder") == "gtvr"),
00144         withPVError(params.getParameter<bool>("usePVError")),
00145         minTrackWeight(params.getParameter<double>("minimumTrackWeight")),
00146         vertexFilter(params.getParameter<edm::ParameterSet>("vertexCuts")),
00147         vertexSorting(params.getParameter<edm::ParameterSet>("vertexSelection"))
00148 {
00149         if (constraint == CONSTRAINT_PV_ERROR_SCALED ||
00150             constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED)
00151                 constraintScaling = params.getParameter<double>("pvErrorScaling");
00152 
00153         if (constraint == CONSTRAINT_PV_BEAMSPOT_SIZE ||
00154             constraint == CONSTRAINT_PV_BS_Z_ERRORS_SCALED ||
00155             constraint == CONSTRAINT_BEAMSPOT ||
00156             constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT )
00157             beamSpotTag = params.getParameter<edm::InputTag>("beamSpotTag");
00158 
00159         produces<SecondaryVertexTagInfoCollection>();
00160 }
00161 
00162 SecondaryVertexProducer::~SecondaryVertexProducer()
00163 {
00164 }
00165 
00166 namespace {
00167         struct SVBuilder :
00168                 public std::unary_function<const Vertex&, SecondaryVertex> {
00169 
00170                 SVBuilder(const reco::Vertex &pv,
00171                           const GlobalVector &direction,
00172                           bool withPVError) :
00173                         pv(pv), direction(direction),
00174                         withPVError(withPVError) {}
00175 
00176                 SecondaryVertex operator () (const reco::Vertex &sv) const
00177                 { return SecondaryVertex(pv, sv, direction, withPVError); }
00178 
00179                 const Vertex            &pv;
00180                 const GlobalVector      &direction;
00181                 bool                    withPVError;
00182         };
00183 
00184         struct SVFilter :
00185                 public std::unary_function<const SecondaryVertex&, bool> {
00186 
00187                 SVFilter(const VertexFilter &filter, const Vertex &pv,
00188                          const GlobalVector &direction) :
00189                         filter(filter), pv(pv), direction(direction) {}
00190 
00191                 inline bool operator () (const SecondaryVertex &sv) const
00192                 { return !filter(pv, sv, direction); }
00193 
00194                 const VertexFilter      &filter;
00195                 const Vertex            &pv;
00196                 const GlobalVector      &direction;
00197         };
00198                         
00199 } // anonynmous namespace
00200 
00201 void SecondaryVertexProducer::produce(edm::Event &event,
00202                                       const edm::EventSetup &es)
00203 {
00204         typedef std::map<TrackBaseRef, TransientTrack,
00205                          RefToBaseLess<Track> > TransientTrackMap;
00206 
00207         edm::ESHandle<TransientTrackBuilder> trackBuilder;
00208         es.get<TransientTrackRecord>().get("TransientTrackBuilder",
00209                                            trackBuilder);
00210 
00211         edm::Handle<TrackIPTagInfoCollection> trackIPTagInfos;
00212         event.getByLabel(trackIPTagInfoLabel, trackIPTagInfos);
00213 
00214         edm::Handle<BeamSpot> beamSpot;
00215         unsigned int bsCovSrc[7] = { 0, };
00216         double sigmaZ = 0.0, beamWidth = 0.0;
00217         switch(constraint) {
00218             case CONSTRAINT_PV_BEAMSPOT_SIZE:
00219                 event.getByLabel(beamSpotTag,beamSpot);
00220                 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = bsCovSrc[6] = 1;
00221                 sigmaZ = beamSpot->sigmaZ();
00222                 beamWidth = beamSpot->BeamWidthX();
00223                 break;
00224 
00225             case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
00226               event.getByLabel(beamSpotTag,beamSpot);
00227                 bsCovSrc[0] = bsCovSrc[1] = 2;
00228                 bsCovSrc[3] = bsCovSrc[4] = bsCovSrc[5] = 1;
00229                 sigmaZ = beamSpot->sigmaZ();
00230                 break;
00231 
00232             case CONSTRAINT_PV_ERROR_SCALED:
00233                 bsCovSrc[0] = bsCovSrc[1] = bsCovSrc[2] = 2;
00234                 break;
00235 
00236             case CONSTRAINT_BEAMSPOT:
00237             case CONSTRAINT_PV_PRIMARIES_IN_FIT:
00238                 event.getByLabel(beamSpotTag,beamSpot);
00239                 break;
00240 
00241             default:
00242                 /* nothing */;
00243         }
00244 
00245         std::auto_ptr<ConfigurableVertexReconstructor> vertexReco;
00246         std::auto_ptr<GhostTrackVertexFinder> vertexRecoGT;
00247         if (useGhostTrack)
00248                 vertexRecoGT.reset(new GhostTrackVertexFinder(
00249                         vtxRecoPSet.getParameter<double>("maxFitChi2"),
00250                         vtxRecoPSet.getParameter<double>("mergeThreshold"),
00251                         vtxRecoPSet.getParameter<double>("primcut"),
00252                         vtxRecoPSet.getParameter<double>("seccut"),
00253                         getGhostTrackFitType(vtxRecoPSet.getParameter<std::string>("fitType"))));
00254         else
00255                 vertexReco.reset(
00256                         new ConfigurableVertexReconstructor(vtxRecoPSet));
00257 
00258         TransientTrackMap primariesMap;
00259 
00260         // result secondary vertices
00261 
00262         std::auto_ptr<SecondaryVertexTagInfoCollection>
00263                         tagInfos(new SecondaryVertexTagInfoCollection);
00264 
00265         for(TrackIPTagInfoCollection::const_iterator iterJets =
00266                 trackIPTagInfos->begin(); iterJets != trackIPTagInfos->end();
00267                 ++iterJets) {
00268 
00269                 std::vector<SecondaryVertexTagInfo::IndexedTrackData> trackData;
00270 
00271                 const Vertex &pv = *iterJets->primaryVertex();
00272 
00273                 std::set<TransientTrack> primaries;
00274                 if (constraint == CONSTRAINT_PV_PRIMARIES_IN_FIT) {
00275                         for(Vertex::trackRef_iterator iter = pv.tracks_begin();
00276                             iter != pv.tracks_end(); ++iter) {
00277                                 TransientTrackMap::iterator pos =
00278                                         primariesMap.lower_bound(*iter);
00279 
00280                                 if (pos != primariesMap.end() &&
00281                                     pos->first == *iter)
00282                                         primaries.insert(pos->second);
00283                                 else {
00284                                         TransientTrack track =
00285                                                 trackBuilder->build(
00286                                                         iter->castTo<TrackRef>());
00287                                         primariesMap.insert(pos,
00288                                                 std::make_pair(*iter, track));
00289                                         primaries.insert(track);
00290                                 }
00291                         }
00292                 }
00293 
00294                 edm::RefToBase<Jet> jetRef = iterJets->jet();
00295 
00296                 GlobalVector jetDir(jetRef->momentum().x(),
00297                                     jetRef->momentum().y(),
00298                                     jetRef->momentum().z());
00299 
00300                 std::vector<std::size_t> indices =
00301                                 iterJets->sortedIndexes(sortCriterium);
00302 
00303                 TrackRefVector trackRefs = iterJets->sortedTracks(indices);
00304 
00305                 const std::vector<TrackIPTagInfo::TrackIPData> &ipData =
00306                                         iterJets->impactParameterData();
00307 
00308                 // build transient tracks used for vertex reconstruction
00309 
00310                 std::vector<TransientTrack> fitTracks;
00311                 std::vector<GhostTrackState> gtStates;
00312                 std::auto_ptr<GhostTrackPrediction> gtPred;
00313                 if (useGhostTrack)
00314                         gtPred.reset(new GhostTrackPrediction(
00315                                                 *iterJets->ghostTrack()));
00316 
00317                 for(unsigned int i = 0; i < indices.size(); i++) {
00318                         typedef SecondaryVertexTagInfo::IndexedTrackData IndexedTrackData;
00319 
00320                         const TrackRef &trackRef = trackRefs[i];
00321 
00322                         trackData.push_back(IndexedTrackData());
00323                         trackData.back().first = indices[i];
00324 
00325                         // select tracks for SV finder
00326 
00327                         if (!trackSelector(*trackRef, ipData[i], *jetRef,
00328                                            RecoVertex::convertPos(
00329                                                         pv.position()))) {
00330                                 trackData.back().second.svStatus =
00331                                         SecondaryVertexTagInfo::TrackData::trackSelected;
00332                                 continue;
00333                         }
00334 
00335                         TransientTrackMap::const_iterator pos =
00336                                         primariesMap.find(
00337                                                 TrackBaseRef(trackRef));
00338                         TransientTrack fitTrack;
00339                         if (pos != primariesMap.end()) {
00340                                 primaries.erase(pos->second);
00341                                 fitTrack = pos->second;
00342                         } else
00343                                 fitTrack = trackBuilder->build(trackRef);
00344                         fitTracks.push_back(fitTrack);
00345 
00346                         trackData.back().second.svStatus =
00347                                 SecondaryVertexTagInfo::TrackData::trackUsedForVertexFit;
00348 
00349                         if (useGhostTrack) {
00350                                 GhostTrackState gtState(fitTrack);
00351                                 GlobalPoint pos =
00352                                         ipData[i].closestToGhostTrack;
00353                                 gtState.linearize(*gtPred, true,
00354                                                   gtPred->lambda(pos));
00355                                 gtState.setWeight(ipData[i].ghostTrackWeight);
00356                                 gtStates.push_back(gtState);
00357                         }
00358                 }
00359 
00360                 std::auto_ptr<GhostTrack> ghostTrack;
00361                 if (useGhostTrack)
00362                         ghostTrack.reset(new GhostTrack(
00363                                 GhostTrackPrediction(
00364                                         RecoVertex::convertPos(pv.position()),
00365                                         RecoVertex::convertError(pv.error()),
00366                                         GlobalVector(
00367                                                 iterJets->ghostTrack()->px(),
00368                                                 iterJets->ghostTrack()->py(),
00369                                                 iterJets->ghostTrack()->pz()),
00370                                         0.05),
00371                                 *gtPred, gtStates,
00372                                 iterJets->ghostTrack()->chi2(),
00373                                 iterJets->ghostTrack()->ndof()));
00374 
00375                 // perform actual vertex finding
00376 
00377                 std::vector<TransientVertex> fittedSVs;
00378                 switch(constraint) {
00379                     case CONSTRAINT_NONE:
00380                         if (useGhostTrack)
00381                                 fittedSVs = vertexRecoGT->vertices(
00382                                                 pv, *ghostTrack);
00383                         else
00384                                 fittedSVs = vertexReco->vertices(fitTracks);
00385                         break;
00386 
00387                     case CONSTRAINT_BEAMSPOT:
00388                         if (useGhostTrack)
00389                                 fittedSVs = vertexRecoGT->vertices(
00390                                                 pv, *beamSpot, *ghostTrack);
00391                         else
00392                                 fittedSVs = vertexReco->vertices(fitTracks,
00393                                                                  *beamSpot);
00394                         break;
00395 
00396                     case CONSTRAINT_PV_BEAMSPOT_SIZE:
00397                     case CONSTRAINT_PV_BS_Z_ERRORS_SCALED:
00398                     case CONSTRAINT_PV_ERROR_SCALED: {
00399                         BeamSpot::CovarianceMatrix cov;
00400                         for(unsigned int i = 0; i < 7; i++) {
00401                                 unsigned int covSrc = bsCovSrc[i];
00402                                 for(unsigned int j = 0; j < 7; j++) {
00403                                         double v;
00404                                         if (!covSrc || bsCovSrc[j] != covSrc)
00405                                                 v = 0.0;
00406                                         else if (covSrc == 1)
00407                                                 v = beamSpot->covariance(i, j);
00408                                         else
00409                                                 v = pv.covariance(i, j) *
00410                                                     constraintScaling;
00411                                         cov(i, j) = v;
00412                                 }
00413                         }
00414 
00415                         BeamSpot bs(pv.position(), sigmaZ,
00416                                     beamSpot.isValid() ? beamSpot->dxdz() : 0.,
00417                                     beamSpot.isValid() ? beamSpot->dydz() : 0.,
00418                                     beamWidth, cov, BeamSpot::Unknown);
00419 
00420                         if (useGhostTrack)
00421                                 fittedSVs = vertexRecoGT->vertices(
00422                                                 pv, bs, *ghostTrack);
00423                         else
00424                                 fittedSVs = vertexReco->vertices(fitTracks, bs);
00425                     }   break;
00426 
00427                     case CONSTRAINT_PV_PRIMARIES_IN_FIT: {
00428                         std::vector<TransientTrack> primaries_(
00429                                         primaries.begin(), primaries.end());
00430                         if (useGhostTrack)
00431                                 fittedSVs = vertexRecoGT->vertices(
00432                                                 pv, *beamSpot, primaries_,
00433                                                 *ghostTrack);
00434                         else
00435                                 fittedSVs = vertexReco->vertices(
00436                                                 primaries_, fitTracks,
00437                                                 *beamSpot);
00438                     }   break;
00439                 }
00440 
00441                 // build combined SV information and filter
00442 
00443                 std::vector<SecondaryVertex> SVs;
00444                 SVBuilder svBuilder(pv, jetDir, withPVError);
00445                 std::remove_copy_if(boost::make_transform_iterator(
00446                                         fittedSVs.begin(), svBuilder),
00447                                     boost::make_transform_iterator(
00448                                         fittedSVs.end(), svBuilder),
00449                                     std::back_inserter(SVs),
00450                                     SVFilter(vertexFilter, pv, jetDir));
00451 
00452                 // clean up now unneeded collections
00453 
00454                 gtPred.reset();
00455                 ghostTrack.reset();
00456                 gtStates.clear();
00457                 fitTracks.clear();
00458                 fittedSVs.clear();
00459 
00460                 // sort SVs by importance
00461 
00462                 std::vector<unsigned int> vtxIndices = vertexSorting(SVs);
00463 
00464                 std::vector<SecondaryVertexTagInfo::VertexData> svData;
00465 
00466                 svData.resize(vtxIndices.size());
00467                 for(unsigned int idx = 0; idx < vtxIndices.size(); idx++) {
00468                         const SecondaryVertex &sv = SVs[vtxIndices[idx]];
00469 
00470                         svData[idx].vertex = sv;
00471                         svData[idx].dist2d = sv.dist2d();
00472                         svData[idx].dist3d = sv.dist3d();
00473                         svData[idx].direction =
00474                                 GlobalVector(sv.x() - pv.x(),
00475                                              sv.y() - pv.y(),
00476                                              sv.z() - pv.z());
00477 
00478                         // mark tracks successfully used in vertex fit
00479 
00480                         for(Vertex::trackRef_iterator iter = sv.tracks_begin();
00481                             iter != sv.tracks_end(); iter++) {
00482                                 if (sv.trackWeight(*iter) < minTrackWeight)
00483                                         continue;
00484 
00485                                 TrackRefVector::const_iterator pos =
00486                                         std::find(trackRefs.begin(), trackRefs.end(),
00487                                                   iter->castTo<TrackRef>());
00488                                 if (pos == trackRefs.end())
00489                                         throw cms::Exception("TrackNotFound")
00490                                                 << "Could not find track from secondary "
00491                                                    "vertex in original tracks."
00492                                                 << std::endl;
00493 
00494                                 unsigned int index = pos - trackRefs.begin();
00495                                 trackData[index].second.svStatus =
00496                                         (SecondaryVertexTagInfo::TrackData::Status)
00497                                         ((unsigned int)SecondaryVertexTagInfo::TrackData::trackAssociatedToVertex + idx);
00498                         }
00499                 }
00500 
00501                 // fill result into tag infos
00502 
00503                 tagInfos->push_back(
00504                         SecondaryVertexTagInfo(
00505                                 trackData, svData, SVs.size(),
00506                                 TrackIPTagInfoRef(trackIPTagInfos,
00507                                         iterJets - trackIPTagInfos->begin())));
00508         }
00509 
00510         event.put(tagInfos);
00511 }
00512 
00513 //define this as a plug-in
00514 DEFINE_FWK_MODULE(SecondaryVertexProducer);
/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoBTag/SecondaryVertex/plugins/SecondaryVertexProducer.cc
