Ask Your Question

Revision history [back]

click to hide/show revision 1
initial version

The functionality of EtaProduct seems to be restricted to creating modular functions. For the purposes of the posted question there is no implemented support. (At least as far as my quick search failed to find any.)

In the following, there will be some lines of code, written and dedicated for the case $$ S_6(\ \Gamma_0(24)\ )\ . $$ (The question was a good occasion for me to learn about the subject, i had to do something in a neighbour room, and was even extermely curious to see the results for $S_6(\Gamma_0(24))$.)

(The code is written, so that it could be reused, at least theoretically.)

As a background for the implementation, it is simplest for me to use the setting (definitions, notations, and results) from [Günther Köhler, Eta Products and Theta Series Identities, Springer Monographs in Mathematics, Springer, 2011], in the sequel [GK2011], for short. For the convenience of the reader, i will try to also illustrate the intermediate steps by the structure found so far.

ONE.

Chapter 2 in [GK2011], conditions for holomorphic eta products, is telling us which data can be used to build relevant eta products.
It is convenient to consider an eta product associated to a symbolic sum of the shape $$ x = \sum_{d|24}x(d)[d]\in{\mathbb Z}[D] $$ in the free abelian group, i.e. free ${\mathbb Z}$--module with the basis $[d]$, $d$ being a pure symbol running in the set $D$ of divisors of $24$. (There will be no $N=24$, for me it is better to fix ideas and use the $24$ explicitly.)
(The additive notation is better suited, the multiplicative notation $\prod_{d|24}[d]^{x(d)}$ represents small the important data $(x(d))_d$, as an upper index.)

