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