| 1 | // Nsubjettiness Package
|
|---|
| 2 | // Questions/Comments? jthaler@jthaler.net
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2011-14
|
|---|
| 5 | // Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
|
|---|
| 6 | //
|
|---|
| 7 | // $Id: ExtraRecombiners.cc 842 2015-08-20 13:44:31Z jthaler $
|
|---|
| 8 | //----------------------------------------------------------------------
|
|---|
| 9 | // This file is part of FastJet contrib.
|
|---|
| 10 | //
|
|---|
| 11 | // It is free software; you can redistribute it and/or modify it under
|
|---|
| 12 | // the terms of the GNU General Public License as published by the
|
|---|
| 13 | // Free Software Foundation; either version 2 of the License, or (at
|
|---|
| 14 | // your option) any later version.
|
|---|
| 15 | //
|
|---|
| 16 | // It is distributed in the hope that it will be useful, but WITHOUT
|
|---|
| 17 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
|
|---|
| 18 | // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
|
|---|
| 19 | // License for more details.
|
|---|
| 20 | //
|
|---|
| 21 | // You should have received a copy of the GNU General Public License
|
|---|
| 22 | // along with this code. If not, see <http://www.gnu.org/licenses/>.
|
|---|
| 23 | //----------------------------------------------------------------------
|
|---|
| 24 |
|
|---|
| 25 | #include "ExtraRecombiners.hh"
|
|---|
| 26 |
|
|---|
| 27 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 28 |
|
|---|
| 29 | namespace contrib{
|
|---|
| 30 |
|
|---|
| 31 | std::string GeneralEtSchemeRecombiner::description() const {
|
|---|
| 32 | return "General Et-scheme recombination";
|
|---|
| 33 | }
|
|---|
| 34 |
|
|---|
| 35 | // recombine pa and pb according to a generalized Et-scheme parameterized by the power delta
|
|---|
| 36 | void GeneralEtSchemeRecombiner::recombine(const fastjet::PseudoJet & pa, const fastjet::PseudoJet & pb, fastjet::PseudoJet & pab) const {
|
|---|
| 37 |
|
|---|
| 38 | // Define new weights for recombination according to delta
|
|---|
| 39 | // definition of ratio done so that we do not encounter issues about numbers being too large for huge values of delta
|
|---|
| 40 | double ratio;
|
|---|
| 41 | if (std::abs(_delta - 1.0) < std::numeric_limits<double>::epsilon()) ratio = pb.perp()/pa.perp(); // save computation time of pow()
|
|---|
| 42 | else ratio = pow(pb.perp()/pa.perp(), _delta);
|
|---|
| 43 | double weighta = 1.0/(1.0 + ratio);
|
|---|
| 44 | double weightb = 1.0/(1.0 + 1.0/ratio);
|
|---|
| 45 |
|
|---|
| 46 | double perp_ab = pa.perp() + pb.perp();
|
|---|
| 47 | // reweight the phi and rap sums according to the weights above
|
|---|
| 48 | if (perp_ab != 0.0) {
|
|---|
| 49 | double y_ab = (weighta * pa.rap() + weightb * pb.rap());
|
|---|
| 50 |
|
|---|
| 51 | double phi_a = pa.phi(), phi_b = pb.phi();
|
|---|
| 52 | if (phi_a - phi_b > pi) phi_b += twopi;
|
|---|
| 53 | if (phi_a - phi_b < -pi) phi_b -= twopi;
|
|---|
| 54 | double phi_ab = (weighta * phi_a + weightb * phi_b);
|
|---|
| 55 |
|
|---|
| 56 | pab.reset_PtYPhiM(perp_ab, y_ab, phi_ab);
|
|---|
| 57 |
|
|---|
| 58 | }
|
|---|
| 59 | else {
|
|---|
| 60 | pab.reset(0.0,0.0,0.0,0.0);
|
|---|
| 61 | }
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 |
|
|---|
| 65 | std::string WinnerTakeAllRecombiner::description() const {
|
|---|
| 66 | return "Winner-Take-All recombination";
|
|---|
| 67 | }
|
|---|
| 68 |
|
|---|
| 69 | // recombine pa and pb by creating pab with energy of the sum of particle energies in the direction of the harder particle
|
|---|
| 70 | // updated recombiner to use more general form of a metric equal to E*(pT/E)^(alpha), which reduces to pT*cosh(rap)^(1-alpha)
|
|---|
| 71 | // alpha is specified by the user. The default is alpha = 1, which is the typical behavior. alpha = 2 provides a metric which more
|
|---|
| 72 | // favors central jets
|
|---|
| 73 | void WinnerTakeAllRecombiner::recombine(const fastjet::PseudoJet & pa, const fastjet::PseudoJet & pb, fastjet::PseudoJet & pab) const {
|
|---|
| 74 | double a_pt = pa.perp(), b_pt = pb.perp(), a_rap = pa.rap(), b_rap = pb.rap();
|
|---|
| 75 |
|
|---|
| 76 | // special case of alpha = 1, everything is just pt (made separate so that pow function isn't called)
|
|---|
| 77 | if (_alpha == 1.0) {
|
|---|
| 78 | if (a_pt >= b_pt) {
|
|---|
| 79 | pab.reset_PtYPhiM(a_pt + b_pt, a_rap, pa.phi());
|
|---|
| 80 | }
|
|---|
| 81 | else if (b_pt > a_pt) {
|
|---|
| 82 | pab.reset_PtYPhiM(a_pt + b_pt, b_rap, pb.phi());
|
|---|
| 83 | }
|
|---|
| 84 | }
|
|---|
| 85 |
|
|---|
| 86 | // every other case uses additional cosh(rap) term
|
|---|
| 87 | else {
|
|---|
| 88 | double a_metric = a_pt*pow(cosh(a_rap), 1.0-_alpha);
|
|---|
| 89 | double b_metric = b_pt*pow(cosh(b_rap), 1.0-_alpha);
|
|---|
| 90 | if (a_metric >= b_metric) {
|
|---|
| 91 | double new_pt = a_pt + b_pt*pow(cosh(b_rap)/cosh(a_rap), 1.0-_alpha);
|
|---|
| 92 | pab.reset_PtYPhiM(new_pt, a_rap, pa.phi());
|
|---|
| 93 | }
|
|---|
| 94 | if (b_metric > a_metric) {
|
|---|
| 95 | double new_pt = b_pt + a_pt*pow(cosh(a_rap)/cosh(b_rap), 1.0-_alpha);
|
|---|
| 96 | pab.reset_PtYPhiM(new_pt, b_rap, pb.phi());
|
|---|
| 97 | }
|
|---|
| 98 | }
|
|---|
| 99 | }
|
|---|
| 100 |
|
|---|
| 101 | } //namespace contrib
|
|---|
| 102 |
|
|---|
| 103 | FASTJET_END_NAMESPACE
|
|---|
