diff --git a/distribsForHNLandBG_byEvent.py b/distribsForHNLandBG_byEvent.py
new file mode 100755
index 0000000..dcca120
--- /dev/null
+++ b/distribsForHNLandBG_byEvent.py
@@ -0,0 +1,721 @@
+import sys
+import ROOT as r
+import numpy as np
+import funcsByBarbara as tools
+
+r.gROOT.SetBatch(r.kTRUE)
+r.gROOT.ProcessLine(".X lhcbStyle.C")
+plots = {}
+pdg = r.TDatabasePDG.Instance()
+
+def addPlot(key, x_axis, tree, string, cuts=None):
+	if not cuts:
+		tree.Draw(string+'>>h1')
+	else:
+		tree.Draw(string+'>>h1', cuts)
+	plots[key] = r.gDirectory.Get('h1')
+	plots[key].GetXaxis().SetTitle(x_axis)
+	plots[key].SetTitle(key)
+	plots[key].SetName(key)
+	c1 = r.TCanvas(); c1.cd()
+	if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+		c1.SetLogy()
+	plots[key].Draw()
+	c1.Print('accPlots/%s.pdf'%key)
+	c1.Print('accPlots/%s.png'%key)
+	key += "-WEIGHTED"
+	if 'bg' in key or ('nu' in key and key != 'mumunu'):
+		if not cuts:
+			tree.Draw(string+'>>h1', "NuWeight")
+		else:
+			cuts = "NuWeight*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+	elif 'cosmics' in key:
+		if not cuts:
+			tree.Draw(string+'>>h1',"HNLw*1.e6/200.")
+		else:
+			cuts = "HNLw*1.e6/200.*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+	else:
+		if not cuts:
+			tree.Draw(string+'>>h1',"HNLw")
+		else:
+			cuts = "HNLw*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+
+
+
+
+fpimu = r.TFile('../DATA/NewPIMU/ShipAna.root','read')#r.TFile('../DATA/MUPI/ShipAna.root','read')
+pimu = fpimu.Get('ShipAna')
+pimu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewPIMU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root')#'../DATA/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-2985556.root')
+
+fmumunu = r.TFile('../DATA/NewMUMUNU/ShipAna.root','read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read')
+mumunu = fmumunu.Get('ShipAna')
+mumunu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewMUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-mumunu.root')#'../DATA/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-TOTAL.root')
+
+ftiny = r.TFile('../DATA/SmallHNLs/ShipAna.root', 'read')
+tiny = ftiny.Get('ShipAna')
+tiny_geo = tools.searchForNodes3_xyz_dict('../DATA/SmallHNLs/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-011.root')
+
+#fbg = r.TFile('../DATA/KLONG/ShipAna.root','read')
+#bg = fbg.Get('ShipAna')
+#bg_geo = tools.searchForNodes2_xyz_dict('../DATA/KLONG/geofile_full.10.0.Genie-TGeant4_70cc30nc.root')
+
+fcosmics = r.TFile('../DATA/Cosmics/ShipAna.root','read')
+cosmics = fcosmics.Get('ShipAna')
+cosmics_geo = tools.searchForNodes2_xyz_dict('../DATA/Cosmics/geofile_full.Signal.10.0.Cosmics3.25.1-8.root')
+
+fnubar = r.TFile('../DATA/NewKLONG/ShipAna_nubar.root','read')#'../DATA/NewKLONG/ShipAna74-75_nubar.root','read')
+nubar = fnubar.Get('ShipAna')
+nubar_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nubar.root')#neutrino74-75.root')
+
+fnu = r.TFile('../DATA/NewKLONG/ShipAna_nu.root','read')#'../DATA/NewKLONG/ShipAna77_nu.root','read')
+nu = fnu.Get('ShipAna')
+nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nu.root')#neutrino77.root')
+
+
+# For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume
+
+studies = {
+	'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-4175},#1000},
+	'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-4125},#7000},
+	'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-4498},
+	#'bg': {'data':bg, 'geo':bg_geo, 'ntot':142856},#99999}#199998
+	'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'seen':[], 'ntot':cosmics.GetEntriesFast()},
+	'nubar': {'data': nubar, 'geo':nubar_geo, 'seen':[], 'ntot':99999*2},#99999*2},
+	'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':99999*2}#99999}
+}
+
+print 'Found %s entries in the cosmics file'%cosmics.GetEntriesFast()
+
+# Cuts:
+# - z_vertex in [vol[0]+5m, straw[0]]
+# - has track converged
+# - chi2/ndof < 10
+# - DOCA < 10 cm
+# - has no veto?
+
+
+
+class Cuts(object):
+    def __init__(self, cutConfigList, study): 
+        self.cutList = []      
+        try:        
+            for (index,cutConfig) in enumerate(cutConfigList):
+                self.cutList.append( self.Cut( self, index, cutConfig['name'], study, cutConfig['params'] ) )               
+        except BaseException,e:
+            print "cutConfigList badly set: %r" %e
+    def process(self, eventList):
+        for event in eventList:
+            for cut in self.cutList:
+                if  not cut.count(event): break
+        return self
+    def select(self, event, imax):
+    	for i in xrange(imax):
+    		if not self.cutList[i].selection(event, self.cutList[i].study, self.cutList[i].params):
+    			return False
+    	return True
+
+    class Cut(object):
+        def __init__(self, cutCollection, index, cutName, study, params):
+        	self.cutCollection = cutCollection
+        	self.index = index
+        	self.name = cutName
+        	self.params = params
+        	self.study = study
+        	self.total = 0 
+        	self.weight = 0
+        	self.eventNumbers = []
+        	self.duplicates = []
+        	self.selection = self.__get(cutName)
+        def __get(self, cutName):        
+            name = "cut_"+cutName
+            if hasattr(self,name):
+                try:
+                    return getattr(self,name)
+                except AttributeError,e:
+                    return None
+            else:
+                raise BaseException ( "bad cut name: %r" %cutName )
+                return None
+        def count(self, event):
+            try:
+                if (self.selection(event,self.study,self.params)):
+                    self.total += 1
+                    if (event.EventNumber not in self.eventNumbers):
+                        self.weight += self.compute_weight(event, self.study)
+                        self.eventNumbers.append(event.EventNumber)
+                    else:
+                    	self.duplicates.append(event.EventNumber)
+                    return True
+                else :
+                    return False
+            except BaseException,e:
+                print "Error performing %r: %r" %(self.name,e)
+        def cut_fit_converged(self, event, study, params):
+        	return fit_converged(event)
+        def cut_chi2(self, event, study, params):
+            return chi2_cut(event,params['chi2max'])
+        def cut_has_muons(self, event, study, params):
+        	for station in params['stations']:
+        		if not has_muon_station(event, study, station):
+        			return False
+        	return True
+        def cut_doca(self, event, study, params):
+        	return doca_cut(event, params['docamax'])
+        def cut_z(self, event, study, params):
+        	return z_cut(event, study)
+        def cut_normalization(self, event, study, params):
+        	return goodTruthVtx(event, study)
+        def cut_good_tracks(self, event, study, params):
+        	return goodTrack(event)
+        def cut_ip(self, event, study, params):
+        	return ip_cut(event, params['ipmax'])
+        def cut_vtx_in_ellipse(self, event, study, params):
+        	return vtxInAcceptance(event, study)
+        def cut_n_candidates(self, event, study, params):
+        	#print self.cutCollection.cutList[self.index-1].name
+        	if (event.EventNumber in self.cutCollection.cutList[self.index-1].duplicates):
+        		#if self.cutCollection.select(event, self.index):
+        		return nParticles_cut(event, params['nmax'])
+       		return True
+        def compute_weight(self, event, study):
+        	return weight(event, study)
+
+
+
+def fit_converged(event):
+	flags = event.DaughtersFitConverged
+	if np.product(flags):
+		return True
+	return False
+
+def chi2_cut(event, chi2max):
+	chi2 = event.DaughtersChi2
+	n = event.DaughtersNPoints
+	for i in xrange(len(n)):
+		if not n[i]:
+			print event.EventNumber
+			return False
+		if chi2[i]/n[i] > chi2max:
+			return False
+	return True
+
+def has_muon_station(event, study, station):
+	zIn = studies[study]['geo']['muondet%s_1'%(station-1)]['z']['pos'] - studies[study]['geo']['muondet%s_1'%(station-1)]['z']['dim']
+	zOut = studies[study]['geo']['muondet%s_1'%(station-1)]['z']['pos'] + studies[study]['geo']['muondet%s_1'%(station-1)]['z']['dim']
+	for z in event.muon_z:
+		if zIn <= z <= zOut:
+			return True
+	return False
+
+def doca_cut(event, doca_max):
+	if event.DOCA < doca_max:
+		return True
+	return False
+
+def z_cut(event, study):
+	z = event.vtxz
+	#try:
+	#	minz = studies[study]['geo']['lidT1I_1']['z']['pos'] + studies[study]['geo']['lidT1I_1']['z']['dim'] + 500.
+	#	maxz = studies[study]['geo']['lidT6I_1']['z']['pos'] - studies[study]['geo']['lidT6I_1']['z']['dim']
+	#except:
+	minz = studies[study]['geo']['Veto_5']['z']['pos'] + studies[study]['geo']['Veto_5']['z']['dim'] + 0.01
+	maxz = studies[study]['geo']['Tr1_1']['z']['pos'] - studies[study]['geo']['Tr1_1']['z']['dim'] - 0.01
+	if minz < z < maxz:
+		return True
+	return False
+
+def goodTruthVtx(event, study):
+	if study in ['pimu', 'mumunu', 'tiny']:
+		if not event.BadTruthVtx:
+			return True
+	return True
+
+def nParticles_cut(event, nmax):
+	if event.NParticles > nmax:
+		return False
+	return True
+
+# This checks that:
+# - there is at least one tracking station before and one after the magnet
+# - hits in these tracking stations are within a 5x10m ellipse
+def goodTrack(event):
+	if event.DaughtersAlwaysIn:
+		return True
+	return False
+
+def ip_cut(event, maxip):
+	if event.IP0 < maxip:
+		return True
+	return False
+
+def nothing_in_veto5(event):
+	if len(event.veto5_x) > 0:
+		return False
+	return True
+
+def nothing_in_liquidscint(event):
+	if len(event.liquidscint_x) > 0:
+		return False
+	return True
+
+def something_in_muon(event):
+	if len(event.muon_x) > 0:
+		return True
+	return False
+
+def something_in_ecal(event):
+	if len(event.ecal_x) > 0:
+		return True
+	return False
+
+def vtxInAcceptance(event, study):
+	x = event.vtxx
+	y = event.vtxy
+	#if 'cosmics' in study or 'bg' in study:
+	#	print study, x, y
+	return tools.inEllipse(x,y,500./2,1000./2)
+
+def weight(event, study):
+	#if event.EventNumber not in studies[study]['seen']:
+		#studies[study]['seen'].append(event.EventNumber)
+	if ('bg' in study) or (study == 'nu' or study == 'nubar'):
+		w = event.NuWeight
+	elif 'cosmics' in study:
+		w = event.HNLw * 1.e6 / 200.
+	elif study in ['pimu', 'mumunu', 'tiny']:
+		w = event.HNLw / studies[study]['ntot']
+	else:
+		print 'WEIGHTING ERROR: ', study
+		sys.exit()
+	#else:
+	#	w = 0.
+		#if study in ['pimu', 'mumunu', 'tiny']:
+		#	print study, event.EventNumber, event.NParticles
+	return w
+
+def computeResidualWeight(filename=None):
+	if not filename:
+		filename = "forElena_nu_notVetoed.txt"
+	fScan = file(filename)
+	lines = fScan.readlines()
+	sample = nu
+	s = 'nu'
+	if ('bar' in filename) or ('anti' in filename):
+		sample = nubar
+		s = 'nubar'
+	lines = lines[3:-1]
+	events = []
+	for line in lines:
+	    events.append(int(line.rsplit("*")[1].replace(" ","")))
+	count_reco = 0
+	count = 0
+	weighted_count = 0.
+	print "unvetoed: ", len(events)
+	print "total Reconstructed: ", sample.GetEntriesFast()
+	for event in sample:
+		if event.EventNumber in events:
+			count_reco +=1
+			if (fit_converged(event) and chi2_cut(event, 5.) and z_cut(event, s) and vtxInAcceptance(event, s)
+				and goodTrack(event) and doca_cut(event, 30) and has_muon_station(event, s, 1)
+				and has_muon_station(event, s, 2) and ip_cut(event, 250.)):
+			#if True:#fit_converged(event): #and doca_cut(event, 50.) and ip_cut(event, 250.) and z_cut(event, 'nu') and vtxInAcceptance(event, 'nu'):
+				#print event.EventNumber, 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+				print event.EventNumber
+				print 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+				#print event.DOCA, event.IP0, event.DaughtersChi2[0]/event.DaughtersNPoints[0], event.DaughtersChi2[1]/event.DaughtersNPoints[1], event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+				print "tracks: ", event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+				print "vtx: ", event.vtxx, event.vtxy, event.vtxz
+				print 'DOCA, IP: ', event.DOCA, event.IP0
+				#print tools.interactionElement(event.DaughtersTruthProdX[0], event.DaughtersTruthProdY[0], event.DaughtersTruthProdZ[0], nu_geo), tools.interactionElement(event.DaughtersTruthProdX[1], event.DaughtersTruthProdY[1], event.DaughtersTruthProdZ[1],nu_geo)
+				count +=1
+				weighted_count += event.NuWeight
+	print "Reconstructed: %s"%count_reco
+	return count, weighted_count
+
+def latex_float(f):
+    float_str = "{0:.3g}".format(f)
+    if "e" in float_str:
+        base, exponent = float_str.split("e")
+        return "${0} \\times 10^{{{1}}}$".format(base, int(exponent))
+    else:
+        return float_str
+
+cutConfigList = [   {'name': 'normalization', 'params':{}},
+					{'name': 'fit_converged', 'params':{}},
+					{'name': 'chi2',          'params':{'chi2max':5.}},
+					{'name': 'z',             'params':{}},
+					{'name': 'vtx_in_ellipse','params':{}},
+					{'name': 'good_tracks',   'params':{}},
+					{'name': 'doca', 		  'params':{'docamax':50.}},
+					{'name': 'ip', 		      'params':{'ipmax':500.}},
+					{'name': 'has_muons',     'params':{'stations':[1,2]}},
+					{'name': 'doca', 		  'params':{'docamax':30.}},
+					{'name': 'ip', 		      'params':{'ipmax':250.}},
+					{'name': 'n_candidates',  'params':{'nmax':1}},
+]
+
+def numbers_v2():
+	mycuts, results = {}, {}
+	for study in studies:
+		print 'Reading %s...'%study
+		mycuts[study] = Cuts(cutConfigList, study).process(studies[study]['data'])
+		for cut in mycuts[study].cutList:
+			print '%s\t%s\t%s\t%s'%(cut.name, cut.total, len(cut.eventNumbers), cut.weight)
+		results[study] = {}
+		results[study]['numbers']=[ cut.total for cut in mycuts[study].cutList[1:] ]
+		results[study]['numbersWithWeights']=[ cut.weight for cut in mycuts[study].cutList[1:] ]
+		print
+		print 'Events with more than one particle:'
+		for entry in mycuts[study].cutList[-2].duplicates: #[ x for x in mycuts[study].cutList[-1].eventNumbers if mycuts[study].cutList[-1].eventNumbers.count(x) > 1 ]:
+			for tr in studies[study]['data']:
+				if tr.EventNumber == entry:
+					if mycuts[study].select(tr, len(mycuts[study].cutList)-1):
+						print tr.EventNumber, tr.vtxz, tr.IP0, [pdg.GetParticle(code).GetName() for code in tr.DaughtersTruthPDG]
+		print
+		print
+	return results
+
+def makeTex(results):
+	print
+	print 'Background:'
+	print '\t\\begin{table}'
+	print '\t\\centering Backgrounds: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['cosmics']['numbersWithWeights'][i]), ' & ', '%s'%results['cosmics']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print
+	print 'Signal:'
+	print '\t\\begin{table}'
+	print '\t\\centering Signal: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['tiny']['numbersWithWeights'][i]), ' & ', '%s'%results['tiny']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print	
+
+
+
+def numbers_for_TP():
+	results = {}
+	for s in ['pimu', 'mumunu', 'tiny', 'nu', 'nubar', 'cosmics']:#studies:
+		print 'Reading %s...'%s
+		t = studies[s]['data']
+		ntot = t.GetEntriesFast()
+		nAfterCuts = [0]*10
+		nAfterCuts_weighted = [0.]*10
+		initialEventNums = []
+		eventNums = []
+		for event in t:
+			w = weight(event,s) # e' importante fare questa cosa a monte, x come e' definita la funzione!
+			if fit_converged(event) and goodTruthVtx(event,s):
+				nAfterCuts[0] +=1
+				nAfterCuts_weighted[0] += w
+				initialEventNums.append(int(event.EventNumber))
+				if chi2_cut(event, 5.):
+					nAfterCuts[1]+=1
+					nAfterCuts_weighted[1]+=w
+					if z_cut(event, s):#chi2_cut(event, 5):
+						nAfterCuts[1+1] += 1
+						nAfterCuts_weighted[1+1] += w
+						if vtxInAcceptance(event, s):#goodTrack(event):
+							nAfterCuts[2+1] +=1
+							nAfterCuts_weighted[2+1] +=w
+							#if 'cosmics' in s:
+							#	print 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+							#	print "tracks: ", event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+							#	print "p candidate: ", event.P, "\t IP: ", event.IP0
+							#	print "pt candidate: ", event.Pt
+							if goodTrack(event):#doca_cut(event, 50):
+								nAfterCuts[3+1] +=1
+								nAfterCuts_weighted[3+1] +=w
+								if doca_cut(event, 50):#vtxInAcceptance(event, s):
+									nAfterCuts[4+1]+=1
+									nAfterCuts_weighted[4+1]+=w
+									if has_muon_station(event, s, 1) and has_muon_station(event, s, 2):#ip_cut(event,500):#something_in_muon(event):
+										nAfterCuts[5+1] +=1
+										nAfterCuts_weighted[5+1] +=w
+										if doca_cut(event, 30.):#something_in_ecal(event):
+											nAfterCuts[6+1] +=1
+											nAfterCuts_weighted[6+1] +=w
+											if ip_cut(event, 500.):
+												nAfterCuts[7+1] +=1
+												nAfterCuts_weighted[7+1] +=w
+												if ip_cut(event, 250.):
+													nAfterCuts[8+1] +=1
+													nAfterCuts_weighted[8+1] +=w
+													eventNums.append(int(event.EventNumber))
+													studies[s]['seen'].append(event.EventNumber)
+
+		print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot
+		print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot
+		results[s] = {}
+		results[s]['ntot']=studies[s]['ntot']
+		results[s]['numbers']=nAfterCuts
+		results[s]['numbersWithWeights']=nAfterCuts_weighted
+		results[s]['listOfSurvivors']=list(set(eventNums))
+		studies[s]['listOfSurvivors']=list(set(eventNums))
+		results[s]['listOfRecoed']=list(set(initialEventNums))
+		studies[s]['listOfRecoed']=list(set(initialEventNums))
+	print
+	for s in ['pimu', 'mumunu', 'tiny']:#studies.keys():#['pimu', 'mumunu', 'bg']:
+	    for lst in results[s]:
+		if not ('Survivors' in lst) and not ('Recoed' in lst):
+		        print s,list,results[s][lst]
+	print
+	print 'Events with surviving nu: %s \t with nubar: %s'%(len(results['nu']['listOfSurvivors']), len(results['nubar']['listOfSurvivors']))
+	print 
+	print 'Background:'
+	print '\t\\begin{table}'
+	print '\t\\centering Backgrounds: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 5 m', 'IP $<$ 2.5 m']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['cosmics']['numbersWithWeights'][i]), ' & ', '%s'%results['cosmics']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print 'Signal:'
+	print '\t\\begin{table}'
+	print '\t\\centering Signal: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 5 m', 'IP $<$ 2.5 m']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['tiny']['numbersWithWeights'][i]), ' & ', '%s'%results['tiny']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print
+	out = {'pimu':[], 'mumunu':[], 'tiny':[]}
+	for s in ['pimu', 'mumunu', 'tiny']:
+		for event in studies[s]['data']:
+			if event.EventNumber in studies[s]['listOfSurvivors']:
+				if event.NParticles > 1:
+					out[s].append(event.EventNumber)
+					#print s, event.EventNumber, event.NParticles, r.gGeoManager.FindNode(event.vtxx, event.vtxy, event.vtxz).GetName()
+		print s, list(set(out[s]))
+
+		
+def badEvents():
+	numbers_for_TP()
+	hnpnu_recoed = r.TH1F('np nu recoed', 'np nu recoed', 60, 0, 59)
+	hnpnu_survivors = r.TH1F('np nu survived', 'np nu survived', 60, 0, 59)
+	seen = []
+	for part in studies['nu']['data']:
+		if not (part.EventNumber in seen):
+			if part.EventNumber in studies['nu']['listOfSurvivors']:
+				hnpnu_survivors.Fill(part.NParticles)
+			if part.EventNumber in studies['nu']['listOfRecoed']:
+				hnpnu_recoed.Fill(part.NParticles)
+			seen.append(part.EventNumber)
+	#c1 = r.TCanvas(); c1.cd()
+	#hnpnu.Draw()
+	#c1.Print('badEvents_nu_nparticles.png')
+	hnpnubar_recoed = r.TH1F('np nubar recoed', 'np nubar recoed', 60, 0, 59)
+	hnpnubar_survivors = r.TH1F('np nubar survived', 'np nubar survived', 60, 0, 59)
+	seen = []
+	for part in studies['nubar']['data']:
+		if not (part.EventNumber in seen):
+			if part.EventNumber in studies['nubar']['listOfSurvivors']:
+				hnpnubar_survivors.Fill(part.NParticles)
+			if part.EventNumber in studies['nubar']['listOfRecoed']:
+				hnpnubar_recoed.Fill(part.NParticles)
+			seen.append(part.EventNumber)
+	#c1 = r.TCanvas(); c1.cd()
+	#hnpnubar.Draw()
+	#c1.Print('badEvents_nubar_nparticles.png')
+	hnpnu_recoed.SetLineColor(r.kBlack); hnpnu_recoed.SetMarkerColor(r.kBlack)
+	hnpnu_survivors.SetLineColor(r.kMagenta); hnpnu_survivors.SetMarkerColor(r.kMagenta)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_recoed.DrawNormalized('p'); hnpnu_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_recoed_and_survived.png')
+	c1.Print('nuAndNubar/nu_recoed_and_survived.pdf')
+	c1.Close()
+	hnpnubar_recoed.SetLineColor(r.kBlack); hnpnubar_recoed.SetMarkerColor(r.kBlack)
+	hnpnubar_survivors.SetLineColor(r.kMagenta); hnpnubar_survivors.SetMarkerColor(r.kMagenta)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnubar_recoed.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nubar_recoed_and_survived.png')
+	c1.Print('nuAndNubar/nubar_recoed_and_survived.pdf')
+	c1.Close()
+	hnpnu_recoed.SetLineColor(r.kBlue); hnpnu_recoed.SetMarkerColor(r.kBlue)
+	hnpnubar_recoed.SetLineColor(r.kGreen); hnpnubar_recoed.SetMarkerColor(r.kGreen)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_recoed.DrawNormalized('p'); hnpnubar_recoed.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_and_nubar_recoed.png')
+	c1.Print('nuAndNubar/nu_and_nubar_recoed.pdf')
+	c1.Close()
+	hnpnu_survivors.SetLineColor(r.kBlue); hnpnu_survivors.SetMarkerColor(r.kBlue)
+	hnpnubar_survivors.SetLineColor(r.kGreen); hnpnubar_survivors.SetMarkerColor(r.kGreen)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_survivors.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_and_nubar_survived.png')
+	c1.Print('nuAndNubar/nu_and_nubar_survived.pdf')
+	c1.Close()
+
+
+
+def do_my_studies():
+	results = {}
+	for s in studies:
+		print 'Reading %s...'%s
+		t = studies[s]['data']
+		ntot = t.GetEntriesFast()
+		nAfterCuts = [0]*10
+		nAfterCuts_weighted = [0.]*10
+		for event in t:
+			w = weight(event,s)
+			if fit_converged(event): 
+				nAfterCuts[0] += 1
+				nAfterCuts_weighted[0] += w
+				if chi2_cut(event,25):
+					nAfterCuts[1] +=1
+					nAfterCuts_weighted[1] +=w
+					#if chi2_cut(event,15):
+					if chi2_cut(event,5):
+						nAfterCuts[2]+=1
+						nAfterCuts_weighted[2]+=w
+						if doca_cut(event, 50.):
+							nAfterCuts[2+1] += 1
+							nAfterCuts_weighted[2+1] += w
+							if doca_cut(event, 30.):
+								nAfterCuts[3+1]+=1
+								nAfterCuts_weighted[3+1]+=w
+								if z_cut(event, s):
+									nAfterCuts[4+1] += 1
+									nAfterCuts_weighted[4+1] += w
+									if ip_cut(event,500.):
+										nAfterCuts[5+1]+=1
+										nAfterCuts_weighted[5+1]+=w
+										if ip_cut(event,250.):
+											nAfterCuts[6+1]+=1
+											nAfterCuts_weighted[6+1]+=w
+											if nothing_in_veto5(event):
+												nAfterCuts[7+1] += 1
+												nAfterCuts_weighted[7+1] += w
+												if nothing_in_liquidscint(event):
+													nAfterCuts[8+1] += 1
+													nAfterCuts_weighted[8+1] += w
+		print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot
+		print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot
+		results[s] = {}
+		results[s]['ntot']=studies[s]['ntot']
+		results[s]['numbers']=nAfterCuts
+		results[s]['numbersWithWeights']=nAfterCuts_weighted
+		addPlot(s+'-0cuts-IP', 'Impact parameter to target [cm]', t, 'IP0')
+		addPlot(s+'-0cuts-Mass', 'Reco inv. mass', t, 'Mass'              )
+		addPlot(s+'-0cuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             )
+		addPlot(s+'-0cuts-P','Reconstructed Momentum', t, 'P')
+		addPlot(s+'-0cuts-Pt', 'Reconstructed Pt', t, 'Pt')
+		addPlot(s+'-0cuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P')
+		addPlot(s+'-0cuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'             )
+		addPlot(s+'-0cuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'             )
+		addPlot(s+'-0cuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' )
+		addPlot(s+'-0cuts-DOCA', 'Daughter tracks DOCA [cm]', t, 'DOCA' )
+		addPlot(s+'-0cuts-dChi2', 'Daughter tracks #chi^{2}/n', t, 'DaughtersChi2/DaughtersNPoints' )
+		#cut = "DaughtersChi2<25*DaughtersNPoints"
+		cut = "DaughtersChi2<5*DaughtersNPoints"
+		addPlot(s+'-chi2cut-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2cut-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2cut-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2cut-P', 'Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2cut-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2cut-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2cut-dMinP', 'Min Tracks P', t, 'DaughterMinP'        , cut     )
+		addPlot(s+'-chi2cut-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'           , cut   )
+		addPlot(s+'-chi2cut-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+		addPlot(s+'-chi2cut-DOCA', 'Daughter tracks DOCA [cm]', t, 'DOCA'  , cut)
+		cut += " && DOCA < 25"
+		addPlot(s+'-chi2andDOCAcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2andDOCAcuts-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2andDOCAcuts-P','Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2andDOCAcuts-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2andDOCAcuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2andDOCAcuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2andDOCAcuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'        , cut      )
+		addPlot(s+'-chi2andDOCAcuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'     , cut         )
+		addPlot(s+'-chi2andDOCAcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+		#try:
+		#	minz = studies[s]['geo']['lidT1I_1']['z']['pos'] + studies[s]['geo']['lidT1I_1']['z']['dim'] + 500.
+		#	maxz = studies[s]['geo']['lidT6I_1']['z']['pos'] - studies[s]['geo']['lidT6I_1']['z']['dim']
+		#except:
+		minz = studies[s]['geo']['Veto_5']['z']['pos'] + studies[s]['geo']['Veto_5']['z']['dim']
+		maxz = studies[s]['geo']['Tr1_1']['z']['pos'] + studies[s]['geo']['Tr1_1']['z']['dim']
+		cut += " && vtxz > %s && vtxz < %s"%(minz, maxz)
+		addPlot(s+'-chi2andDOCAandZcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2andDOCAandZcuts-P','Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2andDOCAandZcuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'         , cut     )
+		addPlot(s+'-chi2andDOCAandZcuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'      , cut        )
+		addPlot(s+'-chi2andDOCAandZcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+	
+	
+	for s in studies.keys():#['pimu', 'mumunu', 'bg']:
+	    for list in results[s]:
+	        print s,list,results[s][list]
diff --git a/distribsForHNLandBG_byParticle.py b/distribsForHNLandBG_byParticle.py
new file mode 100755
index 0000000..dcca120
--- /dev/null
+++ b/distribsForHNLandBG_byParticle.py
@@ -0,0 +1,721 @@
+import sys
+import ROOT as r
+import numpy as np
+import funcsByBarbara as tools
+
+r.gROOT.SetBatch(r.kTRUE)
+r.gROOT.ProcessLine(".X lhcbStyle.C")
+plots = {}
+pdg = r.TDatabasePDG.Instance()
+
+def addPlot(key, x_axis, tree, string, cuts=None):
+	if not cuts:
+		tree.Draw(string+'>>h1')
+	else:
+		tree.Draw(string+'>>h1', cuts)
+	plots[key] = r.gDirectory.Get('h1')
+	plots[key].GetXaxis().SetTitle(x_axis)
+	plots[key].SetTitle(key)
+	plots[key].SetName(key)
+	c1 = r.TCanvas(); c1.cd()
+	if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+		c1.SetLogy()
+	plots[key].Draw()
+	c1.Print('accPlots/%s.pdf'%key)
+	c1.Print('accPlots/%s.png'%key)
+	key += "-WEIGHTED"
+	if 'bg' in key or ('nu' in key and key != 'mumunu'):
+		if not cuts:
+			tree.Draw(string+'>>h1', "NuWeight")
+		else:
+			cuts = "NuWeight*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+	elif 'cosmics' in key:
+		if not cuts:
+			tree.Draw(string+'>>h1',"HNLw*1.e6/200.")
+		else:
+			cuts = "HNLw*1.e6/200.*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+	else:
+		if not cuts:
+			tree.Draw(string+'>>h1',"HNLw")
+		else:
+			cuts = "HNLw*("+cuts+")"
+			tree.Draw(string+'>>h1', cuts)
+		plots[key] = r.gDirectory.Get('h1')
+		plots[key].GetXaxis().SetTitle(x_axis)
+		plots[key].SetTitle(key)
+		plots[key].SetName(key)
+		c1 = r.TCanvas(); c1.cd()
+		if "IP" in key or "Mass" in key or "DOCA" in key or "dChi2" in key:
+			c1.SetLogy()
+		plots[key].Draw()
+		c1.Print('accPlots/%s.pdf'%key)
+		c1.Print('accPlots/%s.png'%key)
+
+
+
+
+fpimu = r.TFile('../DATA/NewPIMU/ShipAna.root','read')#r.TFile('../DATA/MUPI/ShipAna.root','read')
+pimu = fpimu.Get('ShipAna')
+pimu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewPIMU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root')#'../DATA/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-2985556.root')
+
+fmumunu = r.TFile('../DATA/NewMUMUNU/ShipAna.root','read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read')
+mumunu = fmumunu.Get('ShipAna')
+mumunu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewMUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-mumunu.root')#'../DATA/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-TOTAL.root')
+
+ftiny = r.TFile('../DATA/SmallHNLs/ShipAna.root', 'read')
+tiny = ftiny.Get('ShipAna')
+tiny_geo = tools.searchForNodes3_xyz_dict('../DATA/SmallHNLs/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-011.root')
+
+#fbg = r.TFile('../DATA/KLONG/ShipAna.root','read')
+#bg = fbg.Get('ShipAna')
+#bg_geo = tools.searchForNodes2_xyz_dict('../DATA/KLONG/geofile_full.10.0.Genie-TGeant4_70cc30nc.root')
+
+fcosmics = r.TFile('../DATA/Cosmics/ShipAna.root','read')
+cosmics = fcosmics.Get('ShipAna')
+cosmics_geo = tools.searchForNodes2_xyz_dict('../DATA/Cosmics/geofile_full.Signal.10.0.Cosmics3.25.1-8.root')
+
+fnubar = r.TFile('../DATA/NewKLONG/ShipAna_nubar.root','read')#'../DATA/NewKLONG/ShipAna74-75_nubar.root','read')
+nubar = fnubar.Get('ShipAna')
+nubar_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nubar.root')#neutrino74-75.root')
+
+fnu = r.TFile('../DATA/NewKLONG/ShipAna_nu.root','read')#'../DATA/NewKLONG/ShipAna77_nu.root','read')
+nu = fnu.Get('ShipAna')
+nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nu.root')#neutrino77.root')
+
+
+# For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume
+
+studies = {
+	'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-4175},#1000},
+	'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-4125},#7000},
+	'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-4498},
+	#'bg': {'data':bg, 'geo':bg_geo, 'ntot':142856},#99999}#199998
+	'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'seen':[], 'ntot':cosmics.GetEntriesFast()},
+	'nubar': {'data': nubar, 'geo':nubar_geo, 'seen':[], 'ntot':99999*2},#99999*2},
+	'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':99999*2}#99999}
+}
+
+print 'Found %s entries in the cosmics file'%cosmics.GetEntriesFast()
+
+# Cuts:
+# - z_vertex in [vol[0]+5m, straw[0]]
+# - has track converged
+# - chi2/ndof < 10
+# - DOCA < 10 cm
+# - has no veto?
+
+
+
+class Cuts(object):
+    def __init__(self, cutConfigList, study): 
+        self.cutList = []      
+        try:        
+            for (index,cutConfig) in enumerate(cutConfigList):
+                self.cutList.append( self.Cut( self, index, cutConfig['name'], study, cutConfig['params'] ) )               
+        except BaseException,e:
+            print "cutConfigList badly set: %r" %e
+    def process(self, eventList):
+        for event in eventList:
+            for cut in self.cutList:
+                if  not cut.count(event): break
+        return self
+    def select(self, event, imax):
+    	for i in xrange(imax):
+    		if not self.cutList[i].selection(event, self.cutList[i].study, self.cutList[i].params):
+    			return False
+    	return True
+
+    class Cut(object):
+        def __init__(self, cutCollection, index, cutName, study, params):
+        	self.cutCollection = cutCollection
+        	self.index = index
+        	self.name = cutName
+        	self.params = params
+        	self.study = study
+        	self.total = 0 
+        	self.weight = 0
+        	self.eventNumbers = []
+        	self.duplicates = []
+        	self.selection = self.__get(cutName)
+        def __get(self, cutName):        
+            name = "cut_"+cutName
+            if hasattr(self,name):
+                try:
+                    return getattr(self,name)
+                except AttributeError,e:
+                    return None
+            else:
+                raise BaseException ( "bad cut name: %r" %cutName )
+                return None
+        def count(self, event):
+            try:
+                if (self.selection(event,self.study,self.params)):
+                    self.total += 1
+                    if (event.EventNumber not in self.eventNumbers):
+                        self.weight += self.compute_weight(event, self.study)
+                        self.eventNumbers.append(event.EventNumber)
+                    else:
+                    	self.duplicates.append(event.EventNumber)
+                    return True
+                else :
+                    return False
+            except BaseException,e:
+                print "Error performing %r: %r" %(self.name,e)
+        def cut_fit_converged(self, event, study, params):
+        	return fit_converged(event)
+        def cut_chi2(self, event, study, params):
+            return chi2_cut(event,params['chi2max'])
+        def cut_has_muons(self, event, study, params):
+        	for station in params['stations']:
+        		if not has_muon_station(event, study, station):
+        			return False
+        	return True
+        def cut_doca(self, event, study, params):
+        	return doca_cut(event, params['docamax'])
+        def cut_z(self, event, study, params):
+        	return z_cut(event, study)
+        def cut_normalization(self, event, study, params):
+        	return goodTruthVtx(event, study)
+        def cut_good_tracks(self, event, study, params):
+        	return goodTrack(event)
+        def cut_ip(self, event, study, params):
+        	return ip_cut(event, params['ipmax'])
+        def cut_vtx_in_ellipse(self, event, study, params):
+        	return vtxInAcceptance(event, study)
+        def cut_n_candidates(self, event, study, params):
+        	#print self.cutCollection.cutList[self.index-1].name
+        	if (event.EventNumber in self.cutCollection.cutList[self.index-1].duplicates):
+        		#if self.cutCollection.select(event, self.index):
+        		return nParticles_cut(event, params['nmax'])
+       		return True
+        def compute_weight(self, event, study):
+        	return weight(event, study)
+
+
+
+def fit_converged(event):
+	flags = event.DaughtersFitConverged
+	if np.product(flags):
+		return True
+	return False
+
+def chi2_cut(event, chi2max):
+	chi2 = event.DaughtersChi2
+	n = event.DaughtersNPoints
+	for i in xrange(len(n)):
+		if not n[i]:
+			print event.EventNumber
+			return False
+		if chi2[i]/n[i] > chi2max:
+			return False
+	return True
+
+def has_muon_station(event, study, station):
+	zIn = studies[study]['geo']['muondet%s_1'%(station-1)]['z']['pos'] - studies[study]['geo']['muondet%s_1'%(station-1)]['z']['dim']
+	zOut = studies[study]['geo']['muondet%s_1'%(station-1)]['z']['pos'] + studies[study]['geo']['muondet%s_1'%(station-1)]['z']['dim']
+	for z in event.muon_z:
+		if zIn <= z <= zOut:
+			return True
+	return False
+
+def doca_cut(event, doca_max):
+	if event.DOCA < doca_max:
+		return True
+	return False
+
+def z_cut(event, study):
+	z = event.vtxz
+	#try:
+	#	minz = studies[study]['geo']['lidT1I_1']['z']['pos'] + studies[study]['geo']['lidT1I_1']['z']['dim'] + 500.
+	#	maxz = studies[study]['geo']['lidT6I_1']['z']['pos'] - studies[study]['geo']['lidT6I_1']['z']['dim']
+	#except:
+	minz = studies[study]['geo']['Veto_5']['z']['pos'] + studies[study]['geo']['Veto_5']['z']['dim'] + 0.01
+	maxz = studies[study]['geo']['Tr1_1']['z']['pos'] - studies[study]['geo']['Tr1_1']['z']['dim'] - 0.01
+	if minz < z < maxz:
+		return True
+	return False
+
+def goodTruthVtx(event, study):
+	if study in ['pimu', 'mumunu', 'tiny']:
+		if not event.BadTruthVtx:
+			return True
+	return True
+
+def nParticles_cut(event, nmax):
+	if event.NParticles > nmax:
+		return False
+	return True
+
+# This checks that:
+# - there is at least one tracking station before and one after the magnet
+# - hits in these tracking stations are within a 5x10m ellipse
+def goodTrack(event):
+	if event.DaughtersAlwaysIn:
+		return True
+	return False
+
+def ip_cut(event, maxip):
+	if event.IP0 < maxip:
+		return True
+	return False
+
+def nothing_in_veto5(event):
+	if len(event.veto5_x) > 0:
+		return False
+	return True
+
+def nothing_in_liquidscint(event):
+	if len(event.liquidscint_x) > 0:
+		return False
+	return True
+
+def something_in_muon(event):
+	if len(event.muon_x) > 0:
+		return True
+	return False
+
+def something_in_ecal(event):
+	if len(event.ecal_x) > 0:
+		return True
+	return False
+
+def vtxInAcceptance(event, study):
+	x = event.vtxx
+	y = event.vtxy
+	#if 'cosmics' in study or 'bg' in study:
+	#	print study, x, y
+	return tools.inEllipse(x,y,500./2,1000./2)
+
+def weight(event, study):
+	#if event.EventNumber not in studies[study]['seen']:
+		#studies[study]['seen'].append(event.EventNumber)
+	if ('bg' in study) or (study == 'nu' or study == 'nubar'):
+		w = event.NuWeight
+	elif 'cosmics' in study:
+		w = event.HNLw * 1.e6 / 200.
+	elif study in ['pimu', 'mumunu', 'tiny']:
+		w = event.HNLw / studies[study]['ntot']
+	else:
+		print 'WEIGHTING ERROR: ', study
+		sys.exit()
+	#else:
+	#	w = 0.
+		#if study in ['pimu', 'mumunu', 'tiny']:
+		#	print study, event.EventNumber, event.NParticles
+	return w
+
+def computeResidualWeight(filename=None):
+	if not filename:
+		filename = "forElena_nu_notVetoed.txt"
+	fScan = file(filename)
+	lines = fScan.readlines()
+	sample = nu
+	s = 'nu'
+	if ('bar' in filename) or ('anti' in filename):
+		sample = nubar
+		s = 'nubar'
+	lines = lines[3:-1]
+	events = []
+	for line in lines:
+	    events.append(int(line.rsplit("*")[1].replace(" ","")))
+	count_reco = 0
+	count = 0
+	weighted_count = 0.
+	print "unvetoed: ", len(events)
+	print "total Reconstructed: ", sample.GetEntriesFast()
+	for event in sample:
+		if event.EventNumber in events:
+			count_reco +=1
+			if (fit_converged(event) and chi2_cut(event, 5.) and z_cut(event, s) and vtxInAcceptance(event, s)
+				and goodTrack(event) and doca_cut(event, 30) and has_muon_station(event, s, 1)
+				and has_muon_station(event, s, 2) and ip_cut(event, 250.)):
+			#if True:#fit_converged(event): #and doca_cut(event, 50.) and ip_cut(event, 250.) and z_cut(event, 'nu') and vtxInAcceptance(event, 'nu'):
+				#print event.EventNumber, 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+				print event.EventNumber
+				print 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+				#print event.DOCA, event.IP0, event.DaughtersChi2[0]/event.DaughtersNPoints[0], event.DaughtersChi2[1]/event.DaughtersNPoints[1], event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+				print "tracks: ", event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+				print "vtx: ", event.vtxx, event.vtxy, event.vtxz
+				print 'DOCA, IP: ', event.DOCA, event.IP0
+				#print tools.interactionElement(event.DaughtersTruthProdX[0], event.DaughtersTruthProdY[0], event.DaughtersTruthProdZ[0], nu_geo), tools.interactionElement(event.DaughtersTruthProdX[1], event.DaughtersTruthProdY[1], event.DaughtersTruthProdZ[1],nu_geo)
+				count +=1
+				weighted_count += event.NuWeight
+	print "Reconstructed: %s"%count_reco
+	return count, weighted_count
+
+def latex_float(f):
+    float_str = "{0:.3g}".format(f)
+    if "e" in float_str:
+        base, exponent = float_str.split("e")
+        return "${0} \\times 10^{{{1}}}$".format(base, int(exponent))
+    else:
+        return float_str
+
+cutConfigList = [   {'name': 'normalization', 'params':{}},
+					{'name': 'fit_converged', 'params':{}},
+					{'name': 'chi2',          'params':{'chi2max':5.}},
+					{'name': 'z',             'params':{}},
+					{'name': 'vtx_in_ellipse','params':{}},
+					{'name': 'good_tracks',   'params':{}},
+					{'name': 'doca', 		  'params':{'docamax':50.}},
+					{'name': 'ip', 		      'params':{'ipmax':500.}},
+					{'name': 'has_muons',     'params':{'stations':[1,2]}},
+					{'name': 'doca', 		  'params':{'docamax':30.}},
+					{'name': 'ip', 		      'params':{'ipmax':250.}},
+					{'name': 'n_candidates',  'params':{'nmax':1}},
+]
+
+def numbers_v2():
+	mycuts, results = {}, {}
+	for study in studies:
+		print 'Reading %s...'%study
+		mycuts[study] = Cuts(cutConfigList, study).process(studies[study]['data'])
+		for cut in mycuts[study].cutList:
+			print '%s\t%s\t%s\t%s'%(cut.name, cut.total, len(cut.eventNumbers), cut.weight)
+		results[study] = {}
+		results[study]['numbers']=[ cut.total for cut in mycuts[study].cutList[1:] ]
+		results[study]['numbersWithWeights']=[ cut.weight for cut in mycuts[study].cutList[1:] ]
+		print
+		print 'Events with more than one particle:'
+		for entry in mycuts[study].cutList[-2].duplicates: #[ x for x in mycuts[study].cutList[-1].eventNumbers if mycuts[study].cutList[-1].eventNumbers.count(x) > 1 ]:
+			for tr in studies[study]['data']:
+				if tr.EventNumber == entry:
+					if mycuts[study].select(tr, len(mycuts[study].cutList)-1):
+						print tr.EventNumber, tr.vtxz, tr.IP0, [pdg.GetParticle(code).GetName() for code in tr.DaughtersTruthPDG]
+		print
+		print
+	return results
+
+def makeTex(results):
+	print
+	print 'Background:'
+	print '\t\\begin{table}'
+	print '\t\\centering Backgrounds: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['cosmics']['numbersWithWeights'][i]), ' & ', '%s'%results['cosmics']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print
+	print 'Signal:'
+	print '\t\\begin{table}'
+	print '\t\\centering Signal: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'IP $<$ 5 m', 'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 2.5 m', '$N_{candidates} = 1$']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['tiny']['numbersWithWeights'][i]), ' & ', '%s'%results['tiny']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print	
+
+
+
+def numbers_for_TP():
+	results = {}
+	for s in ['pimu', 'mumunu', 'tiny', 'nu', 'nubar', 'cosmics']:#studies:
+		print 'Reading %s...'%s
+		t = studies[s]['data']
+		ntot = t.GetEntriesFast()
+		nAfterCuts = [0]*10
+		nAfterCuts_weighted = [0.]*10
+		initialEventNums = []
+		eventNums = []
+		for event in t:
+			w = weight(event,s) # e' importante fare questa cosa a monte, x come e' definita la funzione!
+			if fit_converged(event) and goodTruthVtx(event,s):
+				nAfterCuts[0] +=1
+				nAfterCuts_weighted[0] += w
+				initialEventNums.append(int(event.EventNumber))
+				if chi2_cut(event, 5.):
+					nAfterCuts[1]+=1
+					nAfterCuts_weighted[1]+=w
+					if z_cut(event, s):#chi2_cut(event, 5):
+						nAfterCuts[1+1] += 1
+						nAfterCuts_weighted[1+1] += w
+						if vtxInAcceptance(event, s):#goodTrack(event):
+							nAfterCuts[2+1] +=1
+							nAfterCuts_weighted[2+1] +=w
+							#if 'cosmics' in s:
+							#	print 'mothers: ', event.DaughtersTruthMotherPDG[0], event.DaughtersTruthMotherPDG[1]
+							#	print "tracks: ", event.DaughtersTruthPDG[0], event.DaughtersTruthPDG[1] 
+							#	print "p candidate: ", event.P, "\t IP: ", event.IP0
+							#	print "pt candidate: ", event.Pt
+							if goodTrack(event):#doca_cut(event, 50):
+								nAfterCuts[3+1] +=1
+								nAfterCuts_weighted[3+1] +=w
+								if doca_cut(event, 50):#vtxInAcceptance(event, s):
+									nAfterCuts[4+1]+=1
+									nAfterCuts_weighted[4+1]+=w
+									if has_muon_station(event, s, 1) and has_muon_station(event, s, 2):#ip_cut(event,500):#something_in_muon(event):
+										nAfterCuts[5+1] +=1
+										nAfterCuts_weighted[5+1] +=w
+										if doca_cut(event, 30.):#something_in_ecal(event):
+											nAfterCuts[6+1] +=1
+											nAfterCuts_weighted[6+1] +=w
+											if ip_cut(event, 500.):
+												nAfterCuts[7+1] +=1
+												nAfterCuts_weighted[7+1] +=w
+												if ip_cut(event, 250.):
+													nAfterCuts[8+1] +=1
+													nAfterCuts_weighted[8+1] +=w
+													eventNums.append(int(event.EventNumber))
+													studies[s]['seen'].append(event.EventNumber)
+
+		print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot
+		print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot
+		results[s] = {}
+		results[s]['ntot']=studies[s]['ntot']
+		results[s]['numbers']=nAfterCuts
+		results[s]['numbersWithWeights']=nAfterCuts_weighted
+		results[s]['listOfSurvivors']=list(set(eventNums))
+		studies[s]['listOfSurvivors']=list(set(eventNums))
+		results[s]['listOfRecoed']=list(set(initialEventNums))
+		studies[s]['listOfRecoed']=list(set(initialEventNums))
+	print
+	for s in ['pimu', 'mumunu', 'tiny']:#studies.keys():#['pimu', 'mumunu', 'bg']:
+	    for lst in results[s]:
+		if not ('Survivors' in lst) and not ('Recoed' in lst):
+		        print s,list,results[s][lst]
+	print
+	print 'Events with surviving nu: %s \t with nubar: %s'%(len(results['nu']['listOfSurvivors']), len(results['nubar']['listOfSurvivors']))
+	print 
+	print 'Background:'
+	print '\t\\begin{table}'
+	print '\t\\centering Backgrounds: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & Weighted $\\nu$ & $\\nu$ entries & Weighted $\\bar{\\nu}$ & $\\bar{\\nu}$ entries & Weighted cosmics & cosmics entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 5 m', 'IP $<$ 2.5 m']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['nu']['numbersWithWeights'][i]), ' & ', '%s'%results['nu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['nubar']['numbersWithWeights'][i]), ' & ', '%s'%results['nubar']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['cosmics']['numbersWithWeights'][i]), ' & ', '%s'%results['cosmics']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print 'Signal:'
+	print '\t\\begin{table}'
+	print '\t\\centering Signal: \\tiny'
+	print '\t\t\\begin{tabular}{|l|c|c|c|c|c|c|}'
+	print '\t\t\t\\hline'
+	print '\t\t\tCut & $\\mathcal{A}(\\pi\\mu)$ & $\\pi\\mu$ entries & $\\mathcal{A}(\\mu\\mu\\nu)$ & $\\mu\\mu\\nu$ entries & $\\mathcal{A}(ee\\nu)$ & $ee\\nu$ entries \\\\'
+	print '\t\t\t\\hline'
+	cuts = ['Track fit converged', '$\chi^2/n < 5$', '$z_{in}+5~\\text{m} < z_{vtx} < z_{out}$', '$z_{vtx} \in $ ellipse', 'Tracks $\in$ ellipse', 'DOCA $<$ 50 cm',
+			'Muon1 \& Muon2', 'DOCA $<$ 30 cm', 'IP $<$ 5 m', 'IP $<$ 2.5 m']
+	for i in xrange(len(results['nu']['numbers'])):
+		print   '\t\t\t', cuts[i], ' & ', '%s'%latex_float(results['pimu']['numbersWithWeights'][i]), ' & ', '%s'%results['pimu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['mumunu']['numbersWithWeights'][i]), ' & ', '%s'%results['mumunu']['numbers'][i], ' & ',\
+				'%s'%latex_float(results['tiny']['numbersWithWeights'][i]), ' & ', '%s'%results['tiny']['numbers'][i], '\\\\ '
+		print '\t\t\t\\hline'
+	print '\t\t\\end{tabular}'
+	print '\t\\end{table}'
+	print
+	out = {'pimu':[], 'mumunu':[], 'tiny':[]}
+	for s in ['pimu', 'mumunu', 'tiny']:
+		for event in studies[s]['data']:
+			if event.EventNumber in studies[s]['listOfSurvivors']:
+				if event.NParticles > 1:
+					out[s].append(event.EventNumber)
+					#print s, event.EventNumber, event.NParticles, r.gGeoManager.FindNode(event.vtxx, event.vtxy, event.vtxz).GetName()
+		print s, list(set(out[s]))
+
+		
+def badEvents():
+	numbers_for_TP()
+	hnpnu_recoed = r.TH1F('np nu recoed', 'np nu recoed', 60, 0, 59)
+	hnpnu_survivors = r.TH1F('np nu survived', 'np nu survived', 60, 0, 59)
+	seen = []
+	for part in studies['nu']['data']:
+		if not (part.EventNumber in seen):
+			if part.EventNumber in studies['nu']['listOfSurvivors']:
+				hnpnu_survivors.Fill(part.NParticles)
+			if part.EventNumber in studies['nu']['listOfRecoed']:
+				hnpnu_recoed.Fill(part.NParticles)
+			seen.append(part.EventNumber)
+	#c1 = r.TCanvas(); c1.cd()
+	#hnpnu.Draw()
+	#c1.Print('badEvents_nu_nparticles.png')
+	hnpnubar_recoed = r.TH1F('np nubar recoed', 'np nubar recoed', 60, 0, 59)
+	hnpnubar_survivors = r.TH1F('np nubar survived', 'np nubar survived', 60, 0, 59)
+	seen = []
+	for part in studies['nubar']['data']:
+		if not (part.EventNumber in seen):
+			if part.EventNumber in studies['nubar']['listOfSurvivors']:
+				hnpnubar_survivors.Fill(part.NParticles)
+			if part.EventNumber in studies['nubar']['listOfRecoed']:
+				hnpnubar_recoed.Fill(part.NParticles)
+			seen.append(part.EventNumber)
+	#c1 = r.TCanvas(); c1.cd()
+	#hnpnubar.Draw()
+	#c1.Print('badEvents_nubar_nparticles.png')
+	hnpnu_recoed.SetLineColor(r.kBlack); hnpnu_recoed.SetMarkerColor(r.kBlack)
+	hnpnu_survivors.SetLineColor(r.kMagenta); hnpnu_survivors.SetMarkerColor(r.kMagenta)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_recoed.DrawNormalized('p'); hnpnu_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_recoed_and_survived.png')
+	c1.Print('nuAndNubar/nu_recoed_and_survived.pdf')
+	c1.Close()
+	hnpnubar_recoed.SetLineColor(r.kBlack); hnpnubar_recoed.SetMarkerColor(r.kBlack)
+	hnpnubar_survivors.SetLineColor(r.kMagenta); hnpnubar_survivors.SetMarkerColor(r.kMagenta)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnubar_recoed.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nubar_recoed_and_survived.png')
+	c1.Print('nuAndNubar/nubar_recoed_and_survived.pdf')
+	c1.Close()
+	hnpnu_recoed.SetLineColor(r.kBlue); hnpnu_recoed.SetMarkerColor(r.kBlue)
+	hnpnubar_recoed.SetLineColor(r.kGreen); hnpnubar_recoed.SetMarkerColor(r.kGreen)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_recoed.DrawNormalized('p'); hnpnubar_recoed.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_and_nubar_recoed.png')
+	c1.Print('nuAndNubar/nu_and_nubar_recoed.pdf')
+	c1.Close()
+	hnpnu_survivors.SetLineColor(r.kBlue); hnpnu_survivors.SetMarkerColor(r.kBlue)
+	hnpnubar_survivors.SetLineColor(r.kGreen); hnpnubar_survivors.SetMarkerColor(r.kGreen)
+	c1 = r.TCanvas(); c1.cd()
+	hnpnu_survivors.DrawNormalized('p'); hnpnubar_survivors.DrawNormalized('p same')
+	c1.Print('nuAndNubar/nu_and_nubar_survived.png')
+	c1.Print('nuAndNubar/nu_and_nubar_survived.pdf')
+	c1.Close()
+
+
+
+def do_my_studies():
+	results = {}
+	for s in studies:
+		print 'Reading %s...'%s
+		t = studies[s]['data']
+		ntot = t.GetEntriesFast()
+		nAfterCuts = [0]*10
+		nAfterCuts_weighted = [0.]*10
+		for event in t:
+			w = weight(event,s)
+			if fit_converged(event): 
+				nAfterCuts[0] += 1
+				nAfterCuts_weighted[0] += w
+				if chi2_cut(event,25):
+					nAfterCuts[1] +=1
+					nAfterCuts_weighted[1] +=w
+					#if chi2_cut(event,15):
+					if chi2_cut(event,5):
+						nAfterCuts[2]+=1
+						nAfterCuts_weighted[2]+=w
+						if doca_cut(event, 50.):
+							nAfterCuts[2+1] += 1
+							nAfterCuts_weighted[2+1] += w
+							if doca_cut(event, 30.):
+								nAfterCuts[3+1]+=1
+								nAfterCuts_weighted[3+1]+=w
+								if z_cut(event, s):
+									nAfterCuts[4+1] += 1
+									nAfterCuts_weighted[4+1] += w
+									if ip_cut(event,500.):
+										nAfterCuts[5+1]+=1
+										nAfterCuts_weighted[5+1]+=w
+										if ip_cut(event,250.):
+											nAfterCuts[6+1]+=1
+											nAfterCuts_weighted[6+1]+=w
+											if nothing_in_veto5(event):
+												nAfterCuts[7+1] += 1
+												nAfterCuts_weighted[7+1] += w
+												if nothing_in_liquidscint(event):
+													nAfterCuts[8+1] += 1
+													nAfterCuts_weighted[8+1] += w
+		print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot
+		print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot
+		results[s] = {}
+		results[s]['ntot']=studies[s]['ntot']
+		results[s]['numbers']=nAfterCuts
+		results[s]['numbersWithWeights']=nAfterCuts_weighted
+		addPlot(s+'-0cuts-IP', 'Impact parameter to target [cm]', t, 'IP0')
+		addPlot(s+'-0cuts-Mass', 'Reco inv. mass', t, 'Mass'              )
+		addPlot(s+'-0cuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             )
+		addPlot(s+'-0cuts-P','Reconstructed Momentum', t, 'P')
+		addPlot(s+'-0cuts-Pt', 'Reconstructed Pt', t, 'Pt')
+		addPlot(s+'-0cuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P')
+		addPlot(s+'-0cuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'             )
+		addPlot(s+'-0cuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'             )
+		addPlot(s+'-0cuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' )
+		addPlot(s+'-0cuts-DOCA', 'Daughter tracks DOCA [cm]', t, 'DOCA' )
+		addPlot(s+'-0cuts-dChi2', 'Daughter tracks #chi^{2}/n', t, 'DaughtersChi2/DaughtersNPoints' )
+		#cut = "DaughtersChi2<25*DaughtersNPoints"
+		cut = "DaughtersChi2<5*DaughtersNPoints"
+		addPlot(s+'-chi2cut-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2cut-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2cut-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2cut-P', 'Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2cut-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2cut-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2cut-dMinP', 'Min Tracks P', t, 'DaughterMinP'        , cut     )
+		addPlot(s+'-chi2cut-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'           , cut   )
+		addPlot(s+'-chi2cut-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+		addPlot(s+'-chi2cut-DOCA', 'Daughter tracks DOCA [cm]', t, 'DOCA'  , cut)
+		cut += " && DOCA < 25"
+		addPlot(s+'-chi2andDOCAcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2andDOCAcuts-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2andDOCAcuts-P','Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2andDOCAcuts-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2andDOCAcuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2andDOCAcuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2andDOCAcuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'        , cut      )
+		addPlot(s+'-chi2andDOCAcuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'     , cut         )
+		addPlot(s+'-chi2andDOCAcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+		#try:
+		#	minz = studies[s]['geo']['lidT1I_1']['z']['pos'] + studies[s]['geo']['lidT1I_1']['z']['dim'] + 500.
+		#	maxz = studies[s]['geo']['lidT6I_1']['z']['pos'] - studies[s]['geo']['lidT6I_1']['z']['dim']
+		#except:
+		minz = studies[s]['geo']['Veto_5']['z']['pos'] + studies[s]['geo']['Veto_5']['z']['dim']
+		maxz = studies[s]['geo']['Tr1_1']['z']['pos'] + studies[s]['geo']['Tr1_1']['z']['dim']
+		cut += " && vtxz > %s && vtxz < %s"%(minz, maxz)
+		addPlot(s+'-chi2andDOCAandZcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-Mass', 'Reco inv. mass', t, 'Mass'              , cut)
+		addPlot(s+'-chi2andDOCAandZcuts-P','Reconstructed Momentum', t, 'P', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-Pt','Reconstructed Pt', t, 'Pt', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-PtOverP', 'Reconstructed Pt/P', t, 'Pt/P', cut)
+		addPlot(s+'-chi2andDOCAandZcuts-dPt', 'Tracks Pt', t, 'DaughtersPt'             , cut)
+		addPlot(s+'-chi2andDOCAandZcuts-dMinP', 'Min Tracks P', t, 'DaughterMinP'         , cut     )
+		addPlot(s+'-chi2andDOCAandZcuts-dMinPt', 'Min Tracks Pt', t, 'DaughterMinPt'      , cut        )
+		addPlot(s+'-chi2andDOCAandZcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut)
+	
+	
+	for s in studies.keys():#['pimu', 'mumunu', 'bg']:
+	    for list in results[s]:
+	        print s,list,results[s][list]