from lookAtGeo import *
import tools
import shipunit as u
from ShipGeoConfig import ConfigRegistry
import shipDet_conf

from operator import mul, add

import sys
sys.path.append('KaterinaLight/')
from StrawHits import StrawHits
## Use it like:
# f = TFile(fileName)
# t = f.Get("cbmsim")
# sh = offline.StrawHits(t, offline.shipDet_conf.configure(offline.__run, r['ShipGeo']), r['ShipGeo'].straw.resol, 0, None, r['ShipGeo'])
# t.GetEntry(58)
# sh.readEvent()
# sh.FitTracks()

#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_nu = ROOT.TFile("histoForWeights_nu.root")
h_GioHans_nu = ff_nu.Get("h_Gio")

ff_antinu = ROOT.TFile("histoForWeights_antinu.root")
h_GioHans_antinu = ff_antinu.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 = 1.#6.98e+11 * 2.e+20 / 5.e+13
        h_GioHans = h_GioHans_antinu
    else: 
        reduction = 1.
        flux = 1.#1.09e+12 * 2.e+20/ 5.e+13
        h_GioHans = h_GioHans_nu
        
    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
NoB_DOCA, NoB_vtxx, NoB_vtxy, NoB_vtxz, NoB_IP0 = None, None, None, None, None
DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 = None, None, None, None, None, None
MaxDaughtersRedChi2, MinDaughtersNdf = None, None
NoB_MaxDaughtersRedChi2, NoB_MinDaughtersNdf = None, None
DaughtersMinP, DaughtersMinPt, Mass, P, Pt = None, None, None, None, None
NoB_DaughtersMinP, NoB_DaughtersMinPt, NoB_Mass, NoB_P, NoB_Pt = None, None, None, None, None

