File indexing completed on 2024-04-06 12:13:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include <memory>
0014 #include <iostream>
0015 #include <set>
0016
0017
0018 #include "Math/Vector4D.h"
0019 #include "Math/Vector4Dfwd.h"
0020
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022 #include "FWCore/Framework/interface/global/EDFilter.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0029
0030
0031
0032
0033 class LHEmttFilter : public edm::global::EDFilter<> {
0034 public:
0035 explicit LHEmttFilter(const edm::ParameterSet&);
0036 ~LHEmttFilter() override;
0037
0038 bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0039
0040 private:
0041
0042
0043 edm::EDGetTokenT<LHEEventProduct> src_;
0044 double ptMin_ = 0;
0045 double MinInvMass_ = -1;
0046 double MaxInvMass_ = -1;
0047 };
0048
0049 using namespace edm;
0050 using namespace std;
0051
0052 LHEmttFilter::LHEmttFilter(const edm::ParameterSet& iConfig)
0053 : ptMin_(iConfig.getParameter<double>("ptMin")),
0054 MinInvMass_(iConfig.getParameter<double>("MinInvMass")),
0055 MaxInvMass_(iConfig.getParameter<double>("MaxInvMass")) {
0056 src_ = consumes<LHEEventProduct>(iConfig.getParameter<edm::InputTag>("src"));
0057 }
0058
0059 LHEmttFilter::~LHEmttFilter() {
0060
0061
0062 }
0063
0064
0065 bool LHEmttFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0066 edm::Handle<LHEEventProduct> EvtHandle;
0067 iEvent.getByToken(src_, EvtHandle);
0068
0069 std::vector<lhef::HEPEUP::FiveVector> lheParticles = EvtHandle->hepeup().PUP;
0070 std::vector<ROOT::Math::PxPyPzMVector> cands;
0071 std::vector<int> pdgId_cands;
0072
0073 for (unsigned int i = 0; i < lheParticles.size(); i++) {
0074 if (EvtHandle->hepeup().ISTUP[i] != 2) {
0075 continue;
0076 }
0077 int pdgId = EvtHandle->hepeup().IDUP[i];
0078 if (std::abs(pdgId) != 6) {
0079 continue;
0080 }
0081 pdgId_cands.push_back(pdgId);
0082 cands.push_back(
0083 ROOT::Math::PxPyPzMVector(lheParticles[i][0], lheParticles[i][1], lheParticles[i][2], lheParticles[i][4]));
0084 }
0085
0086 if (cands.size() != 2) {
0087 edm::LogWarning("LHEmttFilter Error") << "Number of top quarks found in the event != 2" << endl;
0088 }
0089
0090 double ttmass_ = -1;
0091 if (!cands.empty() && (cands.size() == 2)) {
0092 if (pdgId_cands.at(0) + pdgId_cands.at(1) == 0) {
0093 if ((cands.at(0).pt() > ptMin_) && (cands.at(1).pt() > ptMin_)) {
0094 ROOT::Math::PxPyPzMVector tot = cands.at(0);
0095 for (unsigned icand = 1; icand < cands.size(); ++icand) {
0096 tot += cands.at(icand);
0097 }
0098 ttmass_ = tot.mass();
0099 }
0100 } else if (pdgId_cands.at(0) + pdgId_cands.at(1) != 0) {
0101 edm::LogWarning("LHEmttFilter Error") << "Found two t/tbar quarks instead of a ttbar pair" << endl;
0102 }
0103 }
0104
0105 if ((MinInvMass_ > -1 && ttmass_ < MinInvMass_) || (MaxInvMass_ > -1 && ttmass_ > MaxInvMass_)) {
0106 return false;
0107 } else {
0108 return true;
0109 }
0110 }
0111
0112
0113 DEFINE_FWK_MODULE(LHEmttFilter);