LazyTiling9

Collaboration diagram for LazyTiling9:

Public Member Functions
	LazyTiling9 (ClusterSequence &cs)
void	run ()

Protected Member Functions
void	_initialise_tiles ()
int	_tile_index (int ieta, int iphi) const
void	_bj_remove_from_tiles (TiledJet *const jet)
int	_tile_index (const double eta, const double phi) const
void	_tj_set_jetinfo (TiledJet *const jet, const int _jets_index)
void	_print_tiles (TiledJet *briefjets) const
void	_add_neighbours_to_tile_union (const int tile_index, std::vector< int > &tile_union, int &n_near_tiles) const
void	_add_untagged_neighbours_to_tile_union (const int tile_index, std::vector< int > &tile_union, int &n_near_tiles)
void	_add_untagged_neighbours_to_tile_union_using_max_info (const TiledJet *const jet, std::vector< int > &tile_union, int &n_near_tiles)
double	_distance_to_tile (const TiledJet *bj, const Tile2 *) const
void	_update_jetX_jetI_NN (TiledJet *jetX, TiledJet *jetI, std::vector< TiledJet * > &jets_for_minheap)
void	_set_NN (TiledJet *jetI, std::vector< TiledJet * > &jets_for_minheap)
template<class J >
double	_bj_diJ (const J *const jet) const
template<class J >
void	_bj_set_jetinfo (J *const jetA, const int _jets_index) const
template<class J >
double	_bj_dist (const J *const jetA, const J *const jetB) const
template<class J >
double	_bj_dist_not_periodic (const J *const jetA, const J *const jetB) const

Protected Attributes
ClusterSequence &	_cs
const std::vector< PseudoJet > &	_jets
std::vector< Tile2 >	_tiles
double	_Rparam
double	_R2
double	_invR2
double	_tiles_eta_min
double	_tiles_eta_max
double	_tile_size_eta
double	_tile_size_phi
double	_tile_half_size_eta
double	_tile_half_size_phi
int	_n_tiles_phi
int	_tiles_ieta_min
int	_tiles_ieta_max
std::vector< TiledJet * >	_jets_for_minheap

Constructor & Destructor Documentation

LazyTiling9()

FJCORE_BEGIN_NAMESPACE LazyTiling9

(

ClusterSequence &

cs

)

                                             :
  _cs(cs), _jets(cs.jets())
{
#ifdef INSTRUMENT2
  _ncall = 0; // gps tmp
  _ncall_dtt = 0; // gps tmp
#endif // INSTRUMENT2
  _Rparam = cs.jet_def().R();
  _R2 = _Rparam * _Rparam;
  _invR2 = 1.0 / _R2;
  _initialise_tiles();
}

Member Function Documentation

_add_neighbours_to_tile_union()

void _add_neighbours_to_tile_union ( const int  tile_index, std::vector< int > &  tile_union, int &  n_near_tiles  ) const

protected

                                                               {
  for (Tile2 * const * near_tile = _tiles[tile_index].begin_tiles;
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
    n_near_tiles++;
  }
}

_add_untagged_neighbours_to_tile_union()

void _add_untagged_neighbours_to_tile_union ( const int  tile_index, std::vector< int > &  tile_union, int &  n_near_tiles  )

inlineprotected

                                                          {
  for (Tile2 ** near_tile = _tiles[tile_index].begin_tiles;
       near_tile != _tiles[tile_index].end_tiles; near_tile++){
    if (! (*near_tile)->tagged) {
      (*near_tile)->tagged = true;
      tile_union[n_near_tiles] = *near_tile - & _tiles[0];
      n_near_tiles++;
    }
  }
}

_add_untagged_neighbours_to_tile_union_using_max_info()

void _add_untagged_neighbours_to_tile_union_using_max_info ( const TiledJet *const  jet, std::vector< int > &  tile_union, int &  n_near_tiles  )

inlineprotected

                                                          {
  Tile2 & tile = _tiles[jet->tile_index];
  for (Tile2 ** near_tile = tile.begin_tiles; near_tile != tile.end_tiles; near_tile++){
    if ((*near_tile)->tagged) continue;
    double dist = _distance_to_tile(jet, *near_tile) - tile_edge_security_margin;
    if (dist > (*near_tile)->max_NN_dist) continue;
    (*near_tile)->tagged = true;
    tile_union[n_near_tiles] = *near_tile - & _tiles[0];
    n_near_tiles++;
  }
}

_bj_diJ()

double _bj_diJ ( const J *const  jet ) const

