heppy::Hemisphere Class Reference

#include <Hemisphere.h>

Public Member Functions
void	ClearAllNoLists ()
std::vector< float >	getAxis1 ()
std::vector< float >	getAxis2 ()
std::vector< int >	getGrouping ()
int	GetNumLoop ()
	Hemisphere ()
	Hemisphere (std::vector< float > Px_vector, std::vector< float > Py_vector, std::vector< float > Pz_vector, std::vector< float > E_vector, int seed_method, int hemisphere_association_method)
	Hemisphere (std::vector< float > Px_vector, std::vector< float > Py_vector, std::vector< float > Pz_vector, std::vector< float > E_vector)
void	RejectISRDRmax (float drmax)
void	RejectISRPtmax (float ptmax)
void	SetDebug (int debug)
void	SetDRminSeed1 (float rmin)
void	SetMethod (int seed_method, int hemisphere_association_method)
void	SetnItermax (int niter)
void	SetNoAssoc (int object_number)
void	SetNoSeed (int object_number)
	~Hemisphere ()

Private Member Functions
int	Reconstruct ()
int	RejectISR ()

Private Attributes
std::vector< float >	Axis1
std::vector< float >	Axis2
int	dbg
float	dRminSeed1
int	hemi_meth
int	nItermax
int	numLoop
std::vector< float >	Object_E
std::vector< float >	Object_Eta
std::vector< int >	Object_Group
std::vector< int >	Object_Noassoc
std::vector< int >	Object_Noseed
std::vector< float >	Object_P
std::vector< float >	Object_Phi
std::vector< float >	Object_Pt
std::vector< float >	Object_Px
std::vector< float >	Object_Py
std::vector< float >	Object_Pz
int	rejectISR
int	rejectISRDR
float	rejectISRDRmax
int	rejectISRPt
float	rejectISRPtmax
int	seed_meth
int	status

Detailed Description

Definition at line 26 of file Hemisphere.h.

Constructor & Destructor Documentation

Hemisphere() [1/3]

heppy::Hemisphere::Hemisphere ( )

inline

Definition at line 71 of file Hemisphere.h.

{};

Hemisphere() [2/3]

heppy::Hemisphere::Hemisphere	(	std::vector< float >	Px_vector,
		std::vector< float >	Py_vector,
		std::vector< float >	Pz_vector,
		std::vector< float >	E_vector,
		int	seed_method,
		int	hemisphere_association_method
	)

Definition at line 14 of file Hemisphere.cc.

References mps_fire::i, createfilelist::int, numLoop, Object_Noassoc, Object_Noseed, and Object_Px.

      : Object_Px(Px_vector),
        Object_Py(Py_vector),
        Object_Pz(Pz_vector),
        Object_E(E_vector),
        seed_meth(seed_method),
        hemi_meth(hemisphere_association_method),
        status(0),
        dRminSeed1(0.5),
        nItermax(100),
        rejectISR(0),
        rejectISRPt(0),
        rejectISRPtmax(10000.),
        rejectISRDR(0),
        rejectISRDRmax(100.),
        dbg(0) {
    for (int i = 0; i < (int)Object_Px.size(); i++) {
      Object_Noseed.push_back(0);
      Object_Noassoc.push_back(0);
    }
    numLoop = 0;
  }

Hemisphere() [3/3]

heppy::Hemisphere::Hemisphere	(	std::vector< float >	Px_vector,
		std::vector< float >	Py_vector,
		std::vector< float >	Pz_vector,
		std::vector< float >	E_vector
	)

Definition at line 44 of file Hemisphere.cc.

References mps_fire::i, createfilelist::int, numLoop, Object_Noassoc, Object_Noseed, and Object_Px.

      : Object_Px(Px_vector),
        Object_Py(Py_vector),
        Object_Pz(Pz_vector),
        Object_E(E_vector),
        seed_meth(0),
        hemi_meth(0),
        status(0),
        dRminSeed1(0.5),
        nItermax(100),
        rejectISR(0),
        rejectISRPt(0),
        rejectISRPtmax(10000.),
        rejectISRDR(0),
        rejectISRDRmax(100.),
        dbg(0) {
    for (int i = 0; i < (int)Object_Px.size(); i++) {
      Object_Noseed.push_back(0);
      Object_Noassoc.push_back(0);
    }
    numLoop = 0;
  }

~Hemisphere()

heppy::Hemisphere::~Hemisphere ( )

inline

Definition at line 86 of file Hemisphere.h.

{};

Member Function Documentation

ClearAllNoLists()

void heppy::Hemisphere::ClearAllNoLists ( )

inline

Definition at line 125 of file Hemisphere.h.

References mps_fire::i, createfilelist::int, Object_Noassoc, Object_Noseed, and status.

                           {
      for (int i = 0; i < (int)Object_Noseed.size(); ++i) {
        Object_Noassoc[i] = 0;
        Object_Noseed[i] = 0;
        status = 0;
      }
    }

getAxis1()

vector< float > heppy::Hemisphere::getAxis1

(

)

Definition at line 70 of file Hemisphere.cc.

