Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:03

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_P = 0.;
0197     float Sum1_E = 0.;
0198     float Sum2_Px = 0.;
0199     float Sum2_Py = 0.;
0200     float Sum2_Pz = 0.;
0201     float Sum2_P = 0.;
0202     float Sum2_E = 0.;
0203 
0204     if (hemi_meth == 1) {
0205       for (int i = 0; i < vsize; i++) {
0206         float P_Long1 = (*Object[i]).px() * Axis1[0] + (*Object[i]).py() * Axis1[1] + (*Object[i]).pz() * Axis1[2];
0207         float P_Long2 = (*Object[i]).px() * Axis2[0] + (*Object[i]).py() * Axis2[1] + (*Object[i]).pz() * Axis2[2];
0208         if (P_Long1 >= P_Long2) {
0209           if (Object_Group[i] != 1) {
0210             I_Move = true;
0211           }
0212           Object_Group[i] = 1;
0213           Sum1_Px += (*Object[i]).px();
0214           Sum1_Py += (*Object[i]).py();
0215           Sum1_Pz += (*Object[i]).pz();
0216           Sum1_P += (*Object[i]).p();
0217           Sum1_E += (*Object[i]).energy();
0218         } else {
0219           if (Object_Group[i] != 2) {
0220             I_Move = true;
0221           }
0222           Object_Group[i] = 2;
0223           Sum2_Px += (*Object[i]).px();
0224           Sum2_Py += (*Object[i]).py();
0225           Sum2_Pz += (*Object[i]).pz();
0226           Sum2_P += (*Object[i]).p();
0227           Sum2_E += (*Object[i]).energy();
0228         }
0229       }
0230 
0231     } else if (hemi_meth == 2 || hemi_meth == 3) {
0232       for (int i = 0; i < vsize; i++) {
0233         if (i == I_Max) {
0234           Object_Group[i] = 1;
0235           Sum1_Px += (*Object[i]).px();
0236           Sum1_Py += (*Object[i]).py();
0237           Sum1_Pz += (*Object[i]).pz();
0238           Sum1_P += (*Object[i]).p();
0239           Sum1_E += (*Object[i]).energy();
0240         } else if (i == J_Max) {
0241           Object_Group[i] = 2;
0242           Sum2_Px += (*Object[i]).px();
0243           Sum2_Py += (*Object[i]).py();
0244           Sum2_Pz += (*Object[i]).pz();
0245           Sum2_P += (*Object[i]).p();
0246           Sum2_E += (*Object[i]).energy();
0247         } else {
0248           if (!I_Move) {
0249             float NewAxis1_Px = Axis1[0] * Axis1[3];
0250             float NewAxis1_Py = Axis1[1] * Axis1[3];
0251             float NewAxis1_Pz = Axis1[2] * Axis1[3];
0252             float NewAxis1_E = Axis1[4];
0253 
0254             float NewAxis2_Px = Axis2[0] * Axis2[3];
0255             float NewAxis2_Py = Axis2[1] * Axis2[3];
0256             float NewAxis2_Pz = Axis2[2] * Axis2[3];
0257             float NewAxis2_E = Axis2[4];
0258 
0259             if (Object_Group[i] == 1) {
0260               NewAxis1_Px = NewAxis1_Px - (*Object[i]).px();
0261               NewAxis1_Py = NewAxis1_Py - (*Object[i]).py();
0262               NewAxis1_Pz = NewAxis1_Pz - (*Object[i]).pz();
0263               NewAxis1_E = NewAxis1_E - (*Object[i]).energy();
0264 
0265             } else if (Object_Group[i] == 2) {
0266               NewAxis2_Px = NewAxis2_Px - (*Object[i]).px();
0267               NewAxis2_Py = NewAxis2_Py - (*Object[i]).py();
0268               NewAxis2_Pz = NewAxis2_Pz - (*Object[i]).pz();
0269               NewAxis2_E = NewAxis2_E - (*Object[i]).energy();
0270             }
0271 
0272             float mass1 =
0273                 NewAxis1_E -
0274                 (((*Object[i]).px() * NewAxis1_Px + (*Object[i]).py() * NewAxis1_Py + (*Object[i]).pz() * NewAxis1_Pz) /
0275                  (*Object[i]).p());
0276 
0277             float mass2 =
0278                 NewAxis2_E -
0279                 (((*Object[i]).px() * NewAxis2_Px + (*Object[i]).py() * NewAxis2_Py + (*Object[i]).pz() * NewAxis2_Pz) /
0280                  (*Object[i]).p());
0281 
0282             if (hemi_meth == 3) {
0283               mass1 *= NewAxis1_E / ((NewAxis1_E + (*Object[i]).energy()) * (NewAxis1_E + (*Object[i]).energy()));
0284 
0285               mass2 *= NewAxis2_E / ((NewAxis2_E + (*Object[i]).energy()) * (NewAxis2_E + (*Object[i]).energy()));
0286             }
0287 
0288             if (mass1 < mass2) {
0289               if (Object_Group[i] != 1) {
0290                 I_Move = true;
0291               }
0292               Object_Group[i] = 1;
0293 
0294               Sum1_Px += (*Object[i]).px();
0295               Sum1_Py += (*Object[i]).py();
0296               Sum1_Pz += (*Object[i]).pz();
0297               Sum1_P += (*Object[i]).p();
0298               Sum1_E += (*Object[i]).energy();
0299             } else {
0300               if (Object_Group[i] != 2) {
0301                 I_Move = true;
0302               }
0303               Object_Group[i] = 2;
0304               Sum2_Px += (*Object[i]).px();
0305               Sum2_Py += (*Object[i]).py();
0306               Sum2_Pz += (*Object[i]).pz();
0307               Sum2_P += (*Object[i]).p();
0308               Sum2_E += (*Object[i]).energy();
0309             }
0310 
0311           } else {
0312             if (Object_Group[i] == 1) {
0313               Sum1_Px += (*Object[i]).px();
0314               Sum1_Py += (*Object[i]).py();
0315               Sum1_Pz += (*Object[i]).pz();
0316               Sum1_P += (*Object[i]).p();
0317               Sum1_E += (*Object[i]).energy();
0318             }
0319             if (Object_Group[i] == 2) {
0320               Sum2_Px += (*Object[i]).px();
0321               Sum2_Py += (*Object[i]).py();
0322               Sum2_Pz += (*Object[i]).pz();
0323               Sum2_P += (*Object[i]).p();
0324               Sum2_E += (*Object[i]).energy();
0325             }
0326           }
0327         }
0328       }
0329     } else {
0330       throw cms::Exception("Configuration") << "Please give a valid hemisphere association method!";
0331     }
0332 
0333     // recomputing the axes
0334 
0335     Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
0336     if (Axis1[3] < 0.0001) {
0337       edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 1! ";
0338     } else {
0339       Axis1[0] = Sum1_Px / Axis1[3];
0340       Axis1[1] = Sum1_Py / Axis1[3];
0341       Axis1[2] = Sum1_Pz / Axis1[3];
0342       Axis1[4] = Sum1_E;
0343     }
0344 
0345     Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
0346     if (Axis2[3] < 0.0001) {
0347       edm::LogWarning("HemisphereAlgo") << "ZERO objects in group 2!";
0348     } else {
0349       Axis2[0] = Sum2_Px / Axis2[3];
0350       Axis2[1] = Sum2_Py / Axis2[3];
0351       Axis2[2] = Sum2_Pz / Axis2[3];
0352       Axis2[4] = Sum2_E;
0353     }
0354 
0355     LogDebug("HemisphereAlgo") << " Grouping = ";
0356     for (int i = 0; i < vsize; i++) {
0357       LogTrace("HemisphereAlgo") << "  " << Object_Group[i];
0358     }
0359     LogTrace("HemisphereAlgo") << std::endl;
0360   }
0361 
0362   status = 1;
0363   return 1;
0364 }