File indexing completed on 2024-04-06 12:13:56
0001
0002
0003
0004
0005 #include "JetMatchingEWKFxFx.h"
0006 using namespace Pythia8;
0007
0008 JetMatchingEWKFxFx::JetMatchingEWKFxFx(const edm::ParameterSet& iConfig) {}
0009
0010 bool JetMatchingEWKFxFx::initAfterBeams() {
0011
0012
0013
0014
0015
0016 pTfirstSave = -1.;
0017 processSubsetSave.init("(eventProcess)", particleDataPtr);
0018 workEventJetSave.init("(workEventJet)", particleDataPtr);
0019
0020
0021 bool setMad = settingsPtr->flag("JetMatching:setMad");
0022
0023
0024 MadgraphPar par;
0025
0026 string parStr = infoPtr->header("MGRunCard");
0027 if (!parStr.empty()) {
0028 par.parse(parStr);
0029 par.printParams();
0030 }
0031
0032
0033 if (setMad) {
0034 if (par.haveParam("xqcut") && par.haveParam("maxjetflavor") && par.haveParam("alpsfact") &&
0035 par.haveParam("ickkw")) {
0036 settingsPtr->flag("JetMatching:merge", par.getParam("ickkw"));
0037 settingsPtr->parm("JetMatching:qCut", par.getParam("xqcut"));
0038 settingsPtr->mode("JetMatching:nQmatch", par.getParamAsInt("maxjetflavor"));
0039 settingsPtr->parm("JetMatching:clFact", clFact = par.getParam("alpsfact"));
0040 if (par.getParamAsInt("ickkw") == 0)
0041 errorMsg(
0042 "Error in JetMatchingMadgraph:init: "
0043 "Madgraph file parameters are not set for merging");
0044
0045
0046 } else {
0047 errorMsg(
0048 "Warning in JetMatchingMadgraph:init: "
0049 "Madgraph merging parameters not found");
0050 if (!par.haveParam("xqcut"))
0051 errorMsg(
0052 "Warning in "
0053 "JetMatchingMadgraph:init: No xqcut");
0054 if (!par.haveParam("ickkw"))
0055 errorMsg(
0056 "Warning in "
0057 "JetMatchingMadgraph:init: No ickkw");
0058 if (!par.haveParam("maxjetflavor"))
0059 errorMsg(
0060 "Warning in "
0061 "JetMatchingMadgraph:init: No maxjetflavor");
0062 if (!par.haveParam("alpsfact"))
0063 errorMsg(
0064 "Warning in "
0065 "JetMatchingMadgraph:init: No alpsfact");
0066 }
0067 }
0068
0069
0070 doFxFx = settingsPtr->flag("JetMatching:doFxFx");
0071 nPartonsNow = settingsPtr->mode("JetMatching:nPartonsNow");
0072 qCutME = settingsPtr->parm("JetMatching:qCutME");
0073 qCutMESq = pow(qCutME, 2);
0074
0075
0076 doMerge = settingsPtr->flag("JetMatching:merge");
0077 doShowerKt = settingsPtr->flag("JetMatching:doShowerKt");
0078 qCut = settingsPtr->parm("JetMatching:qCut");
0079 nQmatch = settingsPtr->mode("JetMatching:nQmatch");
0080 clFact = settingsPtr->parm("JetMatching:clFact");
0081
0082
0083 jetAlgorithm = settingsPtr->mode("JetMatching:jetAlgorithm");
0084 nJetMax = settingsPtr->mode("JetMatching:nJetMax");
0085 eTjetMin = settingsPtr->parm("JetMatching:eTjetMin");
0086 coneRadius = settingsPtr->parm("JetMatching:coneRadius");
0087 etaJetMax = settingsPtr->parm("JetMatching:etaJetMax");
0088 slowJetPower = settingsPtr->mode("JetMatching:slowJetPower");
0089
0090
0091 jetAllow = settingsPtr->mode("JetMatching:jetAllow");
0092 exclusiveMode = settingsPtr->mode("JetMatching:exclusive");
0093 qCutSq = pow(qCut, 2);
0094 etaJetMaxAlgo = etaJetMax;
0095
0096
0097
0098 if (!doMerge)
0099 return true;
0100
0101
0102 if (exclusiveMode == 2) {
0103
0104 if (nJetMax < 0) {
0105 errorMsg(
0106 "Warning in JetMatchingMadgraph:init: "
0107 "missing jet multiplicity information; running in exclusive mode");
0108 exclusiveMode = 1;
0109 }
0110 }
0111
0112
0113
0114
0115 jetAlgorithm = 2;
0116
0117 slowJet = new SlowJet(slowJetPower, coneRadius, eTjetMin, etaJetMaxAlgo, 2, 2, nullptr, false);
0118
0119
0120 slowJetHard = new SlowJet(slowJetPower, coneRadius, qCutME, etaJetMaxAlgo, 2, 2, nullptr, false);
0121
0122
0123 slowJetDJR = new SlowJet(slowJetPower, coneRadius, qCutME, etaJetMaxAlgo, 2, 2, nullptr, false);
0124
0125
0126 hjSlowJet = new HJSlowJet(slowJetPower, coneRadius, 0.0, 100.0, 1, 2, nullptr, false, true);
0127
0128
0129 eventProcessOrig.init("(eventProcessOrig)", particleDataPtr);
0130 eventProcess.init("(eventProcess)", particleDataPtr);
0131 workEventJet.init("(workEventJet)", particleDataPtr);
0132
0133
0134 if (MATCHINGDEBUG) {
0135 string jetStr = (jetAlgorithm == 1) ? "CellJet"
0136 : (slowJetPower == -1) ? "anti-kT"
0137 : (slowJetPower == 0) ? "C/A"
0138 : (slowJetPower == 1) ? "kT"
0139 : "unknown";
0140 string modeStr = (exclusiveMode) ? "exclusive" : "inclusive";
0141 cout << endl
0142 << " *----- Madgraph matching parameters -----*" << endl
0143 << " | qCut | " << setw(14) << qCut << " |" << endl
0144 << " | nQmatch | " << setw(14) << nQmatch << " |" << endl
0145 << " | clFact | " << setw(14) << clFact << " |" << endl
0146 << " | Jet algorithm | " << setw(14) << jetStr << " |" << endl
0147 << " | eTjetMin | " << setw(14) << eTjetMin << " |" << endl
0148 << " | etaJetMax | " << setw(14) << etaJetMax << " |" << endl
0149 << " | jetAllow | " << setw(14) << jetAllow << " |" << endl
0150 << " | Mode | " << setw(14) << modeStr << " |" << endl
0151 << " *-----------------------------------------*" << endl;
0152 }
0153 return true;
0154 }
0155
0156 bool JetMatchingEWKFxFx::doVetoPartonLevelEarly(const Event& event) {
0157
0158 sortIncomingProcess(event);
0159
0160
0161
0162 if (doShowerKt)
0163 return false;
0164
0165
0166 if (MATCHINGDEBUG) {
0167
0168 cout << endl << "-------- Begin Madgraph Debug --------" << endl;
0169
0170 cout << endl << "Original incoming process:";
0171 eventProcessOrig.list();
0172
0173 cout << endl << "Final-state incoming process:";
0174 eventProcess.list();
0175
0176 for (size_t i = 0; i < typeIdx[0].size(); i++)
0177 cout << ((i == 0) ? "Light jets: " : ", ") << setw(3) << typeIdx[0][i];
0178 if (typeIdx[0].empty())
0179 cout << "Light jets: None";
0180
0181 for (size_t i = 0; i < typeIdx[1].size(); i++)
0182 cout << ((i == 0) ? "\nHeavy jets: " : ", ") << setw(3) << typeIdx[1][i];
0183 for (size_t i = 0; i < typeIdx[2].size(); i++)
0184 cout << ((i == 0) ? "\nOther: " : ", ") << setw(3) << typeIdx[2][i];
0185
0186 cout << endl << endl << "Event:";
0187 event.list();
0188
0189 cout << endl << "Work event:";
0190 workEvent.list();
0191 }
0192
0193
0194 int iTypeEnd = (typeIdx[2].empty()) ? 2 : 3;
0195 for (int iType = 0; iType < iTypeEnd; iType++) {
0196
0197
0198 jetAlgorithmInput(event, iType);
0199
0200
0201 if (MATCHINGDEBUG) {
0202
0203 cout << endl << "Jet algorithm event (iType = " << iType << "):";
0204 workEventJet.list();
0205 }
0206
0207
0208
0209 runJetAlgorithm();
0210
0211
0212 if (matchPartonsToJets(iType) == true) {
0213
0214 if (MATCHINGDEBUG) {
0215 cout << endl
0216 << "Event vetoed"
0217 << "---------- End MLM Debug ----------" << endl;
0218 }
0219 return true;
0220 }
0221 }
0222
0223
0224 if (MATCHINGDEBUG) {
0225 cout << endl
0226 << "Event accepted"
0227 << "---------- End MLM Debug ----------" << endl;
0228 }
0229
0230 return false;
0231 }
0232
0233 void JetMatchingEWKFxFx::sortIncomingProcess(const Event& event) {
0234 omitResonanceDecays(eventProcessOrig, true);
0235 clearDJR();
0236 clear_nMEpartons();
0237
0238
0239 eventProcess = workEvent;
0240
0241 for (int i = 0; i < 3; i++) {
0242 typeIdx[i].clear();
0243 typeSet[i].clear();
0244 origTypeIdx[i].clear();
0245 }
0246
0247 for (int i = 0; i < eventProcess.size(); i++) {
0248
0249 if (!eventProcess[i].isFinal())
0250 continue;
0251 int idx = -1;
0252 int orig_idx = -1;
0253
0254
0255 if (eventProcess[i].isGluon() || (eventProcess[i].idAbs() <= nQmatch)) {
0256 orig_idx = 0;
0257 if (doFxFx) {
0258
0259
0260 idx = (trunc(1000. * eventProcess[i].scale()) == trunc(1000. * infoPtr->scalup())) ? 0 : 2;
0261 } else {
0262
0263
0264 idx = (eventProcess[i].scale() < 1.999 * sqrt(infoPtr->eA() * infoPtr->eB())) ? 0 : 2;
0265 }
0266 }
0267
0268
0269 else if (eventProcess[i].idAbs() > nQmatch && eventProcess[i].idAbs() <= ID_TOP) {
0270 idx = 1;
0271 orig_idx = 1;
0272
0273 } else if (eventProcess[i].colType() != 0 && eventProcess[i].idAbs() > ID_TOP) {
0274 idx = 1;
0275 orig_idx = 1;
0276 }
0277 if (idx < 0)
0278 continue;
0279
0280 typeIdx[idx].push_back(i);
0281 typeSet[idx].insert(eventProcess[i].daughter1());
0282 origTypeIdx[orig_idx].push_back(i);
0283 }
0284
0285 if (exclusiveMode == 2) {
0286
0287 int nParton = origTypeIdx[0].size();
0288 exclusive = (nParton == nJetMax) ? false : true;
0289
0290
0291 } else {
0292 exclusive = (exclusiveMode == 0) ? false : true;
0293 }
0294
0295
0296
0297 subEvent(event);
0298
0299
0300 int nParton = typeIdx[0].size();
0301 processSubsetSave.clear();
0302 for (int i = 0; i < nParton; ++i)
0303 processSubsetSave.append(eventProcess[typeIdx[0][i]]);
0304 }
0305
0306 bool JetMatchingEWKFxFx::doVetoProcessLevel(Event& process) {
0307 eventProcessOrig = process;
0308
0309
0310
0311
0312
0313
0314 sortIncomingProcess(process);
0315
0316
0317 if (!doFxFx && int(typeIdx[0].size()) > nJetMax)
0318 return true;
0319 if (doFxFx && npNLO() < nJetMax && int(typeIdx[0].size()) > nJetMax)
0320 return true;
0321
0322
0323 return false;
0324 }
0325
0326 void JetMatchingEWKFxFx::jetAlgorithmInput(const Event& event, int iType) {
0327
0328 workEventJet = workEvent;
0329
0330
0331 for (int i = 0; i < workEventJet.size(); ++i) {
0332 if (!workEventJet[i].isFinal())
0333 continue;
0334
0335
0336 if (jetAllow == 1) {
0337
0338 if (workEventJet[i].colType() == 0) {
0339 workEventJet[i].statusNeg();
0340 continue;
0341 }
0342 }
0343
0344 int idx = workEventJet[i].daughter1();
0345
0346
0347 while (true) {
0348
0349 if (iType == 0) {
0350
0351 if (typeSet[1].find(idx) != typeSet[1].end() || typeSet[2].find(idx) != typeSet[2].end()) {
0352 workEventJet[i].statusNeg();
0353 break;
0354 }
0355
0356
0357 if (idx == 0)
0358 break;
0359
0360 idx = event[idx].mother1();
0361
0362
0363 } else if (iType == 1) {
0364
0365 if (typeSet[1].find(idx) != typeSet[1].end())
0366 break;
0367
0368
0369
0370 if (idx == 0) {
0371 workEventJet[i].statusNeg();
0372 break;
0373 }
0374
0375
0376 idx = event[idx].mother1();
0377
0378
0379 } else if (iType == 2) {
0380
0381 if (typeSet[2].find(idx) != typeSet[2].end())
0382 break;
0383
0384
0385
0386 if (idx == 0) {
0387 workEventJet[i].statusNeg();
0388 break;
0389 }
0390
0391
0392 idx = event[idx].mother1();
0393 }
0394 }
0395 }
0396 }
0397
0398 void JetMatchingEWKFxFx::runJetAlgorithm() { ; }
0399
0400 bool JetMatchingEWKFxFx::doVetoStep(int iPos, int nISR, int nFSR, const Event& event) {
0401
0402 if (!doShowerKt)
0403 return false;
0404
0405
0406 if (nISR + nFSR > 1)
0407 return false;
0408
0409
0410 if (iPos == 5)
0411 return false;
0412
0413
0414
0415 sortIncomingProcess(event);
0416
0417
0418 double pTfirst = 0.;
0419
0420
0421 vector<int> weakBosons;
0422 for (int i = 0; i < event.size(); i++) {
0423 if (event[i].id() == 22 && event[i].id() == 23 && event[i].idAbs() == 24)
0424 weakBosons.push_back(i);
0425 }
0426
0427 for (int i = workEvent.size() - 1; i > 0; --i) {
0428 if (workEvent[i].isFinal() && workEvent[i].colType() != 0 &&
0429 (workEvent[i].statusAbs() == 43 || workEvent[i].statusAbs() == 51)) {
0430
0431
0432
0433 bool QCDemission = true;
0434
0435
0436 int iPosOld = workEvent[i].daughter1();
0437 for (int j = 0; i < int(weakBosons.size()); ++i)
0438 if (event[iPosOld].isAncestor(j)) {
0439 QCDemission = false;
0440 break;
0441 }
0442
0443 if (QCDemission) {
0444 pTfirst = workEvent[i].pT();
0445 break;
0446 }
0447 }
0448 }
0449
0450
0451 pTfirstSave = pTfirst;
0452
0453
0454
0455
0456 if (doShowerKtVeto(pTfirst))
0457 return true;
0458
0459
0460 return false;
0461 }
0462
0463 void JetMatchingEWKFxFx::setDJR(const Event& event) {
0464
0465 clearDJR();
0466 vector<double> result;
0467
0468
0469 if (!slowJetDJR->setup(event)) {
0470 errorMsg(
0471 "Warning in JetMatchingMadgraph:setDJR"
0472 ": the SlowJet algorithm failed on setup");
0473 return;
0474 }
0475
0476
0477 while (slowJetDJR->sizeAll() - slowJetDJR->sizeJet() > 0) {
0478
0479 result.push_back(sqrt(slowJetDJR->dNext()));
0480
0481 slowJetDJR->doStep();
0482 }
0483
0484
0485 for (int i = int(result.size()) - 1; i >= 0; --i)
0486 DJR.push_back(result[i]);
0487 }
0488
0489 bool JetMatchingEWKFxFx::matchPartonsToJets(int iType) {
0490
0491
0492 if (iType == 0) {
0493
0494
0495 setDJR(workEventJet);
0496 set_nMEpartons(origTypeIdx[0].size(), typeIdx[0].size());
0497
0498 return (matchPartonsToJetsLight() > 0);
0499 } else if (iType == 1) {
0500 return (matchPartonsToJetsHeavy() > 0);
0501 } else {
0502 return (matchPartonsToJetsOther() > 0);
0503 }
0504 }
0505
0506 int JetMatchingEWKFxFx::matchPartonsToJetsLight() {
0507
0508 workEventJetSave = workEventJet;
0509
0510
0511
0512
0513 int nParton = typeIdx[0].size();
0514
0515
0516 if (!slowJet->setup(workEventJet)) {
0517 errorMsg(
0518 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0519 "Light: the SlowJet algorithm failed on setup");
0520 return NONE;
0521 }
0522 double localQcutSq = qCutSq;
0523 double dOld = 0.0;
0524
0525 while (slowJet->sizeAll() - slowJet->sizeJet() > 0) {
0526 if (slowJet->dNext() > localQcutSq)
0527 break;
0528 dOld = slowJet->dNext();
0529 slowJet->doStep();
0530 }
0531 int nJets = slowJet->sizeJet();
0532 int nClus = slowJet->sizeAll();
0533
0534
0535 if (MATCHINGDEBUG)
0536 slowJet->list(true);
0537
0538
0539 int nCLjets = nClus - nJets;
0540
0541
0542
0543
0544
0545 int nRequested = (doFxFx && !(npNLO() == nJetMax && npNLO() == (int)(typeIdx[2].size() - 1)))
0546 ? npNLO() - typeIdx[2].size()
0547 : nParton;
0548
0549
0550
0551
0552 if (doFxFx && npNLO() < nJetMax && !typeIdx[2].empty() && npNLO() == (int)(typeIdx[2].size() - 1)) {
0553 return MORE_JETS;
0554 }
0555
0556
0557
0558
0559
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571
0572
0573
0574 if (nCLjets < nRequested)
0575 return LESS_JETS;
0576
0577
0578 if (exclusive && !doFxFx) {
0579 if (nCLjets > nRequested)
0580 return MORE_JETS;
0581 } else {
0582
0583
0584
0585
0586
0587
0588
0589 if (doFxFx && npNLO() < nJetMax && nCLjets > nRequested)
0590 return MORE_JETS;
0591
0592
0593
0594
0595
0596 if (!slowJet->setup(workEventJet)) {
0597 errorMsg(
0598 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0599 "Light: the SlowJet algorithm failed on setup");
0600 return NONE;
0601 }
0602
0603
0604
0605 if (doFxFx) {
0606 while (slowJet->sizeAll() - slowJet->sizeJet() > 0) {
0607 if (slowJet->dNext() > localQcutSq)
0608 break;
0609 slowJet->doStep();
0610 }
0611
0612
0613 } else {
0614 while (slowJet->sizeAll() - slowJet->sizeJet() > nParton)
0615 slowJet->doStep();
0616 }
0617
0618
0619
0620 localQcutSq = dOld;
0621 if (clFact >= 0. && nParton > 0) {
0622 vector<double> partonPt;
0623 partonPt.reserve(nParton);
0624 for (int i = 0; i < nParton; ++i)
0625 partonPt.push_back(eventProcess[typeIdx[0][i]].pT2());
0626 sort(partonPt.begin(), partonPt.end());
0627 localQcutSq = max(qCutSq, partonPt[0]);
0628 }
0629 nJets = slowJet->sizeJet();
0630 nClus = slowJet->sizeAll();
0631 }
0632
0633 if (clFact != 0.)
0634 localQcutSq *= pow2(clFact);
0635
0636 Event tempEvent;
0637 tempEvent.init("(tempEvent)", particleDataPtr);
0638 int nPass = 0;
0639 double pTminEstimate = -1.;
0640
0641
0642
0643 for (int i = nJets; i < nClus; ++i) {
0644 tempEvent.append(
0645 ID_GLUON, 98, 0, 0, 0, 0, 0, 0, slowJet->p(i).px(), slowJet->p(i).py(), slowJet->p(i).pz(), slowJet->p(i).e());
0646 ++nPass;
0647 pTminEstimate = max(pTminEstimate, slowJet->pT(i));
0648 if (nPass == nRequested)
0649 break;
0650 }
0651
0652 int tempSize = tempEvent.size();
0653
0654 vector<bool> jetAssigned;
0655 jetAssigned.assign(tempSize, false);
0656
0657
0658
0659 vector<vector<bool> > partonMatchesJet;
0660 partonMatchesJet.reserve(nParton);
0661 for (int i = 0; i < nParton; ++i)
0662 partonMatchesJet.push_back(vector<bool>(tempEvent.size(), false));
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676 int iNow = 0;
0677 int nMatched = 0;
0678
0679 while (doFxFx && iNow < tempSize) {
0680
0681 Event tempEventJet;
0682 tempEventJet.init("(tempEventJet)", particleDataPtr);
0683 for (int i = 0; i < nParton; ++i) {
0684
0685
0686
0687
0688
0689 tempEventJet.clear();
0690 tempEventJet.append(ID_GLUON,
0691 98,
0692 0,
0693 0,
0694 0,
0695 0,
0696 0,
0697 0,
0698 tempEvent[iNow].px(),
0699 tempEvent[iNow].py(),
0700 tempEvent[iNow].pz(),
0701 tempEvent[iNow].e());
0702
0703 Vec4 pIn = eventProcess[typeIdx[0][i]].p();
0704 tempEventJet.append(ID_GLUON, 99, 0, 0, 0, 0, 0, 0, pIn.px(), pIn.py(), pIn.pz(), pIn.e());
0705
0706
0707 if (!slowJet->setup(tempEventJet)) {
0708 errorMsg(
0709 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0710 "Light: the SlowJet algorithm failed on setup");
0711 return NONE;
0712 }
0713
0714
0715 if (slowJet->iNext() == tempEventJet.size() - 1 && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq) {
0716 jetAssigned[iNow] = true;
0717 partonMatchesJet[i][iNow] = true;
0718 }
0719 }
0720
0721
0722 if (jetAssigned[iNow])
0723 nMatched++;
0724
0725
0726 ++iNow;
0727 }
0728
0729 if (doFxFx) {
0730
0731
0732
0733
0734
0735
0736
0737 if (npNLO() < nJetMax && nMatched != nRequested)
0738 return UNMATCHED_PARTON;
0739 if (npNLO() == nJetMax && nMatched < nRequested)
0740 return UNMATCHED_PARTON;
0741 }
0742
0743
0744
0745
0746
0747
0748
0749 iNow = 0;
0750 while (!doFxFx && iNow < nParton) {
0751 Event tempEventJet;
0752 tempEventJet.init("(tempEventJet)", particleDataPtr);
0753 for (int i = 0; i < tempSize; ++i) {
0754 if (jetAssigned[i])
0755 continue;
0756 Vec4 pIn = tempEvent[i].p();
0757
0758 tempEventJet.append(ID_GLUON, 98, 0, 0, 0, 0, 0, 0, pIn.px(), pIn.py(), pIn.pz(), pIn.e());
0759 }
0760
0761 Vec4 pIn = eventProcess[typeIdx[0][iNow]].p();
0762
0763 tempEventJet.append(ID_GLUON, 99, 0, 0, 0, 0, 0, 0, pIn.px(), pIn.py(), pIn.pz(), pIn.e());
0764 if (!slowJet->setup(tempEventJet)) {
0765 errorMsg(
0766 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0767 "Light: the SlowJet algorithm failed on setup");
0768 return NONE;
0769 }
0770
0771
0772 if (slowJet->iNext() == tempEventJet.size() - 1 && slowJet->jNext() > -1 && slowJet->dNext() < localQcutSq) {
0773 int iKnt = -1;
0774 for (int i = 0; i != tempSize; ++i) {
0775 if (jetAssigned[i])
0776 continue;
0777 ++iKnt;
0778
0779 if (iKnt == slowJet->jNext())
0780 jetAssigned[i] = true;
0781 }
0782 } else {
0783 return UNMATCHED_PARTON;
0784 }
0785 ++iNow;
0786 }
0787
0788
0789
0790 if (nParton > 0 && pTminEstimate > 0)
0791 eTpTlightMin = pTminEstimate;
0792 else
0793 eTpTlightMin = -1.;
0794
0795
0796 setDJR(workEventJet);
0797
0798
0799 return NONE;
0800 }
0801
0802 int JetMatchingEWKFxFx::matchPartonsToJetsHeavy() {
0803
0804
0805
0806
0807
0808
0809
0810
0811
0812
0813 int nParton = typeIdx[1].size();
0814
0815 Event tempEventJet(workEventJet);
0816
0817 double scaleF(1.0);
0818
0819
0820
0821 for (int i = 0; i < nParton; ++i) {
0822 scaleF = eventProcessOrig[0].e() / workEventJet[typeIdx[1][i]].pT();
0823 tempEventJet[typeIdx[1][i]].rescale5(scaleF);
0824 }
0825
0826 if (!hjSlowJet->setup(tempEventJet)) {
0827 errorMsg(
0828 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0829 "Heavy: the SlowJet algorithm failed on setup");
0830 return NONE;
0831 }
0832
0833 while (hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0) {
0834 if (hjSlowJet->dNext() > qCutSq)
0835 break;
0836 hjSlowJet->doStep();
0837 }
0838
0839 int nCLjets(0);
0840
0841
0842 for (int idx = 0; idx < hjSlowJet->sizeAll(); ++idx) {
0843 if (hjSlowJet->pT(idx) > qCut)
0844 nCLjets++;
0845 }
0846
0847
0848 if (MATCHINGDEBUG)
0849 hjSlowJet->list(true);
0850
0851
0852
0853
0854 int nRequested = nParton;
0855
0856
0857 if (nCLjets < nRequested) {
0858 if (MATCHINGDEBUG)
0859 cout << "veto : hvy LESS_JETS " << endl;
0860 if (MATCHINGDEBUG)
0861 cout << "nCLjets = " << nCLjets << "; nRequest = " << nRequested << endl;
0862 return LESS_JETS;
0863 }
0864
0865
0866 if (exclusive) {
0867 if (nCLjets > nRequested) {
0868 if (MATCHINGDEBUG)
0869 cout << "veto : excl hvy MORE_JETS " << endl;
0870 return MORE_JETS;
0871 }
0872 }
0873
0874
0875 return NONE;
0876 }
0877
0878 int JetMatchingEWKFxFx::matchPartonsToJetsOther() {
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888
0889 int nParton = typeIdx[2].size();
0890
0891 Event tempEventJet(workEventJet);
0892
0893 double scaleF(1.0);
0894
0895
0896
0897 for (int i = 0; i < nParton; ++i) {
0898 scaleF = eventProcessOrig[0].e() / workEventJet[typeIdx[2][i]].pT();
0899 tempEventJet[typeIdx[2][i]].rescale5(scaleF);
0900 }
0901
0902 if (!hjSlowJet->setup(tempEventJet)) {
0903 errorMsg(
0904 "Warning in JetMatchingMadgraph:matchPartonsToJets"
0905 "Heavy: the SlowJet algorithm failed on setup");
0906 return NONE;
0907 }
0908
0909 while (hjSlowJet->sizeAll() - hjSlowJet->sizeJet() > 0) {
0910 if (hjSlowJet->dNext() > qCutSq)
0911 break;
0912 hjSlowJet->doStep();
0913 }
0914
0915 int nCLjets(0);
0916
0917
0918 for (int idx = 0; idx < hjSlowJet->sizeAll(); ++idx) {
0919 if (hjSlowJet->pT(idx) > qCut)
0920 nCLjets++;
0921 }
0922
0923
0924 if (MATCHINGDEBUG)
0925 hjSlowJet->list(true);
0926
0927
0928
0929
0930 int nRequested = nParton;
0931
0932
0933 if (nCLjets < nRequested) {
0934 if (MATCHINGDEBUG)
0935 cout << "veto : other LESS_JETS " << endl;
0936 if (MATCHINGDEBUG)
0937 cout << "nCLjets = " << nCLjets << "; nRequest = " << nRequested << endl;
0938 return LESS_JETS;
0939 }
0940
0941
0942 if (exclusive) {
0943 if (nCLjets > nRequested) {
0944 if (MATCHINGDEBUG)
0945 cout << "veto : excl other MORE_JETS" << endl;
0946 return MORE_JETS;
0947 }
0948 }
0949
0950
0951 return NONE;
0952 }
0953
0954 bool JetMatchingEWKFxFx::doShowerKtVeto(double pTfirst) {
0955
0956 if (!doShowerKt)
0957 return false;
0958
0959
0960 bool doVeto = false;
0961
0962
0963
0964 int nParton = typeIdx[0].size();
0965 double pTminME = 1e10;
0966 for (int i = 0; i < nParton; ++i)
0967 pTminME = min(pTminME, eventProcess[typeIdx[0][i]].pT());
0968
0969
0970 if (nParton > 0 && pow(pTminME, 2) < qCutSq)
0971 doVeto = true;
0972
0973
0974
0975 if (exclusive && pow(pTfirst, 2) > qCutSq) {
0976 doVeto = true;
0977
0978
0979 } else if (!exclusive && nParton > 0 && pTfirst > pTminME) {
0980 doVeto = true;
0981 }
0982
0983
0984 return doVeto;
0985 }
0986
0987
0988
0989
0990 int JetMatchingEWKFxFx::npNLO() {
0991 string npIn = infoPtr->getEventAttribute("npNLO", true);
0992 int np = !npIn.empty() ? std::atoi(npIn.c_str()) : -1;
0993 if (np < 0) {
0994 ;
0995 } else
0996 return np;
0997 return nPartonsNow;
0998 }