import os

# Make debug ntuple
if not os.path.isfile("nuNutple_debug691.root"):
    os.system("python ntuple_nuVeto.py")

# For my output
if not os.path.exists("elena_plots"):
    os.system("mkdir elena_plots")

# Import amenities
from analyseNtupleForVeto import *

# Reminder
print ""
print "Interaction Elements:"
print "0\t Last passive element of Opera"
print "1\t Entrance window of the vacuum tank"
print "2\t Along the vacuum tank"
print "3\t Along Opera"
print "4\t Between windows (in the liquid scintillator)"
print ""

def findDimentionBoxElement(node):
    sh = node.GetVolume().GetShape()
    if sh.IsCylType():
        #print node.GetName(), " is a cylinder", sh.GetRmin(), sh.GetRmax()
        if sh.HasRmin():
            r = (sh.GetRmax() - sh.GetRmin())
        else:
            r = 0.
        return {'x':sh.GetDX(),
                'y':sh.GetDY(),
                'z':sh.GetDZ(),
                'r':r}
    return {'x':sh.GetDX(),
            'y':sh.GetDY(),
            'z':sh.GetDZ()}

def findPositionElement(node):
    pos = node.GetMatrix().GetTranslation()
    sh = node.GetVolume().GetShape()
    if sh.IsCylType():
        if sh.HasRmin():
            r = (sh.GetRmax() - sh.GetRmin())/2.
        else:
            r = sh.GetRmax()
        return {'x':pos[0],
                'y':pos[1],
                'z':pos[2],
                'r':r}
    return {'x':pos[0],
            'y':pos[1],
            'z':pos[2]}

def searchForNodes2_xyz_dict(inputFile, verbose=False):
    d = {}
    r = loadGeometry(inputFile)
    fGeo = r['fGeo']
    ## Get the top volume
    #fGeo = ROOT.gGeoManager
    tv = fGeo.GetTopVolume()
    nodes = tv.GetNodes()
    for (i,n) in enumerate(nodes):
        if verbose:
            print i,n.GetName()
            print "            x: ", findPositionElement(n)['x'],findDimentionBoxElement(n)['x']
            print "            y: ", findPositionElement(n)['y'],findDimentionBoxElement(n)['y']
            print "            z: ", findPositionElement(n)['z'],findDimentionBoxElement(n)['z']
        d[n.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}}
        d[n.GetName()]['x']['pos'] =     findPositionElement(n)['x']
        d[n.GetName()]['x']['dim'] = findDimentionBoxElement(n)['x']
        d[n.GetName()]['y']['pos'] =     findPositionElement(n)['y']
        d[n.GetName()]['y']['dim'] = findDimentionBoxElement(n)['y']
        d[n.GetName()]['z']['pos'] =     findPositionElement(n)['z']
        d[n.GetName()]['z']['dim'] = findDimentionBoxElement(n)['z']
        if n.GetVolume().GetShape().IsCylType():
            d[n.GetName()]['r']['pos'] =     findPositionElement(n)['r']
            d[n.GetName()]['r']['dim'] = findDimentionBoxElement(n)['r']
        else:
            d[n.GetName()]['r']['pos'] = 0.
            d[n.GetName()]['r']['dim'] = 0.
    return d

# Search nodes
nodesDict = searchForNodes2_xyz_dict('geofile_full.10.0.Genie-TGeant4.root')
# Sort nodes by z
from collections import OrderedDict
sortedNodes = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['pos']))

def printByZ():
    # Print nodes by increasing z
    print
    print
    print "Nodes by z:"
    for node in sortedNodes.keys():
        print node, '\t\t', (sortedNodes[node]['z']['pos']-sortedNodes[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodes[node]['z']['dim']/100., ' m'

# Sort nodes by dimension along z (the walls of the vacuum tank should be ~60 m long)
sortedNodesByZdim = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['dim']))

def printByZdim():
    # Print nodes by increasing dimension along z
    print
    print
    print "Nodes by z dimension:"
    for node in sortedNodesByZdim.keys():
        print node, '\t\t', (sortedNodesByZdim[node]['z']['pos']-sortedNodesByZdim[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodesByZdim[node]['z']['dim']/100., ' m'

# Custom sorting
def sortNodesBy(coord='z', what='pos', printdict=True):
    d = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1][coord][what]))
    if printdict:
        print
        print
        print "Nodes by %s %s:"%(coord, what)
        for node in d.keys():
            print node, '\t\t', (d[node][coord]['pos']-d[node][coord]['dim'])/100., ' m', '\t', 2*d[node][coord]['dim']/100., ' m'
    return d



# Load geometry (again?)
r = loadGeometry('geofile_full.10.0.Genie-TGeant4.root')
fGeo = r['fGeo']
#someOtherNodes = ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
#someOtherNodes += ["volLateralS1_d_1", "volLateralS2_d_1", "volLateralS1_c_1", "volLateralS2_c_1"]
#someOtherNodes += ["Target%s_1"%i for i in (0,1,2,3,5)]

# Make plots
from ROOT import *
cSaver = []
cSaver.append(TCanvas("entrance_windows","entrance_windows"))
t2.Draw("startY_nu:startX_nu","isPrimary==1 && (interactionElement==2 || interactionElement==4)")
# Veto
vetoIn = nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim']
vetoOut = nodesDict['Veto_5']['z']['pos'] + nodesDict['Veto_5']['z']['dim']
cSaver.append(TCanvas("veto","veto"))
t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(vetoIn,vetoOut))
# Decay volume
cSaver.append(TCanvas("decay_volume","decay_volume"))
volIn = nodesDict['T2Startplate_1']['z']['pos'] - nodesDict['T2Startplate_1']['z']['dim']
volOut = nodesDict['lidT6I_1']['z']['pos'] + nodesDict['lidT6I_1']['z']['dim']
t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(volIn,volOut))

# Save plots
for c in cSaver:
    c.Print('elena_plots/'+c.GetName()+'.eps')