CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TtHadLRSignalSelObservables.cc
Go to the documentation of this file.
1 #include "TopQuarkAnalysis/TopEventSelection/interface/TtHadLRSignalSelObservables.h"
2 
3 TtHadLRSignalSelObservables::TtHadLRSignalSelObservables()
4 {
5 }
6 
7 TtHadLRSignalSelObservables::~TtHadLRSignalSelObservables()
8 {
9 }
10 
11 void TtHadLRSignalSelObservables::operator() (TtHadEvtSolution &TS)
12 {
13  evtselectVarVal.clear();
14  std::vector<pat::Jet> topJets;
15  topJets.clear();
16  topJets.push_back(TS.getHadp());
17  topJets.push_back(TS.getHadq());
18  topJets.push_back(TS.getHadj());
19  topJets.push_back(TS.getHadk());
20  topJets.push_back(TS.getHadb());
21  topJets.push_back(TS.getHadbbar());
22 
23  //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...
24  TLorentzVector *Hadp = new TLorentzVector();
25  TLorentzVector *Hadq = new TLorentzVector();
26  TLorentzVector *Hadj = new TLorentzVector();
27  TLorentzVector *Hadk = new TLorentzVector();
28  TLorentzVector *Hadb = new TLorentzVector();
29  TLorentzVector *Hadbbar = new TLorentzVector();
30  Hadp->SetPxPyPzE(topJets[0].px(),topJets[0].py(),topJets[0].pz(),topJets[3].energy());
31  Hadq->SetPxPyPzE(topJets[1].px(),topJets[1].py(),topJets[1].pz(),topJets[2].energy());
32  Hadj->SetPxPyPzE(topJets[2].px(),topJets[2].py(),topJets[2].pz(),topJets[2].energy());
33  Hadk->SetPxPyPzE(topJets[3].px(),topJets[3].py(),topJets[3].pz(),topJets[3].energy());
34  Hadb->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy());
35  Hadbbar->SetPxPyPzE(topJets[4].px(),topJets[4].py(),topJets[4].pz(),topJets[1].energy());
36 
37  //sort the topJets in Et
38  std::sort(topJets.begin(),topJets.end(),EtComparator);
39 
40  //Et-Sum of the lightest jets
41  double EtSum = topJets[5].et()+topJets[5].et();
42  double Obs1 = (EtSum>0 ? EtSum : -1);
43  evtselectVarVal.push_back(std::pair<unsigned int,double>(1,Obs1));
44 
45  //Log of the difference in Bdisc between the 2nd and the 3rd jets (ordered in Bdisc) - does that still
46  //make sense for the fully hadronic channel as well?FIXME!!!
47 
48  //sort the jets in Bdiscriminant
49  std::sort(topJets.begin(),topJets.end(),BdiscComparator);
50 
51  double BGap = topJets[1].bDiscriminator("trackCountingJetTags") - topJets[2].bDiscriminator("trackCountingJetTags");
52  double Obs2 = (BGap>0 ? log(BGap) : -1);
53  evtselectVarVal.push_back(std::pair<unsigned int,double>(2,Obs2));
54 
55  //Circularity of the event
56  double N=0,D=0,C_tmp=0,C=1000;
57  double nx,ny,nz;
58 
59  // 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))
60 
61  for(unsigned int i=0;i<6;i++){
62  D += fabs(topJets[i].pt());
63  }
64 
65  if((D>0)){
66  // Loop over all the unit vectors in the transverse plane in order to find the miminum :
67  for(unsigned int i=0; i<360; i++){
68  nx = cos((2*PI/360)*i);
69  ny = sin((2*PI/360)*i);
70  nz = 0;
71  N=0;
72  for(unsigned int i=0;i<4;i++){
73  N += fabs(topJets[i].px()*nx+topJets[i].py()*ny+topJets[i].pz()*nz);
74  }
75  C_tmp = 2*N/D;
76  if(C_tmp<C) C = C_tmp;
77  }
78  }
79 
80  double Obs3 = ( C!=1000 ? C : -1);
81  evtselectVarVal.push_back(std::pair<unsigned int,double>(3,Obs3));
82 
83 
84  //HT variable (Et-sum of the six jets)
85  double HT=0;
86  for(unsigned int i=0;i<6;i++){
87  HT += topJets[i].et();
88  }
89 
90  double Obs4 = ( HT!=0 ? HT : -1);
91  evtselectVarVal.push_back(std::pair<unsigned int,double>(4,Obs4));
92 
93  //Transverse Mass of the system
94  math::XYZTLorentzVector pjets;
95  // for the six jets
96  for(unsigned int i=0;i<6;i++){
97  pjets += topJets[i].p4();
98  }
99  double MT = sqrt(std::pow(pjets.mass(),2));
100  double Obs5 = ( MT>0 ? MT : -1);
101  evtselectVarVal.push_back(std::pair<unsigned int,double>(5,Obs5));
102 
103  //CosTheta(Hadp,Hadq) and CosTheta(Hadj,Hadq)
104  //sort the lightJets in Et
105  std::sort(topJets.begin(),topJets.end(),EtComparator);
106 
107  double px1 = topJets[2].px(); double px2 = topJets[3].px();
108  double py1 = topJets[2].py(); double py2 = topJets[3].py();
109  double pz1 = topJets[2].pz(); double pz2 = topJets[3].pz();
110  double E1 = topJets[2].energy(); double E2 = topJets[3].energy();
111  double px3 = topJets[4].px(); double px4 = topJets[5].px();
112  double py3 = topJets[4].py(); double py4 = topJets[5].py();
113  double pz3 = topJets[4].pz(); double pz4 = topJets[5].pz();
114  double E3 = topJets[4].energy(); double E4 = topJets[5].energy();
115 
116  TLorentzVector *LightJet1 = new TLorentzVector();
117  TLorentzVector *LightJet2 = new TLorentzVector();
118  TLorentzVector *LightJet3 = new TLorentzVector();
119  TLorentzVector *LightJet4 = new TLorentzVector();
120 
121  LightJet1->SetPxPyPzE(px1,py1,pz1,E1);
122  //LightJet1->Boost(BoostBackToCM);
123  LightJet2->SetPxPyPzE(px2,py2,pz2,E2);
124  //LightJet2->Boost(BoostBackToCM);
125  LightJet3->SetPxPyPzE(px3,py3,pz3,E3);
126  //LightJet1->Boost(BoostBackToCM);
127  LightJet4->SetPxPyPzE(px4,py4,pz4,E4);
128  //LightJet2->Boost(BoostBackToCM);
129 
130  double CosTheta1 = cos(LightJet2->Angle(LightJet1->Vect()));
131  double CosTheta2 = cos(LightJet4->Angle(LightJet3->Vect()));
132 
133  double Obs6 = ( -1<CosTheta1 ? CosTheta1 : -2);
134  evtselectVarVal.push_back(std::pair<unsigned int,double>(6,Obs6));
135  double Obs7 = ( -1<CosTheta2 ? CosTheta2 : -2);
136  evtselectVarVal.push_back(std::pair<unsigned int,double>(7,Obs7));
137 
138  delete LightJet1;
139  delete LightJet2;
140  delete LightJet3;
141  delete LightJet4;
142 
143  // try to find out more powerful observables related to the b-disc
144  //sort the jets in Bdiscriminant
145  std::sort(topJets.begin(),topJets.end(),BdiscComparator);
146 
147  double BjetsBdiscSum = topJets[0].bDiscriminator("trackCountingJetTags") + topJets[1].bDiscriminator("trackCountingJetTags");
148  double Ljets1BdiscSum = topJets[2].bDiscriminator("trackCountingJetTags") + topJets[3].bDiscriminator("trackCountingJetTags");
149  double Ljets2BdiscSum = topJets[4].bDiscriminator("trackCountingJetTags") + topJets[5].bDiscriminator("trackCountingJetTags");
150  //std::cout<<"BjetsBdiscSum = "<<BjetsBdiscSum<<std::endl;
151  //std::cout<<"LjetsBdiscSum = "<<LjetsBdiscSum<<std::endl;
152 
153  double Obs8 = (Ljets1BdiscSum !=0 ? log(BjetsBdiscSum/Ljets1BdiscSum) : -1);
154  evtselectVarVal.push_back(std::pair<unsigned int,double>(8,Obs8));
155  double Obs9 = (Ljets2BdiscSum !=0 ? log(BjetsBdiscSum/Ljets2BdiscSum) : -1);
156  evtselectVarVal.push_back(std::pair<unsigned int,double>(9,Obs9));
157  double Obs10 = (BGap>0 ? log(BjetsBdiscSum*BGap) : -1);
158  evtselectVarVal.push_back(std::pair<unsigned int,double>(10,Obs10));
159 
160 
161  // Et-Ratio between the two highest in Et jets and four highest jets
162  double Obs11 = (HT!=0 ? (HT-EtSum)/HT : -1);
163  evtselectVarVal.push_back(std::pair<unsigned int,double>(11,Obs11));
164 
165 
166  //Sphericity and Aplanarity without boosting back the system to CM frame
167 
168  TMatrixDSym Matrix(3);
169  double PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+
170  std::pow(Hadj->Px(),2)+std::pow(Hadk->Px(),2)+std::pow(Hadbbar->Px(),2);
171  double PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+
172  std::pow(Hadj->Py(),2)+std::pow(Hadk->Py(),2)+std::pow(Hadbbar->Py(),2);
173  double PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+
174  std::pow(Hadj->Pz(),2)+std::pow(Hadk->Pz(),2)+std::pow(Hadbbar->Pz(),2);
175  double P2 = PX2+PY2+PZ2;
176 
177  double PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
178  Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
179  double PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
180  Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
181  double PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
182  Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
183 
184  Matrix(0,0) = PX2/P2; Matrix(0,1) = PXY/P2; Matrix(0,2) = PXZ/P2;
185  Matrix(1,0) = PXY/P2; Matrix(1,1) = PY2/P2; Matrix(1,2) = PYZ/P2;
186  Matrix(2,0) = PXZ/P2; Matrix(2,1) = PYZ/P2; Matrix(2,2) = PZ2/P2;
187 
188  TMatrixDSymEigen pTensor(Matrix);
189 
190  std::vector<double> EigValues;
191  EigValues.clear();
192  for(int i=0;i<3;i++){
193  EigValues.push_back(pTensor.GetEigenValues()[i]);
194  }
195 
196  std::sort(EigValues.begin(),EigValues.end(),dComparator);
197 
198  double Sphericity = 1.5*(EigValues[1]+EigValues[2]);
199  double Aplanarity = 1.5*EigValues[2];
200 
201  double Obs12 = (isnan(Sphericity) ? -1 : Sphericity);
202  evtselectVarVal.push_back(std::pair<unsigned int,double>(12,Obs12));
203 
204  double Obs13 = (isnan(Aplanarity) ? -1 : Aplanarity);
205  evtselectVarVal.push_back(std::pair<unsigned int,double>(13,Obs13));
206 
207 
208  //Sphericity and Aplanarity with boosting back the system to CM frame
209  TLorentzVector *TtbarSystem = new TLorentzVector();
210  TtbarSystem->SetPx(Hadp->Px()+Hadq->Px()+Hadb->Px()+Hadj->Px()+Hadk->Px()+Hadbbar->Px());
211  TtbarSystem->SetPy(Hadp->Py()+Hadq->Py()+Hadb->Py()+Hadj->Py()+Hadk->Py()+Hadbbar->Py());
212  TtbarSystem->SetPz(Hadp->Pz()+Hadq->Pz()+Hadb->Pz()+Hadj->Pz()+Hadk->Pz()+Hadbbar->Pz());
213  TtbarSystem->SetE(Hadp->E()+Hadq->E()+Hadb->E()+Hadj->E()+Hadk->E()+Hadbbar->E());
214 
215  TVector3 BoostBackToCM = -(TtbarSystem->BoostVector());
216  Hadp->Boost(BoostBackToCM);
217  Hadq->Boost(BoostBackToCM);
218  Hadb->Boost(BoostBackToCM);
219  Hadj->Boost(BoostBackToCM);
220  Hadk->Boost(BoostBackToCM);
221  Hadbbar->Boost(BoostBackToCM);
222 
223  double BOOST_PX2 = std::pow(Hadp->Px(),2)+std::pow(Hadq->Px(),2)+std::pow(Hadb->Px(),2)+
224  std::pow(Hadj->Px(),2)+std::pow(Hadk->Px(),2)+std::pow(Hadbbar->Px(),2);
225  double BOOST_PY2 = std::pow(Hadp->Py(),2)+std::pow(Hadq->Py(),2)+std::pow(Hadb->Py(),2)+
226  std::pow(Hadj->Py(),2)+std::pow(Hadk->Py(),2)+std::pow(Hadbbar->Py(),2);
227  double BOOST_PZ2 = std::pow(Hadp->Pz(),2)+std::pow(Hadq->Pz(),2)+std::pow(Hadb->Pz(),2)+
228  std::pow(Hadj->Pz(),2)+std::pow(Hadk->Pz(),2)+std::pow(Hadbbar->Pz(),2);
229 
230  double BOOST_P2 = BOOST_PX2+BOOST_PY2+BOOST_PZ2;
231 
232  double BOOST_PXY = Hadp->Px()*Hadp->Py()+Hadq->Px()*Hadq->Py()+Hadb->Px()*Hadb->Py()+
233  Hadj->Px()*Hadj->Py()+Hadk->Px()*Hadk->Py()+Hadbbar->Px()*Hadbbar->Py();
234  double BOOST_PXZ = Hadp->Px()*Hadp->Pz()+Hadq->Px()*Hadq->Pz()+Hadb->Px()*Hadb->Pz()+
235  Hadj->Px()*Hadj->Pz()+Hadk->Px()*Hadk->Pz()+Hadbbar->Px()*Hadbbar->Pz();
236  double BOOST_PYZ = Hadp->Py()*Hadp->Pz()+Hadq->Py()*Hadq->Pz()+Hadb->Py()*Hadb->Pz()+
237  Hadj->Py()*Hadj->Pz()+Hadk->Py()*Hadk->Pz()+Hadbbar->Py()*Hadbbar->Pz();
238 
239  TMatrixDSym BOOST_Matrix(3);
240 
241  BOOST_Matrix(0,0) = BOOST_PX2/BOOST_P2; BOOST_Matrix(0,1) = BOOST_PXY/BOOST_P2; BOOST_Matrix(0,2) = BOOST_PXZ/BOOST_P2;
242  BOOST_Matrix(1,0) = BOOST_PXY/BOOST_P2; BOOST_Matrix(1,1) = BOOST_PY2/BOOST_P2; BOOST_Matrix(1,2) = BOOST_PYZ/BOOST_P2;
243  BOOST_Matrix(2,0) = BOOST_PXZ/BOOST_P2; BOOST_Matrix(2,1) = BOOST_PYZ/BOOST_P2; BOOST_Matrix(2,2) = BOOST_PZ2/BOOST_P2;
244 
245  TMatrixDSymEigen BOOST_pTensor(BOOST_Matrix);
246 
247  std::vector<double> BOOST_EigValues;
248  BOOST_EigValues.clear();
249  for(int i=0;i<3;i++){
250  BOOST_EigValues.push_back(BOOST_pTensor.GetEigenValues()[i]);
251  }
252 
253  std::sort(BOOST_EigValues.begin(),BOOST_EigValues.end(),dComparator);
254 
255  double BOOST_Sphericity = 1.5*(BOOST_EigValues[1]+BOOST_EigValues[2]);
256  double BOOST_Aplanarity = 1.5*BOOST_EigValues[2];
257 
258  double Obs14 = ( isnan(BOOST_Sphericity) ? -1 : BOOST_Sphericity );
259  evtselectVarVal.push_back(std::pair<unsigned int,double>(14,Obs14));
260 
261  double Obs15 = ( isnan(BOOST_Aplanarity) ? -1 : BOOST_Aplanarity );
262  evtselectVarVal.push_back(std::pair<unsigned int,double>(15,Obs15));
263 
264 
265  // Centrality of the six jets
266  double H=0;
267  for(unsigned int i=0;i<6;i++){
268  H += topJets[i].energy();
269  }
270  double Obs16 = ( H != 0 ? HT/H : -1 );
271  evtselectVarVal.push_back(std::pair<unsigned int,double>(16,Obs16));
272 
273 
274  // distance from the origin in the (BjetsBdiscSum, Ljets1BdiscSum) and (BjetsBdiscSum, Ljets2BdiscSum)
275  double Obs17 = ( BjetsBdiscSum != 0 && Ljets1BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets1BdiscSum) : -1 );
276  evtselectVarVal.push_back(std::pair<unsigned int,double>(17,Obs17));
277  double Obs18 = ( BjetsBdiscSum != 0 && Ljets2BdiscSum != 0 ? 0.707*(BjetsBdiscSum-Ljets2BdiscSum) : -1 );
278  evtselectVarVal.push_back(std::pair<unsigned int,double>(18,Obs18));
279 
280  // Ratio of the Et Sum of the two lightest jets over HT=Sum of the Et of the six highest jets in Et
281  double Obs19 = ( HT != 0 ? EtSum/HT : -1 );
282  evtselectVarVal.push_back(std::pair<unsigned int,double>(19,Obs19));
283 
284 
285  // Put the vector in the TtHadEvtSolution
286  TS.setLRSignalEvtObservables(evtselectVarVal);
287 }
i
int i
Definition: DBlmapReader.cc:9
WElectronSkim_cff.MT
string MT
|_ _( _ ) | _ \ | _ \ __ _(_)_ __ ___ | | / _ \/\ |_) | | |_) / _` | | &#39;__/ __| | || (_&gt; &lt; __/ | __/...
Definition: WElectronSkim_cff.py:135
PI
#define PI
Definition: HcalBeamMonitor.cc:14
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TtHadEvtSolution::getHadbbar
pat::Jet getHadbbar() const
Definition: TtHadEvtSolution.cc:86
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1108
TtHadEvtSolution::getHadb
pat::Jet getHadb() const
Definition: TtHadEvtSolution.cc:38
TtHadEvtSolution::getHadq
pat::Jet getHadq() const
Definition: TtHadEvtSolution.cc:70
TtHadEvtSolution::getHadj
pat::Jet getHadj() const
Definition: TtHadEvtSolution.cc:102
hitfit::Matrix
CLHEP::HepMatrix Matrix
Definition: matutil.h:66
math::XYZTLorentzVector
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:30
TtHadEvtSolution::setLRSignalEvtObservables
void setLRSignalEvtObservables(std::vector< std::pair< unsigned int, double > > varval)
Definition: TtHadEvtSolution.cc:263
edm::detail::isnan
bool isnan(float x)
Definition: math.h:13
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
TtHadEvtSolution::getHadp
pat::Jet getHadp() const
Definition: TtHadEvtSolution.cc:54
TtHadEvtSolution::getHadk
pat::Jet getHadk() const
Definition: TtHadEvtSolution.cc:118
N
#define N
Definition: blowfish.cc:9
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
TtHadLRSignalSelObservables::evtselectVarVal
std::vector< std::pair< unsigned int, double > > evtselectVarVal
Definition: TtHadLRSignalSelObservables.h:67
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