from lookAtGeo import *
import tools
import shipunit as u
from ShipGeoConfig import ConfigRegistry
from operator import mul, add

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


def searchForNodes3_xyz_dict(fGeo, verbose=False):
    from tools import findPositionElement, findDimentionBoxElement, findPositionElement2
    d = {}
    #r = loadGeometry(inputFile)
    #fGeo = r['fGeo']
    ## Get the top volume
    #fGeo = ROOT.gGeoManager
    tv = fGeo.GetTopVolume()
    topnodes = tv.GetNodes()
    for (j,topn) in enumerate(topnodes):
        # top volumes
        if verbose:
            print j, topn.GetName()
            print "            x: ", findPositionElement(topn)['x'],findDimentionBoxElement(topn)['x']
            print "            y: ", findPositionElement(topn)['y'],findDimentionBoxElement(topn)['y']
            print "            z: ", findPositionElement(topn)['z'],findDimentionBoxElement(topn)['z']
        d[topn.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}}
        d[topn.GetName()]['x']['pos'] =     findPositionElement(topn)['x']
        d[topn.GetName()]['x']['dim'] = findDimentionBoxElement(topn)['x']
        d[topn.GetName()]['y']['pos'] =     findPositionElement(topn)['y']
        d[topn.GetName()]['y']['dim'] = findDimentionBoxElement(topn)['y']
        d[topn.GetName()]['z']['pos'] =     findPositionElement(topn)['z']
        d[topn.GetName()]['z']['dim'] = findDimentionBoxElement(topn)['z']
        if topn.GetVolume().GetShape().IsCylType():
            d[topn.GetName()]['r']['pos'] =     findPositionElement(topn)['r']
            d[topn.GetName()]['r']['dim'] = findDimentionBoxElement(topn)['r']
        else:
            d[topn.GetName()]['r']['pos'] = 0.
            d[topn.GetName()]['r']['dim'] = 0.
        # First children
        nodes = topn.GetNodes()
        if nodes:
            topPos = topn.GetMatrix().GetTranslation()
            for (i,n) in enumerate(nodes):
                if verbose:
                    print j, topn.GetName(), i, n.GetName()
                    print "            x: ", findPositionElement2(n,topPos)['x'],findDimentionBoxElement(n)['x']
                    print "            y: ", findPositionElement2(n,topPos)['y'],findDimentionBoxElement(n)['y']
                    print "            z: ", findPositionElement2(n,topPos)['z'],findDimentionBoxElement(n)['z']
                d[n.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}}
                d[n.GetName()]['x']['pos'] =     findPositionElement2(n,topPos)['x']
                d[n.GetName()]['x']['dim'] = findDimentionBoxElement(n)['x']
                d[n.GetName()]['y']['pos'] =     findPositionElement2(n,topPos)['y']
                d[n.GetName()]['y']['dim'] = findDimentionBoxElement(n)['y']
                d[n.GetName()]['z']['pos'] =     findPositionElement2(n,topPos)['z']
                d[n.GetName()]['z']['dim'] = findDimentionBoxElement(n)['z']
                if n.GetVolume().GetShape().IsCylType():
                    d[n.GetName()]['r']['pos'] =     findPositionElement2(n,topPos)['r']
                    d[n.GetName()]['r']['dim'] = findDimentionBoxElement(n)['r']
                else:
                    d[n.GetName()]['r']['pos'] = 0.
                    d[n.GetName()]['r']['dim'] = 0.
                # Second children
                cnodes = n.GetNodes()
                if cnodes:
                    localpos = n.GetMatrix().GetTranslation()
                    localToGlobal = []
                    for i in xrange(3):
                        localToGlobal.append(localpos[i]+topPos[i])
                    for (k,cn) in enumerate(cnodes):
                        if verbose:
                            print j, topn.GetName(), i, n.GetName(), k, cn.GetName()
                            print "            x: ", findPositionElement2(cn,localToGlobal)['x'],findDimentionBoxElement(cn)['x']
                            print "            y: ", findPositionElement2(cn,localToGlobal)['y'],findDimentionBoxElement(cn)['y']
                            print "            z: ", findPositionElement2(cn,localToGlobal)['z'],findDimentionBoxElement(cn)['z']
                        d[cn.GetName()] = {'x': {}, 'y':{}, 'z':{}, 'r':{}}
                        d[cn.GetName()]['x']['pos'] =     findPositionElement2(cn,localToGlobal)['x']
                        d[cn.GetName()]['x']['dim'] = findDimentionBoxElement(cn)['x']
                        d[cn.GetName()]['y']['pos'] =     findPositionElement2(cn,localToGlobal)['y']
                        d[cn.GetName()]['y']['dim'] = findDimentionBoxElement(cn)['y']
                        d[cn.GetName()]['z']['pos'] =     findPositionElement2(cn,localToGlobal)['z']
                        d[cn.GetName()]['z']['dim'] = findDimentionBoxElement(cn)['z']
                        if cn.GetVolume().GetShape().IsCylType():
                            d[cn.GetName()]['r']['pos'] =     findPositionElement2(cn,localToGlobal)['r']
                            d[cn.GetName()]['r']['dim'] = findDimentionBoxElement(cn)['r']
                        else:
                            d[cn.GetName()]['r']['pos'] = 0.
                            d[cn.GetName()]['r']['dim'] = 0.
    return d


