JetProducers/plugins/SubEventGenJetProducer.cc

0001
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "RecoJets/JetProducers/plugins/SubEventGenJetProducer.h"
0004
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "FWCore/Utilities/interface/isFinite.h"
0007 #include "RecoJets/JetProducers/interface/JetSpecific.h"
0008 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0009
0010 using namespace std;
0011 using namespace reco;
0012 using namespace edm;
0013 using namespace cms;
0014
0015 namespace {
0016   bool checkHydro(const reco::GenParticle* p) {
0017     const Candidate* m1 = p->mother();
0018     while (m1) {
0019       int pdg = abs(m1->pdgId());
0020       int st = m1->status();
0021       LogDebug("SubEventMothers") << "Pdg ID : " << pdg << endl;
0022       if (st == 3 || pdg < 9 || pdg == 21) {
0023         LogDebug("SubEventMothers") << "Sub-Collision  Found! Pdg ID : " << pdg << endl;
0024         return false;
0025       }
0026       const Candidate* m = m1->mother();
0027       m1 = m;
0028       if (!m1)
0029         LogDebug("SubEventMothers") << "No Mother, particle is : " << pdg << " with status " << st << endl;
0030     }
0031     //      return true;
0032     return true;  // Debugging - to be changed
0033   }
0034 }  // namespace
0035
0036 SubEventGenJetProducer::SubEventGenJetProducer(edm::ParameterSet const& conf) : VirtualJetProducer(conf) {
0037   ignoreHydro_ = conf.getUntrackedParameter<bool>("ignoreHydro", true);
0038
0039   // the subjet collections are set through the config file in the "jetCollInstanceName" field.
0040
0041   input_cand_token_ = consumes<reco::CandidateView>(src_);
0042 }
0043
0044 void SubEventGenJetProducer::inputTowers() {
0045   std::vector<edm::Ptr<reco::Candidate>>::const_iterator inBegin = inputs_.begin(), inEnd = inputs_.end(), i = inBegin;
0046   for (; i != inEnd; ++i) {
0047     reco::CandidatePtr input = inputs_[i - inBegin];
0048     if (edm::isNotFinite(input->pt()))
0049       continue;
0050     if (input->et() < inputEtMin_)
0051       continue;
0052     if (input->energy() < inputEMin_)
0053       continue;
0054     if (isAnomalousTower(input))
0055       continue;
0056
0057     edm::Ptr<reco::Candidate> p = inputs_[i - inBegin];
0058     const GenParticle* pref = dynamic_cast<const GenParticle*>(p.get());
0059     int subevent = pref->collisionId();
0060     LogDebug("SubEventContainers") << "SubEvent is : " << subevent << endl;
0061     LogDebug("SubEventContainers") << "SubSize is : " << subInputs_.size() << endl;
0062
0063     if (subevent >= (int)subInputs_.size()) {
0064       hydroTag_.resize(subevent + 1, -1);
0065       subInputs_.resize(subevent + 1);
0066       LogDebug("SubEventContainers") << "SubSize is : " << subInputs_.size() << endl;
0067       LogDebug("SubEventContainers") << "HydroTagSize is : " << hydroTag_.size() << endl;
0068     }
0069
0070     LogDebug("SubEventContainers") << "HydroTag is : " << hydroTag_[subevent] << endl;
0071     if (hydroTag_[subevent] != 0)
0072       hydroTag_[subevent] = (int)checkHydro(pref);
0073
0074     subInputs_[subevent].push_back(fastjet::PseudoJet(input->px(), input->py(), input->pz(), input->energy()));
0075
0076     subInputs_[subevent].back().set_user_index(i - inBegin);
0077   }
0078 }
0079
0080 void SubEventGenJetProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0081   LogDebug("VirtualJetProducer") << "Entered produce\n";
0082
0083   fjJets_.clear();
0084   subInputs_.clear();
0085   nSubParticles_.clear();
0086   hydroTag_.clear();
0087   inputs_.clear();
0088
0089   // get inputs and convert them to the fastjet format (fastjet::PeudoJet)
0090   edm::Handle<reco::CandidateView> inputsHandle;
0091   iEvent.getByToken(input_cand_token_, inputsHandle);
0092   for (size_t i = 0; i < inputsHandle->size(); ++i) {
0093     inputs_.push_back(inputsHandle->ptrAt(i));
0094   }
0095   LogDebug("VirtualJetProducer") << "Got inputs\n";
0096
0097   inputTowers();
0098   // Convert candidates to fastjet::PseudoJets.
0099   // Also correct to Primary Vertex. Will modify fjInputs_
0100   // and use inputs_
0101
0102   ////////////////
0103
0104   jets_ = std::make_unique<std::vector<GenJet>>();
0105
0106   LogDebug("VirtualJetProducer") << "Inputted towers\n";
0107
0108   size_t nsub = subInputs_.size();
0109
0110   for (size_t isub = 0; isub < nsub; ++isub) {
0111     if (ignoreHydro_ && hydroTag_[isub])
0112       continue;
0113     fjJets_.clear();
0114     fjInputs_.clear();
0115     fjInputs_ = subInputs_[isub];
0116     runAlgorithm(iEvent, iSetup);
0117   }
0118
0119   //Finalize
0120   LogDebug("SubEventJetProducer") << "Wrote jets\n";
0121
0122   iEvent.put(std::move(jets_));
0123   return;
0124 }
0125
0126 void SubEventGenJetProducer::runAlgorithm(edm::Event& iEvent, edm::EventSetup const& iSetup) {
0127   // run algorithm
0128   fjJets_.clear();
0129
0130   fjClusterSeq_ = std::make_shared<fastjet::ClusterSequence>(fjInputs_, *fjJetDefinition_);
0131   fjJets_ = fastjet::sorted_by_pt(fjClusterSeq_->inclusive_jets(jetPtMin_));
0132
0133   using namespace reco;
0134
0135   for (unsigned int ijet = 0; ijet < fjJets_.size(); ++ijet) {
0136     GenJet jet;
0137     const fastjet::PseudoJet& fjJet = fjJets_[ijet];
0138
0139     std::vector<fastjet::PseudoJet> fjConstituents = sorted_by_pt(fjClusterSeq_->constituents(fjJet));
0140
0141     std::vector<CandidatePtr> constituents = getConstituents(fjConstituents);
0142
0143     double px = fjJet.px();
0144     double py = fjJet.py();
0145     double pz = fjJet.pz();
0146     double E = fjJet.E();
0147     double jetArea = 0.0;
0148     double pu = 0.;
0149
0150     writeSpecific(jet, Particle::LorentzVector(px, py, pz, E), vertex_, constituents);
0151
0152     jet.setJetArea(jetArea);
0153     jet.setPileup(pu);
0154
0155     jets_->push_back(jet);
0156   }
0157 }
0158
0159 DEFINE_FWK_MODULE(SubEventGenJetProducer);