File indexing completed on 2024-09-08 23:51:48
0001 #include "GeneratorInterface/Pythia8Interface/plugins/Py8toJetInput.h"
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/LHECommonBlocks.h"
0004 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0005 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0006
0007 const std::vector<fastjet::PseudoJet> Py8toJetInput::fillJetAlgoInput(const Event& event,
0008 const Event& workEvent,
0009 const lhef::LHEEvent* lhee,
0010 const std::vector<int>* typeIdx) {
0011 fJetInput.clear();
0012
0013 Event workEventJet = workEvent;
0014
0015 const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();
0016
0017 std::set<int> typeSet[3];
0018
0019
0020
0021
0022
0023 for (size_t i = 0; i < 3; i++) {
0024 typeSet[i].clear();
0025 for (size_t j = 0; j < typeIdx[i].size(); j++) {
0026
0027
0028
0029
0030
0031
0032 int pos = typeIdx[i][j];
0033 int shift = 3;
0034
0035 for (int ip = 2; ip < pos; ip++) {
0036
0037
0038
0039
0040 if (hepeup.MOTHUP[ip].first == hepeup.MOTHUP[ip].second) {
0041 shift--;
0042 }
0043 }
0044 pos += shift;
0045
0046 typeSet[i].insert(pos);
0047 }
0048 }
0049
0050
0051
0052
0053 int iType = 0;
0054 int jetAllow = 0;
0055
0056
0057
0058
0059 for (int i = 0; i < workEventJet.size(); ++i) {
0060 if (!workEventJet[i].isFinal())
0061 continue;
0062
0063
0064 if (jetAllow == 1) {
0065
0066
0067 int id = workEventJet[i].idAbs();
0068 if ((id >= 11 && id <= 16) || id == ID_TOP || id == ID_PHOTON) {
0069 workEventJet[i].statusNeg();
0070 continue;
0071 }
0072 }
0073
0074
0075 int idx = workEventJet[i].daughter1();
0076
0077
0078 while (true) {
0079
0080 if (iType == 0) {
0081
0082 if (typeSet[1].find(idx) != typeSet[1].end() || typeSet[2].find(idx) != typeSet[2].end()) {
0083 workEventJet[i].statusNeg();
0084 break;
0085 }
0086
0087
0088 if (idx == 0) {
0089 break;
0090 }
0091
0092 idx = event[idx].mother1();
0093
0094
0095 } else if (iType == 1) {
0096
0097 if (typeSet[1].find(idx) != typeSet[1].end())
0098 break;
0099
0100
0101
0102 if (idx == 0) {
0103 workEventJet[i].statusNeg();
0104 break;
0105 }
0106
0107
0108 idx = event[idx].mother1();
0109
0110 }
0111 }
0112 }
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141 for (int i = 0; i < workEventJet.size(); i++) {
0142
0143
0144 if (workEventJet[i].status() < 0)
0145 continue;
0146
0147
0148
0149
0150
0151 if (fabs(workEventJet[i].eta()) > fJetEtaMax)
0152 continue;
0153
0154
0155
0156
0157 fastjet::PseudoJet partTmp = workEventJet[i];
0158 fJetInput.push_back(partTmp);
0159 }
0160
0161 return fJetInput;
0162 }
0163
0164 int Py8toJetInput::getAncestor(int pos, const Event& fullEvent, const Event& workEvent) {
0165 int parentId = fullEvent[pos].mother1();
0166 int parentPrevId = 0;
0167 int counter = pos;
0168
0169 while (parentId > 0) {
0170 if (parentId == fullEvent[counter].mother2())
0171 {
0172 parentPrevId = parentId;
0173 counter = parentId;
0174 parentId = fullEvent[parentPrevId].mother1();
0175 continue;
0176 }
0177
0178
0179
0180
0181
0182 if ((parentId < parentPrevId) || parentId < fullEvent[counter].mother2())
0183 {
0184
0185 if (abs(fullEvent[counter].status()) == 22 || abs(fullEvent[counter].status()) == 23) {
0186
0187
0188 parentId = counter;
0189 break;
0190 } else {
0191 parentPrevId = parentId;
0192 parentId = fullEvent[parentPrevId].mother1();
0193 }
0194 } else if (parentId > parentPrevId ||
0195 parentId > pos)
0196 {
0197 parentId = -1;
0198 break;
0199 }
0200
0201
0202
0203 if (abs(fullEvent[parentId].status()) == 22 || abs(fullEvent[parentId].status()) == 23)
0204 {
0205 break;
0206 }
0207 if (abs(fullEvent[parentId].status()) < 22)
0208 {
0209 parentId = -1;
0210 break;
0211 }
0212 }
0213
0214 return parentId;
0215 }
0216
0217 #include "HepMC/HEPEVT_Wrapper.h"
0218 #include <cassert>
0219
0220 const std::vector<fastjet::PseudoJet> Py8toJetInputHEPEVT::fillJetAlgoInput(const Event& event,
0221 const Event& workEvent,
0222 const lhef::LHEEvent* lhee,
0223 const std::vector<int>*) {
0224 fJetInput.clear();
0225
0226 HepMC::HEPEVT_Wrapper::zero_everything();
0227
0228
0229
0230 std::vector<int> Py8PartonIdx;
0231 Py8PartonIdx.clear();
0232 std::vector<int> HEPEVTPartonIdx;
0233 HEPEVTPartonIdx.clear();
0234
0235
0236
0237 int index = 0;
0238
0239 int Py8PartonCounter = 0;
0240 int HEPEVTPartonCounter = 0;
0241
0242
0243
0244 for (int iprt = 1; iprt < event.size(); iprt++) {
0245 const Particle& part = event[iprt];
0246 if (abs(part.status()) < 22)
0247 continue;
0248
0249
0250 Py8PartonCounter = iprt;
0251 break;
0252 }
0253
0254 const lhef::HEPEUP& hepeup = *lhee->getHEPEUP();
0255
0256 for (int iprt = 2; iprt < hepeup.NUP; iprt++) {
0257 index++;
0258 HepMC::HEPEVT_Wrapper::set_id(index, hepeup.IDUP[iprt]);
0259 HepMC::HEPEVT_Wrapper::set_status(index, 2);
0260 HepMC::HEPEVT_Wrapper::set_momentum(
0261 index, hepeup.PUP[iprt][0], hepeup.PUP[iprt][1], hepeup.PUP[iprt][2], hepeup.PUP[iprt][4]);
0262 HepMC::HEPEVT_Wrapper::set_mass(index, hepeup.PUP[iprt][4]);
0263
0264 HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);
0265 HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);
0266 if (hepeup.MOTHUP[iprt].first > 2 &&
0267 hepeup.MOTHUP[iprt].second > 2)
0268 {
0269 HEPEVTPartonCounter++;
0270 continue;
0271 }
0272 Py8PartonIdx.push_back(Py8PartonCounter);
0273 Py8PartonCounter++;
0274 HEPEVTPartonIdx.push_back(HEPEVTPartonCounter);
0275 HEPEVTPartonCounter++;
0276 }
0277
0278 HepMC::HEPEVT_Wrapper::set_number_entries(index);
0279
0280
0281
0282
0283 for (int iprt = 1; iprt < workEvent.size(); iprt++)
0284 {
0285 const Particle& part = workEvent[iprt];
0286
0287
0288 if (part.status() < 51)
0289 continue;
0290 index++;
0291 HepMC::HEPEVT_Wrapper::set_id(index, part.id());
0292
0293
0294 HepMC::HEPEVT_Wrapper::set_status(index, 1);
0295 HepMC::HEPEVT_Wrapper::set_momentum(index, part.px(), part.py(), part.pz(), part.e());
0296 HepMC::HEPEVT_Wrapper::set_mass(index, part.m());
0297 HepMC::HEPEVT_Wrapper::set_position(index, part.xProd(), part.yProd(), part.zProd(), part.tProd());
0298 HepMC::HEPEVT_Wrapper::set_parents(index, 0, 0);
0299
0300 HepMC::HEPEVT_Wrapper::set_children(index, 0, 0);
0301
0302
0303
0304 int parentId = getAncestor(part.daughter1(), event, workEvent);
0305
0306 if (parentId <= 0)
0307 continue;
0308
0309 for (int idx = 0; idx < (int)Py8PartonIdx.size(); idx++) {
0310 if (parentId == Py8PartonIdx[idx]) {
0311 int idx1 = HEPEVTPartonIdx[idx];
0312 HepMC::HEPEVT_Wrapper::set_parents(index, idx1 + 1, idx1 + 1);
0313 break;
0314 }
0315 }
0316 }
0317
0318 HepMC::HEPEVT_Wrapper::set_number_entries(index);
0319
0320
0321
0322
0323
0324 return fJetInput;
0325 }