particle_halo2.h

Go to the documentation of this file.

 
#ifndef GLAMER_particle_halo2_h
#define GLAMER_particle_halo2_h
#include <thread>
#include <mutex>
 
#include "geometry.h"
#include "quadTree.h"
#include "simpleTree.h"
#include "particle_types.h"
#include "utilities_slsim.h"
#include "lens_halos.h"
#include "oTreeNB.h"
#include "gadget.hh"
 
template<typename DType = float >
class LensHaloParticles : public LensHalo
{
  public:
    LensHaloParticles(
      const std::string& simulation_filename 
      ,SimFileFormat format   
      ,PosType redshift        
      ,double my_inv_area_phys           
      ,PosType mass_particle  
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5             
      ,float Nbucket = 8  
      ,float theta = 0.1      
      ,bool recenter = false 
      ,bool verbose = false
    ) : LensHalo(redshift,cosmo),particle_mass(mass_particle)
       ,inv_area(my_inv_area_phys),Nsmooth(my_Nsmooth),theta2(theta*theta)
       ,Nbucket(Nbucket)
    {
 
      switch (format) {
        case  SimFileFormat::ascii:
          readPositionFileASCII(simulation_filename);
          break;
        default:
          std::cerr << "LensHaloParticles does not accept format of particle data file." << std::endl;
          throw std::invalid_argument("bad format");
      }
      if(verbose) std::cout << "   Read " << pp.size() << " particles from " << simulation_filename << std::endl;
 
      // convert from comoving to physical coordinates
      double scale_factor = 1.0/(1 + redshift);
      for(auto &p : pp){
            p[0] *= scale_factor;
            p[1] *= scale_factor;
            p[2] *= scale_factor;
      }
      setUp(recenter,verbose);
    }
 
    LensHaloParticles(
      std::vector<Point_3d<DType> > &pvector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,PosType mass_particle   
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHalo(redshift,cosmo),particle_mass(mass_particle)
       ,inv_area(my_inv_area_phys),Nsmooth(my_Nsmooth),theta2(theta*theta)
       ,Nbucket(Nbucket)
    {
      std::swap(pp, pvector); // take ownership of the data
      setUp(recenter,verbose);
    }
 
    ~LensHaloParticles(){};
    
    void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
    ){
 
      otree->force2D(xcm,particle_mass,Nsmooth,theta2,inv_area
        ,alpha,kappa,gamma,phi);
    };
 
    size_t getN() const { return pp.size(); };
 
    Point_3d<DType> CenterOfMass(){return mcenter;}
 
    LensHaloParticles(LensHaloParticles &&h):LensHalo(std::move(h)){
      mcenter = h.mcenter;
      pp = std::move(h.pp);
      inv_area = h.inv_area;
      Nsmooth = h.Nsmooth;
      theta2 = h.theta2;
      particle_mass = h.particle_mass;
      Nbucket = h.Nbucket;
      otree = std::move(h.otree);
   }
    LensHaloParticles & operator=(LensHaloParticles &&h){
    if(this == &h) return *this;
      LensHalo::operator=(std::move(h));
      mcenter = h.mcenter;
      pp = std::move(h.pp);
      inv_area = h.inv_area;
      Nsmooth = h.Nsmooth;
      theta2 = h.theta2;
      particle_mass = h.particle_mass;
      Nbucket = h.Nbucket;
      otree = std::move(h.otree);
      return *this;
    }
 
    void rotate(
      PosType theta_x
      ,PosType theta_y
      ,Point_3d<DType> &xo);
    
    // reset the opening angle used in the tree force calculation
    void resetForceAngle(double theta){
      theta2 = theta*theta;
    }
 
    void translate_particles(Point_3d<DType> &xo);
 
