CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
TtSemiLRSignalSelObservables.cc
Go to the documentation of this file.
1 #include "TopQuarkAnalysis/TopEventSelection/interface/TtSemiLRSignalSelObservables.h"
2 #include "FWCore/Utilities/interface/isFinite.h"
3 
4 TtSemiLRSignalSelObservables::TtSemiLRSignalSelObservables() {}
5 
6 TtSemiLRSignalSelObservables::~TtSemiLRSignalSelObservables() {}
7 
8 void TtSemiLRSignalSelObservables::operator()(TtSemiEvtSolution &TS, const std::vector<pat::Jet> &SelectedTopJets) {
9  evtselectVarVal.clear();
10 
11  // choice the B-tag algorithm
12  const char *BtagAlgo = "trackCountingJetTags";
13 
14  // activate the "debug mode"
15  bool DEBUG = false;
16 
17  if (DEBUG)
18  std::cout << "---------- Start calculating the LR observables ----------" << std::endl;
19 
20  std::vector<pat::Jet> TopJets;
21  TopJets.clear();
22 
23  for (unsigned int i = 0; i < SelectedTopJets.size(); i++) {
24  TopJets.push_back(SelectedTopJets[i]);
25  }
26 
27  //sort the TopJets in Et
28  std::sort(TopJets.begin(), TopJets.end(), EtComparator);
29 
30  TLorentzVector *Hadp = new TLorentzVector();
31  Hadp->SetPxPyPzE(TopJets[3].px(), TopJets[3].py(), TopJets[3].pz(), TopJets[3].energy());
32 
33  TLorentzVector *Hadq = new TLorentzVector();
34  Hadq->SetPxPyPzE(TopJets[2].px(), TopJets[2].py(), TopJets[2].pz(), TopJets[2].energy());
35 
36  TLorentzVector *Hadb = new TLorentzVector();
37  Hadb->SetPxPyPzE(TopJets[1].px(), TopJets[1].py(), TopJets[1].pz(), TopJets[1].energy());
38 
39  TLorentzVector *Lepb = new TLorentzVector();
40  Lepb->SetPxPyPzE(TopJets[0].px(), TopJets[0].py(), TopJets[0].pz(), TopJets[0].energy());
41 
42  TLorentzVector *Lept = new TLorentzVector();
43  if (TS.getDecay() == "muon")
44  Lept->SetPxPyPzE(TS.getMuon().px(), TS.getMuon().py(), TS.getMuon().pz(), TS.getMuon().energy());
45  else if (TS.getDecay() == "electron")
46  Lept->SetPxPyPzE(TS.getElectron().px(), TS.getElectron().py(), TS.getElectron().pz(), TS.getElectron().energy());
47 
48  TLorentzVector *Lepn = new TLorentzVector();
49  Lepn->SetPxPyPzE(TS.getNeutrino().px(), TS.getNeutrino().py(), TS.getNeutrino().pz(), TS.getNeutrino().energy());
50 
51  // Calculation of the pz of the neutrino due to W-mass constraint
52 
53  MEzCalculator *Mez = new MEzCalculator();
54  Mez->SetMET(TS.getNeutrino());
55  if (TS.getDecay() == "muon")
56  Mez->SetLepton(TS.getMuon());
57  else
58  Mez->SetLepton(TS.getElectron(), false);
59  double MEZ = Mez->Calculate();
60  Lepn->SetPz(MEZ);
61 
62  // Pt of the lepton
63 
64  double Obs1 = Lept->Pt();
65  evtselectVarVal.push_back(std::pair<unsigned int, double>(1, Obs1));
66  if (DEBUG)
67  std::cout << "------ LR observable 1 " << Obs1 << " calculated ------" << std::endl;
68 
69  // Missing transverse energy
70 
71  double Obs2 = TS.getNeutrino().et();
72  evtselectVarVal.push_back(std::pair<unsigned int, double>(2, Obs2));
73  if (DEBUG)
74  std::cout << "------ LR observable 2 " << Obs2 << " calculated ------" << std::endl;
75 
76  //HT variable (Et-sum of the four jets)
77 
78  double HT = 0;
79  for (unsigned int i = 0; i < 4; i++) {
80  HT += TopJets[i].et();
81  }
82 
83  double Obs3 = HT;
84  evtselectVarVal.push_back(std::pair<unsigned int, double>(3, Obs3));
85  if (DEBUG)
86  std::cout << "------ LR observable 3 " << Obs3 << " calculated ------" << std::endl;
87 
88  //Et-Sum of the lightest jets
89 
90  double EtSum = TopJets[2].et() + TopJets[3].et();
91  double Obs4 = EtSum;
92  evtselectVarVal.push_back(std::pair<unsigned int, double>(4, Obs4));
93  if (DEBUG)
94  std::cout << "------ LR observable 4 " << Obs4 << " calculated ------" << std::endl;
95 
96  // Et-Ratio between the two lowest jets in Et and four highest jets
97 
98  double Obs5 = EtSum / HT;
99  evtselectVarVal.push_back(std::pair<unsigned int, double>(5, Obs5));
100  if (DEBUG)
101  std::cout << "------ LR observable 5 " << Obs5 << " calculated ------" << std::endl;
102 
103  // Et-Ratio between the two highest jets in Et and four highest jets
104 
105  double Obs6 = (HT - EtSum) / HT;
106  evtselectVarVal.push_back(std::pair<unsigned int, double>(6, Obs6));
107  if (DEBUG)
108  std::cout << "------ LR observable 6 " << Obs6 << " calculated ------" << std::endl;
109 
110  // Transverse Mass of the system
111 
112  TLorentzVector TtbarSystem = (*Hadp) + (*Hadq) + (*Hadb) + (*Lepb) + (*Lept) + (*Lepn);
113  double MT = TtbarSystem.Mt();
114  double Obs7 = MT;
115  evtselectVarVal.push_back(std::pair<unsigned int, double>(7, Obs7));
116  if (DEBUG)
117  std::cout << "------ LR observable 7 " << Obs7 << " calculated ------" << std::endl;
118 
119  // Observables related to the b-disc (jets ordered in Bdisc)
120 
121  //sort the TopJets in Bdiscriminant
122  std::sort(TopJets.begin(), TopJets.end(), BdiscComparator);
123 
124  double Obs8;
125  double Obs9;
126  double Obs10;
127  double Obs11;
128 
129  // Difference in bdisc between the 2nd and the 3rd jets
130  double BGap = TopJets[1].bDiscriminator(BtagAlgo) - TopJets[2].bDiscriminator(BtagAlgo);
131 
132  // Sum of the bdisc of the two highest/lowest jets
133  double BjetsBdiscSum = TopJets[0].bDiscriminator(BtagAlgo) + TopJets[1].bDiscriminator(BtagAlgo);
134  double LjetsBdiscSum = TopJets[2].bDiscriminator(BtagAlgo) + TopJets[3].bDiscriminator(BtagAlgo);
135 
136  Obs8 = (TopJets[2].bDiscriminator(BtagAlgo) > -10 ? log(BGap) : -5);
137  if (DEBUG)
138  std::cout << "------ LR observable 8 " << Obs8 << " calculated ------" << std::endl;
139  Obs9 = (BjetsBdiscSum * BGap);
140  if (DEBUG)
141  std::cout << "------ LR observable 9 " << Obs9 << " calculated ------" << std::endl;
142  Obs10 = (BjetsBdiscSum / LjetsBdiscSum);
143  if (DEBUG)
144  std::cout << "------ LR observable 10 " << Obs10 << " calculated ------" << std::endl;
145  // Distance from the origin in the (BjetsBdiscSum, LjetsBdiscSum) plane
146  Obs11 = 0.707 * ((BjetsBdiscSum + LjetsBdiscSum) / 2 + 10);
147  if (DEBUG)
148  std::cout << "------ LR observable 11 " << Obs11 << " calculated ------" << std::endl;
149 
150  evtselectVarVal.push_back(std::pair<unsigned int, double>(8, Obs8));
151  evtselectVarVal.push_back(std::pair<unsigned int, double>(9, Obs9));
152  evtselectVarVal.push_back(std::pair<unsigned int, double>(10, Obs10));
153  evtselectVarVal.push_back(std::pair<unsigned int, double>(11, Obs11));
154 
155  //sort the TopJets in Et
156  std::sort(TopJets.begin(), TopJets.end(), EtComparator);
157 
158  // Circularity of the event
159 
160  double N = 0, D = 0, C_tmp = 0, C = 10000;
161  double N_NoNu = 0, D_NoNu = 0, C_NoNu_tmp = 0, C_NoNu = 10000;
162  double nx, ny, nz;
163 
164  // 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))
165 
166  for (unsigned int i = 0; i < 4; i++) {
167  D += fabs(TopJets[i].pt());
168  }
169  // modified the 26th of September to calculate also the circularity without the neutrino contribution
170  D += fabs(Lept->Pt());
171  D_NoNu = D;
172  D += fabs(Lepn->Pt());
173 
174  if ((D > 0)) {
175  // Loop over all the unit vectors in the transverse plane in order to find the miminum :
176  for (unsigned int i = 0; i < 720; i++) {
177  nx = cos((2 * PI / 720) * i);
178  ny = sin((2 * PI / 720) * i);
179  nz = 0;
180  N = 0;
181 
182  for (unsigned int i = 0; i < 4; i++) {
183  N += fabs(TopJets[i].px() * nx + TopJets[i].py() * ny + TopJets[i].pz() * nz);
184  }
185  // modified the 26th of September to calculate also the circularity without the neutrino contribution
186  N += fabs(Lept->Px() * nx + Lept->Py() * ny + Lept->Pz() * nz);
187  N_NoNu = N;
188  N += fabs(Lepn->Px() * nx + Lepn->Py() * ny + Lepn->Pz() * nz);
189 
190  C_tmp = 2 * N / D;
191 
192  // modified the 26th of September to calculate also the circularity without the neutrino contribution
193  C_NoNu_tmp = 2 * N_NoNu / D_NoNu;
194 
195  if (C_tmp < C)
196  C = C_tmp;
197 
198  // modified the 26th of September to calculate also the circularity without the neutrino contribution
199  if (C_NoNu_tmp < C_NoNu)
200  C_NoNu = C_NoNu_tmp;
201  }
202  }
203 
204  double Obs12 = (C != 10000 ? C : 0);
205  evtselectVarVal.push_back(std::pair<unsigned int, double>(12, Obs12));
206  if (DEBUG)
207  std::cout << "------ LR observable 12 " << Obs12 << " calculated ------" << std::endl;
208 
209  // Circularity of the event without neutrino
210 
211  double Obs13 = (C_NoNu != 10000 ? C_NoNu : 0);
212  evtselectVarVal.push_back(std::pair<unsigned int, double>(13, Obs13));
213  if (DEBUG)
214  std::cout << "------ LR observable 13 " << Obs13 << " calculated ------" << std::endl;
215 
216  // Centrality of the four highest jets
217 
218  double H = 0;
219  for (unsigned int i = 0; i < 4; i++) {
220  H += TopJets[i].energy();
221  }
222 
223  double Obs14 = HT / H;
224  evtselectVarVal.push_back(std::pair<unsigned int, double>(14, Obs14));
225  if (DEBUG)
226  std::cout << "------ LR observable 14 " << Obs14 << " calculated ------" << std::endl;
227 
228  // Sphericity and Aplanarity without boosting back the system to CM frame
229 
230  TMatrixDSym Matrix(3);
231 
232  double PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Lepb->Px(), 2) +
233  std::pow(Lept->Px(), 2) + std::pow(Lepn->Px(), 2);
234  double PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Lepb->Py(), 2) +
235  std::pow(Lept->Py(), 2) + std::pow(Lepn->Py(), 2);
236  double PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Lepb->Pz(), 2) +
237  std::pow(Lept->Pz(), 2) + std::pow(Lepn->Pz(), 2);
238 
239  double P2 = PX2 + PY2 + PZ2;
240 
241  double PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Lepb->Px() * Lepb->Py() +
242  Lept->Px() * Lept->Py() + Lepn->Px() * Lepn->Py();
243  double PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Lepb->Px() * Lepb->Pz() +
244  Lept->Px() * Lept->Pz() + Lepn->Px() * Lepn->Pz();
245  double PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Lepb->Py() * Lepb->Pz() +
246  Lept->Py() * Lept->Pz() + Lepn->Py() * Lepn->Pz();
247 
248  Matrix(0, 0) = PX2 / P2;
249  Matrix(0, 1) = PXY / P2;
250  Matrix(0, 2) = PXZ / P2;
251  Matrix(1, 0) = PXY / P2;
252  Matrix(1, 1) = PY2 / P2;
253  Matrix(1, 2) = PYZ / P2;
254  Matrix(2, 0) = PXZ / P2;
255  Matrix(2, 1) = PYZ / P2;
256  Matrix(2, 2) = PZ2 / P2;
257 
258  TMatrixDSymEigen pTensor(Matrix);
259 
260  std::vector<double> EigValues;
261  EigValues.clear();
262  for (int i = 0; i < 3; i++) {
263  EigValues.push_back(pTensor.GetEigenValues()[i]);
264  }
265 
266  std::sort(EigValues.begin(), EigValues.end(), dComparator);
267 
268  double Sphericity = 1.5 * (EigValues[1] + EigValues[2]);
269  double Aplanarity = 1.5 * EigValues[2];
270 
271  double Obs15 = (edm::isNotFinite(Sphericity) ? 0 : Sphericity);
272  evtselectVarVal.push_back(std::pair<unsigned int, double>(15, Obs15));
273  if (DEBUG)
274  std::cout << "------ LR observable 15 " << Obs15 << " calculated ------" << std::endl;
275 
276  double Obs16 = (edm::isNotFinite(Aplanarity) ? 0 : Aplanarity);
277  evtselectVarVal.push_back(std::pair<unsigned int, double>(16, Obs16));
278  if (DEBUG)
279  std::cout << "------ LR observable 16 " << Obs16 << " calculated ------" << std::endl;
280 
281  // Sphericity and Aplanarity with boosting back the system to CM frame
282 
283  TVector3 BoostBackToCM = -(TtbarSystem.BoostVector());
284  Hadp->Boost(BoostBackToCM);
285  Hadq->Boost(BoostBackToCM);
286  Hadb->Boost(BoostBackToCM);
287  Lepb->Boost(BoostBackToCM);
288  Lept->Boost(BoostBackToCM);
289  Lepn->Boost(BoostBackToCM);
290 
291  double BOOST_PX2 = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
292  std::pow(Lepb->Px(), 2) + std::pow(Lept->Px(), 2) + std::pow(Lepn->Px(), 2);
293  double BOOST_PY2 = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
294  std::pow(Lepb->Py(), 2) + std::pow(Lept->Py(), 2) + std::pow(Lepn->Py(), 2);
295  double BOOST_PZ2 = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
296  std::pow(Lepb->Pz(), 2) + std::pow(Lept->Pz(), 2) + std::pow(Lepn->Pz(), 2);
297 
298  double BOOST_P2 = BOOST_PX2 + BOOST_PY2 + BOOST_PZ2;
299 
300  double BOOST_PXY = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
301  Lepb->Px() * Lepb->Py() + Lept->Px() * Lept->Py() + Lepn->Px() * Lepn->Py();
302  double BOOST_PXZ = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
303  Lepb->Px() * Lepb->Pz() + Lept->Px() * Lept->Pz() + Lepn->Px() * Lepn->Pz();
304  double BOOST_PYZ = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
305  Lepb->Py() * Lepb->Pz() + Lept->Py() * Lept->Pz() + Lepn->Py() * Lepn->Pz();
306 
307  TMatrixDSym BOOST_Matrix(3);
308 
309  BOOST_Matrix(0, 0) = BOOST_PX2 / BOOST_P2;
310  BOOST_Matrix(0, 1) = BOOST_PXY / BOOST_P2;
311  BOOST_Matrix(0, 2) = BOOST_PXZ / BOOST_P2;
312  BOOST_Matrix(1, 0) = BOOST_PXY / BOOST_P2;
313  BOOST_Matrix(1, 1) = BOOST_PY2 / BOOST_P2;
314  BOOST_Matrix(1, 2) = BOOST_PYZ / BOOST_P2;
315  BOOST_Matrix(2, 0) = BOOST_PXZ / BOOST_P2;
316  BOOST_Matrix(2, 1) = BOOST_PYZ / BOOST_P2;
317  BOOST_Matrix(2, 2) = BOOST_PZ2 / BOOST_P2;
318 
319  TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
320 
321  std::vector<double> BOOST_EigValues;
322  BOOST_EigValues.clear();
323  for (int i = 0; i < 3; i++) {
324  BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
325  }
326 
327  std::sort(BOOST_EigValues.begin(), BOOST_EigValues.end(), dComparator);
328 
329  double BOOST_Sphericity = 1.5 * (BOOST_EigValues[1] + BOOST_EigValues[2]);
330  double BOOST_Aplanarity = 1.5 * BOOST_EigValues[2];
331 
332  double Obs17 = (edm::isNotFinite(BOOST_Sphericity) ? 0 : BOOST_Sphericity);
333  evtselectVarVal.push_back(std::pair<unsigned int, double>(17, Obs17));
334  if (DEBUG)
335  std::cout << "------ LR observable 17 " << Obs17 << " calculated ------" << std::endl;
336 
337  double Obs18 = (edm::isNotFinite(BOOST_Aplanarity) ? 0 : BOOST_Aplanarity);
338  evtselectVarVal.push_back(std::pair<unsigned int, double>(18, Obs18));
339  if (DEBUG)
340  std::cout << "------ LR observable 18 " << Obs18 << " calculated ------" << std::endl;
341 
342  //ratio between ET of the fifth jet and HT
343 
344  double Obs19 = (TopJets.size() > 4) ? TopJets[4].et() / HT : 1.0;
345  evtselectVarVal.push_back(std::pair<unsigned int, double>(19, Obs19));
346  if (DEBUG)
347  std::cout << "------ LR observable 19 " << Obs19 << " calculated ------" << std::endl;
348 
349  // HT variable calculated with all the jets in the event.
350 
351  double HT_alljets = 0;
352  double H_alljets = 0;
353  for (unsigned int i = 0; i < TopJets.size(); i++) {
354  HT_alljets += TopJets[i].et();
355  H_alljets += TopJets[i].energy();
356  }
357  double Obs20 = HT_alljets;
358  evtselectVarVal.push_back(std::pair<unsigned int, double>(20, Obs20));
359  if (DEBUG)
360  std::cout << "------ LR observable 20 " << Obs20 << " calculated ------" << std::endl;
361 
362  // HT3 = HT calculated with all jets except the two leading jets
363 
364  double HT3_alljets = 0;
365  for (unsigned int i = 2; i < TopJets.size(); i++) {
366  HT3_alljets += TopJets[i].et();
367  }
368  double Obs21 = HT3_alljets;
369  evtselectVarVal.push_back(std::pair<unsigned int, double>(21, Obs21));
370  if (DEBUG)
371  std::cout << "------ LR observable 21 " << Obs21 << " calculated ------" << std::endl;
372 
373  // ET0, ratio of the Et of the leading and HT_alljets
374 
375  double ET0 = TopJets[0].et() / HT_alljets;
376  double Obs22 = ET0;
377  evtselectVarVal.push_back(std::pair<unsigned int, double>(22, Obs22));
378  if (DEBUG)
379  std::cout << "------ LR observable 22 " << Obs22 << " calculated ------" << std::endl;
380 
381  // Centrality of the event computed with all jets
382 
383  double Obs23 = ((H_alljets > 0) ? HT_alljets / H_alljets : 0);
384  evtselectVarVal.push_back(std::pair<unsigned int, double>(23, Obs23));
385  if (DEBUG)
386  std::cout << "------ LR observable 23 " << Obs23 << " calculated ------" << std::endl;
387 
388  //Fox-Wolfram momenta (1st to 6th), modified for hadron collider and using a Legendre polynomials expansion
389 
390  double FW_momentum_0 = 0, FW_momentum_1 = 0, FW_momentum_2 = 0, FW_momentum_3 = 0, FW_momentum_4 = 0,
391  FW_momentum_5 = 0, FW_momentum_6 = 0;
392 
393  for (unsigned int i = 0; i < TopJets.size(); i++) {
394  for (unsigned int j = 0; j < TopJets.size(); j++) {
395  double ET_ij_over_ETSum2 = TopJets[i].et() * TopJets[j].et() / (std::pow(HT_alljets, 2));
396  double cosTheta_ij =
397  (TopJets[i].px() * TopJets[j].px() + TopJets[i].py() * TopJets[j].py() + TopJets[i].pz() * TopJets[j].pz()) /
398  (TopJets[i].p4().R() * TopJets[j].p4().R());
399  FW_momentum_0 += ET_ij_over_ETSum2;
400  FW_momentum_1 += ET_ij_over_ETSum2 * cosTheta_ij;
401  FW_momentum_2 += ET_ij_over_ETSum2 * 0.5 * (3 * std::pow(cosTheta_ij, 2) - 1);
402  FW_momentum_3 += ET_ij_over_ETSum2 * 0.5 * (5 * std::pow(cosTheta_ij, 3) - 3 * cosTheta_ij);
403  FW_momentum_4 += ET_ij_over_ETSum2 * 0.125 * (35 * std::pow(cosTheta_ij, 4) - 30 * std::pow(cosTheta_ij, 2) + 3);
404  FW_momentum_5 += ET_ij_over_ETSum2 * 0.125 *
405  (63 * std::pow(cosTheta_ij, 5) - 70 * std::pow(cosTheta_ij, 3) + 15 * cosTheta_ij);
406  FW_momentum_6 +=
407  ET_ij_over_ETSum2 * 0.0625 *
408  (231 * std::pow(cosTheta_ij, 6) - 315 * std::pow(cosTheta_ij, 4) + 105 * std::pow(cosTheta_ij, 2) - 5);
409  }
410  }
411 
412  double Obs24 = FW_momentum_0;
413  evtselectVarVal.push_back(std::pair<unsigned int, double>(24, Obs24));
414  if (DEBUG)
415  std::cout << "------ LR observable 24 " << Obs24 << " calculated ------" << std::endl;
416  double Obs25 = FW_momentum_1;
417  evtselectVarVal.push_back(std::pair<unsigned int, double>(25, Obs25));
418  if (DEBUG)
419  std::cout << "------ LR observable 25 " << Obs25 << " calculated ------" << std::endl;
420  double Obs26 = FW_momentum_2;
421  evtselectVarVal.push_back(std::pair<unsigned int, double>(26, Obs26));
422  if (DEBUG)
423  std::cout << "------ LR observable 26 " << Obs26 << " calculated ------" << std::endl;
424  double Obs27 = FW_momentum_3;
425  evtselectVarVal.push_back(std::pair<unsigned int, double>(27, Obs27));
426  if (DEBUG)
427  std::cout << "------ LR observable 27 " << Obs27 << " calculated ------" << std::endl;
428  double Obs28 = FW_momentum_4;
429  evtselectVarVal.push_back(std::pair<unsigned int, double>(28, Obs28));
430  if (DEBUG)
431  std::cout << "------ LR observable 28 " << Obs28 << " calculated ------" << std::endl;
432  double Obs29 = FW_momentum_5;
433  evtselectVarVal.push_back(std::pair<unsigned int, double>(29, Obs29));
434  if (DEBUG)
435  std::cout << "------ LR observable 29 " << Obs29 << " calculated ------" << std::endl;
436  double Obs30 = FW_momentum_6;
437  evtselectVarVal.push_back(std::pair<unsigned int, double>(30, Obs30));
438  if (DEBUG)
439  std::cout << "------ LR observable 30 " << Obs30 << " calculated ------" << std::endl;
440 
441  // Thrust, thrust major and thrust minor
442 
443  TVector3 n(0, 0, 0), n_tmp(0, 0, 0);
444 
445  double Thrust_D = 0, Thrust_N = 0;
446  double Thrust = -1, Thrust_tmp = 0;
447  double Thrust_Major = -1, Thrust_Major_tmp = 0;
448  double Thrust_Minor = -1;
449 
450  for (unsigned int i = 0; i < TopJets.size(); i++) {
451  Thrust_D += TopJets[i].p();
452  }
453 
454  Thrust_D += Lept->P();
455 
456  if (Thrust_D > 0) {
457  // Calculation of the thrust :
458  for (unsigned int i = 0; i < 720; i++) {
459  for (unsigned int j = 0; j < 360; j++) {
460  n_tmp.SetX(cos((2 * PI / 720) * i) * sin((PI / 360) * j));
461  n_tmp.SetY(sin((2 * PI / 720) * i) * sin((PI / 360) * j));
462  n_tmp.SetZ(cos((PI / 360) * j));
463 
464  for (unsigned int k = 0; k < TopJets.size(); k++) {
465  Thrust_N +=
466  fabs(TopJets[k].px() * (n_tmp.x()) + TopJets[k].py() * (n_tmp.y()) + TopJets[k].pz() * (n_tmp.z()));
467  }
468  Thrust_N += fabs(Lept->Px() * (n_tmp.x()) + Lept->Py() * (n_tmp.y()) + Lept->Pz() * (n_tmp.z()));
469 
470  Thrust_tmp = Thrust_N / Thrust_D;
471 
472  Thrust_N = 0;
473  if (Thrust_tmp > Thrust) {
474  Thrust = Thrust_tmp;
475  n.SetXYZ(n_tmp.x(), n_tmp.y(), n_tmp.z());
476  }
477  }
478  }
479 
480  // Calculation of the thrust major :
481  TVector3 nT = n.Orthogonal();
482  nT = nT.Unit();
483  for (unsigned int i = 0; i < 720; i++) {
484  nT.Rotate((2 * PI / 720) * i, n);
485  for (unsigned int j = 0; j < TopJets.size(); j++) {
486  Thrust_N += fabs(TopJets[j].px() * (nT.x()) + TopJets[j].py() * (nT.y()) + TopJets[j].pz() * (nT.z()));
487  }
488  Thrust_N += fabs(Lept->Px() * nT.x() + Lept->Py() * (nT.y()) + Lept->Pz() * (nT.z()));
489 
490  Thrust_Major_tmp = Thrust_N / Thrust_D;
491  Thrust_N = 0;
492 
493  if (Thrust_Major_tmp > Thrust_Major) {
494  Thrust_Major = Thrust_Major_tmp;
495  }
496  }
497 
498  // Calculation of the thrust minor :
499 
500  TVector3 nMinor = nT.Cross(n);
501  nMinor = nMinor.Unit();
502 
503  for (unsigned int i = 0; i < TopJets.size(); i++) {
504  Thrust_N +=
505  fabs(TopJets[i].px() * (nMinor.x()) + TopJets[i].py() * (nMinor.y()) + TopJets[i].pz() * (nMinor.z()));
506  }
507  Thrust_N += fabs(Lept->Px() * nMinor.x() + Lept->Py() * (nMinor.y()) + Lept->Pz() * (nMinor.z()));
508 
509  Thrust_Minor = Thrust_N / Thrust_D;
510  }
511 
512  double Obs31 = Thrust;
513  evtselectVarVal.push_back(std::pair<unsigned int, double>(31, Obs31));
514  if (DEBUG)
515  std::cout << "------ LR observable 31 " << Obs31 << " calculated ------" << std::endl;
516  double Obs32 = Thrust_Major;
517  evtselectVarVal.push_back(std::pair<unsigned int, double>(32, Obs32));
518  if (DEBUG)
519  std::cout << "------ LR observable 32 " << Obs32 << " calculated ------" << std::endl;
520  double Obs33 = Thrust_Minor;
521  evtselectVarVal.push_back(std::pair<unsigned int, double>(33, Obs33));
522  if (DEBUG)
523  std::cout << "------ LR observable 33 " << Obs33 << " calculated ------" << std::endl;
524 
525  // Oblateness
526 
527  double Obs34 = Thrust_Major - Thrust_Minor;
528  evtselectVarVal.push_back(std::pair<unsigned int, double>(34, Obs34));
529  if (DEBUG)
530  std::cout << "------ LR observable 34 " << Obs34 << " calculated ------" << std::endl;
531 
532  // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino
533 
534  TMatrixDSym Matrix_NoNu(3);
535 
536  double PX2_NoNu = std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) +
537  std::pow(Lepb->Px(), 2) + std::pow(Lept->Px(), 2);
538  double PY2_NoNu = std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) +
539  std::pow(Lepb->Py(), 2) + std::pow(Lept->Py(), 2);
540  double PZ2_NoNu = std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) +
541  std::pow(Lepb->Pz(), 2) + std::pow(Lept->Pz(), 2);
542 
543  double P2_NoNu = PX2_NoNu + PY2_NoNu + PZ2_NoNu;
544 
545  double PXY_NoNu = Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() +
546  Lepb->Px() * Lepb->Py() + Lept->Px() * Lept->Py();
547  double PXZ_NoNu = Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() +
548  Lepb->Px() * Lepb->Pz() + Lept->Px() * Lept->Pz();
549  double PYZ_NoNu = Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() +
550  Lepb->Py() * Lepb->Pz() + Lept->Py() * Lept->Pz();
551 
552  Matrix_NoNu(0, 0) = PX2_NoNu / P2_NoNu;
553  Matrix_NoNu(0, 1) = PXY_NoNu / P2_NoNu;
554  Matrix_NoNu(0, 2) = PXZ_NoNu / P2_NoNu;
555  Matrix_NoNu(1, 0) = PXY_NoNu / P2_NoNu;
556  Matrix_NoNu(1, 1) = PY2_NoNu / P2_NoNu;
557  Matrix_NoNu(1, 2) = PYZ_NoNu / P2_NoNu;
558  Matrix_NoNu(2, 0) = PXZ_NoNu / P2_NoNu;
559  Matrix_NoNu(2, 1) = PYZ_NoNu / P2_NoNu;
560  Matrix_NoNu(2, 2) = PZ2_NoNu / P2_NoNu;
561 
562  TMatrixDSymEigen pTensor_NoNu(Matrix_NoNu);
563 
564  std::vector<double> EigValues_NoNu;
565  EigValues_NoNu.clear();
566  for (int i = 0; i < 3; i++) {
567  EigValues_NoNu.push_back(pTensor_NoNu.GetEigenValues()[i]);
568  }
569 
570  std::sort(EigValues_NoNu.begin(), EigValues_NoNu.end(), dComparator);
571 
572  double Sphericity_NoNu = 1.5 * (EigValues_NoNu[1] + EigValues_NoNu[2]);
573  double Aplanarity_NoNu = 1.5 * EigValues_NoNu[2];
574 
575  double Obs35 = (edm::isNotFinite(Sphericity_NoNu) ? 0 : Sphericity_NoNu);
576  evtselectVarVal.push_back(std::pair<unsigned int, double>(35, Obs35));
577  if (DEBUG)
578  std::cout << "------ LR observable 35 " << Obs35 << " calculated ------" << std::endl;
579 
580  double Obs36 = (edm::isNotFinite(Aplanarity_NoNu) ? 0 : Aplanarity_NoNu);
581  evtselectVarVal.push_back(std::pair<unsigned int, double>(36, Obs36));
582  if (DEBUG)
583  std::cout << "------ LR observable 36 " << Obs36 << " calculated ------" << std::endl;
584 
585  // Sphericity and Aplanarity without boosting back the system to CM frame and without neutrino or lepton
586 
587  TMatrixDSym Matrix_NoNuNoLep(3);
588 
589  double PX2_NoNuNoLep =
590  std::pow(Hadp->Px(), 2) + std::pow(Hadq->Px(), 2) + std::pow(Hadb->Px(), 2) + std::pow(Lepb->Px(), 2);
591  double PY2_NoNuNoLep =
592  std::pow(Hadp->Py(), 2) + std::pow(Hadq->Py(), 2) + std::pow(Hadb->Py(), 2) + std::pow(Lepb->Py(), 2);
593  double PZ2_NoNuNoLep =
594  std::pow(Hadp->Pz(), 2) + std::pow(Hadq->Pz(), 2) + std::pow(Hadb->Pz(), 2) + std::pow(Lepb->Pz(), 2);
595 
596  double P2_NoNuNoLep = PX2_NoNuNoLep + PY2_NoNuNoLep + PZ2_NoNuNoLep;
597 
598  double PXY_NoNuNoLep =
599  Hadp->Px() * Hadp->Py() + Hadq->Px() * Hadq->Py() + Hadb->Px() * Hadb->Py() + Lepb->Px() * Lepb->Py();
600  double PXZ_NoNuNoLep =
601  Hadp->Px() * Hadp->Pz() + Hadq->Px() * Hadq->Pz() + Hadb->Px() * Hadb->Pz() + Lepb->Px() * Lepb->Pz();
602  double PYZ_NoNuNoLep =
603  Hadp->Py() * Hadp->Pz() + Hadq->Py() * Hadq->Pz() + Hadb->Py() * Hadb->Pz() + Lepb->Py() * Lepb->Pz();
604 
605  Matrix_NoNuNoLep(0, 0) = PX2_NoNuNoLep / P2_NoNuNoLep;
606  Matrix_NoNuNoLep(0, 1) = PXY_NoNuNoLep / P2_NoNuNoLep;
607  Matrix_NoNuNoLep(0, 2) = PXZ_NoNuNoLep / P2_NoNuNoLep;
608  Matrix_NoNuNoLep(1, 0) = PXY_NoNuNoLep / P2_NoNuNoLep;
609  Matrix_NoNuNoLep(1, 1) = PY2_NoNuNoLep / P2_NoNuNoLep;
610  Matrix_NoNuNoLep(1, 2) = PYZ_NoNuNoLep / P2_NoNuNoLep;
611  Matrix_NoNuNoLep(2, 0) = PXZ_NoNuNoLep / P2_NoNuNoLep;
612  Matrix_NoNuNoLep(2, 1) = PYZ_NoNuNoLep / P2_NoNuNoLep;
613  Matrix_NoNuNoLep(2, 2) = PZ2_NoNuNoLep / P2_NoNuNoLep;
614 
615  TMatrixDSymEigen pTensor_NoNuNoLep(Matrix_NoNuNoLep);
616 
617  std::vector<double> EigValues_NoNuNoLep;
618  EigValues_NoNuNoLep.clear();
619  for (int i = 0; i < 3; i++) {
620  EigValues_NoNuNoLep.push_back(pTensor_NoNuNoLep.GetEigenValues()[i]);
621  }
622 
623  std::sort(EigValues_NoNuNoLep.begin(), EigValues_NoNuNoLep.end(), dComparator);
624 
625  double Sphericity_NoNuNoLep = 1.5 * (EigValues_NoNuNoLep[1] + EigValues_NoNuNoLep[2]);
626  double Aplanarity_NoNuNoLep = 1.5 * EigValues_NoNuNoLep[2];
627 
628  double Obs37 = (edm::isNotFinite(Sphericity_NoNuNoLep) ? 0 : Sphericity_NoNuNoLep);
629  evtselectVarVal.push_back(std::pair<unsigned int, double>(37, Obs37));
630  if (DEBUG)
631  std::cout << "------ LR observable 37 " << Obs37 << " calculated ------" << std::endl;
632 
633  double Obs38 = (edm::isNotFinite(Aplanarity_NoNuNoLep) ? 0 : Aplanarity_NoNuNoLep);
634  evtselectVarVal.push_back(std::pair<unsigned int, double>(38, Obs38));
635  if (DEBUG)
636  std::cout << "------ LR observable 38 " << Obs38 << " calculated ------" << std::endl;
637 
638  // Put the vector in the TtSemiEvtSolution
639  TS.setLRSignalEvtObservables(evtselectVarVal);
640  if (DEBUG)
641  std::cout << "------ Observable values stored in the event ------" << std::endl;
642 
643  delete Hadp;
644  if (DEBUG)
645  std::cout << "------ Pointer to Hadp deleted ------" << std::endl;
646  delete Hadq;
647  if (DEBUG)
648  std::cout << "------ Pointer to Hadq deleted ------" << std::endl;
649  delete Hadb;
650  if (DEBUG)
651  std::cout << "------ Pointer to Hadb deleted ------" << std::endl;
652  delete Lepb;
653  if (DEBUG)
654  std::cout << "------ Pointer to Lepb deleted ------" << std::endl;
655  delete Lept;
656  if (DEBUG)
657  std::cout << "------ Pointer to Lepn deleted ------" << std::endl;
658  delete Lepn;
659  if (DEBUG)
660  std::cout << "------ Pointer to Mez deleted ------" << std::endl;
661  delete Mez;
662 
663  if (DEBUG)
664  std::cout << "------------ LR observables calculated -----------" << std::endl;
665 }
MEzCalculator::SetLepton
void SetLepton(const pat::Particle &lepton, bool isMuon=true)
Set lepton.
Definition: MEzCalculator.h:31
TtSemiLRSignalSelObservables::operator()
void operator()(TtSemiEvtSolution &, const std::vector< pat::Jet > &)
Definition: TtSemiLRSignalSelObservables.cc:8
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
reco::LeafCandidate::pz
double pz() const final
z coordinate of momentum vector
Definition: LeafCandidate.h:144
mps_fire.i
i
Definition: mps_fire.py:428
DiDispStaMuonMonitor_cfi.pt
tuple pt
Definition: DiDispStaMuonMonitor_cfi.py:39
TtSemiEvtSolution::setLRSignalEvtObservables
void setLRSignalEvtObservables(const std::vector< std::pair< unsigned int, double > > &varval)
Definition: TtSemiEvtSolution.cc:211
edm::isNotFinite
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TtSemiEvtSolution::getDecay
std::string getDecay() const
Definition: TtSemiEvtSolution.h:165
TtSemiEvtSolution::getMuon
pat::Muon getMuon() const
Definition: TtSemiEvtSolution.h:51
hitfit::Matrix
CLHEP::HepMatrix Matrix
Definition: matutil.h:62
WZElectronSkims53X_cff.MT
string MT
Definition: WZElectronSkims53X_cff.py:219
MEzCalculator::SetMET
void SetMET(const pat::MET &MET)
Set MET.
Definition: MEzCalculator.h:25
reco::LeafCandidate::px
double px() const final
x coordinate of momentum vector
Definition: LeafCandidate.h:140
MEzCalculator::Calculate
double Calculate(int type=1)
member functions
Definition: MEzCalculator.cc:14
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TtSemiEvtSolution::getElectron
pat::Electron getElectron() const
Definition: TtSemiEvtSolution.h:52
reco::LeafCandidate::py
double py() const final
y coordinate of momentum vector
Definition: LeafCandidate.h:142
PI
#define PI
Definition: QcdUeDQM.h:37
Thrust
Definition: Thrust.h:38
N
#define N
Definition: blowfish.cc:9
funct::D
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
TtSemiEvtSolution::getNeutrino
pat::MET getNeutrino() const
Definition: TtSemiEvtSolution.h:53
reco::LeafCandidate::et
double et() const final
transverse energy
Definition: LeafCandidate.h:127
DEBUG
#define DEBUG
Definition: DMRChecker.cc:120
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
HT
Definition: HT.h:21
TtSemiLRSignalSelObservables::evtselectVarVal
std::vector< std::pair< unsigned int, double > > evtselectVarVal
Definition: TtSemiLRSignalSelObservables.h:37
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
reco::LeafCandidate::energy
double energy() const final
energy
Definition: LeafCandidate.h:125
Generated for CMSSW Reference Manual by doxygen 1.8.5