Heppy/src/Hemisphere.cc

0001 #include "PhysicsTools/Heppy/interface/Hemisphere.h"
0002 #include "DataFormats/Math/interface/deltaPhi.h"
0003 #include "DataFormats/Math/interface/deltaR.h"
0004
0005 using namespace std;
0006
0007 using std::cout;
0008 using std::endl;
0009 using std::vector;
0010
0011 namespace heppy {
0012
0013   // constructor specifying the seed and association methods

0014   Hemisphere::Hemisphere(vector<float> Px_vector,
0015                          vector<float> Py_vector,
0016                          vector<float> Pz_vector,
0017                          vector<float> E_vector,
0018                          int seed_method,
0019                          int hemisphere_association_method)
0020       : Object_Px(Px_vector),
0021         Object_Py(Py_vector),
0022         Object_Pz(Pz_vector),
0023         Object_E(E_vector),
0024         seed_meth(seed_method),
0025         hemi_meth(hemisphere_association_method),
0026         status(0),
0027         dRminSeed1(0.5),
0028         nItermax(100),
0029         rejectISR(0),
0030         rejectISRPt(0),
0031         rejectISRPtmax(10000.),
0032         rejectISRDR(0),
0033         rejectISRDRmax(100.),
0034         dbg(0) {
0035     for (int i = 0; i < (int)Object_Px.size(); i++) {
0036       Object_Noseed.push_back(0);
0037       Object_Noassoc.push_back(0);
0038     }
0039     numLoop = 0;
0040   }
0041
0042   // constructor without specification of the seed and association methods

0043   // in this case, the latter must be given by calling SetMethod before invoking reconstruct()

0044   Hemisphere::Hemisphere(vector<float> Px_vector,
0045                          vector<float> Py_vector,
0046                          vector<float> Pz_vector,
0047                          vector<float> E_vector)
0048       : Object_Px(Px_vector),
0049         Object_Py(Py_vector),
0050         Object_Pz(Pz_vector),
0051         Object_E(E_vector),
0052         seed_meth(0),
0053         hemi_meth(0),
0054         status(0),
0055         dRminSeed1(0.5),
0056         nItermax(100),
0057         rejectISR(0),
0058         rejectISRPt(0),
0059         rejectISRPtmax(10000.),
0060         rejectISRDR(0),
0061         rejectISRDRmax(100.),
0062         dbg(0) {
0063     for (int i = 0; i < (int)Object_Px.size(); i++) {
0064       Object_Noseed.push_back(0);
0065       Object_Noassoc.push_back(0);
0066     }
0067     numLoop = 0;
0068   }
0069
0070   vector<float> Hemisphere::getAxis1() {
0071     if (status != 1) {
0072       if (rejectISR == 0) {
0073         this->Reconstruct();
0074       } else {
0075         this->RejectISR();
0076       }
0077     }
0078     return Axis1;
0079   }
0080   vector<float> Hemisphere::getAxis2() {
0081     if (status != 1) {
0082       if (rejectISR == 0) {
0083         this->Reconstruct();
0084       } else {
0085         this->RejectISR();
0086       }
0087     }
0088     return Axis2;
0089   }
0090
0091   vector<int> Hemisphere::getGrouping() {
0092     if (status != 1) {
0093       if (rejectISR == 0) {
0094         this->Reconstruct();
0095       } else {
0096         this->RejectISR();
0097       }
0098     }
0099     return Object_Group;
0100   }
0101
0102   int Hemisphere::Reconstruct() {
0103     // Performs the actual hemisphere reconstrucion

0104     //

0105     // definition of the vectors used internally:

0106     // Object_xxx :

0107     //        xxx = Px, Py, Pz, E for input values

0108     //        xxx = P, Pt, Eta, Phi, Group for internal use

0109     // Axis1 : final hemisphere axis 1

0110     // Axis2 : final hemisphere axis 2

0111     // Sum1_xxx : hemisphere 1 being updated during the association iterations

0112     // Sum2_xxx : hemisphere 2 being updated during the association iterations

0113     // NewAxis1_xxx, NewAxis1_xxx : temporary axes for calculation in association methods 2 and 3

0114
0115     numLoop = 0;  // initialize numLoop for Zero

0116     int vsize = (int)Object_Px.size();
0117     if ((int)Object_Py.size() != vsize || (int)Object_Pz.size() != vsize) {
0118       cout << "WARNING!!!!! Input vectors have different size! Fix it!" << endl;
0119       return 0;
0120     }
0121     if (dbg > 0) {
0122       //    cout << " Hemisphere method, vsn = " << hemivsn << endl;

0123       cout << " Hemisphere method " << endl;
0124     }
0125
0126     // clear some vectors if method reconstruct() is called again

0127     if (!Object_P.empty()) {
0128       Object_P.clear();
0129       Object_Pt.clear();
0130       Object_Eta.clear();
0131       Object_Phi.clear();
0132       Object_Group.clear();
0133       Axis1.clear();
0134       Axis2.clear();
0135     }
0136     // initialize the vectors

0137     for (int j = 0; j < vsize; ++j) {
0138       Object_P.push_back(0);
0139       Object_Pt.push_back(0);
0140       Object_Eta.push_back(0);
0141       Object_Phi.push_back(0);
0142       Object_Group.push_back(0);
0143     }
0144     for (int j = 0; j < 5; ++j) {
0145       Axis1.push_back(0);
0146       Axis2.push_back(0);
0147     }
0148
0149     // compute additional quantities for vectors Object_xxx

0150     float theta;
0151     for (int i = 0; i < vsize; ++i) {
0152       Object_P[i] = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i] + Object_Pz[i] * Object_Pz[i]);
0153       if (Object_P[i] > Object_E[i] + 0.001) {
0154         cout << "WARNING!!!!! Object " << i << " has E = " << Object_E[i] << " less than P = " << Object_P[i]
0155              << " *** Fix it!" << endl;
0156         return 0;
0157       }
0158       Object_Pt[i] = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]);
0159       // protection for div by 0

