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);