53 double lheHT = 0, lheHTIncoming = 0;
54 unsigned int lheNj = 0, lheNb = 0, lheNc = 0, lheNuds = 0, lheNglu = 0;
58 const auto& hepeup = lheProd.
hepeup();
59 const auto& pup = hepeup.
PUP;
60 int lep = -1, lepBar = -1, nu = -1, nuBar = -1;
61 std::vector<float> vals_pt;
62 std::vector<float> vals_eta;
63 std::vector<float> vals_phi;
64 std::vector<float> vals_mass;
65 std::vector<float> vals_pz;
66 std::vector<int> vals_pid;
67 std::vector<int> vals_status;
68 std::vector<int> vals_spin;
69 alphaS = hepeup.AQCDUP;
70 for (
unsigned int i = 0,
n = pup.size();
i <
n; ++
i) {
74 TLorentzVector
p4(pup[
i][0], pup[
i][1], pup[
i][2], pup[
i][3]);
75 vals_pid.push_back(hepeup.IDUP[
i]);
76 vals_spin.push_back(hepeup.SPINUP[
i]);
77 vals_status.push_back(
status);
80 vals_eta.push_back(0);
81 vals_phi.push_back(0);
82 vals_mass.push_back(0);
83 vals_pz.push_back(
p4.Pz());
85 vals_pt.push_back(
p4.Pt());
86 vals_eta.push_back(
p4.Eta());
87 vals_phi.push_back(
p4.Phi());
88 vals_mass.push_back(
p4.M());
92 if ((
status == 1) && ((idabs == 21) || (idabs > 0 && idabs < 7))) {
101 else if (idabs == 21)
104 double pt = std::hypot(pup[
i][0], pup[
i][1]);
107 hepeup.MOTHUP[
i].first - 1,
109 int mothIdxTwo =
std::max(hepeup.MOTHUP[
i].second - 1, 0);
110 int mothStatus = hepeup.ISTUP[mothIdx];
111 int mothStatusTwo = hepeup.ISTUP[mothIdxTwo];
112 bool hasIncomingAsMother = mothStatus < 0 || mothStatusTwo < 0;
113 if (hasIncomingAsMother)
115 }
else if (idabs == 12 || idabs == 14 || idabs == 16) {
116 (hepeup.IDUP[
i] > 0 ? nu : nuBar) =
i;
117 }
else if (idabs == 11 || idabs == 13 || idabs == 15) {
118 (hepeup.IDUP[
i] > 0 ? lep : lepBar) =
i;
121 std::pair<int, int>
v(0, 0);
122 if (lep != -1 && lepBar != -1)
123 v = std::make_pair(lep, lepBar);
124 else if (lep != -1 && nuBar != -1)
125 v = std::make_pair(lep, nuBar);
126 else if (nu != -1 && lepBar != -1)
127 v = std::make_pair(nu, lepBar);
128 else if (nu != -1 && nuBar != -1)
129 v = std::make_pair(nu, nuBar);
130 if (
v.first != -1 &&
v.second != -1) {
131 lheVpt = std::hypot(pup[
v.first][0] + pup[
v.second][0], pup[
v.first][1] + pup[
v.second][1]);
134 out.addColumnValue<uint8_t>(
"Njets", lheNj,
"Number of jets (partons) at LHE step");
135 out.addColumnValue<uint8_t>(
"Nb", lheNb,
"Number of b partons at LHE step");
136 out.addColumnValue<uint8_t>(
"Nc", lheNc,
"Number of c partons at LHE step");
137 out.addColumnValue<uint8_t>(
"Nuds", lheNuds,
"Number of u,d,s partons at LHE step");
138 out.addColumnValue<uint8_t>(
"Nglu", lheNglu,
"Number of gluon partons at LHE step");
139 out.addColumnValue<
float>(
"HT", lheHT,
"HT, scalar sum of parton pTs at LHE step");
140 out.addColumnValue<
float>(
141 "HTIncoming", lheHTIncoming,
"HT, scalar sum of parton pTs at LHE step, restricted to partons");
142 out.addColumnValue<
float>(
"Vpt", lheVpt,
"pT of the W or Z boson at LHE step");
143 out.addColumnValue<uint8_t>(
"NpNLO", lheProd.
npNLO(),
"number of partons at NLO");
144 out.addColumnValue<uint8_t>(
"NpLO", lheProd.
npLO(),
"number of partons at LO");
145 out.addColumnValue<
float>(
"AlphaS", alphaS,
"Per-event alphaS");
147 auto outPart = std::make_unique<nanoaod::FlatTable>(vals_pt.size(),
"LHEPart",
false);
148 outPart->addColumn<
float>(
"pt", vals_pt,
"Pt of LHE particles", this->
precision_);
149 outPart->addColumn<
float>(
"eta", vals_eta,
"Pseodorapidity of LHE particles", this->
precision_);
150 outPart->addColumn<
float>(
"phi", vals_phi,
"Phi of LHE particles", this->
precision_);
151 outPart->addColumn<
float>(
"mass", vals_mass,
"Mass of LHE particles", this->
precision_);
152 outPart->addColumn<
float>(
"incomingpz", vals_pz,
"Pz of incoming LHE particles", this->
precision_);
153 outPart->addColumn<
int>(
"pdgId", vals_pid,
"PDG ID of LHE particles");
154 outPart->addColumn<
int>(
"status", vals_status,
"LHE particle status; -1:incoming, 1:outgoing");
155 outPart->addColumn<
int>(
"spin", vals_spin,
"Spin of LHE particles");