0160       if (fabs(Object_Pz[i]) > 0.001) {
0161         theta = atan(sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]) / Object_Pz[i]);
0162       } else {
0163         theta = 1.570796327;
0164       }
0165       if (theta < 0.) {
0166         theta = theta + 3.141592654;
0167       }
0168       Object_Eta[i] = -log(tan(0.5 * theta));
0169       Object_Phi[i] = atan2(Object_Py[i], Object_Px[i]);
0170       if (dbg > 0) {
0171         cout << " Object " << i << " Eta = " << Object_Eta[i] << " Phi = " << Object_Phi[i] << endl;
0172       }
0173     }
0174
0175     if (dbg > 0) {
0176       cout << endl;
0177       cout << " Seeding method = " << seed_meth << endl;
0178     }
0179     // I_Max and J_Max are indices of the seeds in the vectors

0180     int I_Max = -1;
0181     int J_Max = -1;
0182
0183     // determine the seeds for seed method 1

0184     if (seed_meth == 1) {
0185       float P_Max = 0.;
0186       float DeltaRP_Max = 0.;
0187
0188       // take highest momentum object as first seed

0189       for (int i = 0; i < vsize; ++i) {
0190         Object_Group[i] = 0;
0191         //cout << "Object_Px[i] = " << Object_Px[i] << ", Object_Py[i] = " << Object_Py[i]

0192         //<< ", Object_Pz[i] = " << Object_Pz[i] << "  << endl;

0193         if (Object_Noseed[i] == 0 && P_Max < Object_P[i]) {
0194           P_Max = Object_P[i];
0195           I_Max = i;
0196         }
0197       }
0198       // if 1st seed is found, save it as initial hemisphere 1 axis

0199       if (I_Max >= 0) {
0200         Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
0201         Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
0202         Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
0203         Axis1[3] = Object_P[I_Max];
0204         Axis1[4] = Object_E[I_Max];
0205       } else {
0206         // cout << " This is an empty event." << endl;

0207         return 0;
0208       }
0209
0210       // take as second seed the object with largest DR*P w.r.t. the first seed

0211       for (int i = 0; i < vsize; ++i) {
0212         //          float DeltaR = sqrt((Object_Eta[i] - Object_Eta[I_Max])*(Object_Eta[i] - Object_Eta[I_Max])

0213         //              + (Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max]))*(Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max]))  );

0214         float DeltaR =
0215             sqrt((Object_Eta[i] - Object_Eta[I_Max]) * (Object_Eta[i] - Object_Eta[I_Max]) +
0216                  (deltaPhi(Object_Phi[i], Object_Phi[I_Max])) * (deltaPhi(Object_Phi[i], Object_Phi[I_Max])));
0217         if (Object_Noseed[i] == 0 && DeltaR > dRminSeed1) {
0218           float DeltaRP = DeltaR * Object_P[i];
0219           if (DeltaRP > DeltaRP_Max) {
0220             DeltaRP_Max = DeltaRP;
0221             J_Max = i;
0222           }
0223         }
0224       }
0225       // if 2nd seed is found, save it as initial hemisphere 2 axis

0226       if (J_Max >= 0) {
0227         Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
0228         Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
0229         Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
0230         Axis2[3] = Object_P[J_Max];
0231         Axis2[4] = Object_E[J_Max];
0232       } else {
0233         // cout << " This is a MONOJET." << endl;

0234         return 0;
0235       }
0236       if (dbg > 0) {
0237         cout << " Axis 1 is Object = " << I_Max << endl;
0238         cout << " Axis 2 is Object = " << J_Max << endl;
0239       }
0240
0241       // determine the seeds for seed methods 2 and 3

0242     } else if (seed_meth == 2 || seed_meth == 3) {
0243       float Mass_Max = 0.;
0244       float InvariantMass = 0.;
0245
0246       // maximize the invariant mass of two objects

0247       for (int i = 0; i < vsize; ++i) {
0248         Object_Group[i] = 0;
0249         if (Object_Noseed[i] == 0) {
0250           for (int j = i + 1; j < vsize; ++j) {
0251             if (Object_Noseed[j] == 0) {
0252               // either the invariant mass

0253               if (seed_meth == 2) {
0254                 InvariantMass = (Object_E[i] + Object_E[j]) * (Object_E[i] + Object_E[j]) -
0255                                 (Object_Px[i] + Object_Px[j]) * (Object_Px[i] + Object_Px[j]) -
0256                                 (Object_Py[i] + Object_Py[j]) * (Object_Py[i] + Object_Py[j]) -
0257                                 (Object_Pz[i] + Object_Pz[j]) * (Object_Pz[i] + Object_Pz[j]);
0258               }
0259               // or the transverse mass

0260               else if (seed_meth == 3) {
0261                 float pti = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]);
0262                 float ptj = sqrt(Object_Px[j] * Object_Px[j] + Object_Py[j] * Object_Py[j]);
0263                 InvariantMass = 2. * (pti * ptj - Object_Px[i] * Object_Px[j] - Object_Py[i] * Object_Py[j]);
0264               }
0265               if (Mass_Max < InvariantMass) {
0266                 Mass_Max = InvariantMass;
0267                 I_Max = i;
0268                 J_Max = j;
0269               }
0270             }
0271           }
0272         }
0273       }
0274
0275       // if both seeds are found, save them as initial hemisphere axes

