QjetsAdder.cc

Go to the documentation of this file.

#include <numeric>
 
#include "RecoJets/JetProducers/interface/QjetsAdder.h"
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequence.hh"
 
#include "FWCore/Framework/interface/MakerMacros.h"
 
using namespace std;
 
void QjetsAdder::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
  // read input collection
  //edm::Handle<edm::View<pat::Jet> > jets;
  edm::Handle<edm::View<reco::Jet>> jets;
  iEvent.getByToken(src_token_, jets);
 
  // prepare room for output
  std::vector<float> QjetsVolatility;
  QjetsVolatility.reserve(jets->size());
 
  for (typename edm::View<reco::Jet>::const_iterator jetIt = jets->begin(); jetIt != jets->end(); ++jetIt) {
    reco::Jet newCand(*jetIt);
 
    if (newCand.pt() < cutoff_) {
      QjetsVolatility.push_back(-1);
      continue;
    }
 
    //refill and recluster
    vector<fastjet::PseudoJet> allconstits;
    //for (unsigned k=0; k < newCand.getPFConstituents().size(); k++){
    for (unsigned k = 0; k < newCand.getJetConstituents().size(); k++) {
      const edm::Ptr<reco::Candidate> thisParticle = newCand.getJetConstituents().at(k);
      allconstits.push_back(
          fastjet::PseudoJet(thisParticle->px(), thisParticle->py(), thisParticle->pz(), thisParticle->energy()));
    }
 
    fastjet::JetDefinition jetDef(fastjet::cambridge_algorithm, jetRad_);
    if (mJetAlgo_ == "AK")
      jetDef.set_jet_algorithm(fastjet::antikt_algorithm);
    else if (mJetAlgo_ == "CA")
      jetDef.set_jet_algorithm(fastjet::cambridge_algorithm);
    else
      throw cms::Exception("GroomedJetFiller") << " unknown jet algorithm " << std::endl;
 
    fastjet::ClusterSequence thisClustering_basic(allconstits, jetDef);
    std::vector<fastjet::PseudoJet> out_jets_basic = sorted_by_pt(thisClustering_basic.inclusive_jets(cutoff_));
    //std::cout << newCand.pt() << " " << out_jets_basic.size() <<std::endl;
    if (out_jets_basic.size() !=
        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
      QjetsVolatility.push_back(-1);
      continue;
    }
 
    // setup objects for qjets computation
    fastjet::JetDefinition qjet_def(&qjetsAlgo_);
 
    std::vector<double> qjetmass;
 
    vector<fastjet::PseudoJet> constits;
    unsigned int nqjetconstits = out_jets_basic.at(0)
                                     .constituents()
                                     .size();  // there should always be exactly one reclsutered jet => always "at(0)"
    if (nqjetconstits < (unsigned int)QjetsPreclustering_)
      constits = out_jets_basic.at(0).constituents();
    else
      constits = out_jets_basic.at(0).associated_cluster_sequence()->exclusive_subjets_up_to(out_jets_basic.at(0),
                                                                                             QjetsPreclustering_);
 
    edm::Service<edm::RandomNumberGenerator> rng;
    CLHEP::HepRandomEngine* engine = &rng->getEngine(iEvent.streamID());
    qjetsAlgo_.SetRNEngine(engine);
    // create probabilistic recusterings
    for (unsigned int ii = 0; ii < (unsigned int)ntrial_; ii++) {
      //qjetsAlgo_.SetRandSeed(iEvent.id().event()*100 + (jetIt - jets->begin())*ntrial_ + ii );// set random seed for reprudcibility. We need a smarted scheme
      fastjet::ClusterSequence qjet_seq(constits, qjet_def);
      vector<fastjet::PseudoJet> inclusive_jets2 = sorted_by_pt(qjet_seq.inclusive_jets(cutoff_));
      if (!inclusive_jets2.empty()) {  // fill the massvalue only if the reclustering was successfull
        qjetmass.push_back(inclusive_jets2[0].m());
      }
 
    }  //end loop on trials
 
    if (qjetmass.empty()) {  //protection against dummy case
      QjetsVolatility.push_back(-1);
      continue;
    }
 
    double mean = std::accumulate(qjetmass.begin(), qjetmass.end(), 0) / qjetmass.size();
    float totalsquared = 0.;
    for (unsigned int i = 0; i < qjetmass.size(); i++) {
      totalsquared += (qjetmass[i] - mean) * (qjetmass[i] - mean);
    }
    float variance = sqrt(totalsquared / qjetmass.size());
 
    QjetsVolatility.push_back(variance / mean);
  }  //end loop on jets
 
  auto outQJV = std::make_unique<edm::ValueMap<float>>();
  edm::ValueMap<float>::Filler fillerQJV(*outQJV);
  fillerQJV.insert(jets, QjetsVolatility.begin(), QjetsVolatility.end());
  fillerQJV.fill();
 
  iEvent.put(std::move(outQJV), "QjetsVolatility");
}
 
DEFINE_FWK_MODULE(QjetsAdder);