File indexing completed on 2024-04-06 12:23:29
0001 #include "PhysicsTools/Heppy/interface/Megajet.h"
0002
0003 #include <vector>
0004 #include <cmath>
0005 #include <TLorentzVector.h>
0006
0007 using namespace std;
0008
0009 namespace heppy {
0010
0011
0012 Megajet::Megajet(vector<float> Px_vector,
0013 vector<float> Py_vector,
0014 vector<float> Pz_vector,
0015 vector<float> E_vector,
0016 int megajet_association_method)
0017 : Object_Px(Px_vector),
0018 Object_Py(Py_vector),
0019 Object_Pz(Pz_vector),
0020 Object_E(E_vector),
0021 megajet_meth(megajet_association_method),
0022 status(0) {
0023 if (Object_Px.size() < 2)
0024 cout << "Error in Megajet: you should provide at least two jets to form Megajets" << endl;
0025 for (int j = 0; j < (int)jIN.size(); ++j) {
0026 jIN.push_back(TLorentzVector(Object_Px[j], Object_Py[j], Object_Pz[j], Object_E[j]));
0027 }
0028 }
0029
0030
0031
0032 Megajet::Megajet(vector<float> Px_vector, vector<float> Py_vector, vector<float> Pz_vector, vector<float> E_vector)
0033 : Object_Px(Px_vector), Object_Py(Py_vector), Object_Pz(Pz_vector), Object_E(E_vector), status(0) {
0034 if (Object_Px.size() < 2)
0035 cout << "Error in Megajet: you should provide at least two jets to form Megajets" << endl;
0036 for (int j = 0; j < (int)jIN.size(); ++j) {
0037 jIN.push_back(TLorentzVector(Object_Px[j], Object_Py[j], Object_Pz[j], Object_E[j]));
0038 }
0039 }
0040
0041 vector<float> Megajet::getAxis1() {
0042 if (status != 1) {
0043 this->Combine();
0044 if (megajet_meth == 1)
0045 this->CombineMinMass();
0046 else if (megajet_meth == 2)
0047 this->CombineMinHT();
0048 else if (megajet_meth == 3)
0049 this->CombineMinEnergyMass();
0050 else if (megajet_meth == 4)
0051 this->CombineGeorgi();
0052 else
0053 this->CombineMinMass();
0054 }
0055 if (jIN.size() > 1) {
0056 Axis1[0] = jOUT[0].Px() / jOUT[0].P();
0057 Axis1[1] = jOUT[0].Py() / jOUT[0].P();
0058 Axis1[2] = jOUT[0].Pz() / jOUT[0].P();
0059 Axis1[3] = jOUT[0].P();
0060 Axis1[4] = jOUT[0].E();
0061 }
0062 return Axis1;
0063 }
0064 vector<float> Megajet::getAxis2() {
0065 if (status != 1) {
0066 this->Combine();
0067 if (megajet_meth == 1)
0068 this->CombineMinMass();
0069 else if (megajet_meth == 2)
0070 this->CombineMinHT();
0071 else if (megajet_meth == 3)
0072 this->CombineMinEnergyMass();
0073 else if (megajet_meth == 4)
0074 this->CombineGeorgi();
0075 else
0076 this->CombineMinMass();
0077 }
0078 if (jIN.size() > 1) {
0079 Axis2[0] = jOUT[1].Px() / jOUT[1].P();
0080 Axis2[1] = jOUT[1].Py() / jOUT[1].P();
0081 Axis2[2] = jOUT[1].Pz() / jOUT[1].P();
0082 Axis2[3] = jOUT[1].P();
0083 Axis2[4] = jOUT[1].E();
0084 }
0085 return Axis2;
0086 }
0087
0088 void Megajet::Combine() {
0089 int N_JETS = (int)jIN.size();
0090
0091 int N_comb = 1;
0092 for (int i = 0; i < N_JETS; i++) {
0093 N_comb *= 2;
0094 }
0095
0096
0097 if (!j1.empty()) {
0098 j1.clear();
0099 j2.clear();
0100 Axis1.clear();
0101 Axis2.clear();
0102 }
0103
0104 for (int j = 0; j < 5; ++j) {
0105 Axis1.push_back(0);
0106 Axis2.push_back(0);
0107 }
0108
0109 int j_count;
0110 for (int i = 1; i < N_comb - 1; i++) {
0111 TLorentzVector j_temp1, j_temp2;
0112 int itemp = i;
0113 j_count = N_comb / 2;
0114 int count = 0;
0115 while (j_count > 0) {
0116 if (itemp / j_count == 1) {
0117 j_temp1 += jIN[count];
0118 } else {
0119 j_temp2 += jIN[count];
0120 }
0121 itemp -= j_count * (itemp / j_count);
0122 j_count /= 2;
0123 count++;
0124 }
0125
0126 j1.push_back(j_temp1);
0127 j2.push_back(j_temp2);
0128 }
0129 }
0130
0131 void Megajet::CombineMinMass() {
0132 double M_min = -1;
0133
0134 TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
0135 TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
0136 for (int i = 0; i < (int)j1.size(); i++) {
0137 double M_temp = j1[i].M2() + j2[i].M2();
0138 if (M_min < 0 || M_temp < M_min) {
0139 M_min = M_temp;
0140 myJ1 = j1[i];
0141 myJ2 = j2[i];
0142 }
0143 }
0144
0145
0146
0147 jOUT.clear();
0148 if (myJ1.Pt() > myJ2.Pt()) {
0149 jOUT.push_back(myJ1);
0150 jOUT.push_back(myJ2);
0151 } else {
0152 jOUT.push_back(myJ2);
0153 jOUT.push_back(myJ1);
0154 }
0155 status = 1;
0156 }
0157
0158 void Megajet::CombineMinEnergyMass() {
0159 double M_min = -1;
0160
0161 TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
0162 TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
0163 for (int i = 0; i < (int)j1.size(); i++) {
0164 double M_temp = j1[i].M2() / j1[i].E() + j2[i].M2() / j2[i].E();
0165 if (M_min < 0 || M_temp < M_min) {
0166 M_min = M_temp;
0167 myJ1 = j1[i];
0168 myJ2 = j2[i];
0169 }
0170 }
0171
0172
0173
0174
0175 jOUT.clear();
0176 if (myJ1.Pt() > myJ2.Pt()) {
0177 jOUT.push_back(myJ1);
0178 jOUT.push_back(myJ2);
0179 } else {
0180 jOUT.push_back(myJ2);
0181 jOUT.push_back(myJ1);
0182 }
0183 status = 1;
0184 }
0185
0186 void Megajet::CombineGeorgi() {
0187 double M_max = -10000;
0188
0189 TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
0190 TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
0191 for (int i = 0; i < (int)j1.size(); i++) {
0192 int myBeta = 2;
0193 double M_temp = (j1[i].E() - myBeta * j1[i].M2() / j1[i].E()) + (j2[i].E() - myBeta * j2[i].M2() / j2[i].E());
0194 if (M_max < -9999 || M_temp > M_max) {
0195 M_max = M_temp;
0196 myJ1 = j1[i];
0197 myJ2 = j2[i];
0198 }
0199 }
0200
0201
0202
0203
0204 jOUT.clear();
0205 if (myJ1.Pt() > myJ2.Pt()) {
0206 jOUT.push_back(myJ1);
0207 jOUT.push_back(myJ2);
0208 } else {
0209 jOUT.push_back(myJ2);
0210 jOUT.push_back(myJ1);
0211 }
0212 status = 1;
0213 }
0214
0215 void Megajet::CombineMinHT() {
0216 double dHT_min = 999999999999999.0;
0217
0218 TLorentzVector myJ1 = TLorentzVector(0, 0, 0, 0);
0219 TLorentzVector myJ2 = TLorentzVector(0, 0, 0, 0);
0220 for (int i = 0; i < (int)j1.size(); i++) {
0221 double dHT_temp = fabs(j1[i].E() - j2[i].E());
0222 if (dHT_temp < dHT_min) {
0223 dHT_min = dHT_temp;
0224 myJ1 = j1[i];
0225 myJ2 = j2[i];
0226 }
0227 }
0228
0229 jOUT.clear();
0230 if (myJ1.Pt() > myJ2.Pt()) {
0231 jOUT.push_back(myJ1);
0232 jOUT.push_back(myJ2);
0233 } else {
0234 jOUT.push_back(myJ2);
0235 jOUT.push_back(myJ1);
0236 }
0237 status = 1;
0238 }
0239
0240 }