SUSYBSM/src/SUSY_HLT_Razor.cc

0001 #include "DataFormats/HLTReco/interface/TriggerObject.h"
0002 #include "FWCore/Common/interface/TriggerNames.h"
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "HLTriggerOffline/SUSYBSM/interface/SUSY_HLT_Razor.h"
0007
0008 SUSY_HLT_Razor::SUSY_HLT_Razor(const edm::ParameterSet &ps) {
0009   edm::LogInfo("SUSY_HLT_Razor") << "Constructor SUSY_HLT_Razor::SUSY_HLT_Razor " << std::endl;
0010   // Get parameters from configuration file
0011   theTrigSummary_ = consumes<trigger::TriggerEvent>(ps.getParameter<edm::InputTag>("trigSummary"));
0012   theMETCollection_ = consumes<edm::View<reco::MET>>(ps.getParameter<edm::InputTag>("METCollection"));
0013   theHemispheres_ = consumes<std::vector<math::XYZTLorentzVector>>(ps.getParameter<edm::InputTag>("hemispheres"));
0014   triggerResults_ = consumes<edm::TriggerResults>(ps.getParameter<edm::InputTag>("TriggerResults"));
0015   triggerPath_ = ps.getParameter<std::string>("TriggerPath");
0016   triggerFilter_ = ps.getParameter<edm::InputTag>("TriggerFilter");
0017   caloFilter_ = ps.getParameter<edm::InputTag>("CaloFilter");
0018   theJetCollection_ = consumes<edm::View<reco::Jet>>(ps.getParameter<edm::InputTag>("jetCollection"));
0019 }
0020
0021 SUSY_HLT_Razor::~SUSY_HLT_Razor() {
0022   edm::LogInfo("SUSY_HLT_Razor") << "Destructor SUSY_HLT_Razor::~SUSY_HLT_Razor " << std::endl;
0023 }
0024
0025 void SUSY_HLT_Razor::bookHistograms(DQMStore::IBooker &ibooker_, edm::Run const &, edm::EventSetup const &) {
0026   edm::LogInfo("SUSY_HLT_Razor") << "SUSY_HLT_Razor::bookHistograms" << std::endl;
0027   // book at beginRun
0028   bookHistos(ibooker_);
0029 }
0030
0031 void SUSY_HLT_Razor::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0032   edm::ParameterSetDescription desc;
0033   desc.add<edm::InputTag>("jetCollection", edm::InputTag("ak4PFJetsCHS"));
0034   desc.add<edm::InputTag>("CaloFilter", edm::InputTag("hltRsqMR200Rsq0p01MR100Calo", "", "HLT"))
0035       ->setComment("Calo HLT filter module used to save razor variable objects");
0036   desc.add<edm::InputTag>("METCollection", edm::InputTag("pfMet"));
0037   desc.add<edm::InputTag>("hemispheres", edm::InputTag("hemispheres"))
0038       ->setComment("hemisphere jets used to compute razor variables");
0039   desc.add<std::string>("TriggerPath", "HLT_RsqMR300_Rsq0p09_MR200_v")->setComment("trigger path name");
0040   desc.add<edm::InputTag>("TriggerFilter", edm::InputTag("hltRsqMR300Rsq0p09MR200", "", "HLT"))
0041       ->setComment("PF HLT filter module used to save razor variable objects");
0042   desc.add<edm::InputTag>("TriggerResults", edm::InputTag("TriggerResults", "", "HLT"));
0043   desc.add<edm::InputTag>("trigSummary", edm::InputTag("hltTriggerSummaryAOD"));
0044   descriptions.add("SUSY_HLT_Razor_Main", desc);
0045 }
0046
0047 void SUSY_HLT_Razor::analyze(edm::Event const &e, edm::EventSetup const &eSetup) {
0048   edm::LogInfo("SUSY_HLT_Razor") << "SUSY_HLT_Razor::analyze" << std::endl;
0049
0050   using namespace std;
0051   using namespace edm;
0052   using namespace reco;
0053
0054   // get hold of collection of objects
0055   Handle<vector<math::XYZTLorentzVector>> hemispheres;
0056   e.getByToken(theHemispheres_, hemispheres);
0057   // get hold of the MET Collection
0058   edm::Handle<edm::View<reco::MET>> inputMet;
0059   e.getByToken(theMETCollection_, inputMet);
0060   if (!inputMet.isValid()) {
0061     edm::LogError("SUSY_HLT_Razor") << "invalid collection: MET"
0062                                     << "\n";
0063     return;
0064   }
0065   //  edm::Handle<reco::PFJetCollection> jetCollection;
0066   edm::Handle<edm::View<reco::Jet>> jetCollection;
0067   e.getByToken(theJetCollection_, jetCollection);
0068   if (!jetCollection.isValid()) {
0069     edm::LogError("SUSY_HLT_Razor") << "invalid collection: jets"
0070                                     << "\n";
0071     return;
0072   }
0073
0074   // check what is in the menu
0075   edm::Handle<edm::TriggerResults> hltresults;
0076   e.getByToken(triggerResults_, hltresults);
0077   if (!hltresults.isValid()) {
0078     edm::LogError("SUSY_HLT_Razor") << "invalid collection: TriggerResults"
0079                                     << "\n";
0080     return;
0081   }
0082
0083   //-------------------------------
0084   //--- Trigger
0085   //-------------------------------
0086   edm::Handle<trigger::TriggerEvent> triggerSummary;
0087   e.getByToken(theTrigSummary_, triggerSummary);
0088   if (!triggerSummary.isValid()) {
0089     edm::LogError("SUSY_HLT_Razor") << "invalid collection: TriggerSummary"
0090                                     << "\n";
0091     return;
0092   }
0093   // get online objects
0094   // HLTriggerOffline/Egamma/python/TriggerTypeDefs.py contains the trigger
0095   // object IDs
0096   double onlineMR = 0, onlineRsq = 0;
0097   double caloMR = 0,
0098          caloRsq = 0;  // online razor variables computed using calo quantities
0099   size_t filterIndex = triggerSummary->filterIndex(triggerFilter_);
0100   size_t caloFilterIndex = triggerSummary->filterIndex(caloFilter_);
0101   // search for online MR and Rsq objects
0102   trigger::TriggerObjectCollection triggerObjects = triggerSummary->getObjects();
0103   if (!(filterIndex >= triggerSummary->sizeFilters())) {
0104     const trigger::Keys &keys = triggerSummary->filterKeys(filterIndex);
0105     for (size_t j = 0; j < keys.size(); ++j) {
0106       trigger::TriggerObject foundObject = triggerObjects[keys[j]];
0107       if (foundObject.id() == 0) {    // the MET object containing MR and Rsq will show up with ID = 0
0108         onlineMR = foundObject.px();  // razor variables stored in dummy reco::MET objects
0109         onlineRsq = foundObject.py();
0110       }
0111     }
0112   }
0113
0114   // search for calo MR and Rsq objects
0115   if (!(caloFilterIndex >= triggerSummary->sizeFilters())) {
0116     const trigger::Keys &keys = triggerSummary->filterKeys(caloFilterIndex);
0117     for (size_t j = 0; j < keys.size(); ++j) {
0118       trigger::TriggerObject foundObject = triggerObjects[keys[j]];
0119       if (foundObject.id() == 0) {
0120         caloMR = foundObject.px();  // razor variables stored in dummy reco::MET objects
0121         caloRsq = foundObject.py();
0122       }
0123     }
0124   }
0125
0126   bool hasFired = false;
0127   const edm::TriggerNames &trigNames = e.triggerNames(*hltresults);
0128   unsigned int numTriggers = trigNames.size();
0129   for (unsigned int hltIndex = 0; hltIndex < numTriggers; ++hltIndex) {
0130     if (trigNames.triggerName(hltIndex).find(triggerPath_) != std::string::npos && hltresults->wasrun(hltIndex) &&
0131         hltresults->accept(hltIndex)) {
0132       hasFired = true;
0133     }
0134   }
0135
0136   float HT = 0.0;
0137
0138   for (unsigned int i = 0; i < jetCollection->size(); i++) {
0139     if (std::abs(jetCollection->at(i).eta()) < 3.0 && jetCollection->at(i).pt() >= 40.0) {
0140       HT += jetCollection->at(i).pt();
0141     }
0142   }
0143   // for (reco::PFJetCollection::const_iterator i_jet = jetCollection->begin();
0144   // i_jet != jetCollection->end(); ++i_jet){
0145   //    if (i_pfjet->pt() < 40) continue;
0146   //    if (fabs(i_pfjet->eta()) > 3.0) continue;
0147   //    pfHT += i_pfjet->pt();
0148   //}
0149   float MET = (inputMet->front()).pt();
0150
0151   // this part is adapted from HLTRFilter.cc
0152
0153   // check that the input collections are available
0154   if (not hemispheres.isValid()) {
0155     // This is happening many times (which is normal) and it's noisy for the
0156     // output, we removed the error message edm::LogError("SUSY_HLT_Razor") <<
0157     // "Hemisphere object is invalid!" << "\n";
0158     return;
0159   }
0160
0161   if (hasFired) {
0162     if (hemispheres->empty()) {  // the Hemisphere Maker will produce an empty
0163                                  // collection of hemispheres if the number of
0164                                  // jets in the
0165       edm::LogError("SUSY_HLT_Razor") << "Cannot calculate M_R and R^2 because there are too many jets! "
0166                                          "(trigger passed automatically without forming the hemispheres)"
0167                                       << endl;
0168       return;
0169     }
0170
0171     //***********************************
0172     // Calculate R
0173
0174     if (!hemispheres->empty() && hemispheres->size() != 2 && hemispheres->size() != 5 && hemispheres->size() != 10) {
0175       edm::LogError("SUSY_HLT_Razor") << "Invalid hemisphere collection!  hemispheres->size() = " << hemispheres->size()
0176                                       << endl;
0177       return;
0178     }
0179
0180     TLorentzVector ja(hemispheres->at(0).x(), hemispheres->at(0).y(), hemispheres->at(0).z(), hemispheres->at(0).t());
0181     TLorentzVector jb(hemispheres->at(1).x(), hemispheres->at(1).y(), hemispheres->at(1).z(), hemispheres->at(1).t());
0182
0183     // dummy vector (this trigger does not care about muons)
0184     std::vector<math::XYZTLorentzVector> muonVec;
0185
0186     double MR = CalcMR(ja, jb);
0187     double R = CalcR(MR, ja, jb, inputMet, muonVec);
0188     double Rsq = R * R;
0189
0190     if (Rsq > 0.15)
0191       h_mr->Fill(MR);
0192     if (MR > 300)
0193       h_rsq->Fill(Rsq);
0194     h_mrRsq->Fill(MR, Rsq);
0195
0196     h_rsq_loose->Fill(Rsq);
0197
0198     if (Rsq > 0.25) {
0199       h_mr_tight->Fill(MR);
0200     }
0201     if (MR > 400) {
0202       h_rsq_tight->Fill(Rsq);
0203     }
0204
0205     h_ht->Fill(HT);
0206     h_met->Fill(MET);
0207     h_htMet->Fill(HT, MET);
0208
0209     h_online_mr_vs_mr->Fill(MR, onlineMR);
0210     h_online_rsq_vs_rsq->Fill(Rsq, onlineRsq);
0211
0212     h_calo_mr_vs_mr->Fill(MR, caloMR);
0213     h_calo_rsq_vs_rsq->Fill(Rsq, caloRsq);
0214   }
0215 }
0216
0217 void SUSY_HLT_Razor::bookHistos(DQMStore::IBooker &ibooker_) {
0218   ibooker_.cd();
0219   ibooker_.setCurrentFolder("HLT/SUSYBSM/" + triggerPath_);
0220
0221   h_mr = ibooker_.book1D("mr", "M_{R} (R^{2} > 0.15) ; GeV", 100, 0.0, 4000);
0222   h_rsq = ibooker_.book1D("rsq", "R^{2} (M_{R} > 300)", 100, 0.0, 1.5);
0223   h_mrRsq = ibooker_.book2D("mrRsq", "R^{2} vs M_{R}; GeV; ", 100, 0.0, 4000.0, 100, 0.0, 1.5);
0224
0225   h_mr_tight = ibooker_.book1D("mr_tight", "M_{R} (R^{2} > 0.25) ; GeV", 100, 0.0, 4000);
0226   h_rsq_tight = ibooker_.book1D("rsq_tight", "R^{2} (M_{R} > 400) ; ", 100, 0.0, 1.5);
0227
0228   h_rsq_loose = ibooker_.book1D("rsq_loose", "R^{2} (M_{R} > 0) ; ", 100, 0.0, 1.5);
0229
0230   h_ht = ibooker_.book1D("ht", "HT; GeV; ", 100, 0.0, 4000.0);
0231   h_met = ibooker_.book1D("met", "MET; GeV; ", 100, 0.0, 1000);
0232   h_htMet = ibooker_.book2D("htMet", "MET vs HT; GeV; ", 100, 0.0, 4000.0, 100, 0.0, 1000);
0233
0234   h_online_mr_vs_mr = ibooker_.book2D("online_mr_vs_mr",
0235                                       "Online M_{R} vs  Offline M_{R} (events passing "
0236                                       "trigger); Offline M_{R} (GeV); Online M_{R} (GeV); ",
0237                                       100,
0238                                       0.0,
0239                                       4000.0,
0240                                       100,
0241                                       0.0,
0242                                       4000.0);
0243   h_calo_mr_vs_mr = ibooker_.book2D("calo_mr_vs_mr",
0244                                     "Calo M_{R} vs  Offline M_{R} (events passing trigger); "
0245                                     "Offline M_{R} (GeV); Calo M_{R} (GeV); ",
0246                                     100,
0247                                     0.0,
0248                                     4000.0,
0249                                     100,
0250                                     0.0,
0251                                     4000.0);
0252   h_online_rsq_vs_rsq = ibooker_.book2D("online_rsq_vs_rsq",
0253                                         "Online R^{2} vs Offline R^{2} (events passing trigger); "
0254                                         "Offline R^{2}; Online R^{2}; ",
0255                                         100,
0256                                         0.0,
0257                                         1.5,
0258                                         100,
0259                                         0.0,
0260                                         1.5);
0261   h_calo_rsq_vs_rsq = ibooker_.book2D("calo_rsq_vs_rsq",
0262                                       "Calo R^{2} vs Offline R^{2} (events passing trigger); "
0263                                       "Offline R^{2}; Calo R^{2}; ",
0264                                       100,
0265                                       0.0,
0266                                       1.5,
0267                                       100,
0268                                       0.0,
0269                                       1.5);
0270
0271   ibooker_.cd();
0272 }
0273
0274 // CalcMR and CalcR borrowed from HLTRFilter.cc
0275 double SUSY_HLT_Razor::CalcMR(TLorentzVector ja, TLorentzVector jb) {
0276   if (ja.Pt() <= 0.1)
0277     return -1;
0278
0279   ja.SetPtEtaPhiM(ja.Pt(), ja.Eta(), ja.Phi(), 0.0);
0280   jb.SetPtEtaPhiM(jb.Pt(), jb.Eta(), jb.Phi(), 0.0);
0281
0282   if (ja.Pt() > jb.Pt()) {
0283     TLorentzVector temp = ja;
0284     ja = jb;
0285     jb = temp;
0286   }
0287
0288   double A = ja.P();
0289   double B = jb.P();
0290   double az = ja.Pz();
0291   double bz = jb.Pz();
0292   TVector3 jaT, jbT;
0293   jaT.SetXYZ(ja.Px(), ja.Py(), 0.0);
0294   jbT.SetXYZ(jb.Px(), jb.Py(), 0.0);
0295   double ATBT = (jaT + jbT).Mag2();
0296
0297   double MR = sqrt((A + B) * (A + B) - (az + bz) * (az + bz) -
0298                    (jbT.Dot(jbT) - jaT.Dot(jaT)) * (jbT.Dot(jbT) - jaT.Dot(jaT)) / (jaT + jbT).Mag2());
0299
0300   double mybeta = (jbT.Dot(jbT) - jaT.Dot(jaT)) / sqrt(ATBT * ((A + B) * (A + B) - (az + bz) * (az + bz)));
0301
0302   double mygamma = 1. / sqrt(1. - mybeta * mybeta);
0303
0304   // use gamma times MRstar
0305   return MR * mygamma;
0306 }
0307
0308 double SUSY_HLT_Razor::CalcR(double MR,
0309                              TLorentzVector ja,
0310                              TLorentzVector jb,
0311                              edm::Handle<edm::View<reco::MET>> inputMet,
0312                              const std::vector<math::XYZTLorentzVector> &muons) {
0313   // now we can calculate MTR
0314   TVector3 met;
0315   met.SetPtEtaPhi((inputMet->front()).pt(), 0.0, (inputMet->front()).phi());
0316
0317   std::vector<math::XYZTLorentzVector>::const_iterator muonIt;
0318   for (muonIt = muons.begin(); muonIt != muons.end(); muonIt++) {
0319     TVector3 tmp;
0320     tmp.SetPtEtaPhi(muonIt->pt(), 0, muonIt->phi());
0321     met -= tmp;
0322   }
0323
0324   double MTR = sqrt(0.5 * (met.Mag() * (ja.Pt() + jb.Pt()) - met.Dot(ja.Vect() + jb.Vect())));
0325
0326   // filter events
0327   return float(MTR) / float(MR);  // R
0328 }
0329
0330 // define this as a plug-in
0331 DEFINE_FWK_MODULE(SUSY_HLT_Razor);