    static std::vector< LensHaloParticles<DType> > make(
      std::string filename     
      ,SimFileFormat format    
      ,float umass             
      ,float ulength           
      ,double inv_area         
      ,const COSMOLOGY& cosmo  
      ,int Nsmooth = 5         
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,DType rmax = 0          
      ,Point_3d<DType> xo = Point_3d<DType>(0,0,0) 
      ,bool verbose = false
    ){
      double r2max= rmax*rmax;
      
        switch (format) {
          case SimFileFormat::gadget2 : 
            {
 
              std::vector<Point_3d<DType> > tmpdata;
              GadgetFile<Point_3d<DType> > gadget_file(filename,tmpdata);
              gadget_file.openFile();
              gadget_file.readPOS();
              gadget_file.closeFile();
 
              if(rmax <= 0){
                r2max = 0; // no excision
                xo[0] = gadget_file.boxl/2; // center of box
                xo[1] = gadget_file.boxl/2;
                xo[2] = gadget_file.boxl/2;
              }
              // convert to physical Mpc
              const double a = ulength*1.0e-6/(1 + gadget_file.redshift)/gadget_file.hubble;
              // split particles into different types
              std::vector< std::vector<Point_3d<DType> > > pos(6);
              size_t n=0;
              int ntypes = 0;
              Point_3d<DType> dx;
              for(int i=0; i<6 ; ++i){
                pos[i].resize(gadget_file.npart[i]);
                if(gadget_file.npart[i] > 0) ++ntypes;
                if(r2max > 0){
                  size_t k = 0;
                  for(size_t j=0; j<gadget_file.npart[i]; ++j){
                      pos[i][k] = tmpdata[n++];
                      dx = pos[i][k] - xo;
                      if(dx.length_sqr() < r2max){ 
                        pos[i][k] -= xo;
                        pos[i][k] *= a; // convert to physical Mpc
                        ++k;
                      }
                  }
                  pos[i].resize(k); // resize to number of particles in sphere
                  gadget_file.npart[i] = k; // update number of particles
                }else{
                  for(size_t j=0; j<gadget_file.npart[i]; ++j){
                      pos[i][j] = (tmpdata[n++]-xo)*a;
                  }
                }
              }
              tmpdata.clear();
 
              std::vector< LensHaloParticles<DType> > halos;
              ntypes = 0;
              for(int i=0; i<6 ; ++i){
                if(gadget_file.npart[i] > 0){
                  halos.push_back( LensHaloParticles<DType>(
                    pos[i]
                    ,gadget_file.redshift
                    ,inv_area
                    ,gadget_file.masstab[i]*umass/gadget_file.hubble
                    ,cosmo
                    ,Nsmooth
                    ,Nbucket
                    ,theta
                    ,false
                    ,verbose)
                  );
                  ++ntypes;
                }
              }
              return halos;
            }
            break;
          default:
            std::cerr << "LensHaloParticles does not accept format of particle data file." << std::endl;
            throw std::invalid_argument("bad format");
        }
      }
 
      void getBoundingBox(
        Point_3d<DType> &p1,  
        Point_3d<DType> &p2   
      ) const {
        otree->getBoundingBox(p1,p2);
      }
 
    void swapout(std::vector<Point_3d<DType> > &pvector){
      std::swap(pp, pvector);
      otree.reset();
    }
protected :
    Point_3d<DType> mcenter;
    //DType *pp;
    std::vector<Point_3d<DType> > pp;
    int Nsmooth;  
    int Nbucket;  
    double theta2;  
    double particle_mass;  
  
    PosType inv_area;
    std::unique_ptr<OTreeNB<Point_3d<DType>, DType > > otree;
 
    void readPositionFileASCII(const std::string &filename);
    // construct tree, particles positions must already by stored in comoving Mpc
    void setUp(
        bool recenter           
        ,bool verbose
    );
 
