Jettoolbox mini aod 706

PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py

@@ -72,7 +72,7 @@ def miniAOD_customizeCommon(process):
     btagDiscriminators = ['jetBProbabilityBJetTags', 'jetProbabilityBJetTags', 'trackCountingHighPurBJetTags', 'trackCountingHighEffBJetTags', 'simpleSecondaryVertexHighEffBJetTags',
                          'simpleSecondaryVertexHighPurBJetTags', 'combinedSecondaryVertexBJetTags' , 'combinedInclusiveSecondaryVertexBJetTags' ],
     )
-    #add CA8   
+    #add AK8   
     from PhysicsTools.PatAlgos.tools.jetTools import addJetCollection
     addJetCollection(process, labelName = 'AK8', jetSource = cms.InputTag('ak8PFJetsCHS'),algo= 'AK', rParam = 0.8, jetCorrections = ('AK7PFchs', cms.vstring(['L1FastJet', 'L2Relative', 'L3Absolute']), 'None') )
     process.patJetsAK8.userData.userFloats.src = [] # start with empty list of user floats
@@ -100,6 +100,26 @@ def miniAOD_customizeCommon(process):
     process.cmsTopTagPFJetsCHSLinksAK8.matched = cms.InputTag("cmsTopTagPFJetsCHS")
     process.patJetsAK8.userData.userFloats.src += ['cmsTopTagPFJetsCHSLinksAK8']
 
+    #QJetsAdder
+    if hasattr(process,"RandomNumberGeneratorService") :
+	    process.RandomNumberGeneratorService.QJetsAdderAK8 = cms.PSet(initialSeed = cms.untracked.uint32(7))
+	    process.RandomNumberGeneratorService.QJetsAdder = cms.PSet(initialSeed = cms.untracked.uint32(7))
+    else:
+	    process.RandomNumberGeneratorService = cms.Service("RandomNumberGeneratorService", QJetsAdderAK8 = cms.PSet(initialSeed = cms.untracked.uint32(7)), QJetsAdder = cms.PSet(initialSeed = cms.untracked.uint32(7)))
+    process.load('RecoJets.JetProducers.qjetsadder_cfi')
+    process.QJetsAdderAK8 = process.QJetsAdder.clone()
+    process.QJetsAdderAK8.src = cms.InputTag("ak8PFJetsCHS")
+    process.QJetsAdderAK8.jetRad = cms.double(0.8)
+    process.QJetsAdderAK8.jetAlgo = cms.string('AK')
+    process.patJetsAK8.userData.userFloats.src += ['QJetsAdderAK8:QjetsVolatility']
+
+    # add Njetiness
+    process.load('RecoJets.JetProducers.nJettinessAdder_cfi')
+    process.NjettinessAK8 = process.Njettiness.clone()
+    process.NjettinessAK8.src = cms.InputTag("ak8PFJetsCHS")
+    process.NjettinessAK8.cone = cms.double(0.8)
+    process.patJetsAK8.userData.userFloats.src += ['NjettinessAK8:tau1','NjettinessAK8:tau2','NjettinessAK8:tau3']
+
     #
     from PhysicsTools.PatAlgos.tools.trigTools import switchOnTriggerStandAlone
     switchOnTriggerStandAlone( process, outputModule = '' )

RecoJets/JetAlgorithms/interface/Qjets.h

