Newer
Older
FairShipTools / lookAtGeo.py
import ROOT,os,sys,getopt,time
import shipunit as u
import shipRoot_conf
from ShipGeoConfig import ConfigRegistry
from ROOT import *

# function to prepare two dictionaries: one that has the node names as key and the related intex 
# in the nodes array as value and the other one with key and value swapped. 
def prepareNodeDictionaries(nodes):
    # key: index, value: name of the node
    nodeDict_index = {}
    # key: name of the node, value: index
    nodeDict_name = {}
    for (i,n) in enumerate(nodes):
        #print i,n.GetName()
        nodeDict_index[i]=n.GetName()
        nodeDict_name[n.GetName()]=i
    return {'nodeDict_index':nodeDict_index, 'nodeDict_name':nodeDict_name}

# function to print the node names and their index. Usefull to search the name of the element you would like
# to study if you don't know it by heart. Then the name can be used with the other functions for the specific studies.
def searchForNodes(inputFile, volName = None):
    r = loadGeometry(inputFile)
    fGeo = r['fGeo']
    ## Get the top volume
    fGeo = ROOT.gGeoManager
    
    if volName is None:
        tv = fGeo.GetTopVolume()
    else: 
         tv = fGeo.GetVolume(volName)
         
    nodes = tv.GetNodes()
    for (i,n) in enumerate(nodes):
        print i,n.GetName()

def searchForNodes2(inputFile,volName=None):
    r = loadGeometry(inputFile)
    fGeo = r['fGeo']
    ## Get the top volume
    #fGeo = ROOT.gGeoManager
    if volName is None:
        tv = fGeo.GetTopVolume()
    else:
        tv = fGeo.GetVolume(volName)
        
    nodes = tv.GetNodes()
    for (i,n) in enumerate(nodes):
        print i,n.GetName(),findPositionElement(n)['z'],findDimentionBoxElement(n)['z']

def searchForNodes2_xyz(inputFile, volName=None):
    r = loadGeometry(inputFile)
    fGeo = r['fGeo']
    ## Get the top volume
    #fGeo = ROOT.gGeoManager
    if volName is None:
        tv = fGeo.GetTopVolume()
    else:
        tv = fGeo.GetVolume(volName)
    nodes = tv.GetNodes()
    for (i,n) in enumerate(nodes):
        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']
        
# basic function to be called to load the geometry from a file        
def loadGeometry(inputFile):
    dy = float( inputFile.split("/")[-1].split(".")[1])

    ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy )
    # init geometry and mag. field
    ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy )
    # -----Create geometry----------------------------------------------
    import shipDet_conf

    tgeom = ROOT.TGeoManager("Geometry", "Geane geometry")
    geofile = inputFile.replace('ship.','geofile_full.').replace('_rec.','.')
    print geofile
    gMan  = tgeom.Import(geofile)
    fGeo = ROOT.gGeoManager
    return {'fGeo':fGeo,'gMan':gMan, 'ShipGeo':ShipGeo}

def getNode(nodeName,fGeo=None, volName=None):
    if fGeo is None:
        r = loadGeometry(inputFile)
        fGeo = r['fGeo']
    if volName is None:
        tv = fGeo.GetTopVolume()
    else:
        tv = fGeo.GetVolume(volName)
        
    nodes = tv.GetNodes()
    dicts = prepareNodeDictionaries(nodes)
    return nodes[dicts['nodeDict_name'][nodeName]]
    
def findDimentionBoxElement(node):
    ## GetDZ() etc gives you half of the dimention
    sh = node.GetVolume().GetShape()
    return {'x':sh.GetDX(),
            'y':sh.GetDY(),
            'z':sh.GetDZ()}

def findPositionElement(node):
    pos = node.GetMatrix().GetTranslation()
    
    return {'x':pos[0],
            'y':pos[1],
            'z':pos[2]}
    
## inputFile: file used from which I would like to retrive the geometry used
## myNodes_name: geometrical elements I would like to analyse (find the position)
def findPositionGeoElement(inputFile, myNodes_name, volName=None):
 
    r = loadGeometry(inputFile)
    fGeo = r['fGeo']
    #volumes = gMan.GetListOfVolumes()
    #for v in volumes:
    #    print v.GetName(), v.GetNumber()
    
    ## Get the top volume
    if volName is None:
        tv = fGeo.GetTopVolume()
    else:
        tv = fGeo.GetVolume(volName)
        
    nodes = tv.GetNodes()

    tmp = prepareNodeDictionaries(nodes)
    nodeDict_name = tmp['nodeDict_name']
    res = {}
    for nd_name in myNodes_name:
        nd_index = nodeDict_name[nd_name]
        nd = nodes[nd_index]
        pos = findPositionElement(nd)
        dims = findDimentionBoxElement(nd)
        res[nd_name]={'z':pos['z'],
                      'dimZ':dims['z'],
                      'node':nd,
                      'dimX':dims['x'],
                      'dimY':dims['y']}
    return res

#inputFile = "../data/neutrino661/ship.10.0.Genie-TGeant4_D.root"
##searchForNodes(inputFile)

### nodes name to be looked at
#myNodes_name = ["volLayer2_%s"%i for i in xrange(0,12)]
#myGeoEl = findPositionGeoElement(inputFile, myNodes_name)
#print myGeoEl

##print myGeoEl["volLayer2_0"]['z'],myGeoEl["volLayer2_11"]['z']
##print myGeoEl["volLayer2_11"]['z']-myGeoEl["volLayer2_0"]['z']