ff = ROOT.TFile("histoForWeights.root")
h_GioHans = ff.Get("h_Gio")
def calcWeightNu(NC, E, w, entries, nuName, ON=True):
	# Only for neutrinos and antineutrinos
    if not ON:
        return 1
    if "bar" in nuName:
        reduction = 0.5
        flux = 6.98e+11 * 2.e+20 / 5.e+13
    else: 
        reduction = 1.
        flux = 1.09e+12 * 2.e+20/ 5.e+13
    crossSec = 6.7e-39*E * reduction
    NA = 6.022e+23
    binN = h_GioHans.GetXaxis().FindBin(E)    
    return crossSec * flux * h_GioHans.GetBinContent(binN) * w * NA / entries


def findWeight(sampleType, NC, E, MCTrack, entries, nuName, ON):
	if sampleType == 'nuBg':
		return calcWeightNu(NC, E, MCTrack.GetWeight(), entries, nuName, ON)
	elif sampleType == 'sig':
		return MCTrack.GetWeight() # for the acceptance, multiply by normalization
	elif sampleType == 'cosmics':
		return MCTrack.GetWeight() # multiply by 1.e6 / 200.



def hasStrawStations(event, trackId, listOfWantedStraws):
    ok = [False]*len(listOfWantedStraws)
    positions = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in listOfWantedStraws ]
    for hit in event.strawtubesPoint:
        if hit.GetTrackID() == trackId:
            for (i,det) in enumerate(listOfWantedStraws):
                if (positions[i][0] < hit.GetZ() < positions[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.):
                    ok[i] = True
    return bool(reduce(mul, ok, 1))

def hasGoodStrawStations(event, trackId):
    #ok = [False]*2
    okupstream = [False]*2
    okdownstream = [False]*2
    upstream = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in ['Tr1_1', 'Tr2_2'] ]
    downstream = [ (nodes[det]['z']['pos'] - nodes[det]['z']['dim'], nodes[det]['z']['pos'] + nodes[det]['z']['dim'] ) for det in ['Tr3_3', 'Tr4_4'] ]
    for hit in event.strawtubesPoint:
        if hit.GetTrackID() == trackId:
            for i in xrange(2):
                if (upstream[i][0] < hit.GetZ() < upstream[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.):
                    okupstream[i] = True
                if (downstream[i][0] < hit.GetZ() < downstream[i][1]) and tools.inEllipse(hit.GetX(), hit.GetY(), 250., 500.):
                    okdownstream[i] = True
    ok = [ bool(reduce(mul, l, 1)) for l in [okupstream, okdownstream] ]
    return bool(reduce(add, ok, 0))

def findHNLvertex(event):
    for t in event.MCTrack:
        if t.GetMotherId() == 1:
            return t.GetStartZ()
    return False

def has_muon_station(event, trackId, station):
	zIn = nodes['muondet%s_1'%(station-1)]['z']['pos'] - nodes['muondet%s_1'%(station-1)]['z']['dim']
	zOut = nodes['muondet%s_1'%(station-1)]['z']['pos'] + nodes['muondet%s_1'%(station-1)]['z']['dim']
	for hit in event.muonPoint:
		if hit.GetTrackID() == trackId:
			if zIn <= hit.GetZ() <= zOut:
				return True
	return False

def hasEcalDeposit(event, trackId, ELossThreshold):
	ELoss = 0.
	for hit in event.EcalPoint:
		if hit.GetTrackID() == trackId:
			ELoss += hit.GetEnergyLoss()
	if ELoss >= ELossThreshold:
		return True
	return False

def hasMuons(event, trackId):
    m1 = 0
    m2 = 0
    m3 = 0
    m4 = 0
    for ahit in event.muonPoint:
    	if ahit.GetTrackID() == trackId:
        	detID = ahit.GetDetectorID()
        	if(detID == 476) :
        	    m1 += 1
        	if(detID == 477) :
        	    m2 += 1
        	if(detID == 478) :
        	    m3 += 1
        	if(detID == 479) :
        	    m4 += 1
    return [bool(m1), bool(m2), bool(m3), bool(m4)]

def myVertex(t1,t2,PosDir):
    # closest distance between two tracks
    # d = |pq . u x v|/|u x v|
    a = ROOT.TVector3(PosDir[t1][0](0) ,PosDir[t1][0](1), PosDir[t1][0](2))
    u = ROOT.TVector3(PosDir[t1][1](0),PosDir[t1][1](1),PosDir[t1][1](2))
    c = ROOT.TVector3(PosDir[t2][0](0) ,PosDir[t2][0](1), PosDir[t2][0](2))
    v = ROOT.TVector3(PosDir[t2][1](0),PosDir[t2][1](1),PosDir[t2][1](2))
    pq = a-c
    uCrossv = u.Cross(v)
    dist  = pq.Dot(uCrossv)/(uCrossv.Mag()+1E-8)
    # u.a - u.c + s*|u|**2 - u.v*t    = 0
    # v.a - v.c + s*v.u    - t*|v|**2 = 0
    E = u.Dot(a) - u.Dot(c) 
    F = v.Dot(a) - v.Dot(c) 
    A,B = u.Mag2(), -u.Dot(v) 
    C,D = u.Dot(v), -v.Mag2()
    t = -(C*E-A*F)/(B*C-A*D)
    X = c.x()+v.x()*t
    Y = c.y()+v.y()*t
    Z = c.z()+v.z()*t
    # sT = ROOT.gROOT.FindAnything('cbmsim')
    #print 'test2 ',X,Y,Z,dist
    #print 'truth',sTree.MCTrack[2].GetStartX(),sTree.MCTrack[2].GetStartY(),sTree.MCTrack[2].GetStartZ()
    return X,Y,Z,abs(dist)

def addFullInfoToTree(elenaTree):
	elenaTree, DaughtersPt = tools.AddVect(elenaTree, 'DaughtersPt', 'float')
	elenaTree, DaughtersChi2 = tools.AddVect(elenaTree, 'DaughtersChi2', 'float')
	elenaTree, DaughtersNPoints = tools.AddVect(elenaTree, 'DaughtersNPoints', 'int')
	elenaTree, DaughtersTruthProdX = tools.AddVect(elenaTree, 'DaughtersTruthProdX', 'float')
	elenaTree, DaughtersTruthProdY = tools.AddVect(elenaTree, 'DaughtersTruthProdY', 'float')
	elenaTree, DaughtersTruthProdZ = tools.AddVect(elenaTree, 'DaughtersTruthProdZ', 'float')
	elenaTree, DaughtersTruthPDG = tools.AddVect(elenaTree, 'DaughtersTruthPDG', 'int')
	elenaTree, DaughtersTruthMotherPDG = tools.AddVect(elenaTree, 'DaughtersTruthMotherPDG', 'int')
	elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int')
	elenaTree, straw_x = tools.AddVect(elenaTree, 'straw_x', 'float')
	elenaTree, straw_y = tools.AddVect(elenaTree, 'straw_y', 'float')
	elenaTree, straw_z = tools.AddVect(elenaTree, 'straw_z', 'float')
	elenaTree, muon_x = tools.AddVect(elenaTree, 'muon_x', 'float')
	elenaTree, muon_y = tools.AddVect(elenaTree, 'muon_y', 'float')
	elenaTree, muon_z = tools.AddVect(elenaTree, 'muon_z', 'float')
	elenaTree, ecal_x = tools.AddVect(elenaTree, 'ecal_x', 'float')
	elenaTree, ecal_y = tools.AddVect(elenaTree, 'ecal_y', 'float')
	elenaTree, ecal_z = tools.AddVect(elenaTree, 'ecal_z', 'float')
	elenaTree, hcal_x = tools.AddVect(elenaTree, 'hcal_x', 'float')
	elenaTree, hcal_y = tools.AddVect(elenaTree, 'hcal_y', 'float')
	elenaTree, hcal_z = tools.AddVect(elenaTree, 'hcal_z', 'float')
	elenaTree, veto5_x = tools.AddVect(elenaTree, 'veto5_x', 'float')
	elenaTree, veto5_y = tools.AddVect(elenaTree, 'veto5_y', 'float')
	elenaTree, veto5_z = tools.AddVect(elenaTree, 'veto5_z', 'float')
	elenaTree, liquidscint_x = tools.AddVect(elenaTree, 'liquidscint_x', 'float')
	elenaTree, liquidscint_y = tools.AddVect(elenaTree, 'liquidscint_y', 'float')
	elenaTree, liquidscint_z = tools.AddVect(elenaTree, 'liquidscint_z', 'float')
	elenaTree, DOCA = tools.AddVar(elenaTree, 'DOCA', 'float')
	elenaTree, vtxx = tools.AddVar(elenaTree, 'vtxx', 'float')
	elenaTree, vtxy = tools.AddVar(elenaTree, 'vtxy', 'float')
	elenaTree, vtxz = tools.AddVar(elenaTree, 'vtxz', 'float')
	elenaTree, IP0 = tools.AddVar(elenaTree, 'IP0', 'float')
	elenaTree, Mass = tools.AddVar(elenaTree, 'Mass', 'float')
	elenaTree, Pt = tools.AddVar(elenaTree, 'Pt', 'float')
	elenaTree, P = tools.AddVar(elenaTree, 'P', 'float')
	elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int')
	elenaTree, HNLw = tools.AddVar(elenaTree, 'HNLw', 'float')
	elenaTree, NuWeight = tools.AddVar(elenaTree, 'NuWeight', 'float')
	elenaTree, EventNumber = tools.AddVar(elenaTree, 'EventNumber', 'int')
	elenaTree, DaughterMinPt = tools.AddVar(elenaTree, 'DaughterMinPt', 'float')
	elenaTree, DaughterMinP = tools.AddVar(elenaTree, 'DaughterMinP', 'float')
	elenaTree, DaughtersAlwaysIn = tools.AddVar(elenaTree, 'DaughtersAlwaysIn', 'int')
	elenaTree, BadTruthVtx = tools.AddVar(elenaTree, 'BadTruthVtx', 'int')

DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal = None, None, None, None, None, None, None
DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 = None, None, None, None, None, None
MaxDaughtersRedChi2, MinDaughtersNdf = None, None

def addOfflineToTree(elenaTree):
    global DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal
    global DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2
    global MaxDaughtersRedChi2, MinDaughtersNdf, HNLw, NuWeight
    elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') #
    elenaTree, DOCA = tools.AddVect(elenaTree, 'DOCA', 'float') #
    elenaTree, vtxx = tools.AddVect(elenaTree, 'vtxx', 'float') #
    elenaTree, vtxy = tools.AddVect(elenaTree, 'vtxy', 'float') #
    elenaTree, vtxz = tools.AddVect(elenaTree, 'vtxz', 'float') #
    elenaTree, IP0 = tools.AddVect(elenaTree, 'IP0', 'float') #
    #elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int') #
    elenaTree, DaughtersAlwaysIn = tools.AddVect(elenaTree, 'DaughtersAlwaysIn', 'int') #
    elenaTree, BadTruthVtx = tools.AddVect(elenaTree, 'BadTruthVtx', 'int') #
    elenaTree, Has1Muon1 = tools.AddVect(elenaTree, 'Has1Muon1', 'int') #
    elenaTree, Has1Muon2 = tools.AddVect(elenaTree, 'Has1Muon2', 'int') #
    elenaTree, Has2Muon1 = tools.AddVect(elenaTree, 'Has2Muon1', 'int') #
    elenaTree, Has2Muon2 = tools.AddVect(elenaTree, 'Has2Muon2', 'int') #
    elenaTree, HasEcal = tools.AddVect(elenaTree, 'HasEcal', 'int') #
    elenaTree, MaxDaughtersRedChi2 = tools.AddVect(elenaTree, 'MaxDaughtersRedChi2', 'float') #
    elenaTree, MinDaughtersNdf = tools.AddVect(elenaTree, 'MinDaughtersNdf', 'int') #
    # Add liquid scintillator segmentation
    tools.makeLSsegments(nodes, elenaTree)

nodes = None
def loadNodes(fGeo):
	global nodes
	nodes = searchForNodes3_xyz_dict(fGeo)

num_bad_z = 0

def signalNormalisationZ(tree, datatype, verbose):
    # By event
    # Uses MC truth!!
    global BadTruthVtx, num_bad_z
    tools.PutToZero(BadTruthVtx)
    z_hnl_vtx = findHNLvertex(tree)
    bad_z = False
    if not z_hnl_vtx:
    	if "sig" in datatype:
    		print 'ERROR: hnl vertex not found!'
    	ii = 0
    	for g in tree.MCTrack:
    		ii +=1
    	if ("sig" in datatype) and ii < 3:
    		pass
    	elif ("sig" in datatype) and ii >= 3:
    		sys.exit()
    if not (nodes['Veto_5']['z']['pos']-nodes['Veto_5']['z']['dim']-500. < z_hnl_vtx < nodes['Tr4_4']['z']['pos']+nodes['Tr4_4']['z']['dim']):
    	bad_z = True
    	num_bad_z += 1
    	if "sig" in datatype:
    		print z_hnl_vtx
    tools.Push(BadTruthVtx, int(bad_z))

def nParticles(tree):
	# By event
    global NParticles
    np = 0
    for HNL in tree.Particles:
        np += 1
    tools.PutToZero(NParticles); tools.Push(NParticles, np)

def hasEcalAndMuons(tree, particle):
	# By particle
	global Has1Muon1, Has1Muon2, Has2Muon1
	global Has2Muon2, HasEcal
	flag2Muon1 = False
	flag2Muon2 = False
	flag1Muon1 = False
	flag1Muon2 = False
	flagEcal = False
	t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) 
	# AND or OR?
	if ( has_muon_station(tree, t1, 1) and has_muon_station(tree, t2, 1) ):
		flag2Muon1 = True
	if ( has_muon_station(tree, t1, 2) and has_muon_station(tree, t2, 2) ):
		flag2Muon2 = True
	if ( has_muon_station(tree, t1, 1) or has_muon_station(tree, t2, 1) ):
		flag1Muon1 = True
	if ( has_muon_station(tree, t1, 2) or has_muon_station(tree, t2, 2) ):
		flag1Muon2 = True
	if ( hasEcalDeposit(tree, t1, 150.*u.MeV) or hasEcalDeposit(tree, t2, 150.*u.MeV)	):
		flagEcal = True
	tools.Push(Has2Muon1, int(flag2Muon1))
	tools.Push(Has2Muon2, int(flag2Muon2))
	tools.Push(Has1Muon1, int(flag1Muon1))
	tools.Push(Has1Muon2, int(flag1Muon2))
	tools.Push(HasEcal, int(flagEcal))