@@ -0,0 +1,64 @@
+#ifndef RecoJets_JetAlgorithms_QJets_h
+#define RecoJets_JetAlgorithms_QJets_h
+#include <queue>
+#include <vector>
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "FWCore/MessageLogger/interface/MessageLogger.h"
+#include "FWCore/Utilities/interface/isFinite.h"
+#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "CLHEP/Random/RandomEngine.h"
+
+struct JetDistance{
+  double dij;
+  int j1;
+  int j2;
+};
+
+class JetDistanceCompare{
+ public:
+  JetDistanceCompare(){};
+  bool operator() (const JetDistance& lhs, const JetDistance&rhs) const {return lhs.dij > rhs.dij;};
+};
+
+class Qjets{
+ private:
+  bool _rand_seed_set;
+  unsigned int _seed;
+  double _zcut, _dcut, _dcut_fctr, _exp_min, _exp_max, _rigidity, _truncation_fctr;
+  std::map<int,bool> _merged_jets;
+  std::priority_queue <JetDistance, std::vector<JetDistance>, JetDistanceCompare> _distances;
+  CLHEP::HepRandomEngine* _rnEngine;
+
+  double d_ij(const fastjet::PseudoJet& v1, const fastjet::PseudoJet& v2) const; 
+  void computeDCut(fastjet::ClusterSequence & cs);
+
+  double Rand();
+  bool Prune(JetDistance& jd,fastjet::ClusterSequence & cs);
+  bool JetsUnmerged(const JetDistance& jd) const;
+  bool JetUnmerged(int num) const;
+  void ComputeNewDistanceMeasures(fastjet::ClusterSequence & cs, unsigned int new_jet);
+  void ComputeAllDistances(const std::vector<fastjet::PseudoJet>& inp);  
+  double ComputeMinimumDistance();
+  double ComputeNormalization(double dmin);
+  JetDistance GetNextDistance();
+  bool Same(const JetDistance& lhs, const JetDistance& rhs);
+
+ public:
+  Qjets(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr, CLHEP::HepRandomEngine* rnEngine)  : _rand_seed_set(false),
+    _zcut(zcut),
+    _dcut(-1.),
+    _dcut_fctr(dcut_fctr),
+    _exp_min(exp_min),
+    _exp_max(exp_max),
+    _rigidity(rigidity),
+    _truncation_fctr(truncation_fctr),
+    _rnEngine(rnEngine)
+    {};
+
+  void Cluster(fastjet::ClusterSequence & cs);
+  void SetRandSeed(unsigned int seed); /* In case you want reproducible behavior */
+};
+#endif

RecoJets/JetAlgorithms/interface/QjetsPlugin.h

@@ -0,0 +1,33 @@
+#ifndef RecoJets_JetAlgorithms_QJETSPLUGIN_h
+#define RecoJets_JetAlgorithms_QJETSPLUGIN_h
+#include "fastjet/JetDefinition.hh"
+#include "fastjet/PseudoJet.hh"
+#include "fastjet/ClusterSequence.hh"
+#include "RecoJets/JetAlgorithms/interface/Qjets.h"
+
+class QjetsPlugin: public fastjet::JetDefinition::Plugin{
+ private:
+  bool _rand_seed_set;
+  unsigned int _seed;
+  int _truncated_length;
+  double _zcut, _dcut_fctr, _exp_min, _exp_max, _rigidity,_truncation_fctr;
+  CLHEP::HepRandomEngine* _rnEngine;
+ public:
+  QjetsPlugin(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr = 0.)  : _rand_seed_set(false),
+    _zcut(zcut),
+    _dcut_fctr(dcut_fctr),
+    _exp_min(exp_min),
+    _exp_max(exp_max),
+    _rigidity(rigidity),
+    _truncation_fctr(truncation_fctr),
+    _rnEngine(0)
+      {};
+  void SetRandSeed(unsigned int seed); /* In case you want reproducible behavior */
+  void SetRNEngine(CLHEP::HepRandomEngine* rnEngine){
+    _rnEngine=rnEngine;
+  };
+  double R() const;
+  std::string description() const;
+  void run_clustering(fastjet::ClusterSequence & cs) const;
+};
+#endif

RecoJets/JetAlgorithms/src/Qjets.cc