Associated to $x$ there is the following eta product: $$ \eta(x)(z)=\prod{d|24}\eta(\; dz\;)^{x(d)}\ . $$ The conditions of holomorphicity at the cusps (of the form $1/c$, where $c$ divides $24$,) can be restated as $$ 24\sum_{d|24}\frac{\text{gcd}(c,d)^2}{cd}\cdot x(d)\ge 0\ , $$ for all $c$ with $c|24$. (It is This isolates a symmetric matrix with entries $$ 24\frac{\text{gcd}(c,d)^2}{cd}\ ,\c,d\in D\ , $$ which are integers.

We fix an order on $D$. In the following code, it corresponds to the entries of the $8\times 1$ matrix d24. (The number $24$ has euler_phi(24)$=8$ divisors.)

The above entries are then arranged corresponding to this order and build a matrix $B=$B24.

We search column vectors $x$ (with entries $x(d)$) so that the holomorphicity conditions are satisfied: $$ Bx\ge 0\ . $$ (The $0$ is the $8\times 1$ zero matrix and the inequality is considered componentwise.) This isolates a cone in ${\mathbb Z}^8$.

The eta function behaves like a modular form of weight $1/2$ for some multiplier system that can be technically described. So $\eta(x)$ behaves like a modular form of weight $$ \frac 12 \sum_{d|24} x(d)\ . $$ In our case, we are searching modular (and even cusp) forms of weight $6$, so the condition is $$ 12 = \sum_{d|24} x(d)\ . $$ We intersect the above cone $Bx\ge 0$ with this hyperplane. The resulted simplex has "known extremal points" $v$, that can be explicitly described for "the same problem", when we replace $24$ by its highest prime power divisors $8$, and respectively $3$. Such extremal points $v$ build correspondingly (quadratic) matrices V8 and V3 with euler_phi(8)$=4$, and respectively euler_phi(3)$=2$ rows and columns. The following routine dBV computes the data $d,B,V$ (order of the divisors, matrix $B$ describing the cone simplex $Bx\ge 0$, and some scaling to integers of the extremal points of the cone $V$). We will soon have:

sage: dBV(8)
(
[1]  [8 4 2 1]  [ 2 -2  0  0]
[2]  [4 8 4 2]  [-1  5 -2  0]
[4]  [2 4 8 4]  [ 0 -2  5 -1]
[8], [1 2 4 8], [ 0  0 -2  2]
)
sage: dBV(3)
(
[1]  [3 1]  [ 3/2 -1/2]
[3], [1 3], [-1/2  3/2]
)

The third matrix is in both cases so that the sum of entries in each column is one.

Then the data for $24$ is collected by building the tensor product of the above matrices. The results are block matrices:

sage: dBV(24)
(
[ 1]  [24  8|12  4| 6  2| 3  1]  [   3   -1|  -3    1|   0    0|   0    0]
[ 3]  [ 8 24| 4 12| 2  6| 1  3]  [  -1    3|   1   -3|   0    0|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 2]  [12  4|24  8|12  4| 6  2]  [-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 6]  [ 4 12| 8 24| 4 12| 2  6]  [ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 4]  [ 6  2|12  4|24  8|12  4]  [   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[12]  [ 2  6| 4 12| 8 24| 4 12]  [   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 8]  [ 3  1| 6  2|12  4|24  8]  [   0    0|   0    0|  -3    1|   3   -1]
[24], [ 1  3| 2  6| 4 12| 8 24], [   0    0|   0    0|   1   -3|  -1    3]
)

(The print shows d24, the chosen order of the divisors of $24$, the matrix B24, and the columns of the next block matrix V24 are extremal points for the cone.) The product of the last two matrices is:

sage: d24, B24, V24 = dBV(24)
sage: B24 * V24
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

An eta product, leading to a (holomorphic) modular form of weight $6$ is thus obtained by taking a linear combination with coefficients (building a vector) $t\ge 0$ of the columns of V24, such that the scalars used in it sum to $12=2\cdot 6$. (Since each column of V24 has entries adding to one.)

We will get an $x$ of the shape $x=Vt$.

Let us type some code first, that collects the matrices $d, B,V$, this will fix (for me) the frame:

def dBV( N ):
    """deliver three matrices, d, B, and V,
    d records the divisors of N, computed in a specific order,
    that will reflect the same order used for B, V.
    B implements the / a positivity condition "pos",
    a solution of "pos" is a column x with Bx >= 0,
    and V is the list of extremal points (as columns) for "pos" whith rational coordinates,
    and such that each column sum is one.

    This will help us to construct eta products as follows.
    The components of d, considered as eta argument multiplicators, 
    and the corresponding components of a solution x
    (x has components in ZZ), considered as corresponding eta powers,
    describe a relevant holomorphic (at cusps) eta product.

    It is a cusp form, iff Bx > 0.
    """
    if N.is_prime_power():
        p, ppower = N.is_prime_power( get_data=True )    # so N = p^ppower

        nv = ppower+1
        d  = matrix( QQ, nv,  1 )
        b  = matrix( QQ, nv, nv )
        v  = matrix( QQ, nv, nv )

        for j in range(nv):
            d[j, 0] = p^j
            for k in range(nv):
                b[j, k] = p^( ppower + 2*min(j,k) - j - k )

        for j in range(nv):
            if   j == 0:
                v[0, 0] = p*(p-1)
                v[1, 0] =    1-p
            elif j == nv-1:
                v[j-1, j] =    1-p
                v[j  , j] = p*(p-1)
            else:
                v[j-1, j] = -p
                v[j  , j] = 1+p^2
                v[j+1, j] = -p
        return d, b, (1/(p-1)^2) * v

    else:
        prime_powers = [ p^r for p, r in N.factor() ]
        d = matrix( QQ, 1, 1, [1,] )
        v = matrix( QQ, 1, 1, [1,] )
        b = matrix( QQ, 1, 1, [1,] )

        for q in prime_powers:
            dq, bq, vq = dBV( q )
            d = d.tensor_product( dq )
            b = b.tensor_product( bq )
            v = v.tensor_product( vq )

        return d, b, v

d24, B24, V24 = dBV( 24 )
dd24 = matrix( QQ, 8, 1, [ 24/d24[k, 0] for k in range(8) ] )

print "d24 is:\n%s\n"  % d24
print "dd24 is:\n%s\n" % dd24
print "B24 is:\n%s\n"  % B24
print "V24 is:\n%s\n"  % V24
print "B24 * V24 is:\n%s\n" % ( B24*V24 )

""" The above calculates, as already mentioned:

d24 is:
[ 1]
[ 3]
[--]
[ 2]
[ 6]
[--]
[ 4]
[12]
[--]
[ 8]
[24]

dd24 is:
[24]
[ 8]
[12]
[ 4]
[ 6]
[ 2]
[ 3]
[ 1]

B24 is:
[24  8|12  4| 6  2| 3  1]
[ 8 24| 4 12| 2  6| 1  3]
[-----+-----+-----+-----]
[12  4|24  8|12  4| 6  2]
[ 4 12| 8 24| 4 12| 2  6]
[-----+-----+-----+-----]
[ 6  2|12  4|24  8|12  4]
[ 2  6| 4 12| 8 24| 4 12]
[-----+-----+-----+-----]
[ 3  1| 6  2|12  4|24  8]
[ 1  3| 2  6| 4 12| 8 24]

V24 is:
[   3   -1|  -3    1|   0    0|   0    0]
[  -1    3|   1   -3|   0    0|   0    0]
[---------+---------+---------+---------]
[-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[---------+---------+---------+---------]
[   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[---------+---------+---------+---------]
[   0    0|   0    0|  -3    1|   3   -1]
[   0    0|   0    0|   1   -3|  -1    3]

B24 * V24 is:
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

TWO:

We want moreover a cusp form resulting from the construction, so for $Bx>0$, $x=Vt$, we need correspondingly $t>0$. In order to get some runnable code, some examples, i will start with the Ansatz, that t has natural entries $\ge 1$. (So there will be no component $1/2$ or $1/3$ or...) This facilitates the enumeration. Then the eight entries of $t$ are at least one each, they sum to $12$, so we have four further points to distribute among the components.

The following simple code collects the candidates:

etaDataList = []

t0 = matrix( QQ, 8, 1, [ 1 for _ in range(8) ] )
tList = []
C = [ [c0, c1, c2, c3]
      for c0 in [ 0..7 ]
      for c1 in [ c0..7 ]
      for c2 in [ c1..7 ]
      for c3 in [ c2..7 ] ]

for c in C:
    t = t0 + matrix( QQ, 8, 1, [ c.count( _ ) for _ in range(8) ] )
    x = V24 * t

    if not all( [ entry in ZZ for entry in (  d24.transpose() * t ) . list() ] ):    continue
    # if not all( [ entry in ZZ for entry in ( dd24.transpose() * t ) . list() ] ):    continue
    if not all( [ entry in ZZ for entry in (                    x ) . list() ] ):    continue

    tList.append( t )
    etaDataList.append( ( [ ZZ(entry) for entry in x.list() ], d24 ) )

(Decommenting the commented condition leads to the same list.) Here, len( etaDataList ) and len( tList ) are both $170$. So we have $170$ candidates to play with.

(We do not need the tList, and also need from etaDataList only the data comming from the $x$-component. But for $170$ entries this is not expensive.)

THREE:

We are in the position to construct some eta products. Let us do this with one entry in the collected etaDataList, for instance with

ed  = etaDataList[ 11 ]

This is

sage: ed
(
                              [ 1]
                              [ 3]
                              [--]
                              [ 2]
                              [ 6]
                              [--]
                              [ 4]
                              [12]
                              [--]
                              [ 8]
[6, -2, -3, 1, 7, 7, -2, -2], [24]
)

Associated to it we have the linear combination $$

\sum x(d)[d]

6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]\ . $$ The eta product for it is $$

\eta(x)(z)

\eta(\; z\; )^6 \eta(\; 3z\; )^{-2} \eta(\; 2z\; )^{-3} \eta(\; 6z\; ) \eta(\; 4z\; )^7 \eta(\; 12z\; )^7 \eta(\; 8z\; )^{-2} \eta(\; 24z\; )^{-2}\ . $$ The corresponding $q$-product representation, $q = \exp( \; 2\pi i\; z\; )$, has $d\cdot z\leftrightarrow q^d$, and $\eta(z)\leftrightarrow q^{1/24}f(q)$, where $f$ "is" f from

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

which gives

sage: f
1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - q^35 - q^40 + q^51 + q^57 - q^70 - q^77 + q^92 + O(q^100)

and $\eta(d\cdot z)\leftrightarrow q^{d/24}f(q^d)$, and in the eta product we have to cumulate the power-of-$q$-prefactor. To the above element $x=6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]$ there corresponds...

ed = etaDataList[ 11 ]
x, d = ed

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

fx = PSR( prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) \
           * \
           q^ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 ) )
fx

This gives a "long" $q$-expansion, and we can / have to truncate it here to fit somehow:

sage: fx + O(q^10)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + O(q^10)

"Is" this element in the space of cusp forms $S_6(\Gamma_0(24))$, at least if we test this for the first $100$ coefficients in the $q$-expansion?! For this, we associate the space CF = CuspForms( Gamma0(24), 6 ) of cusp forms, each cusp form (in the basis of this space) is mapped to its $q$--expansion, truncated modulo $O(q^100)$, then mapped to $\mathbb Q^{100}$ by picking the coefficients with respect to $1,q,q^2,\dots, q^{99}$ modulo $O(q^{100})$ .

This is a vector space $W=$WCF of dimension $16$, let us construct it:

CF = CuspForms( Gamma0(24), 6 ) 
# CF is:
# Cuspidal subspace of dimension 16
# of Modular Forms space of dimension 24
# for Congruence Subgroup Gamma0(24)
# of weight 6 over Rational Field

def fit_to_length( myList, length ):
    lenMyList = len( myList )

    if   lenMyList == length:
        return myList

    elif lenMyList  < length:
        return myList + ( length - lenMyList ) * [0,]

    else:
        return myList[:length]

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

Then the coordinates of the eta product $\eta(x)$ in the sage basis of CF are as follows:

Wf = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wf )
print coord

This gives:

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

In particular, the constructed eta product lives in the space of cusp forms CF$=S_6(\; \Gamma_0(24)\; )$.

FOUR:

Let us print the first coefficients, say working modulo $O(q^{21})$ in the CF basis and in the eta product $\eta(x)$. The $16$ elements in the basis are (starting with)

for cf in CF.basis():
    print cf.qexp( 21 )

This gives:

q + 48*q^19 + O(q^21)
q^2 + 39*q^18 + O(q^21)
q^3 + 30*q^19 + O(q^21)
q^4 + 4*q^16 - 30*q^20 + O(q^21)
q^5 - 20*q^19 + O(q^21)
q^6 - 12*q^18 + O(q^21)
q^7 - 7*q^19 + O(q^21)
q^8 - 6*q^16 + O(q^21)
q^9 - 20*q^19 + O(q^21)
q^10 - 15*q^18 + O(q^21)
(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)
q^11 - 11*q^19 + O(q^21)
q^12 - 4*q^16 + 4*q^20 + O(q^21)
q^13 - 7*q^19 + O(q^21)
q^14 - 6*q^18 + O(q^21)
q^15 - 5*q^19 + O(q^21)
q^17 - q^19 + O(q^21)

And of course, from the $q$--expansion of $\eta(x)$,

sage: fx + O(q^21)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + 70*q^10 
    - 102*q^11 + 108*q^12 - 12*q^13 - 200*q^14 + 114*q^15 + 48*q^16 - 108*q^17 
    + 345*q^18 - 210*q^19 + 72*q^20 + O(q^21)

(Output was manually rearranged.)

The coefficients in

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

are exactly the $16$ coefficients of the $q$--expansion in degrees $1,2,\dots, 14,15,\ 17$. (Which are the "(strating) orders" for the basis elements of CF.)

FIVE:

Let us see that all $170$ eta products found in TWO (variable etaDataList) are also in CF, at least regarding the test of the first $100$ coefficients.

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

count = 0
for ed in etaDataList:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s -> %s" % ( count, str(x), str(coord) )

This gives a long list of coordinates for the found eta products with respect to the standard base of CF...

The first ten entries, reshaped in the print are:

count = 0
for ed in etaDataList[:10]:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s\n    %s" % ( count, str(x), str(coord) )

This gives:

  1 x=[12, -4, -4, 4, 2, 2, 0, 0]   
    (0, 1, -12, 58, -132, 81, 276, -668, 468, 294, -876, 810, 24, -904, 324, -936)
  2 x=[8, 0, -2, 2, 2, 2, 0, 0]     
    (0, 1, -8, 22, -16, -27, 40, -20, 96, -138, -8, 54, -128, 392, -216, -192)
  3 x=[6, -2, 5, 1, -1, 3, 0, 0]    
    (0, 1, -6, 4, 42, -81, -78, 304, -90, -354, 426, -324, -204, 1040, -486, 180)
  4 x=[10, -6, -5, 11, 3, -1, 0, 0] 
    (0, 1, -10, 40, -74, 27, 158, -344, 282, 78, -442, 432, -52, -256, 54, -564)
  5 x=[4, 4, 0, 0, 2, 2, 0, 0]      
    (0, 1, -4, 2, 4, 9, -4, -44, 12, 6, 92, 18, -88, -40, -180, -24)
  6 x=[2, 2, 7, -1, -1, 3, 0, 0]    
    (0, 1, -2, -8, 14, 27, -26, -56, -30, 78, 142, 0, -68, -256, -162, 60)
  7 x=[6, -2, -3, 9, 3, -1, 0, 0]   
    (0, 1, -6, 12, -6, -9, 18, -32, 54, -66, 42, 36, -108, 176, -198, -108)
  8 x=[0, 0, 14, -2, -4, 4, 0, 0]   
    (0, 1, 0, -14, 0, 81, 0, -236, 0, 294, 0, 162, 0, -904, 0, 0)
  9 x=[4, -4, 4, 8, 0, 0, 0, 0]     
    (0, 1, -4, -2, 28, -27, -52, 124, -60, -138, 284, -162, -136, 392, -324, 120)
 10 x=[8, -8, -6, 18, 4, -4, 0, 0]  
    (0, 1, -8, 26, -40, 9, 88, -188, 168, 6, -200, 234, -80, -40, -72, -336)

The functionality of EtaProduct seems to be restricted to creating modular functions. For the purposes of the posted question there is no implemented support. (At least as far as my quick search failed to find any.)

In the following, there will be some lines of code, written and dedicated for the case $$ S_6(\ \Gamma_0(24)\ )\ . $$ (The question was a good occasion for me to learn about the subject, i had to do something in a neighbour room, and was even extermely curious to see the results for $S_6(\Gamma_0(24))$.)

(The code is written, so that it could be reused, at least theoretically.)

As a background for the implementation, it is simplest for me to use the setting (definitions, notations, and results) from [Günther Köhler, Eta Products and Theta Series Identities, Springer Monographs in Mathematics, Springer, 2011], in the sequel [GK2011], for short. For the convenience of the reader, i will try to also illustrate the intermediate steps by the structure found so far.

ONE.

Chapter 2 in [GK2011], conditions for holomorphic eta products, is telling us which data can be used to build relevant eta products.
It is convenient to consider an eta product associated to a symbolic sum of the shape $$ x = \sum_{d|24}x(d)[d]\in{\mathbb Z}[D] $$ in the free abelian group, i.e. free ${\mathbb Z}$--module with the basis $[d]$, $d$ being a pure symbol running in the set $D$ of divisors of $24$. (There will be no $N=24$, for me it is better to fix ideas and use the $24$ explicitly.)
(The additive notation is better suited, the multiplicative notation $\prod_{d|24}[d]^{x(d)}$ represents small the important data $(x(d))_d$, as an upper index.)

Associated to $x$ there is the following eta product: $$ \eta(x)(z)=\prod{d|24}\eta(\; dz\;)^{x(d)}\ . $$ The conditions of holomorphicity at the cusps (of the form $1/c$, where $c$ divides $24$,) can be restated as $$ 24\sum_{d|24}\frac{\text{gcd}(c,d)^2}{cd}\cdot x(d)\ge 0\ , $$ for all $c$ with $c|24$. (It is This isolates a symmetric matrix with entries $$ 24\frac{\text{gcd}(c,d)^2}{cd}\ ,\c,d\in ,\ c,d\in D\ , $$ which are integers.

We fix an order on $D$. In the following code, it corresponds to the entries of the $8\times 1$ matrix d24. (The number $24$ has euler_phi(24)$=8$ divisors.)

The above entries are then arranged corresponding to this order and build a matrix $B=$B24.

We search column vectors $x$ (with entries $x(d)$) so that the holomorphicity conditions are satisfied: $$ Bx\ge 0\ . $$ (The $0$ is the $8\times 1$ zero matrix and the inequality is considered componentwise.) This isolates a cone in ${\mathbb Z}^8$.

The eta function behaves like a modular form of weight $1/2$ for some multiplier system that can be technically described. So $\eta(x)$ behaves like a modular form of weight $$ \frac 12 \sum_{d|24} x(d)\ . $$ In our case, we are searching modular (and even cusp) forms of weight $6$, so the condition is $$ 12 = \sum_{d|24} x(d)\ . $$ We intersect the above cone $Bx\ge 0$ with this hyperplane. The resulted simplex has "known extremal points" $v$, that can be explicitly described for "the same problem", when we replace $24$ by its highest prime power divisors $8$, and respectively $3$. Such extremal points $v$ build correspondingly (quadratic) matrices V8 and V3 with euler_phi(8)$=4$, and respectively euler_phi(3)$=2$ rows and columns. The following routine dBV computes the data $d,B,V$ (order of the divisors, matrix $B$ describing the cone simplex $Bx\ge 0$, and some scaling to integers of the extremal points of the cone $V$). We will soon have:

sage: dBV(8)
(
[1]  [8 4 2 1]  [ 2 -2  0  0]
[2]  [4 8 4 2]  [-1  5 -2  0]
[4]  [2 4 8 4]  [ 0 -2  5 -1]
[8], [1 2 4 8], [ 0  0 -2  2]
)
sage: dBV(3)
(
[1]  [3 1]  [ 3/2 -1/2]
[3], [1 3], [-1/2  3/2]
)

The third matrix is in both cases so that the sum of entries in each column is one.

Then the data for $24$ is collected by building the tensor product of the above matrices. The results are block matrices:

sage: dBV(24)
(
[ 1]  [24  8|12  4| 6  2| 3  1]  [   3   -1|  -3    1|   0    0|   0    0]
[ 3]  [ 8 24| 4 12| 2  6| 1  3]  [  -1    3|   1   -3|   0    0|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 2]  [12  4|24  8|12  4| 6  2]  [-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 6]  [ 4 12| 8 24| 4 12| 2  6]  [ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 4]  [ 6  2|12  4|24  8|12  4]  [   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[12]  [ 2  6| 4 12| 8 24| 4 12]  [   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 8]  [ 3  1| 6  2|12  4|24  8]  [   0    0|   0    0|  -3    1|   3   -1]
[24], [ 1  3| 2  6| 4 12| 8 24], [   0    0|   0    0|   1   -3|  -1    3]
)

(The print shows d24, the chosen order of the divisors of $24$, the matrix B24, and the columns of the next block matrix V24 are extremal points for the cone.) The product of the last two matrices is:

sage: d24, B24, V24 = dBV(24)
sage: B24 * V24
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

An eta product, leading to a (holomorphic) modular form of weight $6$ is thus obtained by taking a linear combination with coefficients (building a vector) $t\ge 0$ of the columns of V24, such that the scalars used in it sum to $12=2\cdot 6$. (Since each column of V24 has entries adding to one.)

We will get an $x$ of the shape $x=Vt$.

Let us type some code first, that collects the matrices $d, B,V$, this will fix (for me) the frame:

def dBV( N ):
    """deliver three matrices, d, B, and V,
    d records the divisors of N, computed in a specific order,
    that will reflect the same order used for B, V.
    B implements the / a positivity condition "pos",
    a solution of "pos" is a column x with Bx >= 0,
    and V is the list of extremal points (as columns) for "pos" whith rational coordinates,
    and such that each column sum is one.

    This will help us to construct eta products as follows.
    The components of d, considered as eta argument multiplicators, 
    and the corresponding components of a solution x
    (x has components in ZZ), considered as corresponding eta powers,
    describe a relevant holomorphic (at cusps) eta product.

    It is a cusp form, iff Bx > 0.
    """
    if N.is_prime_power():
        p, ppower = N.is_prime_power( get_data=True )    # so N = p^ppower

        nv = ppower+1
        d  = matrix( QQ, nv,  1 )
        b  = matrix( QQ, nv, nv )
        v  = matrix( QQ, nv, nv )

        for j in range(nv):
            d[j, 0] = p^j
            for k in range(nv):
                b[j, k] = p^( ppower + 2*min(j,k) - j - k )

        for j in range(nv):
            if   j == 0:
                v[0, 0] = p*(p-1)
                v[1, 0] =    1-p
            elif j == nv-1:
                v[j-1, j] =    1-p
                v[j  , j] = p*(p-1)
            else:
                v[j-1, j] = -p
                v[j  , j] = 1+p^2
                v[j+1, j] = -p
        return d, b, (1/(p-1)^2) * v

    else:
        prime_powers = [ p^r for p, r in N.factor() ]
        d = matrix( QQ, 1, 1, [1,] )
        v = matrix( QQ, 1, 1, [1,] )
        b = matrix( QQ, 1, 1, [1,] )

        for q in prime_powers:
            dq, bq, vq = dBV( q )
            d = d.tensor_product( dq )
            b = b.tensor_product( bq )
            v = v.tensor_product( vq )

        return d, b, v

d24, B24, V24 = dBV( 24 )
dd24 = matrix( QQ, 8, 1, [ 24/d24[k, 0] for k in range(8) ] )

print "d24 is:\n%s\n"  % d24
print "dd24 is:\n%s\n" % dd24
print "B24 is:\n%s\n"  % B24
print "V24 is:\n%s\n"  % V24
print "B24 * V24 is:\n%s\n" % ( B24*V24 )

""" The above calculates, as already mentioned:

d24 is:
[ 1]
[ 3]
[--]
[ 2]
[ 6]
[--]
[ 4]
[12]
[--]
[ 8]
[24]

dd24 is:
[24]
[ 8]
[12]
[ 4]
[ 6]
[ 2]
[ 3]
[ 1]

B24 is:
[24  8|12  4| 6  2| 3  1]
[ 8 24| 4 12| 2  6| 1  3]
[-----+-----+-----+-----]
[12  4|24  8|12  4| 6  2]
[ 4 12| 8 24| 4 12| 2  6]
[-----+-----+-----+-----]
[ 6  2|12  4|24  8|12  4]
[ 2  6| 4 12| 8 24| 4 12]
[-----+-----+-----+-----]
[ 3  1| 6  2|12  4|24  8]
[ 1  3| 2  6| 4 12| 8 24]

V24 is:
[   3   -1|  -3    1|   0    0|   0    0]
[  -1    3|   1   -3|   0    0|   0    0]
[---------+---------+---------+---------]
[-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[---------+---------+---------+---------]
[   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[---------+---------+---------+---------]
[   0    0|   0    0|  -3    1|   3   -1]
[   0    0|   0    0|   1   -3|  -1    3]

B24 * V24 is:
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

TWO:

We want moreover a cusp form resulting from the construction, so for $Bx>0$, $x=Vt$, we need correspondingly $t>0$. In order to get some runnable code, some examples, i will start with the Ansatz, that t has natural entries $\ge 1$. (So there will be no component $1/2$ or $1/3$ or...) This facilitates the enumeration. Then the eight entries of $t$ are at least one each, they sum to $12$, so we have four further points to distribute among the components.

The following simple code collects the candidates:

etaDataList = []

t0 = matrix( QQ, 8, 1, [ 1 for _ in range(8) ] )
tList = []
C = [ [c0, c1, c2, c3]
      for c0 in [ 0..7 ]
      for c1 in [ c0..7 ]
      for c2 in [ c1..7 ]
      for c3 in [ c2..7 ] ]

for c in C:
    t = t0 + matrix( QQ, 8, 1, [ c.count( _ ) for _ in range(8) ] )
    x = V24 * t

    if not all( [ entry in ZZ for entry in (  d24.transpose() * t ) . list() ] ):    continue
    # if not all( [ entry in ZZ for entry in ( dd24.transpose() * t ) . list() ] ):    continue
    if not all( [ entry in ZZ for entry in (                    x ) . list() ] ):    continue

    tList.append( t )
    etaDataList.append( ( [ ZZ(entry) for entry in x.list() ], d24 ) )

(Decommenting the commented condition leads to the same list.) Here, len( etaDataList ) and len( tList ) are both $170$. So we have $170$ candidates to play with.

(We do not need the tList, and also need from etaDataList only the data comming from the $x$-component. But for $170$ entries this is not expensive.)

THREE:

We are in the position to construct some eta products. Let us do this with one entry in the collected etaDataList, for instance with

ed  = etaDataList[ 11 ]

This is

sage: ed
(
                              [ 1]
                              [ 3]
                              [--]
                              [ 2]
                              [ 6]
                              [--]
                              [ 4]
                              [12]
                              [--]
                              [ 8]
[6, -2, -3, 1, 7, 7, -2, -2], [24]
)

Associated to it we have the linear combination $$

$$ \sum x(d)[d]

x(d)[d] = 6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]\ . $$ The eta product for it is $$

\eta(x)(z)

$$ \eta(x)(z) = \eta(\; z\; )^6 \eta(\; 3z\; )^{-2} \eta(\; 2z\; )^{-3} \eta(\; 6z\; ) \eta(\; 4z\; )^7 \eta(\; 12z\; )^7 \eta(\; 8z\; )^{-2} \eta(\; 24z\; )^{-2}\ . $$ The corresponding $q$-product representation, $q = \exp( \; 2\pi i\; z\; )$, has $d\cdot z\leftrightarrow q^d$, and $\eta(z)\leftrightarrow q^{1/24}f(q)$, where $f$ "is" f from

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

which gives

sage: f
1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - q^35 - q^40 + q^51 + q^57 - q^70 - q^77 + q^92 + O(q^100)

and $\eta(d\cdot z)\leftrightarrow q^{d/24}f(q^d)$, and in the eta product we have to cumulate the power-of-$q$-prefactor. To the above element $x=6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]$ there corresponds...

ed = etaDataList[ 11 ]
x, d = ed

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

fx = PSR( prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) \
           * \
           q^ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 ) )
