diff --git a/BTagHelpers/interface/TrackPairInfoBuilder.h b/BTagHelpers/interface/TrackPairInfoBuilder.h new file mode 100644 index 0000000..234d1dd --- /dev/null +++ b/BTagHelpers/interface/TrackPairInfoBuilder.h @@ -0,0 +1,100 @@ +#ifndef BTAGHELPERS_INTERFACE_TRACKPARINFOBUILDER_H_ +#define BTAGHELPERS_INTERFACE_TRACKPAIRINFOBUILDER_H_ + +#include "DataFormats/GeometrySurface/interface/Line.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "TrackingTools/PatternTools/interface/TwoTrackMinimumDistance.h" +#include "DataFormats/PatCandidates/interface/Jet.h" + +namespace deepntuples { + + class TrackPairInfoBuilder { + public: + TrackPairInfoBuilder(); + + void buildTrackPairInfo(const reco::TransientTrack* it, + const reco::TransientTrack* tt, + const reco::Vertex& pv, + const pat::Jet &jet + ); + + const float track_pt() const { return track_pt_; } + const float track_eta() const { return track_eta_; } + const float track_phi() const { return track_phi_; } + const float track_dz() const { return track_dz_; } + const float track_dxy() const { return track_dxy_; } + const float pca_distance() const { return pca_distance_; } + const float pca_significance() const { return pca_significance_; } + const float pcaSeed_x() const { return pcaSeed_x_; } + const float pcaSeed_y() const { return pcaSeed_y_; } + const float pcaSeed_z() const { return pcaSeed_z_; } + const float pcaSeed_xerr() const { return pcaSeed_xerr_; } + const float pcaSeed_yerr() const { return pcaSeed_yerr_; } + const float pcaSeed_zerr() const { return pcaSeed_zerr_; } + const float pcaTrack_x() const { return pcaTrack_x_; } + const float pcaTrack_y() const { return pcaTrack_y_; } + const float pcaTrack_z() const { return pcaTrack_z_; } + const float pcaTrack_xerr() const { return pcaTrack_xerr_; } + const float pcaTrack_yerr() const { return pcaTrack_yerr_; } + const float pcaTrack_zerr() const { return pcaTrack_zerr_; } + const float dotprodTrack() const { return dotprodTrack_; } + const float dotprodSeed() const { return dotprodSeed_; } + const float pcaSeed_dist() const { return pcaSeed_dist_; } + const float pcaTrack_dist() const { return pcaTrack_dist_; } + const float track_candMass() const { return track_candMass_; } + const float track_ip2d() const { return track_ip2d_; } + const float track_ip2dSig() const { return track_ip2dSig_; } + const float track_ip3d() const { return track_ip3d_; } + const float track_ip3dSig() const { return track_ip3dSig_; } + const float dotprodTrackSeed2D() const { return dotprodTrackSeed2D_; } + const float dotprodTrackSeed2DV() const { return dotprodTrackSeed2DV_; } + const float dotprodTrackSeed3D() const { return dotprodTrackSeed3D_; } + const float dotprodTrackSeed3DV() const { return dotprodTrackSeed3DV_; } + const float pca_jetAxis_dist() const { return pca_jetAxis_dist_; } + const float pca_jetAxis_dotprod() const { return pca_jetAxis_dotprod_; } + const float pca_jetAxis_dEta() const { return pca_jetAxis_dEta_; } + const float pca_jetAxis_dPhi() const { return pca_jetAxis_dPhi_; } + + private: + float track_pt_; + float track_eta_; + float track_phi_; + float track_dz_; + float track_dxy_; + float pca_distance_; + float pca_significance_; + float pcaSeed_x_; + float pcaSeed_y_; + float pcaSeed_z_; + float pcaSeed_xerr_; + float pcaSeed_yerr_; + float pcaSeed_zerr_; + float pcaTrack_x_; + float pcaTrack_y_; + float pcaTrack_z_; + float pcaTrack_xerr_; + float pcaTrack_yerr_; + float pcaTrack_zerr_; + float dotprodTrack_; + float dotprodSeed_; + float pcaSeed_dist_; + float pcaTrack_dist_; + float track_candMass_; + float track_ip2d_; + float track_ip2dSig_; + float track_ip3d_; + float track_ip3dSig_; + float dotprodTrackSeed2D_; + float dotprodTrackSeed2DV_; + float dotprodTrackSeed3D_; + float dotprodTrackSeed3DV_; + float pca_jetAxis_dist_; + float pca_jetAxis_dotprod_; + float pca_jetAxis_dEta_; + float pca_jetAxis_dPhi_; + }; + +} // namespace btagbtvdeep + +#endif //RecoBTag_FeatureTools_TrackPairInfoBuilder_h diff --git a/BTagHelpers/src/TrackPairInfoBuilder.cc b/BTagHelpers/src/TrackPairInfoBuilder.cc new file mode 100644 index 0000000..3fbe15a --- /dev/null +++ b/BTagHelpers/src/TrackPairInfoBuilder.cc @@ -0,0 +1,154 @@ +#include "DataFormats/Candidate/interface/Candidate.h" +#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TVector3.h" +#include "DataFormats/PatCandidates/interface/Jet.h" + +#include "DeepNTuples/NtupleCommons/interface/InfinityCatcher.h" +#include "DeepNTuples/BTagHelpers/interface/TrackPairInfoBuilder.h" + +namespace deepntuples { + + TrackPairInfoBuilder::TrackPairInfoBuilder() + : + + track_pt_(0), + track_eta_(0), + track_phi_(0), + track_dz_(0), + track_dxy_(0), + + pca_distance_(0), + pca_significance_(0), + + pcaSeed_x_(0), + pcaSeed_y_(0), + pcaSeed_z_(0), + pcaSeed_xerr_(0), + pcaSeed_yerr_(0), + pcaSeed_zerr_(0), + pcaTrack_x_(0), + pcaTrack_y_(0), + pcaTrack_z_(0), + pcaTrack_xerr_(0), + pcaTrack_yerr_(0), + pcaTrack_zerr_(0), + + dotprodTrack_(0), + dotprodSeed_(0), + pcaSeed_dist_(0), + pcaTrack_dist_(0), + + track_candMass_(0), + track_ip2d_(0), + track_ip2dSig_(0), + track_ip3d_(0), + track_ip3dSig_(0), + + dotprodTrackSeed2D_(0), + dotprodTrackSeed2DV_(0), + dotprodTrackSeed3D_(0), + dotprodTrackSeed3DV_(0), + + pca_jetAxis_dist_(0), + pca_jetAxis_dotprod_(0), + pca_jetAxis_dEta_(0), + pca_jetAxis_dPhi_(0) + + {} + + void TrackPairInfoBuilder::buildTrackPairInfo(const reco::TransientTrack* it, + const reco::TransientTrack* tt, + const reco::Vertex& pv, + const pat::Jet &jet + ) { + + GlobalVector jetdirection(jet.px(), jet.py(), jet.pz()); + GlobalPoint pvp(pv.x(), pv.y(), pv.z()); + + VertexDistance3D distanceComputer; + TwoTrackMinimumDistance dist; + + auto const& iImpactState = it->impactPointState(); + auto const& tImpactState = tt->impactPointState(); + + if (dist.calculate(tImpactState, iImpactState)) { + GlobalPoint ttPoint = dist.points().first; + GlobalError ttPointErr = tImpactState.cartesianError().position(); + GlobalPoint seedPosition = dist.points().second; + GlobalError seedPositionErr = iImpactState.cartesianError().position(); + + Measurement1D m = + distanceComputer.distance(VertexState(seedPosition, seedPositionErr), VertexState(ttPoint, ttPointErr)); + + GlobalPoint cp(dist.crossingPoint()); + + GlobalVector pairMomentum((Basic3DVector)(it->track().momentum() + tt->track().momentum())); + GlobalVector pvToPCA(cp - pvp); + + float pvToPCAseed = (seedPosition - pvp).mag(); + float pvToPCAtrack = (ttPoint - pvp).mag(); + float distance = dist.distance(); + + GlobalVector trackDir2D(tImpactState.globalDirection().x(), tImpactState.globalDirection().y(), 0.); + GlobalVector seedDir2D(iImpactState.globalDirection().x(), iImpactState.globalDirection().y(), 0.); + GlobalVector trackPCADir2D(ttPoint.x() - pvp.x(), ttPoint.y() - pvp.y(), 0.); + GlobalVector seedPCADir2D(seedPosition.x() - pvp.x(), seedPosition.y() - pvp.y(), 0.); + + float dotprodTrack = (ttPoint - pvp).unit().dot(tImpactState.globalDirection().unit()); + float dotprodSeed = (seedPosition - pvp).unit().dot(iImpactState.globalDirection().unit()); + + Line::PositionType pos(pvp); + Line::DirectionType dir(jetdirection); + Line::DirectionType pairMomentumDir(pairMomentum); + Line jetLine(pos, dir); + Line PCAMomentumLine(cp, pairMomentumDir); + +// track_pt_ = tt->track().pt(); +// track_eta_ = tt->track().eta(); +// track_phi_ = tt->track().phi(); +// track_dz_ = tt->track().dz(pv.position()); +// track_dxy_ = tt->track().dxy(pv.position()); + + pca_distance_ = distance; + pca_significance_ = m.significance(); + + pcaSeed_x_ = seedPosition.x(); + pcaSeed_y_ = seedPosition.y(); + pcaSeed_z_ = seedPosition.z(); + pcaSeed_xerr_ = seedPositionErr.cxx(); + pcaSeed_yerr_ = seedPositionErr.cyy(); + pcaSeed_zerr_ = seedPositionErr.czz(); + pcaTrack_x_ = ttPoint.x(); + pcaTrack_y_ = ttPoint.y(); + pcaTrack_z_ = ttPoint.z(); + pcaTrack_xerr_ = ttPointErr.cxx(); + pcaTrack_yerr_ = ttPointErr.cyy(); + pcaTrack_zerr_ = ttPointErr.czz(); + + dotprodTrack_ = dotprodTrack; + dotprodSeed_ = dotprodSeed; + pcaSeed_dist_ = pvToPCAseed; + pcaTrack_dist_ = pvToPCAtrack; + +// track_candMass_ = mass; +// track_ip2d_ = t_ip2d.second.value(); +// track_ip2dSig_ = t_ip2d.second.significance(); +// track_ip3d_ = t_ip.second.value(); +// track_ip3dSig_ = t_ip.second.significance(); + + dotprodTrackSeed2D_ = trackDir2D.unit().dot(seedDir2D.unit()); + dotprodTrackSeed3D_ = iImpactState.globalDirection().unit().dot(tImpactState.globalDirection().unit()); + dotprodTrackSeed2DV_ = trackPCADir2D.unit().dot(seedPCADir2D.unit()); + dotprodTrackSeed3DV_ = (seedPosition - pvp).unit().dot((ttPoint - pvp).unit()); + + pca_jetAxis_dist_ = jetLine.distance(cp).mag(); + pca_jetAxis_dotprod_ = pairMomentum.unit().dot(jetdirection.unit()); + pca_jetAxis_dEta_ = std::fabs(pvToPCA.eta() - jetdirection.eta()); + pca_jetAxis_dPhi_ = std::fabs(pvToPCA.phi() - jetdirection.phi()); + } + } + +} // namespace btagbtvdeep diff --git a/Ntupler/interface/TrackPairFiller.h b/Ntupler/interface/TrackPairFiller.h new file mode 100644 index 0000000..f3c801f --- /dev/null +++ b/Ntupler/interface/TrackPairFiller.h @@ -0,0 +1,43 @@ +#ifndef NTUPLER_INTERFACE_TRACKPAIRFILLER_H_ +#define NTUPLER_INTERFACE_TRACKPAIRFILLER_H_ + +#include "DeepNTuples/BTagHelpers/interface/TrackPairInfoBuilder.h" +#include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h" +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "DeepNTuples/NtupleCommons/interface/NtupleBase.h" + +namespace deepntuples { + +class TrackPairFiller: public NtupleBase { +public: + TrackPairFiller() : TrackPairFiller("", 0.8) {} + TrackPairFiller(std::string branchName, double jetR=0.8) : NtupleBase(branchName, jetR) {} + virtual ~TrackPairFiller() {} + + // get input parameters from the cfg file + virtual void readConfig(const edm::ParameterSet& iConfig, edm::ConsumesCollector && cc) override; + + // read event content or event setup for each event + virtual void readEvent(const edm::Event& iEvent, const edm::EventSetup& iSetup) override; + +protected: + // declare the data branches (name, type, default values) + virtual void book() override; + // fill the branches + virtual bool fill(const pat::Jet &jet, size_t jetidx, const JetHelper &jet_helper) override; + +private: + edm::EDGetTokenT vtxToken_; + edm::Handle vertices; + + edm::EDGetTokenT svToken_; + edm::Handle SVs; + + edm::ESHandle builder_; +}; + +} /* namespace deepntuples */ + +#endif /* NTUPLER_INTERFACE_TRACKPAIRFILLER_H_ */ diff --git a/Ntupler/plugins/DeepNtuplizer.cc b/Ntupler/plugins/DeepNtuplizer.cc index 2c7ced8..39ce92b 100644 --- a/Ntupler/plugins/DeepNtuplizer.cc +++ b/Ntupler/plugins/DeepNtuplizer.cc @@ -20,7 +20,7 @@ #include "DeepNTuples/Ntupler/interface/FatJetInfoFiller.h" #include "DeepNTuples/Ntupler/interface/SVFiller.h" #include "DeepNTuples/Ntupler/interface/PFCompleteFiller.h" - +#include "DeepNTuples/Ntupler/interface/TrackPairFiller.h" using namespace deepntuples; @@ -74,6 +74,9 @@ DeepNtuplizer::DeepNtuplizer(const edm::ParameterSet& iConfig): PFCompleteFiller *parts = new PFCompleteFiller("", jetR); addModule(parts); + + TrackPairFiller *pairmatrix = new TrackPairFiller("", jetR); + addModule(pairmatrix); // read config and init modules for(auto& m: modules_) diff --git a/Ntupler/src/TrackPairFiller.cc b/Ntupler/src/TrackPairFiller.cc new file mode 100644 index 0000000..a4736d6 --- /dev/null +++ b/Ntupler/src/TrackPairFiller.cc @@ -0,0 +1,170 @@ +#include +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" + +#include "DeepNTuples/Ntupler/interface/TrackPairFiller.h" + +namespace deepntuples { + +void TrackPairFiller::readConfig(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& cc) { + vtxToken_ = cc.consumes(iConfig.getParameter("vertices")); + svToken_ = cc.consumes(iConfig.getParameter("SVs")); +} + +void TrackPairFiller::readEvent(const edm::Event& iEvent, const edm::EventSetup& iSetup) { + iEvent.getByToken(vtxToken_, vertices); + iEvent.getByToken(svToken_, SVs); + iSetup.get().get("TransientTrackBuilder", builder_); +} + +void TrackPairFiller::book() { + + //test + data.add("n_pfcandidates", 0); + + data.addMulti("track1_index"); + data.addMulti("track2_index"); + + // pca (1 to 2) + data.addMulti("pca_distance"); + data.addMulti("pca_significance"); + + // pcas (1 and 2) poistions + data.addMulti("pcaSeed_x1"); + data.addMulti("pcaSeed_y1"); + data.addMulti("pcaSeed_z1"); + + data.addMulti("pcaSeed_x2"); + data.addMulti("pcaSeed_y2"); + data.addMulti("pcaSeed_z2"); + + data.addMulti("pcaSeed_xerr1"); + data.addMulti("pcaSeed_yerr1"); + data.addMulti("pcaSeed_zerr1"); + + data.addMulti("pcaSeed_xerr2"); + data.addMulti("pcaSeed_yerr2"); + data.addMulti("pcaSeed_zerr2"); + + // dot prod betweeen track and pca direction + data.addMulti("dotprod1"); + data.addMulti("dotprod2"); + + //pca distance form PV on both tracks + data.addMulti("pca_dist1"); + data.addMulti("pca_dist2"); + + //track track or dir dir dotprod + data.addMulti("dotprod12_2D"); + data.addMulti("dotprod12_2DV"); + data.addMulti("dotprod12_3D"); + data.addMulti("dotprod12_3DV"); + + //jet pca relative + data.addMulti("pca_jetAxis_dist"); + data.addMulti("pca_jetAxis_dotprod"); + data.addMulti("pca_jetAxis_dEta"); + data.addMulti("pca_jetAxis_dPhi_"); + + + +} + +bool TrackPairFiller::fill(const pat::Jet& jet, size_t jetidx, const JetHelper& jet_helper) { + + const auto& pfCands = jet_helper.getJetConstituents(); + + data.fill("n_pfcandidates", pfCands.size()); + + std::vector selectedTracks; + std::vector selectedTracks_pfidx; + int counter = 0; + + + for (const auto& cand : pfCands){ + + //const auto *packed_cand = dynamic_cast(&(*cand)); + + if (cand->bestTrack()){ + selectedTracks.push_back( builder_->build(cand) ); + selectedTracks_pfidx.push_back(counter); + + + } + counter++; + + } + + + for(std::vector::const_iterator it = selectedTracks.begin(); it != selectedTracks.end(); it++){ + + int index1 = it - selectedTracks.begin(); + + for(std::vector::const_iterator tt = selectedTracks.begin(); tt != selectedTracks.end(); tt++){ + + int index2 = tt - selectedTracks.begin(); + + std::cout << index1 << index2 << std::endl; + + if (index1!=index2){ + + data.fillMulti("track1_index", selectedTracks_pfidx[index1]); + data.fillMulti("track2_index", selectedTracks_pfidx[index2]); + + TrackPairInfoBuilder trkpairinfo; + + trkpairinfo.buildTrackPairInfo(&(*it),&(*tt),vertices->at(0),jet); + + data.fillMulti("pca_distance", trkpairinfo.pca_distance()); + data.fillMulti("pca_significance", trkpairinfo.pca_significance()); + + // pcas (1 and 2) poistions + data.fillMulti("pcaSeed_x1", trkpairinfo.pcaSeed_x()); + data.fillMulti("pcaSeed_y1", trkpairinfo.pcaSeed_y()); + data.fillMulti("pcaSeed_z1", trkpairinfo.pcaSeed_z()); + + data.fillMulti("pcaSeed_x2", trkpairinfo.pcaTrack_x()); + data.fillMulti("pcaSeed_y2", trkpairinfo.pcaTrack_y()); + data.fillMulti("pcaSeed_z2", trkpairinfo.pcaTrack_z()); + + data.fillMulti("pcaSeed_xerr1", trkpairinfo.pcaSeed_xerr()); + data.fillMulti("pcaSeed_yerr1", trkpairinfo.pcaSeed_yerr()); + data.fillMulti("pcaSeed_zerr1", trkpairinfo.pcaSeed_zerr()); + + data.fillMulti("pcaSeed_xerr2", trkpairinfo.pcaTrack_xerr()); + data.fillMulti("pcaSeed_yerr2", trkpairinfo.pcaTrack_yerr()); + data.fillMulti("pcaSeed_zerr2", trkpairinfo.pcaTrack_zerr()); + + // dot prod betweeen track and pca direction + data.fillMulti("dotprod1", trkpairinfo.dotprodTrack()); + data.fillMulti("dotprod2", trkpairinfo.dotprodSeed()); + + //pca distance form PV on both tracks + data.fillMulti("pca_dist1", trkpairinfo.pcaSeed_dist()); + data.fillMulti("pca_dist2", trkpairinfo.pcaTrack_dist()); + + //track track or dir dir dotprod + data.fillMulti("dotprod12_2D", trkpairinfo.dotprodTrackSeed2D()); + data.fillMulti("dotprod12_2DV", trkpairinfo.dotprodTrackSeed2DV()); + data.fillMulti("dotprod12_3D", trkpairinfo.dotprodTrackSeed3D()); + data.fillMulti("dotprod12_3DV", trkpairinfo.dotprodTrackSeed3DV()); + + //jet pca relative + data.fillMulti("pca_jetAxis_dist", trkpairinfo.pca_jetAxis_dist()); + data.fillMulti("pca_jetAxis_dotprod", trkpairinfo.pca_jetAxis_dotprod()); + data.fillMulti("pca_jetAxis_dEta", trkpairinfo.pca_jetAxis_dEta()); + data.fillMulti("pca_jetAxis_dPhi_", trkpairinfo.pca_jetAxis_dPhi()); + + } + + + } + + + } + + + return true; +} + +} /* namespace deepntuples */