JetAlgorithms/src/Qjets.cc

0001 #include "RecoJets/JetAlgorithms/interface/Qjets.h"
0002
0003 using namespace std;
0004
0005 void Qjets::SetRandSeed(unsigned int seed) {
0006   _rand_seed_set = true;
0007   _seed = seed;
0008 }
0009
0010 bool Qjets::JetUnmerged(int num) const { return _merged_jets.find(num) == _merged_jets.end(); }
0011
0012 bool Qjets::JetsUnmerged(const JetDistance& jd) const { return JetUnmerged(jd.j1) && JetUnmerged(jd.j2); }
0013
0014 JetDistance Qjets::GetNextDistance() {
0015   vector<pair<JetDistance, double> > popped_distances;
0016   double norm(0.);
0017   JetDistance ret;
0018   ret.j1 = -1;
0019   ret.j2 = -1;
0020   ret.dij = -1.;
0021   bool dmin_set(false);
0022   double dmin(0.);
0023
0024   while (!_distances.empty()) {
0025     JetDistance dist = _distances.top();
0026     _distances.pop();
0027     if (JetsUnmerged(dist)) {
0028       if (!dmin_set) {
0029         dmin = dist.dij;
0030         dmin_set = true;
0031       }
0032       double weight = exp(-_rigidity * (dist.dij - dmin) / dmin);
0033       popped_distances.push_back(make_pair(dist, weight));
0034       norm += weight;
0035       if (weight / norm < _truncation_fctr)
0036         break;
0037     }
0038   }
0039
0040   double rand(Rand()), tot_weight(0.);
0041   for (vector<pair<JetDistance, double> >::iterator it = popped_distances.begin(); it != popped_distances.end(); it++) {
0042     tot_weight += (*it).second / norm;
0043     if (tot_weight >= rand) {
0044       ret = (*it).first;
0045       omega *= ((*it).second);
0046       break;
0047     }
0048   }
0049
0050   // repopulate in reverse (maybe quicker?)
0051   for (vector<pair<JetDistance, double> >::reverse_iterator it = popped_distances.rbegin();
0052        it != popped_distances.rend();
0053        it++)
0054     if (JetsUnmerged((*it).first))
0055       _distances.push((*it).first);
0056
0057   return ret;
0058 }
0059
0060 void Qjets::Cluster(fastjet::ClusterSequence& cs) {
0061   omega = 1.;
0062   QjetsBaseExtras* extras = new QjetsBaseExtras();
0063   computeDCut(cs);
0064
0065   // Populate all the distances
0066   ComputeAllDistances(cs.jets());
0067   JetDistance jd = GetNextDistance();
0068
0069   while (!_distances.empty() && jd.dij != -1.) {
0070     if (!Prune(jd, cs)) {
0071       //      _merged_jets.push_back(jd.j1);
0072       //      _merged_jets.push_back(jd.j2);
0073       _merged_jets[jd.j1] = true;
0074       _merged_jets[jd.j2] = true;
0075
0076       int new_jet = -1;
0077       cs.plugin_record_ij_recombination(jd.j1, jd.j2, 1., new_jet);
0078       if (!JetUnmerged(new_jet))
0079         edm::LogError("QJets Clustering") << "Problem with FastJet::plugin_record_ij_recombination";
0080       ComputeNewDistanceMeasures(cs, new_jet);
0081     } else {
0082       double j1pt = cs.jets()[jd.j1].perp();
0083       double j2pt = cs.jets()[jd.j2].perp();
0084       if (j1pt > j2pt) {
0085         //  _merged_jets.push_back(jd.j2);
0086         _merged_jets[jd.j2] = true;
0087         cs.plugin_record_iB_recombination(jd.j2, 1.);
0088       } else {
0089         //  _merged_jets.push_back(jd.j1);
0090         _merged_jets[jd.j1] = true;
0091         cs.plugin_record_iB_recombination(jd.j1, 1.);
0092       }
0093     }
0094     jd = GetNextDistance();
0095   }
0096
0097   extras->_wij = omega;
0098   cs.plugin_associate_extras(extras);
0099
0100   // merge remaining jets with beam
0101   int num_merged_final(0);
0102   for (unsigned int i = 0; i < cs.jets().size(); i++)
0103     if (JetUnmerged(i)) {
0104       num_merged_final++;
0105       cs.plugin_record_iB_recombination(i, 1.);
0106     }
0107
0108   if (!(num_merged_final < 2))
0109     edm::LogError("QJets Clustering") << "More than 1 jet remains after reclustering";
0110 }
0111
0112 void Qjets::computeDCut(fastjet::ClusterSequence& cs) {
0113   // assume all jets in cs form a single jet.  compute mass and pt
0114   fastjet::PseudoJet sum(0., 0., 0., 0.);
0115   for (vector<fastjet::PseudoJet>::const_iterator it = cs.jets().begin(); it != cs.jets().end(); it++)
0116     sum += (*it);
0117   _dcut = 2. * _dcut_fctr * sum.m() / sum.perp();
0118 }
0119
0120 bool Qjets::Prune(JetDistance& jd, fastjet::ClusterSequence& cs) {
0121   double pt1 = cs.jets()[jd.j1].perp();
0122   double pt2 = cs.jets()[jd.j2].perp();
0123   fastjet::PseudoJet sum_jet = cs.jets()[jd.j1] + cs.jets()[jd.j2];
0124   double sum_pt = sum_jet.perp();
0125   double z = min(pt1, pt2) / sum_pt;
0126   double d2 = cs.jets()[jd.j1].plain_distance(cs.jets()[jd.j2]);
0127   return (d2 > _dcut * _dcut) && (z < _zcut);
0128 }
0129
0130 void Qjets::ComputeAllDistances(const vector<fastjet::PseudoJet>& inp) {
0131   for (unsigned int i = 0; i < inp.size() - 1; i++) {
0132     // jet-jet distances
0133     for (unsigned int j = i + 1; j < inp.size(); j++) {
0134       JetDistance jd;
0135       jd.j1 = i;
0136       jd.j2 = j;
0137       if (jd.j1 != jd.j2) {
0138         jd.dij = d_ij(inp[i], inp[j]);
0139         _distances.push(jd);
0140       }
0141     }
0142   }
0143 }
0144
0145 void Qjets::ComputeNewDistanceMeasures(fastjet::ClusterSequence& cs, unsigned int new_jet) {
0146   // jet-jet distances
0147   for (unsigned int i = 0; i < cs.jets().size(); i++)
0148     if (JetUnmerged(i) && i != new_jet) {
0149       JetDistance jd;
0150       jd.j1 = new_jet;
0151       jd.j2 = i;
0152       jd.dij = d_ij(cs.jets()[jd.j1], cs.jets()[jd.j2]);
0153       _distances.push(jd);
0154     }
0155 }
0156
0157 double Qjets::d_ij(const fastjet::PseudoJet& v1, const fastjet::PseudoJet& v2) const {
0158   double p1 = v1.perp();
0159   double p2 = v2.perp();
0160   double ret = pow(min(p1, p2), _exp_min) * pow(max(p1, p2), _exp_max) * v1.squared_distance(v2);
0161   if (edm::isNotFinite(ret))
0162     edm::LogError("QJets Clustering") << "d_ij is not finite";
0163   return ret;
0164 }
0165
0166 double Qjets::Rand() {
0167   double ret = 0.;
0168   //if(_rand_seed_set)
0169   //  ret = rand_r(&_seed)/(double)RAND_MAX;
0170   //else
0171   //ret = rand()/(double)RAND_MAX;
0172   ret = _rnEngine->flat();
0173   return ret;
0174 }