Page MenuHomeHEPForge
Files F8309690

ModifiedMassDropTagger.cc
No OneTemporary
Actions

ModifiedMassDropTagger.cc
View Options

// $Id: ModifiedMassDropTagger.cc 520 2014-01-24 10:58:19Z gsalam $
//
// Copyright (c) 2014-, Gavin P. Salam
//
//----------------------------------------------------------------------
// This file is part of FastJet contrib.
//
// It is free software; you can redistribute it and/or modify it under
// the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at
// your option) any later version.
//
// It is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
// License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this code. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
#include "ModifiedMassDropTagger.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/ClusterSequenceAreaBase.hh"
#include <algorithm>
#include <cstdlib>
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
namespace contrib{
LimitedWarning ModifiedMassDropTagger::_negative_mass_warning;
LimitedWarning ModifiedMassDropTagger::_mu2_gt1_warning;
LimitedWarning ModifiedMassDropTagger::_nonca_warning;
LimitedWarning ModifiedMassDropTagger::_explicit_ghost_warning;
bool ModifiedMassDropTagger::_verbose = false;
//----------------------------------------------------------------------
PseudoJet ModifiedMassDropTagger::result(const PseudoJet & j) const {
// issue a warning if the jet is not obtained through a C/A
// clustering
if ((! j.has_associated_cluster_sequence()) ||
(j.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm))
_nonca_warning.warn("ModifiedMassDropTagger is designed to be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
if (_subtractor) {
const ClusterSequenceAreaBase * csab =
dynamic_cast<const ClusterSequenceAreaBase *>(j.associated_cs());
if (csab == 0 || (!csab->has_explicit_ghosts()))
_explicit_ghost_warning.warn("ModifiedMassDropTagger: there is no clustering sequence, or it lacks explicit ghosts: subtraction is not guaranteed to function properly");
}
// establish the first subjet and optionally subtract it
PseudoJet subjet = j;
if (_subtractor && (!_input_jet_is_subtracted)) {
subjet = (*_subtractor)(subjet);
}
bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
// variables for tracking what will happen
PseudoJet piece1, piece2;
// now recurse into the jet's structure
while (subjet.has_parents(piece1, piece2)) {
// first sanity check:
// - zero or negative pts are not allowed for the input subjet
// - zero or negative masses are not allowed for configurations
// in which the mass will effectively appear in a denominator
// (The masses will be checked later)
if (subjet.pt2() <= 0) return PseudoJet();
if (_subtractor) {
piece1 = (*_subtractor)(piece1);
piece2 = (*_subtractor)(piece2);
}
// determine the symmetry parameter
double sym;
if (_symmetry_measure == y) {
// the original d_{ij}/m^2 choice from MDT
// first make sure the mass is sensible
if (subjet.m2() <= 0) {
_negative_mass_warning.warn("ModifiedMassDropTagger: cannot calculate y, because (sub)jet mass is negative; bailing out");
return PseudoJet();
}
sym = piece1.kt_distance(piece2) / subjet.m2();
} else if (_symmetry_measure == vector_z) {
// min(pt1, pt2)/(pt), where the denominator is a vector sum
// of the two subjets
sym = min(piece1.pt(), piece2.pt()) / subjet.pt();
} else if (_symmetry_measure == scalar_z) {
// min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum
// of the two subjets
double pt1 = piece1.pt();
double pt2 = piece2.pt();
// make sure denominator is non-zero
sym = pt1 + pt2;
if (sym == 0) return PseudoJet();
sym = min(pt1, pt2) / sym;
} else {
throw Error ("Unrecognized choice of symmetry_measure");
}
// determine the symmetry cut (by default this returns a constant)
double this_symmetry_cut = symmetry_cut_fn(piece1, piece2);
// and make a first tagging decision based on symmetry cut
bool tagged = (sym > this_symmetry_cut);
// if tagged based on symmetry cut, then check the mu cut (if relevant)
// and update the tagging decision. Calculate mu^2 regardless, for cases
// of users not cutting on mu2, but still interested in its value.
double mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
if (tagged && use_mu_cut) {
// first a sanity check -- mu2 won't be sensible if the subjet mass
// is negative, so we can't then trust the mu cut - bail out
if (subjet.m2() <= 0) {
_negative_mass_warning.warn("ModifiedMassDropTagger: cannot trust mu, because (sub)jet mass is negative; bailing out");
return PseudoJet();
}
if (mu2 > 1) _mu2_gt1_warning.warn("ModifiedMassDropTagger encountered mu^2 value > 1");
if (mu2 > pow(_mu_cut,2)) tagged = false;
}
// if we've tagged the splitting, return the jet with its substructure
if (tagged) {
// record relevant information
StructureType * structure = new StructureType(subjet);
structure->_symmetry = sym;
structure->_mu = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
structure->_delta_R = piece1.delta_R(piece2);
subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
return subjet;
}
// otherwise continue unclustering, allowing for the different
// ways of choosing which parent to look into
int choice;
if (_recursion_choice == larger_mt) {
choice = piece1.mt2() > piece2.mt2() ? 1 : 2;
} else if (_recursion_choice == larger_pt) {
choice = piece1.pt2() > piece2.pt2() ? 1 : 2;
} else if (_recursion_choice == larger_m) {
choice = piece1.m2() > piece2.m2() ? 1 : 2;
} else {
throw Error ("Unrecognized value for recursion_choice");
}
if (_verbose) cout << "choice is " << choice << endl;;
subjet = (choice == 1) ? piece1 : piece2;
}
if (_verbose) cout << "reached end; returning null jet " << endl;
return PseudoJet();
}
//----------------------------------------------------------------------
string ModifiedMassDropTagger::symmetry_cut_description() const {
ostringstream ostr;
ostr << _symmetry_cut;
return ostr.str();
}
//----------------------------------------------------------------------
string ModifiedMassDropTagger::description() const {
ostringstream ostr;
ostr << "Recursive Tagger with a symmetry cut ";
switch(_symmetry_measure) {
case y:
ostr << "y"; break;
case scalar_z:
ostr << "scalar_z"; break;
case vector_z:
ostr << "vector_z"; break;
default:
cerr << "failed to interpret symmetry_measure" << endl; exit(-1);
}
ostr << " > " << symmetry_cut_description();
if (_mu_cut != numeric_limits<double>::infinity()) {
ostr << ", mass-drop cut mu=max(m1,m2)/m < " << _mu_cut;
} else {
ostr << ", no mass-drop requirement";
}
ostr << ", recursion into the subjet with larger ";
switch(_recursion_choice) {
case larger_pt:
ostr << "pt"; break;
case larger_mt:
ostr << "mt(=sqrt(m^2+pt^2))"; break;
case larger_m:
ostr << "mass"; break;
default:
cerr << "failed to interpret recursion_choice" << endl; exit(-1);
}
if (_subtractor) {
ostr << " and subtractor: " << _subtractor->description();
if (_input_jet_is_subtracted) {ostr << " (input jet is assumed already subtracted)";}
}
return ostr.str();
}
} // namespace contrib
FASTJET_END_NAMESPACE

File Metadata

Mime Type
text/x-c
Expires
Sat, Dec 21, 4:12 PM (21 h, 47 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023402
Default Alt Text
ModifiedMassDropTagger.cc (7 KB)

Event Timeline