inlineprotected

                                                               {
    double kt2 = jet->kt2;
    if (jet->NN != NULL) {if (jet->NN->kt2 < kt2) {kt2 = jet->NN->kt2;}}
    return jet->NN_dist * kt2;
  }

_bj_dist()

double _bj_dist ( const J *const  jetA, const J *const  jetB  ) const

inlineprotected

    {
    _ncall++; // GPS tmp
#else
    const {
#endif
    double dphi = std::abs(jetA->phi - jetB->phi);
    double deta = (jetA->eta - jetB->eta);
    if (dphi > pi) {dphi = twopi - dphi;}
    return dphi*dphi + deta*deta;
  }

_bj_dist_not_periodic()

double _bj_dist_not_periodic ( const J *const  jetA, const J *const  jetB  ) const

inlineprotected

    {
    _ncall++; // GPS tmp
#else
    const {
#endif
    double dphi = jetA->phi - jetB->phi;
    double deta = (jetA->eta - jetB->eta);
    return dphi*dphi + deta*deta;
  }

_bj_remove_from_tiles()

void _bj_remove_from_tiles ( TiledJet *const  jet )

protected

                                                            {
  Tile2 * tile = & _tiles[jet->tile_index];
  if (jet->previous == NULL) {
    tile->head = jet->next;
  } else {
    jet->previous->next = jet->next;
  }
  if (jet->next != NULL) {
    jet->next->previous = jet->previous;
  }
}

_bj_set_jetinfo()

void _bj_set_jetinfo ( J *const  jetA, const int  _jets_index  ) const

inlineprotected

                                                                         {
    jetA->eta  = _jets[_jets_index].rap();
    jetA->phi  = _jets[_jets_index].phi_02pi();
    jetA->kt2  = _cs.jet_scale_for_algorithm(_jets[_jets_index]);
    jetA->_jets_index = _jets_index;
    jetA->NN_dist = _R2;
    jetA->NN      = NULL;
  }

_distance_to_tile()

double _distance_to_tile ( const TiledJet *  bj, const Tile2 *  tile  ) const

inlineprotected

   {
  _ncall_dtt++; // GPS tmp
#else
  const {
#endif // INSTRUMENT2
  double deta;
  if (_tiles[bj->tile_index].eta_centre == tile->eta_centre) deta = 0;
  else   deta = std::abs(bj->eta - tile->eta_centre) - _tile_half_size_eta;
  double dphi = std::abs(bj->phi - tile->phi_centre);
  if (dphi > pi) dphi = twopi-dphi;
  dphi -= _tile_half_size_phi;
  if (dphi < 0) dphi = 0;
  return dphi*dphi + deta*deta;
}

_initialise_tiles()

void _initialise_tiles ( )

protected

                                    {
  double default_size = max(0.1,_Rparam);
  _tile_size_eta = default_size;
  _n_tiles_phi   = max(3,int(floor(twopi/default_size)));
  _tile_size_phi = twopi / _n_tiles_phi; // >= _Rparam and fits in 2pi
#ifdef _FJCORE_TILING2_USE_TILING_ANALYSIS_
  TilingExtent tiling_analysis(_cs);
  _tiles_eta_min = tiling_analysis.minrap();
  _tiles_eta_max = tiling_analysis.maxrap();
#else
  _tiles_eta_min = 0.0;
  _tiles_eta_max = 0.0;
  const double maxrap = 7.0;
  for(unsigned int i = 0; i < _jets.size(); i++) {
    double eta = _jets[i].rap();
    if (abs(eta) < maxrap) {
      if (eta < _tiles_eta_min) {_tiles_eta_min = eta;}
      if (eta > _tiles_eta_max) {_tiles_eta_max = eta;}
    }
  }
#endif
# define FJCORE_LAZY9_MIN2TILESY
#ifdef FJCORE_LAZY9_MIN2TILESY
   if (_tiles_eta_max - _tiles_eta_min < 2*_tile_size_eta) {
     _tile_size_eta = (_tiles_eta_max - _tiles_eta_min)/2;
     _tiles_ieta_min = 0;
     _tiles_ieta_max = 1;
     _tiles_eta_max -= _tile_size_eta;
   } else {
#endif //FASTJET_LAZY9_MIN2TILESY
  _tiles_ieta_min = int(floor(_tiles_eta_min/_tile_size_eta));
  _tiles_ieta_max = int(floor( _tiles_eta_max/_tile_size_eta));
  _tiles_eta_min = _tiles_ieta_min * _tile_size_eta;
  _tiles_eta_max = _tiles_ieta_max * _tile_size_eta;
#ifdef FJCORE_LAZY9_MIN2TILESY
   }
#endif
  _tile_half_size_eta = _tile_size_eta * 0.5;
  _tile_half_size_phi = _tile_size_phi * 0.5;
  vector<bool> use_periodic_delta_phi(_n_tiles_phi, false);
  if (_n_tiles_phi <= 3) {
    fill(use_periodic_delta_phi.begin(), use_periodic_delta_phi.end(), true);
  } else {
    use_periodic_delta_phi[0] = true;
    use_periodic_delta_phi[_n_tiles_phi-1] = true;
  }
  _tiles.resize((_tiles_ieta_max-_tiles_ieta_min+1)*_n_tiles_phi);
  for (int ieta = _tiles_ieta_min; ieta <= _tiles_ieta_max; ieta++) {
    for (int iphi = 0; iphi < _n_tiles_phi; iphi++) {
      Tile2 * tile = & _tiles[_tile_index(ieta,iphi)];
      tile->head = NULL; // first element of tiles points to itself
      tile->begin_tiles[0] =  tile;
      Tile2 ** pptile = & (tile->begin_tiles[0]);
      pptile++;
      tile->surrounding_tiles = pptile;
      if (ieta > _tiles_ieta_min) {
    // with the itile subroutine, we can safely run tiles from
    // idphi=-1 to idphi=+1, because it takes care of
    // negative and positive boundaries
    for (int idphi = -1; idphi <=+1; idphi++) {
      *pptile = & _tiles[_tile_index(ieta-1,iphi+idphi)];
      pptile++;
    }
      }
      *pptile = & _tiles[_tile_index(ieta,iphi-1)];
      pptile++;
      tile->RH_tiles = pptile;
      *pptile = & _tiles[_tile_index(ieta,iphi+1)];
      pptile++;
      if (ieta < _tiles_ieta_max) {
    for (int idphi = -1; idphi <= +1; idphi++) {
      *pptile = & _tiles[_tile_index(ieta+1,iphi+idphi)];
      pptile++;
    }
      }
      tile->end_tiles = pptile;
      tile->tagged = false;
      tile->use_periodic_delta_phi = use_periodic_delta_phi[iphi];
      tile->max_NN_dist = 0;
      tile->eta_centre = (ieta-_tiles_ieta_min+0.5)*_tile_size_eta + _tiles_eta_min;
      tile->phi_centre = (iphi+0.5)*_tile_size_phi;
    }
  }
}

_print_tiles()

void _print_tiles ( TiledJet *  briefjets ) const

protected

                                                          {
  for (vector<Tile2>::const_iterator tile = _tiles.begin();
       tile < _tiles.end(); tile++) {
    cout << "Tile " << tile - _tiles.begin()<<" = ";
    vector<int> list;
    for (TiledJet * jetI = tile->head; jetI != NULL; jetI = jetI->next) {
      list.push_back(jetI-briefjets);
    }
    sort(list.begin(),list.end());
    for (unsigned int i = 0; i < list.size(); i++) {cout <<" "<<list[i];}
    cout <<"\n";
  }
}

_set_NN()

void _set_NN ( TiledJet *  jetI, std::vector< TiledJet * > &  jets_for_minheap  )

inlineprotected

                                                                     {
  jetI->NN_dist = _R2;
  jetI->NN      = NULL;
  if (!jetI->minheap_update_needed()) {
    jetI->label_minheap_update_needed();
    jets_for_minheap.push_back(jetI);}
  Tile2 * tile_ptr = &_tiles[jetI->tile_index];
    for (Tile2 ** near_tile  = tile_ptr->begin_tiles;
         near_tile != tile_ptr->end_tiles; near_tile++) {
      if (jetI->NN_dist < _distance_to_tile(jetI, *near_tile)) continue;
      for (TiledJet * jetJ  = (*near_tile)->head;
           jetJ != NULL; jetJ = jetJ->next) {
        double dist = _bj_dist(jetI,jetJ);
        if (dist < jetI->NN_dist && jetJ != jetI) {
          jetI->NN_dist = dist; jetI->NN = jetJ;
        }
      }
    }
}

_tile_index() [1/2]

int _tile_index ( const double  eta, const double  phi  ) const

protected

                                                                     {
  int ieta, iphi;
  if      (eta <= _tiles_eta_min) {ieta = 0;}
  else if (eta >= _tiles_eta_max) {ieta = _tiles_ieta_max-_tiles_ieta_min;}
  else {
    ieta = int(((eta - _tiles_eta_min) / _tile_size_eta));
    if (ieta > _tiles_ieta_max-_tiles_ieta_min) {
      ieta = _tiles_ieta_max-_tiles_ieta_min;}
  }
  iphi = int((phi+twopi)/_tile_size_phi) % _n_tiles_phi;
  return (iphi + ieta * _n_tiles_phi);
}

_tile_index() [2/2]

int _tile_index ( int  ieta, int  iphi  ) const

inlineprotected

                                                    {
    return (ieta-_tiles_ieta_min)*_n_tiles_phi
                  + (iphi+_n_tiles_phi) % _n_tiles_phi;
  }

_tj_set_jetinfo()

void _tj_set_jetinfo ( TiledJet *const  jet, const int  _jets_index  )

inlineprotected

                                                 {
  _bj_set_jetinfo<>(jet, _jets_index);
  jet->tile_index = _tile_index(jet->eta, jet->phi);
  Tile2 * tile = &_tiles[jet->tile_index];
  jet->previous   = NULL;
  jet->next       = tile->head;
  if (jet->next != NULL) {jet->next->previous = jet;}
  tile->head      = jet;
}

_update_jetX_jetI_NN()

void _update_jetX_jetI_NN ( TiledJet *  jetX, TiledJet *  jetI, std::vector< TiledJet * > &  jets_for_minheap  )

inlineprotected

                                                                                                                     {
  double dist = _bj_dist(jetI,jetX);
  if (dist < jetI->NN_dist) {
    if (jetI != jetX) {
      jetI->NN_dist = dist;
      jetI->NN = jetX;
      if (!jetI->minheap_update_needed()) {
    jetI->label_minheap_update_needed();
    jets_for_minheap.push_back(jetI);
      }
    }
  }
  if (dist < jetX->NN_dist) {
    if (jetI != jetX) {
      jetX->NN_dist = dist;
      jetX->NN      = jetI;}
  }
}

run()

void run

(

)

                      {
  int n = _jets.size();
  if (n == 0) return;
  TiledJet * briefjets = new TiledJet[n];
  TiledJet * jetA = briefjets, * jetB;
  TiledJet oldB = briefjets[0];
  vector<int> tile_union(3*n_tile_neighbours);
  for (int i = 0; i< n; i++) {
    _tj_set_jetinfo(jetA, i);
    jetA++; // move on to next entry of briefjets
  }
  TiledJet * head = briefjets; // a nicer way of naming start
  vector<Tile2>::iterator tile;
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      for (jetB = tile->head; jetB != jetA; jetB = jetB->next) {
    double dist = _bj_dist_not_periodic(jetA,jetB);
    if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
    if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
      }
    }
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
    }
  }
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    if (tile->use_periodic_delta_phi) {
      for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
          if (relevant_for_jetA || relevant_for_RTile) {
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
              double dist = _bj_dist(jetA,jetB);
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
            }
          }
        }
      }
    } else {
      for (Tile2 ** RTile = tile->RH_tiles; RTile != tile->end_tiles; RTile++) {
        for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
          double dist_to_tile = _distance_to_tile(jetA, *RTile);
          bool relevant_for_jetA  = dist_to_tile <= jetA->NN_dist;
          bool relevant_for_RTile = dist_to_tile <= (*RTile)->max_NN_dist;
          if (relevant_for_jetA || relevant_for_RTile) {
            for (jetB = (*RTile)->head; jetB != NULL; jetB = jetB->next) {
              double dist = _bj_dist_not_periodic(jetA,jetB);
              if (dist < jetA->NN_dist) {jetA->NN_dist = dist; jetA->NN = jetB;}
              if (dist < jetB->NN_dist) {jetB->NN_dist = dist; jetB->NN = jetA;}
            }
          }
        }
      }
    }
  }
  for (tile = _tiles.begin(); tile != _tiles.end(); tile++) {
    tile->max_NN_dist = 0;
    for (jetA = tile->head; jetA != NULL; jetA = jetA->next) {
      if (jetA->NN_dist > tile->max_NN_dist) tile->max_NN_dist = jetA->NN_dist;
    }
  }
