Creating a pyramid around a fractal in VPython

Question

I am trying to create my own fractal using VPython. I currently have this hexagonal shape:

The formula of this fractal is

I have added some height to make it a sort of pyramid without sides:

I would like to add lines from the corners of the hexagons in one layer to corns in the layer beneath. This way I would really generate a pyramid.

Currently, this is my code:

from vpython import *
import math

def draw_hexagonal_snowflake(position, size, depth, height):
    if depth == 0 or size < 1e-6:
        return
    
    # Calculate the vertices of the hexagon
    vertices = []
    for i in range(6):
        theta = math.radians(60 * i)
        px = position.x + size * math.cos(theta)
        py = position.y + size * math.sin(theta)
        vertices.append(vector(px, py, height))
    
    # Draw the hexagon
    for i in range(6):
        if size > 1e-6:
            cylinder(pos=vertices[i], axis=vertices[(i+1)%6]-vertices[i], radius=size/20, color=color.white)
    
    # Recursive function for all sides of the hexagon
    for i in range(6):
        new_position = (vertices[i] + vertices[(i+1)%6]) / 2
        new_height = height - size * math.sin(math.radians(60))
        draw_hexagonal_snowflake(new_position, size / 2, depth - 1, new_height)

scene = canvas(width=800, height=800)

position = vector(0, 0, 0)  
size = 200  
depth = 5 # Amount of iterations
height = 0

draw_hexagonal_snowflake(position, size, depth, height)

I have already tried adding lines (cylinders) on the vertices of the hexagons, but I could only get the lines to be vertical:

Does anyone know how to go about this? Both formula edits or VPython code would be okay solutions for me.