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;
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;
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);
909 }
while ((i1<j1.
n) && (i2<j2.
n));
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++){
965 jet.
v += particles[indices[i]];
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)
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();
1027 case SM_pttilde:
return pt_tilde*pt_tilde;
1028 case SM_Et:
return v.
Et2();
1031 + ptcomparison.SM_scale_name());
double Et2() const
computes transverse energy (squared)
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...
bool operator()(const Cjet &jet1, const Cjet &jet2) const
comparison between 2 jets
int parent_index
particle number in the parent list
base class for dynamic coordinates management
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 merge_collinear_and_remove_soft()
build the list 'p_uncol_hard' from p_remain by clustering collinear particles and removing particles ...
Creference ref
reference number for the vector
double perpmass2() const
transverse mass squared, mt^2 = pt^2+m^2 = E^2 - pz^2
class corresponding to errors that will be thrown by siscone
int init_particles(std::vector< Cmomentum > &_particles)
initialisation function for particle list
int save_contents(FILE *flux)
save final jets
#define EPSILON_SPLITMERGE
The following define enables you to allow for identical protocones to be merged automatically after e...
int init_pleft()
build initial list of left particles
int add_protocones(std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
add a list of protocones
const double twopi
definition of 2*M_PI which is useful a bit everyhere!
static double eta_max
maximal value for eta
double sm_var2
ordering variable used for ordering and overlap in the split–merge.
class for holding a covering range in eta-phi
int full_clear()
full clearance
double phi
particle azimuthal angle
int add_particle(const double eta, const double phi)
add a particle to the range
void build_etaphi()
build eta-phi from 4-momentum info !!! WARNING !!! !!! computing eta and phi is time-consuming !!! !!...
int partial_clear()
partial clearance
double perp2() const
computes pT^2
double pt_tilde
p-scheme pt
int show()
show jets/candidates status
int init(std::vector< Cmomentum > &_particles, std::vector< Cmomentum > *protocones, double R2, double ptmin=0.0)
initialisation function
double perp() const
computes pT
int index
internal particle number
Csplit_merge()
default ctor
~Csplit_merge()
default dtor
#define EPSILON_COLLINEAR
The following parameter controls collinear safety.
static double eta_min
minimal value for eta
double eta
particle pseudo-rapidity
Ceta_phi_range range
covered range in eta-phi
std::vector< int > contents
particle contents (list of indices)
int n
number of particles inside