fx

This gives a "long" $q$-expansion, and we can / have to truncate it here to fit somehow:

sage: fx + O(q^10)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + O(q^10)

"Is" this element in the space of cusp forms $S_6(\Gamma_0(24))$, at least if we test this for the first $100$ coefficients in the $q$-expansion?! For this, we associate the space CF = CuspForms( Gamma0(24), 6 ) of cusp forms, each cusp form (in the basis of this space) is mapped to its $q$--expansion, truncated modulo $O(q^100)$, then mapped to $\mathbb Q^{100}$ by picking the coefficients with respect to $1,q,q^2,\dots, q^{99}$ modulo $O(q^{100})$ .

This is a vector space $W=$WCF of dimension $16$, let us construct it:

CF = CuspForms( Gamma0(24), 6 ) 
# CF is:
# Cuspidal subspace of dimension 16
# of Modular Forms space of dimension 24
# for Congruence Subgroup Gamma0(24)
# of weight 6 over Rational Field

def fit_to_length( myList, length ):
    lenMyList = len( myList )

    if   lenMyList == length:
        return myList

    elif lenMyList  < length:
        return myList + ( length - lenMyList ) * [0,]

    else:
        return myList[:length]

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

Then the coordinates of the eta product $\eta(x)$ in the sage basis of CF are as follows:

