29 #include <siscone/siscone_error.h>
30 #include "split_merge.h"
69 return j1.
v.
E > j2.
v.
E;
107 #ifdef EPSILON_SPLITMERGE
110 #ifdef DEBUG_SPLIT_MERGE
111 cout <<
"Using high-precision ordering tests" << endl;
115 double E_tilde_difference;
116 get_difference(jet1,jet2,&difference,&E_tilde_difference);
125 switch (split_merge_scale){
127 qdiff = E_tilde_sum*E_tilde_difference;
130 qdiff = sum.
E*difference.
E;
138 #endif // EPSILON_SPLITMERGE
146 std::string split_merge_scale_name(Esplit_merge_scale sms) {
149 return "E (IR unsafe for pairs of identical decayed heavy particles)";
151 return "Etilde (sum of E.[1+sin^2(theta_{i,jet})])";
153 return "[SM scale without a name]";
187 (*E_tilde) += p.
E*((norm2_cross_product3(p,jet1_axis)-norm2_cross_product3(p,jet2_axis))/(*particles_norm2)[j1.
contents[i1]]);
193 (*E_tilde) += p.
E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.
contents[i1]];
198 (*E_tilde) -= p.
E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.
contents[i2]];
203 }
while ((i1<j1.
n) && (i2<j2.
n));
209 (*E_tilde) += p.
E*norm2_cross_product3(p,jet1_axis)/(*particles_norm2)[j1.
contents[i1]];
215 (*E_tilde) -= p.
E*norm2_cross_product3(p,jet2_axis)/(*particles_norm2)[j2.
contents[i2]];
231 merge_identical_protocones =
false;
232 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
233 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
234 merge_identical_protocones =
true;
240 ptcomparison.particles = &particles;
241 ptcomparison.particles_norm2 = &particles_norm2;
242 candidates.reset(
new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
245 SM_var2_hardest_cut_off = -1.0;
248 stable_cone_soft_E2_cutoff = -1.0;
251 use_E_weighted_splitting =
false;
270 return add_protocones(protocones, R2, Emin);
283 particles = _particles;
284 n = particles.size();
287 particles_norm2.resize(n);
288 for (
int i=0;i<n;i++){
289 particles_norm2[i] = particles[i].norm2();
293 ptcomparison.particles = &particles;
294 ptcomparison.particles_norm2 = &particles_norm2;
299 indices =
new int[n];
322 particles[i].ref.randomize();
326 p_remain.push_back(particles[i]);
328 p_remain[j].parent_index = i;
335 p_remain[j].index = 1;
339 particles[i].index = 0;
344 n_left = p_remain.size();
347 merge_collinear_and_remove_soft();
364 candidates.reset(
new multiset<CSphjet,CSphsplit_merge_ptcomparison>(ptcomparison));
367 most_ambiguous_split = numeric_limits<double>::max();
370 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
371 if (merge_identical_protocones)
387 if (indices != NULL){
402 vector<CSphmomentum> p_sorted;
406 p_uncol_hard.clear();
409 for (i=0;i<n_left;i++)
410 p_sorted.push_back(p_remain[i]);
411 sort(p_sorted.begin(), p_sorted.end(), momentum_theta_less);
420 if (p_sorted[i].E*p_sorted[i].E<stable_cone_soft_E2_cutoff) {
428 while ((j<n_left) && (fabs(p_sorted[j]._theta-p_sorted[i]._theta)<
EPSILON_COLLINEAR) && (!collinear)){
429 dphi = fabs(p_sorted[j]._phi-p_sorted[i]._phi);
430 if (dphi>M_PI) dphi =
twopi-dphi;
433 #ifdef DEBUG_SPLIT_MERGE
434 cout <<
"# collinear merging at point " << p_sorted[i]._theta <<
", " << p_sorted[j]._phi << endl;
436 p_sorted[j] += p_sorted[i];
438 p_sorted[j].build_norm();
446 p_uncol_hard.push_back(p_sorted[i]);
466 if (protocones->size()==0)
474 #ifdef DEBUG_SPLIT_MERGE
475 cout <<
"particle list: ";
476 for (
int i2=0;i2<n_left;i2++)
477 cout << p_remain[i2].parent_index <<
" "
478 << p_remain[i2].px <<
" " << p_remain[i2].py <<
" "
479 << p_remain[i2].pz <<
" " << p_remain[i2].E << endl;
485 for (vector<CSphmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
493 for (i=0;i<n_left;i++){
495 if (is_closer(v, c, tan2R)){
516 #ifdef DEBUG_SPLIT_MERGE
517 cout <<
"adding protojet: ";
518 for (
int i2=0;i2<jet.
n;i2++)
530 #ifdef DEBUG_SPLIT_MERGE
531 cout <<
"remaining particles: ";
534 for (i=0;i<n_left;i++){
535 if (p_remain[i].index){
537 p_remain[j]=p_remain[i];
538 p_remain[j].parent_index = p_remain[i].parent_index;
541 particles[p_remain[j].parent_index].index = n_pass;
542 #ifdef DEBUG_SPLIT_MERGE
543 cout << p_remain[j].parent_index <<
" ";
548 #ifdef DEBUG_SPLIT_MERGE
554 merge_collinear_and_remove_soft();
576 if (candidates->size()==0)
579 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
580 ostringstream message;
581 message <<
"Illegal value for overlap_tshold, f = " << overlap_tshold;
582 message <<
" (legal values are 0<f<1)";
592 double overlap_tshold2 = overlap_tshold*overlap_tshold;
595 if (candidates->size()>0){
597 j1 = candidates->begin();
602 if (j1->sm_var2<SM_var2_hardest_cut_off) {
break;}
609 while (j2 != candidates->end()){
610 #ifdef DEBUG_SPLIT_MERGE
614 if (get_overlap(*j1, *j2, &overlap2)){
617 #ifdef DEBUG_SPLIT_MERGE
618 cout <<
"overlap between cdt 1 and cdt " << j2_relindex+1 <<
" with overlap "
619 << sqrt(overlap2)/j2->v.E << endl<<endl;
622 if (overlap2<overlap_tshold2*sqr(j2->v.E)){
627 j2 = j1 = candidates->begin();
634 j2 = j1 = candidates->begin();
643 if (j2 != candidates->end()) j2++;
646 if (j1 != candidates->end()) {
651 jets[jets.size()-1].v.build_thetaphi();
652 jets[jets.size()-1].v.build_norm();
655 assert(j1->contents.size() > 0);
656 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
657 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
658 cand_refs.erase(j1->v.ref);
660 candidates->erase(j1);
663 }
while (candidates->size()>0);
666 sort(jets.begin(), jets.end(), jets_E_less);
667 #ifdef DEBUG_SPLIT_MERGE
684 fprintf(flux,
"# %d jets found\n", (
int) jets.size());
685 fprintf(flux,
"# columns are: px, py, pz, E and number of particles for each jet\n");
686 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
688 fprintf(flux,
"%e\t%e\t%e\t%e\t%d\n",
692 fprintf(flux,
"# jet contents\n");
693 fprintf(flux,
"# columns are: px, py, pz, E, particle index and jet number\n");
694 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
696 for (i2=0;i2<j1->
n;i2++)
697 fprintf(flux,
"%e\t%e\t%e\t%e\t%d\t%d\n",
716 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
718 fprintf(stdout,
"jet %2d: %e\t%e\t%e\t%e\t", i1+1,
720 for (i2=0;i2<j->
n;i2++)
721 fprintf(stdout,
"%d ", j->
contents[i2]);
722 fprintf(stdout,
"\n");
725 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
727 fprintf(stdout,
"cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
729 for (i2=0;i2<c->
n;i2++)
730 fprintf(stdout,
"%d ", c->
contents[i2]);
731 fprintf(stdout,
"\n");
734 fprintf(stdout,
"\n");
745 bool CSphsplit_merge::get_overlap(
const CSphjet &j1,
const CSphjet &j2,
double *overlap2){
762 indices[idx_size] = j1.
contents[i1];
765 indices[idx_size] = j2.
contents[i2];
769 indices[idx_size] = j2.
contents[i2];
775 }
while ((i1<j1.
n) && (i2<j2.
n));
781 indices[idx_size] = j1.
contents[i1];
786 indices[idx_size] = j2.
contents[i2];
793 (*overlap2) = sqr(v.
E);
810 bool CSphsplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
813 double E1_weight, E2_weight;
818 const CSphjet & j1 = * it_j1;
819 const CSphjet & j2 = * it_j2;
822 jet2.
v = jet1.v = CSphmomentum();
829 E1_weight = (use_E_weighted_splitting) ? 1.0/j1.v.E/j1.v.E : 1.0;
830 E2_weight = (use_E_weighted_splitting) ? 1.0/j2.v.E/j2.v.E : 1.0;
834 if (j1.contents[i1]<j2.contents[i2]){
836 v = &(particles[j1.contents[i1]]);
837 jet1.contents.push_back(j1.contents[i1]);
841 jet1.range.add_particle(v->_theta,v->_phi);
842 }
else if (j1.contents[i1]>j2.contents[i2]){
844 v = &(particles[j2.contents[i2]]);
845 jet2.contents.push_back(j2.contents[i2]);
849 jet2.range.add_particle(v->_theta,v->_phi);
852 v = &(particles[j1.contents[i1]]);
860 double d1 = get_distance(&(j1.v), v)*E1_weight;
861 double d2 = get_distance(&(j2.v), v)*E2_weight;
863 if (fabs(d1-d2) < most_ambiguous_split)
864 most_ambiguous_split = fabs(d1-d2);
868 jet1.contents.push_back(j1.contents[i1]);
871 jet1.range.add_particle(v->_theta,v->_phi);
874 jet2.contents.push_back(j2.contents[i2]);
877 jet2.range.add_particle(v->_theta,v->_phi);
883 }
while ((i1<j1.n) && (i2<j2.n));
886 v = &(particles[j1.contents[i1]]);
887 jet1.contents.push_back(j1.contents[i1]);
891 jet1.range.add_particle(v->_theta,v->_phi);
894 v = &(particles[j2.contents[i2]]);
895 jet2.contents.push_back(j2.contents[i2]);
899 jet2.range.add_particle(v->_theta,v->_phi);
903 jet1.n = jet1.contents.size();
904 jet2.n = jet2.contents.size();
907 compute_Etilde(jet1);
908 compute_Etilde(jet2);
911 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
912 cand_refs.erase(j1.v.ref);
913 cand_refs.erase(j2.v.ref);
915 candidates->erase(it_j1);
916 candidates->erase(it_j2);
932 bool CSphsplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
938 for (i=0;i<idx_size;i++){
939 jet.contents.push_back(indices[i]);
940 jet.v += particles[indices[i]];
943 jet.n = jet.contents.size();
948 jet.range = range_union(it_j1->range, it_j2->range);
951 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
952 if (merge_identical_protocones){
953 cand_refs.erase(it_j1->v.ref);
954 cand_refs.erase(it_j2->v.ref);
957 candidates->erase(it_j1);
958 candidates->erase(it_j2);
971 bool CSphsplit_merge::insert(CSphjet &jet){
976 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
977 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
986 jet.sm_var2 = get_sm_var2(jet.v, jet.E_tilde);
989 candidates->insert(jet);
1000 double CSphsplit_merge::get_sm_var2(CSphmomentum &v,
double &E_tilde){
1001 switch(ptcomparison.split_merge_scale) {
1002 case SM_E:
return v.E*v.E;
1003 case SM_Etilde:
return E_tilde*E_tilde;
1006 + ptcomparison.SM_scale_name());
1015 void CSphsplit_merge::compute_Etilde(CSphjet &jet){
1018 CSph3vector jet_axis = jet.v;
1026 for (vector<int>::iterator cont_it=jet.contents.begin(); cont_it!=jet.contents.end(); cont_it++){
1027 const CSphmomentum &p = particles[*cont_it];
1028 jet.E_tilde+=p.E*(1.0+norm2_cross_product3(p,jet_axis)/particles_norm2[*cont_it]);
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 n
number of particles inside
CSphsplit_merge()
default ctor
bool operator()(const CSphjet &jet1, const CSphjet &jet2) const
comparison of 2 CSphjet
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
int merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles note that thins in only...
class corresponding to errors that will be thrown by siscone
double _phi
particle phi angle (available ONLY after a call to build_thetaphi)
#define EPSILON_SPLITMERGE
The following define enables you to allow for identical protocones to be merged automatically after e...
int full_clear()
full clearance
siscone::Creference ref
reference number for the vector
int init_particles(std::vector< CSphmomentum > &_particles)
initialisation function for particle list
const double twopi
definition of 2*M_PI which is useful a bit everyhere!
int save_contents(FILE *flux)
save final jets
CSphtheta_phi_range range
covered range in eta-phi
void build_thetaphi()
just a useful tool to store theta and phi locally (in _theta and _phi) in case you need repeated acce...
int add_protocones(std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
add a list of protocones
CSphmomentum v
jet momentum
class for holding a covering range in eta-phi
double E_tilde
sum of E_i [ 1 +|p_i x p_J|^2/(|p_i|^2 E_J^2)]
int parent_index
particle number in the parent list
int init(std::vector< CSphmomentum > &_particles, std::vector< CSphmomentum > *protocones, double R2, double Emin=0.0)
initialisation function
~CSphsplit_merge()
default dtor
std::vector< int > contents
particle contents (list of indices)
double _theta
particle theta angle (available ONLY after a call to build_thetaphi)
int show()
show jets/candidates status
int index
internal particle number
int partial_clear()
partial clearance
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 init_pleft()
build initial list of left particles
base class for dynamic coordinates management
#define EPSILON_COLLINEAR
The following parameter controls collinear safety.
base class for managing the spatial part of Cmomentum (defined after)