CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
TtHadLRSignalSelObservables.cc
Go to the documentation of this file.
1 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
2 #include "FWCore/Utilities/interface/isFinite.h"
3 
4 TtHadLRSignalSelObservables::TtHadLRSignalSelObservables() {}
5 
6 TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables() {}
7 
8 void TtHadLRSignalSelObservables::operator()(TtHadEvtSolution &TS) {
9  evtselectVarVal.clear();
10  std::vector<pat::Jet> topJets;
11  topJets.clear();
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());
17  topJets.push_back(TS.getHadbbar());
18 
19  //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...
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());
32 
33  //sort the topJets in Et
34  std::sort(topJets.begin(), topJets.end(), EtComparator);
35 
36  //Et-Sum of the lightest jets
37  double EtSum = topJets[5].et() + topJets[5].et();
38  double Obs1 = (EtSum > 0 ? EtSum : -1);
39  evtselectVarVal.push_back(std::pair<unsigned int, double>(1, Obs1));
40 
41  //Log of the difference in Bdisc between the 2nd and the 3rd jets (ordered in Bdisc) - does that still
42  //make sense for the fully hadronic channel as well?FIXME!!!
43 
44  //sort the jets in Bdiscriminant
45  std::sort(topJets.begin(), topJets.end(), BdiscComparator);
46 
47  double BGap = topJets[1].bDiscriminator("trackCountingJetTags") - topJets[2].bDiscriminator("trackCountingJetTags");
48  double Obs2 = (BGap > 0 ? log(BGap) : -1);
49  evtselectVarVal.push_back(std::pair<unsigned int, double>(2, Obs2));
50 
51  //Circularity of the event
52  double N = 0, D = 0, C_tmp = 0, C = 1000;
53  double nx, ny, nz;
54 
55  // 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))
56 
57  for (unsigned int i = 0; i < 6; i++) {
58  D += fabs(topJets[i].pt());
59  }
60 
61  if ((D > 0)) {
62  // Loop over all the unit vectors in the transverse plane in order to find the miminum :
63  for (unsigned int i = 0; i < 360; i++) {
64  nx = cos((2 * PI / 360) * i);
65  ny = sin((2 * PI / 360) * i);
66  nz = 0;
67  N = 0;
68  for (unsigned int i = 0; i < 4; i++) {
69  N += fabs(topJets[i].px() * nx + topJets[i].py() * ny + topJets[i].pz() * nz);
70  }
71  C_tmp = 2 * N / D;
72  if (C_tmp < C)
73  C = C_tmp;
74  }
75  }
76 
77  double Obs3 = (C != 1000 ? C : -1);
78  evtselectVarVal.push_back(std::pair<unsigned int, double>(3, Obs3));
79 
80  //HT variable (Et-sum of the six jets)
81  double HT = 0;
82  for (unsigned int i = 0; i < 6; i++) {
83  HT += topJets[i].et();
84  }
85 
86  double Obs4 = (HT != 0 ? HT : -1);
87  evtselectVarVal.push_back(std::pair<unsigned int, double>(4, Obs4));
88 
89  //Transverse Mass of the system
90  math::XYZTLorentzVector pjets;
91  // for the six jets
92  for (unsigned int i = 0; i < 6; i++) {
93  pjets += topJets[i].p4();
94  }
95  double MT = sqrt(std::pow(pjets.mass(), 2));
96  double Obs5 = (MT > 0 ? MT : -1);
97  evtselectVarVal.push_back(std::pair<unsigned int, double>(5, Obs5));
98 
99  //CosTheta(Hadp,Hadq) and CosTheta(Hadj,Hadq)
100  //sort the lightJets in Et
101  std::sort(topJets.begin(), topJets.end(), EtComparator);
102 
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();
119 
120  TLorentzVector *LightJet1 = new TLorentzVector();
121  TLorentzVector *LightJet2 = new TLorentzVector();
122  TLorentzVector *LightJet3 = new TLorentzVector();
123  TLorentzVector *LightJet4 = new TLorentzVector();
124 
125  LightJet1->SetPxPyPzE(px1, py1, pz1, E1);
126  //LightJet1->Boost(BoostBackToCM);
127  LightJet2->SetPxPyPzE(px2, py2, pz2, E2);
128  //LightJet2->Boost(BoostBackToCM);
129  LightJet3->SetPxPyPzE(px3, py3, pz3, E3);
130  //LightJet1->Boost(BoostBackToCM);
131  LightJet4->SetPxPyPzE(px4, py4, pz4, E4);
132  //LightJet2->Boost(BoostBackToCM);
133 
134  double CosTheta1 = cos(LightJet2->Angle(LightJet1->Vect()));
135  double CosTheta2 = cos(LightJet4->Angle(LightJet3->Vect()));
136 
137  double Obs6 = (-1 < CosTheta1 ? CosTheta1 : -2);
138  evtselectVarVal.push_back(std::pair<unsigned int, double>(6, Obs6));
139  double Obs7 = (-1 < CosTheta2 ? CosTheta2 : -2);
140  evtselectVarVal.push_back(std::pair<unsigned int, double>(7, Obs7));
141 
142  delete LightJet1;
143  delete LightJet2;
144  delete LightJet3;
145  delete LightJet4;
146 
147  // try to find out more powerful observables related to the b-disc
148  //sort the jets in Bdiscriminant
149  std::sort(topJets.begin(), topJets.end(), BdiscComparator);
150 
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");
157  //std::cout<<"BjetsBdiscSum = "<<BjetsBdiscSum<<std::endl;
158  //std::cout<<"LjetsBdiscSum = "<<LjetsBdiscSum<<std::endl;
159 
160  double Obs8 = (Ljets1BdiscSum != 0 ? log(BjetsBdiscSum / Ljets1BdiscSum) : -1);
161  evtselectVarVal.push_back(std::pair<unsigned int, double>(8, Obs8));
162  double Obs9 = (Ljets2BdiscSum != 0 ? log(BjetsBdiscSum / Ljets2BdiscSum) : -1);
163  evtselectVarVal.push_back(std::pair<unsigned int, double>(9, Obs9));
164  double Obs10 = (BGap > 0 ? log(BjetsBdiscSum * BGap) : -1);
165  evtselectVarVal.push_back(std::pair<unsigned int, double>(10, Obs10));
166 
167  // Et-Ratio between the two highest in Et jets and four highest jets
168  double Obs11 = (HT != 0 ? (HT - EtSum) / HT : -1);
169  evtselectVarVal.push_back(std::pair<unsigned int, double>(11, Obs11));
170 
171  //Sphericity and Aplanarity without boosting back the system to CM frame
172 
173  TMatrixDSym Matrix(3);
174  double PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Hadj->Px(), 2) +
175  std::pow(Hadk->Px(), 2) + std::pow(Hadbbar->Px(), 2);
176  double PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Hadj->Py(), 2) +
177  std::pow(Hadk->Py(), 2) + std::pow(Hadbbar->Py(), 2);
178  double PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Hadj->Pz(), 2) +
179  std::pow(Hadk->Pz(), 2) + std::pow(Hadbbar->Pz(), 2);
180  double P2 = PX2 + PY2 + PZ2;
181 
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();
188 
189  Matrix(0, 0) = PX2 / P2;
190  Matrix(0, 1) = PXY / P2;
191  Matrix(0, 2) = PXZ / P2;
192  Matrix(1, 0) = PXY / P2;
193  Matrix(1, 1) = PY2 / P2;
194  Matrix(1, 2) = PYZ / P2;
195  Matrix(2, 0) = PXZ / P2;
196  Matrix(2, 1) = PYZ / P2;
197  Matrix(2, 2) = PZ2 / P2;
198 
199  TMatrixDSymEigen pTensor(Matrix);
200 
201  std::vector<double> EigValues;
202  EigValues.clear();
203  for (int i = 0; i < 3; i++) {
204  EigValues.push_back(pTensor.GetEigenValues()[i]);
205  }
206 
207  std::sort(EigValues.begin(), EigValues.end(), dComparator);
208 
209  double Sphericity = 1.5 * (EigValues[1] + EigValues[2]);
210  double Aplanarity = 1.5 * EigValues[2];
211 
212  double Obs12 = (edm::isNotFinite(Sphericity) ? -1 : Sphericity);
213  evtselectVarVal.push_back(std::pair<unsigned int, double>(12, Obs12));
214 
215  double Obs13 = (edm::isNotFinite(Aplanarity) ? -1 : Aplanarity);
216  evtselectVarVal.push_back(std::pair<unsigned int, double>(13, Obs13));
217 
218  //Sphericity and Aplanarity with boosting back the system to CM frame
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());
224 
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);
232 
233  double BOOST_PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
234  std::pow(Hadj->Px(), 2) + std::pow(Hadk->Px(), 2) + std::pow(Hadbbar->Px(), 2);
235  double BOOST_PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
236  std::pow(Hadj->Py(), 2) + std::pow(Hadk->Py(), 2) + std::pow(Hadbbar->Py(), 2);
237  double BOOST_PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
238  std::pow(Hadj->Pz(), 2) + std::pow(Hadk->Pz(), 2) + std::pow(Hadbbar->Pz(), 2);
239 
240  double BOOST_P2 = BOOST_PX2 + BOOST_PY2 + BOOST_PZ2;
241 
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();
248 
249  TMatrixDSym BOOST_Matrix(3);
250 
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;
260 
261  TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
262 
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]);
267  }
268 
269  std::sort(BOOST_EigValues.begin(), BOOST_EigValues.end(), dComparator);
270 
271  double BOOST_Sphericity = 1.5 * (BOOST_EigValues[1] + BOOST_EigValues[2]);
272  double BOOST_Aplanarity = 1.5 * BOOST_EigValues[2];
273 
274  double Obs14 = (edm::isNotFinite(BOOST_Sphericity) ? -1 : BOOST_Sphericity);
275  evtselectVarVal.push_back(std::pair<unsigned int, double>(14, Obs14));
276 
277  double Obs15 = (edm::isNotFinite(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity);
278  evtselectVarVal.push_back(std::pair<unsigned int, double>(15, Obs15));
279 
280  // Centrality of the six jets
281  double H = 0;
282  for (unsigned int i = 0; i < 6; i++) {
283  H += topJets[i].energy();
284  }
285  double Obs16 = (H != 0 ? HT / H : -1);
286  evtselectVarVal.push_back(std::pair<unsigned int, double>(16, Obs16));
287 
288  // distance from the origin in the (BjetsBdiscSum, Ljets1BdiscSum) and (BjetsBdiscSum, Ljets2BdiscSum)
289  double Obs17 = (BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets1BdiscSum) : -1);
290  evtselectVarVal.push_back(std::pair<unsigned int, double>(17, Obs17));
291  double Obs18 = (BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707 * (BjetsBdiscSum - Ljets2BdiscSum) : -1);
292  evtselectVarVal.push_back(std::pair<unsigned int, double>(18, Obs18));
293 
294  // Ratio of the Et Sum of the two lightest jets over HT=Sum of the Et of the six highest jets in Et
295  double Obs19 = (HT != 0 ? EtSum / HT : -1);
296  evtselectVarVal.push_back(std::pair<unsigned int, double>(19, Obs19));
297 
298  // Put the vector in the TtHadEvtSolution
299  TS.setLRSignalEvtObservables(evtselectVarVal);
300 }
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
mps_fire.i
i
Definition: mps_fire.py:428
DiDispStaMuonMonitor_cfi.pt
tuple pt
Definition: DiDispStaMuonMonitor_cfi.py:39
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TtHadEvtSolution::getHadbbar
pat::Jet getHadbbar() const
Definition: TtHadEvtSolution.cc:73
TtHadEvtSolution::setLRSignalEvtObservables
void setLRSignalEvtObservables(const std::vector< std::pair< unsigned int, double > > &varval)
Definition: TtHadEvtSolution.cc:228
TtHadEvtSolution::getHadb
pat::Jet getHadb() const
Definition: TtHadEvtSolution.cc:34
TtHadEvtSolution::getHadq
pat::Jet getHadq() const
Definition: TtHadEvtSolution.cc:60
TtHadEvtSolution::getHadj
pat::Jet getHadj() const
Definition: TtHadEvtSolution.cc:86
TtHadLRSignalSelObservables::evtselectVarVal
std::vector< std::pair< unsigned int, double > > evtselectVarVal
Definition: TtHadLRSignalSelObservables.h:57
hitfit::Matrix
CLHEP::HepMatrix Matrix
Definition: matutil.h:62
WZElectronSkims53X_cff.MT
string MT
Definition: WZElectronSkims53X_cff.py:219
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TtHadEvtSolution::getHadp
pat::Jet getHadp() const
Definition: TtHadEvtSolution.cc:47
PI
#define PI
Definition: QcdUeDQM.h:37
TtHadEvtSolution::getHadk
pat::Jet getHadk() const
Definition: TtHadEvtSolution.cc:99
N
#define N
Definition: blowfish.cc:9
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
HT
Definition: HT.h:21
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Generated for CMSSW Reference Manual by doxygen 1.8.5