Wf = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wf )
print coord

This gives:

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

In particular, the constructed eta product lives in the space of cusp forms CF$=S_6(\; \Gamma_0(24)\; )$.

FOUR:

Let us print the first coefficients, say working modulo $O(q^{21})$ in the CF basis and in the eta product $\eta(x)$. The $16$ elements in the basis are (starting with)

for cf in CF.basis():
    print cf.qexp( 21 )

This gives:

q + 48*q^19 + O(q^21)
q^2 + 39*q^18 + O(q^21)
q^3 + 30*q^19 + O(q^21)
q^4 + 4*q^16 - 30*q^20 + O(q^21)
q^5 - 20*q^19 + O(q^21)
q^6 - 12*q^18 + O(q^21)
q^7 - 7*q^19 + O(q^21)
q^8 - 6*q^16 + O(q^21)
q^9 - 20*q^19 + O(q^21)
q^10 - 15*q^18 + O(q^21)
(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)
q^11 - 11*q^19 + O(q^21)
q^12 - 4*q^16 + 4*q^20 + O(q^21)
q^13 - 7*q^19 + O(q^21)
q^14 - 6*q^18 + O(q^21)
q^15 - 5*q^19 + O(q^21)
q^17 - q^19 + O(q^21)

And of course, from the $q$--expansion of $\eta(x)$,

