Newer
Older
FairShipTools / distribsForHNLandBG.py
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]