@@ -0,0 +1,170 @@
+#include "RecoJets/JetAlgorithms/interface/Qjets.h"
+
+using namespace std;
+
+void Qjets::SetRandSeed(unsigned int seed){
+  _rand_seed_set = true;
+  _seed = seed;
+}
+
+bool Qjets::JetUnmerged(int num) const{
+  return _merged_jets.find(num) == _merged_jets.end();
+}
+
+bool Qjets::JetsUnmerged(const JetDistance& jd) const{
+  return JetUnmerged(jd.j1) && JetUnmerged(jd.j2);
+}
+
+JetDistance Qjets::GetNextDistance(){
+  vector< pair<JetDistance, double> > popped_distances;
+  double norm(0.);
+  JetDistance ret;
+  ret.j1 = -1;
+  ret.j2 = -1;
+  ret.dij = -1.;
+  bool dmin_set(false);
+  double dmin(0.);
+
+  while(!_distances.empty()){
+    JetDistance dist = _distances.top();
+    _distances.pop();
+    if(JetsUnmerged(dist)){
+      if(!dmin_set){
+	dmin = dist.dij;
+	dmin_set = true;
+      }
+      double weight = exp(-_rigidity * (dist.dij-dmin) /dmin);
+      popped_distances.push_back(make_pair(dist,weight));        
+      norm += weight; 
+      if(weight/norm < _truncation_fctr)
+	break;            
+    }
+  }
+
+  double rand(Rand()), tot_weight(0.);
+  for(vector<pair<JetDistance, double> >::iterator it = popped_distances.begin(); it != popped_distances.end(); it++){
+    tot_weight += (*it).second/norm;
+    if(tot_weight >= rand){
+      ret = (*it).first;
+      break;
+    }
+  }
+
+  // repopulate in reverse (maybe quicker?)
+  for(vector<pair<JetDistance, double> >::reverse_iterator it = popped_distances.rbegin(); it != popped_distances.rend(); it++)
+    if(JetsUnmerged((*it).first))
+      _distances.push((*it).first);
+
+  return ret;
+}
+
+void Qjets::Cluster(fastjet::ClusterSequence & cs){
+  computeDCut(cs);
+
+  // Populate all the distances
+  ComputeAllDistances(cs.jets());
+  JetDistance jd = GetNextDistance();
+
+  while(!_distances.empty() && jd.dij != -1.){
+    if(!Prune(jd,cs)){
+      //      _merged_jets.push_back(jd.j1);
+      //      _merged_jets.push_back(jd.j2);
+      _merged_jets[jd.j1] = true;
+      _merged_jets[jd.j2] = true;
+
+      int new_jet=-1;
+      cs.plugin_record_ij_recombination(jd.j1, jd.j2, 1., new_jet);
+      if(!JetUnmerged(new_jet))
+	edm::LogError("QJets Clustering") << "Problem with FastJet::plugin_record_ij_recombination";
+      ComputeNewDistanceMeasures(cs,new_jet);
+    } else {
+      double j1pt = cs.jets()[jd.j1].perp();
+      double j2pt = cs.jets()[jd.j2].perp();
+      if(j1pt>j2pt){
+	//	_merged_jets.push_back(jd.j2);
+	_merged_jets[jd.j2] = true;
+	cs.plugin_record_iB_recombination(jd.j2, 1.);
+      } else {
+	//	_merged_jets.push_back(jd.j1);
+	_merged_jets[jd.j1] = true;
+	cs.plugin_record_iB_recombination(jd.j1, 1.);
+      }
+    }
+    jd = GetNextDistance();
+  } 
+
+  // merge remaining jets with beam
+  int num_merged_final(0);
+  for(unsigned int i = 0 ; i < cs.jets().size(); i++)
+    if(JetUnmerged(i)){
+      num_merged_final++;
+      cs.plugin_record_iB_recombination(i,1.);
+    }
+
+  if(!(num_merged_final < 2))
+    edm::LogError("QJets Clustering") << "More than 1 jet remains after reclustering";
+}
+
+void Qjets::computeDCut(fastjet::ClusterSequence & cs){
+  // assume all jets in cs form a single jet.  compute mass and pt
+  fastjet::PseudoJet sum(0.,0.,0.,0.);
+  for(vector<fastjet::PseudoJet>::const_iterator it = cs.jets().begin(); it != cs.jets().end(); it++)
+    sum += (*it);
+  _dcut = 2. * _dcut_fctr * sum.m()/sum.perp(); 
+}
+
+bool Qjets::Prune(JetDistance& jd,fastjet::ClusterSequence & cs){
+  double pt1 = cs.jets()[jd.j1].perp();
+  double pt2 = cs.jets()[jd.j2].perp();
+  fastjet::PseudoJet sum_jet = cs.jets()[jd.j1]+cs.jets()[jd.j2];
+  double sum_pt = sum_jet.perp();
+  double z = min(pt1,pt2)/sum_pt;
+  double d2 = cs.jets()[jd.j1].plain_distance(cs.jets()[jd.j2]);
+  return (d2 > _dcut*_dcut) && (z < _zcut);
+}
+
+void Qjets::ComputeAllDistances(const vector<fastjet::PseudoJet>& inp){
+  for(unsigned int i = 0 ; i < inp.size()-1; i++){
+    // jet-jet distances
+    for(unsigned int j = i+1 ; j < inp.size(); j++){
+      JetDistance jd;
+      jd.j1 = i;
+      jd.j2 = j;
+      if(jd.j1 != jd.j2){
+	jd.dij = d_ij(inp[i],inp[j]);
+	_distances.push(jd);
+      }
+    }    
+  }
+}
+
+void Qjets::ComputeNewDistanceMeasures(fastjet::ClusterSequence & cs, unsigned int new_jet){
+  // jet-jet distances
+  for(unsigned int i = 0; i < cs.jets().size(); i++)
+    if(JetUnmerged(i) && i != new_jet){
+      JetDistance jd;
+      jd.j1 = new_jet;
+      jd.j2 = i;
+      jd.dij = d_ij(cs.jets()[jd.j1], cs.jets()[jd.j2]);
+      _distances.push(jd);
+    }
+}
+
+double Qjets::d_ij(const fastjet::PseudoJet& v1,const  fastjet::PseudoJet& v2) const{
+  double p1 = v1.perp();
+  double p2 = v2.perp();
+  double ret = pow(min(p1,p2),_exp_min) * pow(max(p1,p2),_exp_max) * v1.squared_distance(v2);
+  if(edm::isNotFinite(ret))
+    edm::LogError("QJets Clustering") << "d_ij is not finite";
+  return ret; 
+}
+
+double Qjets::Rand(){
+  double ret = 0.;
+  //if(_rand_seed_set)
+  //  ret = rand_r(&_seed)/(double)RAND_MAX;
+  //else 
+  //ret = rand()/(double)RAND_MAX;
+  ret = _rnEngine->flat();
+  return ret;
+}