def chi2Ndf(tree, particle):
	# By particle
	global MaxDaughtersRedChi2, MinDaughtersNdf, DaughtersFitConverged
	t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1)
	reducedChi2 = []
	ndfs = []
	converged = []
	for tr in [t1,t2]:
		x = tree.FitTracks[tr]
		ndfs.append( int(round(x.getFitStatus().getNdf())) )
		reducedChi2.append( x.getFitStatus().getChi2()/x.getFitStatus().getNdf() )
		converged.append( x.getFitStatus().isFitConverged() )
	tools.Push(MaxDaughtersRedChi2, max(reducedChi2))
	tools.Push(MinDaughtersNdf, min(ndfs))
	tools.Push( DaughtersFitConverged, int(converged[0]*converged[1]) )

	
def goodBehavedTracks(tree, particle):
	# By particle
	# Uses MC truth!!
	global DaughtersAlwaysIn
	t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) 
	accFlag = True
	for tr in [t1,t2]:
		mctrid = tree.fitTrack2MC[tr]
		if not hasGoodStrawStations(tree, mctrid):
			accFlag = False
	tools.Push(DaughtersAlwaysIn, int(accFlag))

def vertexInfo(tree, particle):
	# By particle
	global vtxx, vtxy, vtxz
	global IP0, DOCA
	t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) 
	PosDir = {} 
	for tr in [t1,t2]:
		xx  = tree.FitTracks[tr].getFittedState()
		PosDir[tr] = [xx.getPos(),xx.getDir()]
	xv,yv,zv,doca = myVertex(t1,t2,PosDir)
	tools.Push(DOCA, doca)
	tools.Push(vtxx, xv); tools.Push(vtxy, yv); tools.Push(vtxz, zv)
	# impact parameter to target
	HNLPos = ROOT.TLorentzVector()
	particle.ProductionVertex(HNLPos)
	HNLMom = ROOT.TLorentzVector()
	particle.Momentum(HNLMom)
	tr = ROOT.TVector3(0,0,ShipGeo.target.z0)
	t = 0
	for i in range(3):   t += HNLMom(i)/HNLMom.P()*(tr(i)-HNLPos(i)) 
	ip = 0
	for i in range(3):   ip += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2
	ip = ROOT.TMath.Sqrt(ip)
	tools.Push(IP0, ip)


