File indexing completed on 2024-04-06 12:31:16
0001 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003
0004 TtHadLRSignalSelObservables::TtHadLRSignalSelObservables() {}
0005
0006 TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables() {}
0007
0008 void TtHadLRSignalSelObservables::operator()(TtHadEvtSolution &TS) {
0009 evtselectVarVal.clear();
0010 std::vector<pat::Jet> topJets;
0011 topJets.clear();
0012 topJets.push_back(TS.getHadp());
0013 topJets.push_back(TS.getHadq());
0014 topJets.push_back(TS.getHadj());
0015 topJets.push_back(TS.getHadk());
0016 topJets.push_back(TS.getHadb());
0017 topJets.push_back(TS.getHadbbar());
0018
0019
0020 TLorentzVector *Hadp = new TLorentzVector();
0021 TLorentzVector *Hadq = new TLorentzVector();
0022 TLorentzVector *Hadj = new TLorentzVector();
0023 TLorentzVector *Hadk = new TLorentzVector();
0024 TLorentzVector *Hadb = new TLorentzVector();
0025 TLorentzVector *Hadbbar = new TLorentzVector();
0026 Hadp->SetPxPyPzE(topJets[0].px(), topJets[0].py(), topJets[0].pz(), topJets[3].energy());
0027 Hadq->SetPxPyPzE(topJets[1].px(), topJets[1].py(), topJets[1].pz(), topJets[2].energy());
0028 Hadj->SetPxPyPzE(topJets[2].px(), topJets[2].py(), topJets[2].pz(), topJets[2].energy());
0029 Hadk->SetPxPyPzE(topJets[3].px(), topJets[3].py(), topJets[3].pz(), topJets[3].energy());
0030 Hadb->SetPxPyPzE(topJets[4].px(), topJets[4].py(), topJets[4].pz(), topJets[1].energy());
0031 Hadbbar->SetPxPyPzE(topJets[4].px(), topJets[4].py(), topJets[4].pz(), topJets[1].energy());
0032
0033
0034 std::sort(topJets.begin(), topJets.end(), EtComparator);
0035
0036
0037 double EtSum = topJets[5].et() + topJets[5].et();
0038 double Obs1 = (EtSum > 0 ? EtSum : -1);
0039 evtselectVarVal.push_back(std::pair<unsigned int, double>(1, Obs1));
0040
0041
0042
0043
0044
0045 std::sort(topJets.begin(), topJets.end(), BdiscComparator);
0046
0047 double BGap = topJets[1].bDiscriminator("trackCountingJetTags") - topJets[2].bDiscriminator("trackCountingJetTags");
0048 double Obs2 = (BGap > 0 ? log(BGap) : -1);
0049 evtselectVarVal.push_back(std::pair<unsigned int, double>(2, Obs2));
0050
0051
0052 double N = 0, D = 0, C_tmp = 0, C = 1000;
0053 double nx, ny, nz;
0054
0055
0056
0057 for (unsigned int i = 0; i < 6; i++) {
0058 D += fabs(topJets[i].pt());
0059 }
0060
0061 if ((D > 0)) {
0062
0063 for (unsigned int i = 0; i < 360; i++) {
0064 nx = cos((2 * PI / 360) * i);
0065 ny = sin((2 * PI / 360) * i);
0066 nz = 0;
0067 N = 0;
0068 for (unsigned int i = 0; i < 4; i++) {
0069 N += fabs(topJets[i].px() * nx + topJets[i].py() * ny + topJets[i].pz() * nz);
0070 }
0071 C_tmp = 2 * N / D;
0072 if (C_tmp < C)
0073 C = C_tmp;
0074 }
0075 }
0076
0077 double Obs3 = (C != 1000 ? C : -1);
0078 evtselectVarVal.push_back(std::pair<unsigned int, double>(3, Obs3));
0079
0080
0081 double HT = 0;
0082 for (unsigned int i = 0; i < 6; i++) {
0083 HT += topJets[i].et();
0084 }
0085
0086 double Obs4 = (HT != 0 ? HT : -1);
0087 evtselectVarVal.push_back(std::pair<unsigned int, double>(4, Obs4));
0088
0089
0090 math::XYZTLorentzVector pjets;
0091
0092 for (unsigned int i = 0; i < 6; i++) {
0093 pjets += topJets[i].p4();
0094 }
0095 double MT = sqrt(std::pow(pjets.mass(), 2));
0096 double Obs5 = (MT > 0 ? MT : -1);
0097 evtselectVarVal.push_back(std::pair<unsigned int, double>(5, Obs5));
0098
0099
0100
0101 std::sort(topJets.begin(), topJets.end(), EtComparator);
0102
0103 double px1 = topJets[2].px();
0104 double px2 = topJets[3].px();
0105 double py1 = topJets[2].py();
0106 double py2 = topJets[3].py();
0107 double pz1 = topJets[2].pz();
0108 double pz2 = topJets[3].pz();
0109 double E1 = topJets[2].energy();
0110 double E2 = topJets[3].energy();
0111 double px3 = topJets[4].px();
0112 double px4 = topJets[5].px();
0113 double py3 = topJets[4].py();
0114 double py4 = topJets[5].py();
0115 double pz3 = topJets[4].pz();
0116 double pz4 = topJets[5].pz();
0117 double E3 = topJets[4].energy();
0118 double E4 = topJets[5].energy();
0119
0120 TLorentzVector *LightJet1 = new TLorentzVector();
0121 TLorentzVector *LightJet2 = new TLorentzVector();
0122 TLorentzVector *LightJet3 = new TLorentzVector();
0123 TLorentzVector *LightJet4 = new TLorentzVector();
0124
0125 LightJet1->SetPxPyPzE(px1, py1, pz1, E1);
0126
0127 LightJet2->SetPxPyPzE(px2, py2, pz2, E2);
0128
0129 LightJet3->SetPxPyPzE(px3, py3, pz3, E3);
0130
0131 LightJet4->SetPxPyPzE(px4, py4, pz4, E4);
0132
0133
0134 double CosTheta1 = cos(LightJet2->Angle(LightJet1->Vect()));
0135 double CosTheta2 = cos(LightJet4->Angle(LightJet3->Vect()));
0136
0137 double Obs6 = (-1 < CosTheta1 ? CosTheta1 : -2);
0138 evtselectVarVal.push_back(std::pair<unsigned int, double>(6, Obs6));
0139 double Obs7 = (-1 < CosTheta2 ? CosTheta2 : -2);
0140 evtselectVarVal.push_back(std::pair<unsigned int, double>(7, Obs7));
0141
0142 delete LightJet1;
0143 delete LightJet2;
0144 delete LightJet3;
0145 delete LightJet4;
0146
0147
0148
0149 std::sort(topJets.begin(), topJets.end(), BdiscComparator);
0150
0151 double BjetsBdiscSum =
0152 topJets[0].bDiscriminator("trackCountingJetTags") + topJets[1].bDiscriminator("trackCountingJetTags");
0153 double Ljets1BdiscSum =
0154 topJets[2].bDiscriminator("trackCountingJetTags") + topJets[3].bDiscriminator("trackCountingJetTags");
0155 double Ljets2BdiscSum =
0156 topJets[4].bDiscriminator("trackCountingJetTags") + topJets[5].bDiscriminator("trackCountingJetTags");
0157
0158
0159
0160 double Obs8 = (Ljets1BdiscSum != 0 ? log(BjetsBdiscSum / Ljets1BdiscSum) : -1);
0161 evtselectVarVal.push_back(std::pair<unsigned int, double>(8, Obs8));
0162 double Obs9 = (Ljets2BdiscSum != 0 ? log(BjetsBdiscSum / Ljets2BdiscSum) : -1);
0163 evtselectVarVal.push_back(std::pair<unsigned int, double>(9, Obs9));
0164 double Obs10 = (BGap > 0 ? log(BjetsBdiscSum * BGap) : -1);
0165 evtselectVarVal.push_back(std::pair<unsigned int, double>(10, Obs10));
0166
0167
0168 double Obs11 = (HT != 0 ? (HT - EtSum) / HT : -1);
0169 evtselectVarVal.push_back(std::pair<unsigned int, double>(11, Obs11));
0170
0171
0172
0173 TMatrixDSym Matrix(3);
0174 double PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Hadj->Px(), 2) +
0175 std::pow(Hadk->Px(), 2) + std::pow(Hadbbar->Px(), 2);
0176 double PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Hadj->Py(), 2) +
0177 std::pow(Hadk->Py(), 2) + std::pow(Hadbbar->Py(), 2);
0178 double PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Hadj->Pz(), 2) +
0179 std::pow(Hadk->Pz(), 2) + std::pow(Hadbbar->Pz(), 2);
0180 double P2 = PX2 + PY2 + PZ2;
0181
0182 double PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Hadj->Px() * Hadj->Py() +
0183 Hadk->Px() * Hadk->Py() + Hadbbar->Px() * Hadbbar->Py();
0184 double PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Hadj->Px() * Hadj->Pz() +
0185 Hadk->Px() * Hadk->Pz() + Hadbbar->Px() * Hadbbar->Pz();
0186 double PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Hadj->Py() * Hadj->Pz() +
0187 Hadk->Py() * Hadk->Pz() + Hadbbar->Py() * Hadbbar->Pz();
0188
0189 Matrix(0, 0) = PX2 / P2;
0190 Matrix(0, 1) = PXY / P2;
0191 Matrix(0, 2) = PXZ / P2;
0192 Matrix(1, 0) = PXY / P2;
0193 Matrix(1, 1) = PY2 / P2;
0194 Matrix(1, 2) = PYZ / P2;
0195 Matrix(2, 0) = PXZ / P2;
0196 Matrix(2, 1) = PYZ / P2;
0197 Matrix(2, 2) = PZ2 / P2;
0198
0199 TMatrixDSymEigen pTensor(Matrix);
0200
0201 std::vector<double> EigValues;
0202 EigValues.clear();
0203 for (int i = 0; i < 3; i++) {
0204 EigValues.push_back(pTensor.GetEigenValues()[i]);
0205 }
0206
0207 std::sort(EigValues.begin(), EigValues.end(), dComparator);
0208
0209 double Sphericity = 1.5 * (EigValues[1] + EigValues[2]);
0210 double Aplanarity = 1.5 * EigValues[2];
0211
0212 double Obs12 = (edm::isNotFinite(Sphericity) ? -1 : Sphericity);
0213 evtselectVarVal.push_back(std::pair<unsigned int, double>(12, Obs12));
0214
0215 double Obs13 = (edm::isNotFinite(Aplanarity) ? -1 : Aplanarity);
0216 evtselectVarVal.push_back(std::pair<unsigned int, double>(13, Obs13));
0217
0218
0219 TLorentzVector *TtbarSystem = new TLorentzVector();
0220 TtbarSystem->SetPx(Hadp->Px() + Hadq->Px() + Hadb->Px() + Hadj->Px() + Hadk->Px() + Hadbbar->Px());
0221 TtbarSystem->SetPy(Hadp->Py() + Hadq->Py() + Hadb->Py() + Hadj->Py() + Hadk->Py() + Hadbbar->Py());
0222 TtbarSystem->SetPz(Hadp->Pz() + Hadq->Pz() + Hadb->Pz() + Hadj->Pz() + Hadk->Pz() + Hadbbar->Pz());
0223 TtbarSystem->SetE(Hadp->E() + Hadq->E() + Hadb->E() + Hadj->E() + Hadk->E() + Hadbbar->E());
0224
0225 TVector3 BoostBackToCM = -(TtbarSystem->BoostVector());
0226 Hadp->Boost(BoostBackToCM);
0227 Hadq->Boost(BoostBackToCM);
0228 Hadb->Boost(BoostBackToCM);
0229 Hadj->Boost(BoostBackToCM);
0230 Hadk->Boost(BoostBackToCM);
0231 Hadbbar->Boost(BoostBackToCM);
0232
0233 double BOOST_PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
0234 std::pow(Hadj->Px(), 2) + std::pow(Hadk->Px(), 2) + std::pow(Hadbbar->Px(), 2);
0235 double BOOST_PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
0236 std::pow(Hadj->Py(), 2) + std::pow(Hadk->Py(), 2) + std::pow(Hadbbar->Py(), 2);
0237 double BOOST_PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
0238 std::pow(Hadj->Pz(), 2) + std::pow(Hadk->Pz(), 2) + std::pow(Hadbbar->Pz(), 2);
0239
0240 double BOOST_P2 = BOOST_PX2 + BOOST_PY2 + BOOST_PZ2;
0241
0242 double BOOST_PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
0243 Hadj->Px() * Hadj->Py() + Hadk->Px() * Hadk->Py() + Hadbbar->Px() * Hadbbar->Py();
0244 double BOOST_PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
0245 Hadj->Px() * Hadj->Pz() + Hadk->Px() * Hadk->Pz() + Hadbbar->Px() * Hadbbar->Pz();
0246 double BOOST_PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
0247 Hadj->Py() * Hadj->Pz() + Hadk->Py() * Hadk->Pz() + Hadbbar->Py() * Hadbbar->Pz();
0248
0249 TMatrixDSym BOOST_Matrix(3);
0250
0251 BOOST_Matrix(0, 0) = BOOST_PX2 / BOOST_P2;
0252 BOOST_Matrix(0, 1) = BOOST_PXY / BOOST_P2;
0253 BOOST_Matrix(0, 2) = BOOST_PXZ / BOOST_P2;
0254 BOOST_Matrix(1, 0) = BOOST_PXY / BOOST_P2;
0255 BOOST_Matrix(1, 1) = BOOST_PY2 / BOOST_P2;
0256 BOOST_Matrix(1, 2) = BOOST_PYZ / BOOST_P2;
0257 BOOST_Matrix(2, 0) = BOOST_PXZ / BOOST_P2;
0258 BOOST_Matrix(2, 1) = BOOST_PYZ / BOOST_P2;
0259 BOOST_Matrix(2, 2) = BOOST_PZ2 / BOOST_P2;
0260
0261 TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
0262
0263 std::vector<double> BOOST_EigValues;
0264 BOOST_EigValues.clear();
0265 for (int i = 0; i < 3; i++) {
0266 BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
0267 }
0268
0269 std::sort(BOOST_EigValues.begin(), BOOST_EigValues.end(), dComparator);
0270
0271 double BOOST_Sphericity = 1.5 * (BOOST_EigValues[1] + BOOST_EigValues[2]);
0272 double BOOST_Aplanarity = 1.5 * BOOST_EigValues[2];
0273
0274 double Obs14 = (edm::isNotFinite(BOOST_Sphericity) ? -1 : BOOST_Sphericity);
0275 evtselectVarVal.push_back(std::pair<unsigned int, double>(14, Obs14));
0276
0277 double Obs15 = (edm::isNotFinite(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity);
0278 evtselectVarVal.push_back(std::pair<unsigned int, double>(15, Obs15));
0279
0280
0281 double H = 0;
0282 for (unsigned int i = 0; i < 6; i++) {
0283 H += topJets[i].energy();
0284 }
0285 double Obs16 = (H != 0 ? HT / H : -1);
0286 evtselectVarVal.push_back(std::pair<unsigned int, double>(16, Obs16));
0287
0288
0289 double Obs17 = (BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets1BdiscSum) : -1);
0290 evtselectVarVal.push_back(std::pair<unsigned int, double>(17, Obs17));
0291 double Obs18 = (BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets2BdiscSum) : -1);
0292 evtselectVarVal.push_back(std::pair<unsigned int, double>(18, Obs18));
0293
0294
0295 double Obs19 = (HT != 0 ? EtSum / HT : -1);
0296 evtselectVarVal.push_back(std::pair<unsigned int, double>(19, Obs19));
0297
0298
0299 TS.setLRSignalEvtObservables(evtselectVarVal);
0300 }