siscone_parton.h

Go to the documentation of this file.

/*
 *  Copyright (C) 2011 Stefan Weinzierl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#ifndef __SISCONE_PARTON_H__
#define __SISCONE_PARTON_H__

namespace siscone_jet_algorithm {

  const double Pi = 3.14159265358979;

  template<class T> class siscone_parton
    {
      // ctors
    public:
      siscone_parton(int n_particles_max);

      // member functions
    public:
      int compute_jets(const std::vector<T> & particles, double R, double f, int n_pass_max=0, double ptmin=0.0);

      void set_flag_total_pt_is_zero(bool flag);

      // member variables 
    public:
      std::vector<unsigned> final_jets;

      std::vector<unsigned> protojets;

      std::vector<T> protojets_axis;
      
      std::vector<double> protojets_pt;

      std::vector<double> protojets_pttilde;

      std::vector<double> protojets_azimuthal_angle;

      std::vector<double> protojets_rapidity;

      std::vector<bool> flag_rapidity;

      bool flag_total_pt_is_zero;
    };

  template<class T> siscone_parton<T>::siscone_parton(int n) : 
  final_jets(), protojets(), protojets_axis(1<<n,T()), protojets_pt(1<<n,double(0.0)), protojets_pttilde(1<<n,double(0.0)),
    protojets_azimuthal_angle(1<<n,double(0.0)), protojets_rapidity(1<<n,double(0.0)), flag_rapidity(1<<n,bool(false)),
    flag_total_pt_is_zero(false)
  {
    final_jets.reserve(n);

    protojets.reserve(1<<n);
  }

  template<class T> void siscone_parton<T>::set_flag_total_pt_is_zero(bool flag)
    {
      flag_total_pt_is_zero = flag;
    }

  template<class T> int siscone_parton<T>::compute_jets(const std::vector<T> & p, double R, double f, int n_pass_max, double ptmin)
  {
    // internally we use CUT_Y = R*R
    double CUT_Y = R*R;

    // ------------------------
    bool flag_new_stable_cone_found;

    int n = p.size();

    // current_set = 2^n - 1
    // bitwise this just all 1's in the lowest n bits
    unsigned mask_all = (1<<n) - 1;
    unsigned current_set = mask_all;
    unsigned mask_remove = 0;

    int i_pass = 0;

    // clean final_jets vector and protojets vector from previous invocations 
    final_jets.clear();
    protojets.clear();

    // pre-calculate all possible cone axis, pt and pttilde values
    //
    // The calculation of the rapidity and the azimuthal angle involves 
    // the log-function and the atan2-function.
    // This is expensive.
    //
    // The following optimisation strategy is used:
    // The rapidity is calculated only when needed. Even with the brute force search strategy
    // for stable cones, which is used here, one can use the information on the azimuthal angles
    // to discard immediately some protocones.
    // For these protocones the rapidity does not need to be known.
    //
    // If the total pt-sum is zero, only half of all possible combinations have to be
    // computed for pt and azimuthal angle, the other half can be obtained from 
    // transverse momentum conservation.

    // pre-calculate cone axis, pt and pttilde
    for (unsigned i=0; i<n; i++)
      {
        for (unsigned j=0; j<(1<<i); j++)
          {
            unsigned index = (1<<i) + j;
            protojets_axis[index] = p[i];
            if (j) protojets_axis[index].sum_up(protojets_axis[j]);

            protojets_pt[index]              = protojets_axis[index].transverse_momentum();

            protojets_pttilde[index]         = protojets_pt[1<<i];
            if (j) protojets_pttilde[index] += protojets_pttilde[j];

            flag_rapidity[index] = false;
          }
      }

    // pre-calculate the azimuthal angles
    if ( flag_total_pt_is_zero )
      {
        for (unsigned index=0; index<(1<<(n-1)); index++) 
          {
            protojets_azimuthal_angle[index]          = protojets_axis[index].azimuthal_angle();
            protojets_azimuthal_angle[index^mask_all] = Pi + protojets_azimuthal_angle[index];
            if ( protojets_azimuthal_angle[index^mask_all] > Pi ) protojets_azimuthal_angle[index^mask_all] -= 2.0*Pi;
          }
      }
    else // !flag_total_pt_is_zero
      {
        for (unsigned index=0; index<(1<<n); index++) protojets_azimuthal_angle[index] = protojets_axis[index].azimuthal_angle();
      }

    // pre-calculate the rapidities for the n final state particles
    for (unsigned i=0; i<n; i++)
      {
        protojets_rapidity[1<<i] = protojets_axis[1<<i].rapidity();
        flag_rapidity[1<<i] = true;
      }

    // ------------------------------------------------
    // search stable cones
    do
      {
        flag_new_stable_cone_found = false;

        // update current_set and mask_remove
        current_set ^= mask_remove;
        mask_remove = 0;

        // there are pow(2,n) possibilities
        // we are not interested in (0,0,...,0) 
        // version 1.0.1: the case (1,1,...,1) is now included
        for (unsigned i_possibilities=1; i_possibilities<=mask_all; i_possibilities++)
          {
            // consider only configurations which are not in contradiction with current_set
            if ( !(i_possibilities & ~current_set) )
              {
                bool flag_stable_cone=true;

                // Pre-check for stable cone: 
                // A possibility can be discarded immediately, if a particle belonging to a protocone 
                // has already an azimuthal-angle distance larger than the cone size.
                // For those protocones the rapidity does not need to be known.
                for (size_t i=0; i<n; i++)
                  {
                    // has to be in current_set
                    if ( flag_stable_cone && (current_set & (1<<i))  )
                      {
                        double delta_phi_angle = fabs( protojets_azimuthal_angle[1<<i] - protojets_azimuthal_angle[i_possibilities] );
                        if (delta_phi_angle > Pi) delta_phi_angle = 2*Pi-delta_phi_angle;

                        double y_R_i = delta_phi_angle*delta_phi_angle;

                        if ( i_possibilities & (1<<i) )
                          {
                            // belongs to cone, but ouside
                            if ( y_R_i>CUT_Y ) flag_stable_cone=false;
                          }
                      }
                  } // i=0; i<n; i++

                // If the pre-check has been passed, we do a full check.
                // At this stage we need to know the rapidity of the protocone.
                if ( flag_stable_cone )
                  {
                    if ( !flag_rapidity[i_possibilities] )
                      {
                        protojets_rapidity[i_possibilities] = protojets_axis[i_possibilities].rapidity();
                        flag_rapidity[i_possibilities] = true;
                      }

                    for (size_t i=0; i<n; i++)
                      {
                        // has to be in current_set
                        if ( flag_stable_cone && (current_set & (1<<i))  )
                          {
                            double delta_rapidity  = protojets_rapidity[1<<i]  - protojets_rapidity[i_possibilities];
                            double delta_phi_angle = fabs( protojets_azimuthal_angle[1<<i] - protojets_azimuthal_angle[i_possibilities] );
                            if (delta_phi_angle > Pi) delta_phi_angle = 2*Pi-delta_phi_angle;

                            double y_R_i = delta_rapidity*delta_rapidity + delta_phi_angle*delta_phi_angle;

                            if ( i_possibilities & (1<<i) )
                              {
                                // belongs to cone, but ouside
                                if ( y_R_i>CUT_Y ) flag_stable_cone=false;
                              }
                            else // !(i_possibilities & (1<<i))
                              {
                                // does not belong to cone, but inside
                                if ( y_R_i<=CUT_Y ) flag_stable_cone=false;
                              }
                          }
                      } // i=0; i<n; i++
                  }

                // stable cone found
                if (flag_stable_cone)
                  {
                    flag_new_stable_cone_found=true;

                    // add to protojets
                    protojets.push_back(i_possibilities);

                    // mark for removal
                    mask_remove |= i_possibilities;
                  }
              } // !(i_possibilities & ~current_set)
          } // i_possibilities

        i_pass++;
      }
    while (flag_new_stable_cone_found && (( i_pass < n_pass_max ) || ( n_pass_max == 0 )) );


    // ------------------------------------------------
    // split-merge procedure

    while ( protojets.size() )
      {
        // remove protojets with pt<ptmin
        for (int i=protojets.size()-1; i>=0; i--)
          {
            if ( protojets_pt[protojets[i]] < ptmin ) protojets.erase(protojets.begin()+i);
          }

        // find the protojet with the largest pttilde
        if ( protojets.size() )
          {
            size_t i_max = 0;
            double pttilde_max = protojets_pttilde[protojets[0]];

            for (size_t i=1; i<protojets.size(); i++)
              {
                if ( protojets_pttilde[protojets[i]] > pttilde_max )
                  {
                    i_max = i;
                    pttilde_max = protojets_pttilde[protojets[i]];
                  }
              }

            // find overlapping protojet with largest pttilde
            bool flag_overlap = false;
            size_t i_overlap;
            double pttilde_overlap;
            for (size_t i=0; i<protojets.size(); i++)
              {
                if ( ( i != i_max ) && (protojets[i] & protojets[i_max]) )
                  {
                    if ( flag_overlap )
                      {
                        // check for largest pttilde_overlap
                        if ( protojets_pttilde[protojets[i]] > pttilde_overlap )
                          {
                            i_overlap = i;
                            pttilde_overlap = protojets_pttilde[protojets[i]];
                          }
                      }
                    else // !flag_overlap
                      {
                        // initialise the variables
                        i_overlap = i;
                        pttilde_overlap = protojets_pttilde[protojets[i]];

                        flag_overlap = true;
                      }
                  }
              }

            if ( flag_overlap )
              {
                unsigned index_i_max     = protojets[i_max];
                unsigned index_i_overlap = protojets[i_overlap];

                unsigned mask_shared = index_i_max & index_i_overlap;

                // determine shared pttilde
                double ptoverlap = protojets_pttilde[index_i_overlap];
                double ptshared  = protojets_pttilde[mask_shared];

                if ( ptshared < f*ptoverlap )
                  {
                    // split
                    unsigned assign_to_i_max = 0;
                    unsigned assign_to_i_overlap = 0;

                    // we need the rapidities
                    if ( !flag_rapidity[index_i_max] )
                      {
                        protojets_rapidity[index_i_max] = protojets_axis[index_i_max].rapidity();
                        flag_rapidity[index_i_max] = true;
                      }

                    if ( !flag_rapidity[index_i_overlap] )
                      {
                        protojets_rapidity[index_i_overlap] = protojets_axis[index_i_overlap].rapidity();
                        flag_rapidity[index_i_overlap] = true;
                      }

                    for (int i=0; i<n; i++)
                      {
                        unsigned mask_shifted = (1<<i);

                        if ( mask_shared & mask_shifted )
                          {
                            double delta_rapidity_i_imax  = protojets_rapidity[1<<i]  - protojets_rapidity[index_i_max];
                            double delta_phi_angle_i_imax = fabs( protojets_azimuthal_angle[1<<i] - protojets_azimuthal_angle[index_i_max] );
                            if (delta_phi_angle_i_imax > Pi) delta_phi_angle_i_imax = 2*Pi-delta_phi_angle_i_imax;


                            double delta_rapidity_i_ioverlap  = protojets_rapidity[1<<i]  - protojets_rapidity[index_i_overlap];
                            double delta_phi_angle_i_ioverlap = fabs( protojets_azimuthal_angle[1<<i] - protojets_azimuthal_angle[index_i_overlap] );
                            if (delta_phi_angle_i_ioverlap > Pi) delta_phi_angle_i_ioverlap = 2*Pi-delta_phi_angle_i_ioverlap;

                            double y_i_imax     = delta_rapidity_i_imax*delta_rapidity_i_imax + delta_phi_angle_i_imax*delta_phi_angle_i_imax;
                            double y_i_ioverlap = delta_rapidity_i_ioverlap*delta_rapidity_i_ioverlap + delta_phi_angle_i_ioverlap*delta_phi_angle_i_ioverlap;

                            if ( y_i_imax < y_i_ioverlap )
                              {
                                // i is assigned to protojet i_max
                                assign_to_i_max |= mask_shifted;
                              }
                            else // y_i_max >= y_i_overlap
                              {
                                // i is assigned to protojet i_overlap
                                assign_to_i_overlap |= mask_shifted;
                              }
                          } // mask_shared & mask_shifted
                      } // i=0; i<n

                    // particles assigned to i_max ? Then update i_overlap
                    if ( assign_to_i_max ) protojets[i_overlap] ^= assign_to_i_max;

                    // particles assigned to i_overlap ? Then update i_max
                    if ( assign_to_i_overlap ) protojets[i_max] ^= assign_to_i_overlap;

                  } // ptshared < f*ptoverlap
                else // ptshared >= f*ptoverlap
                  {
                    // merge
                    protojets[i_max] |= index_i_overlap;
                    protojets.erase(protojets.begin()+i_overlap);
                  }
              }
            else // !flag_overlap
              {
                final_jets.push_back( protojets[i_max] );

                protojets.erase(protojets.begin()+i_max);
              }

          } // if ( protojets.size() )
      } // while ( protojets.size() )

    return final_jets.size();
  }

  template<class T> class siscone_spherical_parton
    {
      // ctors
    public:
      siscone_spherical_parton(int n_particles_max);

      // member functions
    public:
      int compute_jets(const std::vector<T> & particles, double R, double f, int n_pass_max=0, double Emin=0.0);

      // member variables 
    public:
      std::vector<unsigned> final_jets;

      std::vector<unsigned> protojets;

      std::vector<T> protojets_axis;
      
      std::vector<double> protojets_E;

      std::vector<double> protojets_Etilde;

      std::vector<double> protojets_f;
    };

  template<class T> siscone_spherical_parton<T>::siscone_spherical_parton(int n) : 
  final_jets(), protojets(), protojets_axis(1<<n,T()), protojets_E(1<<n,double(0.0)), protojets_Etilde(1<<n,double(0.0)), protojets_f(1<<n,double(0.0))
  {
    final_jets.reserve(n);

    protojets.reserve(1<<n);
  }

  template<class T> int siscone_spherical_parton<T>::compute_jets(const std::vector<T> & p, double R, double f, int n_pass_max, double Emin)
  {
    // the parameter R is the cone half-opening angle in the range [0,Pi/2],
    // internally we use CUT_Y, which is in the range [0,1].
    double CUT_Y = 1.0-cos(R);

    // ------------------------
    bool flag_new_stable_cone_found;

    int n = p.size();

    // current_set = 2^n - 1
    // bitwise this just all 1's in the lowest n bits
    unsigned mask_all = (1<<n) - 1;
    unsigned current_set = mask_all;
    unsigned mask_remove = 0;

    int i_pass = 0;

    // clean final_jets vector and protojets vector from previous invocations 
    final_jets.clear();
    protojets.clear();

    // pre-calculate all possible cone axis and normalise to unit length
    for (unsigned i=0; i<n; i++)
      {
        for (unsigned j=0; j<(1<<i); j++)
          {
            unsigned index = (1<<i) + j;
            protojets_axis[index] = p[i];
            if (j) protojets_axis[index].sum_up(protojets_axis[j]);
          }
      }
    for (unsigned i=1; i<mask_all+1; i++) 
      {
        protojets_E[i] = protojets_axis[i].energy_component();
        protojets_f[i] = sqrt( protojets_axis[i].spatial_norm() );
      }
    for (unsigned i=1; i<mask_all; i++) protojets_axis[i].rescale( 1.0/protojets_f[i] );
    // protojets_f holds the fudge factor for Etilde
    for (unsigned i=1; i<mask_all+1; i++) 
      {
        double fudge = protojets_f[i]/protojets_E[i];
        protojets_f[i] = fudge*fudge;
      }

    // ------------------------------------------------
    // search stable cones
    do
      {
        flag_new_stable_cone_found = false;

        // update current_set and mask_remove
        current_set ^= mask_remove;
        mask_remove = 0;

        // there are pow(2,n) possibilities
        // we are not interested in (0,0,...,0) and (1,1,...,1)
        for (unsigned i_possibilities=1; i_possibilities<mask_all; i_possibilities++)
          {
            // consider only configurations which are not in contradiction with current_set
            if ( !(i_possibilities & ~current_set) )
              {
                // check for stable cone and calculate Etilde
                bool flag_stable_cone=true;
                double Etilde = 0.0;
                for (size_t i=0; i<n; i++)
                  {
                    // has to be in current_set
                    if ( flag_stable_cone && (current_set & (1<<i))  )
                      {
                        double y_R_i = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[i_possibilities]);

                        if ( i_possibilities & (1<<i) )
                          {
                            if ( y_R_i<=CUT_Y ) 
                              {
                                // calculate Etilde
                                Etilde += p[i].energy_component() * (1.0+protojets_f[i_possibilities]*y_R_i*(2.0-y_R_i));
                              }
                            else // y_R_i>CUT_Y
                              {
                                // belongs to cone, but ouside
                                flag_stable_cone=false;
                              }
                          }
                        else // !(i_possibilities & (1<<i))
                          {
                            // does not belong to cone, but inside
                            if ( y_R_i<=CUT_Y ) flag_stable_cone=false;
                          }
                      }
                  } // i=0; i<n; i++

                // stable cone found
                if (flag_stable_cone)
                  {
                    flag_new_stable_cone_found=true;

                    // add to protojets
                    protojets.push_back(i_possibilities);
                    protojets_Etilde[i_possibilities] = Etilde;

                    // mark for removal
                    mask_remove |= i_possibilities;
                  }
              } // !(i_possibilities & ~current_set)
          } // i_possibilities

        i_pass++;
      }
    while (flag_new_stable_cone_found && (( i_pass < n_pass_max ) || ( n_pass_max == 0 )) );


    // ------------------------------------------------
    // split-merge procedure

    while ( protojets.size() )
      {
        // remove protojets with E<Emin
        for (int i=protojets.size()-1; i>=0; i--)
          {
            if ( protojets_E[protojets[i]] < Emin ) protojets.erase(protojets.begin()+i);
          }

        // find the protojet with the largest Etilde
        if ( protojets.size() )
          {
            size_t i_max = 0;
            double Etilde_max = protojets_Etilde[protojets[0]];

            for (size_t i=1; i<protojets.size(); i++)
              {
                if ( protojets_Etilde[protojets[i]] > Etilde_max )
                  {
                    i_max = i;
                    Etilde_max = protojets_Etilde[protojets[i]];
                  }
              }

            // find overlapping protojet with largest Etilde
            bool flag_overlap = false;
            size_t i_overlap;
            double Etilde_overlap;
            for (size_t i=0; i<protojets.size(); i++)
              {
                if ( ( i != i_max ) && (protojets[i] & protojets[i_max]) )
                  {
                    if ( flag_overlap )
                      {
                        // check for largest Etilde_overlap
                        if ( protojets_Etilde[protojets[i]] > Etilde_overlap )
                          {
                            i_overlap = i;
                            Etilde_overlap = protojets_Etilde[protojets[i]];
                          }
                      }
                    else // !flag_overlap
                      {
                        // initialise the variables
                        i_overlap = i;
                        Etilde_overlap = protojets_Etilde[protojets[i]];

                        flag_overlap = true;
                      }
                  }
              }

            if ( flag_overlap )
              {
                unsigned mask_shared = protojets[i_max] & protojets[i_overlap];

                // determine shared energy
                double Eoverlap = protojets_E[protojets[i_overlap]];
                double Eshared  = protojets_E[mask_shared];

                if ( Eshared < f*Eoverlap )
                  {
                    // split
                    unsigned assign_to_i_max = 0;
                    unsigned assign_to_i_overlap = 0;

                    for (int i=0; i<n; i++)
                      {
                        unsigned mask_shifted = (1<<i);

                        if ( mask_shared & mask_shifted )
                          {
                            double y_i_imax     = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[protojets[i_max]]);
                            double y_i_ioverlap = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[protojets[i_overlap]]);

                            if ( y_i_imax < y_i_ioverlap )
                              {
                                // i is assigned to protojet i_max
                                assign_to_i_max |= mask_shifted;
                              }
                            else // y_i_max >= y_i_overlap
                              {
                                // i is assigned to protojet i_overlap
                                assign_to_i_overlap |= mask_shifted;
                              }
                          } // mask_shared & mask_shifted
                      } // i=0; i<n

                    // particles assigned to i_max ? Then update i_overlap
                    if ( assign_to_i_max )
                      {
                        protojets[i_overlap] ^= assign_to_i_max;

                        // calculate new Etilde
                        double Etilde = 0.0;
                        for (int i=0; i<n; i++)
                          {
                            if ( protojets[i_overlap] & (1<<i) )
                              {
                                double y_R_i = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[protojets[i_overlap]]);
                                Etilde += p[i].energy_component() *(1.0+protojets_f[protojets[i_overlap]]*y_R_i*(2.0-y_R_i));
                              }
                          }
                        protojets_Etilde[protojets[i_overlap]] = Etilde;
                      }

                    // particles assigned to i_overlap ? Then update i_max
                    if ( assign_to_i_overlap )
                      {
                        protojets[i_max] ^= assign_to_i_overlap;

                        // calculate new Etilde
                        double Etilde = 0.0;
                        for (int i=0; i<n; i++)
                          {
                            if ( protojets[i_max] & (1<<i) )
                              {
                                double y_R_i = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[protojets[i_max]]);
                                Etilde += p[i].energy_component() *(1.0+protojets_f[protojets[i_max]]*y_R_i*(2.0-y_R_i));
                              }
                          }
                        protojets_Etilde[protojets[i_max]] = Etilde;
                      }

                  } // Eshared < f*Eoverlap
                else // Eshared >= f*Eoverlap
                  {
                    // merge
                    unsigned mask_union = protojets[i_max] | protojets[i_overlap];

                    // calculate new Etilde
                    double Etilde = 0.0;
                    for (int i=0; i<n; i++)
                      {
                        if ( mask_union & (1<<i) )
                          {
                            double y_R_i = 1.0 - protojets_axis[1<<i].spatial_scalar_product(protojets_axis[mask_union]);
                            Etilde += p[i].energy_component() *(1.0+protojets_f[mask_union]*y_R_i*(2.0-y_R_i));
                          }
                      }
                    protojets_Etilde[mask_union] = Etilde;

                    // update protojets
                    protojets[i_max] = mask_union;
                    protojets.erase(protojets.begin()+i_overlap);
                  }
              }
            else // !flag_overlap
              {
                final_jets.push_back( protojets[i_max] );

                protojets.erase(protojets.begin()+i_max);
              }

          } // if ( protojets.size() )
      } // while ( protojets.size() )

    return final_jets.size();
  }

} // namespace siscone_jet_algorithm

#endif // ndef __SISCONE_PARTON_H__

---