#ifdef INSTRUMENT2
  cout << "intermediate ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
#endif // INSTRUMENT2
  vector<double> diJs(n);
  for (int i = 0; i < n; i++) {
    diJs[i] = _bj_diJ(&briefjets[i]);
    briefjets[i].label_minheap_update_done();
  }
  MinHeap minheap(diJs);
  vector<TiledJet *> jets_for_minheap;
  jets_for_minheap.reserve(n);
  int history_location = n-1;
  while (n > 0) {
    double diJ_min = minheap.minval() *_invR2;
    jetA = head + minheap.minloc();
    history_location++;
    jetB = jetA->NN;
    if (jetB != NULL) {
      if (jetA < jetB) {std::swap(jetA,jetB);}
      int nn; // new jet index
      _cs.plugin_record_ij_recombination(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
      _bj_remove_from_tiles(jetA);
      oldB = * jetB;  // take a copy because we will need it...
      _bj_remove_from_tiles(jetB);
      _tj_set_jetinfo(jetB, nn); // cause jetB to become _jets[nn]
    } else {
      _cs.plugin_record_iB_recombination(jetA->_jets_index, diJ_min);
      _bj_remove_from_tiles(jetA);
    }
    minheap.remove(jetA-head);
    int n_near_tiles = 0;
    if (jetB != NULL) {
      Tile2 & jetB_tile = _tiles[jetB->tile_index];
      for (Tile2 ** near_tile  = jetB_tile.begin_tiles;
               near_tile != jetB_tile.end_tiles; near_tile++) {
        double dist_to_tile = _distance_to_tile(jetB, *near_tile);
        bool relevant_for_jetB  = dist_to_tile <= jetB->NN_dist;
        bool relevant_for_near_tile = dist_to_tile <= (*near_tile)->max_NN_dist;
        bool relevant = relevant_for_jetB || relevant_for_near_tile;
        if (! relevant) continue;
        tile_union[n_near_tiles] = *near_tile - & _tiles[0];
        (*near_tile)->tagged = true;
        n_near_tiles++;
        for (TiledJet * jetI = (*near_tile)->head; jetI != NULL; jetI = jetI->next) {
          if (jetI->NN == jetA || jetI->NN == jetB) _set_NN(jetI, jets_for_minheap);
          _update_jetX_jetI_NN(jetB, jetI, jets_for_minheap);
        }
      }
    }
    int n_done_tiles = n_near_tiles;
    _add_untagged_neighbours_to_tile_union_using_max_info(jetA,
                           tile_union, n_near_tiles);
    if (jetB != NULL) {
    _add_untagged_neighbours_to_tile_union_using_max_info(&oldB,
                                  tile_union,n_near_tiles);
      jetB->label_minheap_update_needed();
      jets_for_minheap.push_back(jetB);
    }
    for (int itile = 0; itile < n_done_tiles; itile++) {
      _tiles[tile_union[itile]].tagged = false;
    }
    for (int itile = n_done_tiles; itile < n_near_tiles; itile++) {
      Tile2 * tile_ptr = &_tiles[tile_union[itile]];
      tile_ptr->tagged = false;
      for (TiledJet * jetI = tile_ptr->head; jetI != NULL; jetI = jetI->next) {
        if (jetI->NN == jetA || (jetI->NN == jetB && jetB != NULL)) {
          _set_NN(jetI, jets_for_minheap);
        }
      }
    }
    while (jets_for_minheap.size() > 0) {
      TiledJet * jetI = jets_for_minheap.back();
      jets_for_minheap.pop_back();
      minheap.update(jetI-head, _bj_diJ(jetI));
      jetI->label_minheap_update_done();
      Tile2 & tile_I = _tiles[jetI->tile_index];
      if (tile_I.max_NN_dist < jetI->NN_dist) tile_I.max_NN_dist = jetI->NN_dist;
    }
    n--;
  }
  delete[] briefjets;
#ifdef INSTRUMENT2
  cout << "ncall, dtt = " << _ncall << " " << _ncall_dtt << endl; // GPS tmp
#endif // INSTRUMENT2
}

Member Data Documentation

_cs

ClusterSequence& _cs

protected

_invR2

double _invR2

protected

_jets

const std::vector<PseudoJet>& _jets

protected

_jets_for_minheap

std::vector<TiledJet *> _jets_for_minheap

protected

_n_tiles_phi

int _n_tiles_phi

protected

_R2

double _R2

protected

_Rparam

double _Rparam

protected

_tile_half_size_eta

double _tile_half_size_eta

protected

_tile_half_size_phi

double _tile_half_size_phi

protected

_tile_size_eta

double _tile_size_eta

protected

_tile_size_phi

double _tile_size_phi

protected

_tiles

std::vector<Tile2> _tiles

protected

_tiles_eta_max

double _tiles_eta_max

protected

_tiles_eta_min

double _tiles_eta_min

protected

_tiles_ieta_max

int _tiles_ieta_max

protected

_tiles_ieta_min

int _tiles_ieta_min

protected

---

The documentation for this class was generated from the following file:

• /builds/cms-analysis/general/DasAnalysisSystem/Core/Installer/Core/JetObservables/src/fjcore.cc