fastjet_boosted_higgs.cc File Reference

#include "fastjet/ClusterSequence.hh"
#include <iostream>
#include <sstream>
#include <iomanip>
#include <cmath>

Include dependency graph for fastjet_boosted_higgs.cc:

Go to the source code of this file.

Classes
class	FlavourRecombiner
Typedefs
typedef fj::JetDefinition::DefaultRecombiner	DefRecomb
	set up a class to give standard (by default E-scheme) recombination, with additional tracking of flavour information in the user_index.
Functions
ostream &	operator<< (ostream &, fj::PseudoJet &)
	forward declaration for printing out info about a jet
int	main (int argc, char **argv)

---

Typedef Documentation

typedef fj::JetDefinition::DefaultRecombiner DefRecomb

set up a class to give standard (by default E-scheme) recombination, with additional tracking of flavour information in the user_index.

If you use this, you must explicitly set the user index to 0 for non-flavoured particles (the default value is -1);

This will work for native algorithms, but not for all plugins

Definition at line 73 of file fastjet_boosted_higgs.cc.

---

Function Documentation

int main	(	int	argc,
		char **	argv
	)

now do the subjet decomposition;

when unpeeling a C/A jet, often only a very soft piece may break off; the mass_drop_threshold indicates how much "lighter" the heavier of the two resulting pieces must be in order for us to consider that we've really seen some form of substructure

QCD backgrounds that give larger jet masses have a component where a quite soft gluon is emitted; to eliminate part of this one can place a cut on the asymmetry of the branching;

Here the cut is expressed in terms of y, the kt-distance scaled to the squared jet mass; an easier way to see it is in terms of a requirement on the momentum fraction in the splitting: z/(1-z) and (1-z)/z > rtycut^2 [the correspondence holds only at LO]

Definition at line 98 of file fastjet_boosted_higgs.cc.

References cambridge_algorithm, fastjet::d0::inline_maths::min(), and sorted_by_pt().

                                  {
  
  vector<fj::PseudoJet> particles;

  // read in data in format px py pz E b-tag [last of these is optional]
  string line;
  while (getline(cin,line)) {
    if (line.substr(0,1) == "#") {continue;}
    istringstream linestream(line);
    double px,py,pz,E;
    linestream >> px >> py >> pz >> E;

    // optionally read in btag information
    int    btag;
    if (! (linestream >> btag)) btag = 0;

    // construct the particle
    fj::PseudoJet particle(px,py,pz,E);
    particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
    particles.push_back(particle);
  }


  // set up the jet finding
  double R = 1.2;
  FlavourRecombiner flav_recombiner; // for tracking flavour
  fj::JetDefinition jet_def(fj::cambridge_algorithm, R, &flav_recombiner);

  
  // run the jet finding; find the hardest jet
  fj::ClusterSequence cs(particles, jet_def);
  vector<fj::PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());


  cout << "Ran: " << jet_def.description() << endl << endl;
  cout << "Hardest jet: " << jets[0] << endl << endl;


  //
  double mass_drop_threshold = 0.667; 
  double rtycut              = 0.3;

  fj::PseudoJet this_jet = jets[0], parent1, parent2;
  bool had_parents;

  while ((had_parents = cs.has_parents(this_jet,parent1,parent2))) {
    // make parent1 the more massive jet
    if (parent1.m() < parent2.m()) swap(parent1,parent2);
    //
    // if we pass the conditions on the mass drop and its degree of
    // asymmetry (z/(1-z) \sim kt_dist/m^2 > rtycut), then we've found
    // something interesting, so exit the loop
    if (parent1.m() < mass_drop_threshold * this_jet.m() &&
        parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
      break;
    } else {
      // otherwise try a futher decomposition on the more massive jet
      this_jet = parent1;
    }
  }

  // look to see what we found
  if (had_parents) {
    cout << "Found suitable pair of subjets: " << endl;
    cout << " " << parent1 << endl;
    cout << " " << parent2 << endl;
    cout << "Total = " << endl;
    cout << " " << this_jet << endl << endl;

    // next we "filter" it, to remove UE & pileup contamination
    //
    // [there are two ways of doing this; here we directly use the
    // exsiting cluster sequence and find the exclusive subjets of
    // this_jet (i.e. work backwards within the cs starting from
    // this_jet); alternatively one can recluster just the
    // constituents of the jet]
    //
    // first get separation between the subjets (called Rbb -- assuming it's a Higgs!)
    double   Rbb = sqrt(parent1.squared_distance(parent2));
    double   Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
    unsigned nfilt = 3;               // number of pieces we'll take
    cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;

    double   dcut  = pow(Rfilt/R,2);  // for C/A get a view at Rfilt by
                                    // using a dcut=(Rfilt/R)^2
    vector<fj::PseudoJet> filt_subjets = sorted_by_pt(cs.exclusive_subjets(this_jet, dcut));

    // now print out the filtered jets and reconstruct total 
    // at the same time
    cout << "Filtered pieces are " << endl;
    cout << " " << filt_subjets[0] << endl;
    fj::PseudoJet filtered_total = filt_subjets[0];
    for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
      cout << " " << filt_subjets[i] << endl;
      flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
    }
    cout << "Filtered total is " << endl;
    cout << " " << filtered_total << endl;

  } else {
    cout << "Did not find suitable hard substructure in this event." << endl;
  }
}

ostream & operator<<	(	ostream &	ostr,
		fj::PseudoJet &	jet
	)

forward declaration for printing out info about a jet

does the actual work for printing out a jet

Definition at line 217 of file fastjet_boosted_higgs.cc.

                                                        {
  ostr << "pt, y, phi =" 
       << " " << setw(10) << jet.perp() 
       << " " << setw(6) <<  jet.rap()  
       << " " << setw(6) <<  jet.phi()  
       << ", mass = " << setw(10) << jet.m()
       << ", btag = " << jet.user_index();
  return ostr;
}