0276       if (J_Max > 0) {
0277         Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
0278         Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
0279         Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
0280         Axis1[3] = Object_P[I_Max];
0281         Axis1[4] = Object_E[I_Max];
0282
0283         Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
0284         Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
0285         Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
0286         Axis2[3] = Object_P[J_Max];
0287         Axis2[4] = Object_E[J_Max];
0288       } else {
0289         // cout << " This is a MONOJET." << endl;

0290         return 0;
0291       }
0292       if (dbg > 0) {
0293         cout << " Axis 1 is Object = " << I_Max << endl;
0294         cout << " Axis 2 is Object = " << J_Max << endl;
0295       }
0296
0297     } else if (seed_meth == 4) {
0298       float P_Max1 = 0.;
0299       float P_Max2 = 0.;
0300
0301       // take largest Pt object as first seed

0302       for (int i = 0; i < vsize; ++i) {
0303         Object_Group[i] = 0;
0304         if (Object_Noseed[i] == 0 && P_Max1 < Object_Pt[i]) {
0305           P_Max1 = Object_Pt[i];
0306           I_Max = i;
0307         }
0308       }
0309       if (I_Max < 0)
0310         return 0;
0311
0312       // take second largest Pt object as second seed, but require dR(seed1, seed2) > dRminSeed1

0313       for (int i = 0; i < vsize; ++i) {
0314         if (i == I_Max)
0315           continue;
0316         //          float DeltaR = Util::GetDeltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);

0317         float DeltaR = deltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);
0318         if (Object_Noseed[i] == 0 && P_Max2 < Object_Pt[i] && DeltaR > dRminSeed1) {
0319           P_Max2 = Object_Pt[i];
0320           J_Max = i;
0321         }
0322       }
0323       if (J_Max < 0)
0324         return 0;
0325
0326       // save first seed as initial hemisphere 1 axis