RecoJets/JetAlgorithms/src/QjetsPlugin.cc

@@ -0,0 +1,24 @@
+#include "RecoJets/JetAlgorithms/interface/QjetsPlugin.h"
+
+using namespace std;
+
+void QjetsPlugin::SetRandSeed(unsigned int seed){
+  _rand_seed_set = true;
+  _seed = seed;
+}
+
+double QjetsPlugin::R()const{
+  return 0.;
+}
+
+string QjetsPlugin::description() const{
+  string desc("Qjets pruning plugin");
+  return desc;
+}
+
+void QjetsPlugin::run_clustering(fastjet::ClusterSequence & cs) const{
+  Qjets qjets(_zcut, _dcut_fctr, _exp_min, _exp_max, _rigidity, _truncation_fctr, _rnEngine);
+  if(_rand_seed_set)
+    qjets.SetRandSeed(_seed);
+  qjets.Cluster(cs);
+}

RecoJets/JetProducers/interface/NjettinessAdder.h

@@ -0,0 +1,35 @@
+#ifndef NjettinessAdder_h
+#define NjettinessAdder_h
+
+#include <memory>
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/EDProducer.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+#include "DataFormats/PatCandidates/interface/Jet.h"
+
+class NjettinessAdder : public edm::EDProducer { 
+ public:
+  explicit NjettinessAdder(const edm::ParameterSet& iConfig) :
+    src_(iConfig.getParameter<edm::InputTag>("src")),
+    src_token_(consumes<edm::View<reco::PFJet>>(src_)),
+    cone_(iConfig.getParameter<double>("cone"))
+      {
+	produces<edm::ValueMap<float> >("tau1");
+	produces<edm::ValueMap<float> >("tau2");
+	produces<edm::ValueMap<float> >("tau3");
+      }
+
+    virtual ~NjettinessAdder() {}
+
+    void produce(edm::Event & iEvent, const edm::EventSetup & iSetup) ;
+    float getTau(int num,edm::Ptr<reco::PFJet> object) const;
+
+ private:	
+    edm::InputTag src_ ;
+    edm::EDGetTokenT<edm::View<reco::PFJet>> src_token_;
+    double cone_ ;
+};
+
+#endif