Starting point for finding the enveloping ellipsoid for a point cloud of calorimeter hits

perfectClustering.jl

###
# Compiles a cloud of calorimeter hits for each MCParticle
# Calorimeter cells with multiple hits are assigned to each particle, weighted with the correct energy contribution
###

using LCIO
using DataStructures
using PyPlot
using JuMP

# representation of an ellipsoid x^2/a^2 + y^2/b^2 + z^2/c^2 < r
type Ellipsoid
        x::Float64
        y::Float64
        z::Float64
        a::Float64
        b::Float64
        c::Float64
        r::Float64
end

function find_smallest_ellipsoid(cluster)
        function is_inside_ellipsoid(hit, e)
                if (hit.x-e.x)^2/e.a^2 + (hit.y-e.y)^2/e.b^2 + (hit.z-e.z)^2/e.c^2 < e.r
                        return true
                end
                return false
        end
end

for (idx, i) in enumerate(LCIO.open(ARGS[1]))
        for mcp in getCollection(i, "MCParticle")
                if getMCPGenStatus(mcp) != 1
                        continue
                end
                @printf("%s\t%.2f\n", getMCPDGid(mcp), getMCP4(mcp).t)
        end
        particleHits = DefaultDict(Ptr{Void}, Array, Array{LCIO.CalHit,1})
        for collectionName in ["EcalEndcapHits", "EcalBarrelHits"]#, "HcalEndcapHits", "HcalBarrelHits"]
                collection = getCollection(i, collectionName)
                println(collectionName, " has ", length(collection), " entries")
                for hit in collection
                        pList = getHitMCParticles(hit)
                        energyList = getHitEnergyList(hit)
                        calHit = getSimCaloHit(hit)
                        for (kdx, particle) in enumerate(pList)
                h = LCIO.CalHit(calHit, energyList[kdx])
                push!(particleHits[particle], h)
                end
                end
        end
        xvals = []
        yvals = []
        zvals = []
        for particle in particleHits
                @printf("%s\t%.2f\t%s\t%s\n", getMCPDGid(particle.first), getMCP4(particle.first).t, length(particle.second), sum([x.E for x in particle.second]))
                xs = []
                ys = []
                zs = []
                for hit in particle.second
                        push!(xs, hit.x)
                        push!(xvals, hit.x)
                        push!(ys, hit.y)
                        push!(yvals, hit.y)
                        push!(zs, hit.z)
                        push!(zvals, hit.z)
                end
                scatter3D(xs, ys, zs)
                show()
        end
        scatter3D(xvals, yvals, zvals)
        show()
end