15 std::vector<pat::Jet> topJets;
17 topJets.push_back(TS.
getHadp());
18 topJets.push_back(TS.
getHadq());
19 topJets.push_back(TS.
getHadj());
20 topJets.push_back(TS.
getHadk());
21 topJets.push_back(TS.
getHadb());
25 TLorentzVector *Hadp =
new TLorentzVector();
26 TLorentzVector *Hadq =
new TLorentzVector();
27 TLorentzVector *Hadj =
new TLorentzVector();
28 TLorentzVector *Hadk =
new TLorentzVector();
29 TLorentzVector *Hadb =
new TLorentzVector();
30 TLorentzVector *Hadbbar =
new TLorentzVector();
31 Hadp->SetPxPyPzE(topJets[0].px(),topJets[0].py(),topJets[0].pz(),topJets[3].
energy());
32 Hadq->SetPxPyPzE(topJets[1].px(),topJets[1].py(),topJets[1].pz(),topJets[2].
energy());
33 Hadj->SetPxPyPzE(topJets[2].px(),topJets[2].py(),topJets[2].pz(),topJets[2].
energy());
34 Hadk->SetPxPyPzE(topJets[3].px(),topJets[3].py(),topJets[3].pz(),topJets[3].
energy());
35 Hadb->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].
energy());
36 Hadbbar->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].
energy());
42 double EtSum = topJets[5].et()+topJets[5].et();
43 double Obs1 = (EtSum>0 ? EtSum : -1);
52 double BGap = topJets[1].bDiscriminator(
"trackCountingJetTags") - topJets[2].bDiscriminator(
"trackCountingJetTags");
53 double Obs2 = (BGap>0 ?
log(BGap) : -1);
57 double N=0,
D=0,C_tmp=0,
C=1000;
62 for(
unsigned int i=0;
i<6;
i++){
63 D += fabs(topJets[
i].
pt());
68 for(
unsigned int i=0;
i<360;
i++){
73 for(
unsigned int i=0;
i<4;
i++){
74 N += fabs(topJets[
i].px()*nx+topJets[
i].py()*ny+topJets[
i].pz()*nz);
77 if(C_tmp<
C)
C = C_tmp;
81 double Obs3 = (
C!=1000 ?
C : -1);
87 for(
unsigned int i=0;
i<6;
i++){
88 HT += topJets[
i].et();
91 double Obs4 = ( HT!=0 ? HT : -1);
97 for(
unsigned int i=0;
i<6;
i++){
98 pjets += topJets[
i].p4();
101 double Obs5 = ( MT>0 ? MT : -1);
108 double px1 = topJets[2].px();
double px2 = topJets[3].px();
109 double py1 = topJets[2].py();
double py2 = topJets[3].py();
110 double pz1 = topJets[2].pz();
double pz2 = topJets[3].pz();
111 double E1 = topJets[2].energy();
double E2 = topJets[3].energy();
112 double px3 = topJets[4].px();
double px4 = topJets[5].px();
113 double py3 = topJets[4].py();
double py4 = topJets[5].py();
114 double pz3 = topJets[4].pz();
double pz4 = topJets[5].pz();
115 double E3 = topJets[4].energy();
double E4 = topJets[5].energy();
117 TLorentzVector *LightJet1 =
new TLorentzVector();
118 TLorentzVector *LightJet2 =
new TLorentzVector();
119 TLorentzVector *LightJet3 =
new TLorentzVector();
120 TLorentzVector *LightJet4 =
new TLorentzVector();
122 LightJet1->SetPxPyPzE(px1,py1,pz1,E1);
124 LightJet2->SetPxPyPzE(px2,py2,pz2,E2);
126 LightJet3->SetPxPyPzE(px3,py3,pz3,E3);
128 LightJet4->SetPxPyPzE(px4,py4,pz4,E4);
131 double CosTheta1 =
cos(LightJet2->Angle(LightJet1->Vect()));
132 double CosTheta2 =
cos(LightJet4->Angle(LightJet3->Vect()));
134 double Obs6 = ( -1<CosTheta1 ? CosTheta1 : -2);
136 double Obs7 = ( -1<CosTheta2 ? CosTheta2 : -2);
148 double BjetsBdiscSum = topJets[0].bDiscriminator(
"trackCountingJetTags") + topJets[1].bDiscriminator(
"trackCountingJetTags");
149 double Ljets1BdiscSum = topJets[2].bDiscriminator(
"trackCountingJetTags") + topJets[3].bDiscriminator(
"trackCountingJetTags");
150 double Ljets2BdiscSum = topJets[4].bDiscriminator(
"trackCountingJetTags") + topJets[5].bDiscriminator(
"trackCountingJetTags");
154 double Obs8 = (Ljets1BdiscSum !=0 ?
log(BjetsBdiscSum/Ljets1BdiscSum) : -1);
156 double Obs9 = (Ljets2BdiscSum !=0 ?
log(BjetsBdiscSum/Ljets2BdiscSum) : -1);
158 double Obs10 = (BGap>0 ?
log(BjetsBdiscSum*BGap) : -1);
163 double Obs11 = (HT!=0 ? (HT-EtSum)/HT : -1);
176 double P2 = PX2+PY2+PZ2;
178 double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
179 Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
180 double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
181 Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
182 double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
183 Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
189 TMatrixDSymEigen pTensor(
Matrix);
191 std::vector<double> EigValues;
193 for(
int i=0;
i<3;
i++){
194 EigValues.push_back(pTensor.GetEigenValues()[
i]);
197 std::sort(EigValues.begin(),EigValues.end(),
dComparator);
199 double Sphericity = 1.5*(EigValues[1]+EigValues[2]);
200 double Aplanarity = 1.5*EigValues[2];
210 TLorentzVector *TtbarSystem =
new TLorentzVector();
211 TtbarSystem->SetPx(Hadp->Px()+Hadq->Px()+Hadb->Px()+Hadj->Px()+Hadk->Px()+Hadbbar->Px());
212 TtbarSystem->SetPy(Hadp->Py()+Hadq->Py()+Hadb->Py()+Hadj->Py()+Hadk->Py()+Hadbbar->Py());
213 TtbarSystem->SetPz(Hadp->Pz()+Hadq->Pz()+Hadb->Pz()+Hadj->Pz()+Hadk->Pz()+Hadbbar->Pz());
214 TtbarSystem->SetE(Hadp->E()+Hadq->E()+Hadb->E()+Hadj->E()+Hadk->E()+Hadbbar->E());
216 TVector3 BoostBackToCM = -(TtbarSystem->BoostVector());
217 Hadp->Boost(BoostBackToCM);
218 Hadq->Boost(BoostBackToCM);
219 Hadb->Boost(BoostBackToCM);
220 Hadj->Boost(BoostBackToCM);
221 Hadk->Boost(BoostBackToCM);
222 Hadbbar->Boost(BoostBackToCM);
231 double BOOST_P2 = BOOST_PX2+BOOST_PY2+BOOST_PZ2;
233 double BOOST_PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
234 Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
235 double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
236 Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
237 double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
238 Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
240 TMatrixDSym BOOST_Matrix(3);
242 BOOST_Matrix(0,0) = BOOST_PX2/BOOST_P2; BOOST_Matrix(0,1) = BOOST_PXY/BOOST_P2; BOOST_Matrix(0,2) = BOOST_PXZ/BOOST_P2;
243 BOOST_Matrix(1,0) = BOOST_PXY/BOOST_P2; BOOST_Matrix(1,1) = BOOST_PY2/BOOST_P2; BOOST_Matrix(1,2) = BOOST_PYZ/BOOST_P2;
244 BOOST_Matrix(2,0) = BOOST_PXZ/BOOST_P2; BOOST_Matrix(2,1) = BOOST_PYZ/BOOST_P2; BOOST_Matrix(2,2) = BOOST_PZ2/BOOST_P2;
246 TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
248 std::vector<double> BOOST_EigValues;
249 BOOST_EigValues.clear();
250 for(
int i=0;
i<3;
i++){
251 BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[
i]);
254 std::sort(BOOST_EigValues.begin(),BOOST_EigValues.end(),
dComparator);
256 double BOOST_Sphericity = 1.5*(BOOST_EigValues[1]+BOOST_EigValues[2]);
257 double BOOST_Aplanarity = 1.5*BOOST_EigValues[2];
259 double Obs14 = (
edm::isNotFinite(BOOST_Sphericity) ? -1 : BOOST_Sphericity );
262 double Obs15 = (
edm::isNotFinite(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity );
268 for(
unsigned int i=0;
i<6;
i++){
269 H += topJets[
i].energy();
271 double Obs16 = ( H != 0 ? HT/H : -1 );
276 double Obs17 = ( BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets1BdiscSum) : -1 );
278 double Obs18 = ( BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets2BdiscSum) : -1 );
282 double Obs19 = ( HT != 0 ? EtSum/HT : -1 );
constexpr bool isNotFinite(T x)
Sin< T >::type sin(const T &t)
pat::Jet getHadbbar() const
void setLRSignalEvtObservables(const std::vector< std::pair< unsigned int, double > > &varval)
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
CompareDouble dComparator
std::vector< std::pair< unsigned int, double > > evtselectVarVal
Power< A, B >::type pow(const A &a, const B &b)