    // protected constructor for creating an empty one
    LensHaloParticles(
      PosType redshift        
      ,double my_inv_area_phys      
      ,PosType mass_particle   
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHalo(redshift,cosmo),particle_mass(mass_particle)
       ,inv_area(my_inv_area_phys),Nsmooth(my_Nsmooth),theta2(theta*theta)
       ,Nbucket(Nbucket)
    {
    }
};
 
// construct tree, particles positions must already by stored in comoving Mpc
template<typename DType>
void LensHaloParticles<DType>::setUp(
        bool recenter           
        ,bool verbose
    ){
        long Npoints = pp.size();
  
        Point_3d<double> mcenter(0,0,0);
 
        for(auto &p : pp){
            mcenter[0] += p[0];
            mcenter[1] += p[1];
            mcenter[2] += p[2];
        }
  
        mcenter /= pp.size();
  
        if(verbose) std::cout << "   Center of mass : " << mcenter << std::endl;
  
        if(recenter){
            PosType r2,r2max=0;
            for(auto &p : pp){
                p[0] -= mcenter[0];
                p[1] -= mcenter[1];
                p[2] -= mcenter[2];
            }
            mcenter *= 0;
        }
  
        LensHalo::setTheta(0,0);
        LensHalo::mass = particle_mass * pp.size();
        LensHalo::set_flag_elliptical(false); // shouldn't do anything
 
        if(verbose) std::cout << "   Creating particle tree ... " << std::endl;
        otree.reset( new OTreeNB<Point_3d<DType>,DType >(pp.data(),pp.size()) );
        otree->build(Nbucket);
        if(verbose) std::cout << "done." << std::endl;
};
 
template<typename DType>
void LensHaloParticles<DType>::readPositionFileASCII(const std::string &filename
                                                     ){
  
  int ncoll = Utilities::IO::CountColumns(filename);
  if(ncoll < 3 ){
    std::cerr << filename << " should have three columns!" << std::endl;
  }
  if(ncoll > 4 ){
    std::cerr << filename << " Warning! : LensHaloParticles is not using masses." << std::endl;
  }
  
  std::ifstream myfile(filename);
  long Npoints = 0;
  // find number of particles
  if (myfile.is_open()){
    
    float tmp_mass = 0.0;
    std::string str,label;
    int count =0;
    while(std::getline(myfile, str)){
      std::stringstream ss(str);
      ss >> label;
      if(label == "#"){
        ss >> label;
        if(label == "nparticles"){
          ss >> Npoints;
          ++count;
        }
        {
          if(label == "mass"){
            ss >> tmp_mass;
            ++count;
          }
        }
      }else break;
      if(count == 2 ) break;
    }
    
    if(count == 0){
      std::cerr << "File " << filename << " must have the header lines: " << std::endl
        << "# nparticles   ****" << std::endl;
      throw std::runtime_error("file reading error");
    }
    
    assert(Npoints > 0);
    pp.resize(Npoints);
    size_t row = 0;
    
    // read in particle positions
    do{
        if(str[0] == '#') continue; //for comments
        std::stringstream ss(str);
        
        ss >> pp[row][0];
        if(!(ss >> pp[row][1])) std::cerr << "3 columns are expected in line " << row
          << " of " << filename << std::endl;
        if(!(ss >> pp[row][2])) std::cerr << "3 columns are expected in line " << row
          << " of " << filename << std::endl;
        
        row++;
    }while(std::getline(myfile, str) && row < Npoints);
    
    if(row != Npoints){
      std::cerr << "Number of data rows in " << filename << " does not match expected number of particles."
      << std::endl;
      throw std::runtime_error("file reading error");
    }
  }else{
    std::cerr << "Unable to open file " << filename << std::endl;
    throw std::runtime_error("file reading error");
  }
  
  std::cout << Npoints << " particle positions read from file " << filename << std::endl;
};
 
template<typename DType>
void LensHaloParticles<DType>::rotate(PosType theta_x,PosType theta_y
      ,Point_3d<DType> &xo){
 
  if(theta_x == 0.0 && theta_y == 0.0) return;
  
  PosType coord[3][3];
  PosType cx,cy,sx,sy;
  
  cx = cos(theta_x); sx = sin(theta_x);
  cy = cos(theta_y); sy = sin(theta_y);
  
  coord[0][0] = cy;  coord[1][0] = -sy*sx; coord[2][0] = cx;
  coord[0][1] = 0;   coord[1][1] = cx;     coord[2][1] = sx;
  coord[0][2] = -sy; coord[1][2] = -cy*sx; coord[2][2] = cy*cx;
  
  Point_3d<DType> tmp;
      // rotate particle positions 
  for(auto &p : pp){
    tmp *= 0;
    for(int j=0;j<3;++j){
      double tmp2 = p[j] - xo[j];
      tmp[0] += coord[0][j]*tmp2;
      tmp[1] += coord[1][j]*tmp2;
      tmp[2] += coord[2][j]*tmp2;
    }
    for(int j=0;j<3;++j) p[j]=tmp[j] + xo[j];
  }
 
      // rotate center of mass
  tmp *= 0.0;
  for(int j=0;j<3;++j){
    tmp[0] += coord[0][j]*(mcenter[j] - xo[j]);
    tmp[1] += coord[1][j]*(mcenter[j] - xo[j]);
    tmp[2] += coord[2][j]*(mcenter[j] - xo[j]);
  }
  mcenter = tmp + xo;
 
  otree.reset( new OTreeNB<Point_3d<DType> >(pp.data(),pp.size()) );
  otree->build(Nbucket);
};
 
template<typename DType>
void LensHaloParticles<DType>::translate_particles(Point_3d<DType> &xo){
 
  if(xo[0] == 0.0 && xo[1] == 0.0) return;
  
      // translate particle positions 
  for(auto &p : pp){
    for(int j=0;j<3;++j) p[j] += xo[j];
  }
 
  mcenter += xo;
 
  otree.reset( new OTreeNB<Point_3d<DType> >(pp.data(),pp.size()) );
  otree->build(Nbucket);
};
 
template<typename DType = float >
class LensHaloParticlesM : public LensHaloParticles<DType>
{
public:
 