0327       if (I_Max >= 0) {
0328         Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
0329         Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
0330         Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
0331         Axis1[3] = Object_P[I_Max];
0332         Axis1[4] = Object_E[I_Max];
0333       }
0334
0335       // save second seed as initial hemisphere 2 axis

0336       if (J_Max >= 0) {
0337         Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
0338         Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
0339         Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
0340         Axis2[3] = Object_P[J_Max];
0341         Axis2[4] = Object_E[J_Max];
0342       }
0343
0344       if (dbg > 0) {
0345         cout << " Axis 1 is Object = " << I_Max << " with Pt " << Object_Pt[I_Max] << endl;
0346         cout << " Axis 2 is Object = " << J_Max << " with Pt " << Object_Pt[J_Max] << endl;
0347       }
0348
0349     } else if (!(seed_meth == 0 && (hemi_meth == 8 || hemi_meth == 9))) {
0350       cout << "Please give a valid seeding method!" << endl;
0351       return 0;
0352     }
0353
0354     // seeding done

0355     // now do the hemisphere association

0356
0357     if (dbg > 0) {
0358       cout << endl;
0359       cout << " Association method = " << hemi_meth << endl;
0360     }
0361
0362     bool I_Move = true;
0363
0364     // iterate to associate all objects to hemispheres (methods 1 to 3 only)

0365     //   until no objects are moved from one to the other hemisphere

0366     //   or the maximum number of iterations is reached

