Heppy/src/Megajet.cc

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   // constructor specifying the association methods
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   // constructor without specification of the seed and association methods
0031   // in this case, the latter must be given by calling SetMethod before invoking Combine()
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     // clear some vectors if method Combine() is called again
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     // default value (in case none is found)
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     //  myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
0145     //  myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
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     // default value (in case none is found)
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     //  myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
0173     //  myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
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     // default value (in case none is found)
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     //  myJ1.SetPtEtaPhiM(myJ1.Pt(),myJ1.Eta(),myJ1.Phi(),0.0);
0202     //  myJ2.SetPtEtaPhiM(myJ2.Pt(),myJ2.Eta(),myJ2.Phi(),0.0);
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     // default value (in case none is found)
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 }  // namespace heppy