  LensHaloParticlesM(
      std::vector<Point_3d<DType> > &pvector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,std::vector<DType> Masses   
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<DType>(pvector,redshift,my_inv_area_phys,1,cosmo,my_Nsmooth
      ,Nbucket,theta,recenter,verbose)
    {
      std::swap(masses,Masses);
      this->otree->calcMoments(masses.data());
    };
 
  LensHaloParticlesM(LensHaloParticlesM &&h)
        : LensHaloParticles<DType>(std::move(h)),
          masses(std::move(h.masses)) 
  {
  }
  
  LensHaloParticlesM & operator=(LensHaloParticlesM &&h) {
    if (this == &h) return *this;
    LensHaloParticles<DType>::operator=(std::move(h));
    masses = std::move(h.masses);
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    this->otree->force2D(xcm,masses.data(),this->Nsmooth,this->theta2,this->inv_area
        ,alpha,kappa,gamma,phi);
  };
  void swapoutMasses(std::vector<DType> &Mvector){
    std::swap(masses, Mvector);
    this->otree->calcMoments(masses.data());
  }
 private:
  std::vector<DType> masses;
};
 
 
template<typename DType = float >
class LensHaloParticlesO : public LensHaloParticles<DType>
{
public:
  LensHaloParticlesO(
    const std::string& simulation_filename 
    ,SimFileFormat format    
    ,PosType redshift        
    ,double my_inv_area_phys      
    ,PosType mass_particle   
    ,Point_2d x_hole         
    ,DType hole_radius       
    ,const COSMOLOGY& cosmo  
    ,int my_Nsmooth = 5      
    ,float Nbucket = 8  
    ,float theta = 0.1      
    ,bool recenter = false 
    ,bool verbose = false
  ) : LensHaloParticles<DType>(simulation_filename,format,redshift,my_inv_area_phys,mass_particle,cosmo,my_Nsmooth,Nbucket,theta,recenter,verbose)
  ,xc(x_hole),rhole(hole_radius)
  {};
 
