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