References Axis1, Reconstruct(), RejectISR(), rejectISR, and status.

                                     {
    if (status != 1) {
      if (rejectISR == 0) {
        this->Reconstruct();
      } else {
        this->RejectISR();
      }
    }
    return Axis1;
  }

getAxis2()

vector< float > heppy::Hemisphere::getAxis2

(

)

Definition at line 80 of file Hemisphere.cc.

References Axis2, Reconstruct(), RejectISR(), rejectISR, and status.

                                     {
    if (status != 1) {
      if (rejectISR == 0) {
        this->Reconstruct();
      } else {
        this->RejectISR();
      }
    }
    return Axis2;
  }

getGrouping()

vector< int > heppy::Hemisphere::getGrouping

(

)

Definition at line 91 of file Hemisphere.cc.

References Object_Group, Reconstruct(), RejectISR(), rejectISR, and status.

                                      {
    if (status != 1) {
      if (rejectISR == 0) {
        this->Reconstruct();
      } else {
        this->RejectISR();
      }
    }
    return Object_Group;
  }

GetNumLoop()

int heppy::Hemisphere::GetNumLoop ( )

inline

Definition at line 161 of file Hemisphere.h.

References numLoop.

{ return numLoop; }

Reconstruct()

int heppy::Hemisphere::Reconstruct ( )

private

Definition at line 102 of file Hemisphere.cc.

References Axis1, Axis2, gather_cfg::cout, dbg, SiPixelRawToDigiRegional_cfi::deltaPhi, PbPb_ZMuSkimMuonDPG_cff::deltaR, dRminSeed1, hemi_meth, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, dqm-mbProfile::log, nItermax, numLoop, Object_E, Object_Eta, Object_Group, Object_Noassoc, Object_Noseed, Object_P, Object_Phi, Object_Pt, Object_Px, Object_Py, Object_Pz, seed_meth, mathSSE::sqrt(), status, funct::tan(), and theta().

