PatAlgos/src/HemisphereAlgo.cc

0001 #include "PhysicsTools/PatAlgos/interface/HemisphereAlgo.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "DataFormats/Math/interface/deltaR.h"
0007
0008 using namespace std;
0009 using std::vector;
0010
0011 HemisphereAlgo::HemisphereAlgo(const std::vector<reco::CandidatePtr>& componentPtr_,
0012                                const int seed_method,
0013                                const int hemisphere_association_method)
0014     : Object(componentPtr_),
0015       Object_Group(),
0016       Axis1(),
0017       Axis2(),
0018       seed_meth(seed_method),
0019       hemi_meth(hemisphere_association_method),
0020       status(0) {
0021   for (int i = 0; i < (int)Object.size(); i++) {
0022     Object_Noseed.push_back(0);
0023   }
0024 }
0025
0026 vector<float> HemisphereAlgo::getAxis1() {
0027   if (status != 1) {
0028     this->reconstruct();
0029   }
0030   return Axis1;
0031 }
0032 vector<float> HemisphereAlgo::getAxis2() {
0033   if (status != 1) {
0034     this->reconstruct();
0035   }
0036   return Axis2;
0037 }
0038
0039 vector<int> HemisphereAlgo::getGrouping() {
0040   if (status != 1) {
0041     this->reconstruct();
0042   }
0043   return Object_Group;
0044 }
0045
0046 int HemisphereAlgo::reconstruct() {
0047   int vsize = (int)Object.size();
0048
0049   LogDebug("HemisphereAlgo") << " HemisphereAlgo method ";
0050
0051   Object_Group.clear();
0052   Axis1.clear();
0053   Axis2.clear();
0054
0055   for (int j = 0; j < vsize; j++) {
0056     Object_Group.push_back(0);
0057   }
0058
0059   for (int j = 0; j < 5; j++) {
0060     Axis1.push_back(0);
0061     Axis2.push_back(0);
0062   }
0063
0064   for (int i = 0; i < vsize; i++) {
0065     if ((*(Object)[i]).p() > (*Object[i]).energy() + 0.001) {
0066       edm::LogWarning("HemisphereAlgo") << "Object " << i << " has E = " << (*Object[i]).energy()
0067                                         << " less than P = " << (*Object[i]).p();
0068     }
0069   }
0070
0071   LogDebug("HemisphereAlgo") << " Seeding method = " << seed_meth;
0072   int I_Max = -1;
0073   int J_Max = -1;
0074
0075   if (seed_meth == 1) {
0076     float P_Max = 0.;
0077     float DeltaRP_Max = 0.;
0078
0079     // take highest momentum object as first axis
0080     for (int i = 0; i < vsize; i++) {
0081       Object_Group[i] = 0;
0082       if (Object_Noseed[i] == 0 && P_Max < (*Object[i]).p()) {
0083         P_Max = (*Object[i]).p();
0084         I_Max = i;
0085       }
0086     }
0087
0088     Axis1[0] = (*Object[I_Max]).px() / (*Object[I_Max]).p();
0089     Axis1[1] = (*Object[I_Max]).py() / (*Object[I_Max]).p();
0090     Axis1[2] = (*Object[I_Max]).pz() / (*Object[I_Max]).p();
0091     Axis1[3] = (*Object[I_Max]).p();
0092     Axis1[4] = (*Object[I_Max]).energy();
0093
0094     // take as second axis
0095     for (int i = 0; i < vsize; i++) {
0096       float DeltaR = deltaR((*Object[i]).eta(), (*Object[i]).phi(), (*Object[I_Max]).eta(), (*Object[I_Max]).phi());
0097
0098       if (Object_Noseed[i] == 0 && DeltaR > 0.5) {
0099         float DeltaRP = DeltaR * (*Object[i]).p();
0100         if (DeltaRP > DeltaRP_Max) {
0101           DeltaRP_Max = DeltaRP;
0102           J_Max = i;
0103         }
0104       }
0105     }
0106
0107     if (J_Max >= 0) {
0108       Axis2[0] = (*Object[J_Max]).px() / (*Object[J_Max]).p();
0109       Axis2[1] = (*Object[J_Max]).py() / (*Object[J_Max]).p();
0110       Axis2[2] = (*Object[J_Max]).pz() / (*Object[J_Max]).p();
0111       Axis2[3] = (*Object[J_Max]).p();
0112       Axis2[4] = (*Object[J_Max]).energy();
0113
0114     } else {
0115       return 0;
0116     }
0117     LogDebug("HemisphereAlgo") << " Axis 1 is Object = " << I_Max;
0118     LogDebug("HemisphereAlgo") << " Axis 2 is Object = " << J_Max;
0119
0120   } else if ((seed_meth == 2) | (seed_meth == 3)) {
0121     float Mass_Max = 0.;
0122     float InvariantMass = 0.;
0123
0124     // maximize the invariant mass of two objects
0125     for (int i = 0; i < vsize; i++) {
0126       Object_Group[i] = 0;
0127       if (Object_Noseed[i] == 0) {
0128         for (int j = i + 1; j < vsize; j++) {
0129           if (Object_Noseed[j] == 0) {
0130             // either the invariant mass
0131             if (seed_meth == 2) {
0132               InvariantMass =
0133                   ((*Object[i]).energy() + (*Object[j]).energy()) * ((*Object[i]).energy() + (*Object[j]).energy()) -
0134                   ((*Object[i]).px() + (*Object[j]).px()) * ((*Object[i]).px() + (*Object[j]).px()) -
0135                   ((*Object[i]).py() + (*Object[j]).py()) * ((*Object[i]).py() + (*Object[j]).py()) -
0136                   ((*Object[i]).pz() + (*Object[j]).pz()) * ((*Object[i]).pz() + (*Object[j]).pz());
0137             }
0138             // or the transverse mass
0139             else if (seed_meth == 3) {
0140               float pti = sqrt((*Object[i]).px() * (*Object[i]).px() + (*Object[i]).py() * (*Object[i]).py());
0141               float ptj = sqrt((*Object[j]).px() * (*Object[j]).px() + (*Object[j]).py() * (*Object[j]).py());
0142               InvariantMass =
0143                   2. * (pti * ptj - (*Object[i]).px() * (*Object[j]).px() - (*Object[i]).py() * (*Object[j]).py());
0144             }
0145             if (Mass_Max < InvariantMass) {
0146               Mass_Max = InvariantMass;
0147               I_Max = i;
0148               J_Max = j;
0149             }
0150           }
0151         }
0152       }
0153     }
0154
0155     if (J_Max > 0) {
0156       Axis1[0] = (*Object[I_Max]).px() / (*Object[I_Max]).p();
0157       Axis1[1] = (*Object[I_Max]).py() / (*Object[I_Max]).p();
0158       Axis1[2] = (*Object[I_Max]).pz() / (*Object[I_Max]).p();
0159
0160       Axis1[3] = (*Object[I_Max]).p();
0161       Axis1[4] = (*Object[I_Max]).energy();
0162
0163       Axis2[0] = (*Object[J_Max]).px() / (*Object[J_Max]).p();
0164       Axis2[1] = (*Object[J_Max]).py() / (*Object[J_Max]).p();
0165       Axis2[2] = (*Object[J_Max]).pz() / (*Object[J_Max]).p();
0166
0167       Axis2[3] = (*Object[J_Max]).p();
0168       Axis2[4] = (*Object[J_Max]).energy();
0169
0170     } else {
0171       return 0;
0172     }
0173     LogDebug("HemisphereAlgo") << " Axis 1 is Object = " << I_Max;
0174     LogDebug("HemisphereAlgo") << " Axis 2 is Object = " << J_Max;
0175
0176   } else {
0177     throw cms::Exception("Configuration") << "Please give a valid seeding method!";
0178   }
0179
0180   // seeding done
0181   // now do the hemisphere association
0182
0183   LogDebug("HemisphereAlgo") << " Association method = " << hemi_meth;
0184
0185   int numLoop = 0;
0186   bool I_Move = true;
0187
0188   while (I_Move && (numLoop < 100)) {
0189     I_Move = false;
0190     numLoop++;
0191     LogDebug("HemisphereAlgo") << " Iteration = " << numLoop;
0192
0193     float Sum1_Px = 0.;
0194     float Sum1_Py = 0.;
0195     float Sum1_Pz = 0.;
0196     float Sum1_E = 0.;
0197     float Sum2_Px = 0.;
0198     float Sum2_Py = 0.;
0199     float Sum2_Pz = 0.;
0200     float Sum2_E = 0.;
0201
0202     if (hemi_meth == 1) {
0203       for (int i = 0; i < vsize; i++) {
0204         float P_Long1 = (*Object[i]).px() * Axis1[0] + (*Object[i]).py() * Axis1[1] + (*Object[i]).pz() * Axis1[2];
0205         float P_Long2 = (*Object[i]).px() * Axis2[0] + (*Object[i]).py() * Axis2[1] + (*Object[i]).pz() * Axis2[2];
0206         if (P_Long1 >= P_Long2) {
0207           if (Object_Group[i] != 1) {
0208             I_Move = true;
0209           }
0210           Object_Group[i] = 1;
0211           Sum1_Px += (*Object[i]).px();
0212           Sum1_Py += (*Object[i]).py();
0213           Sum1_Pz += (*Object[i]).pz();
0214           Sum1_E += (*Object[i]).energy();
0215         } else {
0216           if (Object_Group[i] != 2) {
0217             I_Move = true;
0218           }
0219           Object_Group[i] = 2;
0220           Sum2_Px += (*Object[i]).px();
0221           Sum2_Py += (*Object[i]).py();
0222           Sum2_Pz += (*Object[i]).pz();
0223           Sum2_E += (*Object[i]).energy();
0224         }
0225       }
0226
0227     } else if (hemi_meth == 2 || hemi_meth == 3) {
0228       for (int i = 0; i < vsize; i++) {
0229         if (i == I_Max) {
0230           Object_Group[i] = 1;
0231           Sum1_Px += (*Object[i]).px();
0232           Sum1_Py += (*Object[i]).py();
0233           Sum1_Pz += (*Object[i]).pz();
0234           Sum1_E += (*Object[i]).energy();
0235         } else if (i == J_Max) {
0236           Object_Group[i] = 2;
0237           Sum2_Px += (*Object[i]).px();
0238           Sum2_Py += (*Object[i]).py();
0239           Sum2_Pz += (*Object[i]).pz();
0240           Sum2_E += (*Object[i]).energy();
0241         } else {
0242           if (!I_Move) {
0243             float NewAxis1_Px = Axis1[0] * Axis1[3];
0244             float NewAxis1_Py = Axis1[1] * Axis1[3];
0245             float NewAxis1_Pz = Axis1[2] * Axis1[3];
0246             float NewAxis1_E = Axis1[4];
0247
0248             float NewAxis2_Px = Axis2[0] * Axis2[3];
0249             float NewAxis2_Py = Axis2[1] * Axis2[3];
0250             float NewAxis2_Pz = Axis2[2] * Axis2[3];
0251             float NewAxis2_E = Axis2[4];
0252
0253             if (Object_Group[i] == 1) {
0254               NewAxis1_Px = NewAxis1_Px - (*Object[i]).px();
0255               NewAxis1_Py = NewAxis1_Py - (*Object[i]).py();
0256               NewAxis1_Pz = NewAxis1_Pz - (*Object[i]).pz();
0257               NewAxis1_E = NewAxis1_E - (*Object[i]).energy();
0258
0259             } else if (Object_Group[i] == 2) {
0260               NewAxis2_Px = NewAxis2_Px - (*Object[i]).px();
0261               NewAxis2_Py = NewAxis2_Py - (*Object[i]).py();
0262               NewAxis2_Pz = NewAxis2_Pz - (*Object[i]).pz();
0263               NewAxis2_E = NewAxis2_E - (*Object[i]).energy();
0264             }
0265
0266             float mass1 =
0267                 NewAxis1_E -
0268                 (((*Object[i]).px() * NewAxis1_Px + (*Object[i]).py() * NewAxis1_Py + (*Object[i]).pz() * NewAxis1_Pz) /
0269                  (*Object[i]).p());
0270
0271             float mass2 =
0272                 NewAxis2_E -
0273                 (((*Object[i]).px() * NewAxis2_Px + (*Object[i]).py() * NewAxis2_Py + (*Object[i]).pz() * NewAxis2_Pz) /
0274                  (*Object[i]).p());
0275
0276             if (hemi_meth == 3) {
0277               mass1 *= NewAxis1_E / ((NewAxis1_E + (*Object[i]).energy()) * (NewAxis1_E + (*Object[i]).energy()));
0278
0279               mass2 *= NewAxis2_E / ((NewAxis2_E + (*Object[i]).energy()) * (NewAxis2_E + (*Object[i]).energy()));
0280             }
0281
0282             if (mass1 < mass2) {
0283               if (Object_Group[i] != 1) {
0284                 I_Move = true;
0285               }
0286               Object_Group[i] = 1;
0287
0288               Sum1_Px += (*Object[i]).px();
0289               Sum1_Py += (*Object[i]).py();
0290               Sum1_Pz += (*Object[i]).pz();
0291               Sum1_E += (*Object[i]).energy();
0292             } else {
0293               if (Object_Group[i] != 2) {
0294                 I_Move = true;
0295               }
0296               Object_Group[i] = 2;
0297               Sum2_Px += (*Object[i]).px();
0298               Sum2_Py += (*Object[i]).py();
0299               Sum2_Pz += (*Object[i]).pz();
0300               Sum2_E += (*Object[i]).energy();
0301             }
0302
0303           } else {
0304             if (Object_Group[i] == 1) {
0305               Sum1_Px += (*Object[i]).px();
0306               Sum1_Py += (*Object[i]).py();
0307               Sum1_Pz += (*Object[i]).pz();
0308               Sum1_E += (*Object[i]).energy();
0309             }
0310             if (Object_Group[i] == 2) {
0311               Sum2_Px += (*Object[i]).px();
0312               Sum2_Py += (*Object[i]).py();
0313               Sum2_Pz += (*Object[i]).pz();
0314               Sum2_E += (*Object[i]).energy();
0315             }
0316           }
0317         }
0318       }
0319     } else {
0320       throw cms::Exception("Configuration") << "Please give a valid hemisphere association method!";
0321     }
0322
0323     // recomputing the axes
0324
0325     Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
0326     if (Axis1[3] < 0.0001) {
0327       edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 1! ";
0328     } else {
0329       Axis1[0] = Sum1_Px / Axis1[3];
0330       Axis1[1] = Sum1_Py / Axis1[3];
0331       Axis1[2] = Sum1_Pz / Axis1[3];
0332       Axis1[4] = Sum1_E;
0333     }
0334
0335     Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
0336     if (Axis2[3] < 0.0001) {
0337       edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 2!";
0338     } else {
0339       Axis2[0] = Sum2_Px / Axis2[3];
0340       Axis2[1] = Sum2_Py / Axis2[3];
0341       Axis2[2] = Sum2_Pz / Axis2[3];
0342       Axis2[4] = Sum2_E;
0343     }
0344
0345     LogDebug("HemisphereAlgo") << " Grouping = ";
0346     for (int i = 0; i < vsize; i++) {
0347       LogTrace("HemisphereAlgo") << "  " << Object_Group[i];
0348     }
0349     LogTrace("HemisphereAlgo") << std::endl;
0350   }
0351
0352   status = 1;
0353   return 1;
0354 }