JetDefinition.cc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/JetDefinition.hh"
00032 #include "fastjet/Error.hh"
00033 #include<sstream>
00034
00035 FASTJET_BEGIN_NAMESPACE
00036
00037 using namespace std;
00038
00039
00040
00041
00042 JetDefinition::JetDefinition(JetAlgorithm jet_algorithm,
00043 double R,
00044 Strategy strategy,
00045 RecombinationScheme recomb_scheme,
00046 int nparameters) :
00047 _jet_algorithm(jet_algorithm), _Rparam(R), _strategy(strategy) {
00048
00049
00050 if (jet_algorithm == ee_kt_algorithm) {
00051 _Rparam = 4.0;
00052
00053
00054
00055 } else if (jet_algorithm != ee_genkt_algorithm) {
00056 assert(_Rparam <= 0.5*pi);
00057 }
00058
00059
00060
00061 ostringstream oss;
00062 switch (jet_algorithm) {
00063 case ee_kt_algorithm:
00064 if (nparameters != 0) oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with "
00065 << nparameters << " parameter(s)\n";
00066 break;
00067 case genkt_algorithm:
00068 case ee_genkt_algorithm:
00069 if (nparameters != 2) oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with "
00070 << nparameters << " parameter(s)\n";
00071 break;
00072 default:
00073 if (nparameters != 1)
00074 oss << "The jet algorithm you requested ("
00075 << jet_algorithm << ") should be constructed with 1 parameter but was called with "
00076 << nparameters << " parameter(s)\n";
00077 }
00078 if (oss.str() != "") throw Error(oss.str());
00079
00080
00081 assert (_strategy != plugin_strategy);
00082
00083 _plugin = NULL;
00084 set_recombination_scheme(recomb_scheme);
00085 set_extra_param(0.0);
00086 }
00087
00088
00089
00090 string JetDefinition::description() const {
00091 ostringstream name;
00092 if (jet_algorithm() == plugin_algorithm) {
00093 return plugin()->description();
00094 } else if (jet_algorithm() == kt_algorithm) {
00095 name << "Longitudinally invariant kt algorithm with R = " << R();
00096 name << " and " << recombiner()->description();
00097 } else if (jet_algorithm() == cambridge_algorithm) {
00098 name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
00099 << R() ;
00100 name << " and " << recombiner()->description();
00101 } else if (jet_algorithm() == antikt_algorithm) {
00102 name << "Longitudinally invariant anti-kt algorithm with R = "
00103 << R() ;
00104 name << " and " << recombiner()->description();
00105 } else if (jet_algorithm() == genkt_algorithm) {
00106 name << "Longitudinally invariant generalised kt algorithm with R = "
00107 << R() << ", p = " << extra_param();
00108 name << " and " << recombiner()->description();
00109 } else if (jet_algorithm() == cambridge_for_passive_algorithm) {
00110 name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
00111 << R() << "and a special hack whereby particles with kt < "
00112 << extra_param() << "are treated as passive ghosts";
00113 } else if (jet_algorithm() == ee_kt_algorithm) {
00114 name << "e+e- kt (Durham) algorithm (NB: no R)";
00115 name << " with " << recombiner()->description();
00116 } else if (jet_algorithm() == ee_genkt_algorithm) {
00117 name << "e+e- generalised kt algorithm with R = "
00118 << R() << ", p = " << extra_param();
00119 name << " and " << recombiner()->description();
00120 } else {
00121 throw Error("JetDefinition::description(): unrecognized jet_finder");
00122 }
00123 return name.str();
00124 }
00125
00126
00127 void JetDefinition::set_recombination_scheme(
00128 RecombinationScheme recomb_scheme) {
00129 _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
00130 _recombiner = 0;
00131 }
00132
00133
00134 string JetDefinition::DefaultRecombiner::description() const {
00135 switch(_recomb_scheme) {
00136 case E_scheme:
00137 return "E scheme recombination";
00138 case pt_scheme:
00139 return "pt scheme recombination";
00140 case pt2_scheme:
00141 return "pt2 scheme recombination";
00142 case Et_scheme:
00143 return "Et scheme recombination";
00144 case Et2_scheme:
00145 return "Et2 scheme recombination";
00146 case BIpt_scheme:
00147 return "boost-invariant pt scheme recombination";
00148 case BIpt2_scheme:
00149 return "boost-invariant pt2 scheme recombination";
00150 default:
00151 ostringstream err;
00152 err << "DefaultRecombiner: unrecognized recombination scheme "
00153 << _recomb_scheme;
00154 throw Error(err.str());
00155 }
00156 }
00157
00158
00159 void JetDefinition::DefaultRecombiner::recombine(
00160 const PseudoJet & pa, const PseudoJet & pb,
00161 PseudoJet & pab) const {
00162
00163 double weighta, weightb;
00164
00165 switch(_recomb_scheme) {
00166 case E_scheme:
00167 pab = pa + pb;
00168 pab.set_user_index(0);
00169 return;
00170
00171
00172 case pt_scheme:
00173 case Et_scheme:
00174 case BIpt_scheme:
00175 weighta = pa.perp();
00176 weightb = pb.perp();
00177 break;
00178 case pt2_scheme:
00179 case Et2_scheme:
00180 case BIpt2_scheme:
00181 weighta = pa.perp2();
00182 weightb = pb.perp2();
00183 break;
00184 default:
00185 ostringstream err;
00186 err << "DefaultRecombiner: unrecognized recombination scheme "
00187 << _recomb_scheme;
00188 throw Error(err.str());
00189 }
00190
00191 double perp_ab = pa.perp() + pb.perp();
00192 if (perp_ab != 0.0) {
00193 double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
00194
00195
00196 double phi_a = pa.phi(), phi_b = pb.phi();
00197 if (phi_a - phi_b > pi) phi_b += twopi;
00198 if (phi_a - phi_b < -pi) phi_b -= twopi;
00199 double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
00200
00201 pab = PseudoJet(perp_ab*cos(phi_ab),
00202 perp_ab*sin(phi_ab),
00203 perp_ab*sinh(y_ab),
00204 perp_ab*cosh(y_ab));
00205 } else {
00206 pab = PseudoJet(0.0,0.0,0.0,0.0);
00207 }
00208 pab.set_user_index(0);
00209 }
00210
00211
00212 void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const {
00213 switch(_recomb_scheme) {
00214 case E_scheme:
00215 case BIpt_scheme:
00216 case BIpt2_scheme:
00217 break;
00218 case pt_scheme:
00219 case pt2_scheme:
00220 {
00221
00222
00223 double newE = sqrt(p.perp2()+p.pz()*p.pz());
00224 int user_index = p.user_index();
00225 p = PseudoJet(p.px(), p.py(), p.pz(), newE);
00226 p.set_user_index(user_index);
00227 }
00228 break;
00229 case Et_scheme:
00230 case Et2_scheme:
00231 {
00232
00233
00234 double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
00235 int user_index = p.user_index();
00236 p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
00237 p.set_user_index(user_index);
00238 }
00239 break;
00240 default:
00241 ostringstream err;
00242 err << "DefaultRecombiner: unrecognized recombination scheme "
00243 << _recomb_scheme;
00244 throw Error(err.str());
00245 }
00246 }
00247
00248 void JetDefinition::Plugin::set_ghost_separation_scale(double scale) const {
00249 throw Error("set_ghost_separation_scale not supported");
00250 }
00251
00252
00253 FASTJET_END_NAMESPACE