0367     while (I_Move && (numLoop < nItermax) && hemi_meth != 8 && hemi_meth != 9) {
0368       I_Move = false;
0369       numLoop++;
0370       if (dbg > 0) {
0371         cout << " Iteration = " << numLoop << endl;
0372       }
0373       if (numLoop == nItermax - 1) {
0374         cout << " Hemishpere: warning - reaching max number of iterations " << endl;
0375       }
0376
0377       // initialize the current sums of Px, Py, Pz, E for the two hemispheres

0378       float Sum1_Px = 0.;
0379       float Sum1_Py = 0.;
0380       float Sum1_Pz = 0.;
0381       float Sum1_E = 0.;
0382       float Sum2_Px = 0.;
0383       float Sum2_Py = 0.;
0384       float Sum2_Pz = 0.;
0385       float Sum2_E = 0.;
0386
0387       // associate the objects for method 1

0388       if (hemi_meth == 1) {
0389         for (int i = 0; i < vsize; ++i) {
0390           if (Object_Noassoc[i] == 0) {
0391             float P_Long1 = Object_Px[i] * Axis1[0] + Object_Py[i] * Axis1[1] + Object_Pz[i] * Axis1[2];
0392             float P_Long2 = Object_Px[i] * Axis2[0] + Object_Py[i] * Axis2[1] + Object_Pz[i] * Axis2[2];
0393             if (P_Long1 >= P_Long2) {
0394               if (Object_Group[i] != 1) {
0395                 I_Move = true;
0396               }
0397               Object_Group[i] = 1;
0398               Sum1_Px += Object_Px[i];
0399               Sum1_Py += Object_Py[i];
0400               Sum1_Pz += Object_Pz[i];
0401               Sum1_E += Object_E[i];
0402             } else {
0403               if (Object_Group[i] != 2) {
0404                 I_Move = true;
0405               }
0406               Object_Group[i] = 2;
0407               Sum2_Px += Object_Px[i];
0408               Sum2_Py += Object_Py[i];
0409               Sum2_Pz += Object_Pz[i];
0410               Sum2_E += Object_E[i];
0411             }
0412           }
0413         }
0414
0415         // associate the objects for methods 2 and 3

0416       } else if (hemi_meth == 2 || hemi_meth == 3) {
0417         for (int i = 0; i < vsize; ++i) {
0418           // add the seeds to the sums, as they remain fixed

0419           if (i == I_Max) {
0420             Object_Group[i] = 1;
0421             Sum1_Px += Object_Px[i];
0422             Sum1_Py += Object_Py[i];
0423             Sum1_Pz += Object_Pz[i];
0424             Sum1_E += Object_E[i];
0425           } else if (i == J_Max) {
0426             Object_Group[i] = 2;
0427             Sum2_Px += Object_Px[i];
0428             Sum2_Py += Object_Py[i];
0429             Sum2_Pz += Object_Pz[i];
0430             Sum2_E += Object_E[i];
0431
0432             // for the other objects

0433           } else {
0434             if (Object_Noassoc[i] == 0) {
0435               // only 1 object maximum is moved in a given iteration

0436               if (!I_Move) {
0437                 // initialize the new hemispheres as the current ones

0438                 float NewAxis1_Px = Axis1[0] * Axis1[3];
0439                 float NewAxis1_Py = Axis1[1] * Axis1[3];
0440                 float NewAxis1_Pz = Axis1[2] * Axis1[3];
0441                 float NewAxis1_E = Axis1[4];
0442                 float NewAxis2_Px = Axis2[0] * Axis2[3];
0443                 float NewAxis2_Py = Axis2[1] * Axis2[3];
0444                 float NewAxis2_Pz = Axis2[2] * Axis2[3];
0445                 float NewAxis2_E = Axis2[4];
0446                 // subtract the object from its hemisphere

0447                 if (Object_Group[i] == 1) {
0448                   NewAxis1_Px = NewAxis1_Px - Object_Px[i];
0449                   NewAxis1_Py = NewAxis1_Py - Object_Py[i];
0450                   NewAxis1_Pz = NewAxis1_Pz - Object_Pz[i];
0451                   NewAxis1_E = NewAxis1_E - Object_E[i];
0452                 } else if (Object_Group[i] == 2) {
0453                   NewAxis2_Px = NewAxis2_Px - Object_Px[i];
0454                   NewAxis2_Py = NewAxis2_Py - Object_Py[i];
0455                   NewAxis2_Pz = NewAxis2_Pz - Object_Pz[i];
0456                   NewAxis2_E = NewAxis2_E - Object_E[i];
0457                 }
0458
0459                 // compute the invariant mass squared with each hemisphere (method 2)

0460                 float mass1 = NewAxis1_E -
0461                               ((Object_Px[i] * NewAxis1_Px + Object_Py[i] * NewAxis1_Py + Object_Pz[i] * NewAxis1_Pz) /
0462                                Object_P[i]);
0463                 float mass2 = NewAxis2_E -
0464                               ((Object_Px[i] * NewAxis2_Px + Object_Py[i] * NewAxis2_Py + Object_Pz[i] * NewAxis2_Pz) /
0465                                Object_P[i]);
0466                 // or the Lund distance (method 3)

0467                 if (hemi_meth == 3) {
0468                   mass1 *= NewAxis1_E / ((NewAxis1_E + Object_E[i]) * (NewAxis1_E + Object_E[i]));
0469                   mass2 *= NewAxis2_E / ((NewAxis2_E + Object_E[i]) * (NewAxis2_E + Object_E[i]));
0470                 }
0471                 // and associate the object to the best hemisphere and add it to the sum

0472                 if (mass1 < mass2) {
0473                   //if (Object_Group[i] != 1){

0474                   if (Object_Group[i] != 1 && Object_Group[i] != 0) {
0475                     I_Move = true;
0476                   }
0477                   Object_Group[i] = 1;
0478                   Sum1_Px += Object_Px[i];
0479                   Sum1_Py += Object_Py[i];
0480                   Sum1_Pz += Object_Pz[i];
0481                   Sum1_E += Object_E[i];
0482                 } else {
0483                   //if (Object_Group[i] != 2){

0484                   if (Object_Group[i] != 2 && Object_Group[i] != 0) {
0485                     I_Move = true;
0486                   }
0487                   Object_Group[i] = 2;
0488                   Sum2_Px += Object_Px[i];
0489                   Sum2_Py += Object_Py[i];
0490                   Sum2_Pz += Object_Pz[i];
0491                   Sum2_E += Object_E[i];
0492                 }
0493
0494                 // but if a previous object was moved, add all other associated objects to the sum

0495               } else {
0496                 if (Object_Group[i] == 1) {
0497                   Sum1_Px += Object_Px[i];
0498                   Sum1_Py += Object_Py[i];
0499                   Sum1_Pz += Object_Pz[i];
0500                   Sum1_E += Object_E[i];
0501                 } else if (Object_Group[i] == 2) {
0502                   Sum2_Px += Object_Px[i];
0503                   Sum2_Py += Object_Py[i];
0504                   Sum2_Pz += Object_Pz[i];
0505                   Sum2_E += Object_E[i];
0506                 }
0507               }
0508             }
0509           }  // end loop over objects, Sum1_ and Sum2_ are now the updated hemispheres

0510         }
0511
0512       } else {
0513         cout << "Please give a valid hemisphere association method!" << endl;
0514         return 0;
0515       }
0516
0517       // recomputing the axes for next iteration

0518
0519       Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
0520       if (Axis1[3] < 0.0001) {
0521         cout << "ZERO objects in group 1! " << endl;
0522       } else {
0523         Axis1[0] = Sum1_Px / Axis1[3];
0524         Axis1[1] = Sum1_Py / Axis1[3];
0525         Axis1[2] = Sum1_Pz / Axis1[3];
0526         Axis1[4] = Sum1_E;
0527       }
0528       Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
0529       if (Axis2[3] < 0.0001) {
0530         cout << " ZERO objects in group 2! " << endl;
0531       } else {
0532         Axis2[0] = Sum2_Px / Axis2[3];
0533         Axis2[1] = Sum2_Py / Axis2[3];
0534         Axis2[2] = Sum2_Pz / Axis2[3];
0535         Axis2[4] = Sum2_E;
0536       }
0537
0538       if (dbg > 0) {
0539         cout << " Grouping = ";
0540         for (int i = 0; i < vsize; i++) {
0541           cout << "  " << Object_Group[i];
0542         }
0543         cout << endl;
0544       }
0545
0546       if (numLoop <= 1)
0547         I_Move = true;
0548
0549     }  // end of iteration

