special/plugins/HLTPPSJetComparisonFilter.cc

0001 // Author: Mariana Araujo
0002 // Created: 2019-12-10
0003 /*
0004 Description:
0005 HLT filter module to select events according to matching of central (jets) and PPS (RP tracks) kinematics
0006
0007 Implementation:
0008 Matching can be done on the xi and/or mass+rapidity variables, using the do_xi and do_my booleans. If both are set to true, both matching conditions must be met
0009 */
0010
0011 // include files
0012 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/global/EDFilter.h"
0017
0018 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0019 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0020 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0021 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0022
0023 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0024 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0025
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027
0028 // class declaration
0029 //
0030 class HLTPPSJetComparisonFilter : public edm::global::EDFilter<> {
0031 public:
0032   explicit HLTPPSJetComparisonFilter(const edm::ParameterSet &);
0033   ~HLTPPSJetComparisonFilter() override;
0034
0035   static void fillDescriptions(edm::ConfigurationDescriptions &);
0036   bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0037
0038 private:
0039   // ----------member data ---------------------------
0040   edm::ESGetToken<LHCInfo, LHCInfoRcd> const lhcInfoRcdToken_;
0041
0042   edm::ParameterSet param_;
0043
0044   edm::InputTag jetInputTag_;  // Input tag identifying the jet track
0045   edm::EDGetTokenT<reco::PFJetCollection> jet_token_;
0046
0047   edm::InputTag forwardProtonInputTag_;  // Input tag identifying the forward proton collection
0048   edm::EDGetTokenT<std::vector<reco::ForwardProton>> recoProtonSingleRPToken_;
0049
0050   std::string lhcInfoLabel_;
0051
0052   double maxDiffxi_;
0053   double maxDiffm_;
0054   double maxDiffy_;
0055
0056   unsigned int n_jets_;
0057
0058   bool do_xi_;
0059   bool do_my_;
0060 };
0061
0062 // fill descriptions
0063 //
0064 void HLTPPSJetComparisonFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0065   edm::ParameterSetDescription desc;
0066
0067   desc.add<edm::InputTag>("jetInputTag", edm::InputTag("hltAK4PFJetsCorrected"))
0068       ->setComment("input tag of the jet track collection");
0069   desc.add<edm::InputTag>("forwardProtonInputTag", edm::InputTag("ctppsProtons", "singleRP"))
0070       ->setComment("input tag of the forward proton collection");
0071
0072   desc.add<std::string>("lhcInfoLabel", std::string(""))->setComment("label used for LHCInfo");
0073
0074   desc.add<double>("maxDiffxi", 1.)
0075       ->setComment("maximum relative deviation of RP xi from dijet xi. Used with do_xi option");
0076   desc.add<double>("maxDiffm", 1.)
0077       ->setComment("maximum relative deviation of RP m from dijet m- Used with do_my option");
0078   desc.add<double>("maxDiffy", 1.)
0079       ->setComment("maximum absolute deviation of RP y from dijet y. Used with do_my option");
0080
0081   desc.add<unsigned int>("nJets", 2)->setComment("number of jets to be used");
0082
0083   desc.add<bool>("do_xi", true)->setComment("flag to require xi matching");
0084   desc.add<bool>("do_my", false)->setComment("flag to require m,y matching");
0085
0086   descriptions.addWithDefaultLabel(desc);
0087   return;
0088 }
0089
0090 // destructor and constructor
0091 //
0092 HLTPPSJetComparisonFilter::~HLTPPSJetComparisonFilter() = default;
0093
0094 HLTPPSJetComparisonFilter::HLTPPSJetComparisonFilter(const edm::ParameterSet &iConfig)
0095     : lhcInfoRcdToken_(esConsumes()),
0096       jetInputTag_(iConfig.getParameter<edm::InputTag>("jetInputTag")),
0097       jet_token_(consumes<reco::PFJetCollection>(jetInputTag_)),
0098
0099       forwardProtonInputTag_(iConfig.getParameter<edm::InputTag>("forwardProtonInputTag")),
0100       recoProtonSingleRPToken_(consumes<std::vector<reco::ForwardProton>>(forwardProtonInputTag_)),
0101
0102       lhcInfoLabel_(iConfig.getParameter<std::string>("lhcInfoLabel")),
0103
0104       maxDiffxi_(iConfig.getParameter<double>("maxDiffxi")),
0105       maxDiffm_(iConfig.getParameter<double>("maxDiffm")),
0106       maxDiffy_(iConfig.getParameter<double>("maxDiffy")),
0107
0108       n_jets_(iConfig.getParameter<unsigned int>("nJets")),
0109
0110       do_xi_(iConfig.getParameter<bool>("do_xi")),
0111       do_my_(iConfig.getParameter<bool>("do_my")) {}
0112
0113 // member functions
0114 //
0115 bool HLTPPSJetComparisonFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0116   auto const &hLHCInfo = iSetup.getHandle(lhcInfoRcdToken_);
0117   float sqs = 2. * hLHCInfo->energy();  // get sqrt(s)
0118
0119   edm::Handle<reco::PFJetCollection> jets;
0120   iEvent.getByToken(jet_token_, jets);  // get jet collection
0121
0122   edm::Handle<std::vector<reco::ForwardProton>> recoSingleRPProtons;
0123   iEvent.getByToken(recoProtonSingleRPToken_, recoSingleRPProtons);  // get RP proton collection
0124
0125   if (jets->size() < n_jets_)
0126     return false;  // test for nr jets
0127
0128   if (do_xi_ && maxDiffxi_ > 0) {  // xi matching bloc
0129
0130     float sum45 = 0, sum56 = 0;
0131
0132     for (unsigned int i = 0; i < n_jets_; i++) {
0133       sum45 += (*jets)[i].energy() + (*jets)[i].pz();
0134       sum56 += (*jets)[i].energy() - (*jets)[i].pz();
0135     }
0136
0137     const float xi45 = sum45 / sqs;  // get arm45 xi for n leading-pT jets
0138     const float xi56 = sum56 / sqs;  // get arm56 xi for n leading-pT jets
0139
0140     float min45 = 1000., min56 = 1000.;
0141
0142     for (const auto &proton : *recoSingleRPProtons)  // cycle over proton tracks
0143     {
0144       if (proton.validFit())  // Check that the track fit is valid
0145       {
0146         const auto &xi = proton.xi();  // Get the proton xi
0147
0148         CTPPSDetId rpId(
0149             (*proton.contributingLocalTracks().begin())->rpId());  // get RP ID (rpId.arm() is 0 for 45 and 1 for 56)
0150
0151         if (rpId.arm() == 0 && std::abs(xi - xi45) < min45)
0152           min45 = std::abs(xi - xi45);
0153         if (rpId.arm() == 1 && std::abs(xi - xi56) < min56)
0154           min56 = std::abs(xi - xi56);
0155       }
0156     }
0157
0158     if (min56 / xi56 > maxDiffxi_ || min45 / xi45 > maxDiffxi_)
0159       return false;  // fail cond for xi matching
0160   }
0161
0162   if (do_my_) {  // m, y matching bloc
0163
0164     // get the mass and rap of the n jets
0165     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float>> j_sum;
0166     for (unsigned int i = 0; i < n_jets_; i++)
0167       j_sum = j_sum + (*jets)[i].p4();
0168
0169     const auto &mjet = j_sum.M();
0170     const auto &yjet = j_sum.Rapidity();
0171
0172     for (const auto &proton1 : *recoSingleRPProtons)  // cycle over first RP (only arm45)
0173     {
0174       if (proton1.validFit()) {
0175         CTPPSDetId rpId1(
0176             (*proton1.contributingLocalTracks().begin())->rpId());  // get RP ID (rpId.arm() is 0 for 45 and 1 for 56)
0177         if (rpId1.arm() == 0) {
0178           const auto &xi_45 = proton1.xi();
0179
0180           for (const auto &proton2 : *recoSingleRPProtons)  // cycle over second RP (only arm56)
0181           {
0182             if (proton2.validFit()) {
0183               CTPPSDetId rpId2((*proton2.contributingLocalTracks().begin())->rpId());
0184               if (rpId2.arm() == 1) {
0185                 const auto &xi_56 = proton2.xi();
0186
0187                 // m, y matching tests
0188                 const auto &m = sqs * sqrt(xi_45 * xi_56);
0189                 const auto &y = 0.5 * log(xi_45 / xi_56);
0190                 if ((std::abs(m - mjet) / mjet < maxDiffm_ || maxDiffm_ <= 0) &&
0191                     (std::abs(y - yjet) < maxDiffy_ || maxDiffy_ <= 0))
0192                   return true;  // pass cond, immediately return true
0193               }
0194             }
0195           }
0196         }
0197       }
0198     }
0199     return false;  // fail cond for m,y matching (pass cond never met in cycle)
0200   }
0201
0202   return true;  // if none of the fail conds are met, event has passed the trigger
0203 }
0204
0205 DEFINE_FWK_MODULE(HLTPPSJetComparisonFilter);