27#include "split_merge.h"
28#include "siscone_error.h"
56 pass = CJET_INEXISTENT_PASS;
69 return j1.v.perp2() >
j2.v.perp2();
107#ifdef EPSILON_SPLITMERGE
156std::string split_merge_scale_name(Esplit_merge_scale
sms) {
159 return "pt (IR unsafe)";
161 return "Et (boost dep.)";
163 return "mt (IR safe except for pairs of identical decayed heavy particles)";
165 return "pttilde (scalar sum of pt's)";
167 return "[SM scale without a name]";
190 if (
j1.contents[
i1]==
j2.contents[
i2]) {
193 }
else if (
j1.contents[
i1]<
j2.contents[
i2]){
194 (*v) += (*particles)[
j1.contents[
i1]];
195 (*pt_tilde) += (*pt)[
j1.contents[
i1]];
197 }
else if (
j1.contents[
i1]>
j2.contents[
i2]){
198 (*v) -= (*particles)[
j2.contents[
i2]];
199 (*pt_tilde) -= (*pt)[
j2.contents[
i2]];
202 throw Csiscone_error(
"get_non_overlap reached part it should never have seen...");
208 (*v) += (*particles)[
j1.contents[
i1]];
209 (*pt_tilde) += (*pt)[
j1.contents[
i1]];
213 (*v) -= (*particles)[
j2.contents[
i2]];
214 (*pt_tilde) -= (*pt)[
j2.contents[
i2]];
228#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
229#ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
248 use_pt_weighted_splitting =
false;
285 for (
int i=0;i<
n;i++)
350 epr.eta_min = eta_min-0.01;
351 epr.eta_max = eta_max+0.01;
376#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
436 if (dphi>M_PI) dphi =
twopi-dphi;
469 if (protocones->size()==0)
472 pt_min2 = ptmin*ptmin;
489 jet.contents.clear();
506 jet.n=
jet.contents.size();
517#ifdef DEBUG_SPLIT_MERGE
518 cout <<
"adding jet: ";
531#ifdef DEBUG_SPLIT_MERGE
532 cout <<
"remaining particles: ";
543#ifdef DEBUG_SPLIT_MERGE
549#ifdef DEBUG_SPLIT_MERGE
582 if (protocones->size()==0)
585 pt_min2 = ptmin*ptmin;
660 jets[
jets.size()-1].v.build_etaphi();
662#ifdef DEBUG_SPLIT_MERGE
663 cout <<
"PR-Jet " <<
jets.size() <<
" [size " <<
jet.contents.size() <<
"]:";
670 for (
int index=0;index<
n_left;index++){
676#ifdef DEBUG_SPLIT_MERGE
695#ifdef DEBUG_SPLIT_MERGE
721 pt_min2 = ptmin*ptmin;
723 if (candidates->size()==0)
728 message <<
"Illegal value for overlap_tshold, f = " <<
overlap_tshold;
729 message <<
" (legal values are 0<f<1)";
742 if (candidates->size()>0){
744 j1 = candidates->begin();
755 while (
j2 != candidates->end()){
756#ifdef DEBUG_SPLIT_MERGE
763#ifdef DEBUG_SPLIT_MERGE
764 cout <<
"overlap between cdt 1 and cdt " <<
j2_relindex+1 <<
" with overlap "
772 j2 =
j1 = candidates->begin();
779 j2 =
j1 = candidates->begin();
788 if (
j2 != candidates->end())
j2++;
791 if (
j1 != candidates->end()) {
796 jets[
jets.size()-1].v.build_etaphi();
801#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
804 candidates->erase(
j1);
813 }
while (candidates->size()>0);
817#ifdef DEBUG_SPLIT_MERGE
835 fprintf(
flux,
"# columns are: eta, phi, pt and number of particles for each jet\n");
838 j1->v.build_etaphi();
840 j1->v.eta,
j1->v.phi,
j1->v.perp(),
j1->n);
844 fprintf(
flux,
"# columns are: eta, phi, pt, particle index and jet number\n");
875 for (
it_c = candidates->begin(),
i1=0 ;
it_c != candidates->end() ;
it_c++,
i1++){
878 c->v.px,
c->v.py,
c->v.pz,
c->v.E,
sqrt(
c->sm_var2));
897 if (!is_range_overlap(
j1.range,
j2.range))
912 if (
j1.contents[
i1]<
j2.contents[
i2]){
927 }
while ((i1<j1.
n) && (i2<j2.
n));
945 (*overlap2) = get_sm_var2(v, pt_tilde);
962bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
965 double eta1, phi1, pt1_weight, eta2, phi2, pt2_weight;
966 double dx1, dy1, dx2, dy2;
971 const Cjet & j1 = * it_j1;
972 const Cjet & j2 = * it_j2;
975 jet2.
v = jet1.v = Cmomentum();
976 jet2.pt_tilde = jet1.pt_tilde = 0.0;
987 pt1_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
993 pt2_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
995 jet1.v = jet2.v = Cmomentum();
999 if (j1.contents[i1]<j2.contents[i2]){
1002 jet1.contents.push_back(j1.contents[i1]);
1004 jet1.pt_tilde +=
pt[j1.contents[i1]];
1006 jet1.range.add_particle(v->eta,v->phi);
1007 }
else if (j1.contents[i1]>j2.contents[i2]){
1010 jet2.contents.push_back(j2.contents[i2]);
1012 jet2.pt_tilde +=
pt[j2.contents[i2]];
1014 jet2.range.add_particle(v->eta,v->phi);
1020 dx1 = eta1 - v->eta;
1021 dy1 = fabs(phi1 - v->phi);
1026 dx2 = eta2 - v->eta;
1027 dy2 = fabs(phi2 - v->phi);
1035 double d1sq = (dx1*dx1+dy1*dy1)*pt1_weight;
1036 double d2sq = (dx2*dx2+dy2*dy2)*pt2_weight;
1043 jet1.contents.push_back(j1.contents[i1]);
1045 jet1.pt_tilde +=
pt[j1.contents[i1]];
1046 jet1.range.add_particle(v->eta,v->phi);
1049 jet2.contents.push_back(j2.contents[i2]);
1051 jet2.pt_tilde +=
pt[j2.contents[i2]];
1052 jet2.range.add_particle(v->eta,v->phi);
1058 }
while ((i1<j1.n) && (i2<j2.n));
1062 jet1.contents.push_back(j1.contents[i1]);
1064 jet1.pt_tilde +=
pt[j1.contents[i1]];
1066 jet1.range.add_particle(v->eta,v->phi);
1070 jet2.contents.push_back(j2.contents[i2]);
1072 jet2.pt_tilde +=
pt[j2.contents[i2]];
1074 jet2.range.add_particle(v->eta,v->phi);
1078 jet1.n = jet1.contents.size();
1079 jet2.n = jet2.contents.size();
1085#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1086 cand_refs.erase(j1.v.ref);
1087 cand_refs.erase(j2.v.ref);
1089 candidates->erase(it_j1);
1090 candidates->erase(it_j2);
1106bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
1113 jet.contents.push_back(
indices[i]);
1117 jet.n = jet.contents.size();
1120 jet.range = range_union(it_j1->range, it_j2->range);
1123#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1125 cand_refs.erase(it_j1->v.ref);
1126 cand_refs.erase(it_j2->v.ref);
1129 candidates->erase(it_j1);
1130 candidates->erase(it_j2);
1143bool Csplit_merge::insert(Cjet &jet){
1148#ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1154 if (jet.v.perp2()<pt_min2)
1158 jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
1161 candidates->insert(jet);
1172double Csplit_merge::get_sm_var2(Cmomentum &v,
double &pt_tilde){
1174 case SM_pt:
return v.perp2();
1175 case SM_mt:
return v.perpmass2();
1176 case SM_pttilde:
return pt_tilde*pt_tilde;
1177 case SM_Et:
return v.Et2();
1179 throw Csiscone_error(
"Unsupported split-merge scale choice: "
class for holding a covering range in eta-phi
double pt_tilde
p-scheme pt
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
int n
number of particles inside
int pass
pass at which the jet has been found It starts at 0 (first pass), -1 means infinite rapidity (it will...
std::vector< int > contents
particle contents (list of indices)
base class for dynamic coordinates management
int index
internal particle number
double eta
particle pseudo-rapidity
int parent_index
particle number in the parent list
double phi
particle azimuthal angle
class corresponding to errors that will be thrown by siscone
std::vector< double > * pt
pointer to the pt of the particles
bool operator()(const Cjet &jet1, const Cjet &jet2) const
comparison between 2 jets
void get_difference(const Cjet &j1, const Cjet &j2, Cmomentum *v, double *pt_tilde) const
get the difference between 2 jets, calculated such that rounding errors will not affect the result ev...
std::vector< Cmomentum > * particles
pointer to the list of particles
std::string SM_scale_name() const
return the name corresponding to the SM scale variable
Esplit_merge_scale split_merge_scale
the following parameter controls the variable we're using for the split-merge process i....
int init(std::vector< Cmomentum > &_particles, std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
initialisation function
int n_pass
index of the run
int idx_size
number of elements in indices1
double SM_var2_hardest_cut_off
stop split–merge or progressive-removal when the squared SM_var of the hardest protojet is below this...
int n_left
numer of particles that does not belong to any jet
int * indices
maximal size array for indices works
int add_hardest_protocone_to_jets(std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
remove the hardest protocone and declare it a jet
bool merge_identical_protocones
The following flag indicates that identical protocones are to be merged automatically each time aroun...
std::vector< double > pt
list of particles' pt
int partial_clear()
partial clearance
int show()
show jets/candidates status
int init_pleft()
build initial list of left particles
Csplit_merge_ptcomparison ptcomparison
member used for detailed comparisons of pt's
int merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles and removing particles ...
int perform(double overlap_tshold, double ptmin=0.0)
really do the splitting and merging At the end, the vector jets is filled with the jets found.
std::vector< Cmomentum > particles
list of particles
double most_ambiguous_split
minimal difference in squared distance between a particle and two overlapping protojets when doing a ...
int add_protocones(std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
add a list of protocones
~Csplit_merge()
default dtor
int save_contents(FILE *flux)
save final jets
double stable_cone_soft_pt2_cutoff
pt cutoff for the particles to put in p_uncol_hard this is meant to allow removing soft particles in ...
std::vector< Cmomentum > p_uncol_hard
list of particles remaining with collinear clustering
std::vector< Cjet > jets
list of jets
Csplit_merge()
default ctor
int full_clear()
full clearance
std::vector< Cmomentum > p_remain
list of particles remaining to deal with
int init_particles(std::vector< Cmomentum > &_particles)
initialisation function for particle list
a circulator that is foreseen to take as template member either a pointer or an iterator;
#define EPSILON_SPLITMERGE
The following define enables you to allow for identical protocones to be merged automatically after e...
#define EPSILON_COLLINEAR
The following parameter controls collinear safety.
const double twopi
definition of 2*M_PI which is useful a bit everyhere!