Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:16 

0001 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelObservables.h"
0002 #include "FWCore/Utilities/interface/isFinite.h"
0003 
0004 TtSemiLRSignalSelObservables::TtSemiLRSignalSelObservables() {}
0005 
0006 TtSemiLRSignalSelObservables::~TtSemiLRSignalSelObservables() {}
0007 
0008 void TtSemiLRSignalSelObservables::operator()(TtSemiEvtSolution &TS, const std::vector<pat::Jet> &SelectedTopJets) {
0009   evtselectVarVal.clear();
0010 
0011   // choice the B-tag algorithm
0012   const char *BtagAlgo = "trackCountingJetTags";
0013 
0014   // activate the "debug mode"
0015   bool DEBUG = false;
0016 
0017   if (DEBUG)
0018     std::cout << "---------- Start calculating the LR observables ----------" << std::endl;
0019 
0020   std::vector<pat::Jet> TopJets;
0021   TopJets.clear();
0022 
0023   for (unsigned int i = 0; i < SelectedTopJets.size(); i++) {
0024     TopJets.push_back(SelectedTopJets[i]);
0025   }
0026 
0027   //sort the TopJets in Et
0028   std::sort(TopJets.begin(), TopJets.end(), EtComparator);
0029 
0030   TLorentzVector *Hadp = new TLorentzVector();
0031   Hadp->SetPxPyPzE(TopJets[3].px(), TopJets[3].py(), TopJets[3].pz(), TopJets[3].energy());
0032 
0033   TLorentzVector *Hadq = new TLorentzVector();
0034   Hadq->SetPxPyPzE(TopJets[2].px(), TopJets[2].py(), TopJets[2].pz(), TopJets[2].energy());
0035 
0036   TLorentzVector *Hadb = new TLorentzVector();
0037   Hadb->SetPxPyPzE(TopJets[1].px(), TopJets[1].py(), TopJets[1].pz(), TopJets[1].energy());
0038 
0039   TLorentzVector *Lepb = new TLorentzVector();
0040   Lepb->SetPxPyPzE(TopJets[0].px(), TopJets[0].py(), TopJets[0].pz(), TopJets[0].energy());
0041 
0042   TLorentzVector *Lept = new TLorentzVector();
0043   if (TS.getDecay() == "muon")
0044     Lept->SetPxPyPzE(TS.getMuon().px(), TS.getMuon().py(), TS.getMuon().pz(), TS.getMuon().energy());
0045   else if (TS.getDecay() == "electron")
0046     Lept->SetPxPyPzE(TS.getElectron().px(), TS.getElectron().py(), TS.getElectron().pz(), TS.getElectron().energy());
0047 
0048   TLorentzVector *Lepn = new TLorentzVector();
0049   Lepn->SetPxPyPzE(TS.getNeutrino().px(), TS.getNeutrino().py(), TS.getNeutrino().pz(), TS.getNeutrino().energy());
0050 
0051   // Calculation of the pz of the neutrino due to W-mass constraint
0052 
0053   MEzCalculator *Mez = new MEzCalculator();
0054   Mez->SetMET(TS.getNeutrino());
0055   if (TS.getDecay() == "muon")
0056     Mez->SetLepton(TS.getMuon());
0057   else
0058     Mez->SetLepton(TS.getElectron(), false);
0059   double MEZ = Mez->Calculate();
0060   Lepn->SetPz(MEZ);
0061 
0062   // Pt of the lepton
0063 
0064   double Obs1 = Lept->Pt();
0065   evtselectVarVal.push_back(std::pair<unsigned int, double>(1, Obs1));
0066   if (DEBUG)
0067     std::cout << "------ LR observable 1 " << Obs1 << " calculated ------" << std::endl;
0068 
0069   // Missing transverse energy
0070 
0071   double Obs2 = TS.getNeutrino().et();
0072   evtselectVarVal.push_back(std::pair<unsigned int, double>(2, Obs2));
0073   if (DEBUG)
0074     std::cout << "------ LR observable 2 " << Obs2 << " calculated ------" << std::endl;
0075 
0076   //HT variable (Et-sum of the four jets)
0077 
0078   double HT = 0;
0079   for (unsigned int i = 0; i < 4; i++) {
0080     HT += TopJets[i].et();
0081   }
0082 
0083   double Obs3 = HT;
0084   evtselectVarVal.push_back(std::pair<unsigned int, double>(3, Obs3));
0085   if (DEBUG)
0086     std::cout << "------ LR observable 3 " << Obs3 << " calculated ------" << std::endl;
0087 
0088   //Et-Sum of the lightest jets
0089 
0090   double EtSum = TopJets[2].et() + TopJets[3].et();
0091   double Obs4 = EtSum;
0092   evtselectVarVal.push_back(std::pair<unsigned int, double>(4, Obs4));
0093   if (DEBUG)
0094     std::cout << "------ LR observable 4 " << Obs4 << " calculated ------" << std::endl;
0095 
0096   // Et-Ratio between the two lowest jets in Et and four highest jets
0097 
0098   double Obs5 = EtSum / HT;
0099   evtselectVarVal.push_back(std::pair<unsigned int, double>(5, Obs5));
0100   if (DEBUG)
0101     std::cout << "------ LR observable 5 " << Obs5 << " calculated ------" << std::endl;
0102 
0103   // Et-Ratio between the two highest jets in Et and four highest jets
0104 
0105   double Obs6 = (HT - EtSum) / HT;
0106   evtselectVarVal.push_back(std::pair<unsigned int, double>(6, Obs6));
0107   if (DEBUG)
0108     std::cout << "------ LR observable 6 " << Obs6 << " calculated ------" << std::endl;
0109 
0110   // Transverse Mass of the system
0111 
0112   TLorentzVector TtbarSystem = (*Hadp) + (*Hadq) + (*Hadb) + (*Lepb) + (*Lept) + (*Lepn);
0113   double MT = TtbarSystem.Mt();
0114   double Obs7 = MT;
0115   evtselectVarVal.push_back(std::pair<unsigned int, double>(7, Obs7));
0116   if (DEBUG)
0117     std::cout << "------ LR observable 7 " << Obs7 << " calculated ------" << std::endl;
0118 
0119   // Observables related to the b-disc (jets ordered in Bdisc)
0120 
0121   //sort the TopJets in Bdiscriminant
0122   std::sort(TopJets.begin(), TopJets.end(), BdiscComparator);
0123 
0124   double Obs8;
0125   double Obs9;
0126   double Obs10;
0127   double Obs11;
0128 
0129   // Difference in bdisc between the 2nd and the 3rd jets
0130   double BGap = TopJets[1].bDiscriminator(BtagAlgo) - TopJets[2].bDiscriminator(BtagAlgo);
0131 
0132   // Sum of the bdisc of the two highest/lowest jets
0133   double BjetsBdiscSum = TopJets[0].bDiscriminator(BtagAlgo) + TopJets[1].bDiscriminator(BtagAlgo);
0134   double LjetsBdiscSum = TopJets[2].bDiscriminator(BtagAlgo) + TopJets[3].bDiscriminator(BtagAlgo);
0135 
0136   Obs8 = (TopJets[2].bDiscriminator(BtagAlgo) > -10 ? log(BGap) : -5);
0137   if (DEBUG)
0138     std::cout << "------ LR observable 8 " << Obs8 << " calculated ------" << std::endl;
0139   Obs9 = (BjetsBdiscSum * BGap);
0140   if (DEBUG)
0141     std::cout << "------ LR observable 9 " << Obs9 << " calculated ------" << std::endl;
0142   Obs10 = (BjetsBdiscSum / LjetsBdiscSum);
0143   if (DEBUG)
0144     std::cout << "------ LR observable 10 " << Obs10 << " calculated ------" << std::endl;
0145   // Distance from the origin in the (BjetsBdiscSum, LjetsBdiscSum) plane
0146   Obs11 = 0.707 * ((BjetsBdiscSum + LjetsBdiscSum) / 2 + 10);
0147   if (DEBUG)
0148     std::cout << "------ LR observable 11 " << Obs11 << " calculated ------" << std::endl;
0149 
0150   evtselectVarVal.push_back(std::pair<unsigned int, double>(8, Obs8));
0151   evtselectVarVal.push_back(std::pair<unsigned int, double>(9, Obs9));
0152   evtselectVarVal.push_back(std::pair<unsigned int, double>(10, Obs10));
0153   evtselectVarVal.push_back(std::pair<unsigned int, double>(11, Obs11));
0154 
0155   //sort the TopJets in Et
0156   std::sort(TopJets.begin(), TopJets.end(), EtComparator);
0157 
0158   // Circularity of the event
0159 
0160   double N = 0, D = 0, C_tmp = 0, C = 10000;
0161   double N_NoNu = 0, D_NoNu = 0, C_NoNu_tmp = 0, C_NoNu = 10000;
0162   double nx, ny, nz;
0163 
0164   // 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))
0165 
0166   for (unsigned int i = 0; i < 4; i++) {
0167     D += fabs(TopJets[i].pt());
0168   }
0169   // modified the 26th of September to calculate also the circularity without the neutrino contribution
0170   D += fabs(Lept->Pt());
0171   D_NoNu = D;
0172   D += fabs(Lepn->Pt());
0173 
0174   if ((D > 0)) {
0175     // Loop over all the unit vectors in the transverse plane in order to find the miminum :
0176     for (unsigned int i = 0; i < 720; i++) {
0177       nx = cos((2 * PI / 720) * i);
0178       ny = sin((2 * PI / 720) * i);
0179       nz = 0;
0180       N = 0;
0181 
0182       for (unsigned int i = 0; i < 4; i++) {
0183         N += fabs(TopJets[i].px() * nx + TopJets[i].py() * ny + TopJets[i].pz() * nz);
0184       }
0185       // modified the 26th of September to calculate also the circularity without the neutrino contribution
0186       N += fabs(Lept->Px() * nx + Lept->Py() * ny + Lept->Pz() * nz);
0187       N_NoNu = N;
0188       N += fabs(Lepn->Px() * nx + Lepn->Py() * ny + Lepn->Pz() * nz);
0189 
0190       C_tmp = 2 * N / D;
0191 
0192       // modified the 26th of September to calculate also the circularity without the neutrino contribution
0193       C_NoNu_tmp = 2 * N_NoNu / D_NoNu;
0194 
0195       if (C_tmp < C)
0196         C = C_tmp;
0197 
0198       // modified the 26th of September to calculate also the circularity without the neutrino contribution
0199       if (C_NoNu_tmp < C_NoNu)
0200         C_NoNu = C_NoNu_tmp;
0201     }
0202   }
0203 
0204   double Obs12 = (C != 10000 ? C : 0);
0205   evtselectVarVal.push_back(std::pair<unsigned int, double>(12, Obs12));
0206   if (DEBUG)
0207     std::cout << "------ LR observable 12 " << Obs12 << " calculated ------" << std::endl;
0208 
0209   // Circularity of the event without neutrino
0210 
0211   double Obs13 = (C_NoNu != 10000 ? C_NoNu : 0);
0212   evtselectVarVal.push_back(std::pair<unsigned int, double>(13, Obs13));
0213   if (DEBUG)
0214     std::cout << "------ LR observable 13 " << Obs13 << " calculated ------" << std::endl;
0215 
0216   // Centrality of the four highest jets
0217 
0218   double H = 0;
0219   for (unsigned int i = 0; i < 4; i++) {
0220     H += TopJets[i].energy();
0221   }
0222 
0223   double Obs14 = HT / H;
0224   evtselectVarVal.push_back(std::pair<unsigned int, double>(14, Obs14));
0225   if (DEBUG)
0226     std::cout << "------ LR observable 14 " << Obs14 << " calculated ------" << std::endl;
0227 
0228   // Sphericity and Aplanarity without boosting back the system to CM frame
0229 
0230   TMatrixDSym Matrix(3);
0231 
0232   double PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Lepb->Px(), 2) +
0233                std::pow(Lept->Px(), 2) + std::pow(Lepn->Px(), 2);
0234   double PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Lepb->Py(), 2) +
0235                std::pow(Lept->Py(), 2) + std::pow(Lepn->Py(), 2);
0236   double PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Lepb->Pz(), 2) +
0237                std::pow(Lept->Pz(), 2) + std::pow(Lepn->Pz(), 2);
0238 
0239   double P2 = PX2 + PY2 + PZ2;
0240 
0241   double PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Lepb->Px() * Lepb->Py() +
0242                Lept->Px() * Lept->Py() + Lepn->Px() * Lepn->Py();
0243   double PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Lepb->Px() * Lepb->Pz() +
0244                Lept->Px() * Lept->Pz() + Lepn->Px() * Lepn->Pz();
0245   double PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Lepb->Py() * Lepb->Pz() +
0246                Lept->Py() * Lept->Pz() + Lepn->Py() * Lepn->Pz();
0247 
0248   Matrix(0, 0) = PX2 / P2;
0249   Matrix(0, 1) = PXY / P2;
0250   Matrix(0, 2) = PXZ / P2;
0251   Matrix(1, 0) = PXY / P2;
0252   Matrix(1, 1) = PY2 / P2;
0253   Matrix(1, 2) = PYZ / P2;
0254   Matrix(2, 0) = PXZ / P2;
0255   Matrix(2, 1) = PYZ / P2;
0256   Matrix(2, 2) = PZ2 / P2;
0257 
0258   TMatrixDSymEigen pTensor(Matrix);
0259 
0260   std::vector<double> EigValues;
0261   EigValues.clear();
0262   for (int i = 0; i < 3; i++) {
0263     EigValues.push_back(pTensor.GetEigenValues()[i]);
0264   }
0265 
0266   std::sort(EigValues.begin(), EigValues.end(), dComparator);
0267 
0268   double Sphericity = 1.5 * (EigValues[1] + EigValues[2]);
0269   double Aplanarity = 1.5 * EigValues[2];
0270 
0271   double Obs15 = (edm::isNotFinite(Sphericity) ? 0 : Sphericity);
0272   evtselectVarVal.push_back(std::pair<unsigned int, double>(15, Obs15));
0273   if (DEBUG)
0274     std::cout << "------ LR observable 15 " << Obs15 << " calculated ------" << std::endl;
0275 
0276   double Obs16 = (edm::isNotFinite(Aplanarity) ? 0 : Aplanarity);
0277   evtselectVarVal.push_back(std::pair<unsigned int, double>(16, Obs16));
0278   if (DEBUG)
0279     std::cout << "------ LR observable 16 " << Obs16 << " calculated ------" << std::endl;
0280 
0281   // Sphericity and Aplanarity with boosting back the system to CM frame
0282 
0283   TVector3 BoostBackToCM = -(TtbarSystem.BoostVector());
0284   Hadp->Boost(BoostBackToCM);
0285   Hadq->Boost(BoostBackToCM);
0286   Hadb->Boost(BoostBackToCM);
0287   Lepb->Boost(BoostBackToCM);
0288   Lept->Boost(BoostBackToCM);
0289   Lepn->Boost(BoostBackToCM);
0290 
0291   double BOOST_PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
0292                      std::pow(Lepb->Px(), 2) + std::pow(Lept->Px(), 2) + std::pow(Lepn->Px(), 2);
0293   double BOOST_PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
0294                      std::pow(Lepb->Py(), 2) + std::pow(Lept->Py(), 2) + std::pow(Lepn->Py(), 2);
0295   double BOOST_PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
0296                      std::pow(Lepb->Pz(), 2) + std::pow(Lept->Pz(), 2) + std::pow(Lepn->Pz(), 2);
0297 
0298   double BOOST_P2 = BOOST_PX2 + BOOST_PY2 + BOOST_PZ2;
0299 
0300   double BOOST_PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
0301                      Lepb->Px() * Lepb->Py() + Lept->Px() * Lept->Py() + Lepn->Px() * Lepn->Py();
0302   double BOOST_PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
0303                      Lepb->Px() * Lepb->Pz() + Lept->Px() * Lept->Pz() + Lepn->Px() * Lepn->Pz();
0304   double BOOST_PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
0305                      Lepb->Py() * Lepb->Pz() + Lept->Py() * Lept->Pz() + Lepn->Py() * Lepn->Pz();
0306 
0307   TMatrixDSym BOOST_Matrix(3);
0308 
0309   BOOST_Matrix(0, 0) = BOOST_PX2 / BOOST_P2;
0310   BOOST_Matrix(0, 1) = BOOST_PXY / BOOST_P2;
0311   BOOST_Matrix(0, 2) = BOOST_PXZ / BOOST_P2;
0312   BOOST_Matrix(1, 0) = BOOST_PXY / BOOST_P2;
0313   BOOST_Matrix(1, 1) = BOOST_PY2 / BOOST_P2;
0314   BOOST_Matrix(1, 2) = BOOST_PYZ / BOOST_P2;
0315   BOOST_Matrix(2, 0) = BOOST_PXZ / BOOST_P2;
0316   BOOST_Matrix(2, 1) = BOOST_PYZ / BOOST_P2;
0317   BOOST_Matrix(2, 2) = BOOST_PZ2 / BOOST_P2;
0318 
0319   TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
0320 
0321   std::vector<double> BOOST_EigValues;
0322   BOOST_EigValues.clear();
0323   for (int i = 0; i < 3; i++) {
0324     BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
0325   }
0326 
0327   std::sort(BOOST_EigValues.begin(), BOOST_EigValues.end(), dComparator);
0328 
0329   double BOOST_Sphericity = 1.5 * (BOOST_EigValues[1] + BOOST_EigValues[2]);
0330   double BOOST_Aplanarity = 1.5 * BOOST_EigValues[2];
0331 
0332   double Obs17 = (edm::isNotFinite(BOOST_Sphericity) ? 0 : BOOST_Sphericity);
0333   evtselectVarVal.push_back(std::pair<unsigned int, double>(17, Obs17));
0334   if (DEBUG)
0335     std::cout << "------ LR observable 17 " << Obs17 << " calculated ------" << std::endl;
0336 
0337   double Obs18 = (edm::isNotFinite(BOOST_Aplanarity) ? 0 : BOOST_Aplanarity);
0338   evtselectVarVal.push_back(std::pair<unsigned int, double>(18, Obs18));
0339   if (DEBUG)
0340     std::cout << "------ LR observable 18 " << Obs18 << " calculated ------" << std::endl;
0341 
0342   //ratio between ET of the fifth jet and HT
0343 
0344   double Obs19 = (TopJets.size() > 4) ? TopJets[4].et() / HT : 1.0;
0345   evtselectVarVal.push_back(std::pair<unsigned int, double>(19, Obs19));
0346   if (DEBUG)
0347     std::cout << "------ LR observable 19 " << Obs19 << " calculated ------" << std::endl;
0348 
0349   // HT variable calculated with all the jets in the event.
0350 
0351   double HT_alljets = 0;
0352   double H_alljets = 0;
0353   for (unsigned int i = 0; i < TopJets.size(); i++) {
0354     HT_alljets += TopJets[i].et();
0355     H_alljets += TopJets[i].energy();
0356   }
0357   double Obs20 = HT_alljets;
0358   evtselectVarVal.push_back(std::pair<unsigned int, double>(20, Obs20));
0359   if (DEBUG)
0360     std::cout << "------ LR observable 20 " << Obs20 << " calculated ------" << std::endl;
0361 
0362   // HT3 = HT calculated with all jets except the two leading jets
0363 
0364   double HT3_alljets = 0;
0365   for (unsigned int i = 2; i < TopJets.size(); i++) {
0366     HT3_alljets += TopJets[i].et();
0367   }
0368   double Obs21 = HT3_alljets;
0369   evtselectVarVal.push_back(std::pair<unsigned int, double>(21, Obs21));
0370   if (DEBUG)
0371     std::cout << "------ LR observable 21 " << Obs21 << " calculated ------" << std::endl;
0372 
0373   // ET0, ratio of the Et of the leading and HT_alljets
0374 
0375   double ET0 = TopJets[0].et() / HT_alljets;
0376   double Obs22 = ET0;
0377   evtselectVarVal.push_back(std::pair<unsigned int, double>(22, Obs22));
0378   if (DEBUG)
0379     std::cout << "------ LR observable 22 " << Obs22 << " calculated ------" << std::endl;
0380 
0381   // Centrality of the event computed with all jets
0382 
0383   double Obs23 = ((H_alljets > 0) ? HT_alljets / H_alljets : 0);
0384   evtselectVarVal.push_back(std::pair<unsigned int, double>(23, Obs23));
0385   if (DEBUG)
0386     std::cout << "------ LR observable 23 " << Obs23 << " calculated ------" << std::endl;
0387 
0388   //Fox-Wolfram momenta (1st to 6th), modified for hadron collider and using a Legendre polynomials expansion
0389 
0390   double FW_momentum_0 = 0, FW_momentum_1 = 0, FW_momentum_2 = 0, FW_momentum_3 = 0, FW_momentum_4 = 0,
0391          FW_momentum_5 = 0, FW_momentum_6 = 0;
0392 
0393   for (unsigned int i = 0; i < TopJets.size(); i++) {
0394     for (unsigned int j = 0; j < TopJets.size(); j++) {
0395       double ET_ij_over_ETSum2 = TopJets[i].et() * TopJets[j].et() / (std::pow(HT_alljets, 2));
0396       double cosTheta_ij =
0397           (TopJets[i].px() * TopJets[j].px() + TopJets[i].py() * TopJets[j].py() + TopJets[i].pz() * TopJets[j].pz()) /
0398           (TopJets[i].p4().R() * TopJets[j].p4().R());
0399       FW_momentum_0 += ET_ij_over_ETSum2;
0400       FW_momentum_1 += ET_ij_over_ETSum2 * cosTheta_ij;
0401       FW_momentum_2 += ET_ij_over_ETSum2 * 0.5 * (3 * std::pow(cosTheta_ij, 2) - 1);
0402       FW_momentum_3 += ET_ij_over_ETSum2 * 0.5 * (5 * std::pow(cosTheta_ij, 3) - 3 * cosTheta_ij);
0403       FW_momentum_4 += ET_ij_over_ETSum2 * 0.125 * (35 * std::pow(cosTheta_ij, 4) - 30 * std::pow(cosTheta_ij, 2) + 3);
0404       FW_momentum_5 += ET_ij_over_ETSum2 * 0.125 *
0405                        (63 * std::pow(cosTheta_ij, 5) - 70 * std::pow(cosTheta_ij, 3) + 15 * cosTheta_ij);
0406       FW_momentum_6 +=
0407           ET_ij_over_ETSum2 * 0.0625 *
0408           (231 * std::pow(cosTheta_ij, 6) - 315 * std::pow(cosTheta_ij, 4) + 105 * std::pow(cosTheta_ij, 2) - 5);
0409     }
0410   }
0411 
0412   double Obs24 = FW_momentum_0;
0413   evtselectVarVal.push_back(std::pair<unsigned int, double>(24, Obs24));
0414   if (DEBUG)
0415     std::cout << "------ LR observable 24 " << Obs24 << " calculated ------" << std::endl;
0416   double Obs25 = FW_momentum_1;
0417   evtselectVarVal.push_back(std::pair<unsigned int, double>(25, Obs25));
0418   if (DEBUG)
0419     std::cout << "------ LR observable 25 " << Obs25 << " calculated ------" << std::endl;
0420   double Obs26 = FW_momentum_2;
0421   evtselectVarVal.push_back(std::pair<unsigned int, double>(26, Obs26));
0422   if (DEBUG)
0423     std::cout << "------ LR observable 26 " << Obs26 << " calculated ------" << std::endl;
0424   double Obs27 = FW_momentum_3;
0425   evtselectVarVal.push_back(std::pair<unsigned int, double>(27, Obs27));
0426   if (DEBUG)
0427     std::cout << "------ LR observable 27 " << Obs27 << " calculated ------" << std::endl;
0428   double Obs28 = FW_momentum_4;
0429   evtselectVarVal.push_back(std::pair<unsigned int, double>(28, Obs28));
0430   if (DEBUG)
0431     std::cout << "------ LR observable 28 " << Obs28 << " calculated ------" << std::endl;
0432   double Obs29 = FW_momentum_5;
0433   evtselectVarVal.push_back(std::pair<unsigned int, double>(29, Obs29));
0434   if (DEBUG)
0435     std::cout << "------ LR observable 29 " << Obs29 << " calculated ------" << std::endl;
0436   double Obs30 = FW_momentum_6;
0437   evtselectVarVal.push_back(std::pair<unsigned int, double>(30, Obs30));
0438   if (DEBUG)
0439     std::cout << "------ LR observable 30 " << Obs30 << " calculated ------" << std::endl;
0440 
0441   // Thrust, thrust major and thrust minor
0442 
0443   TVector3 n(0, 0, 0), n_tmp(0, 0, 0);
0444 
0445   double Thrust_D = 0, Thrust_N = 0;
0446   double Thrust = -1, Thrust_tmp = 0;
0447   double Thrust_Major = -1, Thrust_Major_tmp = 0;
0448   double Thrust_Minor = -1;
0449 
0450   for (unsigned int i = 0; i < TopJets.size(); i++) {
0451     Thrust_D += TopJets[i].p();
0452   }
0453 
0454   Thrust_D += Lept->P();
0455 
0456   if (Thrust_D > 0) {
0457     // Calculation of the thrust :
0458     for (unsigned int i = 0; i < 720; i++) {
0459       for (unsigned int j = 0; j < 360; j++) {
0460         n_tmp.SetX(cos((2 * PI / 720) * i) * sin((PI / 360) * j));
0461         n_tmp.SetY(sin((2 * PI / 720) * i) * sin((PI / 360) * j));
0462         n_tmp.SetZ(cos((PI / 360) * j));
0463 
0464         for (unsigned int k = 0; k < TopJets.size(); k++) {
0465           Thrust_N +=
0466               fabs(TopJets[k].px() * (n_tmp.x()) + TopJets[k].py() * (n_tmp.y()) + TopJets[k].pz() * (n_tmp.z()));
0467         }
0468         Thrust_N += fabs(Lept->Px() * (n_tmp.x()) + Lept->Py() * (n_tmp.y()) + Lept->Pz() * (n_tmp.z()));
0469 
0470         Thrust_tmp = Thrust_N / Thrust_D;
0471 
0472         Thrust_N = 0;
0473         if (Thrust_tmp > Thrust) {
0474           Thrust = Thrust_tmp;
0475           n.SetXYZ(n_tmp.x(), n_tmp.y(), n_tmp.z());
0476         }
0477       }
0478     }
0479 
0480     // Calculation of the thrust major :
0481     TVector3 nT = n.Orthogonal();
0482     nT = nT.Unit();
0483     for (unsigned int i = 0; i < 720; i++) {
0484       nT.Rotate((2 * PI / 720) * i, n);
0485       for (unsigned int j = 0; j < TopJets.size(); j++) {
0486         Thrust_N += fabs(TopJets[j].px() * (nT.x()) + TopJets[j].py() * (nT.y()) + TopJets[j].pz() * (nT.z()));
0487       }
0488       Thrust_N += fabs(Lept->Px() * nT.x() + Lept->Py() * (nT.y()) + Lept->Pz() * (nT.z()));
0489 
0490       Thrust_Major_tmp = Thrust_N / Thrust_D;
0491       Thrust_N = 0;
0492 
0493       if (Thrust_Major_tmp > Thrust_Major) {
0494         Thrust_Major = Thrust_Major_tmp;
0495       }
0496     }
0497 
0498     // Calculation of the thrust minor :
0499 
0500     TVector3 nMinor = nT.Cross(n);
0501     nMinor = nMinor.Unit();
0502 
0503     for (unsigned int i = 0; i < TopJets.size(); i++) {
0504       Thrust_N +=
0505           fabs(TopJets[i].px() * (nMinor.x()) + TopJets[i].py() * (nMinor.y()) + TopJets[i].pz() * (nMinor.z()));
0506     }
0507     Thrust_N += fabs(Lept->Px() * nMinor.x() + Lept->Py() * (nMinor.y()) + Lept->Pz() * (nMinor.z()));
0508 
0509     Thrust_Minor = Thrust_N / Thrust_D;
0510   }
0511 
0512   double Obs31 = Thrust;
0513   evtselectVarVal.push_back(std::pair<unsigned int, double>(31, Obs31));
0514   if (DEBUG)
0515     std::cout << "------ LR observable 31 " << Obs31 << " calculated ------" << std::endl;
0516   double Obs32 = Thrust_Major;
0517   evtselectVarVal.push_back(std::pair<unsigned int, double>(32, Obs32));
0518   if (DEBUG)
0519     std::cout << "------ LR observable 32 " << Obs32 << " calculated ------" << std::endl;
0520   double Obs33 = Thrust_Minor;
0521   evtselectVarVal.push_back(std::pair<unsigned int, double>(33, Obs33));
0522   if (DEBUG)
0523     std::cout << "------ LR observable 33 " << Obs33 << " calculated ------" << std::endl;
0524 
0525   // Oblateness
0526 
0527   double Obs34 = Thrust_Major - Thrust_Minor;
0528   evtselectVarVal.push_back(std::pair<unsigned int, double>(34, Obs34));
0529   if (DEBUG)
0530     std::cout << "------ LR observable 34 " << Obs34 << " calculated ------" << std::endl;
0531 
0532   // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino
0533 
0534   TMatrixDSym Matrix_NoNu(3);
0535 
0536   double PX2_NoNu = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
0537                     std::pow(Lepb->Px(), 2) + std::pow(Lept->Px(), 2);
0538   double PY2_NoNu = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
0539                     std::pow(Lepb->Py(), 2) + std::pow(Lept->Py(), 2);
0540   double PZ2_NoNu = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
0541                     std::pow(Lepb->Pz(), 2) + std::pow(Lept->Pz(), 2);
0542 
0543   double P2_NoNu = PX2_NoNu + PY2_NoNu + PZ2_NoNu;
0544 
0545   double PXY_NoNu = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
0546                     Lepb->Px() * Lepb->Py() + Lept->Px() * Lept->Py();
0547   double PXZ_NoNu = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
0548                     Lepb->Px() * Lepb->Pz() + Lept->Px() * Lept->Pz();
0549   double PYZ_NoNu = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
0550                     Lepb->Py() * Lepb->Pz() + Lept->Py() * Lept->Pz();
0551 
0552   Matrix_NoNu(0, 0) = PX2_NoNu / P2_NoNu;
0553   Matrix_NoNu(0, 1) = PXY_NoNu / P2_NoNu;
0554   Matrix_NoNu(0, 2) = PXZ_NoNu / P2_NoNu;
0555   Matrix_NoNu(1, 0) = PXY_NoNu / P2_NoNu;
0556   Matrix_NoNu(1, 1) = PY2_NoNu / P2_NoNu;
0557   Matrix_NoNu(1, 2) = PYZ_NoNu / P2_NoNu;
0558   Matrix_NoNu(2, 0) = PXZ_NoNu / P2_NoNu;
0559   Matrix_NoNu(2, 1) = PYZ_NoNu / P2_NoNu;
0560   Matrix_NoNu(2, 2) = PZ2_NoNu / P2_NoNu;
0561 
0562   TMatrixDSymEigen pTensor_NoNu(Matrix_NoNu);
0563 
0564   std::vector<double> EigValues_NoNu;
0565   EigValues_NoNu.clear();
0566   for (int i = 0; i < 3; i++) {
0567     EigValues_NoNu.push_back(pTensor_NoNu.GetEigenValues()[i]);
0568   }
0569 
0570   std::sort(EigValues_NoNu.begin(), EigValues_NoNu.end(), dComparator);
0571 
0572   double Sphericity_NoNu = 1.5 * (EigValues_NoNu[1] + EigValues_NoNu[2]);
0573   double Aplanarity_NoNu = 1.5 * EigValues_NoNu[2];
0574 
0575   double Obs35 = (edm::isNotFinite(Sphericity_NoNu) ? 0 : Sphericity_NoNu);
0576   evtselectVarVal.push_back(std::pair<unsigned int, double>(35, Obs35));
0577   if (DEBUG)
0578     std::cout << "------ LR observable 35 " << Obs35 << " calculated ------" << std::endl;
0579 
0580   double Obs36 = (edm::isNotFinite(Aplanarity_NoNu) ? 0 : Aplanarity_NoNu);
0581   evtselectVarVal.push_back(std::pair<unsigned int, double>(36, Obs36));
0582   if (DEBUG)
0583     std::cout << "------ LR observable 36 " << Obs36 << " calculated ------" << std::endl;
0584 
0585   // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino or lepton
0586 
0587   TMatrixDSym Matrix_NoNuNoLep(3);
0588 
0589   double PX2_NoNuNoLep =
0590       std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Lepb->Px(), 2);
0591   double PY2_NoNuNoLep =
0592       std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Lepb->Py(), 2);
0593   double PZ2_NoNuNoLep =
0594       std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Lepb->Pz(), 2);
0595 
0596   double P2_NoNuNoLep = PX2_NoNuNoLep + PY2_NoNuNoLep + PZ2_NoNuNoLep;
0597 
0598   double PXY_NoNuNoLep =
0599       Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Lepb->Px() * Lepb->Py();
0600   double PXZ_NoNuNoLep =
0601       Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Lepb->Px() * Lepb->Pz();
0602   double PYZ_NoNuNoLep =
0603       Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Lepb->Py() * Lepb->Pz();
0604 
0605   Matrix_NoNuNoLep(0, 0) = PX2_NoNuNoLep / P2_NoNuNoLep;
0606   Matrix_NoNuNoLep(0, 1) = PXY_NoNuNoLep / P2_NoNuNoLep;
0607   Matrix_NoNuNoLep(0, 2) = PXZ_NoNuNoLep / P2_NoNuNoLep;
0608   Matrix_NoNuNoLep(1, 0) = PXY_NoNuNoLep / P2_NoNuNoLep;
0609   Matrix_NoNuNoLep(1, 1) = PY2_NoNuNoLep / P2_NoNuNoLep;
0610   Matrix_NoNuNoLep(1, 2) = PYZ_NoNuNoLep / P2_NoNuNoLep;
0611   Matrix_NoNuNoLep(2, 0) = PXZ_NoNuNoLep / P2_NoNuNoLep;
0612   Matrix_NoNuNoLep(2, 1) = PYZ_NoNuNoLep / P2_NoNuNoLep;
0613   Matrix_NoNuNoLep(2, 2) = PZ2_NoNuNoLep / P2_NoNuNoLep;
0614 
0615   TMatrixDSymEigen pTensor_NoNuNoLep(Matrix_NoNuNoLep);
0616 
0617   std::vector<double> EigValues_NoNuNoLep;
0618   EigValues_NoNuNoLep.clear();
0619   for (int i = 0; i < 3; i++) {
0620     EigValues_NoNuNoLep.push_back(pTensor_NoNuNoLep.GetEigenValues()[i]);
0621   }
0622 
0623   std::sort(EigValues_NoNuNoLep.begin(), EigValues_NoNuNoLep.end(), dComparator);
0624 
0625   double Sphericity_NoNuNoLep = 1.5 * (EigValues_NoNuNoLep[1] + EigValues_NoNuNoLep[2]);
0626   double Aplanarity_NoNuNoLep = 1.5 * EigValues_NoNuNoLep[2];
0627 
0628   double Obs37 = (edm::isNotFinite(Sphericity_NoNuNoLep) ? 0 : Sphericity_NoNuNoLep);
0629   evtselectVarVal.push_back(std::pair<unsigned int, double>(37, Obs37));
0630   if (DEBUG)
0631     std::cout << "------ LR observable 37 " << Obs37 << " calculated ------" << std::endl;
0632 
0633   double Obs38 = (edm::isNotFinite(Aplanarity_NoNuNoLep) ? 0 : Aplanarity_NoNuNoLep);
0634   evtselectVarVal.push_back(std::pair<unsigned int, double>(38, Obs38));
0635   if (DEBUG)
0636     std::cout << "------ LR observable 38 " << Obs38 << " calculated ------" << std::endl;
0637 
0638   // Put the vector in the TtSemiEvtSolution
0639   TS.setLRSignalEvtObservables(evtselectVarVal);
0640   if (DEBUG)
0641     std::cout << "------  Observable values stored in the event  ------" << std::endl;
0642 
0643   delete Hadp;
0644   if (DEBUG)
0645     std::cout << "------     Pointer to Hadp deleted             ------" << std::endl;
0646   delete Hadq;
0647   if (DEBUG)
0648     std::cout << "------     Pointer to Hadq deleted             ------" << std::endl;
0649   delete Hadb;
0650   if (DEBUG)
0651     std::cout << "------     Pointer to Hadb deleted             ------" << std::endl;
0652   delete Lepb;
0653   if (DEBUG)
0654     std::cout << "------     Pointer to Lepb deleted      ------" << std::endl;
0655   delete Lept;
0656   if (DEBUG)
0657     std::cout << "------     Pointer to Lepn deleted      ------" << std::endl;
0658   delete Lepn;
0659   if (DEBUG)
0660     std::cout << "------     Pointer to Mez deleted       ------" << std::endl;
0661   delete Mez;
0662 
0663   if (DEBUG)
0664     std::cout << "------------ LR observables calculated -----------" << std::endl;
0665 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!