New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Addition of code for nSubjettiness and QJets #2429
Changes from 4 commits
a637315
e3db711
968e9b9
b5f2ad2
68372e2
3770d1d
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,52 @@ | ||
#ifndef _QJETS_ | ||
#define _QJETS_ | ||
#include <queue> | ||
#include <vector> | ||
#include <list> | ||
#include <algorithm> | ||
#include "fastjet/JetDefinition.hh" | ||
#include "fastjet/PseudoJet.hh" | ||
#include "fastjet/ClusterSequence.hh" | ||
|
||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Overall, placing most of it in the reco namespace may be beneficial |
||
using namespace std; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. using directives are forbidden in the header files, please remove |
||
|
||
struct jet_distance{ | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Consider starting with a capital letter (CMS style for non-trivial types) |
||
double dij; | ||
int j1; | ||
int j2; | ||
}; | ||
|
||
class JetDistanceCompare{ | ||
public: | ||
JetDistanceCompare(){}; | ||
bool operator() (const jet_distance& lhs, const jet_distance&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; | ||
map<int,bool> _merged_jets; | ||
priority_queue <jet_distance, vector<jet_distance>, JetDistanceCompare> _distances; | ||
|
||
double d_ij(const fastjet::PseudoJet& v1, const fastjet::PseudoJet& v2) const; | ||
void ComputeDCut(fastjet::ClusterSequence & cs); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. methods start with lower case (CMS style) |
||
|
||
double Rand(); | ||
bool Prune(jet_distance& jd,fastjet::ClusterSequence & cs); | ||
bool JetsUnmerged(const jet_distance& jd) const; | ||
bool JetUnmerged(int num) const; | ||
void ComputeNewDistanceMeasures(fastjet::ClusterSequence & cs, unsigned int new_jet); | ||
void ComputeAllDistances(const vector<fastjet::PseudoJet>& inp); | ||
double ComputeMinimumDistance(); | ||
double ComputeNormalization(double dmin); | ||
jet_distance GetNextDistance(); | ||
bool Same(const jet_distance& lhs, const jet_distance& rhs); | ||
|
||
public: | ||
Qjets(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr); | ||
void Cluster(fastjet::ClusterSequence & cs); | ||
void SetRandSeed(unsigned int seed); /* In case you want reproducible behavior */ | ||
}; | ||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,21 @@ | ||
#ifndef _QJETSPLUGIN_ | ||
#define _QJETSPLUGIN_ | ||
#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; | ||
public: | ||
QjetsPlugin(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr = 0.); | ||
void SetRandSeed(unsigned int seed); /* In case you want reproducible behavior */ | ||
double R() const; | ||
string description() const; | ||
void run_clustering(fastjet::ClusterSequence & cs) const; | ||
}; | ||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,176 @@ | ||
#include "RecoJets/JetAlgorithms/interface/Qjets.h" | ||
|
||
Qjets::Qjets(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. this could've been inlined, but ok |
||
: _rand_seed_set(false), | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Defaults for RECO should have reproducible results (better be so in most of the code) There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. For reproducibility, the random seed is set in RecoJets/JetProducers/plugins/QjetsAdder.cc |
||
_zcut(zcut), | ||
_dcut(-1.), | ||
_dcut_fctr(dcut_fctr), | ||
_exp_min(exp_min), | ||
_exp_max(exp_max), | ||
_rigidity(rigidity), | ||
_truncation_fctr(truncation_fctr) | ||
{ | ||
} | ||
|
||
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 jet_distance& jd) const{ | ||
return JetUnmerged(jd.j1) && JetUnmerged(jd.j2); | ||
} | ||
|
||
jet_distance Qjets::GetNextDistance(){ | ||
vector< pair<jet_distance, double> > popped_distances; | ||
double norm(0.); | ||
jet_distance ret; | ||
ret.j1 = -1; | ||
ret.j2 = -1; | ||
ret.dij = -1.; | ||
bool dmin_set(false); | ||
double dmin(0.); | ||
|
||
while(!_distances.empty()){ | ||
jet_distance 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<jet_distance, 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<jet_distance, 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()); | ||
jet_distance 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; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. here and elsewhere, initialize a variable before passing it to another function |
||
cs.plugin_record_ij_recombination(jd.j1, jd.j2, 1., new_jet); | ||
assert(JetUnmerged(new_jet)); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Does it have to be an assert here (kill the job just because this condition should never happen)? |
||
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.); | ||
} | ||
|
||
assert(num_merged_final < 2); | ||
} | ||
|
||
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(jet_distance& 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 d = sqrt(cs.jets()[jd.j1].plain_distance(cs.jets()[jd.j2])); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. no need for sqrt here, just use comparison with _dcut**2 |
||
return (d > _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++){ | ||
jet_distance 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){ | ||
jet_distance 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); | ||
assert(!std::isnan(ret)); | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more.
|
||
return ret; | ||
} | ||
|
||
double Qjets::Rand(){ | ||
double ret = 0.; | ||
if(_rand_seed_set) | ||
ret = rand_r(&_seed)/(double)RAND_MAX; | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. don't use system random function, use CMS random engines |
||
else | ||
ret = rand()/(double)RAND_MAX; | ||
return ret; | ||
} |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,33 @@ | ||
#include "RecoJets/JetAlgorithms/interface/QjetsPlugin.h" | ||
|
||
QjetsPlugin::QjetsPlugin(double zcut, double dcut_fctr, double exp_min, double exp_max, double rigidity, double truncation_fctr) | ||
: _rand_seed_set(false), | ||
_zcut(zcut), | ||
_dcut_fctr(dcut_fctr), | ||
_exp_min(exp_min), | ||
_exp_max(exp_max), | ||
_rigidity(rigidity), | ||
_truncation_fctr(truncation_fctr) | ||
{ | ||
} | ||
|
||
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); | ||
if(_rand_seed_set) | ||
qjets.SetRandSeed(_seed); | ||
qjets.Cluster(cs); | ||
} |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,33 @@ | ||
#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")), | ||
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_ ; | ||
double cone_ ; | ||
}; | ||
|
||
#endif |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,46 @@ | ||
#ifndef QjetsAdder_h | ||
#define QjetsAdder_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" | ||
#include "RecoJets/JetAlgorithms/interface/QjetsPlugin.h" | ||
|
||
class QjetsAdder : public edm::EDProducer { | ||
public: | ||
explicit QjetsAdder(const edm::ParameterSet& iConfig) : | ||
src_(iConfig.getParameter<edm::InputTag>("src")), | ||
qjetsAlgo_( iConfig.getParameter<double>("zcut"), | ||
iConfig.getParameter<double>("dcutfctr"), | ||
iConfig.getParameter<double>("expmin"), | ||
iConfig.getParameter<double>("expmax"), | ||
iConfig.getParameter<double>("rigidity")), | ||
ntrial_(iConfig.getParameter<int>("ntrial")), | ||
cutoff_(iConfig.getParameter<double>("cutoff")), | ||
jetRad_(iConfig.getParameter<double>("jetRad")), | ||
mJetAlgo_(iConfig.getParameter<std::string>("jetAlgo")) , | ||
QjetsPreclustering_(iConfig.getParameter<int>("preclustering")) | ||
{ | ||
produces<edm::ValueMap<float> >("QjetsVolatility"); | ||
} | ||
|
||
virtual ~QjetsAdder() {} | ||
|
||
void produce(edm::Event & iEvent, const edm::EventSetup & iSetup) ; | ||
|
||
private: | ||
edm::InputTag src_ ; | ||
QjetsPlugin qjetsAlgo_ ; | ||
int ntrial_; | ||
double cutoff_; | ||
double jetRad_; | ||
std::string mJetAlgo_; | ||
int QjetsPreclustering_; | ||
}; | ||
|
||
|
||
#endif |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Please use more unique include guards