0550
0551     // associate all objects to hemispheres for method 8

0552     if (hemi_meth == 8 || hemi_meth == 9) {
0553       float sumtot = 0.;
0554       for (int i = 0; i < vsize; ++i) {
0555         if (Object_Noassoc[i] != 0)
0556           continue;
0557         if (hemi_meth == 8) {
0558           sumtot += Object_E[i];
0559         } else
0560           sumtot += Object_Pt[i];
0561       }
0562       float sumdiff = fabs(sumtot - 2 * Object_E[0]);
0563       float sum1strt = 0.;
0564       int ibest = 0, jbest = 0;
0565
0566       // start by choosing object 0 in hemi 1, all others in hemi 2

0567       // then add each of the other objects one by one to hemi 1

0568       // then add object 1 to hemi one and add each of the others in turn, etc

0569       if (vsize > 2) {
0570         for (int i = 0; i < vsize - 1; ++i) {
0571           if (Object_Noassoc[i] != 0)
0572             continue;
0573           if (hemi_meth == 8) {
0574             sum1strt += Object_E[i];
0575           } else {
0576             sum1strt += Object_Pt[i];
0577           }
0578           for (int j = i + 1; j < vsize; ++j) {
0579             if (Object_Noassoc[j] != 0)
0580               continue;
0581             float sum1_E = 0;
0582             if (hemi_meth == 8) {
0583               sum1_E = sum1strt + Object_E[j];
0584             } else {
0585               sum1_E = sum1strt + Object_Pt[j];
0586             }
0587             float sum2_E = sumtot - sum1_E;
0588             if (sumdiff >= fabs(sum1_E - sum2_E)) {
0589               sumdiff = fabs(sum1_E - sum2_E);
0590               ibest = i;
0591               jbest = j;
0592             }
0593           }
0594         }
0595       }
0596
0597       // then store the best combination into the hemisphere axes

0598       float Sum1_Px = 0, Sum1_Py = 0, Sum1_Pz = 0, Sum1_E = 0;
0599       float Sum2_Px = 0, Sum2_Py = 0, Sum2_Pz = 0, Sum2_E = 0;
0600       for (int i = 0; i < vsize; ++i) {
0601         if (Object_Noassoc[i] != 0)
0602           Object_Group[i] = 0;
0603         else if (i <= ibest || i == jbest) {
0604           Sum1_Px += Object_Px[i];
0605           Sum1_Py += Object_Py[i];
0606           Sum1_Pz += Object_Pz[i];
0607           Sum1_E += Object_E[i];
0608           Object_Group[i] = 1;
0609         } else {
0610           Sum2_Px += Object_Px[i];
0611           Sum2_Py += Object_Py[i];
0612           Sum2_Pz += Object_Pz[i];
0613           Sum2_E += Object_E[i];
0614           Object_Group[i] = 2;
0615         }
0616       }
0617       Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
0618       Axis1[0] = Sum1_Px / Axis1[3];
0619       Axis1[1] = Sum1_Py / Axis1[3];
0620       Axis1[2] = Sum1_Pz / Axis1[3];
0621       Axis1[4] = Sum1_E;
0622       Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
0623       Axis2[0] = Sum2_Px / Axis2[3];
0624       Axis2[1] = Sum2_Py / Axis2[3];
0625       Axis2[2] = Sum2_Pz / Axis2[3];
0626       Axis2[4] = Sum2_E;
0627     }
0628
0629     status = 1;
0630     return 1;
0631   }
0632
0633   int Hemisphere::RejectISR() {
0634     // tries to avoid including ISR jets into hemispheres

0635     //

0636
0637     // iterate to remove all ISR objects from the hemispheres

0638     //   until no ISR objects are found

0639     //   cout << " entered RejectISR() with rejectISRDR = " << rejectISRDR << endl;

0640     bool I_Move = true;
0641     while (I_Move) {
0642       I_Move = false;
0643       float valmax = 0.;
0644       int imax = -1;
0645
0646       // start by making a hemisphere reconstruction

0647       int hemiOK = Reconstruct();
0648       if (hemiOK == 0) {
0649         return 0;
0650       }
0651
0652       // convert the axes into Px, Py, Pz, E vectors

0653       float newAxis1_Px = Axis1[0] * Axis1[3];
0654       float newAxis1_Py = Axis1[1] * Axis1[3];
0655       float newAxis1_Pz = Axis1[2] * Axis1[3];
0656       //float newAxis1_E = Axis1[4];

0657       float newAxis2_Px = Axis2[0] * Axis2[3];
0658       float newAxis2_Py = Axis2[1] * Axis2[3];
0659       float newAxis2_Pz = Axis2[2] * Axis2[3];
0660       //float newAxis2_E = Axis2[4];

0661
0662       // loop over all objects associated to a hemisphere

0663       int vsize = (int)Object_Px.size();
0664       for (int i = 0; i < vsize; ++i) {
0665         if (Object_Group[i] == 1 || Object_Group[i] == 2) {
0666           //         cout << "  Object = " << i << ", Object_Group = " << Object_Group[i] << endl;

0667
0668           // collect the hemisphere data

0669           float newPx = 0.;
0670           float newPy = 0.;
0671           float newPz = 0.;
0672           //float newE = 0.;

0673           if (Object_Group[i] == 1) {
0674             newPx = newAxis1_Px;
0675             newPy = newAxis1_Py;
0676             newPz = newAxis1_Pz;
0677             //newE = newAxis1_E;

0678           } else if (Object_Group[i] == 2) {
0679             newPx = newAxis2_Px;
0680             newPy = newAxis2_Py;
0681             newPz = newAxis2_Pz;
0682             //newE = newAxis2_E;

0683           }
0684
0685           // compute the quantities to test whether the object is ISR

0686           float ptHemi = 0.;
0687           float hemiEta = 0.;
0688           float hemiPhi = 0.;
0689           if (rejectISRPt == 1) {
0690             float plHemi = (Object_Px[i] * newPx + Object_Py[i] * newPy + Object_Pz[i] * newPz) /
0691                            sqrt(newPx * newPx + newPy * newPy + newPz * newPz);
0692             float pObjsq = Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i] + Object_Pz[i] * Object_Pz[i];
0693             ptHemi = sqrt(pObjsq - plHemi * plHemi);
0694             if (ptHemi > valmax) {
0695               valmax = ptHemi;
0696               imax = i;
0697             }
0698           } else if (rejectISRDR == 1) {
0699             float theta = 1.570796327;
0700             // compute the new hemisphere eta, phi and DeltaR

0701             float pdiff = fabs(newPx - Object_Px[i]) + fabs(newPy - Object_Py[i]) + fabs(newPz - Object_Pz[i]);
0702             if (pdiff > 0.001) {
0703               if (fabs(newPz) > 0.001) {
0704                 theta = atan(sqrt(newPx * newPx + newPy * newPy) / newPz);
0705               }
0706               if (theta < 0.) {
0707                 theta = theta + 3.141592654;
0708               }
0709               hemiEta = -log(tan(0.5 * theta));
0710               hemiPhi = atan2(newPy, newPx);
0711               //                        float DeltaR = sqrt((Object_Eta[i] - hemiEta)*(Object_Eta[i] - hemiEta)

0712               //                            + (Util::DeltaPhi(Object_Phi[i], hemiPhi))*(Util::DeltaPhi(Object_Phi[i], hemiPhi)) );

0713               float DeltaR = sqrt((Object_Eta[i] - hemiEta) * (Object_Eta[i] - hemiEta) +
0714                                   (deltaPhi(Object_Phi[i], hemiPhi)) * (deltaPhi(Object_Phi[i], hemiPhi)));
0715               //             cout << "  Object = " << i << ", DeltaR = " << DeltaR << endl;

0716               if (DeltaR > valmax) {
0717                 valmax = DeltaR;
0718                 imax = i;
0719               }
0720             }
0721           }
0722         }
0723       }  // end loop over objects

0724
0725       // verify whether the ISR tests are fulfilled

0726       if (imax < 0) {
0727         I_Move = false;
0728       } else if (rejectISRPt == 1 && valmax > rejectISRPtmax) {
0729         SetNoAssoc(imax);
0730         I_Move = true;
0731       } else if (rejectISRDR == 1 && valmax > rejectISRDRmax) {
0732         SetNoAssoc(imax);
0733         I_Move = true;
0734       }
0735
0736     }  // end iteration over ISR objects

0737
0738     status = 1;
0739     return 1;
0740   }
0741
0742 }  // namespace heppy