def prepareFillingsByParticle():
	# By event
	global DaughtersAlwaysIn, DaughtersFitConverged, MinDaughtersNdf, MaxDaughtersRedChi2
	global Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2, HasEcal
	global vtxx, vtxy, vtxz, IP0, DOCA
	tools.PutToZero(DaughtersAlwaysIn)
	tools.PutToZero(Has2Muon1); tools.PutToZero(Has2Muon2); tools.PutToZero(HasEcal)
	tools.PutToZero(Has1Muon1); tools.PutToZero(Has1Muon2)
	tools.PutToZero(DOCA)
	tools.PutToZero(vtxx); tools.PutToZero(vtxy); tools.PutToZero(vtxz)
	tools.PutToZero(IP0)
	tools.PutToZero(MinDaughtersNdf); tools.PutToZero(MaxDaughtersRedChi2)
	tools.PutToZero(DaughtersFitConverged)


def pushOfflineByEvent(tree, datatype, verbose):
    # True HNL decay vertex (only for signal normalisation)
    signalNormalisationZ(tree, datatype, verbose)
    ## Number of particles
    #nParticles(tree)
    # Empties arrays filled by particle
    prepareFillingsByParticle()
    # Liquid scintillator segments
    global nodes
    tools.hitSegments(tree, nodes)

def pushOfflineByParticle(tree, particle):
	hasEcalAndMuons(tree, particle)
	goodBehavedTracks(tree, particle)
	chi2Ndf(tree, particle)
	vertexInfo(tree, particle)