  LensHaloParticlesO(
      std::vector<Point_3d<DType> > &pvector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,PosType mass_particle   
      ,Point_2d x_hole         
      ,DType hole_radius       
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<DType>(pvector,redshift,my_inv_area_phys,mass_particle,cosmo,my_Nsmooth
      ,Nbucket,theta,recenter,verbose),xc(x_hole),rhole(hole_radius)
    {
    };
 
  LensHaloParticlesO(LensHaloParticlesO &&h)
        : LensHaloParticles<DType>(std::move(h)),
          xc(std::move(h.xc)),
          rhole(h.rhole)
  {}
  
  LensHaloParticlesO & operator=(LensHaloParticlesO &&h) {
    if (this == &h) return *this;
    LensHaloParticles<DType>::operator=(std::move(h));
    xc = std::move(h.xc);
    rhole = h.rhole;
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    this->otree->force2D_hole(xcm,this->particle_mass,this->Nsmooth,this->theta2,this->inv_area
        ,rhole,xc.data()
        ,alpha,kappa,gamma,phi);
  };
 
  void setHole(
    Point_2d &x_hole,DType hole_radius
  ){
    xc = x_hole;
    rhole = hole_radius;
  };
  
 private:
  Point_2d xc;
  PosType rhole;
};
 
template<typename DType = float >
class LensHaloParticlesP : public LensHaloParticles<DType>
{
public:
  LensHaloParticlesP(
    const std::string& simulation_filename 
    ,SimFileFormat format    
    ,PosType redshift        
    ,double my_inv_area_phys      
    ,PosType mass_particle   
    ,Point_2d lower_left     
    ,Point_2d upper_right    
    ,const COSMOLOGY& cosmo  
    ,int my_Nsmooth = 5      
    ,float Nbucket = 8  
    ,float theta = 0.1      
    ,bool recenter = false 
    ,bool verbose = false
  ) : LensHaloParticles<DType>(simulation_filename,format,redshift,my_inv_area_phys,mass_particle,cosmo,my_Nsmooth,Nbucket,theta,recenter,verbose)
  ,length(upper_right-lower_left)
  {// check that particles are within bounding box
      Point_3d<float> p1,p2;
      this->otree->getBoundingBox(p1,p2);
      double Dl = cosmo.angDist(redshift);
      length *= Dl;
      p1 /= Dl;
      p2 /= Dl;
      if(  lower_left[0] > p1[0] 
        || lower_left[1] > p1[1] 
        || upper_right[0] < p2[0] 
        || upper_right[1] < p2[1] 
      ){
        
        std::cerr << "Error : LensParticlesP bounding box does not contain all the particles" 
        << std::endl;
        std::cerr << "   angular field " << lower_left << " " << upper_right << std::endl;
        std::cerr << "   bounding box " << p1 << " " << p2 << std::endl;
        throw std::runtime_error("Bad bounding box."); 
      }
  };
 
  LensHaloParticlesP(
      std::vector<Point_3d<DType> > &pvector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,PosType mass_particle   
      ,Point_2d lower_left     
      ,Point_2d upper_right   
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<DType>(pvector,redshift,my_inv_area_phys,mass_particle,cosmo,my_Nsmooth
      ,Nbucket,theta,recenter,verbose),length(upper_right-lower_left)
    {
      // check that particles are within bounding box
      Point_3d<float> p1,p2;
      this->otree->getBoundingBox(p1,p2);
      double Dl = cosmo.angDist(redshift);
      length *= Dl;
      p1 /= Dl;
      p2 /= Dl;
      if(  lower_left[0] > p1[0] 
        || lower_left[1] > p1[1] 
        || upper_right[0] < p2[0] 
        || upper_right[1] < p2[1] 
      ){
        
        std::cerr << "Error : LensParticlesP bounding box does not contain all the particles" 
        << std::endl;
        std::cerr << "   angular field " << lower_left << " " << upper_right << std::endl;
        std::cerr << "   bounding box " << p1 << " " << p2 << std::endl;
        throw std::runtime_error("Bad bounding box."); 
      }
    };
 
  LensHaloParticlesP(LensHaloParticlesP &&h)
        : LensHaloParticles<DType>(std::move(h)),
          length(h.length)
  {}
  
  LensHaloParticlesP & operator=(LensHaloParticlesP &&h) {
    if (this == &h) return *this;
    LensHaloParticles<DType>::operator=(std::move(h));
    length = h.length;
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    Point_2d xx;
    for(int i = -1 ; i <= 1 ; ++i ){
      xx[0] = xcm[0] + i*length[0];
      for(int j = -1 ; j <= 1 ; ++j ){
        xx[1] = xcm[1] + j*length[1];
        this->otree->force2D(xx.x,this->particle_mass,this->Nsmooth,this->theta2,this->inv_area
        ,alpha,kappa,gamma,phi);
      }
    }
  };
 
 private:
  Point_2d length;
};
 
template<typename DType = float >
class LensHaloParticlesMO : public LensHaloParticles<DType>
{
public:
 
