27 #include "split_merge.h"
28 #include "siscone_error.h"
66 bool jets_pt_less(
const Cjet &j1,
const Cjet &j2){
105 #ifdef EPSILON_SPLITMERGE
110 double pt_tilde_difference;
111 get_difference(jet1,jet2,&difference,&pt_tilde_difference);
120 switch (split_merge_scale){
122 qdiff = sum.
E*difference.
E - sum.
pz*difference.
pz;
125 qdiff = sum.
px*difference.
px + sum.
py*difference.
py;
128 qdiff = pt_tilde_sum*pt_tilde_difference;
135 qdiff = jet1.
v.
E*jet1.
v.
E*
136 ((sum.
px*difference.
px + sum.
py*difference.
py)*jet1.
v.
pz*jet1.
v.
pz
146 #endif // EPSILON_SPLITMERGE
154 std::string split_merge_scale_name(Esplit_merge_scale sms) {
157 return "pt (IR unsafe)";
159 return "Et (boost dep.)";
161 return "mt (IR safe except for pairs of identical decayed heavy particles)";
163 return "pttilde (scalar sum of pt's)";
165 return "[SM scale without a name]";
192 (*v) += (*particles)[j1.
contents[i1]];
193 (*pt_tilde) += (*pt)[j1.
contents[i1]];
196 (*v) -= (*particles)[j2.
contents[i2]];
197 (*pt_tilde) -= (*pt)[j2.
contents[i2]];
200 throw Csiscone_error(
"get_non_overlap reached part it should never have seen...");
202 }
while ((i1<j1.
n) && (i2<j2.
n));
206 (*v) += (*particles)[j1.
contents[i1]];
207 (*pt_tilde) += (*pt)[j1.
contents[i1]];
211 (*v) -= (*particles)[j2.
contents[i2]];
212 (*pt_tilde) -= (*pt)[j2.
contents[i2]];
225 merge_identical_protocones =
false;
226 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
227 #ifdef MERGE_IDENTICAL_PROTOCONES_DEFAULT_TRUE
228 merge_identical_protocones =
true;
234 ptcomparison.particles = &particles;
235 ptcomparison.pt = &pt;
236 candidates.reset(
new multiset<Cjet,Csplit_merge_ptcomparison>(ptcomparison));
239 SM_var2_hardest_cut_off = -1.0;
242 stable_cone_soft_pt2_cutoff = -1.0;
245 use_pt_weighted_splitting =
false;
262 int Csplit_merge::init(vector<Cmomentum> & , vector<Cmomentum> *protocones,
double R2,
double ptmin){
264 return add_protocones(protocones, R2, ptmin);
277 particles = _particles;
278 n = particles.size();
282 for (
int i=0;i<n;i++)
283 pt[i] = particles[i].perp();
286 ptcomparison.particles = &particles;
287 ptcomparison.pt = &pt;
292 indices =
new int[n];
318 particles[i].ref.randomize();
321 if (fabs(particles[i].pz) < (particles[i].E)){
322 p_remain.push_back(particles[i]);
324 p_remain[j].parent_index = i;
331 p_remain[j].index = 1;
335 particles[i].index = 0;
337 eta_min = min(eta_min, particles[i].eta);
338 eta_max = max(eta_max, particles[i].eta);
340 particles[i].index = -1;
343 n_left = p_remain.size();
350 merge_collinear_and_remove_soft();
367 candidates.reset(
new multiset<Cjet,Csplit_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<Cmomentum> 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_eta_less);
423 if (p_sorted[i].perp2()<stable_cone_soft_pt2_cutoff) {
431 while ((j<n_left) && (fabs(p_sorted[j].eta-p_sorted[i].eta)<
EPSILON_COLLINEAR) && (!collinear)){
432 dphi = fabs(p_sorted[j].phi-p_sorted[i].phi);
433 if (dphi>M_PI) dphi =
twopi-dphi;
436 p_sorted[j] += p_sorted[i];
444 p_uncol_hard.push_back(p_sorted[i]);
466 if (protocones->size()==0)
469 pt_min2 = ptmin*ptmin;
474 for (vector<Cmomentum>::iterator p_it = protocones->begin();p_it != protocones->end();p_it++){
487 for (i=0;i<n_left;i++){
493 dy = fabs(phi - v->
phi);
514 #ifdef DEBUG_SPLIT_MERGE
515 cout <<
"adding jet: ";
516 for (
int i2=0;i2<jet.
n;i2++)
528 #ifdef DEBUG_SPLIT_MERGE
529 cout <<
"remaining particles: ";
532 for (i=0;i<n_left;i++){
533 if (p_remain[i].index){
535 p_remain[j]=p_remain[i];
536 p_remain[j].parent_index = p_remain[i].parent_index;
539 particles[p_remain[j].parent_index].index = n_pass;
540 #ifdef DEBUG_SPLIT_MERGE
541 cout << p_remain[j].parent_index <<
" ";
546 #ifdef DEBUG_SPLIT_MERGE
552 merge_collinear_and_remove_soft();
572 pt_min2 = ptmin*ptmin;
574 if (candidates->size()==0)
577 if (overlap_tshold>=1.0 || overlap_tshold <= 0) {
578 ostringstream message;
579 message <<
"Illegal value for overlap_tshold, f = " << overlap_tshold;
580 message <<
" (legal values are 0<f<1)";
590 double overlap_tshold2 = overlap_tshold*overlap_tshold;
593 if (candidates->size()>0){
595 j1 = candidates->begin();
599 if (j1->sm_var2<SM_var2_hardest_cut_off) {
break;}
606 while (j2 != candidates->end()){
607 #ifdef DEBUG_SPLIT_MERGE
611 if (get_overlap(*j1, *j2, &overlap2)){
614 #ifdef DEBUG_SPLIT_MERGE
615 cout <<
"overlap between cdt 1 and cdt " << j2_relindex+1 <<
" with overlap "
616 << sqrt(overlap2/j2->sm_var2) << endl<<endl;
618 if (overlap2<overlap_tshold2*j2->sm_var2){
623 j2 = j1 = candidates->begin();
630 j2 = j1 = candidates->begin();
639 if (j2 != candidates->end()) j2++;
642 if (j1 != candidates->end()) {
647 jets[jets.size()-1].v.build_etaphi();
650 assert(j1->contents.size() > 0);
651 jets[jets.size()-1].pass = particles[j1->contents[0]].index;
652 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
653 cand_refs.erase(j1->v.ref);
655 candidates->erase(j1);
664 }
while (candidates->size()>0);
667 sort(jets.begin(), jets.end(), jets_pt_less);
668 #ifdef DEBUG_SPLIT_MERGE
685 fprintf(flux,
"# %d jets found\n", (
int) jets.size());
686 fprintf(flux,
"# columns are: eta, phi, pt and number of particles for each jet\n");
687 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
690 fprintf(flux,
"%f\t%f\t%e\t%d\n",
694 fprintf(flux,
"# jet contents\n");
695 fprintf(flux,
"# columns are: eta, phi, pt, particle index and jet number\n");
696 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
698 for (i2=0;i2<j1->
n;i2++)
699 fprintf(flux,
"%f\t%f\t%e\t%d\t%d\n",
717 for (it_j = jets.begin(), i1=0 ; it_j != jets.end() ; it_j++, i1++){
719 fprintf(stdout,
"jet %2d: %e\t%e\t%e\t%e\t", i1+1,
721 for (i2=0;i2<j->
n;i2++)
722 fprintf(stdout,
"%d ", j->
contents[i2]);
723 fprintf(stdout,
"\n");
726 for (it_c = candidates->begin(), i1=0 ; it_c != candidates->end() ; it_c++, i1++){
728 fprintf(stdout,
"cdt %2d: %e\t%e\t%e\t%e\t%e\t", i1+1,
730 for (i2=0;i2<c->
n;i2++)
731 fprintf(stdout,
"%d ", c->
contents[i2]);
732 fprintf(stdout,
"\n");
735 fprintf(stdout,
"\n");
746 bool Csplit_merge::get_overlap(
const Cjet &j1,
const Cjet &j2,
double *overlap2){
764 indices[idx_size] = j1.
contents[i1];
767 indices[idx_size] = j2.
contents[i2];
772 indices[idx_size] = j1.
contents[i1];
778 }
while ((i1<j1.
n) && (i2<j2.
n));
784 indices[idx_size] = j1.
contents[i1];
789 indices[idx_size] = j2.
contents[i2];
796 (*overlap2) = get_sm_var2(v, pt_tilde);
813 bool Csplit_merge::split(cjet_iterator &it_j1, cjet_iterator &it_j2){
816 double eta1, phi1, pt1_weight, eta2, phi2, pt2_weight;
817 double dx1, dy1, dx2, dy2;
822 const Cjet & j1 = * it_j1;
823 const Cjet & j2 = * it_j2;
826 jet2.
v = jet1.v = Cmomentum();
827 jet2.pt_tilde = jet1.pt_tilde = 0.0;
838 pt1_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
844 pt2_weight = (use_pt_weighted_splitting) ? 1.0/tmp.perp2() : 1.0;
846 jet1.v = jet2.v = Cmomentum();
850 if (j1.contents[i1]<j2.contents[i2]){
852 v = &(particles[j1.contents[i1]]);
853 jet1.contents.push_back(j1.contents[i1]);
855 jet1.pt_tilde += pt[j1.contents[i1]];
857 jet1.range.add_particle(v->eta,v->phi);
858 }
else if (j1.contents[i1]>j2.contents[i2]){
860 v = &(particles[j2.contents[i2]]);
861 jet2.contents.push_back(j2.contents[i2]);
863 jet2.pt_tilde += pt[j2.contents[i2]];
865 jet2.range.add_particle(v->eta,v->phi);
868 v = &(particles[j1.contents[i1]]);
872 dy1 = fabs(phi1 - v->phi);
878 dy2 = fabs(phi2 - v->phi);
886 double d1sq = (dx1*dx1+dy1*dy1)*pt1_weight;
887 double d2sq = (dx2*dx2+dy2*dy2)*pt2_weight;
889 if (fabs(d1sq-d2sq) < most_ambiguous_split)
890 most_ambiguous_split = fabs(d1sq-d2sq);
894 jet1.contents.push_back(j1.contents[i1]);
896 jet1.pt_tilde += pt[j1.contents[i1]];
897 jet1.range.add_particle(v->eta,v->phi);
900 jet2.contents.push_back(j2.contents[i2]);
902 jet2.pt_tilde += pt[j2.contents[i2]];
903 jet2.range.add_particle(v->eta,v->phi);
909 }
while ((i1<j1.n) && (i2<j2.n));
912 v = &(particles[j1.contents[i1]]);
913 jet1.contents.push_back(j1.contents[i1]);
915 jet1.pt_tilde += pt[j1.contents[i1]];
917 jet1.range.add_particle(v->eta,v->phi);
920 v = &(particles[j2.contents[i2]]);
921 jet2.contents.push_back(j2.contents[i2]);
923 jet2.pt_tilde += pt[j2.contents[i2]];
925 jet2.range.add_particle(v->eta,v->phi);
929 jet1.n = jet1.contents.size();
930 jet2.n = jet2.contents.size();
936 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
937 cand_refs.erase(j1.v.ref);
938 cand_refs.erase(j2.v.ref);
940 candidates->erase(it_j1);
941 candidates->erase(it_j2);
957 bool Csplit_merge::merge(cjet_iterator &it_j1, cjet_iterator &it_j2){
963 for (i=0;i<idx_size;i++){
964 jet.contents.push_back(indices[i]);
965 jet.v += particles[indices[i]];
966 jet.pt_tilde += pt[indices[i]];
968 jet.n = jet.contents.size();
971 jet.range = range_union(it_j1->range, it_j2->range);
974 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
975 if (merge_identical_protocones){
976 cand_refs.erase(it_j1->v.ref);
977 cand_refs.erase(it_j2->v.ref);
980 candidates->erase(it_j1);
981 candidates->erase(it_j2);
994 bool Csplit_merge::insert(Cjet &jet){
999 #ifdef ALLOW_MERGE_IDENTICAL_PROTOCONES
1000 if ((merge_identical_protocones) && (!cand_refs.insert(jet.v.ref).second))
1005 if (jet.v.perp2()<pt_min2)
1009 jet.sm_var2 = get_sm_var2(jet.v, jet.pt_tilde);
1012 candidates->insert(jet);
1023 double Csplit_merge::get_sm_var2(Cmomentum &v,
double &pt_tilde){
1024 switch(ptcomparison.split_merge_scale) {
1025 case SM_pt:
return v.perp2();
1026 case SM_mt:
return v.perpmass2();
1027 case SM_pttilde:
return pt_tilde*pt_tilde;
1028 case SM_Et:
return v.Et2();
1030 throw Csiscone_error(
"Unsupported split-merge scale choice: "
1031 + ptcomparison.SM_scale_name());