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 "CondTools/RunInfo/interface/LHCInfoCombined.h"
0013 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0014 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0015 #include "DataFormats/JetReco/interface/PFJetCollection.h"
0016 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Framework/interface/global/EDFilter.h"
0021 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/isFinite.h"
0024
0025 // class declaration
0026 //
0027 class HLTPPSJetComparisonFilter : public edm::global::EDFilter<edm::LuminosityBlockCache<float>> {
0028 public:
0029   explicit HLTPPSJetComparisonFilter(const edm::ParameterSet &);
0030   ~HLTPPSJetComparisonFilter() override = default;
0031
0032   static void fillDescriptions(edm::ConfigurationDescriptions &);
0033
0034   std::shared_ptr<float> globalBeginLuminosityBlock(const edm::LuminosityBlock &,
0035                                                     const edm::EventSetup &) const override;
0036
0037   void globalEndLuminosityBlock(const edm::LuminosityBlock &, const edm::EventSetup &) const override;
0038
0039   bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override;
0040
0041 private:
0042   // ----------member data ---------------------------
0043   const edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoToken_;
0044   const edm::ESGetToken<LHCInfoPerLS, LHCInfoPerLSRcd> lhcInfoPerLSToken_;
0045   const edm::ESGetToken<LHCInfoPerFill, LHCInfoPerFillRcd> lhcInfoPerFillToken_;
0046
0047   const bool useNewLHCInfo_;
0048
0049   const edm::InputTag jetInputTag_;  // Input tag identifying the jet track
0050   const edm::EDGetTokenT<reco::PFJetCollection> jet_token_;
0051
0052   const edm::InputTag forwardProtonInputTag_;  // Input tag identifying the forward proton collection
0053   const edm::EDGetTokenT<std::vector<reco::ForwardProton>> recoProtonSingleRPToken_;
0054
0055   const double maxDiffxi_;
0056   const double maxDiffm_;
0057   const double maxDiffy_;
0058
0059   const unsigned int n_jets_;
0060
0061   const bool do_xi_;
0062   const bool do_my_;
0063 };
0064
0065 // fill descriptions
0066 //
0067 void HLTPPSJetComparisonFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0068   edm::ParameterSetDescription desc;
0069
0070   desc.add<edm::InputTag>("jetInputTag", edm::InputTag("hltAK4PFJetsCorrected"))
0071       ->setComment("input tag of the jet track collection");
0072   desc.add<edm::InputTag>("forwardProtonInputTag", edm::InputTag("ctppsProtons", "singleRP"))
0073       ->setComment("input tag of the forward proton collection");
0074
0075   desc.add<std::string>("lhcInfoLabel", "")->setComment("label used for LHCInfo");
0076   desc.add<std::string>("lhcInfoPerLSLabel", "")->setComment("label of the LHCInfoPerLS record");
0077   desc.add<std::string>("lhcInfoPerFillLabel", "")->setComment("label of the LHCInfoPerFill record");
0078   desc.add<bool>("useNewLHCInfo", true)->setComment("flag whether to use new LHCInfoPer* records or old LHCInfo");
0079
0080   desc.add<double>("maxDiffxi", 1.)
0081       ->setComment("maximum relative deviation of RP xi from dijet xi. Used with do_xi option");
0082   desc.add<double>("maxDiffm", 1.)
0083       ->setComment("maximum relative deviation of RP m from dijet m- Used with do_my option");
0084   desc.add<double>("maxDiffy", 1.)
0085       ->setComment("maximum absolute deviation of RP y from dijet y. Used with do_my option");
0086
0087   desc.add<unsigned int>("nJets", 2)->setComment("number of jets to be used");
0088
0089   desc.add<bool>("do_xi", true)->setComment("flag to require xi matching");
0090   desc.add<bool>("do_my", false)->setComment("flag to require m,y matching");
0091
0092   descriptions.addWithDefaultLabel(desc);
0093   return;
0094 }
0095
0096 std::shared_ptr<float> HLTPPSJetComparisonFilter::globalBeginLuminosityBlock(const edm::LuminosityBlock &,
0097                                                                              const edm::EventSetup &iSetup) const {
0098   auto cache = std::make_shared<float>();
0099
0100   LHCInfoCombined lhcInfoCombined(iSetup, lhcInfoPerLSToken_, lhcInfoPerFillToken_, lhcInfoToken_, useNewLHCInfo_);
0101   float sqs = 2. * lhcInfoCombined.energy;
0102
0103   if (sqs == 0.f || !edm::isFinite(sqs)) {
0104     edm::LogError("HLTPPSJetComparisonFilter")
0105         << "LHC energy is zero (sqrt(s) = 0). All events in this IOV will be rejected.";
0106   }
0107
0108   *cache = sqs;
0109   return cache;
0110 }
0111
0112 void HLTPPSJetComparisonFilter::globalEndLuminosityBlock(const edm::LuminosityBlock &, const edm::EventSetup &) const {}
0113
0114 HLTPPSJetComparisonFilter::HLTPPSJetComparisonFilter(const edm::ParameterSet &iConfig)
0115     : lhcInfoToken_(esConsumes<LHCInfo, LHCInfoRcd, edm::Transition::BeginLuminosityBlock>(
0116           edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0117       lhcInfoPerLSToken_(esConsumes<LHCInfoPerLS, LHCInfoPerLSRcd, edm::Transition::BeginLuminosityBlock>(
0118           edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerLSLabel")))),
0119       lhcInfoPerFillToken_(esConsumes<LHCInfoPerFill, LHCInfoPerFillRcd, edm::Transition::BeginLuminosityBlock>(
0120           edm::ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerFillLabel")))),
0121       useNewLHCInfo_(iConfig.getParameter<bool>("useNewLHCInfo")),
0122
0123       jetInputTag_(iConfig.getParameter<edm::InputTag>("jetInputTag")),
0124       jet_token_(consumes<reco::PFJetCollection>(jetInputTag_)),
0125
0126       forwardProtonInputTag_(iConfig.getParameter<edm::InputTag>("forwardProtonInputTag")),
0127       recoProtonSingleRPToken_(consumes<std::vector<reco::ForwardProton>>(forwardProtonInputTag_)),
0128
0129       maxDiffxi_(iConfig.getParameter<double>("maxDiffxi")),
0130       maxDiffm_(iConfig.getParameter<double>("maxDiffm")),
0131       maxDiffy_(iConfig.getParameter<double>("maxDiffy")),
0132
0133       n_jets_(iConfig.getParameter<unsigned int>("nJets")),
0134
0135       do_xi_(iConfig.getParameter<bool>("do_xi")),
0136       do_my_(iConfig.getParameter<bool>("do_my")) {}
0137
0138 // member functions
0139 //
0140 bool HLTPPSJetComparisonFilter::filter(edm::StreamID, edm::Event &iEvent, const edm::EventSetup &iSetup) const {
0141   // get the cached value of the LHC energy from the cache
0142   const float sqs = *luminosityBlockCache(iEvent.getLuminosityBlock().index());
0143
0144   // early return in case of not physical energy
0145   if (sqs == 0.f || !edm::isFinite(sqs))
0146     return false;
0147
0148   edm::Handle<reco::PFJetCollection> jets;
0149   iEvent.getByToken(jet_token_, jets);  // get jet collection
0150
0151   edm::Handle<std::vector<reco::ForwardProton>> recoSingleRPProtons;
0152   iEvent.getByToken(recoProtonSingleRPToken_, recoSingleRPProtons);  // get RP proton collection
0153
0154   if (jets->size() < n_jets_)
0155     return false;  // test for nr jets
0156
0157   if (do_xi_ && maxDiffxi_ > 0) {  // xi matching bloc
0158
0159     float sum45 = 0, sum56 = 0;
0160
0161     for (unsigned int i = 0; i < n_jets_; i++) {
0162       sum45 += (*jets)[i].energy() + (*jets)[i].pz();
0163       sum56 += (*jets)[i].energy() - (*jets)[i].pz();
0164     }
0165
0166     const float xi45 = sum45 / sqs;  // get arm45 xi for n leading-pT jets
0167     const float xi56 = sum56 / sqs;  // get arm56 xi for n leading-pT jets
0168
0169     float min45 = 1000., min56 = 1000.;
0170
0171     for (const auto &proton : *recoSingleRPProtons)  // cycle over proton tracks
0172     {
0173       if (proton.validFit())  // Check that the track fit is valid
0174       {
0175         const auto &xi = proton.xi();  // Get the proton xi
0176
0177         CTPPSDetId rpId(
0178             (*proton.contributingLocalTracks().begin())->rpId());  // get RP ID (rpId.arm() is 0 for 45 and 1 for 56)
0179
0180         if (rpId.arm() == 0 && std::abs(xi - xi45) < min45)
0181           min45 = std::abs(xi - xi45);
0182         if (rpId.arm() == 1 && std::abs(xi - xi56) < min56)
0183           min56 = std::abs(xi - xi56);
0184       }
0185     }
0186
0187     if (min56 / xi56 > maxDiffxi_ || min45 / xi45 > maxDiffxi_)
0188       return false;  // fail cond for xi matching
0189   }
0190
0191   if (do_my_) {  // m, y matching bloc
0192
0193     // get the mass and rap of the n jets
0194     ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float>> j_sum;
0195     for (unsigned int i = 0; i < n_jets_; i++)
0196       j_sum = j_sum + (*jets)[i].p4();
0197
0198     const auto &mjet = j_sum.M();
0199     const auto &yjet = j_sum.Rapidity();
0200
0201     for (const auto &proton1 : *recoSingleRPProtons)  // cycle over first RP (only arm45)
0202     {
0203       if (proton1.validFit()) {
0204         CTPPSDetId rpId1(
0205             (*proton1.contributingLocalTracks().begin())->rpId());  // get RP ID (rpId.arm() is 0 for 45 and 1 for 56)
0206         if (rpId1.arm() == 0) {
0207           const auto &xi_45 = proton1.xi();
0208
0209           for (const auto &proton2 : *recoSingleRPProtons)  // cycle over second RP (only arm56)
0210           {
0211             if (proton2.validFit()) {
0212               CTPPSDetId rpId2((*proton2.contributingLocalTracks().begin())->rpId());
0213               if (rpId2.arm() == 1) {
0214                 const auto &xi_56 = proton2.xi();
0215
0216                 // m, y matching tests
0217                 const auto &m = sqs * sqrt(xi_45 * xi_56);
0218                 const auto &y = 0.5 * log(xi_45 / xi_56);
0219                 if ((std::abs(m - mjet) / mjet < maxDiffm_ || maxDiffm_ <= 0) &&
0220                     (std::abs(y - yjet) < maxDiffy_ || maxDiffy_ <= 0))
0221                   return true;  // pass cond, immediately return true
0222               }
0223             }
0224           }
0225         }
0226       }
0227     }
0228     return false;  // fail cond for m,y matching (pass cond never met in cycle)
0229   }
0230
0231   return true;  // if none of the fail conds are met, event has passed the trigger
0232 }
0233
0234 DEFINE_FWK_MODULE(HLTPPSJetComparisonFilter);