Referenced by getAxis1(), getAxis2(), getGrouping(), and RejectISR().

                              {
    // Performs the actual hemisphere reconstrucion
    //
    // definition of the vectors used internally:
    // Object_xxx :
    //        xxx = Px, Py, Pz, E for input values
    //        xxx = P, Pt, Eta, Phi, Group for internal use
    // Axis1 : final hemisphere axis 1
    // Axis2 : final hemisphere axis 2
    // Sum1_xxx : hemisphere 1 being updated during the association iterations
    // Sum2_xxx : hemisphere 2 being updated during the association iterations
    // NewAxis1_xxx, NewAxis1_xxx : temporary axes for calculation in association methods 2 and 3

    numLoop = 0;  // initialize numLoop for Zero
    int vsize = (int)Object_Px.size();
    if ((int)Object_Py.size() != vsize || (int)Object_Pz.size() != vsize) {
      cout << "WARNING!!!!! Input vectors have different size! Fix it!" << endl;
      return 0;
    }
    if (dbg > 0) {
      //    cout << " Hemisphere method, vsn = " << hemivsn << endl;
      cout << " Hemisphere method " << endl;
    }

    // clear some vectors if method reconstruct() is called again
    if (!Object_P.empty()) {
      Object_P.clear();
      Object_Pt.clear();
      Object_Eta.clear();
      Object_Phi.clear();
      Object_Group.clear();
      Axis1.clear();
      Axis2.clear();
    }
    // initialize the vectors
    for (int j = 0; j < vsize; ++j) {
      Object_P.push_back(0);
      Object_Pt.push_back(0);
      Object_Eta.push_back(0);
      Object_Phi.push_back(0);
      Object_Group.push_back(0);
    }
    for (int j = 0; j < 5; ++j) {
      Axis1.push_back(0);
      Axis2.push_back(0);
    }

    // compute additional quantities for vectors Object_xxx
    float theta;
    for (int i = 0; i < vsize; ++i) {
      Object_P[i] = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i] + Object_Pz[i] * Object_Pz[i]);
      if (Object_P[i] > Object_E[i] + 0.001) {
        cout << "WARNING!!!!! Object " << i << " has E = " << Object_E[i] << " less than P = " << Object_P[i]
             << " *** Fix it!" << endl;
        return 0;
      }
      Object_Pt[i] = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]);
      // protection for div by 0
      if (fabs(Object_Pz[i]) > 0.001) {
        theta = atan(sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]) / Object_Pz[i]);
      } else {
        theta = 1.570796327;
      }
      if (theta < 0.) {
        theta = theta + 3.141592654;
      }
      Object_Eta[i] = -log(tan(0.5 * theta));
      Object_Phi[i] = atan2(Object_Py[i], Object_Px[i]);
      if (dbg > 0) {
        cout << " Object " << i << " Eta = " << Object_Eta[i] << " Phi = " << Object_Phi[i] << endl;
      }
    }

    if (dbg > 0) {
      cout << endl;
      cout << " Seeding method = " << seed_meth << endl;
    }
    // I_Max and J_Max are indices of the seeds in the vectors
    int I_Max = -1;
    int J_Max = -1;

    // determine the seeds for seed method 1
    if (seed_meth == 1) {
      float P_Max = 0.;
      float DeltaRP_Max = 0.;

      // take highest momentum object as first seed
      for (int i = 0; i < vsize; ++i) {
        Object_Group[i] = 0;
        //cout << "Object_Px[i] = " << Object_Px[i] << ", Object_Py[i] = " << Object_Py[i]
        //<< ", Object_Pz[i] = " << Object_Pz[i] << "  << endl;
        if (Object_Noseed[i] == 0 && P_Max < Object_P[i]) {
          P_Max = Object_P[i];
          I_Max = i;
        }
      }
      // if 1st seed is found, save it as initial hemisphere 1 axis
      if (I_Max >= 0) {
        Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
        Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
        Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
        Axis1[3] = Object_P[I_Max];
        Axis1[4] = Object_E[I_Max];
      } else {
        // cout << " This is an empty event." << endl;
        return 0;
      }

      // take as second seed the object with largest DR*P w.r.t. the first seed
      for (int i = 0; i < vsize; ++i) {
        //          float DeltaR = sqrt((Object_Eta[i] - Object_Eta[I_Max])*(Object_Eta[i] - Object_Eta[I_Max])
        //              + (Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max]))*(Util::DeltaPhi(Object_Phi[i], Object_Phi[I_Max]))  );
        float DeltaR =
            sqrt((Object_Eta[i] - Object_Eta[I_Max]) * (Object_Eta[i] - Object_Eta[I_Max]) +
                 (deltaPhi(Object_Phi[i], Object_Phi[I_Max])) * (deltaPhi(Object_Phi[i], Object_Phi[I_Max])));
        if (Object_Noseed[i] == 0 && DeltaR > dRminSeed1) {
          float DeltaRP = DeltaR * Object_P[i];
          if (DeltaRP > DeltaRP_Max) {
            DeltaRP_Max = DeltaRP;
            J_Max = i;
          }
        }
      }
      // if 2nd seed is found, save it as initial hemisphere 2 axis
      if (J_Max >= 0) {
        Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
        Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
        Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
        Axis2[3] = Object_P[J_Max];
        Axis2[4] = Object_E[J_Max];
      } else {
        // cout << " This is a MONOJET." << endl;
        return 0;
      }
      if (dbg > 0) {
        cout << " Axis 1 is Object = " << I_Max << endl;
        cout << " Axis 2 is Object = " << J_Max << endl;
      }

      // determine the seeds for seed methods 2 and 3
    } else if (seed_meth == 2 || seed_meth == 3) {
      float Mass_Max = 0.;
      float InvariantMass = 0.;

      // maximize the invariant mass of two objects
      for (int i = 0; i < vsize; ++i) {
        Object_Group[i] = 0;
        if (Object_Noseed[i] == 0) {
          for (int j = i + 1; j < vsize; ++j) {
            if (Object_Noseed[j] == 0) {
              // either the invariant mass
              if (seed_meth == 2) {
                InvariantMass = (Object_E[i] + Object_E[j]) * (Object_E[i] + Object_E[j]) -
                                (Object_Px[i] + Object_Px[j]) * (Object_Px[i] + Object_Px[j]) -
                                (Object_Py[i] + Object_Py[j]) * (Object_Py[i] + Object_Py[j]) -
                                (Object_Pz[i] + Object_Pz[j]) * (Object_Pz[i] + Object_Pz[j]);
              }
              // or the transverse mass
              else if (seed_meth == 3) {
                float pti = sqrt(Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i]);
                float ptj = sqrt(Object_Px[j] * Object_Px[j] + Object_Py[j] * Object_Py[j]);
                InvariantMass = 2. * (pti * ptj - Object_Px[i] * Object_Px[j] - Object_Py[i] * Object_Py[j]);
              }
              if (Mass_Max < InvariantMass) {
                Mass_Max = InvariantMass;
                I_Max = i;
                J_Max = j;
              }
            }
          }
        }
      }

      // if both seeds are found, save them as initial hemisphere axes
      if (J_Max > 0) {
        Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
        Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
        Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
        Axis1[3] = Object_P[I_Max];
        Axis1[4] = Object_E[I_Max];

        Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
        Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
        Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
        Axis2[3] = Object_P[J_Max];
        Axis2[4] = Object_E[J_Max];
      } else {
        // cout << " This is a MONOJET." << endl;
        return 0;
      }
      if (dbg > 0) {
        cout << " Axis 1 is Object = " << I_Max << endl;
        cout << " Axis 2 is Object = " << J_Max << endl;
      }

    } else if (seed_meth == 4) {
      float P_Max1 = 0.;
      float P_Max2 = 0.;

      // take largest Pt object as first seed
      for (int i = 0; i < vsize; ++i) {
        Object_Group[i] = 0;
        if (Object_Noseed[i] == 0 && P_Max1 < Object_Pt[i]) {
          P_Max1 = Object_Pt[i];
          I_Max = i;
        }
      }
      if (I_Max < 0)
        return 0;

      // take second largest Pt object as second seed, but require dR(seed1, seed2) > dRminSeed1
      for (int i = 0; i < vsize; ++i) {
        if (i == I_Max)
          continue;
        //          float DeltaR = Util::GetDeltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);
        float DeltaR = deltaR(Object_Eta[i], Object_Eta[I_Max], Object_Phi[i], Object_Phi[I_Max]);
        if (Object_Noseed[i] == 0 && P_Max2 < Object_Pt[i] && DeltaR > dRminSeed1) {
          P_Max2 = Object_Pt[i];
          J_Max = i;
        }
      }
      if (J_Max < 0)
        return 0;

      // save first seed as initial hemisphere 1 axis
      if (I_Max >= 0) {
        Axis1[0] = Object_Px[I_Max] / Object_P[I_Max];
        Axis1[1] = Object_Py[I_Max] / Object_P[I_Max];
        Axis1[2] = Object_Pz[I_Max] / Object_P[I_Max];
        Axis1[3] = Object_P[I_Max];
        Axis1[4] = Object_E[I_Max];
      }

      // save second seed as initial hemisphere 2 axis
      if (J_Max >= 0) {
        Axis2[0] = Object_Px[J_Max] / Object_P[J_Max];
        Axis2[1] = Object_Py[J_Max] / Object_P[J_Max];
        Axis2[2] = Object_Pz[J_Max] / Object_P[J_Max];
        Axis2[3] = Object_P[J_Max];
        Axis2[4] = Object_E[J_Max];
      }

      if (dbg > 0) {
        cout << " Axis 1 is Object = " << I_Max << " with Pt " << Object_Pt[I_Max] << endl;
        cout << " Axis 2 is Object = " << J_Max << " with Pt " << Object_Pt[J_Max] << endl;
      }

    } else if (!(seed_meth == 0 && (hemi_meth == 8 || hemi_meth == 9))) {
      cout << "Please give a valid seeding method!" << endl;
      return 0;
    }

    // seeding done
    // now do the hemisphere association

    if (dbg > 0) {
      cout << endl;
      cout << " Association method = " << hemi_meth << endl;
    }

    bool I_Move = true;

    // iterate to associate all objects to hemispheres (methods 1 to 3 only)
    //   until no objects are moved from one to the other hemisphere
    //   or the maximum number of iterations is reached
    while (I_Move && (numLoop < nItermax) && hemi_meth != 8 && hemi_meth != 9) {
      I_Move = false;
      numLoop++;
      if (dbg > 0) {
        cout << " Iteration = " << numLoop << endl;
      }
      if (numLoop == nItermax - 1) {
        cout << " Hemishpere: warning - reaching max number of iterations " << endl;
      }

      // initialize the current sums of Px, Py, Pz, E for the two hemispheres
      float Sum1_Px = 0.;
      float Sum1_Py = 0.;
      float Sum1_Pz = 0.;
      float Sum1_E = 0.;
      float Sum2_Px = 0.;
      float Sum2_Py = 0.;
      float Sum2_Pz = 0.;
      float Sum2_E = 0.;

      // associate the objects for method 1
      if (hemi_meth == 1) {
        for (int i = 0; i < vsize; ++i) {
          if (Object_Noassoc[i] == 0) {
            float P_Long1 = Object_Px[i] * Axis1[0] + Object_Py[i] * Axis1[1] + Object_Pz[i] * Axis1[2];
            float P_Long2 = Object_Px[i] * Axis2[0] + Object_Py[i] * Axis2[1] + Object_Pz[i] * Axis2[2];
            if (P_Long1 >= P_Long2) {
              if (Object_Group[i] != 1) {
                I_Move = true;
              }
              Object_Group[i] = 1;
              Sum1_Px += Object_Px[i];
              Sum1_Py += Object_Py[i];
              Sum1_Pz += Object_Pz[i];
              Sum1_E += Object_E[i];
            } else {
              if (Object_Group[i] != 2) {
                I_Move = true;
              }
              Object_Group[i] = 2;
              Sum2_Px += Object_Px[i];
              Sum2_Py += Object_Py[i];
              Sum2_Pz += Object_Pz[i];
              Sum2_E += Object_E[i];
            }
          }
        }

        // associate the objects for methods 2 and 3
      } else if (hemi_meth == 2 || hemi_meth == 3) {
        for (int i = 0; i < vsize; ++i) {
          // add the seeds to the sums, as they remain fixed
          if (i == I_Max) {
            Object_Group[i] = 1;
            Sum1_Px += Object_Px[i];
            Sum1_Py += Object_Py[i];
            Sum1_Pz += Object_Pz[i];
            Sum1_E += Object_E[i];
          } else if (i == J_Max) {
            Object_Group[i] = 2;
            Sum2_Px += Object_Px[i];
            Sum2_Py += Object_Py[i];
            Sum2_Pz += Object_Pz[i];
            Sum2_E += Object_E[i];

            // for the other objects
          } else {
            if (Object_Noassoc[i] == 0) {
              // only 1 object maximum is moved in a given iteration
              if (!I_Move) {
                // initialize the new hemispheres as the current ones
                float NewAxis1_Px = Axis1[0] * Axis1[3];
                float NewAxis1_Py = Axis1[1] * Axis1[3];
                float NewAxis1_Pz = Axis1[2] * Axis1[3];
                float NewAxis1_E = Axis1[4];
                float NewAxis2_Px = Axis2[0] * Axis2[3];
                float NewAxis2_Py = Axis2[1] * Axis2[3];
                float NewAxis2_Pz = Axis2[2] * Axis2[3];
                float NewAxis2_E = Axis2[4];
                // subtract the object from its hemisphere
                if (Object_Group[i] == 1) {
                  NewAxis1_Px = NewAxis1_Px - Object_Px[i];
                  NewAxis1_Py = NewAxis1_Py - Object_Py[i];
                  NewAxis1_Pz = NewAxis1_Pz - Object_Pz[i];
                  NewAxis1_E = NewAxis1_E - Object_E[i];
                } else if (Object_Group[i] == 2) {
                  NewAxis2_Px = NewAxis2_Px - Object_Px[i];
                  NewAxis2_Py = NewAxis2_Py - Object_Py[i];
                  NewAxis2_Pz = NewAxis2_Pz - Object_Pz[i];
                  NewAxis2_E = NewAxis2_E - Object_E[i];
                }

                // compute the invariant mass squared with each hemisphere (method 2)
                float mass1 = NewAxis1_E -
                              ((Object_Px[i] * NewAxis1_Px + Object_Py[i] * NewAxis1_Py + Object_Pz[i] * NewAxis1_Pz) /
                               Object_P[i]);
                float mass2 = NewAxis2_E -
                              ((Object_Px[i] * NewAxis2_Px + Object_Py[i] * NewAxis2_Py + Object_Pz[i] * NewAxis2_Pz) /
                               Object_P[i]);
                // or the Lund distance (method 3)
                if (hemi_meth == 3) {
                  mass1 *= NewAxis1_E / ((NewAxis1_E + Object_E[i]) * (NewAxis1_E + Object_E[i]));
                  mass2 *= NewAxis2_E / ((NewAxis2_E + Object_E[i]) * (NewAxis2_E + Object_E[i]));
                }
                // and associate the object to the best hemisphere and add it to the sum
                if (mass1 < mass2) {
                  //if (Object_Group[i] != 1){
                  if (Object_Group[i] != 1 && Object_Group[i] != 0) {
                    I_Move = true;
                  }
                  Object_Group[i] = 1;
                  Sum1_Px += Object_Px[i];
                  Sum1_Py += Object_Py[i];
                  Sum1_Pz += Object_Pz[i];
                  Sum1_E += Object_E[i];
                } else {
                  //if (Object_Group[i] != 2){
                  if (Object_Group[i] != 2 && Object_Group[i] != 0) {
                    I_Move = true;
                  }
                  Object_Group[i] = 2;
                  Sum2_Px += Object_Px[i];
                  Sum2_Py += Object_Py[i];
                  Sum2_Pz += Object_Pz[i];
                  Sum2_E += Object_E[i];
                }

                // but if a previous object was moved, add all other associated objects to the sum
              } else {
                if (Object_Group[i] == 1) {
                  Sum1_Px += Object_Px[i];
                  Sum1_Py += Object_Py[i];
                  Sum1_Pz += Object_Pz[i];
                  Sum1_E += Object_E[i];
                } else if (Object_Group[i] == 2) {
                  Sum2_Px += Object_Px[i];
                  Sum2_Py += Object_Py[i];
                  Sum2_Pz += Object_Pz[i];
                  Sum2_E += Object_E[i];
                }
              }
            }
          }  // end loop over objects, Sum1_ and Sum2_ are now the updated hemispheres
        }

      } else {
        cout << "Please give a valid hemisphere association method!" << endl;
        return 0;
      }

      // recomputing the axes for next iteration

      Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
      if (Axis1[3] < 0.0001) {
        cout << "ZERO objects in group 1! " << endl;
      } else {
        Axis1[0] = Sum1_Px / Axis1[3];
        Axis1[1] = Sum1_Py / Axis1[3];
        Axis1[2] = Sum1_Pz / Axis1[3];
        Axis1[4] = Sum1_E;
      }
      Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
      if (Axis2[3] < 0.0001) {
        cout << " ZERO objects in group 2! " << endl;
      } else {
        Axis2[0] = Sum2_Px / Axis2[3];
        Axis2[1] = Sum2_Py / Axis2[3];
        Axis2[2] = Sum2_Pz / Axis2[3];
        Axis2[4] = Sum2_E;
      }

      if (dbg > 0) {
        cout << " Grouping = ";
        for (int i = 0; i < vsize; i++) {
          cout << "  " << Object_Group[i];
        }
        cout << endl;
      }

      if (numLoop <= 1)
        I_Move = true;

    }  // end of iteration

    // associate all objects to hemispheres for method 8
    if (hemi_meth == 8 || hemi_meth == 9) {
      float sumtot = 0.;
      for (int i = 0; i < vsize; ++i) {
        if (Object_Noassoc[i] != 0)
          continue;
        if (hemi_meth == 8) {
          sumtot += Object_E[i];
        } else
          sumtot += Object_Pt[i];
      }
      float sumdiff = fabs(sumtot - 2 * Object_E[0]);
      float sum1strt = 0.;
      int ibest = 0, jbest = 0;

      // start by choosing object 0 in hemi 1, all others in hemi 2
      // then add each of the other objects one by one to hemi 1
      // then add object 1 to hemi one and add each of the others in turn, etc
      if (vsize > 2) {
        for (int i = 0; i < vsize - 1; ++i) {
          if (Object_Noassoc[i] != 0)
            continue;
          if (hemi_meth == 8) {
            sum1strt += Object_E[i];
          } else {
            sum1strt += Object_Pt[i];
          }
          for (int j = i + 1; j < vsize; ++j) {
            if (Object_Noassoc[j] != 0)
              continue;
            float sum1_E = 0;
            if (hemi_meth == 8) {
              sum1_E = sum1strt + Object_E[j];
            } else {
              sum1_E = sum1strt + Object_Pt[j];
            }
            float sum2_E = sumtot - sum1_E;
            if (sumdiff >= fabs(sum1_E - sum2_E)) {
              sumdiff = fabs(sum1_E - sum2_E);
              ibest = i;
              jbest = j;
            }
          }
        }
      }

      // then store the best combination into the hemisphere axes
      float Sum1_Px = 0, Sum1_Py = 0, Sum1_Pz = 0, Sum1_E = 0;
      float Sum2_Px = 0, Sum2_Py = 0, Sum2_Pz = 0, Sum2_E = 0;
      for (int i = 0; i < vsize; ++i) {
        if (Object_Noassoc[i] != 0)
          Object_Group[i] = 0;
        else if (i <= ibest || i == jbest) {
          Sum1_Px += Object_Px[i];
          Sum1_Py += Object_Py[i];
          Sum1_Pz += Object_Pz[i];
          Sum1_E += Object_E[i];
          Object_Group[i] = 1;
        } else {
          Sum2_Px += Object_Px[i];
          Sum2_Py += Object_Py[i];
          Sum2_Pz += Object_Pz[i];
          Sum2_E += Object_E[i];
          Object_Group[i] = 2;
        }
      }
      Axis1[3] = sqrt(Sum1_Px * Sum1_Px + Sum1_Py * Sum1_Py + Sum1_Pz * Sum1_Pz);
      Axis1[0] = Sum1_Px / Axis1[3];
      Axis1[1] = Sum1_Py / Axis1[3];
      Axis1[2] = Sum1_Pz / Axis1[3];
      Axis1[4] = Sum1_E;
      Axis2[3] = sqrt(Sum2_Px * Sum2_Px + Sum2_Py * Sum2_Py + Sum2_Pz * Sum2_Pz);
      Axis2[0] = Sum2_Px / Axis2[3];
      Axis2[1] = Sum2_Py / Axis2[3];
      Axis2[2] = Sum2_Pz / Axis2[3];
      Axis2[4] = Sum2_E;
    }

    status = 1;
    return 1;
  }