  LensHaloParticlesMO(
      std::vector<Point_3d<DType> > &pvector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,std::vector<DType> Masses   
      ,Point_2d x_hole         
      ,DType hole_radius       
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 5      
      ,float Nbucket = 8       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<DType>(pvector,redshift,my_inv_area_phys,1,cosmo,my_Nsmooth
      ,Nbucket,theta,recenter,verbose),xc(x_hole),rhole(hole_radius)
    {
      std::swap(masses,Masses);
      this->otree->calcMoments(masses.data());
    };
 
  LensHaloParticlesMO(LensHaloParticlesMO &&h)
        : LensHaloParticles<DType>(std::move(h)),
          xc(std::move(h.xc)),
          rhole(h.rhole),
          masses(std::move(h.masses)) 
  {
    
  }
  
  LensHaloParticlesMO & operator=(LensHaloParticlesMO &&h) {
    if (this == &h) return *this;
    LensHaloParticles<DType>::operator=(std::move(h));
    xc = std::move(h.xc);
    rhole = h.rhole;
    masses = std::move(h.masses);
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    this->otree->force2D_hole(xcm,masses.data(),this->Nsmooth,this->theta2,this->inv_area
        ,rhole,xc.data()
        ,alpha,kappa,gamma,phi);
  };
 
  void setHole(
    Point_2d &x_hole,DType hole_radius
  ){
    xc[0] = x_hole[0];
    xc[1] = x_hole[1];
    rhole = hole_radius;
  };
 
 private:
  Point_2d xc;
  PosType rhole;
  std::vector<DType> masses;
};
 
template<typename HType >
class LensHaloH : public LensHaloParticles<float>
{
public:
 
  LensHaloH(
      std::vector<HType> &halo_vector 
      ,PosType redshift        
      ,double my_inv_area_phys      
      ,float hole_radius_angle 
      ,const COSMOLOGY& cosmo  
      ,float Nbucket = 4       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<float>(redshift,my_inv_area_phys,1,cosmo,0 
      ,Nbucket,theta,recenter,verbose),rhole(hole_radius_angle*cosmo.angDist(redshift))
    {
      std::swap(halo_vector,halos);
      size_t N = halos.size();
      this->pp.resize(N);
      double Dl = cosmo.angDist(redshift);
      Point_2d p;
      for(size_t i=0 ; i<N ; ++i){
        p = halos[i].getTheta();
        this->pp[i][0] = p[0]*Dl;
        this->pp[i][1] = p[1]*Dl;
        this->pp[i][2] = 0;
        halos[i].setRsize(1.0e3);
      }
 
      this->setUp(recenter,verbose);
      this->otree->calcMoments(halos);
    };
 
  LensHaloH(LensHaloH &&h)
        : LensHaloParticles<float>(std::move(h)),
          rhole(h.rhole),
          halos(std::move(h.halos)) 
  {
  }
  
  LensHaloH & operator=(LensHaloH &&h) {
    if (this == &h) return *this;
    LensHaloParticles<float>::operator=(std::move(h));
    rhole = h.rhole;
    halos = std::move(h.halos);
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    Point_2d xx;
    this->otree->force2D_halo(xcm,halos.data()
    ,this->theta2,this->inv_area
        ,rhole,xcm 
        ,alpha,kappa,gamma,phi);
  };
 
  void swapoutHalos(std::vector<HType> &hvector){
    std::swap(halos, hvector);
    this->otree.reset();
  }
 private:
  PosType rhole;
  std::vector<HType> halos;
};
 
template<typename HType >
class LensHaloHP : public LensHaloParticles<float>
{
public:
 
