29 #include <siscone/siscone_error.h>
30 #include "split_merge.h"
39 namespace siscone_spherical{
58 pass = CJET_INEXISTENT_PASS;
71 return j1.
v.
E > j2.
v.
E;
109 #ifdef EPSILON_SPLITMERGE
112 #ifdef DEBUG_SPLIT_MERGE
113 cout <<
"Using high-precision ordering tests" << endl;
117 double E_tilde_difference;
118 get_difference(jet1,jet2,&difference,&E_tilde_difference);
127 switch (split_merge_scale){
129 qdiff = E_tilde_sum*E_tilde_difference;
132 qdiff = sum.
E*difference.
E;
148 std::string split_merge_scale_name(Esplit_merge_scale sms) {
151 return "E (IR unsafe for pairs of identical decayed heavy particles)";
153 return "Etilde (sum of E.[1+sin^2(theta_{i,jet})])";
155 return "[SM scale without a name]";
189 (*E_tilde) += p.
E*((norm2_cross_product3(p,jet1_axis)-norm2_cross_product3(p,jet2_axis))/(*particles_norm2)[j1.
contents[i1]]);
195 (*E_tilde) += p.
E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.
contents[i1]];
200 (*E_tilde) -= p.
E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.
contents[i2]];
205 }
while ((i1<j1.
n) && (i2<j2.
n));
211 (*E_tilde) += p.
E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.
contents[i1]];
217 (*E_tilde) -= p.
E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.
contents[i2]];
233 merge_identical_protocones =
false;
234 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
235 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
236 merge_identical_protocones =
true;
243 ptcomparison.particles = &particles;
244 ptcomparison.particles_norm2 = &particles_norm2;
245 candidates.reset(
new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
248 SM_var2_hardest_cut_off = -numeric_limits<double>::max();
251 stable_cone_soft_E2_cutoff = -1.0;
254 use_E_weighted_splitting =
false;
273 return add_protocones(protocones, R2, Emin);
286 particles = _particles;
287 n = particles.size();
290 particles_norm2.resize(n);
291 for (
int i=0;i<n;i++){
292 particles_norm2[i] = particles[i].norm2();
296 ptcomparison.particles = &particles;
297 ptcomparison.particles_norm2 = &particles_norm2;
302 indices =
new int[n];
325 particles[i].ref.randomize();
329 p_remain.push_back(particles[i]);
331 p_remain[j].parent_index = i;
338 p_remain[j].index = 1;
342 particles[i].index = 0;
347 n_left = p_remain.size();
350 merge_collinear_and_remove_soft();
367 candidates.reset(
new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
370 most_ambiguous_split = numeric_limits<double>::max();
373 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
374 if (merge_identical_protocones)
390 if (indices != NULL){
405 vector<CSphmomentum> p_sorted;
409 p_uncol_hard.clear();
412 for (i=0;i<n_left;i++)
413 p_sorted.push_back(p_remain[i]);
414 sort(p_sorted.begin(), p_sorted.end(), momentum_theta_less);
423 if (p_sorted[i].E*p_sorted[i].E<stable_cone_soft_E2_cutoff) {
431 while ((j<n_left) && (fabs(p_sorted[j]._theta-p_sorted[i]._theta)<
EPSILON_COLLINEAR) && (!collinear)){
432 dphi = fabs(p_sorted[j]._phi-p_sorted[i]._phi);
433 if (dphi>M_PI) dphi =
twopi-dphi;
436 #ifdef DEBUG_SPLIT_MERGE
437 cout <<
"# collinear merging at point " << p_sorted[i]._theta <<
", " << p_sorted[j]._phi << endl;
439 p_sorted[j] += p_sorted[i];
441 p_sorted[j].build_norm();
449 p_uncol_hard.push_back(p_sorted[i]);
469 if (protocones->size()==0)
477 #ifdef DEBUG_SPLIT_MERGE
478 cout <<
"particle list: ";
479 for (
int i2=0;i2<n_left;i2++)
480 cout << p_remain[i2].parent_index <<
" "
481 << p_remain[i2].px <<
" " << p_remain[i2].py <<
" "
482 << p_remain[i2].pz <<
" " << p_remain[i2].E << endl;
488 for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
496 for (i=0;i<n_left;i++){
498 if (is_closer(v, c, tan2R)){
519 #ifdef DEBUG_SPLIT_MERGE
520 cout <<
"adding protojet: ";
523 for (
unsigned int i2=0;i2<32;i2++) fprintf(stdout,
"%d", (phirange&(1<<i2)) >> i2 );
524 fprintf(stdout,
"\t");
526 for (
unsigned int i2=0;i2<32;i2++) fprintf(stdout,
"%d", (thetarange&(1<<i2)) >> i2);
527 fprintf(stdout,
"\t");
529 for (
int i2=0;i2<jet.
n;i2++)
541 #ifdef DEBUG_SPLIT_MERGE
542 cout <<
"remaining particles: ";
545 for (i=0;i<n_left;i++){
546 if (p_remain[i].index){
548 p_remain[j]=p_remain[i];
549 p_remain[j].parent_index = p_remain[i].parent_index;
552 particles[p_remain[j].parent_index].index = n_pass;
553 #ifdef DEBUG_SPLIT_MERGE
554 cout << p_remain[j].parent_index <<
" ";
559 #ifdef DEBUG_SPLIT_MERGE
565 merge_collinear_and_remove_soft();
587 bool found_jet =
false;
589 if (protocones->size()==0)
599 for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
607 for (i=0;i<n_left;i++){
609 if (is_closer(v, c, tan2R)){
611 jet_candidate.
v+= *v;
615 jet_candidate.
n=jet_candidate.
contents.size();
620 compute_Etilde(jet_candidate);
624 *c = jet_candidate.
v;
631 if (jet_candidate.
v.
E<E_min)
638 jet_candidate.
sm_var2 = (*_user_scale)(jet_candidate);
641 jet_candidate.
sm_var2 = get_sm_var2(jet_candidate.
v, jet_candidate.
E_tilde);
646 (_user_scale ? _user_scale->is_larger(jet_candidate, jet)
647 : ptcomparison(jet_candidate, jet))){
654 if (!found_jet)
return 1;
659 jets[jets.size()-1].v.build_norm();
661 #ifdef DEBUG_SPLIT_MERGE
662 cout <<
"PR-Jet " << jets.size() <<
" [size " << jet.
contents.size() <<
"]:";
666 int p_remain_index = 0;
667 int contents_index = 0;
669 for (
int index=0;index<n_left;index++){
670 if ((contents_index<(
int) jet.
contents.size()) &&
671 (p_remain[index].parent_index == jet.
contents[contents_index])){
674 particles[p_remain[index].parent_index].index = n_pass;
675 #ifdef DEBUG_SPLIT_MERGE
676 cout <<
" " << jet.
contents[contents_index];
681 p_remain[p_remain_index] = p_remain[index];
682 p_remain[p_remain_index].parent_index = p_remain[index].parent_index;
683 p_remain[p_remain_index].index=1;
687 p_remain.resize(n_left-jet.
contents.size());
688 n_left = p_remain.size();
689 jets[jets.size()-1].pass = particles[jet.
contents[0]].index;
694 #ifdef DEBUG_SPLIT_MERGE
701 merge_collinear_and_remove_soft();
722 if (candidates->size()==0)
725 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
726 ostringstream message;
727 message <<
"Illegal value for overlap_tshold, f = " << overlap_tshold;
728 message <<
" (legal values are 0<f<1)";
738 double overlap_tshold2 = overlap_tshold*overlap_tshold;
741 if (candidates->size()>0){
743 j1 = candidates->begin();
748 if (j1->sm_var2<SM_var2_hardest_cut_off) {
break;}
755 while (j2 != candidates->end()){
756 #ifdef DEBUG_SPLIT_MERGE
757 if (j2_relindex==1) show();
758 cout <<
"check overlap between cdt 1 and cdt " << j2_relindex+1 <<
" with overlap " << endl;
761 if (get_overlap(*j1, *j2, &overlap2)){
764 #ifdef DEBUG_SPLIT_MERGE
765 cout <<
"overlap between cdt 1 and cdt " << j2_relindex+1 <<
" with overlap "
766 << sqrt(overlap2)/j2->v.E << endl<<endl;
769 if (overlap2<overlap_tshold2*sqr(j2->v.E)){
770 #ifdef DEBUG_SPLIT_MERGE
771 cout <<
" --> split" << endl<<endl;
777 j2 = j1 = candidates->begin();
780 #ifdef DEBUG_SPLIT_MERGE
781 cout <<
" --> merge" << endl<<endl;
787 j2 = j1 = candidates->begin();
796 if (j2 != candidates->end()) j2++;
799 if (j1 != candidates->end()) {
804 jets[jets.size()-1].v.build_thetaphi();
805 jets[jets.size()-1].v.build_norm();
808 assert(j1->contents.size() > 0);
809 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
810 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
811 cand_refs.erase(j1->v.ref);
813 candidates->erase(j1);
816 }
while (candidates->size()>0);
819 sort(jets.begin(), jets.end(), jets_E_less);
820 #ifdef DEBUG_SPLIT_MERGE
837 fprintf(flux,
"# %d jets found\n", (
int) jets.size());
838 fprintf(flux,
"# columns are: px, py, pz, E and number of particles for each jet\n");
839 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
841 fprintf(flux,
"%e\t%e\t%e\t%e\t%d\n",
845 fprintf(flux,
"# jet contents\n");
846 fprintf(flux,
"# columns are: px, py, pz, E, particle index and jet number\n");
847 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
849 for (i2=0;i2<j1->
n;i2++)
850 fprintf(flux,
"%e\t%e\t%e\t%e\t%d\t%d\n",
869 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
871 fprintf(stdout,
"jet %2d: %e\t%e\t%e\t%e\t", i1+1,
875 for (i2=0;i2<32;i2++) fprintf(stdout,
"%d", (phirange&(1<<i2)) >> i2 );
876 fprintf(stdout,
"\t");
878 for (i2=0;i2<32;i2++) fprintf(stdout,
"%d", (thetarange&(1<<i2)) >> i2);
879 fprintf(stdout,
"\t");
881 for (i2=0;i2<j->
n;i2++)
882 fprintf(stdout,
"%d ", j->
contents[i2]);
883 fprintf(stdout,
"\n");
886 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
888 fprintf(stdout,
"cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
892 for (i2=0;i2<32;i2++) fprintf(stdout,
"%d", (phirange&(1<<i2)) >> i2 );
893 fprintf(stdout,
"\t");
895 for (i2=0;i2<32;i2++) fprintf(stdout,
"%d", (thetarange&(1<<i2)) >> i2);
896 fprintf(stdout,
"\t");
898 for (i2=0;i2<c->
n;i2++)
899 fprintf(stdout,
"%d ", c->
contents[i2]);
900 fprintf(stdout,
"\n");
903 fprintf(stdout,
"\n");
914 bool CSphsplit_merge::get_overlap(
const CSphjet &j1,
const CSphjet &j2,
double *overlap2){
931 indices[idx_size] = j1.
contents[i1];
934 indices[idx_size] = j2.
contents[i2];
938 indices[idx_size] = j2.
contents[i2];
944 }
while ((i1<j1.
n) && (i2<j2.
n));
950 indices[idx_size] = j1.
contents[i1];
955 indices[idx_size] = j2.
contents[i2];
962 (*overlap2) = sqr(v.
E);
979 bool CSphsplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
982 double E1_weight, E2_weight;
987 const CSphjet & j1 = * it_j1;
988 const CSphjet & j2 = * it_j2;
991 jet2.
v = jet1.v = CSphmomentum();
998 E1_weight = (use_E_weighted_splitting) ? 1.0/j1.v.E/j1.v.E : 1.0;
999 E2_weight = (use_E_weighted_splitting) ? 1.0/j2.v.E/j2.v.E : 1.0;
1003 if (j1.contents[i1]<j2.contents[i2]){
1005 v = &(particles[j1.contents[i1]]);
1006 jet1.contents.push_back(j1.contents[i1]);
1010 jet1.range.add_particle(v->_theta,v->_phi);
1011 }
else if (j1.contents[i1]>j2.contents[i2]){
1013 v = &(particles[j2.contents[i2]]);
1014 jet2.contents.push_back(j2.contents[i2]);
1018 jet2.range.add_particle(v->_theta,v->_phi);
1021 v = &(particles[j1.contents[i1]]);
1029 double d1 = get_distance(&(j1.v), v)*E1_weight;
1030 double d2 = get_distance(&(j2.v), v)*E2_weight;
1032 if (fabs(d1-d2) < most_ambiguous_split)
1033 most_ambiguous_split = fabs(d1-d2);
1037 jet1.contents.push_back(j1.contents[i1]);
1040 jet1.range.add_particle(v->_theta,v->_phi);
1043 jet2.contents.push_back(j2.contents[i2]);
1046 jet2.range.add_particle(v->_theta,v->_phi);
1052 }
while ((i1<j1.n) && (i2<j2.n));
1055 v = &(particles[j1.contents[i1]]);
1056 jet1.contents.push_back(j1.contents[i1]);
1060 jet1.range.add_particle(v->_theta,v->_phi);
1063 v = &(particles[j2.contents[i2]]);
1064 jet2.contents.push_back(j2.contents[i2]);
1068 jet2.range.add_particle(v->_theta,v->_phi);
1072 jet1.n = jet1.contents.size();
1073 jet2.n = jet2.contents.size();
1076 compute_Etilde(jet1);
1077 compute_Etilde(jet2);
1080 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1081 cand_refs.erase(j1.v.ref);
1082 cand_refs.erase(j2.v.ref);
1084 candidates->erase(it_j1);
1085 candidates->erase(it_j2);
1101 bool CSphsplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
1107 for (i=0;i<idx_size;i++){
1108 jet.contents.push_back(indices[i]);
1109 jet.v += particles[indices[i]];
1112 jet.n = jet.contents.size();
1114 compute_Etilde(jet);
1117 jet.range = range_union(it_j1->range, it_j2->range);
1120 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1121 if (merge_identical_protocones){
1122 cand_refs.erase(it_j1->v.ref);
1123 cand_refs.erase(it_j2->v.ref);
1126 candidates->erase(it_j1);
1127 candidates->erase(it_j2);
1140 bool CSphsplit_merge::insert(CSphjet &jet){
1145 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1146 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1155 jet.sm_var2 = get_sm_var2(jet.v, jet.E_tilde);
1158 candidates->insert(jet);
1169 double CSphsplit_merge::get_sm_var2(CSphmomentum &v,
double &E_tilde){
1170 switch(ptcomparison.split_merge_scale) {
1171 case SM_E:
return v.E*v.E;
1172 case SM_Etilde:
return E_tilde*E_tilde;
1175 + ptcomparison.SM_scale_name());
1184 void CSphsplit_merge::compute_Etilde(CSphjet &jet){
1187 CSph3vector jet_axis = jet.v;
1195 for (vector<int>::iterator cont_it=jet.contents.begin(); cont_it!=jet.contents.end(); cont_it++){
1196 const CSphmomentum &p = particles[*cont_it];
1197 jet.E_tilde+=p.E*(1.0+norm2_cross_product3(p,jet_axis)/particles_norm2[*cont_it]);
class corresponding to errors that will be thrown by siscone
base class for managing the spatial part of Cmomentum (defined after)
double _theta
particle theta angle (available ONLY after a call to build_thetaphi)
double _phi
particle phi angle (available ONLY after a call to build_thetaphi)
siscone::Creference ref
reference number for the vector
void build_thetaphi()
just a useful tool to store theta and phi locally (in _theta and _phi) in case you need repeated acce...
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
double E_tilde
sum of E_i [ 1 +|p_i x p_J|^2/(|p_i|^2 E_J^2)]
CSphmomentum v
jet momentum
std::vector< int > contents
particle contents (list of indices)
CSphtheta_phi_range range
covered range in eta-phi
int n
number of particles inside
base class for dynamic coordinates management
int index
internal particle number
int parent_index
particle number in the parent list
bool operator()(const CSphjet &jet1, const CSphjet &jet2) const
comparison of 2 CSphjet
void get_difference(const CSphjet &j1, const CSphjet &j2, CSphmomentum *v, double *E_tilde) const
get the difference between 2 jets, calculated such that rounding errors will not affect the result ev...
int show()
show jets/candidates status
int perform(double overlap_tshold, double Emin=0.0)
really do the splitting and merging At the end, the vector jets is filled with the jets found.
int partial_clear()
partial clearance
int add_hardest_protocone_to_jets(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
remove the hardest protocone and declare it a a jet
int full_clear()
full clearance
CSphsplit_merge()
default ctor
int init_pleft()
build initial list of left particles
int init(std::vector< CSphmomentum > &_particles, std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
initialisation function
~CSphsplit_merge()
default dtor
int save_contents(FILE *flux)
save final jets
int merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles note that thins in only...
int init_particles(std::vector< CSphmomentum > &_particles)
initialisation function for particle list
int add_protocones(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
add a list of protocones
class for holding a covering range in eta-phi
unsigned int theta_range
theta range as a binary coding of covered cells
unsigned int phi_range
phi range as a binary coding of covered cells
#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!