sage: fx + O(q^21)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + 70*q^10 
    - 102*q^11 + 108*q^12 - 12*q^13 - 200*q^14 + 114*q^15 + 48*q^16 - 108*q^17 
    + 345*q^18 - 210*q^19 + 72*q^20 + O(q^21)

(Output was manually rearranged.)

The coefficients in

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

are exactly the $16$ coefficients of the $q$--expansion in degrees $1,2,\dots, 14,15,\ 17$. (Which are the "(strating) orders" for the basis elements of CF.)

FIVE:

Let us see that all $170$ eta products found in TWO (variable etaDataList) are also in CF, at least regarding the test of the first $100$ coefficients.

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

count = 0
for ed in etaDataList:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s -> %s" % ( count, str(x), str(coord) )

This gives a long list of coordinates for the found eta products with respect to the standard base of CF...

The first ten entries, reshaped in the print are:

count = 0
for ed in etaDataList[:10]:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s\n    %s" % ( count, str(x), str(coord) )

This gives:

  1 x=[12, -4, -4, 4, 2, 2, 0, 0]   
    (0, 1, -12, 58, -132, 81, 276, -668, 468, 294, -876, 810, 24, -904, 324, -936)
  2 x=[8, 0, -2, 2, 2, 2, 0, 0]     
    (0, 1, -8, 22, -16, -27, 40, -20, 96, -138, -8, 54, -128, 392, -216, -192)
  3 x=[6, -2, 5, 1, -1, 3, 0, 0]    
    (0, 1, -6, 4, 42, -81, -78, 304, -90, -354, 426, -324, -204, 1040, -486, 180)
  4 x=[10, -6, -5, 11, 3, -1, 0, 0] 
    (0, 1, -10, 40, -74, 27, 158, -344, 282, 78, -442, 432, -52, -256, 54, -564)
  5 x=[4, 4, 0, 0, 2, 2, 0, 0]      
    (0, 1, -4, 2, 4, 9, -4, -44, 12, 6, 92, 18, -88, -40, -180, -24)
  6 x=[2, 2, 7, -1, -1, 3, 0, 0]    
    (0, 1, -2, -8, 14, 27, -26, -56, -30, 78, 142, 0, -68, -256, -162, 60)
  7 x=[6, -2, -3, 9, 3, -1, 0, 0]   
    (0, 1, -6, 12, -6, -9, 18, -32, 54, -66, 42, 36, -108, 176, -198, -108)
  8 x=[0, 0, 14, -2, -4, 4, 0, 0]   
    (0, 1, 0, -14, 0, 81, 0, -236, 0, 294, 0, 162, 0, -904, 0, 0)
  9 x=[4, -4, 4, 8, 0, 0, 0, 0]     
    (0, 1, -4, -2, 28, -27, -52, 124, -60, -138, 284, -162, -136, 392, -324, 120)
 10 x=[8, -8, -6, 18, 4, -4, 0, 0]  
    (0, 1, -8, 26, -40, 9, 88, -188, 168, 6, -200, 234, -80, -40, -72, -336)