  LensHaloHP(
      std::vector<HType> &halo_vector 
      ,PosType redshift        
      ,Point_2d lower_left     
      ,Point_2d upper_right   
      ,float hole_radius_angle 
      ,const COSMOLOGY& cosmo  
      ,int my_Nsmooth = 0      
      ,float Nbucket = 4       
      ,float theta = 0.1       
      ,bool recenter = false   
      ,bool verbose = false
    ) : LensHaloParticles<float>(redshift
      ,0,1,cosmo,my_Nsmooth
      ,Nbucket,theta,recenter,verbose),rhole(hole_radius_angle*cosmo.angDist(redshift))
      ,length(upper_right-lower_left)
    {
      std::swap(halo_vector,halos);
      size_t N = halos.size();
      this->pp.resize(N);
      double Dl = cosmo.angDist(redshift);
      
      Point_2d p;
      for(size_t i=0 ; i<N ; ++i){
        p = halos[i].getTheta();
        this->pp[i][0] = p[0]*Dl;
        this->pp[i][1] = p[1]*Dl;
        this->pp[i][2] = 0;
        halos[i].setRsize(1.0e3);
      }
      length *= Dl;
      LensHaloParticles<float>::inv_area = 1.0/length[0]/length[1];
      this->setUp(recenter,verbose);
      this->otree->calcMoments(halos);
 
      // check that particles are within bounding box
      Point_3d<float> p1,p2;
      this->otree->getBoundingBox(p1,p2);
      p1 /= Dl;
      p2 /= Dl;
      if(  lower_left[0] > p1[0] 
        || lower_left[1] > p1[1] 
        || upper_right[0] < p2[0] 
        || upper_right[1] < p2[1] 
      ){
        
        std::cerr << "Error : LensHaloHP bounding box does not contain all the halos" 
        << std::endl;
        std::cerr << "   angular field " << lower_left << " " << upper_right << std::endl;
        std::cerr << "   bounding box " << p1 << " " << p2 << std::endl;
        throw std::runtime_error("Bad bounding box."); 
      }
    };
 
  LensHaloHP(LensHaloHP &&h)
        : LensHaloParticles<float>(std::move(h)),
          rhole(h.rhole),
          halos(std::move(h.halos)),
          length(h.length)
  {
  }
  
  LensHaloHP & operator=(LensHaloHP &&h) {
    if (this == &h) return *this;
    LensHaloParticles<float>::operator=(std::move(h));
    rhole = h.rhole;
    halos = std::move(h.halos);
    length = h.length;
    return *this;
  }
  void force_halo(double *alpha
        ,KappaType *kappa
        ,KappaType *gamma
        ,KappaType *phi
        ,double const *xcm 
        ,bool subtract_point=false
        ,PosType screening = 1     // here so that it overrides the LensHalo::force_halo                               
  ){
 
    Point_2d xx,zero(0,0);
    for(int i = -1 ; i <= 1 ; ++i ){
      xx[0] = xcm[0] + i*length[0];
      for(int j = -1 ; j <= 1 ; ++j ){
        xx[1] = xcm[1] + j*length[1];;
        //this->otree->force2D_halo(xx.x,halos.data(),this->Nsmooth,this->theta2,this->inv_area
        //  ,rhole,alpha,kappa,gamma,phi);
        this->otree->force2D_halo(xx.x,halos.data(),this->theta2,this->inv_area
           ,rhole,xx.x
           ,alpha,kappa,gamma,phi);
      }
    }
  };
 
 private:
  PosType rhole;
  std::vector<HType> halos;
  Point_2d length;
};
 
#endif