File indexing completed on 2021-02-14 13:07:31
0001 #ifndef TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC
0002 #define TRUTHRIVETTOOLS_HIGGSTEMPLATECROSSSECTIONS_CC
0003
0004
0005 #include "Rivet/Analysis.hh"
0006 #include "Rivet/Particle.hh"
0007 #include "Rivet/Projections/FastJets.hh"
0008
0009
0010
0011 #include "SimDataFormats/HTXS/interface/HiggsTemplateCrossSections.h" //
0012
0013 #include <atomic>
0014
0015 namespace Rivet {
0016
0017
0018
0019
0020
0021 class HiggsTemplateCrossSections : public Analysis {
0022 public:
0023
0024 HiggsTemplateCrossSections() : Analysis("HiggsTemplateCrossSections"), m_HiggsProdMode(HTXS::UNKNOWN) {}
0025
0026 public:
0027
0028
0029
0030
0031
0032
0033 Particle getLastInstance(Particle ptcl) {
0034 if (ptcl.genParticle()->end_vertex()) {
0035 if (!hasChild(ptcl.genParticle(), ptcl.pid()))
0036 return ptcl;
0037 else
0038 return getLastInstance(ptcl.children()[0]);
0039 }
0040 return ptcl;
0041 }
0042
0043
0044 bool originateFrom(const Particle &p, const Particles &ptcls) {
0045 const ConstGenVertexPtr prodVtx = p.genParticle()->production_vertex();
0046 if (prodVtx == nullptr)
0047 return false;
0048
0049 for (const auto &ancestor : HepMCUtils::particles(prodVtx, HepMC::ancestors)) {
0050 for (const auto &part : ptcls)
0051 if (ancestor == part.genParticle())
0052 return true;
0053 }
0054
0055 return false;
0056 }
0057
0058
0059 bool originateFrom(const Particle &p, const Particle &p2) {
0060 Particles ptcls = {p2};
0061 return originateFrom(p, ptcls);
0062 }
0063
0064
0065 bool hasChild(const ConstGenParticlePtr &ptcl, int pdgID) {
0066 for (const auto &child : Particle(*ptcl).children()) {
0067 if (child.pid() == pdgID) {
0068 return true;
0069 }
0070 }
0071 return false;
0072 }
0073
0074
0075 bool hasParent(const ConstGenParticlePtr &ptcl, int pdgID) {
0076 for (auto parent : HepMCUtils::particles(ptcl->production_vertex(), HepMC::parents))
0077 if (parent->pdg_id() == pdgID)
0078 return true;
0079 return false;
0080 }
0081
0082
0083 bool quarkDecay(const Particle &p) {
0084 for (const auto &child : p.children()) {
0085 if (PID::isQuark(child.pid())) {
0086 return true;
0087 }
0088 }
0089 return false;
0090 }
0091
0092
0093 bool ChLeptonDecay(const Particle &p) {
0094 for (const auto &child : p.children()) {
0095 if (PID::isChLepton(child.pid())) {
0096 return true;
0097 }
0098 }
0099 return false;
0100 }
0101
0102
0103
0104 HiggsClassification error(HiggsClassification &cat,
0105 HTXS::ErrorCode err,
0106 std::string msg = "",
0107 int NmaxWarnings = 20) {
0108
0109 cat.errorCode = err;
0110 ++m_errorCount[err];
0111
0112
0113 static std::atomic<int> Nwarnings{0};
0114 if (!msg.empty() && ++Nwarnings < NmaxWarnings)
0115 MSG_WARNING(msg);
0116
0117 return cat;
0118 }
0119
0120
0121
0122 HiggsClassification classifyEvent(const Event &event, const HTXS::HiggsProdMode prodMode) {
0123 if (m_HiggsProdMode == HTXS::UNKNOWN)
0124 m_HiggsProdMode = prodMode;
0125
0126
0127 HiggsClassification cat;
0128 cat.prodMode = prodMode;
0129 cat.errorCode = HTXS::UNDEFINED;
0130 cat.stage0_cat = HTXS::Stage0::UNKNOWN;
0131 cat.stage1_cat_pTjet25GeV = HTXS::Stage1::UNKNOWN;
0132 cat.stage1_cat_pTjet30GeV = HTXS::Stage1::UNKNOWN;
0133 cat.stage1_1_cat_pTjet25GeV = HTXS::Stage1_1::UNKNOWN;
0134 cat.stage1_1_cat_pTjet30GeV = HTXS::Stage1_1::UNKNOWN;
0135 cat.stage1_1_fine_cat_pTjet25GeV = HTXS::Stage1_1_Fine::UNKNOWN;
0136 cat.stage1_1_fine_cat_pTjet30GeV = HTXS::Stage1_1_Fine::UNKNOWN;
0137 cat.stage1_2_cat_pTjet25GeV = HTXS::Stage1_2::UNKNOWN;
0138 cat.stage1_2_cat_pTjet30GeV = HTXS::Stage1_2::UNKNOWN;
0139 cat.stage1_2_fine_cat_pTjet25GeV = HTXS::Stage1_2_Fine::UNKNOWN;
0140 cat.stage1_2_fine_cat_pTjet30GeV = HTXS::Stage1_2_Fine::UNKNOWN;
0141
0142 if (prodMode == HTXS::UNKNOWN)
0143 return error(cat,
0144 HTXS::PRODMODE_DEFINED,
0145 "Unkown Higgs production mechanism. Cannot classify event."
0146 " Classification for all events will most likely fail.");
0147
0148
0149
0150
0151
0152
0153
0154 ConstGenVertexPtr HSvtx = event.genEvent()->signal_process_vertex();
0155 int Nhiggs = 0;
0156 for (const ConstGenParticlePtr ptcl : HepMCUtils::particles(event.genEvent())) {
0157
0158 if (!PID::isHiggs(ptcl->pdg_id()))
0159 continue;
0160
0161 if (ptcl->end_vertex() && !hasChild(ptcl, PID::HIGGS)) {
0162 cat.higgs = Particle(ptcl);
0163 ++Nhiggs;
0164 }
0165
0166
0167 if (HSvtx == nullptr && ptcl->production_vertex() && !hasParent(ptcl, PID::HIGGS))
0168 HSvtx = ptcl->production_vertex();
0169 }
0170
0171
0172 if (Nhiggs != 1)
0173 return error(cat,
0174 HTXS::HIGGS_IDENTIFICATION,
0175 "Current event has " + std::to_string(Nhiggs) + " Higgs bosons. There must be only one.");
0176 if (cat.higgs.children().size() < 2)
0177 return error(cat, HTXS::HIGGS_DECAY_IDENTIFICATION, "Could not identify Higgs boson decay products.");
0178
0179 if (HSvtx == nullptr)
0180 return error(cat, HTXS::HS_VTX_IDENTIFICATION, "Cannot find hard-scatter vertex of current event.");
0181
0182
0183
0184
0185
0186
0187
0188 bool is_uncatdV = false;
0189 Particles uncatV_decays;
0190 FourMomentum uncatV_p4(0, 0, 0, 0);
0191 FourVector uncatV_v4(0, 0, 0, 0);
0192 int nWs = 0, nZs = 0, nTop = 0;
0193 if (isVH(prodMode)) {
0194 for (auto ptcl : HepMCUtils::particles(HSvtx, HepMC::children)) {
0195 if (PID::isW(ptcl->pdg_id())) {
0196 ++nWs;
0197 cat.V = Particle(ptcl);
0198 }
0199 if (PID::isZ(ptcl->pdg_id())) {
0200 ++nZs;
0201 cat.V = Particle(ptcl);
0202 }
0203 }
0204 if (nWs + nZs > 0)
0205 cat.V = getLastInstance(cat.V);
0206 else {
0207 for (auto ptcl : HepMCUtils::particles(HSvtx, HepMC::children)) {
0208 if (!PID::isHiggs(ptcl->pdg_id())) {
0209 uncatV_decays += Particle(ptcl);
0210 uncatV_p4 += Particle(ptcl).momentum();
0211
0212 }
0213 }
0214
0215 is_uncatdV = true;
0216 cat.V = Particle(24, uncatV_p4);
0217 }
0218 }
0219
0220 if (!is_uncatdV) {
0221 if (isVH(prodMode) && !cat.V.genParticle()->end_vertex())
0222 return error(cat, HTXS::VH_DECAY_IDENTIFICATION, "Vector boson does not decay!");
0223
0224 if (isVH(prodMode) && cat.V.children().size() < 2)
0225 return error(cat, HTXS::VH_DECAY_IDENTIFICATION, "Vector boson does not decay!");
0226
0227 if ((prodMode == HTXS::WH && (nZs > 0 || nWs != 1)) ||
0228 ((prodMode == HTXS::QQ2ZH || prodMode == HTXS::GG2ZH) && (nZs != 1 || nWs > 0)))
0229 return error(cat,
0230 HTXS::VH_IDENTIFICATION,
0231 "Found " + std::to_string(nWs) + " W-bosons and " + std::to_string(nZs) +
0232 " Z-bosons. Inconsitent with VH expectation.");
0233 }
0234
0235
0236 Particles Ws;
0237 if (prodMode == HTXS::TTH || prodMode == HTXS::TH) {
0238
0239 for (auto ptcl : HepMCUtils::particles(HSvtx, HepMC::children)) {
0240 if (!PID::isTop(ptcl->pdg_id()))
0241 continue;
0242 ++nTop;
0243 Particle top = getLastInstance(Particle(ptcl));
0244 if (top.genParticle()->end_vertex())
0245 for (const auto &child : top.children())
0246 if (PID::isW(child.pid()))
0247 Ws += getLastInstance(child);
0248 }
0249 }
0250
0251
0252 if ((prodMode == HTXS::TTH && Ws.size() < 2) || (prodMode == HTXS::TH && Ws.empty()))
0253 return error(cat, HTXS::TOP_W_IDENTIFICATION, "Failed to identify W-boson(s) from t-decay!");
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264 Particles leptonicVs;
0265 if (!is_uncatdV) {
0266 if (isVH(prodMode) && !quarkDecay(cat.V))
0267 leptonicVs += cat.V;
0268 } else
0269 leptonicVs = uncatV_decays;
0270 for (const auto &W : Ws)
0271 if (W.genParticle()->end_vertex() && !quarkDecay(W))
0272 leptonicVs += W;
0273
0274
0275 const Particles FS = apply<FinalState>(event, "FS").particles();
0276 Particles hadrons;
0277 FourMomentum sum(0, 0, 0, 0), vSum(0, 0, 0, 0), hSum(0, 0, 0, 0);
0278 for (const Particle &p : FS) {
0279
0280 sum += p.momentum();
0281
0282 if (originateFrom(p, cat.higgs)) {
0283 hSum += p.momentum();
0284 continue;
0285 }
0286
0287 if (isVH(prodMode) && !is_uncatdV && originateFrom(p, Ws))
0288 vSum += p.momentum();
0289
0290 if (!leptonicVs.empty() && originateFrom(p, leptonicVs))
0291 continue;
0292
0293 hadrons += p;
0294 }
0295
0296 cat.p4decay_higgs = hSum;
0297 cat.p4decay_V = vSum;
0298
0299 FinalState fps_temp;
0300 FastJets jets(fps_temp, FastJets::ANTIKT, 0.4);
0301 jets.calc(hadrons);
0302
0303 cat.jets25 = jets.jetsByPt(Cuts::pT > 25.0);
0304 cat.jets30 = jets.jetsByPt(Cuts::pT > 30.0);
0305
0306
0307
0308
0309
0310
0311
0312
0313
0314 if (is_uncatdV)
0315 return error(cat, HTXS::VH_IDENTIFICATION, "Failed to identify associated V-boson!");
0316
0317
0318
0319
0320
0321
0322
0323 cat.isZ2vvDecay = false;
0324 if ((prodMode == HTXS::GG2ZH || prodMode == HTXS::QQ2ZH) && !quarkDecay(cat.V) && !ChLeptonDecay(cat.V))
0325 cat.isZ2vvDecay = true;
0326 cat.stage0_cat = getStage0Category(prodMode, cat.higgs, cat.V);
0327 cat.stage1_cat_pTjet25GeV = getStage1Category(prodMode, cat.higgs, cat.jets25, cat.V);
0328 cat.stage1_cat_pTjet30GeV = getStage1Category(prodMode, cat.higgs, cat.jets30, cat.V);
0329 cat.stage1_1_cat_pTjet25GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0330 cat.stage1_1_cat_pTjet30GeV = getStage1_1_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0331 cat.stage1_1_fine_cat_pTjet25GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0332 cat.stage1_1_fine_cat_pTjet30GeV = getStage1_1_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0333 cat.stage1_2_cat_pTjet25GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0334 cat.stage1_2_cat_pTjet30GeV = getStage1_2_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0335 cat.stage1_2_fine_cat_pTjet25GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets25, cat.V);
0336 cat.stage1_2_fine_cat_pTjet30GeV = getStage1_2_Fine_Category(prodMode, cat.higgs, cat.jets30, cat.V);
0337
0338 cat.errorCode = HTXS::SUCCESS;
0339 ++m_errorCount[HTXS::SUCCESS];
0340 ++m_sumevents;
0341
0342 return cat;
0343 }
0344
0345
0346
0347
0348
0349
0350
0351
0352 int getBin(double x, std::vector<double> bins) {
0353 for (size_t i = 1; i < bins.size(); ++i)
0354 if (x < bins[i])
0355 return i - 1;
0356 return bins.size() - 1;
0357 }
0358
0359
0360
0361
0362 int vbfTopology(const Jets &jets, const Particle &higgs) {
0363 if (jets.size() < 2)
0364 return 0;
0365 const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0366 bool VBFtopo = (j1 + j2).mass() > 400.0 && std::abs(j1.rapidity() - j2.rapidity()) > 2.8;
0367 return VBFtopo ? (j1 + j2 + higgs.momentum()).pt() < 25 ? 2 : 1 : 0;
0368 }
0369
0370
0371
0372
0373
0374 int vbfTopology_Stage1_X(const Jets &jets, const Particle &higgs) {
0375 if (jets.size() < 2)
0376 return 0;
0377 const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0378 double mjj = (j1 + j2).mass();
0379 if (mjj > 350 && mjj <= 700)
0380 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0381 else if (mjj > 700)
0382 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0383 else
0384 return 0;
0385 }
0386
0387
0388
0389
0390
0391
0392
0393 int vbfTopology_Stage1_X_Fine(const Jets &jets, const Particle &higgs) {
0394 if (jets.size() < 2)
0395 return 0;
0396 const FourMomentum &j1 = jets[0].momentum(), &j2 = jets[1].momentum();
0397 double mjj = (j1 + j2).mass();
0398 if (mjj > 350 && mjj <= 700)
0399 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 1 : 2;
0400 else if (mjj > 700 && mjj <= 1000)
0401 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 3 : 4;
0402 else if (mjj > 1000 && mjj <= 1500)
0403 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 5 : 6;
0404 else if (mjj > 1500)
0405 return (j1 + j2 + higgs.momentum()).pt() < 25 ? 7 : 8;
0406 else
0407 return 0;
0408 }
0409
0410
0411 bool isVH(HTXS::HiggsProdMode p) { return p == HTXS::WH || p == HTXS::QQ2ZH || p == HTXS::GG2ZH; }
0412
0413
0414 HTXS::Stage0::Category getStage0Category(const HTXS::HiggsProdMode prodMode,
0415 const Particle &higgs,
0416 const Particle &V) {
0417 using namespace HTXS::Stage0;
0418 int ctrlHiggs = std::abs(higgs.rapidity()) < 2.5;
0419
0420 if ((prodMode == HTXS::WH || prodMode == HTXS::QQ2ZH) && quarkDecay(V)) {
0421 return ctrlHiggs ? VH2HQQ : VH2HQQ_FWDH;
0422 } else if (prodMode == HTXS::GG2ZH && quarkDecay(V)) {
0423 return Category(HTXS::GGF * 10 + ctrlHiggs);
0424 }
0425
0426 return Category(prodMode * 10 + ctrlHiggs);
0427 }
0428
0429
0430 HTXS::Stage1::Category getStage1Category(const HTXS::HiggsProdMode prodMode,
0431 const Particle &higgs,
0432 const Jets &jets,
0433 const Particle &V) {
0434 using namespace HTXS::Stage1;
0435 int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0436 double pTj1 = !jets.empty() ? jets[0].momentum().pt() : 0;
0437 int vbfTopo = vbfTopology(jets, higgs);
0438
0439
0440
0441 if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0442 if (fwdHiggs)
0443 return GG2H_FWDH;
0444 if (Njets == 0)
0445 return GG2H_0J;
0446 else if (Njets == 1)
0447 return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0448 else if (Njets >= 2) {
0449
0450 if (higgs.pt() <= 200) {
0451 if (vbfTopo == 2)
0452 return GG2H_VBFTOPO_JET3VETO;
0453 else if (vbfTopo == 1)
0454 return GG2H_VBFTOPO_JET3;
0455 }
0456
0457 return Category(GG2H_GE2J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0458 }
0459 }
0460
0461 else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0462 if (std::abs(higgs.rapidity()) > 2.5)
0463 return QQ2HQQ_FWDH;
0464 if (pTj1 > 200)
0465 return QQ2HQQ_PTJET1_GT200;
0466 if (vbfTopo == 2)
0467 return QQ2HQQ_VBFTOPO_JET3VETO;
0468 if (vbfTopo == 1)
0469 return QQ2HQQ_VBFTOPO_JET3;
0470 double mjj = jets.size() > 1 ? (jets[0].mom() + jets[1].mom()).mass() : 0;
0471 if (60 < mjj && mjj < 120)
0472 return QQ2HQQ_VH2JET;
0473 return QQ2HQQ_REST;
0474 }
0475
0476 else if (prodMode == HTXS::WH) {
0477 if (fwdHiggs)
0478 return QQ2HLNU_FWDH;
0479 else if (V.pt() < 150)
0480 return QQ2HLNU_PTV_0_150;
0481 else if (V.pt() > 250)
0482 return QQ2HLNU_PTV_GT250;
0483
0484 return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0485 }
0486
0487 else if (prodMode == HTXS::QQ2ZH) {
0488 if (fwdHiggs)
0489 return QQ2HLL_FWDH;
0490 else if (V.pt() < 150)
0491 return QQ2HLL_PTV_0_150;
0492 else if (V.pt() > 250)
0493 return QQ2HLL_PTV_GT250;
0494
0495 return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0496 }
0497
0498 else if (prodMode == HTXS::GG2ZH) {
0499 if (fwdHiggs)
0500 return GG2HLL_FWDH;
0501 if (V.pt() < 150)
0502 return GG2HLL_PTV_0_150;
0503 else if (jets.empty())
0504 return GG2HLL_PTV_GT150_0J;
0505 return GG2HLL_PTV_GT150_GE1J;
0506 }
0507
0508 else if (prodMode == HTXS::TTH)
0509 return Category(TTH_FWDH + ctrlHiggs);
0510 else if (prodMode == HTXS::BBH)
0511 return Category(BBH_FWDH + ctrlHiggs);
0512 else if (prodMode == HTXS::TH)
0513 return Category(TH_FWDH + ctrlHiggs);
0514 return UNKNOWN;
0515 }
0516
0517
0518 HTXS::Stage1_1::Category getStage1_1_Category(const HTXS::HiggsProdMode prodMode,
0519 const Particle &higgs,
0520 const Jets &jets,
0521 const Particle &V) {
0522 using namespace HTXS::Stage1_1;
0523 int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0524 int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0525
0526
0527
0528 if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0529 if (fwdHiggs)
0530 return GG2H_FWDH;
0531 if (higgs.pt() > 200)
0532 return GG2H_PTH_GT200;
0533 if (Njets == 0)
0534 return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0535 if (Njets == 1)
0536 return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0537 if (Njets > 1) {
0538
0539 if (vbfTopo)
0540 return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0541
0542 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0543 }
0544 }
0545
0546
0547 else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0548 if (std::abs(higgs.rapidity()) > 2.5)
0549 return QQ2HQQ_FWDH;
0550 int Njets = jets.size();
0551 if (Njets == 0)
0552 return QQ2HQQ_0J;
0553 else if (Njets == 1)
0554 return QQ2HQQ_1J;
0555 else if (Njets >= 2) {
0556 double mjj = (jets[0].mom() + jets[1].mom()).mass();
0557 if (mjj < 60)
0558 return QQ2HQQ_MJJ_0_60;
0559 else if (60 < mjj && mjj < 120)
0560 return QQ2HQQ_MJJ_60_120;
0561 else if (120 < mjj && mjj < 350)
0562 return QQ2HQQ_MJJ_120_350;
0563 else if (mjj > 350) {
0564 if (higgs.pt() > 200)
0565 return QQ2HQQ_MJJ_GT350_PTH_GT200;
0566 if (vbfTopo)
0567 return Category(QQ2HQQ_MJJ_GT350_PTH_GT200 + vbfTopo);
0568 }
0569 }
0570 }
0571
0572 else if (prodMode == HTXS::WH) {
0573 if (fwdHiggs)
0574 return QQ2HLNU_FWDH;
0575 else if (V.pt() < 75)
0576 return QQ2HLNU_PTV_0_75;
0577 else if (V.pt() < 150)
0578 return QQ2HLNU_PTV_75_150;
0579 else if (V.pt() > 250)
0580 return QQ2HLNU_PTV_GT250;
0581
0582 return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0583 }
0584
0585 else if (prodMode == HTXS::QQ2ZH) {
0586 if (fwdHiggs)
0587 return QQ2HLL_FWDH;
0588 else if (V.pt() < 75)
0589 return QQ2HLL_PTV_0_75;
0590 else if (V.pt() < 150)
0591 return QQ2HLL_PTV_75_150;
0592 else if (V.pt() > 250)
0593 return QQ2HLL_PTV_GT250;
0594
0595 return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0596 }
0597
0598 else if (prodMode == HTXS::GG2ZH) {
0599 if (fwdHiggs)
0600 return GG2HLL_FWDH;
0601 else if (V.pt() < 75)
0602 return GG2HLL_PTV_0_75;
0603 else if (V.pt() < 150)
0604 return GG2HLL_PTV_75_150;
0605 else if (V.pt() > 250)
0606 return GG2HLL_PTV_GT250;
0607 return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0608 }
0609
0610 else if (prodMode == HTXS::TTH)
0611 return Category(TTH_FWDH + ctrlHiggs);
0612 else if (prodMode == HTXS::BBH)
0613 return Category(BBH_FWDH + ctrlHiggs);
0614 else if (prodMode == HTXS::TH)
0615 return Category(TH_FWDH + ctrlHiggs);
0616 return UNKNOWN;
0617 }
0618
0619
0620 HTXS::Stage1_1_Fine::Category getStage1_1_Fine_Category(const HTXS::HiggsProdMode prodMode,
0621 const Particle &higgs,
0622 const Jets &jets,
0623 const Particle &V) {
0624 using namespace HTXS::Stage1_1_Fine;
0625 int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0626 int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0627
0628
0629
0630 if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0631 if (fwdHiggs)
0632 return GG2H_FWDH;
0633 if (higgs.pt() > 200)
0634 return GG2H_PTH_GT200;
0635 if (Njets == 0)
0636 return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0637 if (Njets == 1)
0638 return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0639 if (Njets > 1) {
0640
0641 double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0642
0643 if (vbfTopo)
0644 return Category(GG2H_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0645
0646 if (pTHjj < 25)
0647 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0648 else
0649 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0650 }
0651 }
0652
0653
0654 else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0655 if (std::abs(higgs.rapidity()) > 2.5)
0656 return QQ2HQQ_FWDH;
0657 int Njets = jets.size();
0658 if (Njets == 0)
0659 return QQ2HQQ_0J;
0660 else if (Njets == 1)
0661 return QQ2HQQ_1J;
0662 else if (Njets >= 2) {
0663 double mjj = (jets[0].mom() + jets[1].mom()).mass();
0664 double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0665 if (mjj < 350) {
0666 if (pTHjj < 25)
0667 return Category(QQ2HQQ_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0668 else
0669 return Category(QQ2HQQ_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0670 } else {
0671 if (higgs.pt() < 200) {
0672 return Category(QQ2HQQ_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0673 } else {
0674 return Category(QQ2HQQ_PTH_GT200_MJJ_350_700_PTHJJ_0_25 + vbfTopo - 1);
0675 }
0676 }
0677 }
0678 }
0679
0680 else if (prodMode == HTXS::WH) {
0681 if (fwdHiggs)
0682 return QQ2HLNU_FWDH;
0683 int Njets = jets.size();
0684 if (Njets == 0)
0685 return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0686 if (Njets == 1)
0687 return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0688 return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0689 }
0690
0691 else if (prodMode == HTXS::QQ2ZH) {
0692 if (fwdHiggs)
0693 return QQ2HLL_FWDH;
0694 int Njets = jets.size();
0695 if (Njets == 0)
0696 return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0697 if (Njets == 1)
0698 return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0699 return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0700 }
0701
0702 else if (prodMode == HTXS::GG2ZH) {
0703 if (fwdHiggs)
0704 return GG2HLL_FWDH;
0705 int Njets = jets.size();
0706 if (Njets == 0)
0707 return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0708 if (Njets == 1)
0709 return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0710 return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0711 }
0712
0713 else if (prodMode == HTXS::TTH)
0714 return Category(TTH_FWDH + ctrlHiggs);
0715 else if (prodMode == HTXS::BBH)
0716 return Category(BBH_FWDH + ctrlHiggs);
0717 else if (prodMode == HTXS::TH)
0718 return Category(TH_FWDH + ctrlHiggs);
0719 return UNKNOWN;
0720 }
0721
0722
0723 HTXS::Stage1_2::Category getStage1_2_Category(const HTXS::HiggsProdMode prodMode,
0724 const Particle &higgs,
0725 const Jets &jets,
0726 const Particle &V) {
0727 using namespace HTXS::Stage1_2;
0728 int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0729 int vbfTopo = vbfTopology_Stage1_X(jets, higgs);
0730
0731
0732
0733 if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0734 if (fwdHiggs)
0735 return GG2H_FWDH;
0736 if (higgs.pt() > 200)
0737 return Category(GG2H_PTH_200_300 + getBin(higgs.pt(), {200, 300, 450, 650}));
0738 if (Njets == 0)
0739 return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0740 if (Njets == 1)
0741 return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0742 if (Njets > 1) {
0743
0744 if (vbfTopo)
0745 return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0746
0747 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0748 }
0749 }
0750
0751
0752 else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0753 if (std::abs(higgs.rapidity()) > 2.5)
0754 return QQ2HQQ_FWDH;
0755 int Njets = jets.size();
0756 if (Njets == 0)
0757 return QQ2HQQ_0J;
0758 else if (Njets == 1)
0759 return QQ2HQQ_1J;
0760 else if (Njets >= 2) {
0761 double mjj = (jets[0].mom() + jets[1].mom()).mass();
0762 if (mjj < 60)
0763 return QQ2HQQ_GE2J_MJJ_0_60;
0764 else if (60 < mjj && mjj < 120)
0765 return QQ2HQQ_GE2J_MJJ_60_120;
0766 else if (120 < mjj && mjj < 350)
0767 return QQ2HQQ_GE2J_MJJ_120_350;
0768 else if (mjj > 350) {
0769 if (higgs.pt() > 200)
0770 return QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200;
0771 if (vbfTopo)
0772 return Category(QQ2HQQ_GE2J_MJJ_GT350_PTH_GT200 + vbfTopo);
0773 }
0774 }
0775 }
0776
0777 else if (prodMode == HTXS::WH) {
0778 if (fwdHiggs)
0779 return QQ2HLNU_FWDH;
0780 else if (V.pt() < 75)
0781 return QQ2HLNU_PTV_0_75;
0782 else if (V.pt() < 150)
0783 return QQ2HLNU_PTV_75_150;
0784 else if (V.pt() > 250)
0785 return QQ2HLNU_PTV_GT250;
0786
0787 return jets.empty() ? QQ2HLNU_PTV_150_250_0J : QQ2HLNU_PTV_150_250_GE1J;
0788 }
0789
0790 else if (prodMode == HTXS::QQ2ZH) {
0791 if (fwdHiggs)
0792 return QQ2HLL_FWDH;
0793 else if (V.pt() < 75)
0794 return QQ2HLL_PTV_0_75;
0795 else if (V.pt() < 150)
0796 return QQ2HLL_PTV_75_150;
0797 else if (V.pt() > 250)
0798 return QQ2HLL_PTV_GT250;
0799
0800 return jets.empty() ? QQ2HLL_PTV_150_250_0J : QQ2HLL_PTV_150_250_GE1J;
0801 }
0802
0803 else if (prodMode == HTXS::GG2ZH) {
0804 if (fwdHiggs)
0805 return GG2HLL_FWDH;
0806 else if (V.pt() < 75)
0807 return GG2HLL_PTV_0_75;
0808 else if (V.pt() < 150)
0809 return GG2HLL_PTV_75_150;
0810 else if (V.pt() > 250)
0811 return GG2HLL_PTV_GT250;
0812 return jets.empty() ? GG2HLL_PTV_150_250_0J : GG2HLL_PTV_150_250_GE1J;
0813 }
0814
0815 else if (prodMode == HTXS::TTH) {
0816 if (fwdHiggs)
0817 return TTH_FWDH;
0818 else
0819 return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300}));
0820 } else if (prodMode == HTXS::BBH)
0821 return Category(BBH_FWDH + ctrlHiggs);
0822 else if (prodMode == HTXS::TH)
0823 return Category(TH_FWDH + ctrlHiggs);
0824 return UNKNOWN;
0825 }
0826
0827
0828 HTXS::Stage1_2_Fine::Category getStage1_2_Fine_Category(const HTXS::HiggsProdMode prodMode,
0829 const Particle &higgs,
0830 const Jets &jets,
0831 const Particle &V) {
0832 using namespace HTXS::Stage1_2_Fine;
0833 int Njets = jets.size(), ctrlHiggs = std::abs(higgs.rapidity()) < 2.5, fwdHiggs = !ctrlHiggs;
0834 int vbfTopo = vbfTopology_Stage1_X_Fine(jets, higgs);
0835
0836
0837
0838 if (prodMode == HTXS::GGF || (prodMode == HTXS::GG2ZH && quarkDecay(V))) {
0839 if (fwdHiggs)
0840 return GG2H_FWDH;
0841 if (higgs.pt() > 200) {
0842 if (Njets > 0) {
0843 double pTHj = (jets[0].momentum() + higgs.momentum()).pt();
0844 if (pTHj / higgs.pt() > 0.15)
0845 return Category(GG2H_PTH_200_300_PTHJoverPTH_GT15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0846 else
0847 return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0848 } else
0849 return Category(GG2H_PTH_200_300_PTHJoverPTH_0_15 + getBin(higgs.pt(), {200, 300, 450, 650}));
0850 }
0851 if (Njets == 0)
0852 return higgs.pt() < 10 ? GG2H_0J_PTH_0_10 : GG2H_0J_PTH_GT10;
0853 if (Njets == 1)
0854 return Category(GG2H_1J_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200}));
0855 if (Njets > 1) {
0856
0857 double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0858
0859 if (vbfTopo)
0860 return Category(GG2H_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0861
0862 if (pTHjj < 25)
0863 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_0_25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0864 else
0865 return Category(GG2H_GE2J_MJJ_0_350_PTH_0_60_PTHJJ_GT25 + getBin(higgs.pt(), {0, 60, 120, 200}));
0866 }
0867 }
0868
0869
0870 else if (prodMode == HTXS::VBF || (isVH(prodMode) && quarkDecay(V))) {
0871 if (std::abs(higgs.rapidity()) > 2.5)
0872 return QQ2HQQ_FWDH;
0873 int Njets = jets.size();
0874 if (Njets == 0)
0875 return QQ2HQQ_0J;
0876 else if (Njets == 1)
0877 return QQ2HQQ_1J;
0878 else if (Njets >= 2) {
0879 double mjj = (jets[0].mom() + jets[1].mom()).mass();
0880 double pTHjj = (jets[0].momentum() + jets[1].momentum() + higgs.momentum()).pt();
0881 if (mjj < 350) {
0882 if (pTHjj < 25)
0883 return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_0_25 + getBin(mjj, {0, 60, 120, 350}));
0884 else
0885 return Category(QQ2HQQ_GE2J_MJJ_0_60_PTHJJ_GT25 + getBin(mjj, {0, 60, 120, 350}));
0886 } else {
0887 if (higgs.pt() < 200) {
0888 return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_0_200_PTHJJ_0_25 + vbfTopo - 1);
0889 } else {
0890 return Category(QQ2HQQ_GE2J_MJJ_350_700_PTH_GT200_PTHJJ_0_25 + vbfTopo - 1);
0891 }
0892 }
0893 }
0894 }
0895
0896 else if (prodMode == HTXS::WH) {
0897 if (fwdHiggs)
0898 return QQ2HLNU_FWDH;
0899 int Njets = jets.size();
0900 if (Njets == 0)
0901 return Category(QQ2HLNU_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0902 if (Njets == 1)
0903 return Category(QQ2HLNU_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0904 return Category(QQ2HLNU_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0905 }
0906
0907 else if (prodMode == HTXS::QQ2ZH) {
0908 if (fwdHiggs)
0909 return QQ2HLL_FWDH;
0910 int Njets = jets.size();
0911 if (Njets == 0)
0912 return Category(QQ2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0913 if (Njets == 1)
0914 return Category(QQ2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0915 return Category(QQ2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0916 }
0917
0918 else if (prodMode == HTXS::GG2ZH) {
0919 if (fwdHiggs)
0920 return GG2HLL_FWDH;
0921 int Njets = jets.size();
0922 if (Njets == 0)
0923 return Category(GG2HLL_PTV_0_75_0J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0924 if (Njets == 1)
0925 return Category(GG2HLL_PTV_0_75_1J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0926 return Category(GG2HLL_PTV_0_75_GE2J + getBin(V.pt(), {0, 75, 150, 250, 400}));
0927 }
0928
0929 else if (prodMode == HTXS::TTH) {
0930 if (fwdHiggs)
0931 return TTH_FWDH;
0932 else
0933 return Category(TTH_PTH_0_60 + getBin(higgs.pt(), {0, 60, 120, 200, 300, 450}));
0934 } else if (prodMode == HTXS::BBH)
0935 return Category(BBH_FWDH + ctrlHiggs);
0936 else if (prodMode == HTXS::TH)
0937 return Category(TH_FWDH + ctrlHiggs);
0938 return UNKNOWN;
0939 }
0940
0941
0942
0943
0944
0945
0946
0947 void setHiggsProdMode(HTXS::HiggsProdMode prodMode) { m_HiggsProdMode = prodMode; }
0948
0949
0950
0951
0952 void init() override {
0953 printf("==============================================================\n");
0954 printf("======== HiggsTemplateCrossSections Initialization =========\n");
0955 printf("==============================================================\n");
0956
0957
0958 if (m_HiggsProdMode == HTXS::UNKNOWN) {
0959 char *pm_env = std::getenv("HIGGSPRODMODE");
0960 string pm(pm_env == nullptr ? "" : pm_env);
0961 if (pm == "GGF")
0962 m_HiggsProdMode = HTXS::GGF;
0963 else if (pm == "VBF")
0964 m_HiggsProdMode = HTXS::VBF;
0965 else if (pm == "WH")
0966 m_HiggsProdMode = HTXS::WH;
0967 else if (pm == "ZH")
0968 m_HiggsProdMode = HTXS::QQ2ZH;
0969 else if (pm == "QQ2ZH")
0970 m_HiggsProdMode = HTXS::QQ2ZH;
0971 else if (pm == "GG2ZH")
0972 m_HiggsProdMode = HTXS::GG2ZH;
0973 else if (pm == "TTH")
0974 m_HiggsProdMode = HTXS::TTH;
0975 else if (pm == "BBH")
0976 m_HiggsProdMode = HTXS::BBH;
0977 else if (pm == "TH")
0978 m_HiggsProdMode = HTXS::TH;
0979 else {
0980 MSG_WARNING("No HIGGSPRODMODE shell variable found. Needed when running Rivet stand-alone.");
0981 }
0982 }
0983
0984
0985 const FinalState FS;
0986 declare(FS, "FS");
0987
0988
0989 initializeHistos();
0990 m_sumw = 0.0;
0991 m_sumevents = 0;
0992 printf("==============================================================\n");
0993 printf("======== Higgs prod mode %d =========\n", m_HiggsProdMode);
0994 printf("======== Sucessful Initialization =========\n");
0995 printf("==============================================================\n");
0996 }
0997
0998
0999 void analyze(const Event &event) override {
1000
1001 HiggsClassification cat = classifyEvent(event, m_HiggsProdMode);
1002
1003
1004 const double weight = 1.0;
1005 m_sumw += weight;
1006
1007 int F = cat.stage0_cat % 10, P = cat.stage1_cat_pTjet30GeV / 100;
1008 hist_stage0->fill(cat.stage0_cat / 10 * 2 + F, weight);
1009
1010
1011 vector<int> offset({0, 1, 13, 19, 24, 29, 33, 35, 37, 39});
1012 int off = offset[P];
1013 hist_stage1_pTjet25->fill(cat.stage1_cat_pTjet25GeV % 100 + off, weight);
1014 hist_stage1_pTjet30->fill(cat.stage1_cat_pTjet30GeV % 100 + off, weight);
1015
1016
1017 static vector<int> offset1_2({0, 1, 18, 29, 35, 41, 47, 53, 55, 57});
1018 int off1_2 = offset1_2[P];
1019
1020 static vector<int> offset1_2f({0, 1, 29, 54, 70, 86, 102, 109, 111, 113});
1021 int off1_2f = offset1_2f[P];
1022 hist_stage1_2_pTjet25->fill(cat.stage1_2_cat_pTjet25GeV % 100 + off1_2, weight);
1023 hist_stage1_2_pTjet30->fill(cat.stage1_2_cat_pTjet30GeV % 100 + off1_2, weight);
1024 hist_stage1_2_fine_pTjet25->fill(cat.stage1_2_fine_cat_pTjet25GeV % 100 + off1_2f, weight);
1025 hist_stage1_2_fine_pTjet30->fill(cat.stage1_2_fine_cat_pTjet30GeV % 100 + off1_2f, weight);
1026
1027
1028 hist_pT_Higgs->fill(cat.higgs.pT(), weight);
1029 hist_y_Higgs->fill(cat.higgs.rapidity(), weight);
1030 hist_pT_V->fill(cat.V.pT(), weight);
1031
1032 hist_Njets25->fill(cat.jets25.size(), weight);
1033 hist_Njets30->fill(cat.jets30.size(), weight);
1034
1035 hist_isZ2vv->fill(cat.isZ2vvDecay, weight);
1036
1037
1038 if (!cat.jets30.empty())
1039 hist_pT_jet1->fill(cat.jets30[0].pt(), weight);
1040 if (cat.jets30.size() >= 2) {
1041 const FourMomentum &j1 = cat.jets30[0].momentum(), &j2 = cat.jets30[1].momentum();
1042 hist_deltay_jj->fill(std::abs(j1.rapidity() - j2.rapidity()), weight);
1043 hist_dijet_mass->fill((j1 + j2).mass(), weight);
1044 hist_pT_Hjj->fill((j1 + j2 + cat.higgs.momentum()).pt(), weight);
1045 }
1046 }
1047
1048 void printClassificationSummary() {
1049 MSG_INFO(" ====================================================== ");
1050 MSG_INFO(" Higgs Template X-Sec Categorization Tool ");
1051 MSG_INFO(" Status Code Summary ");
1052 MSG_INFO(" ====================================================== ");
1053 bool allSuccess = (m_sumevents == m_errorCount[HTXS::SUCCESS]);
1054 if (allSuccess)
1055 MSG_INFO(" >>>> All " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized!");
1056 else {
1057 MSG_INFO(" >>>> " << m_errorCount[HTXS::SUCCESS] << " events successfully categorized");
1058 MSG_INFO(" >>>> --> the following errors occured:");
1059 MSG_INFO(" >>>> " << m_errorCount[HTXS::PRODMODE_DEFINED] << " had an undefined Higgs production mode.");
1060 MSG_INFO(" >>>> " << m_errorCount[HTXS::MOMENTUM_CONSERVATION] << " failed momentum conservation.");
1061 MSG_INFO(" >>>> " << m_errorCount[HTXS::HIGGS_IDENTIFICATION]
1062 << " failed to identify a valid Higgs boson.");
1063 MSG_INFO(" >>>> " << m_errorCount[HTXS::HS_VTX_IDENTIFICATION]
1064 << " failed to identify the hard scatter vertex.");
1065 MSG_INFO(" >>>> " << m_errorCount[HTXS::VH_IDENTIFICATION] << " VH: to identify a valid V-boson.");
1066 MSG_INFO(" >>>> " << m_errorCount[HTXS::TOP_W_IDENTIFICATION]
1067 << " failed to identify valid Ws from top decay.");
1068 }
1069 MSG_INFO(" ====================================================== ");
1070 MSG_INFO(" ====================================================== ");
1071 }
1072
1073 void finalize() override {
1074 printClassificationSummary();
1075 double sf = m_sumw > 0 ? 1.0 / m_sumw : 1.0;
1076 for (const auto &hist : {hist_stage0,
1077 hist_stage1_pTjet25,
1078 hist_stage1_pTjet30,
1079 hist_stage1_2_pTjet25,
1080 hist_stage1_2_pTjet30,
1081 hist_stage1_2_fine_pTjet25,
1082 hist_stage1_2_fine_pTjet30,
1083 hist_Njets25,
1084 hist_Njets30,
1085 hist_pT_Higgs,
1086 hist_y_Higgs,
1087 hist_pT_V,
1088 hist_pT_jet1,
1089 hist_deltay_jj,
1090 hist_dijet_mass,
1091 hist_pT_Hjj})
1092 scale(hist, sf);
1093 }
1094
1095
1096
1097
1098
1099 void initializeHistos() {
1100 book(hist_stage0, "HTXS_stage0", 20, 0, 20);
1101 book(hist_stage1_pTjet25, "HTXS_stage1_pTjet25", 40, 0, 40);
1102 book(hist_stage1_pTjet30, "HTXS_stage1_pTjet30", 40, 0, 40);
1103 book(hist_stage1_2_pTjet25, "HTXS_stage1_2_pTjet25", 57, 0, 57);
1104 book(hist_stage1_2_pTjet30, "HTXS_stage1_2_pTjet30", 57, 0, 57);
1105 book(hist_stage1_2_fine_pTjet25, "HTXS_stage1_2_fine_pTjet25", 113, 0, 113);
1106 book(hist_stage1_2_fine_pTjet30, "HTXS_stage1_2_fine_pTjet30", 113, 0, 113);
1107 book(hist_pT_Higgs, "pT_Higgs", 80, 0, 400);
1108 book(hist_y_Higgs, "y_Higgs", 80, -4, 4);
1109 book(hist_pT_V, "pT_V", 80, 0, 400);
1110 book(hist_pT_jet1, "pT_jet1", 80, 0, 400);
1111 book(hist_deltay_jj, "deltay_jj", 50, 0, 10);
1112 book(hist_dijet_mass, "m_jj", 50, 0, 2000);
1113 book(hist_pT_Hjj, "pT_Hjj", 50, 0, 250);
1114 book(hist_Njets25, "Njets25", 10, 0, 10);
1115 book(hist_Njets30, "Njets30", 10, 0, 10);
1116 book(hist_isZ2vv, "isZ2vv", 2, 0, 2);
1117 }
1118
1119
1120
1121
1122
1123
1124 private:
1125 double m_sumw;
1126 size_t m_sumevents;
1127 HTXS::HiggsProdMode m_HiggsProdMode;
1128 std::map<HTXS::ErrorCode, size_t> m_errorCount;
1129 Histo1DPtr hist_stage0;
1130 Histo1DPtr hist_stage1_pTjet25, hist_stage1_pTjet30;
1131 Histo1DPtr hist_stage1_2_pTjet25, hist_stage1_2_pTjet30;
1132 Histo1DPtr hist_stage1_2_fine_pTjet25, hist_stage1_2_fine_pTjet30;
1133 Histo1DPtr hist_pT_Higgs, hist_y_Higgs;
1134 Histo1DPtr hist_pT_V, hist_pT_jet1;
1135 Histo1DPtr hist_deltay_jj, hist_dijet_mass, hist_pT_Hjj;
1136 Histo1DPtr hist_Njets25, hist_Njets30;
1137 Histo1DPtr hist_isZ2vv;
1138 };
1139
1140
1141
1142
1143 #ifdef RIVET_ANALYSIS_PATH
1144
1145 DECLARE_RIVET_PLUGIN(HiggsTemplateCrossSections);
1146 #endif
1147
1148 }
1149
1150 #endif