PatAlgos/plugins/PATHemisphereProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    PhysicsTools/PatAlgos
0004 // Class:      PATHemisphereProducer
0005 //
0006 /**\class PATHemisphereProducer PATHemisphereProducer.cc "PhysicsTools/PatAlgos/plugins/PATHemisphereProducer.cc"
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Authors:  Christian Autermann, Tanja Rommerskirchen
0015 //          Created:  Sat Mar 22 12:58:04 CET 2008
0016 //
0017 //
0018
0019 #include "DataFormats/PatCandidates/interface/Jet.h"
0020 #include "DataFormats/PatCandidates/interface/MET.h"
0021 #include "DataFormats/PatCandidates/interface/Muon.h"
0022 #include "DataFormats/PatCandidates/interface/Electron.h"
0023 #include "DataFormats/PatCandidates/interface/Photon.h"
0024 #include "DataFormats/PatCandidates/interface/Tau.h"
0025 #include "DataFormats/PatCandidates/interface/Hemisphere.h"
0026 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0027 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0028 #include "DataFormats/TrackReco/interface/Track.h"
0029 #include "DataFormats/Candidate/interface/Particle.h"
0030 #include "FWCore/Framework/interface/global/EDProducer.h"
0031 #include "FWCore/Framework/interface/Event.h"
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0034 #include "PhysicsTools/UtilAlgos/interface/ParameterAdapter.h"
0035 #include "PhysicsTools/PatAlgos/interface/HemisphereAlgo.h"
0036
0037 #include <map>
0038 #include <memory>
0039 #include <utility>
0040 #include <vector>
0041
0042 class PATHemisphereProducer : public edm::global::EDProducer<> {
0043 public:
0044   explicit PATHemisphereProducer(const edm::ParameterSet&);
0045   ~PATHemisphereProducer() override;
0046
0047   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0048
0049 private:
0050   // ----------member data ---------------------------
0051   /// Input: All PAT objects that are to cross-clean  or needed for that
0052   const edm::EDGetTokenT<reco::CandidateView> _patJetsToken;
0053   //       edm::EDGetTokenT<reco::CandidateView> _patMetsToken;
0054   const edm::EDGetTokenT<reco::CandidateView> _patMuonsToken;
0055   const edm::EDGetTokenT<reco::CandidateView> _patElectronsToken;
0056   const edm::EDGetTokenT<reco::CandidateView> _patPhotonsToken;
0057   const edm::EDGetTokenT<reco::CandidateView> _patTausToken;
0058
0059   const float _minJetEt;
0060   const float _minMuonEt;
0061   const float _minElectronEt;
0062   const float _minTauEt;
0063   const float _minPhotonEt;
0064
0065   const float _maxJetEta;
0066   const float _maxMuonEta;
0067   const float _maxElectronEta;
0068   const float _maxTauEta;
0069   const float _maxPhotonEta;
0070
0071   const int _seedMethod;
0072   const int _combinationMethod;
0073
0074   typedef std::vector<float> HemiAxis;
0075 };
0076
0077 using namespace pat;
0078
0079 //
0080 // constants, enums and typedefs
0081 //
0082
0083 //
0084 // static data member definitions
0085 //
0086
0087 //
0088 // constructors and destructor
0089 //
0090 PATHemisphereProducer::PATHemisphereProducer(const edm::ParameterSet& iConfig)
0091     : _patJetsToken(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("patJets"))),
0092       _patMuonsToken(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("patMuons"))),
0093       _patElectronsToken(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("patElectrons"))),
0094       _patPhotonsToken(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("patPhotons"))),
0095       _patTausToken(consumes<reco::CandidateView>(iConfig.getParameter<edm::InputTag>("patTaus"))),
0096
0097       _minJetEt(iConfig.getParameter<double>("minJetEt")),
0098       _minMuonEt(iConfig.getParameter<double>("minMuonEt")),
0099       _minElectronEt(iConfig.getParameter<double>("minElectronEt")),
0100       _minTauEt(iConfig.getParameter<double>("minTauEt")),
0101       _minPhotonEt(iConfig.getParameter<double>("minPhotonEt")),
0102
0103       _maxJetEta(iConfig.getParameter<double>("maxJetEta")),
0104       _maxMuonEta(iConfig.getParameter<double>("maxMuonEta")),
0105       _maxElectronEta(iConfig.getParameter<double>("maxElectronEta")),
0106       _maxTauEta(iConfig.getParameter<double>("maxTauEta")),
0107       _maxPhotonEta(iConfig.getParameter<double>("maxPhotonEta")),
0108
0109       _seedMethod(iConfig.getParameter<int>("seedMethod")),
0110       _combinationMethod(iConfig.getParameter<int>("combinationMethod"))
0111
0112 {
0113   produces<std::vector<pat::Hemisphere>>();
0114 }
0115
0116 PATHemisphereProducer::~PATHemisphereProducer() {
0117   // do anything here that needs to be done at desctruction time
0118   // (e.g. close files, deallocate resources etc.)
0119 }
0120
0121 //
0122 // member functions
0123 //
0124
0125 // ------------ method called to produce the data  ------------
0126 void PATHemisphereProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0127   using namespace edm;
0128   using namespace std;
0129
0130   std::vector<float> vPx, vPy, vPz, vE;
0131   std::vector<float> vA1, vA2;
0132   std::vector<int> vgroups;
0133   std::vector<reco::CandidatePtr> componentPtrs;
0134
0135   //Jets
0136   Handle<reco::CandidateView> pJets;
0137   iEvent.getByToken(_patJetsToken, pJets);
0138
0139   //Muons
0140   Handle<reco::CandidateView> pMuons;
0141   iEvent.getByToken(_patMuonsToken, pMuons);
0142
0143   //Electrons
0144   Handle<reco::CandidateView> pElectrons;
0145   iEvent.getByToken(_patElectronsToken, pElectrons);
0146
0147   //Photons
0148   Handle<reco::CandidateView> pPhotons;
0149   iEvent.getByToken(_patPhotonsToken, pPhotons);
0150
0151   //Taus
0152   Handle<reco::CandidateView> pTaus;
0153   iEvent.getByToken(_patTausToken, pTaus);
0154
0155   //fill e,p vector with information from all objects (hopefully cleaned before)
0156   for (int i = 0; i < (int)(*pJets).size(); i++) {
0157     if ((*pJets)[i].pt() < _minJetEt || fabs((*pJets)[i].eta()) > _maxJetEta)
0158       continue;
0159
0160     componentPtrs.push_back(pJets->ptrAt(i));
0161   }
0162
0163   for (int i = 0; i < (int)(*pMuons).size(); i++) {
0164     if ((*pMuons)[i].pt() < _minMuonEt || fabs((*pMuons)[i].eta()) > _maxMuonEta)
0165       continue;
0166
0167     componentPtrs.push_back(pMuons->ptrAt(i));
0168   }
0169
0170   for (int i = 0; i < (int)(*pElectrons).size(); i++) {
0171     if ((*pElectrons)[i].pt() < _minElectronEt || fabs((*pElectrons)[i].eta()) > _maxElectronEta)
0172       continue;
0173
0174     componentPtrs.push_back(pElectrons->ptrAt(i));
0175   }
0176
0177   for (int i = 0; i < (int)(*pPhotons).size(); i++) {
0178     if ((*pPhotons)[i].pt() < _minPhotonEt || fabs((*pPhotons)[i].eta()) > _maxPhotonEta)
0179       continue;
0180
0181     componentPtrs.push_back(pPhotons->ptrAt(i));
0182   }
0183
0184   //aren't taus included in jets?
0185   for (int i = 0; i < (int)(*pTaus).size(); i++) {
0186     if ((*pTaus)[i].pt() < _minTauEt || fabs((*pTaus)[i].eta()) > _maxTauEta)
0187       continue;
0188
0189     componentPtrs.push_back(pTaus->ptrAt(i));
0190   }
0191
0192   // create product
0193   auto hemispheres = std::make_unique<std::vector<Hemisphere>>();
0194   hemispheres->reserve(2);
0195
0196   //calls HemiAlgorithm for seed method 3 (transv. inv. Mass) and association method 3 (Lund algo)
0197   HemisphereAlgo myHemi(componentPtrs, _seedMethod, _combinationMethod);
0198
0199   //get Hemisphere Axis
0200   vA1 = myHemi.getAxis1();
0201   vA2 = myHemi.getAxis2();
0202
0203   reco::Particle::LorentzVector p1(vA1[0] * vA1[3], vA1[1] * vA1[3], vA1[2] * vA1[3], vA1[4]);
0204   hemispheres->push_back(Hemisphere(p1));
0205
0206   reco::Particle::LorentzVector p2(vA2[0] * vA2[3], vA2[1] * vA2[3], vA2[2] * vA2[3], vA2[4]);
0207   hemispheres->push_back(Hemisphere(p2));
0208
0209   //get information to which Hemisphere each object belongs
0210   vgroups = myHemi.getGrouping();
0211
0212   for (unsigned int i = 0; i < vgroups.size(); ++i) {
0213     if (vgroups[i] == 1) {
0214       (*hemispheres)[0].addDaughter(componentPtrs[i]);
0215     } else {
0216       (*hemispheres)[1].addDaughter(componentPtrs[i]);
0217     }
0218   }
0219
0220   iEvent.put(std::move(hemispheres));
0221 }
0222
0223 //define this as a plug-in
0224 #include "FWCore/Framework/interface/MakerMacros.h"
0225 DEFINE_FWK_MODULE(PATHemisphereProducer);