The functionality of EtaProduct seems to be restricted to creating modular functions. For the purposes of the posted question there is no implemented support. (At least as far as my quick search failed to find any.)

In the following, there will be some lines of code, written and dedicated for the case $$ S_6(\ \Gamma_0(24)\ )\ . $$ (The question was a good occasion for me to learn about the subject, i had to do something in a neighbour room, and was even extermely curious to see the results for $S_6(\Gamma_0(24))$.)

(The code is written, so that it could be reused, at least theoretically.)

As a background for the implementation, it is simplest for me to use the setting (definitions, notations, and results) from [Günther Köhler, Eta Products and Theta Series Identities, Springer Monographs in Mathematics, Springer, 2011], in the sequel [GK2011], for short. For the convenience of the reader, i will try to also illustrate the intermediate steps by the structure found so far.

ONE.

Chapter 2 in [GK2011], conditions for holomorphic eta products, is telling us which data can be used to build relevant eta products.
It is convenient to consider an eta product associated to a symbolic sum of the shape $$ x = \sum_{d|24}x(d)[d]\in{\mathbb Z}[D] $$ in the free abelian group, i.e. free ${\mathbb Z}$--module with the basis $[d]$, $d$ being a pure symbol running in the set $D$ of divisors of $24$. (There will be no $N=24$, for me it is better to fix ideas and use the $24$ explicitly.)
(The additive notation is better suited, the multiplicative notation $\prod_{d|24}[d]^{x(d)}$ represents small the important data $(x(d))_d$, as an upper index.)

