Index: contrib/contribs/LundPlane/trunk/LundEEGenerator.cc =================================================================== --- contrib/contribs/LundPlane/trunk/LundEEGenerator.cc (revision 1289) +++ contrib/contribs/LundPlane/trunk/LundEEGenerator.cc (revision 1290) @@ -1,222 +0,0 @@ -// $Id: LundEEGenerator.cc 1153 2018-08-22 11:56:35Z frdreyer $ -// -// Copyright (c) 2018-, Frederic A. Dreyer, A. Karlberg, Gavin P. Salam, -// Ludovic Scyboz, Gregory Soyez -// -//---------------------------------------------------------------------- -// 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 . -//---------------------------------------------------------------------- - -#include -#include "LundEEGenerator.hh" - -using namespace std; - -FASTJET_BEGIN_NAMESPACE - -namespace contrib{ - -LundEEDeclustering::LundEEDeclustering(const PseudoJet& pair, - const PseudoJet& j1, const PseudoJet& j2, - int iplane, double psi, double psibar, - int depth, int leaf_iplane, int sign_s) - : iplane_(iplane), psi_(psi), psibar_(psibar), m_(pair.m()), pair_(pair), depth_(depth), leaf_iplane_(leaf_iplane), sign_s_(sign_s) { - - double omc = one_minus_costheta(j1,j2); - - if (omc > sqrt(numeric_limits::epsilon())) { - double cos_theta = 1.0 - omc; - double theta = acos(cos_theta); - sin_theta_ = sin(theta); - eta_ = -log(tan(theta/2.0)); - } else { - // we are at small angles, so use small-angle formulas - double theta = sqrt(2. * omc); - sin_theta_ = theta; - eta_ = -log(theta/2); - } - - // establish which of j1 and j2 is softer - if (j1.modp2() > j2.modp2()) { - harder_ = j1; - softer_ = j2; - } else { - harder_ = j2; - softer_ = j1; - } - - // now work out the various LundEE declustering variables - double softer_modp = softer_.modp(); - z_ = softer_modp / (softer_modp + harder_.modp()); - kt_ = softer_modp * sin_theta_; - lnkt_ = log(kt_); - kappa_ = z_ * sin_theta_; -} - -//---------------------------------------------------------------------- -/// return a vector of the declusterings of the primary plane of the -/// jet, registering information about the plane from which this -/// comes -std::vector LundEEGenerator::result_iplane(const PseudoJet& jet, int iplane) const { - std::vector result; - PseudoJet j = jet; - PseudoJet pair, j1, j2, cross; - pair = j; - while (pair.has_parents(j1, j2)) { - // order the jets in their 3-mom norm - if (j1.modp2() < j2.modp2()) std::swap(j1,j2); - // on the first invocation cross has its default value of 0; - result.push_back(single_declustering(pair, j1, j2, cross, iplane)); - pair = j1; - } - return result; -} - -//---------------------------------------------------------------------- -/// get declustering corresponding to a given pseudojet and its parents. -/// -/// If cross is a null vector, then the psi value for the declustering -/// will be set to zero and the cross vector will be defined as -/// j1 x j2. -/// -/// If cross is non-zero, then the psi value is determined as the -/// angle between this j1 x j2 and the existing cross. -LundEEDeclustering LundEEGenerator::single_declustering(const PseudoJet& pair, - const PseudoJet& j1, - const PseudoJet& j2, - PseudoJet& cross, - int iplane) const { - // get the angle between current declustering plane and previous one - PseudoJet j1xj2 = cross_product(j1,j2,false); - // the first time one encounters the next line cross == 0 - // (default PJ initialisation), so one simply defines the direction - // meant by psi=0 - double psi; - if (cross == double(0)) { - psi = 0; - cross = j1xj2; - // normalise it - cross /= cross.modp(); - } else { - psi = theta(cross, j1xj2); - } - - // prepare for next iteration - // return declustering element - return LundEEDeclustering(pair, j1, j2, iplane, psi); -} - -//---------------------------------------------------------------------- -/// create an ordered declustering sequence for one jet -std::vector LundEEGeneratorOrdered::result_iplane(const PseudoJet& jet, - int iplane) const { - std::priority_queue< std::pair > pq_declusts; - this->fill_pq(jet, iplane, pq_declusts); - return this->pq_to_vector(pq_declusts); -} - -//---------------------------------------------------------------------- -/// create an ordered declustering sequence for two jets -std::vector LundEEGeneratorOrdered::result_twojets(const PseudoJet& jet1, - const PseudoJet& jet2) const { - std::priority_queue< std::pair > pq_declusts; - this->fill_pq(jet1, 0, pq_declusts); - this->fill_pq(jet2, 1, pq_declusts); - return this->pq_to_vector(pq_declusts); -} - -//---------------------------------------------------------------------- -/// create an ordered declustering sequence for multiple jets, along with n secondary planes -std::vector LundEEGeneratorOrdered::result_jets(const std::vector& jets, - int n) const{ - std::priority_queue< std::pair > pq_declusts; - unsigned int i=0; - for(; i < jets.size(); ++i) - this->fill_pq(jets[i], i, pq_declusts); - std::vector< std::pair > tmp; - // get the n elements from which we want secondary planes - for(int isec=0; isec < n; ++isec) { - tmp.push_back(pq_declusts.top()); - pq_declusts.pop(); - } - // get the n secondary declustering sequences and insert them into priority queue - for(int isec=0; isec < n; ++isec) { - this->fill_pq(tmp[isec].second.softer(), isec+i, pq_declusts); - pq_declusts.push(tmp[isec]); - } - // return sorted vector of all the declusterings - return this->pq_to_vector(pq_declusts); -} - -//---------------------------------------------------------------------- -/// fill a priority queue with ordered declusterings -void LundEEGeneratorOrdered::fill_pq(const PseudoJet& jet, int iplane, - std::priority_queue< std::pair >& pq) const { - - PseudoJet j = jet; - PseudoJet pair, j1, j2, cross; - pair = j; - while (pair.has_parents(j1, j2)) { - // order the jets in their 3-mom norm - if (j1.modp2() < j2.modp2()) std::swap(j1,j2); - LundEEDeclustering d = single_declustering(pair, j1, j2, cross, iplane); - pq.push(std::make_pair(this->ordering_measure(d), d)); - pair = j1; - } -} - -//---------------------------------------------------------------------- -/// turn a priority queue to an ordered vector -std::vector LundEEGeneratorOrdered::pq_to_vector(std::priority_queue< std::pair >& pq) const { - std::vector result; - PseudoJet cross; - while (!pq.empty()) { - LundEEDeclustering d = pq.top().second; - // update the psi angle according to the new ordering - PseudoJet newcross = PseudoJet(d.harder().py()*d.softer().pz() - - d.harder().pz()*d.softer().py(), - d.harder().pz()*d.softer().px() - - d.harder().px()*d.softer().pz(), - d.harder().px()*d.softer().py() - - d.harder().py()*d.softer().px(), 0.0); - if (cross != double(0)) d.set_psi(theta(cross, newcross)); - cross = newcross; - // add declustering to result vector - result.push_back(d); - pq.pop(); - } - return result; -} - -} // namespace contrib - - -FASTJET_END_NAMESPACE - -std::ostream & operator<<(std::ostream & ostr, const fastjet::contrib::LundEEDeclustering & d) { - ostr << "kt = " << d.kt() - << " z = " << d.z() - << " eta = " << d.eta() - << " psi = " << d.psi() - << " psibar = " << d.psibar() - << " m = " << d.m() - << " iplane = " << d.iplane() - << " depth = " << d.depth() - << " leaf_iplane = " << d.leaf_iplane(); - return ostr; -} - Index: contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.cc =================================================================== --- contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.cc (revision 0) +++ contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.cc (revision 1290) @@ -0,0 +1,86 @@ +// $Id$ +// +// Copyright (c) 2018-, Frederic A. Dreyer, Keith Hamilton, Alexander Karlberg, +// Gavin P. Salam, Ludovic Scyboz, Gregory Soyez, Rob Verheyen +// +//---------------------------------------------------------------------- +// 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 . +//---------------------------------------------------------------------- + +#include +#include "RecursiveLundEEGenerator.hh" + +using namespace std; + +FASTJET_BEGIN_NAMESPACE + +namespace contrib{ + +LundEEDeclustering::LundEEDeclustering(const PseudoJet& pair, + const PseudoJet& j1, const PseudoJet& j2, + int iplane, double psi, double psibar, + int depth, int leaf_iplane, int sign_s) + : iplane_(iplane), psi_(psi), psibar_(psibar), m_(pair.m()), pair_(pair), depth_(depth), leaf_iplane_(leaf_iplane), sign_s_(sign_s) { + + double omc = one_minus_costheta(j1,j2); + + if (omc > sqrt(numeric_limits::epsilon())) { + double cos_theta = 1.0 - omc; + double theta = acos(cos_theta); + sin_theta_ = sin(theta); + eta_ = -log(tan(theta/2.0)); + } else { + // we are at small angles, so use small-angle formulas + double theta = sqrt(2. * omc); + sin_theta_ = theta; + eta_ = -log(theta/2); + } + + // establish which of j1 and j2 is softer + if (j1.modp2() > j2.modp2()) { + harder_ = j1; + softer_ = j2; + } else { + harder_ = j2; + softer_ = j1; + } + + // now work out the various LundEE declustering variables + double softer_modp = softer_.modp(); + z_ = softer_modp / (softer_modp + harder_.modp()); + kt_ = softer_modp * sin_theta_; + lnkt_ = log(kt_); + kappa_ = z_ * sin_theta_; +} + +} // namespace contrib + + +FASTJET_END_NAMESPACE + +std::ostream & operator<<(std::ostream & ostr, const fastjet::contrib::LundEEDeclustering & d) { + ostr << "kt = " << d.kt() + << " z = " << d.z() + << " eta = " << d.eta() + << " psi = " << d.psi() + << " psibar = " << d.psibar() + << " m = " << d.m() + << " iplane = " << d.iplane() + << " depth = " << d.depth() + << " leaf_iplane = " << d.leaf_iplane(); + return ostr; +} + Property changes on: contrib/contribs/LundPlane/trunk/RecursiveLundEEGenerator.cc ___________________________________________________________________ Added: svn:keywords ## -0,0 +1 ## +Id \ No newline at end of property