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