def addOfflineToTree(elenaTree):
    global DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal
    global NoB_DOCA, NoB_vtxx, NoB_vtxy, NoB_vtxz, NoB_IP0
    global DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2
    global MaxDaughtersRedChi2, MinDaughtersNdf, HNLw, NuWeight, NoB_MaxDaughtersRedChi2, NoB_MinDaughtersNdf
    global DaughtersMinP, DaughtersMinPt, Mass, P, Pt
    global NoB_DaughtersMinP, NoB_DaughtersMinPt, NoB_Mass, NoB_P, NoB_Pt
    elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') #
    elenaTree, DOCA = tools.AddVect(elenaTree, 'DOCA', 'float') #
    elenaTree, NoB_DOCA = tools.AddVect(elenaTree, 'NoB_DOCA', 'float') #
    elenaTree, vtxx = tools.AddVect(elenaTree, 'vtxxSqr', 'float') #
    elenaTree, vtxy = tools.AddVect(elenaTree, 'vtxySqr', 'float') #
    elenaTree, vtxz = tools.AddVect(elenaTree, 'vtxz', 'float') #
    elenaTree, NoB_vtxx = tools.AddVect(elenaTree, 'NoB_vtxxSqr', 'float') #
    elenaTree, NoB_vtxy = tools.AddVect(elenaTree, 'NoB_vtxySqr', 'float') #
    elenaTree, NoB_vtxz = tools.AddVect(elenaTree, 'NoB_vtxz', 'float') #
    elenaTree, IP0 = tools.AddVect(elenaTree, 'IP0', 'float') #
    elenaTree, NoB_IP0 = tools.AddVect(elenaTree, 'NoB_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') #
    elenaTree, NoB_MaxDaughtersRedChi2 = tools.AddVect(elenaTree, 'NoB_MaxDaughtersRedChi2', 'float') #
    elenaTree, NoB_MinDaughtersNdf = tools.AddVect(elenaTree, 'NoB_MinDaughtersNdf', 'int') #
    elenaTree, DaughtersMinP = tools.AddVect(elenaTree, 'DaughtersMinP', 'float')
    elenaTree, DaughtersMinPt = tools.AddVect(elenaTree, 'DaughtersMinPt', 'float')
    elenaTree, P = tools.AddVect(elenaTree, 'P', 'float')
    elenaTree, Pt = tools.AddVect(elenaTree, 'Pt', 'float')
    elenaTree, Mass = tools.AddVect(elenaTree, 'Mass', 'float')
    elenaTree, NoB_DaughtersMinP = tools.AddVect(elenaTree, 'NoB_DaughtersMinP', 'float')
    elenaTree, NoB_DaughtersMinPt = tools.AddVect(elenaTree, 'NoB_DaughtersMinPt', 'float')
    elenaTree, NoB_P = tools.AddVect(elenaTree, 'NoB_P', 'float')
    elenaTree, NoB_Pt = tools.AddVect(elenaTree, 'NoB_Pt', 'float')
    elenaTree, NoB_Mass = tools.AddVect(elenaTree, 'NoB_Mass', 'float')
    # 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 = tree.fitTrack2MC[particle.GetDaughter(0)], tree.fitTrack2MC[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
    # This also work, but may be slower
    #muons1 = hasMuons(tree, t1)
    #muons2 = hasMuons(tree, t2)
    #if muons1[0] or muons2[0]: flag1Muon1 = True
    #if muons1[1] or muons2[1]: flag1Muon2 = True
    #if muons1[0] and muons2[0]: flag2Muon1 = True
    #if muons1[1] and muons2[1]: flag2Muon2 = 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, ntr, nref):
    # By particle
    global MaxDaughtersRedChi2, MinDaughtersNdf
    t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1)
    if ntr>1 and nref==2:#nf>1
        t1r,t2r = sh.getReFitTrIDs()[0], sh.getReFitTrIDs()[1]
        chi2red_1 = sh.getReFitChi2Ndf(t1r)
        ndf_1 = int(round(sh.getReFitNdf(t1r)))
        chi2red_2 = sh.getReFitChi2Ndf(t2r)
        ndf_2 = int(round(sh.getReFitNdf(t2r)))
        reducedChi2 = [chi2red_1, chi2red_2]
        ndfs = [ndf_1, ndf_2]
    # if the refit didn't work
    if (ntr<2) or (nref!=2) or (not ndf_1) or (not ndf_2) or (not chi2red_1) or (not chi2red_2):
        reducedChi2 = []
        ndfs = []
        for tr in [t1,t2]:
            x = tree.FitTracks[tr]
            ndfs.append( int(round(x.getFitStatus().getNdf())) )
            reducedChi2.append( x.getFitStatus().getChi2()/x.getFitStatus().getNdf() )
    tools.Push(MaxDaughtersRedChi2, max(reducedChi2))
    tools.Push(MinDaughtersNdf, min(ndfs))

    
def NoB_chi2Ndf(tree, particle):
    # By particle
    global NoB_MaxDaughtersRedChi2, NoB_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(NoB_MaxDaughtersRedChi2, max(reducedChi2))
    tools.Push(NoB_MinDaughtersNdf, min(ndfs))
    tools.Push( DaughtersFitConverged, int(converged[0]*converged[1]) )

def NoB_kinematics(tree, particle):
    global NoB_DaughtersMinP, NoB_DaughtersMinPt, NoB_P, NoB_Pt, NoB_Mass
    t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1)
    dp, dpt = [], []
    for tr in [t1, t2]:
        x = tree.FitTracks[tr]
        xx  = x.getFittedState()
        dp.append(xx.getMom().Mag()); dpt.append(xx.getMom().Pt())
    tools.Push(NoB_DaughtersMinP, min(dp))
    tools.Push(NoB_DaughtersMinPt, min(dpt))
    HNLMom = ROOT.TLorentzVector()
    particle.Momentum(HNLMom)
    tools.Push(NoB_Mass, HNLMom.M())
    tools.Push(NoB_Pt, HNLMom.Pt())
    tools.Push(NoB_P, HNLMom.P())
    
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 NoB_vertexInfo(tree, particle):
    # By particle
    global NoB_vtxx, NoB_vtxy, NoB_vtxz
    global NoB_IP0, NoB_DOCA
    HNLPos = ROOT.TLorentzVector()
    particle.ProductionVertex(HNLPos)
    xv, yv, zv, doca = HNLPos.X(),HNLPos.Y(),HNLPos.Z(),HNLPos.T()
    tools.Push(NoB_DOCA, doca)
    tools.Push(NoB_vtxx, xv*xv); tools.Push(NoB_vtxy, yv*yv); tools.Push(NoB_vtxz, zv)
    # impact parameter to target
    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(NoB_IP0, ip)
    """
    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(NoB_DOCA, doca)
    #tools.Push(NoB_vtxx, xv); tools.Push(NoB_vtxy, yv); tools.Push(NoB_vtxz, zv)
    tools.Push(NoB_vtxx, xv*xv); tools.Push(NoB_vtxy, yv*yv); tools.Push(NoB_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(NoB_IP0, ip)
    """