RejectISR()

int heppy::Hemisphere::RejectISR ( )

private

Definition at line 633 of file Hemisphere.cc.

References Axis1, Axis2, SiPixelRawToDigiRegional_cfi::deltaPhi, electronAnalyzer_cfi::DeltaR, mps_fire::i, createfilelist::int, dqm-mbProfile::log, Object_Eta, Object_Group, Object_Phi, Object_Px, Object_Py, Object_Pz, Reconstruct(), rejectISRDR, rejectISRDRmax, rejectISRPt, rejectISRPtmax, SetNoAssoc(), mathSSE::sqrt(), status, funct::tan(), and theta().

Referenced by getAxis1(), getAxis2(), and getGrouping().

                            {
    // tries to avoid including ISR jets into hemispheres
    //

    // iterate to remove all ISR objects from the hemispheres
    //   until no ISR objects are found
    //   cout << " entered RejectISR() with rejectISRDR = " << rejectISRDR << endl;
    bool I_Move = true;
    while (I_Move) {
      I_Move = false;
      float valmax = 0.;
      int imax = -1;

      // start by making a hemisphere reconstruction
      int hemiOK = Reconstruct();
      if (hemiOK == 0) {
        return 0;
      }

      // convert the axes into Px, Py, Pz, E vectors
      float newAxis1_Px = Axis1[0] * Axis1[3];
      float newAxis1_Py = Axis1[1] * Axis1[3];
      float newAxis1_Pz = Axis1[2] * Axis1[3];
      //float newAxis1_E = Axis1[4];
      float newAxis2_Px = Axis2[0] * Axis2[3];
      float newAxis2_Py = Axis2[1] * Axis2[3];
      float newAxis2_Pz = Axis2[2] * Axis2[3];
      //float newAxis2_E = Axis2[4];

      // loop over all objects associated to a hemisphere
      int vsize = (int)Object_Px.size();
      for (int i = 0; i < vsize; ++i) {
        if (Object_Group[i] == 1 || Object_Group[i] == 2) {
          //         cout << "  Object = " << i << ", Object_Group = " << Object_Group[i] << endl;

          // collect the hemisphere data
          float newPx = 0.;
          float newPy = 0.;
          float newPz = 0.;
          //float newE = 0.;
          if (Object_Group[i] == 1) {
            newPx = newAxis1_Px;
            newPy = newAxis1_Py;
            newPz = newAxis1_Pz;
            //newE = newAxis1_E;
          } else if (Object_Group[i] == 2) {
            newPx = newAxis2_Px;
            newPy = newAxis2_Py;
            newPz = newAxis2_Pz;
            //newE = newAxis2_E;
          }

          // compute the quantities to test whether the object is ISR
          float ptHemi = 0.;
          float hemiEta = 0.;
          float hemiPhi = 0.;
          if (rejectISRPt == 1) {
            float plHemi = (Object_Px[i] * newPx + Object_Py[i] * newPy + Object_Pz[i] * newPz) /
                           sqrt(newPx * newPx + newPy * newPy + newPz * newPz);
            float pObjsq = Object_Px[i] * Object_Px[i] + Object_Py[i] * Object_Py[i] + Object_Pz[i] * Object_Pz[i];
            ptHemi = sqrt(pObjsq - plHemi * plHemi);
            if (ptHemi > valmax) {
              valmax = ptHemi;
              imax = i;
            }
          } else if (rejectISRDR == 1) {
            float theta = 1.570796327;
            // compute the new hemisphere eta, phi and DeltaR
            float pdiff = fabs(newPx - Object_Px[i]) + fabs(newPy - Object_Py[i]) + fabs(newPz - Object_Pz[i]);
            if (pdiff > 0.001) {
              if (fabs(newPz) > 0.001) {
                theta = atan(sqrt(newPx * newPx + newPy * newPy) / newPz);
              }
              if (theta < 0.) {
                theta = theta + 3.141592654;
              }
              hemiEta = -log(tan(0.5 * theta));
              hemiPhi = atan2(newPy, newPx);
              //                        float DeltaR = sqrt((Object_Eta[i] - hemiEta)*(Object_Eta[i] - hemiEta)
              //                            + (Util::DeltaPhi(Object_Phi[i], hemiPhi))*(Util::DeltaPhi(Object_Phi[i], hemiPhi)) );
              float DeltaR = sqrt((Object_Eta[i] - hemiEta) * (Object_Eta[i] - hemiEta) +
                                  (deltaPhi(Object_Phi[i], hemiPhi)) * (deltaPhi(Object_Phi[i], hemiPhi)));
              //             cout << "  Object = " << i << ", DeltaR = " << DeltaR << endl;
              if (DeltaR > valmax) {
                valmax = DeltaR;
                imax = i;
              }
            }
          }
        }
      }  // end loop over objects

      // verify whether the ISR tests are fulfilled
      if (imax < 0) {
        I_Move = false;
      } else if (rejectISRPt == 1 && valmax > rejectISRPtmax) {
        SetNoAssoc(imax);
        I_Move = true;
      } else if (rejectISRDR == 1 && valmax > rejectISRDRmax) {
        SetNoAssoc(imax);
        I_Move = true;
      }

    }  // end iteration over ISR objects

    status = 1;
    return 1;
  }

