CMS 3D CMS Logo
Megajet.cc
Go to the documentation of this file.
1 #include "PhysicsTools/Heppy/interface/Megajet.h"
2 
3 #include <vector>
4 #include <cmath>
5 #include <TLorentzVector.h>
6 
7 using namespace std;
8 
9 namespace heppy {
10 
11  // constructor specifying the association methods
12  Megajet::Megajet(vector<float> Px_vector,
13  vector<float> Py_vector,
14  vector<float> Pz_vector,
15  vector<float> E_vector,
16  int megajet_association_method)
17  : Object_Px(Px_vector),
18  Object_Py(Py_vector),
19  Object_Pz(Pz_vector),
20  Object_E(E_vector),
21  megajet_meth(megajet_association_method),
22  status(0) {
23  if (Object_Px.size() < 2)
24  cout << "Error in Megajet: you should provide at least two jets to form Megajets" << endl;
25  for (int j = 0; j < (int)jIN.size(); ++j) {
26  jIN.push_back(TLorentzVector(Object_Px[j], Object_Py[j], Object_Pz[j], Object_E[j]));
27  }
28  }
29 
30  // constructor without specification of the seed and association methods
31  // in this case, the latter must be given by calling SetMethod before invoking Combine()
32  Megajet::Megajet(vector<float> Px_vector, vector<float> Py_vector, vector<float> Pz_vector, vector<float> E_vector)
33  : Object_Px(Px_vector), Object_Py(Py_vector), Object_Pz(Pz_vector), Object_E(E_vector), status(0) {
34  if (Object_Px.size() < 2)
35  cout << "Error in Megajet: you should provide at least two jets to form Megajets" << endl;
36  for (int j = 0; j < (int)jIN.size(); ++j) {
37  jIN.push_back(TLorentzVector(Object_Px[j], Object_Py[j], Object_Pz[j], Object_E[j]));
38  }
39  }
40 
41  vector<float> Megajet::getAxis1() {
42  if (status != 1) {
43  this->Combine();
44  if (megajet_meth == 1)
45  this->CombineMinMass();
46  else if (megajet_meth == 2)
47  this->CombineMinHT();
48  else if (megajet_meth == 3)
49  this->CombineMinEnergyMass();
50  else if (megajet_meth == 4)
51  this->CombineGeorgi();
52  else
53  this->CombineMinMass();
54  }
55  if (jIN.size() > 1) {
56  Axis1[0] = jOUT[0].Px() / jOUT[0].P();
57  Axis1[1] = jOUT[0].Py() / jOUT[0].P();
58  Axis1[2] = jOUT[0].Pz() / jOUT[0].P();
59  Axis1[3] = jOUT[0].P();
60  Axis1[4] = jOUT[0].E();
61  }
62  return Axis1;
63  }
64  vector<float> Megajet::getAxis2() {
65  if (status != 1) {
66  this->Combine();
67  if (megajet_meth == 1)
68  this->CombineMinMass();
69  else if (megajet_meth == 2)
70  this->CombineMinHT();
71  else if (megajet_meth == 3)
72  this->CombineMinEnergyMass();
73  else if (megajet_meth == 4)
74  this->CombineGeorgi();
75  else
76  this->CombineMinMass();
77  }
78  if (jIN.size() > 1) {
79  Axis2[0] = jOUT[1].Px() / jOUT[1].P();
80  Axis2[1] = jOUT[1].Py() / jOUT[1].P();
81  Axis2[2] = jOUT[1].Pz() / jOUT[1].P();
82  Axis2[3] = jOUT[1].P();
83  Axis2[4] = jOUT[1].E();
84  }
85  return Axis2;
86  }
87 
88  void Megajet::Combine() {
89  int N_JETS = (int)jIN.size();
90 
91  int N_comb = 1;
92  for (int i = 0; i < N_JETS; i++) {
93  N_comb *= 2;
94  }
95 
96  // clear some vectors if method Combine() is called again
97  if (!j1.empty()) {
98  j1.clear();
99  j2.clear();
100  Axis1.clear();
101  Axis2.clear();
102  }
103 
104  for (int j = 0; j < 5; ++j) {
105  Axis1.push_back(0);
106  Axis2.push_back(0);
107  }
108 
109  int j_count;
110  for (int i = 1; i < N_comb - 1; i++) {
111  TLorentzVector j_temp1, j_temp2;
112  int itemp = i;
113  j_count = N_comb / 2;
114  int count = 0;
115  while (j_count > 0) {
116  if (itemp / j_count == 1) {
117  j_temp1 += jIN[count];
118  } else {
119  j_temp2 += jIN[count];
120  }
121  itemp -= j_count * (itemp / j_count);
122  j_count /= 2;
123  count++;
124  }
125 
126  j1.push_back(j_temp1);
127  j2.push_back(j_temp2);
128  }
129  }
130 
131  void Megajet::CombineMinMass() {
132  double M_min = -1;
133  // default value (in case none is found)
134  TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
135  TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
136  for (int i = 0; i < (int)j1.size(); i++) {
137  double M_temp = j1[i].M2() + j2[i].M2();
138  if (M_min < 0 || M_temp < M_min) {
139  M_min = M_temp;
140  myJ1 = j1[i];
141  myJ2 = j2[i];
142  }
143  }
144  // myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
145  // myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
146 
147  jOUT.clear();
148  if (myJ1.Pt() > myJ2.Pt()) {
149  jOUT.push_back(myJ1);
150  jOUT.push_back(myJ2);
151  } else {
152  jOUT.push_back(myJ2);
153  jOUT.push_back(myJ1);
154  }
155  status = 1;
156  }
157 
158  void Megajet::CombineMinEnergyMass() {
159  double M_min = -1;
160  // default value (in case none is found)
161  TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
162  TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
163  for (int i = 0; i < (int)j1.size(); i++) {
164  double M_temp = j1[i].M2() / j1[i].E() + j2[i].M2() / j2[i].E();
165  if (M_min < 0 || M_temp < M_min) {
166  M_min = M_temp;
167  myJ1 = j1[i];
168  myJ2 = j2[i];
169  }
170  }
171 
172  // myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
173  // myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
174 
175  jOUT.clear();
176  if (myJ1.Pt() > myJ2.Pt()) {
177  jOUT.push_back(myJ1);
178  jOUT.push_back(myJ2);
179  } else {
180  jOUT.push_back(myJ2);
181  jOUT.push_back(myJ1);
182  }
183  status = 1;
184  }
185 
186  void Megajet::CombineGeorgi() {
187  double M_max = -10000;
188  // default value (in case none is found)
189  TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
190  TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
191  for (int i = 0; i < (int)j1.size(); i++) {
192  int myBeta = 2;
193  double M_temp = (j1[i].E() - myBeta * j1[i].M2() / j1[i].E()) + (j2[i].E() - myBeta * j2[i].M2() / j2[i].E());
194  if (M_max < -9999 || M_temp > M_max) {
195  M_max = M_temp;
196  myJ1 = j1[i];
197  myJ2 = j2[i];
198  }
199  }
200 
201  // myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
202  // myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
203 
204  jOUT.clear();
205  if (myJ1.Pt() > myJ2.Pt()) {
206  jOUT.push_back(myJ1);
207  jOUT.push_back(myJ2);
208  } else {
209  jOUT.push_back(myJ2);
210  jOUT.push_back(myJ1);
211  }
212  status = 1;
213  }
214 
215  void Megajet::CombineMinHT() {
216  double dHT_min = 999999999999999.0;
217  // default value (in case none is found)
218  TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
219  TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
220  for (int i = 0; i < (int)j1.size(); i++) {
221  double dHT_temp = fabs(j1[i].E() - j2[i].E());
222  if (dHT_temp < dHT_min) {
223  dHT_min = dHT_temp;
224  myJ1 = j1[i];
225  myJ2 = j2[i];
226  }
227  }
228 
229  jOUT.clear();
230  if (myJ1.Pt() > myJ2.Pt()) {
231  jOUT.push_back(myJ1);
232  jOUT.push_back(myJ2);
233  } else {
234  jOUT.push_back(myJ2);
235  jOUT.push_back(myJ1);
236  }
237  status = 1;
238  }
239 
240 } // namespace heppy
heppy::Megajet::j1
std::vector< TLorentzVector > j1
Definition: Megajet.h:96
mps_fire.i
i
Definition: mps_fire.py:428
heppy::Megajet::jOUT
std::vector< TLorentzVector > jOUT
Definition: Megajet.h:95
Megajet.h
mps_update.status
status
Definition: mps_update.py:69
gather_cfg.cout
cout
Definition: gather_cfg.py:144
heppy::Megajet::Object_Px
std::vector< float > Object_Px
Definition: Megajet.h:83
heppy::Megajet::CombineGeorgi
void CombineGeorgi()
Combining the jets in two hemispheres by maximizing (E1-Beta*m1^2/E1 + E2-Beta*m1^2/E2)
Definition: Megajet.cc:186
heppy::Megajet::jIN
std::vector< TLorentzVector > jIN
Definition: Megajet.h:94
heppy::Megajet::getAxis2
std::vector< float > getAxis2()
Definition: Megajet.cc:64
heppy::Megajet::Megajet
Megajet()
Definition: Megajet.h:35
submitPVResolutionJobs.count
count
Definition: submitPVResolutionJobs.py:352
heppy::Megajet::Combine
void Combine()
Combine the jets in all the possible pairs of hemispheres.
Definition: Megajet.cc:88
heppy::Megajet::Axis2
std::vector< float > Axis2
Definition: Megajet.h:89
heppy::Megajet::Axis1
std::vector< float > Axis1
Definition: Megajet.h:88
heppy
TAKEN FROM http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/ElectroWeakAnalysis/Utilities/src/PdfWeig...
Definition: AlphaT.h:16
heppy::Megajet::CombineMinEnergyMass
void CombineMinEnergyMass()
Combining the jets in two hemispheres by minimizing m1^2/E1 + m2^2/E2.
Definition: Megajet.cc:158
heppy::Megajet::Object_Pz
std::vector< float > Object_Pz
Definition: Megajet.h:85
createfilelist.int
int
Definition: createfilelist.py:10
heppy::Megajet::j2
std::vector< TLorentzVector > j2
Definition: Megajet.h:97
heppy::Megajet::Object_E
std::vector< float > Object_E
Definition: Megajet.h:86
std
Definition: JetResolutionObject.h:76
heppy::Megajet::getAxis1
std::vector< float > getAxis1()
Definition: Megajet.cc:41
heppy::Megajet::CombineMinHT
void CombineMinHT()
Definition: Megajet.cc:215
heppy::Megajet::CombineMinMass
void CombineMinMass()
Definition: Megajet.cc:131
heppy::Megajet::Object_Py
std::vector< float > Object_Py
Definition: Megajet.h:84
heppy::Megajet::megajet_meth
int megajet_meth
Definition: Megajet.h:91
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
heppy::Megajet::status
int status
Definition: Megajet.h:92
Generated for CMSSW Reference Manual by doxygen 1.8.16