def kinematics(tree, particle, ntr, nref):
    global DaughtersMinP, DaughtersMinPt, P, Pt, Mass
    t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1)
    dminpt, dminp = 0., 0.

    if ntr>1 and nref==2:
        t1r,t2r = sh.getReFitTrIDs()[0], sh.getReFitTrIDs()[1]
        Pos1, Dir1, Mom1= sh.getReFitPosDirPval(t1r)
        Pos2, Dir2, Mom2= sh.getReFitPosDirPval(t2r)
        mass1 = pdg.GetParticle(tree.FitTracks[t1].getFittedState().getPDG()).Mass()
        mass2 = pdg.GetParticle(tree.FitTracks[t2].getFittedState().getPDG()).Mass()
        LV1 = ROOT.TLorentzVector(Mom1*Dir1, ROOT.TMath.Sqrt( mass1*mass1 + Mom1*Mom1 ))
        LV2 = ROOT.TLorentzVector(Mom2*Dir2, ROOT.TMath.Sqrt( mass2*mass2 + Mom2*Mom2 ))
        HNLMom = LV1+LV2
        if LV1 and LV2:
            dminpt = min([LV1.Pt(), LV2.Pt()])
            dminp = min([LV1.P(), LV2.P()])

    if (ntr<2) or (nref!=2) or (not dminp) or (not dminpt) or (not HNLMom):
        dp, dpt = [], []
        for tr in [t1, t2]:
            x = tree.FitTracks[tr]
            xx  = x.getFittedState()
            dp.append(xx.getMom().Mag()); dpt.append(xx.getMom().Pt())
        dminpt = min(dpt)
        dminp = min(dp)
        HNLMom = ROOT.TLorentzVector()
        particle.Momentum(HNLMom)
    tools.Push(DaughtersMinP, dminp)
    tools.Push(DaughtersMinPt, dminpt)
    tools.Push(Mass, HNLMom.M())
    tools.Push(Pt, HNLMom.Pt())
    tools.Push(P, HNLMom.P())
    