RejectISRDRmax()

void heppy::Hemisphere::RejectISRDRmax ( float  drmax )

inline

Definition at line 151 of file Hemisphere.h.

References rejectISR, rejectISRDR, rejectISRDRmax, rejectISRPt, and status.

                                     {
      rejectISRDRmax = drmax;
      rejectISRDR = 1;
      rejectISRPt = 0;
      rejectISR = 1;
      status = 0;
    }

RejectISRPtmax()

void heppy::Hemisphere::RejectISRPtmax ( float  ptmax )

inline

Definition at line 142 of file Hemisphere.h.

References muonTiming_cfi::ptmax, rejectISR, rejectISRDR, rejectISRPt, rejectISRPtmax, and status.

                                     {
      rejectISRPtmax = ptmax;
      rejectISRPt = 1;
      rejectISRDR = 0;
      rejectISR = 1;
      status = 0;
    }

SetDebug()

void heppy::Hemisphere::SetDebug ( int  debug )

inline

Definition at line 160 of file Hemisphere.h.

References dbg, and debug.

{ dbg = debug; }

SetDRminSeed1()

void heppy::Hemisphere::SetDRminSeed1 ( float  rmin )

inline

Definition at line 135 of file Hemisphere.h.

References dRminSeed1.

{ dRminSeed1 = rmin; }

