TopEventSelection/src/TtHadLRSignalSelObservables.cc

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   //Assume the highest Et jets are the b-jets, the others the jets from the W - but how to work out which two are from which W? FIXME!!! for now don't sort jets in Et before after this...
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   //sort the topJets in Et
0034   std::sort(topJets.begin(), topJets.end(), EtComparator);
0035
0036   //Et-Sum of the lightest jets
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   //Log of the difference in Bdisc between the 2nd and the 3rd jets (ordered in Bdisc) - does that still
0042   //make sense for the fully hadronic channel as well?FIXME!!!
0043
0044   //sort the jets in Bdiscriminant
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   //Circularity of the event
0052   double N = 0, D = 0, C_tmp = 0, C = 1000;
0053   double nx, ny, nz;
0054
0055   // C = 2min(E(pt.n)^2/E(pt)^2) = 2*N/D but it is theorically preferable to use C'=PI/2*min(E|pt.n|/E|pt|), sum over all jets+lepton+MET (cf PhysRevD 48 R3953(Nov 1993))
0056
0057   for (unsigned int i = 0; i < 6; i++) {
0058     D += fabs(topJets[i].pt());
0059   }
0060
0061   if ((D > 0)) {
0062     // Loop over all the unit vectors in the transverse plane in order to find the miminum :
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   //HT variable (Et-sum of the six jets)
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   //Transverse Mass of the system
0090   math::XYZTLorentzVector pjets;
0091   // for the six jets
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   //CosTheta(Hadp,Hadq) and CosTheta(Hadj,Hadq)
0100   //sort the lightJets in Et
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   //LightJet1->Boost(BoostBackToCM);
0127   LightJet2->SetPxPyPzE(px2, py2, pz2, E2);
0128   //LightJet2->Boost(BoostBackToCM);
0129   LightJet3->SetPxPyPzE(px3, py3, pz3, E3);
0130   //LightJet1->Boost(BoostBackToCM);
0131   LightJet4->SetPxPyPzE(px4, py4, pz4, E4);
0132   //LightJet2->Boost(BoostBackToCM);
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   // try to find out more powerful observables related to the b-disc
0148   //sort the jets in Bdiscriminant
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   //std::cout<<"BjetsBdiscSum = "<<BjetsBdiscSum<<std::endl;
0158   //std::cout<<"LjetsBdiscSum = "<<LjetsBdiscSum<<std::endl;
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   // Et-Ratio between the two highest in Et jets and four highest jets
0168   double Obs11 = (HT != 0 ? (HT - EtSum) / HT : -1);
0169   evtselectVarVal.push_back(std::pair<unsigned int, double>(11, Obs11));
0170
0171   //Sphericity and Aplanarity without boosting back the system to CM frame
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   //Sphericity and Aplanarity with boosting back the system to CM frame
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   // Centrality of the six jets
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   // distance from the origin in the (BjetsBdiscSum, Ljets1BdiscSum) and  (BjetsBdiscSum, Ljets2BdiscSum)
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   // Ratio of the Et Sum of the two lightest jets over HT=Sum of the Et of the six highest jets in Et
0295   double Obs19 = (HT != 0 ? EtSum / HT : -1);
0296   evtselectVarVal.push_back(std::pair<unsigned int, double>(19, Obs19));
0297
0298   // Put the vector in the TtHadEvtSolution
0299   TS.setLRSignalEvtObservables(evtselectVarVal);
0300 }