def vertexInfo(tree, particle, ntr, nref):
    # By particle
    global vtxx, vtxy, vtxz
    global IP0, DOCA
    t1,t2 = particle.GetDaughter(0),particle.GetDaughter(1) 

    if ntr>1 and nref==2:#nf>1
        assert( len(sh.getReFitTrIDs())==2 )
        t1r,t2r = sh.getReFitTrIDs()[0], sh.getReFitTrIDs()[1]
        #print tree.fitTrack2MC[t1], t1r, tree.fitTrack2MC[t2], t2r
        #print ntr, nref, len(sh._StrawHits__docaEval)
        doca = sh.getDoca()#sh._StrawHits__docaEval[-1]#getDoca()
        v = sh.getReFitVertex()
        if v and doca:
            xv = v.X(); yv = v.Y(); zv = v.Z()
            Pos1, Dir1, Mom1= sh.getReFitPosDirPval(t1r)
            Pos2, Dir2, Mom2= sh.getReFitPosDirPval(t2r)
            mass1 = pdg.GetParticle(tree.FitTracks[t1].getFittedState().getPDG()).Mass()
            mass2 = pdg.GetParticle(tree.FitTracks[t2].getFittedState().getPDG()).Mass()
            LV1 = ROOT.TLorentzVector(Mom1*Dir1, ROOT.TMath.Sqrt( mass1*mass1 + Mom1*Mom1 ))
            LV2 = ROOT.TLorentzVector(Mom2*Dir2, ROOT.TMath.Sqrt( mass2*mass2 + Mom2*Mom2 ))
            HNLMom = LV1+LV2

    # If something went wrong, take the standard values
    if (ntr<2) or (nref!=2) or (not v) or (not doca) or (not HNLMom):#(nf<2)
        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)
        HNLMom = ROOT.TLorentzVector()
        particle.Momentum(HNLMom)

    tools.Push(DOCA, doca)
    #tools.Push(vtxx, xv); tools.Push(vtxy, yv); tools.Push(vtxz, zv)
    tools.Push(vtxx, xv*xv); tools.Push(vtxy, yv*yv); tools.Push(vtxz, zv)

    # impact parameter to target
    #HNLPos = ROOT.TLorentzVector()
    #particle.ProductionVertex(HNLPos)
    HNLPos = ROOT.TVector3(xv, yv, zv)
    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 NoB_MinDaughtersNdf, NoB_MaxDaughtersRedChi2
    global Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2, HasEcal
    global vtxx, vtxy, vtxz, IP0, DOCA
    global NoB_vtxx, NoB_vtxy, NoB_vtxz, NoB_IP0, NoB_DOCA
    global DaughtersMinP, DaughtersMinPt, Mass, P, Pt
    global NoB_DaughtersMinP, NoB_DaughtersMinPt, NoB_Mass, NoB_P, NoB_Pt
    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(NoB_DOCA)
    tools.PutToZero(NoB_vtxx); tools.PutToZero(NoB_vtxy); tools.PutToZero(NoB_vtxz)
    tools.PutToZero(NoB_IP0)
    tools.PutToZero(MinDaughtersNdf); tools.PutToZero(MaxDaughtersRedChi2)
    tools.PutToZero(NoB_MinDaughtersNdf); tools.PutToZero(NoB_MaxDaughtersRedChi2)
    tools.PutToZero(DaughtersFitConverged)
    tools.PutToZero(DaughtersMinP); tools.PutToZero(DaughtersMinPt)
    tools.PutToZero(P); tools.PutToZero(Pt); tools.PutToZero(Mass)
    tools.PutToZero(NoB_DaughtersMinP); tools.PutToZero(NoB_DaughtersMinPt)
    tools.PutToZero(NoB_P); tools.PutToZero(NoB_Pt); tools.PutToZero(NoB_Mass)
    ntr = sh.readEvent()
    nref = 0
    if ntr>1:
        nref = sh.FitTracks()
        #print ntr, nref
    return ntr, nref 


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

def pushOfflineByParticle(tree, particle, ntr, nref):
    hasEcalAndMuons(tree, particle)
    goodBehavedTracks(tree, particle)
    NoB_chi2Ndf(tree, particle)
    chi2Ndf(tree, particle, ntr, nref)
    NoB_vertexInfo(tree, particle)
    vertexInfo(tree, particle, ntr, nref)
    NoB_kinematics(tree, particle)
    kinematics(tree, particle, ntr, nref)

fM, tgeom, gMan, geoMat, matEff, modules, run = None, None, None, None, None, None, None

def initBField(fileNameGeo):
    global fM, tgeom, gMan, geoMat, matEff, modules, run, sh
    run     = ROOT.FairRunSim()
    modules = shipDet_conf.configure(run,ShipGeo)
    tgeom = ROOT.TGeoManager("Geometry", "Geane geometry")
    gMan  = tgeom.Import(fileNameGeo)
    geoMat =  ROOT.genfit.TGeoMaterialInterface()
    matEff = ROOT.genfit.MaterialEffects.getInstance()
    matEff.init(geoMat)
    bfield = ROOT.genfit.BellField(ShipGeo.Bfield.max, ShipGeo.Bfield.z, 2, ShipGeo.Yheight/2.)
    fM = ROOT.genfit.FieldManager.getInstance()
    fM.init(bfield)

pdg, sh = None, None