SetMethod()

void heppy::Hemisphere::SetMethod ( int  seed_method, int  hemisphere_association_method  )

inline

Definition at line 102 of file Hemisphere.h.

References hemi_meth, seed_meth, and status.

                                                                       {
      seed_meth = seed_method;
      hemi_meth = hemisphere_association_method;
      status = 0;
    }

SetnItermax()

void heppy::Hemisphere::SetnItermax ( int  niter )

inline

Definition at line 139 of file Hemisphere.h.

References nItermax.

{ nItermax = niter; }

SetNoAssoc()

void heppy::Hemisphere::SetNoAssoc ( int  object_number )

inline

Definition at line 117 of file Hemisphere.h.

References Object_Noassoc, Object_Noseed, and status.

Referenced by RejectISR().

                                       {
      Object_Noassoc[object_number] = 1;
      Object_Noseed[object_number] = 1;
      status = 0;
    }

SetNoSeed()

void heppy::Hemisphere::SetNoSeed ( int  object_number )

inline

Definition at line 110 of file Hemisphere.h.

References Object_Noseed, and status.

                                      {
      Object_Noseed[object_number] = 1;
      status = 0;
    }

Member Data Documentation

Axis1

std::vector<float> heppy::Hemisphere::Axis1

private

Definition at line 180 of file Hemisphere.h.

