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]
Ubuntu