diff --git a/t/hej_test.cc b/t/hej_test.cc index 389b7ed..9dca54e 100644 --- a/t/hej_test.cc +++ b/t/hej_test.cc @@ -1,625 +1,641 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #include #include #include #include #include #include #include #include "HEJ/EWConstants.hh" #include "HEJ/Particle.hh" #include "HEJ/PDG_codes.hh" namespace { const HEJ::ParticleProperties Wprop{80.385, 2.085}; const HEJ::ParticleProperties Zprop{91.187, 2.495}; const HEJ::ParticleProperties Hprop{125, 0.004165}; constexpr double vev = 246.2196508; const HEJ::EWConstants ew_parameters{vev, Wprop, Zprop, Hprop}; } HEJ::Event::EventData get_process(int const njet, int const pos_boson){ using namespace HEJ::pid; HEJ::Event::EventData ev; if(njet == 0){ // jet idx: -1 -1 ev.outgoing.push_back({gluon, { -24, 12, -57, 63}, {}}); ev.outgoing.push_back({gluon, { 24, -12, 41, 49}, {}}); ev.incoming[0] = {gluon, { 0, 0, -64, 64}, {}}; ev.incoming[1] = {gluon, { 0, 0, 48, 48}, {}}; return ev; } else if(njet == 1){ // jet idx: 0 -1 -1 ev.outgoing.push_back({gluon, { 23, 28, -44, 57}, {}}); ev.outgoing.push_back({gluon, { -11, -24, -12, 29}, {}}); ev.outgoing.push_back({gluon, { -12, -4, 39, 41}, {}}); ev.incoming[0] = {gluon, { 0, 0, -72, 72}, {}}; ev.incoming[1] = {gluon, { 0, 0, 55, 55}, {}}; return ev; } else if(njet == 2){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 198, 33, -170, 291}, {}}); ev.outgoing.push_back({gluon, {-154, 68, 44, 174}, {}}); ev.outgoing.push_back({gluon, { -44, -101, 88, 141}, {}}); ev.incoming[0] = {gluon, { 0, 0, -322, 322}, {}}; ev.incoming[1] = {gluon, { 0, 0, 284, 284}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { -6, 82, -159, 179}, {}}); ev.outgoing.push_back({higgs, { 195, -106, 74, 265}, {}}); ev.outgoing.push_back({gluon, {-189, 24, 108, 219}, {}}); ev.incoming[0] = {gluon, { 0, 0, -320, 320}, {}}; ev.incoming[1] = {gluon, { 0, 0, 343, 343}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { -80, -80, -140, 180}, {}}); ev.outgoing.push_back({gluon, { -60, -32, 0, 68}, {}}); ev.outgoing.push_back({higgs, { 140, 112, 177, 281}, {}}); ev.incoming[0] = {gluon, { 0, 0, -246, 246}, {}}; ev.incoming[1] = {gluon, { 0, 0, 283, 283}, {}}; return ev; default: ev.outgoing.push_back({gluon, { -72, 24, 18, 78}, {}}); ev.outgoing.push_back({gluon, { 72, -24, 74, 106}, {}}); ev.incoming[0] = {gluon, { 0, 0, -46, 46}, {}}; ev.incoming[1] = {gluon, { 0, 0, 138, 138}, {}}; return ev; } } if(njet == 3){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 152, -117, -88, 245}, {}}); ev.outgoing.push_back({gluon, {-146, 62, -11, 159}, {}}); ev.outgoing.push_back({gluon, { 126, -72, 96, 174}, {}}); ev.outgoing.push_back({gluon, {-132, 127, 144, 233}, {}}); ev.incoming[0] = {gluon, { 0, 0, -335, 335}, {}}; ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, {-191, 188, -128, 297}, {}}); ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, { 184, -172, -8, 252}, {}}); ev.outgoing.push_back({gluon, {-192, -88, 54, 218}, {}}); ev.incoming[0] = {gluon, { 0, 0, -591, 591}, {}}; ev.incoming[1] = {gluon, { 0, 0, 433, 433}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { -42, 18, -49, 67}, {}}); ev.outgoing.push_back({gluon, { -12, -54, -28, 62}, {}}); ev.outgoing.push_back({higgs, { 99, 32, -16, 163}, {}}); ev.outgoing.push_back({gluon, { -45, 4, 72, 85}, {}}); ev.incoming[0] = {gluon, { 0, 0, -199, 199}, {}}; ev.incoming[1] = {gluon, { 0, 0, 178, 178}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { -65, -32, -76, 105}, {}}); ev.outgoing.push_back({gluon, { -22, 31, -34, 51}, {}}); ev.outgoing.push_back({gluon, { -12, -67, -36, 77}, {}}); ev.outgoing.push_back({higgs, { 99, 68, -4, 173}, {}}); ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}}; ev.incoming[1] = {gluon, { 0, 0, 128, 128}, {}}; return ev; default: ev.outgoing.push_back({gluon, { -90, -135, 30, 165}, {}}); ev.outgoing.push_back({gluon, {-108, 198, 76, 238}, {}}); ev.outgoing.push_back({gluon, { 198, -63, 126, 243}, {}}); ev.incoming[0] = {gluon, { 0, 0, -207, 207}, {}}; ev.incoming[1] = {gluon, { 0, 0, 439, 439}, {}}; return ev; } } if(njet == 4){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, {-200, -155, -64, 261}, {}}); ev.outgoing.push_back({gluon, { 198, 194, 57, 283}, {}}); ev.outgoing.push_back({gluon, { 1, 32, 8, 33}, {}}); ev.outgoing.push_back({gluon, {-198, -143, 186, 307}, {}}); ev.incoming[0] = {gluon, { 0, 0, -515, 515}, {}}; ev.incoming[1] = {gluon, { 0, 0, 626, 626}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { 198, 61, -162, 263}, {}}); ev.outgoing.push_back({higgs, { 199, 72, -76, 257}, {}}); ev.outgoing.push_back({gluon, {-200, 135, 144, 281}, {}}); ev.outgoing.push_back({gluon, {-198, -186, 171, 321}, {}}); ev.outgoing.push_back({gluon, { 1, -82, 122, 147}, {}}); ev.incoming[0] = {gluon, { 0, 0, -535, 535}, {}}; ev.incoming[1] = {gluon, { 0, 0, 734, 734}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, {-180, -27, -164, 245}, {}}); ev.outgoing.push_back({gluon, {-108, 78, -36, 138}, {}}); ev.outgoing.push_back({higgs, { 196, -189, 68, 307}, {}}); ev.outgoing.push_back({gluon, {-107, 136, 76, 189}, {}}); ev.outgoing.push_back({gluon, { 199, 2, 178, 267}, {}}); ev.incoming[0] = {gluon, { 0, 0, -512, 512}, {}}; ev.incoming[1] = {gluon, { 0, 0, 634, 634}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { -12, -30, -84, 90}, {}}); ev.outgoing.push_back({gluon, { -72, 22, -96, 122}, {}}); ev.outgoing.push_back({gluon, { 68, 0, -51, 85}, {}}); ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}}); ev.outgoing.push_back({gluon, { -48, -64, 84, 116}, {}}); ev.incoming[0] = {gluon, { 0, 0, -409, 409}, {}}; ev.incoming[1] = {gluon, { 0, 0, 181, 181}, {}}; return ev; case 4: ev.outgoing.push_back({gluon, { -72, -49, -72, 113}, {}}); ev.outgoing.push_back({gluon, { -48, 0, -36, 60}, {}}); ev.outgoing.push_back({gluon, { -12, 54, -36, 66}, {}}); ev.outgoing.push_back({gluon, { 68, -77, -56, 117}, {}}); ev.outgoing.push_back({higgs, { 64, 72, -81, 177}, {}}); ev.incoming[0] = {gluon, { 0, 0, -407, 407}, {}}; ev.incoming[1] = {gluon, { 0, 0, 126, 126}, {}}; return ev; default: ev.outgoing.push_back({gluon, { 248, -56, -122, 282}, {}}); ev.outgoing.push_back({gluon, { 249, 30, -10, 251}, {}}); ev.outgoing.push_back({gluon, {-249, -18, 26, 251}, {}}); ev.outgoing.push_back({gluon, {-248, 44, 199, 321}, {}}); ev.incoming[0] = {gluon, { 0, 0, -506, 506}, {}}; ev.incoming[1] = {gluon, { 0, 0, 599, 599}, {}}; return ev; } } if(njet == 6){ switch(pos_boson){ case 0: ev.outgoing.push_back({higgs, { 349, 330, -94, 505}, {}}); ev.outgoing.push_back({gluon, {-315, -300, 0, 435}, {}}); ev.outgoing.push_back({gluon, { 347, 306, 18, 463}, {}}); ev.outgoing.push_back({gluon, {-249, -342, 162, 453}, {}}); ev.outgoing.push_back({gluon, { 345, 312, 284, 545}, {}}); ev.outgoing.push_back({gluon, {-324, -126, 292, 454}, {}}); ev.outgoing.push_back({gluon, {-153, -180, 304, 385}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1137, 1137}, {}}; ev.incoming[1] = {gluon, { 0, 0, 2103, 2103}, {}}; return ev; case 1: ev.outgoing.push_back({gluon, { 242, 241, -182, 387}, {}}); ev.outgoing.push_back({higgs, { 243, 238, -190, 409}, {}}); ev.outgoing.push_back({gluon, {-218, -215, -74, 315}, {}}); ev.outgoing.push_back({gluon, {-224, -224, 112, 336}, {}}); ev.outgoing.push_back({gluon, { 241, 182, 154, 339}, {}}); ev.outgoing.push_back({gluon, { -53, -234, 126, 271}, {}}); ev.outgoing.push_back({gluon, {-231, 12, 156, 279}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1117, 1117}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1219, 1219}, {}}; return ev; case 2: ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}}); ev.outgoing.push_back({gluon, { -86, -46, -17, 99}, {}}); ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}}); ev.outgoing.push_back({gluon, { -60, -135, 64, 161}, {}}); ev.outgoing.push_back({gluon, { 150, 123, 110, 223}, {}}); ev.outgoing.push_back({gluon, {-154, -49, 98, 189}, {}}); ev.outgoing.push_back({gluon, {-153, -148, 144, 257}, {}}); ev.incoming[0] = {gluon, { 0, 0, -504, 504}, {}}; ev.incoming[1] = {gluon, { 0, 0, 861, 861}, {}}; return ev; case 3: ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}}); ev.outgoing.push_back({gluon, {-198, -189, -54, 279}, {}}); ev.outgoing.push_back({gluon, {-200, -64, 2, 210}, {}}); ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}}); ev.outgoing.push_back({gluon, {-199, -184, 172, 321}, {}}); ev.outgoing.push_back({gluon, { 196, 168, 177, 313}, {}}); ev.outgoing.push_back({gluon, { 4, -86, 92, 126}, {}}); ev.incoming[0] = {gluon, { 0, 0, -745, 745}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1074, 1074}, {}}; return ev; case 4: ev.outgoing.push_back({gluon, { 151, 102, -42, 187}, {}}); ev.outgoing.push_back({gluon, { -86, -133, -14, 159}, {}}); ev.outgoing.push_back({gluon, {-154, -104, -8, 186}, {}}); ev.outgoing.push_back({gluon, { -60, 11, 0, 61}, {}}); ev.outgoing.push_back({higgs, { 152, 153, 0, 249}, {}}); ev.outgoing.push_back({gluon, { 150, 125, 90, 215}, {}}); ev.outgoing.push_back({gluon, {-153, -154, 126, 251}, {}}); ev.incoming[0] = {gluon, { 0, 0, -578, 578}, {}}; ev.incoming[1] = {gluon, { 0, 0, 730, 730}, {}}; return ev; case 5: ev.outgoing.push_back({gluon, { -15, -90, -94, 131}, {}}); ev.outgoing.push_back({gluon, { -11, 82, -74, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -80, -64, 105}, {}}); ev.outgoing.push_back({gluon, { -48, -25, -36, 65}, {}}); ev.outgoing.push_back({gluon, { -12, 99, -16, 101}, {}}); ev.outgoing.push_back({higgs, { 68, 92, -18, 170}, {}}); ev.outgoing.push_back({gluon, { -5, -78, 54, 95}, {}}); ev.incoming[0] = {gluon, { 0, 0, -513, 513}, {}}; ev.incoming[1] = {gluon, { 0, 0, 265, 265}, {}}; return ev; case 6: ev.outgoing.push_back({gluon, { 198, 197, -66, 287}, {}}); ev.outgoing.push_back({gluon, { 4, -84, -18, 86}, {}}); ev.outgoing.push_back({gluon, {-198, -60, -36, 210}, {}}); ev.outgoing.push_back({gluon, { 196, -78, -36, 214}, {}}); ev.outgoing.push_back({gluon, {-200, 45, 0, 205}, {}}); ev.outgoing.push_back({gluon, {-199, -178, 2, 267}, {}}); ev.outgoing.push_back({higgs, { 199, 158, 6, 283}, {}}); ev.incoming[0] = {gluon, { 0, 0, -850, 850}, {}}; ev.incoming[1] = {gluon, { 0, 0, 702, 702}, {}}; return ev; default: ev.outgoing.push_back({gluon, {-350, -112, -280, 462}, {}}); ev.outgoing.push_back({gluon, { 347, 266, -322, 543}, {}}); ev.outgoing.push_back({gluon, {-349, -314, -38, 471}, {}}); ev.outgoing.push_back({gluon, { 349, 348, 12, 493}, {}}); ev.outgoing.push_back({gluon, {-342, -54, 23, 347}, {}}); ev.outgoing.push_back({gluon, { 345, -134, 138, 395}, {}}); ev.incoming[0] = {gluon, { 0, 0, -1589, 1589}, {}}; ev.incoming[1] = {gluon, { 0, 0, 1122, 1122}, {}}; return ev; } } if(njet == 7){ switch(pos_boson){ case -1: // jet idx: -1 0 1 2 3 4 5 ev.outgoing.push_back({gluon, { -15, -18, -54, 59}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { 23, -100, -64, 121}, {}}); ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -92, -12, 93}, {}}); ev.outgoing.push_back({gluon, { -48, -76, -2, 90}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -439, 439}, {}}; ev.incoming[1] = {gluon, { 0, 0, 273, 273}, {}}; return ev; case -2: // jet idx: 0 1 2 3 4 -1 -1 ev.outgoing.push_back({gluon, { -5, -86, -82, 119}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -14, 20, 54}, {}}); ev.outgoing.push_back({gluon, { 23, -50, 26, 61}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -15, -18, 54, 59}, {}}); ev.outgoing.push_back({gluon, { -11, -20, 88, 91}, {}}); ev.incoming[0] = {gluon, { 0, 0, -215, 215}, {}}; ev.incoming[1] = {gluon, { 0, 0, 397, 397}, {}}; return ev; case -3: // jet idx: 0 0 1 2 2 3 4 // jet pt fraction: 0.6 0.38 1 0.49 0.51 1 1 ev.outgoing.push_back({gluon, { 23, -94, -62, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -5, -62, -34, 71}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 1.1e+02}, {}}); ev.outgoing.push_back({gluon, { -11, 98, 70, 1.2e+02}, {}}); ev.outgoing.push_back({gluon, { -48, -1e+02, 82, 1.4e+02}, {}}); ev.outgoing.push_back({gluon, { -15, -30, 78, 85}, {}}); ev.incoming[0] = {gluon, { 0, 0, -2.7e+02, 2.7e+02}, {}}; ev.incoming[1] = {gluon, { 0, 0, 4.8e+02, 4.8e+02}, {}}; return ev; case -4: // jet idx: 0 1 1 2 3 4 4 // jet pt fraction: 1 0.51 0.49 1 1 0.25 0.75 ev.outgoing.push_back({gluon, { -5, -88, -64, 109}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { 23, -70, 22, 77}, {}}); ev.outgoing.push_back({gluon, { -15, -32, 16, 39}, {}}); ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}}); ev.incoming[0] = {gluon, { 0, 0, -381, 381}, {}}; ev.incoming[1] = {gluon, { 0, 0, 324, 324}, {}}; return ev; case -5: // jet idx: 0 1 -1 -1 2 3 4 ev.outgoing.push_back({gluon, { -15, -26, -62, 69}, {}}); ev.outgoing.push_back({gluon, { -48, -60, -54, 94}, {}}); ev.outgoing.push_back({gluon, { 23, 10, -14, 29}, {}}); ev.outgoing.push_back({gluon, { -5, -20, 0, 21}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -11, -92, 40, 101}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -278, 278}, {}}; ev.incoming[1] = {gluon, { 0, 0, 264, 264}, {}}; return ev; case -6: // jet idx: 0 1 1 2 -1 2 3 // jet pt fraction: 1 0.63 0.36 0.49 1 0.51 1 ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -100, 26, 114}, {}}); ev.outgoing.push_back({gluon, { -15, -62, 26, 69}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}}); ev.outgoing.push_back({gluon, { -11, 98, 70, 121}, {}}); ev.outgoing.push_back({gluon, { 23, -96, 92, 135}, {}}); ev.incoming[0] = {gluon, { 0, 0, -216, 216}, {}}; ev.incoming[1] = {gluon, { 0, 0, 486, 486}, {}}; return ev; case -7: // jet idx: 0 1 2 2 3 3 4 // jet pt fraction: 1 1 0.51 0.49 0.18 0.82 1 ev.outgoing.push_back({gluon, { -15, -80, -100, 129}, {}}); ev.outgoing.push_back({gluon, { 23, -96, -92, 135}, {}}); ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -22, -10, 25}, {}}); ev.outgoing.push_back({gluon, { -48, -88, -31, 105}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.incoming[0] = {gluon, { 0, 0, -541, 541}, {}}; ev.incoming[1] = {gluon, { 0, 0, 202, 202}, {}}; return ev; case -8: // jet idx: 0 1 2 2 2 3 4 // jet pt fraction: 1 1 0.21 0.37 0.41 1 1 ev.outgoing.push_back({gluon, { -48, -44, -62, 90}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -50, -22, 55}, {}}); ev.outgoing.push_back({gluon, { 23, -90, -34, 99}, {}}); ev.outgoing.push_back({gluon, { -15, -100, -28, 105}, {}}); ev.outgoing.push_back({gluon, { 68, 93, -20, 117}, {}}); ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}}); ev.incoming[0] = {gluon, { 0, 0, -423, 423}, {}}; ev.incoming[1] = {gluon, { 0, 0, 277, 277}, {}}; return ev; case -9: // jet idx: 0 1 2 1 3 0 4 // jet pt fraction: 0.72 0.51 1 0.49 1 0.28 1 ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -76, -40, 89}, {}}); ev.outgoing.push_back({gluon, { -11, 92, -40, 101}, {}}); ev.outgoing.push_back({gluon, { -48, -68, -34, 90}, {}}); ev.outgoing.push_back({gluon, { -5, -38, -14, 41}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.incoming[0] = {gluon, { 0, 0, -446, 446}, {}}; ev.incoming[1] = {gluon, { 0, 0, 220, 220}, {}}; return ev; case -10: // jet idx: 0 1 3 2 4 3 1 // jet pt fraction: 1 0.33 0.51 1 1 0.49 0.67 ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -48, 16, 51}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { 23, -76, 52, 95}, {}}); ev.outgoing.push_back({gluon, { -48, -60, 54, 94}, {}}); ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}}); ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}}); ev.incoming[0] = {gluon, { 0, 0, -183, 183}, {}}; ev.incoming[1] = {gluon, { 0, 0, 521, 521}, {}}; return ev; case -11: // jet idx: 0 1 2 3 4 -1 5 // jet pt fraction: 1 1 1 1 1 1 1 ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { 23, -90, -2, 93}, {}}); ev.outgoing.push_back({gluon, { -48, -76, 2, 90}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -5, -22, 10, 25}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.incoming[0] = {gluon, { 0, 0, -360, 360}, {}}; ev.incoming[1] = {gluon, { 0, 0, 314, 314}, {}}; return ev; case -12: // jet idx: 0 1 -1 2 3 4 3 // jet pt fraction: 1 1 1 1 0.35 1 0.65 ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 4, 29}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -15, -58, 34, 69}, {}}); ev.outgoing.push_back({gluon, { -11, 92, 68, 115}, {}}); ev.outgoing.push_back({gluon, { -48, -100, 82, 138}, {}}); ev.incoming[0] = {gluon, { 0, 0, -302, 302}, {}}; ev.incoming[1] = {gluon, { 0, 0, 392, 392}, {}}; return ev; case -13: // jet idx: 0 1 2 3 3 4 2 // jet pt fraction: 1 1 0.5 0.35 0.65 1 0.5 ev.outgoing.push_back({gluon, { -15, -98, -62, 117}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -5, -28, 20, 35}, {}}); ev.outgoing.push_back({gluon, { -48, -96, 75, 131}, {}}); ev.outgoing.push_back({gluon, { 23, -62, 50, 83}, {}}); ev.outgoing.push_back({gluon, { -11, 96, 76, 123}, {}}); ev.incoming[0] = {gluon, { 0, 0, -241, 241}, {}}; ev.incoming[1] = {gluon, { 0, 0, 476, 476}, {}}; return ev; case -14: // jet idx: 0 1 2 3 3 4 2 // jet pt fraction: 1 1 0.52 0.35 0.65 1 0.48 ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -12, 95, 56, 111}, {}}); ev.outgoing.push_back({gluon, { -15, -96, 72, 121}, {}}); ev.outgoing.push_back({gluon, { -5, -42, 38, 57}, {}}); ev.outgoing.push_back({gluon, { -48, -44, 62, 90}, {}}); ev.outgoing.push_back({gluon, { -11, 88, 88, 125}, {}}); ev.incoming[0] = {gluon, { 0, 0, -231, 231}, {}}; ev.incoming[1] = {gluon, { 0, 0, 505, 505}, {}}; return ev; case -15: // jet idx: 0 -1 1 0 2 3 4 // jet pt fraction: 0.51 1 1 0.49 1 1 1 ev.outgoing.push_back({gluon, { -11, 98, -70, 121}, {}}); ev.outgoing.push_back({gluon, { -5, -16, -12, 21}, {}}); ev.outgoing.push_back({gluon, { 23, -94, -62, 115}, {}}); ev.outgoing.push_back({gluon, { -12, 95, -56, 111}, {}}); ev.outgoing.push_back({gluon, { 68, 93, 20, 117}, {}}); ev.outgoing.push_back({gluon, { -48, -76, 70, 114}, {}}); ev.outgoing.push_back({gluon, { -15, -100, 80, 129}, {}}); ev.incoming[0] = {gluon, { 0, 0, -379, 379}, {}}; ev.incoming[1] = {gluon, { 0, 0, 349, 349}, {}}; return ev; } } throw HEJ::unknown_option{"unknown process"}; } HEJ::Event::EventData parse_configuration( std::array const & in, std::vector const & out, int const overwrite_boson ){ auto boson = std::find_if(out.cbegin(), out.cend(), [](std::string id){ return !HEJ::is_parton(HEJ::to_ParticleID(id)); }); int const pos_boson = (overwrite_boson!=0)?overwrite_boson: ((boson == out.cend())?-1:std::distance(out.cbegin(), boson)); std::size_t njets = out.size(); if( (overwrite_boson == 0) && boson != out.cend()) --njets; HEJ::Event::EventData ev{get_process(njets, pos_boson)}; ASSERT((pos_boson<0) || (ev.outgoing[pos_boson].type == HEJ::ParticleID::higgs)); for(std::size_t i=0; i const & in, std::vector const & out, bool reconstruct /* = false */, std::unordered_map< size_t, std::vector > decays /* = {} */ ) { using namespace HEJ::pid; HEJ::Event::EventData evd; const size_t n_out {out.size()}; // Parameters double pT {40.}; std::array cs_phi {1, 0}; std::array cs_dphi {cos(2.*M_PI/n_out), sin(2.*M_PI/n_out)}; double y {-4.5}, dy {std::max(-2.*y/(n_out-1), 0.5)}; double sum_pE {0}, sum_pz {0}; // Outgoing evd.outgoing.reserve(n_out); for(size_t i = 0; i < n_out; ++i) { auto pid = HEJ::to_ParticleID(out[i]); // Kinematics double mT2 = pT*pT; if(HEJ::is_AWZH_boson(pid)) { auto const ref = ew_parameters.prop(pid); mT2 += ref.mass*ref.mass; } double pz {sqrt(mT2)*sinh(y)}, pE {sqrt(mT2)*cosh(y)}; evd.outgoing.push_back({pid, {pT*cs_phi[0], pT*cs_phi[1], pz, pE}, {}}); sum_pE+=pE; sum_pz+=pz; // Update cos,sin, y cs_phi = {cs_phi[0] * cs_dphi[0] - cs_phi[1] * cs_dphi[1], cs_phi[1] * cs_dphi[0] + cs_phi[0] * cs_dphi[1]}; y += dy; } // Specified decays evd.decays.reserve(decays.size()); for(auto const & decay : decays) { auto const parent = evd.outgoing.at(decay.first); std::vector progeny = decay_kinematics(parent); for(std::size_t i = 0; i < decay.second.size(); ++i) { progeny[i].type = HEJ::to_ParticleID(decay.second[i]); } evd.decays.emplace(decay.first, std::move(progeny)); } // Reconstruct decays if(reconstruct) { evd.reconstruct_intermediate(ew_parameters); } // Incoming double p0 {0.5*(sum_pE + sum_pz)}; evd.incoming[0] = {HEJ::to_ParticleID(in[0]),{0,0,-p0,p0},{}}; double p1 {0.5*(sum_pE - sum_pz)}; evd.incoming[1] = {HEJ::to_ParticleID(in[1]),{0,0,+p1,p1},{}}; return evd; } namespace { static std::mt19937_64 ran{0}; } void shuffle_particles(HEJ::Event::EventData & ev) { // incoming std::shuffle(ev.incoming.begin(), ev.incoming.end(), ran); // outgoing (through index) auto old_outgoing = std::move(ev).outgoing; std::vector idx(old_outgoing.size()); std::iota(idx.begin(), idx.end(), 0); std::shuffle(idx.begin(), idx.end(), ran); ev.outgoing.clear(); ev.outgoing.reserve(old_outgoing.size()); for(std::size_t i: idx) { ev.outgoing.emplace_back(std::move(old_outgoing[i])); } // find decays again if(!ev.decays.empty()){ auto old_decays = std::move(ev).decays; ev.decays.clear(); for(std::size_t i=0; isecond)); } for(auto & decay: ev.decays){ std::shuffle(decay.second.begin(), decay.second.end(), ran); } } } +bool couple_quark(std::string const & boson, std::string & quark){ + if(std::abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){ + auto qflav{ HEJ::to_ParticleID(quark) }; + if(!HEJ::is_anyquark(qflav)) return false; + const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1; + if(W_charge*qflav < 0 && !(std::abs(qflav)%2)) return false; // not anti-down + if(W_charge*qflav > 0 && (std::abs(qflav)%2)) return false; // not up + quark=std::to_string(qflav-W_charge); + } + if(HEJ::to_ParticleID(boson) == HEJ::ParticleID::Z_photon_mix){ + auto qflav{ HEJ::to_ParticleID(quark) }; + if(!HEJ::is_anyquark(qflav)) return false; + } + return true; +} + std::vector decay_kinematics( HEJ::Particle const & parent ) { std::vector decay_products(2); const double E = parent.m()/2; const double theta = 2.*M_PI*ran()/static_cast(ran.max()); const double cos_phi = 2.*ran()/static_cast(ran.max())-1.; const double sin_phi = std::sqrt(1. - cos_phi*cos_phi); // Know 0 < phi < pi const double px = E*std::cos(theta)*sin_phi; const double py = E*std::sin(theta)*sin_phi; const double pz = E*cos_phi; decay_products[0].p.reset(px, py, pz, E); decay_products[1].p.reset(-px, -py, -pz, E); for(auto & particle: decay_products) particle.p.boost(parent.p); return decay_products; } std::vector decay_W( HEJ::Particle const & parent ){ if(parent.m() == 0.) // we can't decay massless partons return {}; std::array decays; if(parent.type==HEJ::ParticleID::Wp){ // order matters: first particle, second anti decays[0] = HEJ::ParticleID::nu_e; decays[1] = HEJ::ParticleID::e_bar; } else { // this function is for testing: we don't check that parent==W boson decays[0] = HEJ::ParticleID::e; decays[1] = HEJ::ParticleID::nu_e_bar; } std::vector decay_products = decay_kinematics(parent); for(std::size_t i = 0; i < decay_products.size(); ++i){ decay_products[i].type = decays[i]; } return decay_products; } std::vector decay_Z( HEJ::Particle const & parent ){ if(parent.m() == 0.) // we can't decay massless partons return {}; std::array decays; // order matters: first particle, second anti decays[0] = HEJ::ParticleID::electron; decays[1] = HEJ::ParticleID::positron; std::vector decay_products = decay_kinematics(parent); for(std::size_t i = 0; i < decay_products.size(); ++i){ decay_products[i].type = decays[i]; } return decay_products; } diff --git a/t/hej_test.hh b/t/hej_test.hh index 0270050..bd53558 100644 --- a/t/hej_test.hh +++ b/t/hej_test.hh @@ -1,76 +1,79 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #pragma once #include #include #include #include "HEJ/Event.hh" #include "HEJ/exceptions.hh" namespace HEJ { struct Particle; } //! throw error if condition not fulfilled #define ASSERT(x) if(!(x)) { \ throw std::logic_error("Assertion '" #x "' failed."); \ } //! throw error if prop is different between ev1 and ev2 #define ASSERT_PROPERTY(ev1,ev2,prop) ASSERT(ev1.prop == ev2.prop) //! throw error if condition not fulfilled #define ASSERT_THROW(x, exception) try { \ x; \ std::cerr << "'" #x "' did not throw an exception.\n"; \ throw; \ } catch(exception const &){} \ catch (...) { \ std::cerr << "Unexpected exception thrown for '" #x "'.\n"; \ throw; \ } /** @brief get specific Phase Space Points for njets with boson at pos_boson * * if pos_boson = -1 (or not implemented) -> no boson * * njet==7 is special: has less jets, i.e. multiple parton in one jet, * all partons are massive (4 GeV) -> can be boson/decay * pos_boson < 0 to select process (see list for details) */ HEJ::Event::EventData get_process(int njet, int pos_boson); //! select process from string input (see also get_process) //! //! overwrite_boson to force a specific boson position, indepentent from input //! (useful for njet == 7) HEJ::Event::EventData parse_configuration( std::array const & in, std::vector const & out, int overwrite_boson = 0 ); /** @brief generate a rapidity ordered event */ HEJ::Event::EventData rapidity_order_ps( std::array const & in, std::vector const & out, bool reconstruct = false, std::unordered_map > decays = {} ); //! shuffle particles around void shuffle_particles(HEJ::Event::EventData & ev); +//! Helper function to couple quarksfor flavour-changing bosons +bool couple_quark(std::string const & boson, std::string & quark); + //! Decay kinematics for 1->2 std::vector decay_kinematics( HEJ::Particle const & parent ); //! Decay W boson to lepton & neutrino std::vector decay_W( HEJ::Particle const & parent ); //! Decay Z to electron-positron std::vector decay_Z( HEJ::Particle const & parent ); diff --git a/t/test_classify.cc b/t/test_classify.cc index bfc252a..117fe6d 100644 --- a/t/test_classify.cc +++ b/t/test_classify.cc @@ -1,596 +1,580 @@ /** * \authors The HEJ collaboration (see AUTHORS for details) * \date 2019-2020 * \copyright GPLv2 or later */ #include "hej_test.hh" #include #include #include #include #include #include #include "fastjet/JetDefinition.hh" #include "HEJ/Event.hh" #include "HEJ/event_types.hh" #include "HEJ/exceptions.hh" #include "HEJ/PDG_codes.hh" namespace { const fastjet::JetDefinition jet_def{fastjet::JetAlgorithm::antikt_algorithm, 0.4}; const double min_jet_pt{30.}; const std::vector all_quarks{"-4","-1","1","2","3","4"}; const std::vector all_partons{"g","-2","-1","1","2","3","4"}; const std::vector all_bosons{"h", "Wp", "Wm", "Z_photon_mix"}; const std::vector all_gZ{"photon", "Z"}; const std::vector all_w{"W+", "W-"}; static std::mt19937_64 ran{0}; - bool couple_quark(std::string const & boson, std::string & quark){ - if(std::abs(HEJ::to_ParticleID(boson)) == HEJ::ParticleID::Wp){ - auto qflav{ HEJ::to_ParticleID(quark) }; - if(!HEJ::is_anyquark(qflav)) return false; - const int W_charge = HEJ::to_ParticleID(boson)>0?1:-1; - if(W_charge*qflav < 0 && !(std::abs(qflav)%2)) return false; // not anti-down - if(W_charge*qflav > 0 && (std::abs(qflav)%2)) return false; // not up - quark=std::to_string(qflav-W_charge); - } - if(HEJ::to_ParticleID(boson) == HEJ::ParticleID::Z_photon_mix){ - auto qflav{ HEJ::to_ParticleID(quark) }; - if(!HEJ::is_anyquark(qflav)) return false; - } - return true; - } - //! Upon clustering does EventData match the expected classification bool match_expectation(HEJ::event_type::EventType expected, HEJ::Event::EventData evd ){ HEJ::Event ev { evd.cluster(jet_def, min_jet_pt) }; if(ev.type() != expected){ std::cerr << "Expected type " << name(expected) << " but found " << name(ev.type()) << "\n" << ev; auto jet_idx{ ev.particle_jet_indices() }; std::cout << "Particle Jet indices: "; for(int const i: jet_idx) std::cout << i << " "; std::cout << std::endl; return false; } return true; } //! Does match_expectation for EventData from parse_configuration bool match_expectation( HEJ::event_type::EventType expected, std::array const & in, std::vector const & out, int const overwrite_boson = 0 ){ return match_expectation(expected, parse_configuration( in,out,overwrite_boson )); } //! test FKL configurations //! if implemented==false : check processes that are not in HEJ yet bool check_fkl( bool const implemented=true ){ using namespace HEJ; auto const type{ implemented?event_type::FKL:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = all_bosons; else { bosons = all_gZ; } for(std::string const & first: all_partons) // all quark flavours for(std::string const & last: all_partons){ for(int njet=2; njet<=6; ++njet){ // all multiplicities if(njet==5) continue; std::array base_in{first,last}; std::vector base_out(njet, "g"); base_out.front() = first; base_out.back() = last; if(implemented && !match_expectation(type, base_in, base_out)) return false; for(auto const & boson: bosons) // any boson for(int pos=0; pos<=njet; ++pos){ // at any position auto in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(ran); if(!couple_quark(boson, couple_idx?out.back():out.front())) continue; out.insert(out.begin()+pos, boson); if(!match_expectation(type, in, out)) return false; } } } return true; } //! test unordered configurations //! if implemented==false : check processes that are not in HEJ yet bool check_uno( bool const implemented=true ){ using namespace HEJ; auto const b{ implemented?event_type::unob:event_type::non_resummable }; auto const f{ implemented?event_type::unof:event_type::non_resummable }; std::vector bosons; if(implemented) { bosons = all_bosons; } else { bosons = all_gZ; } for(std::string const & uno: all_quarks) // all quark flavours for(std::string const & fkl: all_partons){ for(int njet=3; njet<=6; ++njet){ // all multiplicities >2 if(njet==5) continue; for(int i=0; i<2; ++i){ // forward & backwards std::array base_in; std::vector base_out(njet, "g"); const int uno_pos = i?1:(njet-2); const int fkl_pos = i?(njet-1):0; base_in[i?0:1] = uno; base_in[i?1:0] = fkl; base_out[uno_pos] = uno; base_out[fkl_pos] = fkl; auto expectation{ i?b:f }; if( implemented && !match_expectation(expectation, base_in, base_out) ) return false; for(auto const & boson: bosons){ // any boson // at any position (higgs only inside FKL chain) int start = 0; int end = njet; if(to_ParticleID(boson) == pid::higgs){ start = i?(uno_pos+1):fkl_pos; end = i?(fkl_pos+1):uno_pos; } for(int pos=start; pos<=end; ++pos){ auto in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(ran); if(!couple_quark(boson, couple_idx?out[fkl_pos]:out[uno_pos])) continue; out.insert(out.begin()+pos, boson); if(!match_expectation(expectation, in, out)) return false; } } } } } return true; } //! test extremal qqx configurations //! if implemented==false : check processes that are not in HEJ yet bool check_extremal_qqx( bool const implemented=true ){ using namespace HEJ; auto const b{ implemented?event_type::qqxexb:event_type::non_resummable }; auto const f{ implemented?event_type::qqxexf:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = all_w; else { bosons = all_gZ; bosons.push_back("h"); bosons.push_back("Z_photon_mix"); } for(std::string const & qqx: all_quarks) // all quark flavours for(std::string const & fkl: all_partons){ std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) }; for(int njet=3; njet<=6; ++njet){ // all multiplicities >2 if(njet==5) continue; for(int i=0; i<2; ++i){ // forward & backwards std::array base_in; std::vector base_out(njet, "g"); const int qqx_pos = i?0:(njet-2); const int fkl_pos = i?(njet-1):0; base_in[i?0:1] = "g"; base_in[i?1:0] = fkl; base_out[fkl_pos] = fkl; base_out[qqx_pos] = qqx; base_out[qqx_pos+1] = qqx2; auto expectation{ i?b:f }; if( implemented && !match_expectation(expectation, base_in, base_out) ) return false; for(auto const & boson: bosons){ // all bosons // at any position (higgs only inside FKL chain) int start = 0; int end = njet; if(to_ParticleID(boson) == pid::higgs){ start = i?(qqx_pos+2):fkl_pos; end = i?(fkl_pos+1):qqx_pos; } for(int pos=start; pos<=end; ++pos){ auto in{base_in}; auto out{base_out}; // change quark flavours for W const bool couple_idx = std::uniform_int_distribution{0,1}(ran); if(couple_idx || !couple_quark(boson, out[fkl_pos]) ){ // (randomly) try couple to FKL, else fall-back to qqx if(!couple_quark(boson, out[qqx_pos])) couple_quark(boson, out[qqx_pos+1]); } out.insert(out.begin()+pos, boson); if(!match_expectation(expectation, in, out)) return false; } } } } // test allowed jet configurations if( implemented){ if( !( match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -3) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -4) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -5) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -5) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -6) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -7) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -7) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -8) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -8) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -9) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -10) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -11) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -11) && match_expectation(f,{fkl,"g"},{fkl,"g","g","g","g",qqx,qqx2}, -12) && match_expectation(b,{"g",fkl},{qqx,qqx2,"g","g","g","g",fkl}, -12) )) return false; if (fkl == "2") { if( !( match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -3) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -4) && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -5) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -5) && match_expectation(f,{"2","g"},{"1","g","Wp","g","g",qqx,qqx2}, -6) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -7) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -7) && match_expectation(f,{"2","g"},{"1","Wp","g","g","g",qqx,qqx2}, -8) && match_expectation(b,{"g","2"},{qqx,qqx2,"Wp","g","g","g","1"}, -8) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -9) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -10) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -11) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","g","g","Wp","1"}, -11) && match_expectation(f,{"2","g"},{"1","g","g","g","Wp",qqx,qqx2}, -12) && match_expectation(b,{"g","2"},{qqx,qqx2,"g","Wp","g","g","1"}, -12) )) return false; } } } return true; } //! test central qqx configurations //! if implemented==false : check processes that are not in HEJ yet bool check_central_qqx(bool const implemented=true){ using namespace HEJ; auto const t{ implemented?event_type::qqxmid:event_type::non_resummable }; std::vector bosons; if(implemented) bosons = all_w; else { bosons = all_gZ; bosons.push_back("h"); bosons.push_back("Z_photon_mix"); } for(std::string const & qqx: all_quarks) // all quark flavours for(std::string const & fkl1: all_partons) for(std::string const & fkl2: all_partons){ std::string const qqx2{ std::to_string(HEJ::to_ParticleID(qqx)*-1) }; for(int njet=4; njet<=6; ++njet){ // all multiplicities >3 if(njet==5) continue; for(int qqx_pos=1; qqx_pos base_in; std::vector base_out(njet, "g"); base_in[0] = fkl1; base_in[1] = fkl2; base_out.front() = fkl1; base_out.back() = fkl2; base_out[qqx_pos] = qqx; base_out[qqx_pos+1] = qqx2; if( implemented && !match_expectation(t, base_in, base_out) ) return false; for(auto const & boson: bosons) // any boson for(int pos=0; pos<=njet; ++pos){ // at any position if( to_ParticleID(boson) == pid::higgs && (pos==qqx_pos || pos==qqx_pos+1) ) continue; auto in{base_in}; auto out{base_out}; // change quark flavours for W const int couple_idx{ std::uniform_int_distribution{0,2}(ran) }; // (randomly) try couple to FKL, else fall-back to qqx if( couple_idx == 0 && couple_quark(boson, out.front()) ){} else if( couple_idx == 1 && couple_quark(boson, out.back()) ){} else { if(!couple_quark(boson, out[qqx_pos])) couple_quark(boson, out[qqx_pos+1]); } out.insert(out.begin()+pos, boson); if(!match_expectation(t, in, out)) return false; } } } } return true; } // this checks a (non excessive) list of non-resummable states bool check_non_resummable(){ auto type{ HEJ::event_type::non_resummable}; return // 2j - crossing lines match_expectation(type, {"g","2"}, {"2","g"}) && match_expectation(type, {"-1","g"}, {"g","-1"}) && match_expectation(type, {"1","-1"}, {"-1","1"}) && match_expectation(type, {"g","2"}, {"2","g","h"}) && match_expectation(type, {"1","2"}, {"2","h","1"}) && match_expectation(type, {"1","-1"}, {"h","-1","1"}) && match_expectation(type, {"g","2"}, {"Wp","1","g"}) && match_expectation(type, {"1","-1"}, {"-2","Wp","1"}) && match_expectation(type, {"4","g"}, {"g","3","Wp"}) && match_expectation(type, {"1","-2"}, {"-1","Wm","1"}) && match_expectation(type, {"g","3"}, {"4","g","Wm"}) && match_expectation(type, {"1","3"}, {"Wm","4","1"}) && match_expectation(type, {"g","2"}, {"Z_photon_mix","2","g"}) && match_expectation(type, {"1","-1"}, {"-1","Z_photon_mix","1"}) && match_expectation(type, {"4","g"}, {"g","4","Z_photon_mix"}) // 2j - qqx && match_expectation(type, {"g","g"}, {"1","-1"}) && match_expectation(type, {"g","g"}, {"-2","2","h"}) && match_expectation(type, {"g","g"}, {"-4","Wp","3"}) && match_expectation(type, {"g","g"}, {"Wm","-1","2"}) && match_expectation(type, {"g","g"}, {"-3","Z_photon_mix","3"}) // 3j - crossing lines && match_expectation(type, {"g","4"}, {"4","g","g"}) && match_expectation(type, {"-1","g"}, {"g","g","-1"}) && match_expectation(type, {"1","3"}, {"3","g","1"}) && match_expectation(type, {"-2","2"}, {"2","g","-2","h"}) && match_expectation(type, {"-3","g"}, {"g","g","Wp","-4"}) && match_expectation(type, {"1","-2"}, {"Wm","-1","g","1"}) && match_expectation(type, {"-1","g"}, {"1","-1","-1"}) && match_expectation(type, {"1","-4"}, {"Z_photon_mix","-4","g","1"}) // higgs inside uno && match_expectation(type, {"-1","g"}, {"g","h","-1","g"}) && match_expectation(type, {"-1","1"}, {"g","h","-1","1"}) && match_expectation(type, {"g","2"}, {"g","2","h","g"}) && match_expectation(type, {"-1","1"}, {"-1","1","h","g"}) // higgs outside uno && match_expectation(type, {"-1","g"}, {"h","g","-1","g"}) && match_expectation(type, {"-1","1"}, {"-1","1","g","h"}) // higgs inside qqx && match_expectation(type, {"g","g"}, {"-1","h","1","g","g"}) && match_expectation(type, {"g","g"}, {"g","-1","h","1","g"}) && match_expectation(type, {"g","g"}, {"g","g","2","h","-2"}) // higgs outside qqx && match_expectation(type, {"g","g"}, {"h","-1","1","g","g"}) && match_expectation(type, {"g","g"}, {"g","g","2","-2","h"}) // 4j - two uno && match_expectation(type, {"-2","2"}, {"g","-2","2","g"}) && match_expectation(type, {"1","3"}, {"g","1","h","3","g"}) && match_expectation(type, {"1","2"}, {"g","1","3","Wp","g"}) && match_expectation(type, {"1","-2"}, {"g","Wm","1","-1","g"}) && match_expectation(type, {"3","2"}, {"g","3","Z_photon_mix","2","g"}) // 4j - two gluon outside && match_expectation(type, {"g","4"}, {"g","4","g","g"}) && match_expectation(type, {"1","3"}, {"1","3","h","g","g"}) && match_expectation(type, {"1","2"}, {"1","3","g","Wp","g"}) && match_expectation(type, {"1","-2"}, {"1","Wm","-1","g","g"}) && match_expectation(type, {"-1","g"}, {"g","g","-1","g"}) && match_expectation(type, {"1","3"}, {"g","g","1","3","h"}) && match_expectation(type, {"1","2"}, {"g","g","1","Wp","3"}) && match_expectation(type, {"1","-2"}, {"Wm","g","g","1","-1"}) && match_expectation(type, {"-1","2"}, {"g","g","-1","Z_photon_mix","2"}) // 4j - ggx+uno && match_expectation(type, {"g","4"}, {"1","-1","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","-3","3"}) && match_expectation(type, {"g","4"}, {"1","-1","h","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","-3","3","h"}) && match_expectation(type, {"g","4"}, {"Wp","1","-1","3","g"}) && match_expectation(type, {"2","g"}, {"g","2","-4","Wp","3"}) && match_expectation(type, {"g","4"}, {"2","Wm","-1","4","g"}) && match_expectation(type, {"2","g"}, {"g","2","Wp","-3","4"}) && match_expectation(type, {"-4","g"}, {"g","-4","-3","3","Z_photon_mix"}) // 3j - crossing+uno && match_expectation(type, {"1","4"}, {"g","4","1"}) && match_expectation(type, {"1","4"}, {"4","1","g"}) && match_expectation(type, {"1","4"}, {"g","h","4","1"}) && match_expectation(type, {"-1","-3"},{"Wm","g","-4","-1"}) && match_expectation(type, {"1","4"}, {"3","1","Wp","g"}) && match_expectation(type, {"1","4"}, {"3","1","g","h"}) && match_expectation(type, {"2","3"}, {"3","2","Z_photon_mix","g"}) // 3j - crossing+qqx && match_expectation(type, {"1","g"}, {"-1","1","g","1"}) && match_expectation(type, {"1","g"}, {"-1","1","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","1","-1"}) && match_expectation(type, {"g","1"}, {"g","1","1","-1"}) && match_expectation(type, {"1","g"}, {"2","-2","g","1"}) && match_expectation(type, {"1","g"}, {"2","-2","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","-2","2"}) && match_expectation(type, {"g","1"}, {"g","1","-2","2"}) && match_expectation(type, {"1","g"}, {"-1","1","h","g","1"}) && match_expectation(type, {"1","g"}, {"-1","h","1","1","g"}) && match_expectation(type, {"g","1"}, {"1","g","1","h","-1"}) && match_expectation(type, {"g","1"}, {"h","g","1","1","-1"}) && match_expectation(type, {"1","g"}, {"2","-2","1","g","h"}) && match_expectation(type, {"g","1"}, {"g","h","1","-2","2"}) && match_expectation(type, {"1","g"}, {"Wp","3","-4","g","1"}) && match_expectation(type, {"3","g"}, {"-2","Wm","1","3","g"}) && match_expectation(type, {"g","1"}, {"1","g","Wm","-3","4"}) && match_expectation(type, {"g","-3"}, {"g","-3","-1","Wp","2"}) && match_expectation(type, {"g","2"}, {"2","g","Z_photon_mix","4","-4"}) // 4j- gluon in qqx && match_expectation(type, {"g","1"}, {"1","g","-1","1"}) && match_expectation(type, {"1","g"}, {"1","-1","g","1"}) && match_expectation(type, {"g","1"}, {"1","g","Wm","-2","1"}) && match_expectation(type, {"2","g"}, {"2","-2","g","Wp","1"}) && match_expectation(type, {"g","g"}, {"Wp","3","g","-4","g"}) && match_expectation(type, {"1","g"}, {"1","h","-1","g","1"}) && match_expectation(type, {"3","g"}, {"3","1","g","Z_photon_mix","-1"}) // 6j - two qqx && match_expectation(type, {"g","g"}, {"1","-1","g","g","1","-1"}) && match_expectation(type, {"g","g"}, {"1","-1","g","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","-1","g","1","-1"}) && match_expectation(type, {"g","g"}, {"g","1","-1","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","1","-1","-1","g"}) && match_expectation(type, {"g","g"}, {"h","1","-1","g","g","1","-1"}) && match_expectation(type, {"g","g"}, {"1","Wp","-2","g","1","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","Wp","-1","g","1","-2"}) && match_expectation(type, {"g","g"}, {"g","1","-1","Wm","2","-1","g"}) && match_expectation(type, {"g","g"}, {"g","1","2","-1","Wm","-1","g"}) && match_expectation(type, {"g","g"}, {"2","-2","g","-1","1","Z_photon_mix","g"}) // random stuff (can be non-physical) && match_expectation(type, {"g","g"}, {"1","-2","2","-1"}) // != 2 qqx && match_expectation(type, {"g","g"}, {"1","-2","2","g"}) // could be qqx && match_expectation(type, {"e+","e-"},{"1","-1"}) // bad initial state && match_expectation(type, {"1","e-"}, {"g","1","Wm"}) // bad initial state && match_expectation(type, {"h","g"}, {"g","g"}) // bad initial state && match_expectation(type, {"-1","g"}, {"-1","1","1"}) // bad qqx && match_expectation(type, {"-1","g"}, {"1","1","-1"}) // crossing in bad qqx && match_expectation(type, {"-1","g"}, {"-2","1","1","Wp"}) // bad qqx && match_expectation(type, {"1","2"}, {"1","-1","g","g","g","2"}) // bad qqx && match_expectation(type, {"1","2"}, {"1","-1","-2","g","g","2"}) // gluon in bad qqx && match_expectation(type, {"g","g"}, {"-1","2","g","g"}) // wrong back qqx && match_expectation(type, {"g","g"}, {"g","g","2","1"}) // wrong forward qqx && match_expectation(type, {"g","g"}, {"g","-2","1","g"}) // wrong central qqx && match_expectation(type, {"1","g"}, {"1","-2","g","g","Wp"}) // extra quark && match_expectation(type, {"g","1"}, {"g","g","-2","1","Wp"}) // extra quark && match_expectation(type, {"g","1"}, {"g","g","Wp","-2","1"}) // extra quark && match_expectation(type, {"g","1"}, {"g","-2","1","g","Wp"}) // extra quark && match_expectation(type, {"g","g"}, {"g","g","g","-2","1","-1","Wp"}) // extra quark && match_expectation(type, {"1","g"}, {"g","Wp","1","-2","g"}) // extra quark && match_expectation(type, {"g","g"}, {"1","-1","-2","g","g","g","Wp"}) // extra quark ; } // Two boson states, that are currently not implemented // Check for supported final states bool check_bad_FS(){ auto type{ HEJ::event_type::bad_final_state}; return match_expectation(type, {"g","g"}, {"g","p","p","g"}) // protons && match_expectation(type, {"-4","-1"},{"-4","g","11","-11","-2"}) // leptons should be in decay && match_expectation(type, {"-4","-1"},{"-4","g","-13","g","-2"}) // leptons should be in decay ; } // not 2 jets bool check_not_2_jets(){ auto type{ HEJ::event_type::no_2_jets}; return match_expectation(type, {"g","g"}, {}) && match_expectation(type, {"1","-1"}, {}) && match_expectation(type, {"g","-1"}, {"-1"}) && match_expectation(type, {"g","g"}, {"g"}) ; } // 2 boson final-states bool check_2_boson(bool const implemented=true){ auto type = (implemented) ? HEJ::event_type::FKL : HEJ::event_type::non_resummable; // Implemented if ( implemented ) { // Wp+Wp for(std::string const & q1 : all_quarks){ for(std::string const & q2 : all_quarks){ for(size_t njet = 2; njet < 6; ++njet){ // Event setup std::array in {q1, q2}; std::vector out(njet, "g"); out.front() = q1; out.back() = q2; // Adjust quark flavours for the boson, if possible if( !couple_quark("Wp", out.front()) || !couple_quark("Wp", out.back() ) ) { continue; } // skip incompatible std::unordered_map< size_t, std::vector > decays; for(size_t i = 0; i < 2; ++i){ out.emplace_back("Wp"); decays.emplace( out.size()-1, std::initializer_list{"e+", "nu_e"} ); } // Check classification if( !match_expectation(type, rapidity_order_ps(in, out, false, decays)) ) return false; } } } return true; } // Not Implemented (non_resummable) else { return // h h match_expectation(type, {"g","g"}, {"h","g","h","g"}) // Wp Wm && match_expectation(type, rapidity_order_ps( {"u","d"}, {"e+","nu_e","d","mu-","nu_mu_bar","u"}, true) ) // Wm Wm && match_expectation(type, rapidity_order_ps( {"d","d"}, {"e-","nu_e_bar","u","mu-","nu_mu_bar","u"}, true) ) // Wm h && match_expectation(type, rapidity_order_ps( {"u","d"}, {"e-","nu_e_bar","u","h","u"}, true) ) // Wp h && match_expectation(type, rapidity_order_ps( {"u","d"}, {"e+","nu_e","d","h","u"}, true) ) ; } } // not implemented processes bool check_not_implemented(){ return check_fkl(false) && check_uno(false) && check_extremal_qqx(false) && check_central_qqx(false) && check_2_boson(false); } } int main() { // tests for "no false negatives" // i.e. all HEJ-configurations get classified correctly if(!check_fkl()) return EXIT_FAILURE; if(!check_uno()) return EXIT_FAILURE; if(!check_extremal_qqx()) return EXIT_FAILURE; if(!check_central_qqx()) return EXIT_FAILURE; if(!check_2_boson()) return EXIT_FAILURE; // test for "no false positive" // i.e. non-resummable gives non-resummable if(!check_non_resummable()) return EXIT_FAILURE; if(!check_bad_FS()) return EXIT_FAILURE; if(!check_not_2_jets()) return EXIT_FAILURE; if(!check_not_implemented()) return EXIT_FAILURE; return EXIT_SUCCESS; }