Referenced by getAxis1(), Reconstruct(), and RejectISR().

Axis2

std::vector<float> heppy::Hemisphere::Axis2

private

Definition at line 181 of file Hemisphere.h.

Referenced by getAxis2(), Reconstruct(), and RejectISR().

dbg

int heppy::Hemisphere::dbg

private

Definition at line 194 of file Hemisphere.h.

Referenced by Reconstruct(), and SetDebug().

dRminSeed1

float heppy::Hemisphere::dRminSeed1

private

Definition at line 187 of file Hemisphere.h.

Referenced by Reconstruct(), and SetDRminSeed1().

hemi_meth

int heppy::Hemisphere::hemi_meth

private

Definition at line 185 of file Hemisphere.h.

Referenced by Reconstruct(), and SetMethod().

nItermax

int heppy::Hemisphere::nItermax

private

Definition at line 188 of file Hemisphere.h.

Referenced by Reconstruct(), and SetnItermax().

numLoop

int heppy::Hemisphere::numLoop

private

Definition at line 195 of file Hemisphere.h.

Referenced by GetNumLoop(), Hemisphere(), and Reconstruct().

Object_E

std::vector<float> heppy::Hemisphere::Object_E

private

Definition at line 173 of file Hemisphere.h.

Referenced by Reconstruct().