Associated to $x$ there is the following eta product: $$ \eta(x)(z)=\prod{d|24}\eta(\; dz\;)^{x(d)}\ . $$ The conditions of holomorphicity at the cusps (of the form $1/c$, where $c$ divides $24$,) can be restated as $$ 24\sum_{d|24}\frac{\text{gcd}(c,d)^2}{cd}\cdot x(d)\ge 0\ , $$ for all $c$ with $c|24$. (It is This isolates a symmetric matrix with entries $$ 24\frac{\text{gcd}(c,d)^2}{cd}\ ,\ c,d\in D\ , $$ which are integers.

We fix an order on $D$. In the following code, it corresponds to the entries of the $8\times 1$ matrix d24. (The number $24$ has euler_phi(24)$=8$ divisors.)

The above entries are then arranged corresponding to this order and build a matrix $B=$B24.

We search column vectors $x$ (with entries $x(d)$) so that the holomorphicity conditions are satisfied: $$ Bx\ge 0\ . $$ (The $0$ is the $8\times 1$ zero matrix and the inequality is considered componentwise.) This isolates a cone in ${\mathbb Z}^8$.

The eta function behaves like a modular form of weight $1/2$ for some multiplier system that can be technically described. So $\eta(x)$ behaves like a modular form of weight $$ \frac 12 \sum_{d|24} x(d)\ . $$ In our case, we are searching modular (and even cusp) forms of weight $6$, so the condition is $$ 12 = \sum_{d|24} x(d)\ . $$ We intersect the above cone $Bx\ge 0$ with this hyperplane. The resulted simplex has "known extremal points" $v$, that can be explicitly described for "the same problem", when we replace $24$ by its highest prime power divisors $8$, and respectively $3$. Such extremal points $v$ build correspondingly (quadratic) matrices V8 and V3 with euler_phi(8)$=4$, and respectively euler_phi(3)$=2$ rows and columns. The following routine dBV computes the data $d,B,V$ (order of the divisors, matrix $B$ describing the cone simplex $Bx\ge 0$, and some scaling to integers of the extremal points of the cone $V$). We will soon have:

sage: dBV(8)
(
[1]  [8 4 2 1]  [ 2 -2  0  0]
[2]  [4 8 4 2]  [-1  5 -2  0]
[4]  [2 4 8 4]  [ 0 -2  5 -1]
[8], [1 2 4 8], [ 0  0 -2  2]
)
sage: dBV(3)
(
[1]  [3 1]  [ 3/2 -1/2]
[3], [1 3], [-1/2  3/2]
)

The third matrix is in both cases so that the sum of entries in each column is one.

Then the data for $24$ is collected by building the tensor product of the above matrices. The results are block matrices:

sage: dBV(24)
(
[ 1]  [24  8|12  4| 6  2| 3  1]  [   3   -1|  -3    1|   0    0|   0    0]
[ 3]  [ 8 24| 4 12| 2  6| 1  3]  [  -1    3|   1   -3|   0    0|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 2]  [12  4|24  8|12  4| 6  2]  [-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 6]  [ 4 12| 8 24| 4 12| 2  6]  [ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 4]  [ 6  2|12  4|24  8|12  4]  [   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[12]  [ 2  6| 4 12| 8 24| 4 12]  [   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[--]  [-----+-----+-----+-----]  [---------+---------+---------+---------]
[ 8]  [ 3  1| 6  2|12  4|24  8]  [   0    0|   0    0|  -3    1|   3   -1]
[24], [ 1  3| 2  6| 4 12| 8 24], [   0    0|   0    0|   1   -3|  -1    3]
)

(The print shows d24, the chosen order of the divisors of $24$, the matrix B24, and the columns of the next block matrix V24 are extremal points for the cone.) The product of the last two matrices is:

sage: d24, B24, V24 = dBV(24)
sage: B24 * V24
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

An eta product, leading to a (holomorphic) modular form of weight $6$ is thus obtained by taking a linear combination with coefficients (building a vector) $t\ge 0$ of the columns of V24, such that the scalars used in it sum to $12=2\cdot 6$. (Since each column of V24 has entries adding to one.)

We will get an $x$ of the shape $x=Vt$.

Let us type some code first, that collects the matrices $d, B,V$, this will fix (for me) the frame:

def dBV( N ):
    """deliver three matrices, d, B, and V,
    d records the divisors of N, computed in a specific order,
    that will reflect the same order used for B, V.
    B implements the / a positivity condition "pos",
    a solution of "pos" is a column x with Bx >= 0,
    and V is the list of extremal points (as columns) for "pos" whith rational coordinates,
    and such that each column sum is one.

    This will help us to construct eta products as follows.
    The components of d, considered as eta argument multiplicators, 
    and the corresponding components of a solution x
    (x has components in ZZ), considered as corresponding eta powers,
    describe a relevant holomorphic (at cusps) eta product.

    It is a cusp form, iff Bx > 0.
    """
    if N.is_prime_power():
        p, ppower = N.is_prime_power( get_data=True )    # so N = p^ppower

        nv = ppower+1
        d  = matrix( QQ, nv,  1 )
        b  = matrix( QQ, nv, nv )
        v  = matrix( QQ, nv, nv )

        for j in range(nv):
            d[j, 0] = p^j
            for k in range(nv):
                b[j, k] = p^( ppower + 2*min(j,k) - j - k )

        for j in range(nv):
            if   j == 0:
                v[0, 0] = p*(p-1)
                v[1, 0] =    1-p
            elif j == nv-1:
                v[j-1, j] =    1-p
                v[j  , j] = p*(p-1)
            else:
                v[j-1, j] = -p
                v[j  , j] = 1+p^2
                v[j+1, j] = -p
        return d, b, (1/(p-1)^2) * v

    else:
        prime_powers = [ p^r for p, r in N.factor() ]
        d = matrix( QQ, 1, 1, [1,] )
        v = matrix( QQ, 1, 1, [1,] )
        b = matrix( QQ, 1, 1, [1,] )

        for q in prime_powers:
            dq, bq, vq = dBV( q )
            d = d.tensor_product( dq )
            b = b.tensor_product( bq )
            v = v.tensor_product( vq )

        return d, b, v

d24, B24, V24 = dBV( 24 )
dd24 = matrix( QQ, 8, 1, [ 24/d24[k, 0] for k in range(8) ] )

print "d24 is:\n%s\n"  % d24
print "dd24 is:\n%s\n" % dd24
print "B24 is:\n%s\n"  % B24
print "V24 is:\n%s\n"  % V24
print "B24 * V24 is:\n%s\n" % ( B24*V24 )

The above calculates, as already mentioned:

d24 is:
[ 1]
[ 3]
[--]
[ 2]
[ 6]
[--]
[ 4]
[12]
[--]
[ 8]
[24]

dd24 is:
[24]
[ 8]
[12]
[ 4]
[ 6]
[ 2]
[ 3]
[ 1]

B24 is:
[24  8|12  4| 6  2| 3  1]
[ 8 24| 4 12| 2  6| 1  3]
[-----+-----+-----+-----]
[12  4|24  8|12  4| 6  2]
[ 4 12| 8 24| 4 12| 2  6]
[-----+-----+-----+-----]
[ 6  2|12  4|24  8|12  4]
[ 2  6| 4 12| 8 24| 4 12]
[-----+-----+-----+-----]
[ 3  1| 6  2|12  4|24  8]
[ 1  3| 2  6| 4 12| 8 24]

V24 is:
[   3   -1|  -3    1|   0    0|   0    0]
[  -1    3|   1   -3|   0    0|   0    0]
[---------+---------+---------+---------]
[-3/2  1/2|15/2 -5/2|  -3    1|   0    0]
[ 1/2 -3/2|-5/2 15/2|   1   -3|   0    0]
[---------+---------+---------+---------]
[   0    0|  -3    1|15/2 -5/2|-3/2  1/2]
[   0    0|   1   -3|-5/2 15/2| 1/2 -3/2]
[---------+---------+---------+---------]
[   0    0|   0    0|  -3    1|   3   -1]
[   0    0|   0    0|   1   -3|  -1    3]

B24 * V24 is:
[48  0  0  0  0  0  0  0]
[ 0 48  0  0  0  0  0  0]
[ 0  0 96  0  0  0  0  0]
[ 0  0  0 96  0  0  0  0]
[ 0  0  0  0 96  0  0  0]
[ 0  0  0  0  0 96  0  0]
[ 0  0  0  0  0  0 48  0]
[ 0  0  0  0  0  0  0 48]

TWO:

We want moreover a cusp form resulting from the construction, so for $Bx>0$, $x=Vt$, we need correspondingly $t>0$. In order to get some runnable code, some examples, i will start with the Ansatz, that t has natural entries $\ge 1$. (So there will be no component $1/2$ or $1/3$ or...) This facilitates the enumeration. Then the eight entries of $t$ are at least one each, they sum to $12$, so we have four further points to distribute among the components.

The following simple code collects the candidates:

etaDataList = []

t0 = matrix( QQ, 8, 1, [ 1 for _ in range(8) ] )
tList = []
C = [ [c0, c1, c2, c3]
      for c0 in [ 0..7 ]
      for c1 in [ c0..7 ]
      for c2 in [ c1..7 ]
      for c3 in [ c2..7 ] ]

for c in C:
    t = t0 + matrix( QQ, 8, 1, [ c.count( _ ) for _ in range(8) ] )
    x = V24 * t

    if not all( [ entry in ZZ for entry in (  d24.transpose() * t ) . list() ] ):    continue
    # if not all( [ entry in ZZ for entry in ( dd24.transpose() * t ) . list() ] ):    continue
    if not all( [ entry in ZZ for entry in (                    x ) . list() ] ):    continue

    tList.append( t )
    etaDataList.append( ( [ ZZ(entry) for entry in x.list() ], d24 ) )

(Decommenting the commented condition leads to the same list.) Here, len( etaDataList ) and len( tList ) are both $170$. So we have $170$ candidates to play with.

(We do not need the tList, and also need from etaDataList only the data comming from the $x$-component. But for $170$ entries this is not expensive.)

THREE:

We are in the position to construct some eta products. Let us do this with one entry in the collected etaDataList, for instance with

ed  = etaDataList[ 11 ]

This is

sage: ed
(
                              [ 1]
                              [ 3]
                              [--]
                              [ 2]
                              [ 6]
                              [--]
                              [ 4]
                              [12]
                              [--]
                              [ 8]
[6, -2, -3, 1, 7, 7, -2, -2], [24]
)

Associated to it we have the linear combination $$ \sum x(d)[d] = 6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]\ . $$ The eta product for it is $$ \eta(x)(z) = \eta(\; z\; )^6 \eta(\; 3z\; )^{-2} \eta(\; 2z\; )^{-3} \eta(\; 6z\; ) \eta(\; 4z\; )^7 \eta(\; 12z\; )^7 \eta(\; 8z\; )^{-2} \eta(\; 24z\; )^{-2}\ . $$ The corresponding $q$-product representation, $q = \exp( \; 2\pi i\; z\; )$, has $d\cdot z\leftrightarrow q^d$, and $\eta(z)\leftrightarrow q^{1/24}f(q)$, where $f$ "is" f from

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

which gives

sage: f
1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - q^35 - q^40 + q^51 + q^57 - q^70 - q^77 + q^92 + O(q^100)

and $\eta(d\cdot z)\leftrightarrow q^{d/24}f(q^d)$, and in the eta product we have to cumulate the power-of-$q$-prefactor. To the above element $x=6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24] \in{\mathbb Z}[D]$ there corresponds...

ed = etaDataList[ 11 ]
x, d = ed

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

fx = PSR( prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) \
           * \
           q^ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 ) )
