27 #include "split_merge.h"
28 #include "siscone_error.h"
56 pass = CJET_INEXISTENT_PASS;
68 bool jets_pt_less(
const Cjet &j1,
const Cjet &j2){
107 #ifdef EPSILON_SPLITMERGE
112 double pt_tilde_difference;
113 get_difference(jet1,jet2,&difference,&pt_tilde_difference);
122 switch (split_merge_scale){
124 qdiff = sum.
E*difference.
E - sum.
pz*difference.
pz;
127 qdiff = sum.
px*difference.
px + sum.
py*difference.
py;
130 qdiff = pt_tilde_sum*pt_tilde_difference;
137 qdiff = jet1.
v.
E*jet1.
v.
E*
138 ((sum.
px*difference.
px + sum.
py*difference.
py)*jet1.
v.
pz*jet1.
v.
pz
156 std::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]";
194 (*v) += (*particles)[j1.
contents[i1]];
195 (*pt_tilde) += (*pt)[j1.
contents[i1]];
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...");
204 }
while ((i1<j1.
n) && (i2<j2.
n));
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]];
227 merge_identical_protocones =
false;
228 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
229 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
230 merge_identical_protocones =
true;
237 ptcomparison.particles = &particles;
238 ptcomparison.pt = &pt;
239 candidates.reset(
new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
242 SM_var2_hardest_cut_off = -numeric_limits<double>::max();
245 stable_cone_soft_pt2_cutoff = -1.0;
248 use_pt_weighted_splitting =
false;
265 int Csplit_merge::init(vector<Cmomentum> & , vector<Cmomentum> *protocones,
double R2,
double ptmin){
267 return add_protocones(protocones, R2, ptmin);
280 particles = _particles;
281 n = particles.size();
285 for (
int i=0;i<n;i++)
286 pt[i] = particles[i].perp();
289 ptcomparison.particles = &particles;
290 ptcomparison.pt = &pt;
295 indices =
new int[n];
321 particles[i].ref.randomize();
324 if (fabs(particles[i].pz) < (particles[i].E)){
325 p_remain.push_back(particles[i]);
327 p_remain[j].parent_index = i;
334 p_remain[j].index = 1;
338 particles[i].index = 0;
340 eta_min = min(eta_min, particles[i].eta);
341 eta_max = max(eta_max, particles[i].eta);
343 particles[i].index = -1;
346 n_left = p_remain.size();
353 merge_collinear_and_remove_soft();
370 candidates.reset(
new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
373 most_ambiguous_split = numeric_limits<double>::max();
376 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
377 if (merge_identical_protocones)
393 if (indices != NULL){
408 vector<Cmomentum> p_sorted;
412 p_uncol_hard.clear();
415 for (i=0;i<n_left;i++)
416 p_sorted.push_back(p_remain[i]);
417 sort(p_sorted.begin(), p_sorted.end(), momentum_eta_less);
426 if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) {
434 while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<
EPSILON_COLLINEAR) && (!collinear)){
435 dphi = fabs(p_sorted[j].phi-p_sorted[i].phi);
436 if (dphi>M_PI) dphi =
twopi-dphi;
439 p_sorted[j] += p_sorted[i];
447 p_uncol_hard.push_back(p_sorted[i]);
469 if (protocones->size()==0)
472 pt_min2 = ptmin*ptmin;
477 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
490 for (i=0;i<n_left;i++){
496 dy = fabs(phi - v->
phi);
517 #ifdef DEBUG_SPLIT_MERGE
518 cout <<
"adding jet: ";
519 for (
int i2=0;i2<jet.
n;i2++)
531 #ifdef DEBUG_SPLIT_MERGE
532 cout <<
"remaining particles: ";
535 for (i=0;i<n_left;i++){
536 if (p_remain[i].index){
538 p_remain[j]=p_remain[i];
539 p_remain[j].parent_index = p_remain[i].parent_index;
542 particles[p_remain[j].parent_index].index = n_pass;
543 #ifdef DEBUG_SPLIT_MERGE
544 cout << p_remain[j].parent_index <<
" ";
549 #ifdef DEBUG_SPLIT_MERGE
555 merge_collinear_and_remove_soft();
579 Cjet jet, jet_candidate;
580 bool found_jet =
false;
582 if (protocones->size()==0)
585 pt_min2 = ptmin*ptmin;
590 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
605 for (i=0;i<n_left;i++){
611 dy = fabs(phi - v->
phi);
616 jet_candidate.
v+= *v;
621 jet_candidate.
n=jet_candidate.
contents.size();
625 *c = jet_candidate.
v;
633 if (jet_candidate.
v.
perp2()<pt_min2)
640 jet_candidate.
sm_var2 = (*_user_scale)(jet_candidate);
643 jet_candidate.
sm_var2 = get_sm_var2(jet_candidate.
v, jet_candidate.
pt_tilde);
648 (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
649 : ptcomparison(jet_candidate, jet))){
656 if (!found_jet)
return 1;
662 #ifdef DEBUG_SPLIT_MERGE
663 cout <<
"PR-Jet " << jets.size() <<
" [size " << jet.
contents.size() <<
"]:";
667 int p_remain_index = 0;
668 int contents_index = 0;
670 for (
int index=0;index<n_left;index++){
671 if ((contents_index<(
int) jet.
contents.size()) &&
672 (p_remain[index].parent_index == jet.
contents[contents_index])){
675 particles[p_remain[index].parent_index].index = n_pass;
676 #ifdef DEBUG_SPLIT_MERGE
677 cout <<
" " << jet.
contents[contents_index];
682 p_remain[p_remain_index] = p_remain[index];
683 p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
684 p_remain[p_remain_index].index=1;
688 p_remain.resize(n_left-jet.
contents.size());
689 n_left = p_remain.size();
690 jets[jets.size()-1].pass = particles[jet.
contents[0]].index;
695 #ifdef DEBUG_SPLIT_MERGE
702 merge_collinear_and_remove_soft();
721 pt_min2 = ptmin*ptmin;
723 if (candidates->size()==0)
726 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
727 ostringstream message;
728 message <<
"Illegal value for overlap_tshold, f = " << overlap_tshold;
729 message <<
" (legal values are 0<f<1)";
739 double overlap_tshold2 = overlap_tshold*overlap_tshold;
742 if (candidates->size()>0){
744 j1 = candidates->begin();
748 if (j1->sm_var2<SM_var2_hardest_cut_off) {
break;}
755 while (j2 != candidates->end()){
756 #ifdef DEBUG_SPLIT_MERGE
760 if (get_overlap(*j1, *j2, &overlap2)){
763 #ifdef DEBUG_SPLIT_MERGE
764 cout <<
"overlap between cdt 1 and cdt " << j2_relindex+1 <<
" with overlap "
765 << sqrt(overlap2/j2->sm_var2) << endl<<endl;
767 if (overlap2<overlap_tshold2*j2->sm_var2){
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();
799 assert(j1->contents.size() > 0);
800 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
801 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
802 cand_refs.erase(j1->v.ref);
804 candidates->erase(j1);
813 }
while (candidates->size()>0);
816 sort(jets.begin(), jets.end(), jets_pt_less);
817 #ifdef DEBUG_SPLIT_MERGE
834 fprintf(flux,
"# %d jets found\n", (
int) jets.size());
835 fprintf(flux,
"# columns are: eta, phi, pt and number of particles for each jet\n");
836 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
839 fprintf(flux,
"%f\t%f\t%e\t%d\n",
843 fprintf(flux,
"# jet contents\n");
844 fprintf(flux,
"# columns are: eta, phi, pt, particle index and jet number\n");
845 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
847 for (i2=0;i2<j1->
n;i2++)
848 fprintf(flux,
"%f\t%f\t%e\t%d\t%d\n",
866 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
868 fprintf(stdout,
"jet %2d: %e\t%e\t%e\t%e\t", i1+1,
870 for (i2=0;i2<j->
n;i2++)
871 fprintf(stdout,
"%d ", j->
contents[i2]);
872 fprintf(stdout,
"\n");
875 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
877 fprintf(stdout,
"cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
879 for (i2=0;i2<c->
n;i2++)
880 fprintf(stdout,
"%d ", c->
contents[i2]);
881 fprintf(stdout,
"\n");
884 fprintf(stdout,
"\n");
895 bool Csplit_merge::get_overlap(
const Cjet &j1,
const Cjet &j2,
double *overlap2){
913 indices[idx_size] = j1.
contents[i1];
916 indices[idx_size] = j2.
contents[i2];
921 indices[idx_size] = j1.
contents[i1];
927 }
while ((i1<j1.
n) && (i2<j2.
n));
933 indices[idx_size] = j1.
contents[i1];
938 indices[idx_size] = j2.
contents[i2];
945 (*overlap2) = get_sm_var2(v, pt_tilde);
962 bool 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]){
1001 v = &(particles[j1.contents[i1]]);
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]){
1009 v = &(particles[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);
1017 v = &(particles[j1.contents[i1]]);
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;
1038 if (fabs(d1sq-d2sq) < most_ambiguous_split)
1039 most_ambiguous_split = fabs(d1sq-d2sq);
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));
1061 v = &(particles[j1.contents[i1]]);
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);
1069 v = &(particles[j2.contents[i2]]);
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);
1106 bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
1112 for (i=0;i<idx_size;i++){
1113 jet.contents.push_back(indices[i]);
1114 jet.v += particles[indices[i]];
1115 jet.pt_tilde += pt[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
1124 if (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);
1143 bool Csplit_merge::insert(Cjet &jet){
1148 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1149 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1154 if (jet.v.perp2()<pt_min2)
1158 jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
1161 candidates->insert(jet);
1172 double Csplit_merge::get_sm_var2(Cmomentum &v,
double &pt_tilde){
1173 switch(ptcomparison.split_merge_scale) {
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: "
1180 + ptcomparison.SM_scale_name());
class for holding a covering range in eta-phi
static double eta_max
maximal value for eta
static double eta_min
minimal value for eta
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
Ceta_phi_range range
covered range in eta-phi
std::vector< int > contents
particle contents (list of indices)
base class for dynamic coordinates management
double perp2() const
computes pT^2
Creference ref
reference number for the vector
int index
internal particle number
double eta
particle pseudo-rapidity
void build_etaphi()
build eta-phi from 4-momentum info !!! WARNING !!! !!! computing eta and phi is time-consuming !...
int parent_index
particle number in the parent list
double perp() const
computes pT
double phi
particle azimuthal angle
class corresponding to errors that will be thrown by siscone
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...
int init(std::vector< Cmomentum > &_particles, std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
initialisation function
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
int partial_clear()
partial clearance
int show()
show jets/candidates status
int init_pleft()
build initial list of left particles
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.
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
Csplit_merge()
default ctor
int full_clear()
full clearance
int init_particles(std::vector< Cmomentum > &_particles)
initialisation function for particle list
#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!