Fork me on GitHub

source: git/external/fastjet/contribs/Nsubjettiness/ExtraRecombiners.cc@ 7e6c201

Last change on this file since 7e6c201 was 973b92a, checked in by Pavel Demin <pavel.demin@…>, 9 years ago

update FastJet library to 3.1.3 and Nsubjettiness library to 2.2.1

  • Property mode set to 100644
File size: 4.0 KB
RevLine 
[9687203]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//
[973b92a]7// $Id: ExtraRecombiners.cc 842 2015-08-20 13:44:31Z jthaler $
[9687203]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
[973b92a]25#include "ExtraRecombiners.hh"
[9687203]26
27FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
28
29namespace contrib{
[973b92a]30
31std::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
36void 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
[9687203]64
65std::string WinnerTakeAllRecombiner::description() const {
[973b92a]66 return "Winner-Take-All recombination";
[9687203]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
73void 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
103FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.