fx

This gives a "long" $q$-expansion, and we can / have to truncate it here to fit somehow:

sage: fx + O(q^10)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + O(q^10)

"Is" this element in the space of cusp forms $S_6(\Gamma_0(24))$, at least if we test this for the first $100$ coefficients in the $q$-expansion?! For this, we associate the space CF = CuspForms( Gamma0(24), 6 ) of cusp forms, each cusp form (in the basis of this space) is mapped to its $q$--expansion, truncated modulo $O(q^100)$, then mapped to $\mathbb Q^{100}$ by picking the coefficients with respect to $1,q,q^2,\dots, q^{99}$ modulo $O(q^{100})$ .

This is a vector space $W=$WCF of dimension $16$, let us construct it:

CF = CuspForms( Gamma0(24), 6 ) 
# CF is:
# Cuspidal subspace of dimension 16
# of Modular Forms space of dimension 24
# for Congruence Subgroup Gamma0(24)
# of weight 6 over Rational Field

def fit_to_length( myList, length ):
    lenMyList = len( myList )

    if   lenMyList == length:
        return myList

    elif lenMyList  < length:
        return myList + ( length - lenMyList ) * [0,]

    else:
        return myList[:length]

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

Then the coordinates of the eta product $\eta(x)$ in the sage basis of CF are as follows:

Wf = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wf )
print coord

This gives:

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

In particular, the constructed eta product lives in the space of cusp forms CF$=S_6(\; \Gamma_0(24)\; )$.

FOUR:

Let us print the first coefficients, say working modulo $O(q^{21})$ in the CF basis and in the eta product $\eta(x)$. The $16$ elements in the basis are (starting with)

for cf in CF.basis():
    print cf.qexp( 21 )

This gives:

q + 48*q^19 + O(q^21)
q^2 + 39*q^18 + O(q^21)
q^3 + 30*q^19 + O(q^21)
q^4 + 4*q^16 - 30*q^20 + O(q^21)
q^5 - 20*q^19 + O(q^21)
q^6 - 12*q^18 + O(q^21)
q^7 - 7*q^19 + O(q^21)
q^8 - 6*q^16 + O(q^21)
q^9 - 20*q^19 + O(q^21)
q^10 - 15*q^18 + O(q^21)
(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)
q^11 - 11*q^19 + O(q^21)
q^12 - 4*q^16 + 4*q^20 + O(q^21)
q^13 - 7*q^19 + O(q^21)
q^14 - 6*q^18 + O(q^21)
q^15 - 5*q^19 + O(q^21)
q^17 - q^19 + O(q^21)

And of course, from the $q$--expansion of $\eta(x)$,

sage: fx + O(q^21)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + 70*q^10 
    - 102*q^11 + 108*q^12 - 12*q^13 - 200*q^14 + 114*q^15 + 48*q^16 - 108*q^17 
    + 345*q^18 - 210*q^19 + 72*q^20 + O(q^21)

(Output was manually rearranged.)

The coefficients in

(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)

are exactly the $16$ coefficients of the $q$--expansion in degrees $1,2,\dots, 14,15,\ 17$. (Which are the "(strating) orders" for the basis elements of CF.)

FIVE:

Let us see that all $170$ eta products found in TWO (variable etaDataList) are also in CF, at least regarding the test of the first $100$ coefficients.

PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )

VS = VectorSpace( QQ, 100 )
W  = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
                               for cf in CF.basis() ] )

count = 0
for ed in etaDataList:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s -> %s" % ( count, str(x), str(coord) )

This gives a long list of coordinates for the found eta products with respect to the standard base of CF...

The first ten entries, reshaped in the print are:

count = 0
for ed in etaDataList[:10]:
    count += 1
    x, d = ed
    qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
    fx  = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
    Wfx = VS( fit_to_length( fx.list(), 100 ) )
    coord = W.coordinate_vector( Wfx )
    print "%3s x=%-30s\n    %s" % ( count, str(x), str(coord) )

This gives:

  1 x=[12, -4, -4, 4, 2, 2, 0, 0]   
    (0, 1, -12, 58, -132, 81, 276, -668, 468, 294, -876, 810, 24, -904, 324, -936)
  2 x=[8, 0, -2, 2, 2, 2, 0, 0]     
    (0, 1, -8, 22, -16, -27, 40, -20, 96, -138, -8, 54, -128, 392, -216, -192)
  3 x=[6, -2, 5, 1, -1, 3, 0, 0]    
    (0, 1, -6, 4, 42, -81, -78, 304, -90, -354, 426, -324, -204, 1040, -486, 180)
  4 x=[10, -6, -5, 11, 3, -1, 0, 0] 
    (0, 1, -10, 40, -74, 27, 158, -344, 282, 78, -442, 432, -52, -256, 54, -564)
  5 x=[4, 4, 0, 0, 2, 2, 0, 0]      
    (0, 1, -4, 2, 4, 9, -4, -44, 12, 6, 92, 18, -88, -40, -180, -24)
  6 x=[2, 2, 7, -1, -1, 3, 0, 0]    
    (0, 1, -2, -8, 14, 27, -26, -56, -30, 78, 142, 0, -68, -256, -162, 60)
  7 x=[6, -2, -3, 9, 3, -1, 0, 0]   
    (0, 1, -6, 12, -6, -9, 18, -32, 54, -66, 42, 36, -108, 176, -198, -108)
  8 x=[0, 0, 14, -2, -4, 4, 0, 0]   
    (0, 1, 0, -14, 0, 81, 0, -236, 0, 294, 0, 162, 0, -904, 0, 0)
  9 x=[4, -4, 4, 8, 0, 0, 0, 0]     
    (0, 1, -4, -2, 28, -27, -52, 124, -60, -138, 284, -162, -136, 392, -324, 120)
 10 x=[8, -8, -6, 18, 4, -4, 0, 0]  
    (0, 1, -8, 26, -40, 9, 88, -188, 168, 6, -200, 234, -80, -40, -72, -336)