JetProducers/plugins/QjetsAdder.cc

0001 #include <numeric>
0002
0003 #include "RecoJets/JetProducers/interface/QjetsAdder.h"
0004 #include "fastjet/PseudoJet.hh"
0005 #include "fastjet/ClusterSequence.hh"
0006
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008
0009 using namespace std;
0010
0011 void QjetsAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0012   // read input collection
0013   //edm::Handle<edm::View<pat::Jet> > jets;
0014   edm::Handle<edm::View<reco::Jet>> jets;
0015   iEvent.getByToken(src_token_, jets);
0016
0017   // prepare room for output
0018   std::vector<float> QjetsVolatility;
0019   QjetsVolatility.reserve(jets->size());
0020
0021   for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
0022     const reco::Jet& newCand(*jetIt);
0023
0024     if (newCand.pt() < cutoff_) {
0025       QjetsVolatility.push_back(-1);
0026       continue;
0027     }
0028
0029     //refill and recluster
0030     vector<fastjet::PseudoJet> allconstits;
0031     //for (unsigned k=0; k < newCand.getPFConstituents().size(); k++){
0032     for (unsigned k = 0; k < newCand.getJetConstituents().size(); k++) {
0033       const edm::Ptr<reco::Candidate> thisParticle = newCand.getJetConstituents().at(k);
0034       allconstits.push_back(
0035           fastjet::PseudoJet(thisParticle->px(), thisParticle->py(), thisParticle->pz(), thisParticle->energy()));
0036     }
0037
0038     fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, jetRad_);
0039     if (mJetAlgo_ == "AK")
0040       jetDef.set_jet_algorithm(fastjet::antikt_algorithm);
0041     else if (mJetAlgo_ == "CA")
0042       jetDef.set_jet_algorithm(fastjet::cambridge_algorithm);
0043     else
0044       throw cms::Exception("GroomedJetFiller") << " unknown jet algorithm " << std::endl;
0045
0046     fastjet::ClusterSequence thisClustering_basic(allconstits, jetDef);
0047     std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(cutoff_));
0048     //std::cout << newCand.pt() << " " << out_jets_basic.size() <<std::endl;
0049     if (out_jets_basic.size() !=
0050         1) {  // jet reclustering did not return exactly 1 jet, most likely due to the higher cutoff or large cone size. Use a recognizeable default value for this jet
0051       QjetsVolatility.push_back(-1);
0052       continue;
0053     }
0054
0055     // setup objects for qjets computation
0056     fastjet::JetDefinition qjet_def(&qjetsAlgo_);
0057
0058     std::vector<double> qjetmass;
0059
0060     vector<fastjet::PseudoJet> constits;
0061     unsigned int nqjetconstits = out_jets_basic.at(0)
0062                                      .constituents()
0063                                      .size();  // there should always be exactly one reclsutered jet => always "at(0)"
0064     if (nqjetconstits < (unsigned int)QjetsPreclustering_)
0065       constits = out_jets_basic.at(0).constituents();
0066     else
0067       constits = out_jets_basic.at(0).associated_cluster_sequence()->exclusive_subjets_up_to(out_jets_basic.at(0),
0068                                                                                              QjetsPreclustering_);
0069
0070     edm::Service<edm::RandomNumberGenerator> rng;
0071     CLHEP::HepRandomEngine* engine = &rng->getEngine(iEvent.streamID());
0072     qjetsAlgo_.SetRNEngine(engine);
0073     // create probabilistic recusterings
0074     for (unsigned int ii = 0; ii < (unsigned int)ntrial_; ii++) {
0075       //qjetsAlgo_.SetRandSeed(iEvent.id().event()*100 + (jetIt - jets->begin())*ntrial_ + ii );// set random seed for reprudcibility. We need a smarted scheme
0076       fastjet::ClusterSequence qjet_seq(constits, qjet_def);
0077       vector<fastjet::PseudoJet> inclusive_jets2 = sorted_by_pt(qjet_seq.inclusive_jets(cutoff_));
0078       if (!inclusive_jets2.empty()) {  // fill the massvalue only if the reclustering was successfull
0079         qjetmass.push_back(inclusive_jets2[0].m());
0080       }
0081
0082     }  //end loop on trials
0083
0084     if (qjetmass.empty()) {  //protection against dummy case
0085       QjetsVolatility.push_back(-1);
0086       continue;
0087     }
0088
0089     double mean = std::accumulate(qjetmass.begin(), qjetmass.end(), 0) / qjetmass.size();
0090     float totalsquared = 0.;
0091     for (unsigned int i = 0; i < qjetmass.size(); i++) {
0092       totalsquared += (qjetmass[i] - mean) * (qjetmass[i] - mean);
0093     }
0094     float variance = sqrt(totalsquared / qjetmass.size());
0095
0096     QjetsVolatility.push_back(variance / mean);
0097   }  //end loop on jets
0098
0099   auto outQJV = std::make_unique<edm::ValueMap<float>>();
0100   edm::ValueMap<float>::Filler fillerQJV(*outQJV);
0101   fillerQJV.insert(jets, QjetsVolatility.begin(), QjetsVolatility.end());
0102   fillerQJV.fill();
0103
0104   iEvent.put(std::move(outQJV), "QjetsVolatility");
0105 }
0106
0107 DEFINE_FWK_MODULE(QjetsAdder);