Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:24:33 

0001 // -*- C++ -*-
0002 //
0003 // Package:    ​RecoBTag/​SecondaryVertex
0004 // Class:      BoostedDoubleSVProducer
0005 //
0006 /**\class BoostedDoubleSVProducer BoostedDoubleSVProducer.cc ​RecoBTag/​SecondaryVertex/plugins/BoostedDoubleSVProducer.cc
0007   *
0008   * Description: EDProducer that produces collection of BoostedDoubleSVTagInfos
0009   *
0010   * Implementation:
0011   *    A collection of SecondaryVertexTagInfos is taken as input and a collection of BoostedDoubleSVTagInfos
0012   *    is produced as output.
0013   */
0014 //
0015 // Original Author:  Dinko Ferencek
0016 //         Created:  Thu, 06 Oct 2016 14:02:30 GMT
0017 //
0018 //
0019 
0020 // system include files
0021 #include <memory>
0022 
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/StreamID.h"
0032 #include "FWCore/Utilities/interface/isFinite.h"
0033 #include "FWCore/Utilities/interface/ESGetToken.h"
0034 
0035 #include "DataFormats/BTauReco/interface/CandIPTagInfo.h"
0036 #include "DataFormats/BTauReco/interface/CandSecondaryVertexTagInfo.h"
0037 #include "DataFormats/BTauReco/interface/BoostedDoubleSVTagInfo.h"
0038 #include "DataFormats/BTauReco/interface/TaggingVariable.h"
0039 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0040 #include "DataFormats/PatCandidates/interface/PackedCandidate.h"
0041 #include "DataFormats/VertexReco/interface/Vertex.h"
0042 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0043 #include "DataFormats/Candidate/interface/VertexCompositePtrCandidate.h"
0044 
0045 #include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h"
0046 #include "RecoBTag/SecondaryVertex/interface/V0Filter.h"
0047 #include "RecoBTag/SecondaryVertex/interface/TrackSelector.h"
0048 #include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h"
0049 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0050 #include "TrackingTools/IPTools/interface/IPTools.h"
0051 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0052 
0053 #include "fastjet/PseudoJet.hh"
0054 #include "fastjet/contrib/Njettiness.hh"
0055 
0056 #include <map>
0057 
0058 //
0059 // class declaration
0060 //
0061 
0062 class BoostedDoubleSVProducer : public edm::stream::EDProducer<> {
0063 public:
0064   explicit BoostedDoubleSVProducer(const edm::ParameterSet&);
0065   ~BoostedDoubleSVProducer() override;
0066 
0067   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0068 
0069 private:
0070   void beginStream(edm::StreamID) override;
0071   void produce(edm::Event&, const edm::EventSetup&) override;
0072   void endStream() override;
0073 
0074   void calcNsubjettiness(const reco::JetBaseRef& jet,
0075                          float& tau1,
0076                          float& tau2,
0077                          std::vector<fastjet::PseudoJet>& currentAxes) const;
0078   void setTracksPVBase(const reco::TrackRef& trackRef, const reco::VertexRef& vertexRef, float& PVweight) const;
0079   void setTracksPV(const reco::CandidatePtr& trackRef, const reco::VertexRef& vertexRef, float& PVweight) const;
0080   void etaRelToTauAxis(const reco::VertexCompositePtrCandidate& vertex,
0081                        const fastjet::PseudoJet& tauAxis,
0082                        std::vector<float>& tau_trackEtaRel) const;
0083 
0084   // ----------member data ---------------------------
0085   const edm::EDGetTokenT<std::vector<reco::CandSecondaryVertexTagInfo>> svTagInfos_;
0086 
0087   const double beta_;
0088   const double R0_;
0089 
0090   const double maxSVDeltaRToJet_;
0091   const double maxDistToAxis_;
0092   const double maxDecayLen_;
0093   reco::V0Filter trackPairV0Filter;
0094   reco::TrackSelector trackSelector;
0095 
0096   edm::EDGetTokenT<edm::ValueMap<float>> weightsToken_;
0097   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> trackBuilderToken_;
0098   edm::Handle<edm::ValueMap<float>> weightsHandle_;
0099 
0100   // static variables
0101   static constexpr float dummyZ_ratio = -3.0f;
0102   static constexpr float dummyTrackSip3dSig = -50.0f;
0103   static constexpr float dummyTrackSip2dSigAbove = -19.0f;
0104   static constexpr float dummyTrackEtaRel = -1.0f;
0105   static constexpr float dummyVertexMass = -1.0f;
0106   static constexpr float dummyVertexEnergyRatio = -1.0f;
0107   static constexpr float dummyVertexDeltaR = -1.0f;
0108   static constexpr float dummyFlightDistance2dSig = -1.0f;
0109 
0110   static constexpr float charmThreshold = 1.5f;
0111   static constexpr float bottomThreshold = 5.2f;
0112 };
0113 
0114 //
0115 // constants, enums and typedefs
0116 //
0117 
0118 //
0119 // static data member definitions
0120 //
0121 
0122 //
0123 // constructors and destructor
0124 //
0125 BoostedDoubleSVProducer::BoostedDoubleSVProducer(const edm::ParameterSet& iConfig)
0126     : svTagInfos_(
0127           consumes<std::vector<reco::CandSecondaryVertexTagInfo>>(iConfig.getParameter<edm::InputTag>("svTagInfos"))),
0128       beta_(iConfig.getParameter<double>("beta")),
0129       R0_(iConfig.getParameter<double>("R0")),
0130       maxSVDeltaRToJet_(iConfig.getParameter<double>("maxSVDeltaRToJet")),
0131       maxDistToAxis_(iConfig.getParameter<edm::ParameterSet>("trackSelection").getParameter<double>("maxDistToAxis")),
0132       maxDecayLen_(iConfig.getParameter<edm::ParameterSet>("trackSelection").getParameter<double>("maxDecayLen")),
0133       trackPairV0Filter(iConfig.getParameter<edm::ParameterSet>("trackPairV0Filter")),
0134       trackSelector(iConfig.getParameter<edm::ParameterSet>("trackSelection")) {
0135   edm::InputTag srcWeights = iConfig.getParameter<edm::InputTag>("weights");
0136   trackBuilderToken_ =
0137       esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0138   if (!srcWeights.label().empty())
0139     weightsToken_ = consumes<edm::ValueMap<float>>(srcWeights);
0140   produces<std::vector<reco::BoostedDoubleSVTagInfo>>();
0141 }
0142 
0143 BoostedDoubleSVProducer::~BoostedDoubleSVProducer() {
0144   // do anything here that needs to be done at destruction time
0145   // (e.g. close files, deallocate resources etc.)
0146 }
0147 
0148 //
0149 // member functions
0150 //
0151 
0152 // ------------ method called to produce the data  ------------
0153 void BoostedDoubleSVProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0154   // get the track builder
0155   edm::ESHandle<TransientTrackBuilder> trackBuilder = iSetup.getHandle(trackBuilderToken_);
0156 
0157   // get input secondary vertex TagInfos
0158   edm::Handle<std::vector<reco::CandSecondaryVertexTagInfo>> svTagInfos;
0159   iEvent.getByToken(svTagInfos_, svTagInfos);
0160 
0161   if (!weightsToken_.isUninitialized())
0162     iEvent.getByToken(weightsToken_, weightsHandle_);
0163 
0164   // create the output collection
0165   auto tagInfos = std::make_unique<std::vector<reco::BoostedDoubleSVTagInfo>>();
0166 
0167   // loop over TagInfos
0168   for (std::vector<reco::CandSecondaryVertexTagInfo>::const_iterator iterTI = svTagInfos->begin();
0169        iterTI != svTagInfos->end();
0170        ++iterTI) {
0171     // get TagInfos
0172     const reco::CandIPTagInfo& ipTagInfo = *(iterTI->trackIPTagInfoRef().get());
0173     const reco::CandSecondaryVertexTagInfo& svTagInfo = *(iterTI);
0174 
0175     // default variable values
0176     float z_ratio = dummyZ_ratio;
0177     float trackSip3dSig_3 = dummyTrackSip3dSig, trackSip3dSig_2 = dummyTrackSip3dSig,
0178           trackSip3dSig_1 = dummyTrackSip3dSig, trackSip3dSig_0 = dummyTrackSip3dSig;
0179     float tau2_trackSip3dSig_0 = dummyTrackSip3dSig, tau1_trackSip3dSig_0 = dummyTrackSip3dSig,
0180           tau2_trackSip3dSig_1 = dummyTrackSip3dSig, tau1_trackSip3dSig_1 = dummyTrackSip3dSig;
0181     float trackSip2dSigAboveCharm_0 = dummyTrackSip2dSigAbove, trackSip2dSigAboveBottom_0 = dummyTrackSip2dSigAbove,
0182           trackSip2dSigAboveBottom_1 = dummyTrackSip2dSigAbove;
0183     float tau1_trackEtaRel_0 = dummyTrackEtaRel, tau1_trackEtaRel_1 = dummyTrackEtaRel,
0184           tau1_trackEtaRel_2 = dummyTrackEtaRel;
0185     float tau2_trackEtaRel_0 = dummyTrackEtaRel, tau2_trackEtaRel_1 = dummyTrackEtaRel,
0186           tau2_trackEtaRel_2 = dummyTrackEtaRel;
0187     float tau1_vertexMass = dummyVertexMass, tau1_vertexEnergyRatio = dummyVertexEnergyRatio,
0188           tau1_vertexDeltaR = dummyVertexDeltaR, tau1_flightDistance2dSig = dummyFlightDistance2dSig;
0189     float tau2_vertexMass = dummyVertexMass, tau2_vertexEnergyRatio = dummyVertexEnergyRatio,
0190           tau2_vertexDeltaR = dummyVertexDeltaR, tau2_flightDistance2dSig = dummyFlightDistance2dSig;
0191     float jetNTracks = 0, nSV = 0, tau1_nSecondaryVertices = 0, tau2_nSecondaryVertices = 0;
0192 
0193     // get the jet reference
0194     const reco::JetBaseRef jet = svTagInfo.jet();
0195 
0196     std::vector<fastjet::PseudoJet> currentAxes;
0197     float tau2, tau1;
0198     // calculate N-subjettiness
0199     calcNsubjettiness(jet, tau1, tau2, currentAxes);
0200 
0201     const reco::VertexRef& vertexRef = ipTagInfo.primaryVertex();
0202     GlobalPoint pv(0., 0., 0.);
0203     if (ipTagInfo.primaryVertex().isNonnull())
0204       pv = GlobalPoint(vertexRef->x(), vertexRef->y(), vertexRef->z());
0205 
0206     const std::vector<reco::CandidatePtr>& selectedTracks = ipTagInfo.selectedTracks();
0207     const std::vector<reco::btag::TrackIPData>& ipData = ipTagInfo.impactParameterData();
0208     size_t trackSize = selectedTracks.size();
0209 
0210     reco::TrackKinematics allKinematics;
0211     std::vector<float> IP3Ds, IP3Ds_1, IP3Ds_2;
0212     int contTrk = 0;
0213 
0214     // loop over tracks associated to the jet
0215     for (size_t itt = 0; itt < trackSize; ++itt) {
0216       const reco::CandidatePtr trackRef = selectedTracks[itt];
0217 
0218       float track_PVweight = 0.;
0219       setTracksPV(trackRef, vertexRef, track_PVweight);
0220       if (track_PVweight > 0.5)
0221         allKinematics.add(trackRef);
0222 
0223       const reco::btag::TrackIPData& data = ipData[itt];
0224       bool isSelected = false;
0225       if (trackSelector(trackRef, data, *jet, pv))
0226         isSelected = true;
0227 
0228       // check if the track is from V0
0229       bool isfromV0 = false, isfromV0Tight = false;
0230       std::vector<reco::CandidatePtr> trackPairV0Test(2);
0231 
0232       trackPairV0Test[0] = trackRef;
0233 
0234       for (size_t jtt = 0; jtt < trackSize; ++jtt) {
0235         if (itt == jtt)
0236           continue;
0237 
0238         const reco::btag::TrackIPData& pairTrackData = ipData[jtt];
0239         const reco::CandidatePtr pairTrackRef = selectedTracks[jtt];
0240 
0241         trackPairV0Test[1] = pairTrackRef;
0242 
0243         if (!trackPairV0Filter(trackPairV0Test)) {
0244           isfromV0 = true;
0245 
0246           if (trackSelector(pairTrackRef, pairTrackData, *jet, pv))
0247             isfromV0Tight = true;
0248         }
0249 
0250         if (isfromV0 && isfromV0Tight)
0251           break;
0252       }
0253 
0254       if (isSelected && !isfromV0Tight)
0255         jetNTracks += 1.;
0256 
0257       reco::TransientTrack transientTrack = trackBuilder->build(trackRef);
0258       GlobalVector direction(jet->px(), jet->py(), jet->pz());
0259 
0260       int index = 0;
0261       if (currentAxes.size() > 1 &&
0262           reco::deltaR2(trackRef->momentum(), currentAxes[1]) < reco::deltaR2(trackRef->momentum(), currentAxes[0]))
0263         index = 1;
0264       direction = GlobalVector(currentAxes[index].px(), currentAxes[index].py(), currentAxes[index].pz());
0265 
0266       // decay distance and track distance wrt to the closest tau axis
0267       float decayLengthTau = -1;
0268       float distTauAxis = -1;
0269 
0270       TrajectoryStateOnSurface closest = IPTools::closestApproachToJet(
0271           transientTrack.impactPointState(), *vertexRef, direction, transientTrack.field());
0272       if (closest.isValid())
0273         decayLengthTau = (closest.globalPosition() - RecoVertex::convertPos(vertexRef->position())).mag();
0274 
0275       distTauAxis = std::abs(IPTools::jetTrackDistance(transientTrack, direction, *vertexRef).second.value());
0276 
0277       float IP3Dsig = ipTagInfo.impactParameterData()[itt].ip3d.significance();
0278 
0279       if (!isfromV0 && decayLengthTau < maxDecayLen_ && distTauAxis < maxDistToAxis_) {
0280         IP3Ds.push_back(IP3Dsig < -50. ? -50. : IP3Dsig);
0281         ++contTrk;
0282         if (currentAxes.size() > 1) {
0283           if (reco::deltaR2(trackRef->momentum(), currentAxes[0]) < reco::deltaR2(trackRef->momentum(), currentAxes[1]))
0284             IP3Ds_1.push_back(IP3Dsig < -50. ? -50. : IP3Dsig);
0285           else
0286             IP3Ds_2.push_back(IP3Dsig < -50. ? -50. : IP3Dsig);
0287         } else
0288           IP3Ds_1.push_back(IP3Dsig < -50. ? -50. : IP3Dsig);
0289       }
0290     }
0291 
0292     std::vector<size_t> indices = ipTagInfo.sortedIndexes(reco::btag::IP2DSig);
0293     bool charmThreshSet = false;
0294 
0295     reco::TrackKinematics kin;
0296     for (size_t i = 0; i < indices.size(); ++i) {
0297       size_t idx = indices[i];
0298       const reco::btag::TrackIPData& data = ipData[idx];
0299       const reco::CandidatePtr trackRef = selectedTracks[idx];
0300 
0301       kin.add(trackRef);
0302 
0303       if (kin.vectorSum().M() > charmThreshold  // charm cut
0304           && !charmThreshSet) {
0305         trackSip2dSigAboveCharm_0 = data.ip2d.significance();
0306 
0307         charmThreshSet = true;
0308       }
0309 
0310       if (kin.vectorSum().M() > bottomThreshold)  // bottom cut
0311       {
0312         trackSip2dSigAboveBottom_0 = data.ip2d.significance();
0313         if ((i + 1) < indices.size())
0314           trackSip2dSigAboveBottom_1 = (ipData[indices[i + 1]]).ip2d.significance();
0315 
0316         break;
0317       }
0318     }
0319 
0320     float dummyTrack = -50.;
0321 
0322     std::sort(IP3Ds.begin(), IP3Ds.end(), std::greater<float>());
0323     std::sort(IP3Ds_1.begin(), IP3Ds_1.end(), std::greater<float>());
0324     std::sort(IP3Ds_2.begin(), IP3Ds_2.end(), std::greater<float>());
0325     int num_1 = IP3Ds_1.size();
0326     int num_2 = IP3Ds_2.size();
0327 
0328     switch (contTrk) {
0329       case 0:
0330 
0331         trackSip3dSig_0 = dummyTrack;
0332         trackSip3dSig_1 = dummyTrack;
0333         trackSip3dSig_2 = dummyTrack;
0334         trackSip3dSig_3 = dummyTrack;
0335 
0336         break;
0337 
0338       case 1:
0339 
0340         trackSip3dSig_0 = IP3Ds.at(0);
0341         trackSip3dSig_1 = dummyTrack;
0342         trackSip3dSig_2 = dummyTrack;
0343         trackSip3dSig_3 = dummyTrack;
0344 
0345         break;
0346 
0347       case 2:
0348 
0349         trackSip3dSig_0 = IP3Ds.at(0);
0350         trackSip3dSig_1 = IP3Ds.at(1);
0351         trackSip3dSig_2 = dummyTrack;
0352         trackSip3dSig_3 = dummyTrack;
0353 
0354         break;
0355 
0356       case 3:
0357 
0358         trackSip3dSig_0 = IP3Ds.at(0);
0359         trackSip3dSig_1 = IP3Ds.at(1);
0360         trackSip3dSig_2 = IP3Ds.at(2);
0361         trackSip3dSig_3 = dummyTrack;
0362 
0363         break;
0364 
0365       default:
0366 
0367         trackSip3dSig_0 = IP3Ds.at(0);
0368         trackSip3dSig_1 = IP3Ds.at(1);
0369         trackSip3dSig_2 = IP3Ds.at(2);
0370         trackSip3dSig_3 = IP3Ds.at(3);
0371     }
0372 
0373     switch (num_1) {
0374       case 0:
0375 
0376         tau1_trackSip3dSig_0 = dummyTrack;
0377         tau1_trackSip3dSig_1 = dummyTrack;
0378 
0379         break;
0380 
0381       case 1:
0382 
0383         tau1_trackSip3dSig_0 = IP3Ds_1.at(0);
0384         tau1_trackSip3dSig_1 = dummyTrack;
0385 
0386         break;
0387 
0388       default:
0389 
0390         tau1_trackSip3dSig_0 = IP3Ds_1.at(0);
0391         tau1_trackSip3dSig_1 = IP3Ds_1.at(1);
0392     }
0393 
0394     switch (num_2) {
0395       case 0:
0396 
0397         tau2_trackSip3dSig_0 = dummyTrack;
0398         tau2_trackSip3dSig_1 = dummyTrack;
0399 
0400         break;
0401 
0402       case 1:
0403         tau2_trackSip3dSig_0 = IP3Ds_2.at(0);
0404         tau2_trackSip3dSig_1 = dummyTrack;
0405 
0406         break;
0407 
0408       default:
0409 
0410         tau2_trackSip3dSig_0 = IP3Ds_2.at(0);
0411         tau2_trackSip3dSig_1 = IP3Ds_2.at(1);
0412     }
0413 
0414     math::XYZVector jetDir = jet->momentum().Unit();
0415     reco::TrackKinematics tau1Kinematics;
0416     reco::TrackKinematics tau2Kinematics;
0417     std::vector<float> tau1_trackEtaRels, tau2_trackEtaRels;
0418 
0419     std::map<double, size_t> VTXmap;
0420     for (size_t vtx = 0; vtx < svTagInfo.nVertices(); ++vtx) {
0421       const reco::VertexCompositePtrCandidate& vertex = svTagInfo.secondaryVertex(vtx);
0422       // get the vertex kinematics
0423       reco::TrackKinematics vertexKinematic(vertex);
0424 
0425       if (currentAxes.size() > 1) {
0426         if (reco::deltaR2(svTagInfo.flightDirection(vtx), currentAxes[1]) <
0427             reco::deltaR2(svTagInfo.flightDirection(vtx), currentAxes[0])) {
0428           tau2Kinematics = tau2Kinematics + vertexKinematic;
0429           if (tau2_flightDistance2dSig < 0) {
0430             tau2_flightDistance2dSig = svTagInfo.flightDistance(vtx, true).significance();
0431             tau2_vertexDeltaR = reco::deltaR(svTagInfo.flightDirection(vtx), currentAxes[1]);
0432           }
0433           etaRelToTauAxis(vertex, currentAxes[1], tau2_trackEtaRels);
0434           tau2_nSecondaryVertices += 1.;
0435         } else {
0436           tau1Kinematics = tau1Kinematics + vertexKinematic;
0437           if (tau1_flightDistance2dSig < 0) {
0438             tau1_flightDistance2dSig = svTagInfo.flightDistance(vtx, true).significance();
0439             tau1_vertexDeltaR = reco::deltaR(svTagInfo.flightDirection(vtx), currentAxes[0]);
0440           }
0441           etaRelToTauAxis(vertex, currentAxes[0], tau1_trackEtaRels);
0442           tau1_nSecondaryVertices += 1.;
0443         }
0444 
0445       } else if (!currentAxes.empty()) {
0446         tau1Kinematics = tau1Kinematics + vertexKinematic;
0447         if (tau1_flightDistance2dSig < 0) {
0448           tau1_flightDistance2dSig = svTagInfo.flightDistance(vtx, true).significance();
0449           tau1_vertexDeltaR = reco::deltaR(svTagInfo.flightDirection(vtx), currentAxes[0]);
0450         }
0451         etaRelToTauAxis(vertex, currentAxes[0], tau1_trackEtaRels);
0452         tau1_nSecondaryVertices += 1.;
0453       }
0454 
0455       const GlobalVector& flightDir = svTagInfo.flightDirection(vtx);
0456       if (reco::deltaR2(flightDir, jetDir) < (maxSVDeltaRToJet_ * maxSVDeltaRToJet_))
0457         VTXmap[svTagInfo.flightDistance(vtx).error()] = vtx;
0458     }
0459     nSV = VTXmap.size();
0460 
0461     math::XYZTLorentzVector allSum = allKinematics.weightedVectorSum();
0462     if (tau1_nSecondaryVertices > 0.) {
0463       const math::XYZTLorentzVector& tau1_vertexSum = tau1Kinematics.weightedVectorSum();
0464       if (allSum.E() > 0.)
0465         tau1_vertexEnergyRatio = tau1_vertexSum.E() / allSum.E();
0466       if (tau1_vertexEnergyRatio > 50.)
0467         tau1_vertexEnergyRatio = 50.;
0468 
0469       tau1_vertexMass = tau1_vertexSum.M();
0470     }
0471 
0472     if (tau2_nSecondaryVertices > 0.) {
0473       const math::XYZTLorentzVector& tau2_vertexSum = tau2Kinematics.weightedVectorSum();
0474       if (allSum.E() > 0.)
0475         tau2_vertexEnergyRatio = tau2_vertexSum.E() / allSum.E();
0476       if (tau2_vertexEnergyRatio > 50.)
0477         tau2_vertexEnergyRatio = 50.;
0478 
0479       tau2_vertexMass = tau2_vertexSum.M();
0480     }
0481 
0482     float dummyEtaRel = -1.;
0483 
0484     std::sort(tau1_trackEtaRels.begin(), tau1_trackEtaRels.end());
0485     std::sort(tau2_trackEtaRels.begin(), tau2_trackEtaRels.end());
0486 
0487     switch (tau2_trackEtaRels.size()) {
0488       case 0:
0489 
0490         tau2_trackEtaRel_0 = dummyEtaRel;
0491         tau2_trackEtaRel_1 = dummyEtaRel;
0492         tau2_trackEtaRel_2 = dummyEtaRel;
0493 
0494         break;
0495 
0496       case 1:
0497 
0498         tau2_trackEtaRel_0 = tau2_trackEtaRels.at(0);
0499         tau2_trackEtaRel_1 = dummyEtaRel;
0500         tau2_trackEtaRel_2 = dummyEtaRel;
0501 
0502         break;
0503 
0504       case 2:
0505 
0506         tau2_trackEtaRel_0 = tau2_trackEtaRels.at(0);
0507         tau2_trackEtaRel_1 = tau2_trackEtaRels.at(1);
0508         tau2_trackEtaRel_2 = dummyEtaRel;
0509 
0510         break;
0511 
0512       default:
0513 
0514         tau2_trackEtaRel_0 = tau2_trackEtaRels.at(0);
0515         tau2_trackEtaRel_1 = tau2_trackEtaRels.at(1);
0516         tau2_trackEtaRel_2 = tau2_trackEtaRels.at(2);
0517     }
0518 
0519     switch (tau1_trackEtaRels.size()) {
0520       case 0:
0521 
0522         tau1_trackEtaRel_0 = dummyEtaRel;
0523         tau1_trackEtaRel_1 = dummyEtaRel;
0524         tau1_trackEtaRel_2 = dummyEtaRel;
0525 
0526         break;
0527 
0528       case 1:
0529 
0530         tau1_trackEtaRel_0 = tau1_trackEtaRels.at(0);
0531         tau1_trackEtaRel_1 = dummyEtaRel;
0532         tau1_trackEtaRel_2 = dummyEtaRel;
0533 
0534         break;
0535 
0536       case 2:
0537 
0538         tau1_trackEtaRel_0 = tau1_trackEtaRels.at(0);
0539         tau1_trackEtaRel_1 = tau1_trackEtaRels.at(1);
0540         tau1_trackEtaRel_2 = dummyEtaRel;
0541 
0542         break;
0543 
0544       default:
0545 
0546         tau1_trackEtaRel_0 = tau1_trackEtaRels.at(0);
0547         tau1_trackEtaRel_1 = tau1_trackEtaRels.at(1);
0548         tau1_trackEtaRel_2 = tau1_trackEtaRels.at(2);
0549     }
0550 
0551     int cont = 0;
0552     GlobalVector flightDir_0, flightDir_1;
0553     reco::Candidate::LorentzVector SV_p4_0, SV_p4_1;
0554     double vtxMass = 0.;
0555 
0556     for (std::map<double, size_t>::iterator iVtx = VTXmap.begin(); iVtx != VTXmap.end(); ++iVtx) {
0557       ++cont;
0558       const reco::VertexCompositePtrCandidate& vertex = svTagInfo.secondaryVertex(iVtx->second);
0559       if (cont == 1) {
0560         flightDir_0 = svTagInfo.flightDirection(iVtx->second);
0561         SV_p4_0 = vertex.p4();
0562         vtxMass = SV_p4_0.mass();
0563 
0564         if (vtxMass > 0.)
0565           z_ratio = reco::deltaR(currentAxes[1], currentAxes[0]) * SV_p4_0.pt() / vtxMass;
0566       }
0567       if (cont == 2) {
0568         flightDir_1 = svTagInfo.flightDirection(iVtx->second);
0569         SV_p4_1 = vertex.p4();
0570         vtxMass = (SV_p4_1 + SV_p4_0).mass();
0571 
0572         if (vtxMass > 0.)
0573           z_ratio = reco::deltaR(flightDir_0, flightDir_1) * SV_p4_1.pt() / vtxMass;
0574 
0575         break;
0576       }
0577     }
0578 
0579     // when only one tau axis has SVs assigned, they are all assigned to the 1st tau axis
0580     // in the special case below need to swap values
0581     if ((tau1_vertexMass < 0 && tau2_vertexMass > 0)) {
0582       float temp = tau1_trackEtaRel_0;
0583       tau1_trackEtaRel_0 = tau2_trackEtaRel_0;
0584       tau2_trackEtaRel_0 = temp;
0585 
0586       temp = tau1_trackEtaRel_1;
0587       tau1_trackEtaRel_1 = tau2_trackEtaRel_1;
0588       tau2_trackEtaRel_1 = temp;
0589 
0590       temp = tau1_trackEtaRel_2;
0591       tau1_trackEtaRel_2 = tau2_trackEtaRel_2;
0592       tau2_trackEtaRel_2 = temp;
0593 
0594       temp = tau1_flightDistance2dSig;
0595       tau1_flightDistance2dSig = tau2_flightDistance2dSig;
0596       tau2_flightDistance2dSig = temp;
0597 
0598       tau1_vertexDeltaR = tau2_vertexDeltaR;
0599 
0600       temp = tau1_vertexEnergyRatio;
0601       tau1_vertexEnergyRatio = tau2_vertexEnergyRatio;
0602       tau2_vertexEnergyRatio = temp;
0603 
0604       temp = tau1_vertexMass;
0605       tau1_vertexMass = tau2_vertexMass;
0606       tau2_vertexMass = temp;
0607     }
0608 
0609     reco::TaggingVariableList vars;
0610 
0611     vars.insert(reco::btau::jetNTracks, jetNTracks, true);
0612     vars.insert(reco::btau::jetNSecondaryVertices, nSV, true);
0613     vars.insert(reco::btau::trackSip3dSig_0, trackSip3dSig_0, true);
0614     vars.insert(reco::btau::trackSip3dSig_1, trackSip3dSig_1, true);
0615     vars.insert(reco::btau::trackSip3dSig_2, trackSip3dSig_2, true);
0616     vars.insert(reco::btau::trackSip3dSig_3, trackSip3dSig_3, true);
0617     vars.insert(reco::btau::tau1_trackSip3dSig_0, tau1_trackSip3dSig_0, true);
0618     vars.insert(reco::btau::tau1_trackSip3dSig_1, tau1_trackSip3dSig_1, true);
0619     vars.insert(reco::btau::tau2_trackSip3dSig_0, tau2_trackSip3dSig_0, true);
0620     vars.insert(reco::btau::tau2_trackSip3dSig_1, tau2_trackSip3dSig_1, true);
0621     vars.insert(reco::btau::trackSip2dSigAboveCharm, trackSip2dSigAboveCharm_0, true);
0622     vars.insert(reco::btau::trackSip2dSigAboveBottom_0, trackSip2dSigAboveBottom_0, true);
0623     vars.insert(reco::btau::trackSip2dSigAboveBottom_1, trackSip2dSigAboveBottom_1, true);
0624     vars.insert(reco::btau::tau1_trackEtaRel_0, tau1_trackEtaRel_0, true);
0625     vars.insert(reco::btau::tau1_trackEtaRel_1, tau1_trackEtaRel_1, true);
0626     vars.insert(reco::btau::tau1_trackEtaRel_2, tau1_trackEtaRel_2, true);
0627     vars.insert(reco::btau::tau2_trackEtaRel_0, tau2_trackEtaRel_0, true);
0628     vars.insert(reco::btau::tau2_trackEtaRel_1, tau2_trackEtaRel_1, true);
0629     vars.insert(reco::btau::tau2_trackEtaRel_2, tau2_trackEtaRel_2, true);
0630     vars.insert(reco::btau::tau1_vertexMass, tau1_vertexMass, true);
0631     vars.insert(reco::btau::tau1_vertexEnergyRatio, tau1_vertexEnergyRatio, true);
0632     vars.insert(reco::btau::tau1_flightDistance2dSig, tau1_flightDistance2dSig, true);
0633     vars.insert(reco::btau::tau1_vertexDeltaR, tau1_vertexDeltaR, true);
0634     vars.insert(reco::btau::tau2_vertexMass, tau2_vertexMass, true);
0635     vars.insert(reco::btau::tau2_vertexEnergyRatio, tau2_vertexEnergyRatio, true);
0636     vars.insert(reco::btau::tau2_flightDistance2dSig, tau2_flightDistance2dSig, true);
0637     vars.insert(reco::btau::z_ratio, z_ratio, true);
0638 
0639     vars.finalize();
0640 
0641     tagInfos->push_back(reco::BoostedDoubleSVTagInfo(
0642         vars, edm::Ref<std::vector<reco::CandSecondaryVertexTagInfo>>(svTagInfos, iterTI - svTagInfos->begin())));
0643   }
0644 
0645   // put the output in the event
0646   iEvent.put(std::move(tagInfos));
0647 }
0648 
0649 void BoostedDoubleSVProducer::calcNsubjettiness(const reco::JetBaseRef& jet,
0650                                                 float& tau1,
0651                                                 float& tau2,
0652                                                 std::vector<fastjet::PseudoJet>& currentAxes) const {
0653   std::vector<fastjet::PseudoJet> fjParticles;
0654 
0655   // loop over jet constituents and push them in the vector of FastJet constituents
0656   for (const reco::CandidatePtr& daughter : jet->daughterPtrVector()) {
0657     if (daughter.isNonnull() && daughter.isAvailable()) {
0658       const reco::Jet* subjet = dynamic_cast<const reco::Jet*>(daughter.get());
0659       // if the daughter is actually a subjet
0660       if (subjet && daughter->numberOfDaughters() > 1) {
0661         // loop over subjet constituents and push them in the vector of FastJet constituents
0662         for (size_t i = 0; i < daughter->numberOfDaughters(); ++i) {
0663           const reco::CandidatePtr& constit = subjet->daughterPtr(i);
0664 
0665           if (constit.isNonnull() && constit->pt() > std::numeric_limits<double>::epsilon()) {
0666             // Check if any values were nan or inf
0667             float valcheck = constit->px() + constit->py() + constit->pz() + constit->energy();
0668             if (edm::isNotFinite(valcheck)) {
0669               edm::LogWarning("FaultyJetConstituent")
0670                   << "Jet constituent required for N-subjettiness computation contains Nan/Inf values!";
0671               continue;
0672             }
0673             if (subjet->isWeighted()) {
0674               float w = 0.0;
0675               if (!weightsToken_.isUninitialized())
0676                 w = (*weightsHandle_)[constit];
0677               else {
0678                 throw cms::Exception("MissingConstituentWeight")
0679                     << "BoostedDoubleSVProducer: No weights (e.g. PUPPI) given for weighted jet collection"
0680                     << std::endl;
0681               }
0682               if (w > 0) {
0683                 fjParticles.push_back(
0684                     fastjet::PseudoJet(constit->px() * w, constit->py() * w, constit->pz() * w, constit->energy() * w));
0685               }
0686             } else {
0687               fjParticles.push_back(fastjet::PseudoJet(constit->px(), constit->py(), constit->pz(), constit->energy()));
0688             }
0689           } else
0690             edm::LogWarning("MissingJetConstituent")
0691                 << "Jet constituent required for N-subjettiness computation is missing!";
0692         }
0693       } else {
0694         // Check if any values were nan or inf
0695         float valcheck = daughter->px() + daughter->py() + daughter->pz() + daughter->energy();
0696         if (edm::isNotFinite(valcheck)) {
0697           edm::LogWarning("FaultyJetConstituent")
0698               << "Jet constituent required for N-subjettiness computation contains Nan/Inf values!";
0699           continue;
0700         }
0701         if (jet->isWeighted()) {
0702           float w = 0.0;
0703           if (!weightsToken_.isUninitialized())
0704             w = (*weightsHandle_)[daughter];
0705           else {
0706             throw cms::Exception("MissingConstituentWeight")
0707                 << "BoostedDoubleSVProducer: No weights (e.g. PUPPI) given for weighted jet collection" << std::endl;
0708           }
0709           if (w > 0 && daughter->pt() > std::numeric_limits<double>::epsilon()) {
0710             fjParticles.push_back(
0711                 fastjet::PseudoJet(daughter->px() * w, daughter->py() * w, daughter->pz() * w, daughter->energy() * w));
0712           }
0713         } else {
0714           fjParticles.push_back(fastjet::PseudoJet(daughter->px(), daughter->py(), daughter->pz(), daughter->energy()));
0715         }
0716       }
0717     } else
0718       edm::LogWarning("MissingJetConstituent") << "Jet constituent required for N-subjettiness computation is missing!";
0719   }
0720 
0721   // N-subjettiness calculator
0722   fastjet::contrib::Njettiness njettiness(fastjet::contrib::OnePass_KT_Axes(),
0723                                           fastjet::contrib::NormalizedMeasure(beta_, R0_));
0724 
0725   // calculate N-subjettiness
0726   tau1 = njettiness.getTau(1, fjParticles);
0727   tau2 = njettiness.getTau(2, fjParticles);
0728   currentAxes = njettiness.currentAxes();
0729 }
0730 
0731 void BoostedDoubleSVProducer::setTracksPVBase(const reco::TrackRef& trackRef,
0732                                               const reco::VertexRef& vertexRef,
0733                                               float& PVweight) const {
0734   PVweight = 0.;
0735 
0736   const reco::TrackBaseRef trackBaseRef(trackRef);
0737 
0738   typedef reco::Vertex::trackRef_iterator IT;
0739 
0740   const reco::Vertex& vtx = *(vertexRef);
0741   // loop over tracks in vertices
0742   for (IT it = vtx.tracks_begin(); it != vtx.tracks_end(); ++it) {
0743     const reco::TrackBaseRef& baseRef = *it;
0744     // one of the tracks in the vertex is the same as the track considered in the function
0745     if (baseRef == trackBaseRef) {
0746       PVweight = vtx.trackWeight(baseRef);
0747       break;
0748     }
0749   }
0750 }
0751 
0752 void BoostedDoubleSVProducer::setTracksPV(const reco::CandidatePtr& trackRef,
0753                                           const reco::VertexRef& vertexRef,
0754                                           float& PVweight) const {
0755   PVweight = 0.;
0756 
0757   const pat::PackedCandidate* pcand = dynamic_cast<const pat::PackedCandidate*>(trackRef.get());
0758 
0759   if (pcand)  // MiniAOD case
0760   {
0761     if (pcand->fromPV() == pat::PackedCandidate::PVUsedInFit) {
0762       PVweight = 1.;
0763     }
0764   } else {
0765     const reco::PFCandidate* pfcand = dynamic_cast<const reco::PFCandidate*>(trackRef.get());
0766 
0767     setTracksPVBase(pfcand->trackRef(), vertexRef, PVweight);
0768   }
0769 }
0770 
0771 void BoostedDoubleSVProducer::etaRelToTauAxis(const reco::VertexCompositePtrCandidate& vertex,
0772                                               const fastjet::PseudoJet& tauAxis,
0773                                               std::vector<float>& tau_trackEtaRel) const {
0774   math::XYZVector direction(tauAxis.px(), tauAxis.py(), tauAxis.pz());
0775   const std::vector<reco::CandidatePtr>& tracks = vertex.daughterPtrVector();
0776 
0777   for (std::vector<reco::CandidatePtr>::const_iterator track = tracks.begin(); track != tracks.end(); ++track)
0778     tau_trackEtaRel.push_back(std::abs(reco::btau::etaRel(direction.Unit(), (*track)->momentum())));
0779 }
0780 
0781 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0782 void BoostedDoubleSVProducer::beginStream(edm::StreamID) {}
0783 
0784 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0785 void BoostedDoubleSVProducer::endStream() {}
0786 
0787 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0788 void BoostedDoubleSVProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0789   edm::ParameterSetDescription desc;
0790   desc.add<double>("beta", 1.0);
0791   desc.add<double>("R0", 0.8);
0792   desc.add<double>("maxSVDeltaRToJet", 0.7);
0793   {
0794     edm::ParameterSetDescription trackSelection;
0795     trackSelection.setAllowAnything();
0796     desc.add<edm::ParameterSetDescription>("trackSelection", trackSelection);
0797   }
0798   {
0799     edm::ParameterSetDescription trackPairV0Filter;
0800     trackPairV0Filter.add<double>("k0sMassWindow", 0.03);
0801     desc.add<edm::ParameterSetDescription>("trackPairV0Filter", trackPairV0Filter);
0802   }
0803   desc.add<edm::InputTag>("svTagInfos", edm::InputTag("pfInclusiveSecondaryVertexFinderAK8TagInfos"));
0804   desc.add<edm::InputTag>("weights", edm::InputTag(""));
0805   descriptions.addDefault(desc);
0806 }
0807 
0808 //define this as a plug-in
0809 DEFINE_FWK_MODULE(BoostedDoubleSVProducer);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!