Object_Eta

std::vector<float> heppy::Hemisphere::Object_Eta

private

Definition at line 175 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

Object_Group

std::vector<int> heppy::Hemisphere::Object_Group

private

Definition at line 176 of file Hemisphere.h.

Referenced by getGrouping(), Reconstruct(), and RejectISR().

Object_Noassoc

std::vector<int> heppy::Hemisphere::Object_Noassoc

private

Definition at line 178 of file Hemisphere.h.

Referenced by ClearAllNoLists(), Hemisphere(), Reconstruct(), and SetNoAssoc().

Object_Noseed

std::vector<int> heppy::Hemisphere::Object_Noseed

private

Definition at line 177 of file Hemisphere.h.

Referenced by ClearAllNoLists(), Hemisphere(), Reconstruct(), SetNoAssoc(), and SetNoSeed().

Object_P

std::vector<float> heppy::Hemisphere::Object_P

private

Definition at line 171 of file Hemisphere.h.

Referenced by Reconstruct().

Object_Phi

std::vector<float> heppy::Hemisphere::Object_Phi

private

Definition at line 174 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

Object_Pt

std::vector<float> heppy::Hemisphere::Object_Pt

private

Definition at line 172 of file Hemisphere.h.

Referenced by Reconstruct().

Object_Px

std::vector<float> heppy::Hemisphere::Object_Px

private

Definition at line 168 of file Hemisphere.h.

Referenced by Hemisphere(), Reconstruct(), and RejectISR().

Object_Py

std::vector<float> heppy::Hemisphere::Object_Py

private

Definition at line 169 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

Object_Pz

std::vector<float> heppy::Hemisphere::Object_Pz

private

Definition at line 170 of file Hemisphere.h.

Referenced by Reconstruct(), and RejectISR().

rejectISR

int heppy::Hemisphere::rejectISR

private

Definition at line 189 of file Hemisphere.h.

Referenced by getAxis1(), getAxis2(), getGrouping(), RejectISRDRmax(), and RejectISRPtmax().

rejectISRDR

int heppy::Hemisphere::rejectISRDR

private

Definition at line 192 of file Hemisphere.h.

Referenced by RejectISR(), RejectISRDRmax(), and RejectISRPtmax().

rejectISRDRmax

float heppy::Hemisphere::rejectISRDRmax

private

Definition at line 193 of file Hemisphere.h.

Referenced by RejectISR(), and RejectISRDRmax().

rejectISRPt

int heppy::Hemisphere::rejectISRPt

private

Definition at line 190 of file Hemisphere.h.

Referenced by RejectISR(), RejectISRDRmax(), and RejectISRPtmax().

rejectISRPtmax

float heppy::Hemisphere::rejectISRPtmax

private

Definition at line 191 of file Hemisphere.h.

Referenced by RejectISR(), and RejectISRPtmax().

seed_meth

int heppy::Hemisphere::seed_meth

private

Definition at line 184 of file Hemisphere.h.

Referenced by Reconstruct(), and SetMethod().

status

int heppy::Hemisphere::status

private

Definition at line 186 of file Hemisphere.h.

Referenced by dirstructure.Comparison::__make_image(), dirstructure.Comparison::__repr__(), ClearAllNoLists(), getAxis1(), getAxis2(), getGrouping(), Reconstruct(), RejectISR(), RejectISRDRmax(), RejectISRPtmax(), SetMethod(), SetNoAssoc(), and SetNoSeed().