diff --git a/newGen/ntuple_nuVeto_small_newVar.py b/newGen/ntuple_nuVeto_small_newVar.py
new file mode 100644
index 0000000..f443c3a
--- /dev/null
+++ b/newGen/ntuple_nuVeto_small_newVar.py
@@ -0,0 +1,717 @@
+import sys, re, getopt
+from array import array
+from ROOT import *
+
+fileName, fileNameGeo, outputFileName, sampleType, jobID = None, None, None, None, None
+maxevents = 0
+verbose = True
+# Energy deposit threshold for the liquid scintillator
+threshold = 0.045
+
+# Parse commandline inputs
+try:
+    opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:j:", ['Ethr='])
+except getopt.GetoptError:
+    # print help information and exit:
+    print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents -v verbose (0 or 1) -j jobID (0 to ...) '
+    sys.exit()
+for o, a in opts:
+    if o in ("-f",):
+        fileName = a
+    if o in ("-o",):
+        outputFileName = a
+    if o in ("-n",):
+        print "-n",a
+        maxevents = int(a)
+    if o in ("-j",):
+        print "-j",a
+        jobID = long(a)
+    if o in ("-t",):
+        sampleType = str(a)
+    if o in ("-v",):
+        verbose = bool(int(a))
+    if o in ("--Ethr=",):
+        threshold = float(a)
+    
+fileNameGeo = fileName.replace('ship', 'geofile_full').replace("_rec","")
+# fileNameGeo = "geofile_full.10.0.newGen_15340000_nuBg81.root"#fileName.replace('ship', 'geofile_full').replace("_rec","")
+
+def namestr(obj, namespace=globals()):
+    return [name for name in namespace if namespace[name] is obj]
+
+for item in [fileName, fileNameGeo, outputFileName, sampleType]:
+    if not item:
+        print '\tFATAL! %s not defined!'%namestr(item)[0]
+        sys.exit()
+
+for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose, threshold]:
+    print '\tINFO: %s set to %s'%(namestr(item)[0], item)
+
+# Load SHiP (LHCb) style
+gROOT.ProcessLine(".x mystyle.C")
+# Obsolete debug flags
+oldGeo = False
+ON = True
+# Load PDG database
+pdg = TDatabasePDG.Instance()
+# Add elements to PDG database
+import pythia8_conf
+pythia8_conf.addHNLtoROOT()
+# TODO: MAKE THIS A DICTIONARY AND DE-VERBOSIZE IT!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
+if not(pdg.GetParticle(1000010020)):
+    pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020)
+    pdg.AddAntiParticle("AntiDeuteron", - 1000010020)
+if not(pdg.GetParticle(1000010030)):
+    pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030);
+    pdg.AddAntiParticle("AntiTriton", - 1000010030);
+if not(pdg.GetParticle(1000020040) ):
+    pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040);
+    pdg.AddAntiParticle("AntiAlpha", - 1000020040);
+if not(pdg.GetParticle(1000020030) ):
+    pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030);
+    pdg.AddAntiParticle("AntiHE3", -1000020030);
+if not(pdg.GetParticle(1000030070) ):
+    print "TO BE CHECKED the data insert for Li7Nucleus"
+    pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070);
+if not(pdg.GetParticle(1000060120) ):
+    print "ERROR: random values insert for C12Nucleus"
+    pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120);
+if not(pdg.GetParticle(1000030060) ):
+    print "ERROR: random values insert for Li6Nucleus"
+    pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060);       
+if not(pdg.GetParticle(1000070140) ):
+    print "ERROR: random values insert for for N14"
+    pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140);       
+if not(pdg.GetParticle(1000050100) ):
+    print "TO BE CHECKED the data insert for B10"
+    pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100);
+if not(pdg.GetParticle(1000020060) ):
+    print "TO BE CHECKED the data insert for He6"
+    pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060);
+if not(pdg.GetParticle(1000040080) ):
+    print "TO BE CHECKED the data insert for Be8"
+    pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080);
+if not(pdg.GetParticle(1000030080) ):
+    print "TO BE CHECKED the data insert for Li8"
+    pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080);    
+if not(pdg.GetParticle(1000040170) ):
+    print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!"
+    pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170);    
+if not(pdg.GetParticle(1000040100) ):
+    print "TO BE CHECKED the data insert for Be10"
+    pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100);    
+if not(pdg.GetParticle(1000040070) ):
+    print "TO BE CHECKED the data insert for Be7"
+    pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070);    
+if not(pdg.GetParticle(1000230470) ):
+    print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!"
+    pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470);    
+if not(pdg.GetParticle(1000080170) ):
+    print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!"
+    pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170);    
+if not(pdg.GetParticle(1000240500) ):
+    print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!"
+    pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500);    
+if not(pdg.GetParticle(1000210450) ):
+    print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!"
+    pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450);    
+if not(pdg.GetParticle(1000040090) ):
+    print "TO BE CHECKED the data insert for Be9"
+    pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090);    
+if not(pdg.GetParticle(1000080160) ):
+    print "TO BE CHECKED the data insert for O16"
+    pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160);    
+if not(pdg.GetParticle(1000220460) ):
+    print "TO BE CHECKED the data insert for Ar40"
+    pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460);   
+if not(pdg.GetParticle(1000050110) ):
+    print "TO BE CHECKED the data insert for B11"
+    pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110);   
+
+def PrintEventPart(particles,pdg):
+    print "Particles in the event:"
+    for (pid,part) in enumerate(particles):
+        print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId()
+    print
+
+def PrintRPChits(rpc,pdg):
+    print "Hits in:"
+    for (ix,RPCpt) in enumerate(rpc):
+        tmpName = pdg.GetParticle(RPCpt.PdgCode())
+        print ix,RPCpt.GetZ(), tmpName, RPCpt.GetTrackID() , fGeo.FindNode(RPCpt.GetX(),RPCpt.GetY(),RPCpt.GetZ()).GetVolume().GetName()
+    print
+
+# weight computation moved to offlineForBarbara.py
+
+def addVect(t,name,vectType):
+    vect =  vector(vectType)()
+    t.Branch( name, vect )
+    return t, vect
+
+def addVar(t,name,varType):
+    var = array(varType,[-999])
+    t.Branch(name,var,name+"/"+varType.upper())
+    return t, var
+
+def putToZero(var):
+    var[0] = 0
+    
+def getPartName(partId):
+    return pdg.GetParticle(partId).GetName()
+
+def lookingForDecay(particles):
+    decayMotherList = []
+    for p in particles:
+        mID = p.GetMotherId()
+        if mID>=0 and not (mID in decayMotherList): 
+            decayMotherList.append(mID)
+            decayPartIndex.push_back(mID)
+            decayPartID.push_back(particles[mID].GetPdgCode())
+            decayPartStrID.push_back(pdg.GetParticle(particles[mID].GetPdgCode()).GetName())
+            decayPos_x.push_back(p.GetStartX())
+            decayPos_y.push_back(p.GetStartY())
+            decayPos_z.push_back(p.GetStartZ())
+            decayMother.push_back(particles[mID].GetMotherId())
+            decayP.push_back(p.GetP())
+
+def wasFired(indexKids, detPoints, detPos, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):
+    #
+    def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId):
+        charges = []
+        trackids = []
+        if hitCharges is None:
+            assert(hitTrackId is None)
+        for pos in detPos:
+            foundStat = False
+            foundStat_eff = False
+            for hit in detPoints:
+                if (indexKids is None) or (hit.GetTrackID() in indexKids):
+                    # check if it is in one of the considered active layer
+                    if pos[0]<=hit.GetZ()<=pos[1]:
+                        Eloss += hit.GetEnergyLoss()
+                        hitID = hit.GetTrackID()
+                        if not hitID in trackids and not hitCharges is None:
+                            if  hit.PdgCode()>100000:
+                                charges.append(9)
+                            else:
+                                charges.append(int(pdg.GetParticle(hit.PdgCode()).Charge()))
+                            trackids.append(hitID)
+                            hitCharges.push_back(charges[-1])
+                            hitTrackId.push_back(trackids[-1])
+                        if pointVects is not None:
+                            pointVects[0].push_back(hit.GetX())
+                            pointVects[1].push_back(hit.GetY())
+                            pointVects[2].push_back(hit.GetZ())
+                        flag = True
+                        foundStat = True
+                        eff_val = gRandom.Uniform(0.,1.)
+                        if eff_val<0.9:
+                            flag_eff = flag_eff or True
+                            foundStat_eff = True
+                        else:
+                            flag_eff = flag_eff or False
+            if foundStat:
+                nstat+=1
+            if foundStat_eff:
+                nstat_eff+=1
+        particles = 0
+        flag_Ethr = Eloss>=Ethr
+        return flag, flag_eff, nstat, nstat_eff, Eloss, flag_Ethr,particles
+    #
+    # Now in partKidTrackID I should have the trackID connected to my charged particles
+    flag_eff = False
+    flag = False
+    nstat = 0
+    nstat_eff = 0
+    Eloss = 0
+    flag,flag_eff,nstat,nstat_eff,Eloss,flag_Ethr,particles = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,indexKids,Eloss,Ethr,hitCharges, hitTrackId)
+    if flag==False and checkOn:
+        print "To be study event %s"%entry        
+    return flag, flag_eff, nstat, nstat_eff, flag_Ethr, particles
+
+def convertion(tmpName):
+    if "Rib" in tmpName: tmpName = "Rib"
+    elif "LiSc" in tmpName: tmpName= "LiSc"
+    elif "Tr" in tmpName: tmpName = "Tr"
+    elif "Startplate" in tmpName: tmpName = "Startplate"
+    elif re.search("T\d+O", tmpName): tmpName = "TO"
+    elif re.search("T\d+I", tmpName): tmpName = "TI"
+    elif "Endplate" in tmpName: tmpName = "Endplate"
+    elif "volCoil" in tmpName: tmpName = "volCoil"
+    elif "volFeYoke" in tmpName: tmpName = "volFeYoke"
+    elif "gas" in tmpName: tmpName = "gas"
+    elif "wire" in tmpName: tmpName == "wire"
+    elif "VetoTimeDet" in tmpName: tmpName = "upstreamVeto"
+    elif "Veto" in tmpName: tmpName = "strawVeto"
+    return tmpName
+
+def wasFired_node(indexKids, detPoints, detVolNames, hitCharges, hitTrackId, checkOn, pointVects=None, Ethr=0):
+    #
+    def lookingForHits(detPoints,flag,hitCharges, hitTrackId):
+        if hitCharges is None:
+            assert(hitTrackId is None)
+        for hit in detPoints:
+            if (indexKids is None) or (hit.GetTrackID() in indexKids):
+                tmpName = fGeo.FindNode(hit.GetX(),hit.GetY(),hit.GetZ()).GetVolume().GetName()
+                tmpName = convertion(tmpName)
+                if hit.GetZ()>8000:
+                    continue
+                if tmpName in detVolNames:
+                    ## trick to remove the case of the tracking system the hits in the gas or straw from the straw-veto:
+                    if 'trackingSystem' in detVolNames and hit.GetZ()<0:
+                        continue
+                    if "strawVetoSystem" in detVolNames and hit.GetZ()>0:
+                        continue
+                    if "upstreamVeto" in detVolNames and not(hit.GetZ()>= upstreamVetoPos[0][0] and hit.GetZ()<=upstreamVetoPos[0][1]):    
+                        continue
+                    # check if it is in one of the considered active layer
+                    if True: #pos[0]<=hit.GetZ()<=pos[1]:
+                        flag = True
+        return flag
+    #
+    # Now in partKidTrackID I should have the trackID connected to my charged particles
+    flag = False
+    flag = lookingForHits(detPoints,flag,hitCharges, hitTrackId)
+    if flag==False and checkOn:
+        print "To be study event %s"%entry        
+    return flag
+
+
+""" Main program """
+
+# Open file
+f = TFile(fileName)
+t = f.Get("cbmsim")
+
+totentries = t.GetEntries()
+if verbose:
+    print
+    print
+    print "Program started, reading %s from %s"%(t.GetName(), fileName)
+if (maxevents>0) and (maxevents < totentries):
+    entries = maxevents
+else:
+    entries = totentries
+if verbose:
+    print "Requested to process %s events out of %s in tree"%(entries, totentries)
+
+# Load geometry
+from lookAtGeo import *
+r = loadGeometry(fileNameGeo)
+fGeo = r['fGeo']
+
+# Functions by Elena for offline selection
+# (Also includes another dictionary of all geometry nodes)
+import offlineForBarbara as offline
+offline.loadNodes(fGeo)
+offline.ShipGeo = r['ShipGeo']
+offline.initBField(fileNameGeo)
+offline.sh = offline.StrawHits(t, offline.modules, offline.ShipGeo.straw.resol, 0, None, offline.ShipGeo)
+offline.pdg = pdg
+
+# Handle geometry
+# Names of useful geometry nodes
+myNodes_name = ["VetoTimeDet_1"]
+myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
+myNodes_name += ["Veto_5"]
+# Positions of selected nodes
+myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name,None)
+# Places where a neutrino can interact
+lastPassive_nodeName = ["volIron_23"]
+OPERA_nodeName = ["volIron", "volRPC", "volHPT","volArm2MS","volFeYokes","volCoil","volFeYoke"]#["volIron_%s"%i for i in xrange(12,24)]+["volRpc_%s"%i for i in xrange(11,22)]+["volArm2MS_1"] 
+entrance_nodeName = ["T1Lid"]#['T1lid_1']
+volumeIn_nodeName = ["TI"]#['T%sI'%i for i in [1,2,3,5]]
+volumeOut_nodeName = ["TO"]#['T%sO'%i for i in [1,2,3,5]]
+OPERA_volNames = ["volIron","volRpc","volHPT"]
+strawVeto_volNames = ["Veto"]
+# Z positions of the VETO and tracking systems
+Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']]
+upstreamVeto = [myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]
+trackStationsPos = [[myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ'],myGeoEl["Tr%s_%s"%(i,i)]['z']+myGeoEl["Tr%s_%s"%(i,i)]['dimZ']] for i in xrange(1,5)]
+strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]]
+vetoWall = [[myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+0.001,myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ']]]#myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']+6000.]]
+upstreamVetoPos = [[myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ'],myGeoEl["VetoTimeDet_1"]['z']+myGeoEl["VetoTimeDet_1"]['dimZ']]]
+RPCstationsPos = [[-3500,myGeoEl["VetoTimeDet_1"]['z']-myGeoEl["VetoTimeDet_1"]['dimZ']-0.5]]
+# Places where a neutrino can interact (some duplicates?)
+trackStationsPos_node = ["Tr", "gas", "straw", 'trackingSystem']
+## should be that this does not work since it could be that hits are also assigned to gas or straw
+strawVetoPos_node = ["strawVeto", "gas", "straw", "strawVetoSystem"]
+vetoWall_node = ["LiSc","cave","Rib","TI","TO","Tr","strawVeto"]
+upstreamVetoPos_node = ["upstreamVeto","cave"]
+RPCstationsPos_node = ["volRpc","volMagneticSpectrometer","volIron","volHPT"]
+# Printout positions
+print "OPERApos: ",RPCstationsPos
+print "trackingPos: ",trackStationsPos
+print "strawVetoPos: ",strawVetoPos
+print "scintPos: ",vetoWall
+print "upstreamVetoPos: ",upstreamVetoPos
+
+# Prepare output ntuple
+nf = TFile(outputFileName,"RECREATE")
+nt = TTree("t","t")
+# Event number --------------------------------------------
+nt, event = addVar(nt, 'event','i')
+# Neutrino information ------------------------------------
+nt, nNu = addVar(nt, 'nNu', 'i')
+nt, isPrimary = addVar(nt, 'isPrimary','i')
+nt, startZ_nu = addVar(nt, 'startZ_nu', 'f')
+nt, startY_nu = addVar(nt, 'startY_nu', 'f')
+nt, startX_nu = addVar(nt, 'startX_nu', 'f')
+nt, nuE = addVar(nt, 'nuE', 'f')
+nt, weight = addVar(nt, 'weight', 'f')
+# Interaction element code --------------------------------
+## 0: last passive material opera-mu-system
+## 1: two windows around liquid scintillator
+## 2: along vaccum tank
+## 3: along all OPERA
+## 4: between the two entrance windows
+## 5: along vacuum tank outer window
+## 6: along vacuum tank inner window
+## 999: wrong OPERA place ## needed until we have the new production
+## -1: anything else, but what???? #it should never been present
+nt, interactionElement = addVar(nt, 'interactionElement','i')
+# Neutrino daughters --------------------------------------
+# Do the particles come from a neutrino?
+nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i')
+nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int')
+#nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string')
+# Do charged particles come from a neutrino?
+nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i')
+nt, idChargedPart_fromNu =  addVect(nt, 'idChargedPart_fromNu', 'int')
+#nt, idStrChargedPart_fromNu =  addVect(nt, 'idStrChargedPart_fromNu', 'string')
+# Do neutral particles come from a neutrino?
+nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i')
+nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int')
+#nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string')
+# RPC -----------------------------------------------------
+# In case something (no matter if it comes from the nu-interaction kids) fired the RPC
+nt, RPCany = addVar(nt,'RPCany', 'i' )
+# accounting for an eff. of each station of 90%
+nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' )
+# upstreamVeto --------------------------------------------
+nt, upstreamVetoany = addVar(nt,'upstreamVetoAny', 'i' )
+# accounting for an eff. of each station of 90%
+nt, upstreamVetoany_eff = addVar(nt,'upstreamVetoAny_eff', 'i' )
+# scintVeto -----------------------------------------------
+## In case something (no matter if it comes from the nu-interaction kids) fired the surroundin walls
+nt, scintVetoAny = addVar(nt,'scintVetoAny', 'i' )
+nt, scintVetoAny_eff = addVar(nt,'scintVetoAny_eff', 'i' )
+# strawVeto -----------------------------------------------
+# In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto
+nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' )
+nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' )
+# TrackSyst -----------------------------------------------
+# only the case of anything fired it is considered, obviously
+nt, TrackSyst = addVar(nt, "TrackSyst", 'i')
+nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i')
+# VETO and tracking systems with thresholds ---------------
+nt, strawVetoAny_Ethr = addVar(nt,"strawVetoAny_Ethr","i")
+nt, scintVetoAny_Ethr = addVar(nt,"scintVetoAny_Ethr","i")
+nt, upstreamVetoAny_Ethr = addVar(nt,"upstreamVetoAny_Ethr","i")
+nt, RPCany_Ethr = addVar(nt,"RPCany_Ethr","i")
+nt, TrackSyst_Ethr = addVar(nt, "TrackSyst_Ethr", "i")
+# Decay information ---------------------------------------
+nt, decayPartIndex = addVect(nt, "decayPartIndex", "float")
+nt, decayPartID = addVect(nt, "decayPartID", "float")
+nt, decayPartStrID = addVect(nt, "decayPartStrID", "string")
+nt, decayPos_x = addVect(nt, "decayPos_x", "float")
+nt, decayPos_y = addVect(nt, "decayPos_y", "float")
+nt, decayPos_z = addVect(nt, "decayPos_z", "float")
+nt, decayMother = addVect(nt, "decayMother", "float")
+nt, decayP = addVect(nt,"decayP","float")
+# Is this a NC interaction? -------------------------------
+nt, NC = addVar(nt, "NC", "i")
+# Store where the neutrino interacted ---------------------
+#nt, nuInteractionNode = addVect(nt, "nuInteractionNode", "string")
+nt, nuIntNumSimpl = addVar(nt,"nuIntNumSimpl","i")
+dictNodeNames = {'volIron':0, 'cave':1, 'LiSc':2, 'Startplate':3, 'TI':4, 'rockD':5, 'Endplate':6, 'Rib':7, 'volFeYoke':8, 'volHPT':9, 'TO':10, 'volRpc':11, 'volCoil':12, 'T1Lid':13, 'straw':14, 'strawVeto':15, 'volBase':16, 'Tr':17, 'gas':18, 'wire':19, 'others':20}
+# Was the event reconstructed? ----------------------------
+nt, recoed = addVar(nt, "recoed","i")
+# How many candidates in the reco event? ------------------
+nt, nRecoed = addVar(nt, "nRecoed","i")
+# Add stuff for online selection --------------------------
+offline.addOfflineToTree(nt)
+# Leaf for job ID
+nt, jID = addVar(nt, "jobID", "l")
+
+####### New Variables after TP
+# original interaction element
+nt, originalIntEl = addVect(nt, "originalIntEl","string")
+# PDG Code of Particles with Hits in Tracking system
+nt, PdgPartTrackSyst  = addVect(nt, "PdgPartTrackSyst","int")
+# It's Momentum
+nt, PPartTrackSyst = addVect(nt, "PPartTrackSyst","float")
+# It's Energy
+nt, EPartTrackSyst = addVect(nt, "EPartTrackSyst","float")
+# Where it interacted
+nt, intElPartTrackSyst = addVect(nt, "intElPartTrackSyst","string")
+# number of Hits left by the particle in the Tracking system per particle
+nt, nHitsPartTrackSyst = addVect(nt, "nHitsPartTrackSyst","int")
+# number of Hits per Event
+nt, nHitsTrackSystEvent = addVar(nt, "nHitsTrackSystEvent","i")
+# PDG Code of mother of these particles
+nt, PdgMotherTrackSyst = addVect(nt, "PdgMotherTrackSyst","int")
+# interaction Element of the mother
+nt, intElMotherTrackSyst = addVect(nt, "intElMotherTrackSyst", "string")
+# Energy mother
+nt, EMotherTrackSyst = addVect(nt,"EMotherTrackSyst","float")
+
+
+
+
+# Now loop on the tree
+# t = original tree
+# nt = new tree
+for entry in xrange(entries):
+    if not (entry%1000):
+        print "\tProcessing entry %s of %s..."%(entry,entries)
+    t.GetEntry(entry)
+    putToZero(jID)
+    jID[0] = jobID
+    event[0] = entry #+(jobID*entries)
+
+    ####################
+    particles = t.MCTrack
+    rpc = t.ShipRpcPoint
+    scint = t.vetoPoint
+    strawPoints = t.strawtubesPoint
+    vetoScint = t.vetoPoint
+    recoParts = t.Particles
+    trackFit = t.FitTracks
+    ####################
+    
+    # initialize containers
+    putToZero(nNu)
+    putToZero(nPart_fromNu)
+    putToZero(nChargedPart_fromNu)
+    putToZero(nNeutrPart_fromNu)
+    idPart_fromNu.clear()
+    idChargedPart_fromNu.clear()
+    idNeutrPart_fromNu.clear()
+    #idStrPart_fromNu.clear()
+    #idStrChargedPart_fromNu.clear()
+    #idStrNeutrPart_fromNu.clear()
+    decayPartIndex.clear()
+    decayPartID.clear()
+    decayPartStrID.clear()
+    decayPos_x.clear()
+    decayPos_y.clear()
+    decayPos_z.clear()
+    decayMother.clear()
+    decayP.clear()
+
+    ## new vars
+    originalIntEl.clear()
+    PdgPartTrackSyst.clear()
+    PPartTrackSyst.clear()
+    EPartTrackSyst.clear()
+    intElPartTrackSyst.clear()
+    nHitsPartTrackSyst.clear()
+    putToZero(nHitsTrackSystEvent)
+    PdgMotherTrackSyst.clear()
+    intElMotherTrackSyst.clear()
+    EMotherTrackSyst.clear()
+
+    
+    #nuInteractionNode.clear()
+    primaryDone=False
+    isPrimary[0] = int(False)
+    lookingForDecay(particles)
+    nRecoed[0] = 0
+    recoed[0] = 0
+    # Which one is the "primary" particle?
+    if sampleType=="nuBg" or sampleType == "cosmics":
+        startPartRange = [0] 
+        primaryMum = -1
+    elif sampleType=="sig":
+        startPartRange = [1] 
+        primaryMum = 0
+    # Look at the primary particle
+    ip = startPartRange[0]
+    skipEvent=False
+    # Were there any hit in the tracking stations?
+    TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, None, None, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] )
+    # If no, skip this event
+    if TrackSyst[0]==0:
+        skipEvent=True
+        #print "i am jumping this one"
+        continue
+        # exit from loop on particles
+        #break
+    # Select the primary MC particle
+    part = particles[ip]
+    pdgPart = pdg.GetParticle(part.GetPdgCode())
+    if sampleType=="nuBg":
+        assert("nu" in pdgPart.GetName())        
+    ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother
+    #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena
+    ## (and de-indented the following block)
+    ## Starting the counter of how many particles were produced by the interaction of this specific nu
+    for recoP in recoParts:
+        nRecoed[0] += 1
+        #offline.pushOfflineByParticle(t, recoP)
+    if nRecoed[0]>0:
+        recoed[0]=1
+    else:
+        skipEvent=True
+        continue
+        #break
+    nNu[0]+=1
+    if part.GetMotherId()==primaryMum:
+        assert(primaryDone==False)
+        primaryDone=True
+        isPrimary[0] = int(True) 
+        assert(len(idPart_fromNu)==0)
+        assert(len(idNeutrPart_fromNu)==0)
+    else: 
+        continue
+    # Where did this guy interact?
+    tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName()
+    #nuInteractionNode.push_back(tmpName)
+    tmpName = convertion(tmpName)
+    if not tmpName in dictNodeNames:
+        tmpName = "others"
+    nuIntNumSimpl[0] = dictNodeNames[tmpName]
+    originalIntEl.push_back(tmpName)
+    # Store interaction point coordinates
+    startZ_nu[0] = part.GetStartZ()
+    startY_nu[0] = part.GetStartY()
+    startX_nu[0] = part.GetStartX()
+    # Primary particle information
+    nuE[0] = part.GetEnergy()
+    # Looking for particles produced by the neutrino interaction
+    interacted = False
+    partKidsId = []
+    NC[0] = 0
+
+    # Add information for offline selection
+    # (mainly signal normalisation and zeroing of particle-wise arrays)
+    ntr, nref = offline.pushOfflineByEvent(t, vetoScint, sampleType, verbose, threshold)
+    for recoP in recoParts:
+       offline.pushOfflineByParticle(t, recoP, ntr, nref)
+
+    # Loop on MCTracks
+    for ip2 in xrange(0,len(particles)):
+        part2 = particles[ip2]
+        # exit if we have reached the empty part of the array
+        if not (type(part2)==type(ShipMCTrack())):
+            break
+        if part2.GetMotherId()==ip:
+            interacted = True 
+            part2Id = part2.GetPdgCode()
+            partKidsId.append(part2Id)
+            # Was this a neutral current interaction?
+            if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0):
+                NC[0] = 1
+            nPart_fromNu[0]+=1
+            idPart_fromNu.push_back(part2Id)
+            if pdg.GetParticle(part2.GetPdgCode()) == None:
+                part2Name = str(part2Id)
+            else:
+                part2Name = getPartName(part2Id)
+            #idStrPart_fromNu.push_back(part2Name)
+            # Counting how many particles produced by the nu-interaction are charged or neutral
+            if (pdg.GetParticle(part2.GetPdgCode())) == None or (int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0)):
+                nNeutrPart_fromNu[0]+=1
+                idNeutrPart_fromNu.push_back(part2Id)
+                #idStrNeutrPart_fromNu.push_back(part2Name)
+            else:
+                nChargedPart_fromNu[0]+=1
+                idChargedPart_fromNu.push_back(part2Id)
+                #idStrChargedPart_fromNu.push_back(part2Name)
+            # How many are produced by the interaction in the last passive element of the "opera-mu system"?
+            # Finding out the "interaction element code"
+            part2Z = part2.GetStartZ()
+            part2X = part2.GetStartX()
+            part2Y = part2.GetStartY()
+            somewhere = False
+            nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName()
+            if  nodeName in lastPassive_nodeName:
+                somewhere = True
+                interactionElement[0] = 0 
+                intElName = 'OPERA'
+            # To know if it was in the full OPERA-system (excluded last passive)
+            if tmpName in OPERA_nodeName and somewhere==False:
+                somewhere = True
+                interactionElement[0] = 3
+                intElName = 'OPERA'
+            # To know if it was between the two windows
+            if tmpName in entrance_nodeName:
+                assert(somewhere==False)
+                somewhere = True
+                interactionElement[0] = 1
+            # Vacuum tank outer window
+            if tmpName in volumeIn_nodeName and somewhere==False:
+                interactionElement[0] = 5
+                somewhere = True 
+            # Vacuum tank inner window
+            elif nodeName in volumeOut_nodeName and somewhere==False:
+                interactionElement[0] = 6
+                somewhere = True 
+            # Vacuum tank ribs
+            if "Rib" in tmpName:
+                interactionElement[0] = 7
+            # Somewhere else
+            if somewhere==False:
+                interactionElement[0] = -1
+        # End loop on MCTracks
+    # Should be changed to avoid entering in the loop if there isn't anything in the tracking
+    if skipEvent:
+        continue
+    strawVetoAny[0],strawVetoAny_eff[0],dummy,dummy,strawVetoAny_Ethr[0],dummy  = wasFired(None, strawPoints, strawVetoPos,None,None, checkOn=False, pointVects=None)#[strawVetoPoint_x, strawVetoPoint_y, strawVetoPoint_z])
+    upstreamVetoany[0], upstreamVetoany_eff[0], dummy,dummy, upstreamVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, upstreamVetoPos, None, None, checkOn=False, pointVects=None)#[upstreamVetoPoint_x, upstreamVetoPoint_y, upstreamVetoPoint_z])
+    RPCany[0], RPCany_eff[0],dummy,dummy, RPCany_Ethr[0],dummy = wasFired(None, rpc, RPCstationsPos, None, None, checkOn=False, pointVects=None)#[rpcPoint_x, rpcPoint_y, rpcPoint_z])
+    scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=threshold)
+    #assert(len(idStrPart_fromNu)==len(idPart_fromNu))
+    #assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu))
+    #assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu))
+    #assert(len(idStrPart_fromNu)==nChargedPart_fromNu[0]+nNeutrPart_fromNu[0])
+    #assert(len(idStrChargedPart_fromNu)==nChargedPart_fromNu[0])
+    #assert(len(idStrNeutrPart_fromNu)==nNeutrPart_fromNu[0])
+    if not primaryDone:
+        print "It looks like there is no primary nu interacting"
+        PrintEventPart(particles,pdg)
+        sys.exit()
+    weight[0] = offline.findWeight(sampleType, NC[0], nuE[0], part, entries, pdgPart.GetName(), ON)#calcWeight(NC[0], nuE[0], part.GetWeight(), entries, pdgPart.GetName(), ON)
+
+
+    # Fill Particle in Tracking System Variables
+    # helping variables: hitlist and tracklist
+    hitlist={}
+    for hit in strawPoints:
+        trId = hit.GetTrackID()
+        if (trId in hitlist.keys()):
+            hitlist[trId] += 1
+        else:
+            hitlist[trId] = 1
+    # Save the trackIDs of the FitTracks in a list
+    tracklist = []
+    for i in xrange(len(trackFit)):
+        tracklist.append(t.fitTrack2MC[i])
+    # Fill the information to the Tree
+    for trId in tracklist:
+        part = particles[trId]
+        mother = particles[part.GetMotherId()]
+        PdgPartTrackSyst.push_back(part.GetPdgCode())
+        PdgMotherTrackSyst.push_back(mother.GetPdgCode())
+        EPartTrackSyst.push_back(part.GetEnergy())
+        EMotherTrackSyst.push_back(mother.GetEnergy())
+        tmpPart = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName()
+        tmpPart = convertion(tmpPart)
+        tmpMother = fGeo.FindNode(mother.GetStartX(),mother.GetStartY(),mother.GetStartZ()).GetVolume().GetName()
+        tmpMother = convertion(tmpMother)
+        intElPartTrackSyst.push_back(tmpPart)
+        intElMotherTrackSyst.push_back(tmpMother)
+        PPartTrackSyst.push_back(part.GetP())
+        nHitsPartTrackSyst.push_back(hitlist[trId])
+        nHitsTrackSystEvent[0] += hitlist[trId]
+
+    
+    nt.Fill()
+    # End loop on tree        
+nt.Write()
+nf.Save()
+nf.Close()
+f.Close()
+print "Output wrote to %s"%outputFileName
+print "Number of events with vertex outside Veto_5-500 - Tr4_4: %s"%offline.num_bad_z
diff --git a/newGen/offlineForBarbara.py b/newGen/offlineForBarbara.py
index 56468eb..e855618 100644
--- a/newGen/offlineForBarbara.py
+++ b/newGen/offlineForBarbara.py
@@ -1,669 +1,670 @@
-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.EcalPointLite:
-        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
+    return 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.EcalPointLite:
+        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