diff --git a/ShipStyle.py b/ShipStyle.py
new file mode 100644
index 0000000..0bb389e
--- /dev/null
+++ b/ShipStyle.py
@@ -0,0 +1,189 @@
+# Global scope variables
+from ROOT import *
+from ROOT import Double
+lhcbstyle = TStyle()     # general lhcb style
+lhcbName  = TPaveText()  # standard lhcb text for plot
+lhcbLabel = TText  # style for Ttext
+lhcbLatex = TLatex #style for TLatex
+
+# define names for colours
+black=1
+red=2
+green=3
+blue=4
+yellow=5 
+magenta=6
+cyan=7
+purple=9
+
+def lhcbstyleSetup():
+
+##################################
+# PURPOSE:
+#
+# This macro defines a standard style for (black-and-white) 
+# "publication quality" LHCb ROOT plots. 
+#
+# USAGE:
+#
+# Include the lines
+#   gROOT.ProcessLine(".L lhcbstyle.C")
+#   lhcbstyle()
+# at the beginning of your root macro.
+#
+# Example usage is given in myPlot.C
+#
+# COMMENTS:
+#
+# Font:
+# 
+# The font is chosen to be 62, this is helvetica-bold-r-normal with
+# precision 2.
+#
+# "Landscape histograms":
+#
+# The style here is designed for more or less square plots.
+# For longer histograms, or canvas with many pads, adjustements are needed. 
+# For instance, for a canvas with 1x5 histograms:
+#  TCanvas* c1 = new TCanvas("c1", "L0 muons", 600, 800)
+#  c1.Divide(1,5)
+#  Adaptions like the following will be needed:
+#  lhcbstyle.SetTickLength(0.05,"x")
+#  lhcbstyle.SetTickLength(0.01,"y")
+#  lhcbstyle.SetLabelSize(0.15,"x")
+#  lhcbstyle.SetLabelSize(0.1,"y")
+#  lhcbstyle.SetStatW(0.15)
+#  lhcbstyle.SetStatH(0.5)
+#
+# Authors: Thomas Schietinger, Andrew Powell, Chris Parkes
+# Maintained by Editorial board member (currently Chris)
+#################################/
+
+ lhcbstyle=TStyle("lhcbstyle","Standard LHCb plots style")
+
+# use helvetica-bold-r-normal, precision 2 (rotatable)
+ lhcbFont = 62
+# line thickness
+ lhcbWidth = int(3.00)
+
+# use plain black on white colors
+ lhcbstyle.SetFrameBorderMode(0)
+ lhcbstyle.SetCanvasBorderMode(0)
+ lhcbstyle.SetPadBorderMode(0)
+ lhcbstyle.SetPadColor(0)
+ lhcbstyle.SetCanvasColor(0)
+ lhcbstyle.SetStatColor(0)
+ lhcbstyle.SetPalette(1)
+
+# set the paper & margin sizes
+ lhcbstyle.SetPaperSize(20,26)
+ lhcbstyle.SetPadTopMargin(0.05)
+ lhcbstyle.SetPadRightMargin(0.05) # increase for colz plots
+ lhcbstyle.SetPadBottomMargin(0.16)
+ lhcbstyle.SetPadLeftMargin(0.14)
+
+# use large fonts
+ lhcbstyle.SetTextFont(lhcbFont)
+ lhcbstyle.SetTextSize(0.08)
+ lhcbstyle.SetLabelFont(lhcbFont,"x")
+ lhcbstyle.SetLabelFont(lhcbFont,"y")
+ lhcbstyle.SetLabelFont(lhcbFont,"z")
+ lhcbstyle.SetLabelSize(0.05,"x")
+ lhcbstyle.SetLabelSize(0.05,"y")
+ lhcbstyle.SetLabelSize(0.05,"z")
+ lhcbstyle.SetTitleFont(lhcbFont)
+ lhcbstyle.SetTitleSize(0.06,"x")
+ lhcbstyle.SetTitleSize(0.06,"y")
+ lhcbstyle.SetTitleSize(0.06,"z")
+
+# use bold lines and markers
+ lhcbstyle.SetLineWidth(lhcbWidth)
+ lhcbstyle.SetFrameLineWidth(lhcbWidth)
+ lhcbstyle.SetHistLineWidth(lhcbWidth)
+ lhcbstyle.SetFuncWidth(lhcbWidth)
+ lhcbstyle.SetGridWidth(lhcbWidth)
+ lhcbstyle.SetLineStyleString(2,"[12 12]") # postscript dashes
+ lhcbstyle.SetMarkerStyle(20)
+ lhcbstyle.SetMarkerSize(1.5)
+
+# label offsets
+ lhcbstyle.SetLabelOffset(0.015)
+
+# by default, do not display histogram decorations:
+ lhcbstyle.SetOptStat(0)  
+ lhcbstyle.SetOptStat("emr")  # show only nent -e , mean - m , rms -r
+# full opts at http:#root.cern.ch/root/html/TStyle.html#TStyle:SetOptStat
+ lhcbstyle.SetStatFormat("6.3g") # specified as c printf options
+ lhcbstyle.SetOptTitle(0)
+ lhcbstyle.SetOptFit(0)
+#lhcbstyle.SetOptFit(1011) # order is probability, Chi2, errors, parameters
+
+# look of the statistics box:
+ lhcbstyle.SetStatBorderSize(0)
+ lhcbstyle.SetStatFont(lhcbFont)
+ lhcbstyle.SetStatFontSize(0.05)
+ lhcbstyle.SetStatX(0.9)
+ lhcbstyle.SetStatY(0.9)
+ lhcbstyle.SetStatW(0.25)
+ lhcbstyle.SetStatH(0.15)
+# put tick marks on top and RHS of plots
+ lhcbstyle.SetPadTickX(1)
+ lhcbstyle.SetPadTickY(1)
+
+# histogram divisions: only 5 in x to avoid label overlaps
+ lhcbstyle.SetNdivisions(505,"x")
+ lhcbstyle.SetNdivisions(510,"y")
+
+
+#define style for text
+ lhcbLabel = TText()
+ lhcbLabel.SetTextFont(lhcbFont)
+ lhcbLabel.SetTextColor(1)
+ lhcbLabel.SetTextSize(0.04)
+ lhcbLabel.SetTextAlign(12)
+
+# define style of latex text
+ lhcbLatex = TLatex()
+ lhcbLatex.SetTextFont(lhcbFont)
+ lhcbLatex.SetTextColor(1)
+ lhcbLatex.SetTextSize(0.04)
+ lhcbLatex.SetTextAlign(12)
+
+# set this style
+ gROOT.SetStyle("lhcbstyle")
+ gROOT.ForceStyle()
+
+
+def printLHCb(optLR="L", optPrelim="Final", optText=""):
+#####################################
+# routine to print 'LHCb', 'LHCb Preliminary' on plots 
+# options: optLR=L (top left) / R (top right) of plots
+#          optPrelim= Final (LHCb), Prelim (LHCb Preliminary), Other
+#          optText= text printed if 'Other' specified
+##################################
+  if optLR=="R" :   
+    lhcbName = TPaveText(0.70 - lhcbstyle.GetPadRightMargin(),
+					   0.75 - lhcbstyle.SetPadTopMargin(0.05),
+					   0.95 - lhcbstyle.GetPadRightMargin(),
+					   0.85 - lhcbstyle.SetPadTopMargin(0.05),
+					   "BRNDC")
+  elif optLR=="L":
+    lhcbName = TPaveText(lhcbstyle.GetPadLeftMargin() + 0.05,
+				       0.87 - lhcbstyle.GetPadTopMargin(),
+				       lhcbstyle.GetPadLeftMargin() + 0.30,
+				       0.95 - lhcbstyle.GetPadTopMargin(),
+                                        "BRNDC")
+  else :
+   print "printLHCb: option unknown" , optLR  
+  if optPrelim=="Final":
+    lhcbName.AddText("LHCb")
+  elif optPrelim=="Prelim":
+    lhcbName.AddText("#splitline{LHCb}{#scale[1.0]{Preliminary}}")  
+  elif optPrelim=="Other":
+    lhcbName.AddText(optText)
+  else :
+    print "printLHCb: option unknown " , optPrelim
+  lhcbName.SetFillColor(0)
+  lhcbName.SetTextAlign(12)
+  lhcbName.SetBorderSize(0)
+  lhcbName.Draw()
diff --git a/newGen.tar.gz b/newGen.tar.gz
new file mode 100644
index 0000000..6286cc5
--- /dev/null
+++ b/newGen.tar.gz
Binary files differ
diff --git a/newGen/offlineForBarbara.py b/newGen/offlineForBarbara.py
index 954d478..0731538 100644
--- a/newGen/offlineForBarbara.py
+++ b/newGen/offlineForBarbara.py
@@ -1,669 +1,669 @@
-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
+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
diff --git a/newGen/offlineForBarbara_maybeOld.py b/newGen/offlineForBarbara_maybeOld.py
new file mode 100644
index 0000000..954d478
--- /dev/null
+++ b/newGen/offlineForBarbara_maybeOld.py
@@ -0,0 +1,669 @@
+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