| 1 | //STARTHEADER
|
|---|
| 2 | // $Id: JetDefinition.cc,v 1.1 2008-11-06 14:32:15 ovyn Exp $
|
|---|
| 3 | //
|
|---|
| 4 | // Copyright (c) 2005-2007, Matteo Cacciari and Gavin Salam
|
|---|
| 5 | //
|
|---|
| 6 | //----------------------------------------------------------------------
|
|---|
| 7 | // This file is part of FastJet.
|
|---|
| 8 | //
|
|---|
| 9 | // FastJet is free software; you can redistribute it and/or modify
|
|---|
| 10 | // it under the terms of the GNU General Public License as published by
|
|---|
| 11 | // the Free Software Foundation; either version 2 of the License, or
|
|---|
| 12 | // (at your option) any later version.
|
|---|
| 13 | //
|
|---|
| 14 | // The algorithms that underlie FastJet have required considerable
|
|---|
| 15 | // development and are described in hep-ph/0512210. If you use
|
|---|
| 16 | // FastJet as part of work towards a scientific publication, please
|
|---|
| 17 | // include a citation to the FastJet paper.
|
|---|
| 18 | //
|
|---|
| 19 | // FastJet is distributed in the hope that it will be useful,
|
|---|
| 20 | // but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 21 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 22 | // GNU General Public License for more details.
|
|---|
| 23 | //
|
|---|
| 24 | // You should have received a copy of the GNU General Public License
|
|---|
| 25 | // along with FastJet; if not, write to the Free Software
|
|---|
| 26 | // Foundation, Inc.:
|
|---|
| 27 | // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
|
|---|
| 28 | //----------------------------------------------------------------------
|
|---|
| 29 | //ENDHEADER
|
|---|
| 30 |
|
|---|
| 31 | #include "../include/fastjet/JetDefinition.hh"
|
|---|
| 32 | #include "../include/fastjet/Error.hh"
|
|---|
| 33 | #include<sstream>
|
|---|
| 34 |
|
|---|
| 35 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 36 |
|
|---|
| 37 | using namespace std;
|
|---|
| 38 |
|
|---|
| 39 | string JetDefinition::description() const {
|
|---|
| 40 | ostringstream name;
|
|---|
| 41 | if (jet_algorithm() == plugin_algorithm) {
|
|---|
| 42 | return plugin()->description();
|
|---|
| 43 | } else if (jet_algorithm() == kt_algorithm) {
|
|---|
| 44 | name << "Longitudinally invariant kt algorithm with R = " << R();
|
|---|
| 45 | name << " and " << recombiner()->description();
|
|---|
| 46 | } else if (jet_algorithm() == cambridge_algorithm) {
|
|---|
| 47 | name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
|
|---|
| 48 | << R() ;
|
|---|
| 49 | name << " and " << recombiner()->description();
|
|---|
| 50 | } else if (jet_algorithm() == antikt_algorithm) {
|
|---|
| 51 | name << "Longitudinally invariant anti-kt algorithm with R = "
|
|---|
| 52 | << R() ;
|
|---|
| 53 | name << " and " << recombiner()->description();
|
|---|
| 54 | } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
|
|---|
| 55 | name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
|
|---|
| 56 | << R() << "and a special hack whereby particles with kt < "
|
|---|
| 57 | << extra_param() << "are treated as passive ghosts";
|
|---|
| 58 | } else {
|
|---|
| 59 | throw Error("Unrecognized jet_finder");
|
|---|
| 60 | }
|
|---|
| 61 | return name.str();
|
|---|
| 62 | }
|
|---|
| 63 |
|
|---|
| 64 |
|
|---|
| 65 | void JetDefinition::set_recombination_scheme(
|
|---|
| 66 | RecombinationScheme recomb_scheme) {
|
|---|
| 67 | _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
|
|---|
| 68 | _recombiner = 0;
|
|---|
| 69 | }
|
|---|
| 70 |
|
|---|
| 71 |
|
|---|
| 72 | string JetDefinition::DefaultRecombiner::description() const {
|
|---|
| 73 | switch(_recomb_scheme) {
|
|---|
| 74 | case E_scheme:
|
|---|
| 75 | return "E scheme recombination";
|
|---|
| 76 | case pt_scheme:
|
|---|
| 77 | return "pt scheme recombination";
|
|---|
| 78 | case pt2_scheme:
|
|---|
| 79 | return "pt2 scheme recombination";
|
|---|
| 80 | case Et_scheme:
|
|---|
| 81 | return "Et scheme recombination";
|
|---|
| 82 | case Et2_scheme:
|
|---|
| 83 | return "Et2 scheme recombination";
|
|---|
| 84 | case BIpt_scheme:
|
|---|
| 85 | return "boost-invariant pt scheme recombination";
|
|---|
| 86 | case BIpt2_scheme:
|
|---|
| 87 | return "boost-invariant pt2 scheme recombination";
|
|---|
| 88 | default:
|
|---|
| 89 | ostringstream err;
|
|---|
| 90 | err << "DefaultRecombiner: unrecognized recombination scheme "
|
|---|
| 91 | << _recomb_scheme;
|
|---|
| 92 | throw Error(err.str());
|
|---|
| 93 | }
|
|---|
| 94 | }
|
|---|
| 95 |
|
|---|
| 96 |
|
|---|
| 97 | void JetDefinition::DefaultRecombiner::recombine(
|
|---|
| 98 | const PseudoJet & pa, const PseudoJet & pb,
|
|---|
| 99 | PseudoJet & pab) const {
|
|---|
| 100 |
|
|---|
| 101 | double weighta, weightb;
|
|---|
| 102 |
|
|---|
| 103 | switch(_recomb_scheme) {
|
|---|
| 104 | case E_scheme:
|
|---|
| 105 | pab = pa + pb;
|
|---|
| 106 | pab.set_user_index(0);
|
|---|
| 107 | return;
|
|---|
| 108 | // all remaining schemes are massless recombinations and locally
|
|---|
| 109 | // we just set weights, while the hard work is done below...
|
|---|
| 110 | case pt_scheme:
|
|---|
| 111 | case Et_scheme:
|
|---|
| 112 | case BIpt_scheme:
|
|---|
| 113 | weighta = pa.perp();
|
|---|
| 114 | weightb = pb.perp();
|
|---|
| 115 | break;
|
|---|
| 116 | case pt2_scheme:
|
|---|
| 117 | case Et2_scheme:
|
|---|
| 118 | case BIpt2_scheme:
|
|---|
| 119 | weighta = pa.perp2();
|
|---|
| 120 | weightb = pb.perp2();
|
|---|
| 121 | break;
|
|---|
| 122 | default:
|
|---|
| 123 | ostringstream err;
|
|---|
| 124 | err << "DefaultRecombiner: unrecognized recombination scheme "
|
|---|
| 125 | << _recomb_scheme;
|
|---|
| 126 | throw Error(err.str());
|
|---|
| 127 | }
|
|---|
| 128 |
|
|---|
| 129 | double perp_ab = pa.perp() + pb.perp();
|
|---|
| 130 | if (perp_ab != 0.0) { // weights also non-zero...
|
|---|
| 131 | double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
|
|---|
| 132 |
|
|---|
| 133 | // take care with periodicity in phi...
|
|---|
| 134 | double phi_a = pa.phi(), phi_b = pb.phi();
|
|---|
| 135 | if (phi_a - phi_b > pi) phi_b += twopi;
|
|---|
| 136 | if (phi_a - phi_b < -pi) phi_b -= twopi;
|
|---|
| 137 | double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
|
|---|
| 138 |
|
|---|
| 139 | pab = PseudoJet(perp_ab*cos(phi_ab),
|
|---|
| 140 | perp_ab*sin(phi_ab),
|
|---|
| 141 | perp_ab*sinh(y_ab),
|
|---|
| 142 | perp_ab*cosh(y_ab));
|
|---|
| 143 | } else { // weights are zero
|
|---|
| 144 | pab = PseudoJet(0.0,0.0,0.0,0.0);
|
|---|
| 145 | }
|
|---|
| 146 | pab.set_user_index(0);
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 |
|
|---|
| 150 | void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const {
|
|---|
| 151 | switch(_recomb_scheme) {
|
|---|
| 152 | case E_scheme:
|
|---|
| 153 | case BIpt_scheme:
|
|---|
| 154 | case BIpt2_scheme:
|
|---|
| 155 | break;
|
|---|
| 156 | case pt_scheme:
|
|---|
| 157 | case pt2_scheme:
|
|---|
| 158 | {
|
|---|
| 159 | // these schemes (as in the ktjet implementation) need massless
|
|---|
| 160 | // initial 4-vectors with essentially E=|p|.
|
|---|
| 161 | double newE = sqrt(p.perp2()+p.pz()*p.pz());
|
|---|
| 162 | int user_index = p.user_index();
|
|---|
| 163 | p = PseudoJet(p.px(), p.py(), p.pz(), newE);
|
|---|
| 164 | p.set_user_index(user_index);
|
|---|
| 165 | }
|
|---|
| 166 | break;
|
|---|
| 167 | case Et_scheme:
|
|---|
| 168 | case Et2_scheme:
|
|---|
| 169 | {
|
|---|
| 170 | // these schemes (as in the ktjet implementation) need massless
|
|---|
| 171 | // initial 4-vectors with essentially E=|p|.
|
|---|
| 172 | double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
|
|---|
| 173 | int user_index = p.user_index();
|
|---|
| 174 | p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
|
|---|
| 175 | p.set_user_index(user_index);
|
|---|
| 176 | }
|
|---|
| 177 | break;
|
|---|
| 178 | default:
|
|---|
| 179 | ostringstream err;
|
|---|
| 180 | err << "DefaultRecombiner: unrecognized recombination scheme "
|
|---|
| 181 | << _recomb_scheme;
|
|---|
| 182 | throw Error(err.str());
|
|---|
| 183 | }
|
|---|
| 184 | }
|
|---|
| 185 |
|
|---|
| 186 | void JetDefinition::Plugin::set_ghost_separation_scale(double scale) const {
|
|---|
| 187 | throw Error("set_ghost_separation_scale not supported");
|
|---|
| 188 | }
|
|---|
| 189 |
|
|---|
| 190 |
|
|---|
| 191 | FASTJET_END_NAMESPACE
|
|---|
