10 std::vector<pat::Jet> topJets;
12 topJets.push_back(TS.
getHadp());
13 topJets.push_back(TS.
getHadq());
14 topJets.push_back(TS.
getHadj());
15 topJets.push_back(TS.
getHadk());
16 topJets.push_back(TS.
getHadb());
20 TLorentzVector *Hadp =
new TLorentzVector();
21 TLorentzVector *Hadq =
new TLorentzVector();
22 TLorentzVector *Hadj =
new TLorentzVector();
23 TLorentzVector *Hadk =
new TLorentzVector();
24 TLorentzVector *Hadb =
new TLorentzVector();
25 TLorentzVector *Hadbbar =
new TLorentzVector();
26 Hadp->SetPxPyPzE(topJets[0].
px(), topJets[0].
py(), topJets[0].pz(), topJets[3].
energy());
27 Hadq->SetPxPyPzE(topJets[1].
px(), topJets[1].
py(), topJets[1].pz(), topJets[2].
energy());
28 Hadj->SetPxPyPzE(topJets[2].
px(), topJets[2].
py(), topJets[2].pz(), topJets[2].
energy());
29 Hadk->SetPxPyPzE(topJets[3].
px(), topJets[3].
py(), topJets[3].pz(), topJets[3].
energy());
30 Hadb->SetPxPyPzE(topJets[4].
px(), topJets[4].
py(), topJets[4].pz(), topJets[1].
energy());
31 Hadbbar->SetPxPyPzE(topJets[4].
px(), topJets[4].
py(), topJets[4].pz(), topJets[1].
energy());
37 double EtSum = topJets[5].et() + topJets[5].et();
38 double Obs1 = (EtSum > 0 ? EtSum : -1);
47 double BGap = topJets[1].bDiscriminator(
"trackCountingJetTags") - topJets[2].bDiscriminator(
"trackCountingJetTags");
48 double Obs2 = (BGap > 0 ?
log(BGap) : -1);
52 double N = 0,
D = 0, C_tmp = 0,
C = 1000;
57 for (
unsigned int i = 0;
i < 6;
i++) {
58 D += fabs(topJets[
i].
pt());
63 for (
unsigned int i = 0;
i < 360;
i++) {
64 nx =
cos((2 *
PI / 360) *
i);
65 ny =
sin((2 *
PI / 360) *
i);
68 for (
unsigned int i = 0;
i < 4;
i++) {
69 N += fabs(topJets[
i].
px() * nx + topJets[
i].
py() * ny + topJets[
i].pz() * nz);
77 double Obs3 = (
C != 1000 ?
C : -1);
82 for (
unsigned int i = 0;
i < 6;
i++) {
83 HT += topJets[
i].et();
86 double Obs4 = (
HT != 0 ?
HT : -1);
92 for (
unsigned int i = 0;
i < 6;
i++) {
93 pjets += topJets[
i].p4();
96 double Obs5 = (
MT > 0 ?
MT : -1);
103 double px1 = topJets[2].px();
104 double px2 = topJets[3].px();
105 double py1 = topJets[2].py();
106 double py2 = topJets[3].py();
107 double pz1 = topJets[2].pz();
108 double pz2 = topJets[3].pz();
109 double E1 = topJets[2].energy();
110 double E2 = topJets[3].energy();
111 double px3 = topJets[4].px();
112 double px4 = topJets[5].px();
113 double py3 = topJets[4].py();
114 double py4 = topJets[5].py();
115 double pz3 = topJets[4].pz();
116 double pz4 = topJets[5].pz();
117 double E3 = topJets[4].energy();
118 double E4 = topJets[5].energy();
120 TLorentzVector *LightJet1 =
new TLorentzVector();
121 TLorentzVector *LightJet2 =
new TLorentzVector();
122 TLorentzVector *LightJet3 =
new TLorentzVector();
123 TLorentzVector *LightJet4 =
new TLorentzVector();
125 LightJet1->SetPxPyPzE(px1, py1, pz1, E1);
127 LightJet2->SetPxPyPzE(px2, py2, pz2, E2);
129 LightJet3->SetPxPyPzE(px3, py3, pz3, E3);
131 LightJet4->SetPxPyPzE(px4, py4, pz4, E4);
134 double CosTheta1 =
cos(LightJet2->Angle(LightJet1->Vect()));
135 double CosTheta2 =
cos(LightJet4->Angle(LightJet3->Vect()));
137 double Obs6 = (-1 < CosTheta1 ? CosTheta1 : -2);
139 double Obs7 = (-1 < CosTheta2 ? CosTheta2 : -2);
151 double BjetsBdiscSum =
152 topJets[0].bDiscriminator(
"trackCountingJetTags") + topJets[1].bDiscriminator(
"trackCountingJetTags");
153 double Ljets1BdiscSum =
154 topJets[2].bDiscriminator(
"trackCountingJetTags") + topJets[3].bDiscriminator(
"trackCountingJetTags");
155 double Ljets2BdiscSum =
156 topJets[4].bDiscriminator(
"trackCountingJetTags") + topJets[5].bDiscriminator(
"trackCountingJetTags");
160 double Obs8 = (Ljets1BdiscSum != 0 ?
log(BjetsBdiscSum / Ljets1BdiscSum) : -1);
162 double Obs9 = (Ljets2BdiscSum != 0 ?
log(BjetsBdiscSum / Ljets2BdiscSum) : -1);
164 double Obs10 = (BGap > 0 ?
log(BjetsBdiscSum * BGap) : -1);
168 double Obs11 = (
HT != 0 ? (
HT - EtSum) /
HT : -1);
180 double P2 = PX2 + PY2 + PZ2;
182 double PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Hadj->Px() * Hadj->Py() +
183 Hadk->Px() * Hadk->Py() + Hadbbar->Px() * Hadbbar->Py();
184 double PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Hadj->Px() * Hadj->Pz() +
185 Hadk->Px() * Hadk->Pz() + Hadbbar->Px() * Hadbbar->Pz();
186 double PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Hadj->Py() * Hadj->Pz() +
187 Hadk->Py() * Hadk->Pz() + Hadbbar->Py() * Hadbbar->Pz();
199 TMatrixDSymEigen pTensor(
Matrix);
201 std::vector<double> EigValues;
203 for (
int i = 0;
i < 3;
i++) {
204 EigValues.push_back(pTensor.GetEigenValues()[
i]);
209 double Sphericity = 1.5 * (EigValues[1] + EigValues[2]);
210 double Aplanarity = 1.5 * EigValues[2];
219 TLorentzVector *TtbarSystem =
new TLorentzVector();
220 TtbarSystem->SetPx(Hadp->Px() + Hadq->Px() + Hadb->Px() + Hadj->Px() + Hadk->Px() + Hadbbar->Px());
221 TtbarSystem->SetPy(Hadp->Py() + Hadq->Py() + Hadb->Py() + Hadj->Py() + Hadk->Py() + Hadbbar->Py());
222 TtbarSystem->SetPz(Hadp->Pz() + Hadq->Pz() + Hadb->Pz() + Hadj->Pz() + Hadk->Pz() + Hadbbar->Pz());
223 TtbarSystem->SetE(Hadp->E() + Hadq->E() + Hadb->E() + Hadj->E() + Hadk->E() + Hadbbar->E());
225 TVector3 BoostBackToCM = -(TtbarSystem->BoostVector());
226 Hadp->Boost(BoostBackToCM);
227 Hadq->Boost(BoostBackToCM);
228 Hadb->Boost(BoostBackToCM);
229 Hadj->Boost(BoostBackToCM);
230 Hadk->Boost(BoostBackToCM);
231 Hadbbar->Boost(BoostBackToCM);
240 double BOOST_P2 = BOOST_PX2 + BOOST_PY2 + BOOST_PZ2;
242 double BOOST_PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
243 Hadj->Px() * Hadj->Py() + Hadk->Px() * Hadk->Py() + Hadbbar->Px() * Hadbbar->Py();
244 double BOOST_PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
245 Hadj->Px() * Hadj->Pz() + Hadk->Px() * Hadk->Pz() + Hadbbar->Px() * Hadbbar->Pz();
246 double BOOST_PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
247 Hadj->Py() * Hadj->Pz() + Hadk->Py() * Hadk->Pz() + Hadbbar->Py() * Hadbbar->Pz();
249 TMatrixDSym BOOST_Matrix(3);
251 BOOST_Matrix(0, 0) = BOOST_PX2 / BOOST_P2;
252 BOOST_Matrix(0, 1) = BOOST_PXY / BOOST_P2;
253 BOOST_Matrix(0, 2) = BOOST_PXZ / BOOST_P2;
254 BOOST_Matrix(1, 0) = BOOST_PXY / BOOST_P2;
255 BOOST_Matrix(1, 1) = BOOST_PY2 / BOOST_P2;
256 BOOST_Matrix(1, 2) = BOOST_PYZ / BOOST_P2;
257 BOOST_Matrix(2, 0) = BOOST_PXZ / BOOST_P2;
258 BOOST_Matrix(2, 1) = BOOST_PYZ / BOOST_P2;
259 BOOST_Matrix(2, 2) = BOOST_PZ2 / BOOST_P2;
261 TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
263 std::vector<double> BOOST_EigValues;
264 BOOST_EigValues.clear();
265 for (
int i = 0;
i < 3;
i++) {
266 BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[
i]);
271 double BOOST_Sphericity = 1.5 * (BOOST_EigValues[1] + BOOST_EigValues[2]);
272 double BOOST_Aplanarity = 1.5 * BOOST_EigValues[2];
274 double Obs14 = (
edm::isNotFinite(BOOST_Sphericity) ? -1 : BOOST_Sphericity);
277 double Obs15 = (
edm::isNotFinite(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity);
282 for (
unsigned int i = 0;
i < 6;
i++) {
283 H += topJets[
i].energy();
285 double Obs16 = (
H != 0 ?
HT /
H : -1);
289 double Obs17 = (BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets1BdiscSum) : -1);
291 double Obs18 = (BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets2BdiscSum) : -1);
295 double Obs19 = (
HT != 0 ? EtSum /
HT : -1);
constexpr bool isNotFinite(T x)
Sin< T >::type sin(const T &t)
void setLRSignalEvtObservables(const std::vector< std::pair< unsigned int, double > > &varval)
std::vector< std::pair< unsigned int, double > > evtselectVarVal
CompareBdisc BdiscComparator
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Cos< T >::type cos(const T &t)
DecomposeProduct< arg, typename Div::arg > D
pat::Jet getHadbbar() const
CompareDouble dComparator
portabletest::Matrix Matrix
Power< A, B >::type pow(const A &a, const B &b)