import ROOT as r import numpy as np import funcsByBarbara as tools r.gROOT.SetBatch(r.kTRUE) r.gROOT.ProcessLine(".X lhcbStyle.C") plots = {} 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-009-eenu.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/ShipAna74-75_nubar.root','read') nubar = fnubar.Get('ShipAna') nubar_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_neutrino74-75.root') fnu = r.TFile('../DATA/NewKLONG/ShipAna77_nu.root','read') nu = fnu.Get('ShipAna') nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_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, 'ntot':20000-4175},#1000}, 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'ntot':20000-4125},#7000}, 'tiny': {'data': tiny, 'geo':tiny_geo, 'ntot':20000-4452}, 'bg': {'data':bg, 'geo':bg_geo, 'ntot':142856},#99999}#199998 'cosmics': {'data':cosmics, 'geo':cosmics_geo, 'ntot':cosmics.GetEntriesFast()}, 'nubar': {'data': nubar, 'geo':nubar_geo, 'ntot':99999*2}, 'nu': {'data': nu, 'geo':nu_geo, 'ntot':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? 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'] maxz = studies[study]['geo']['Tr1_1']['z']['pos'] + studies[study]['geo']['Tr1_1']['z']['dim'] 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 # 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): w = event.HNLw / studies[study]['ntot'] if ('bg' in study) or (study == 'nu' or 'study' == 'nubar'): w = event.NuWeight if 'cosmics' in study: w = event.HNLw * 1.e6 / 200. return w def computeResidualWeight(filename): if not filename: filename = "forElena_nu_notVetoed.txt" fScan = file(filename) lines = fScan.readlines() 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: ", nu.GetEntriesFast() #for entry in events: # bg.GetEntry(entry) # particles = bg.Particles # for p in particles: # count+=1 # weighted_count += bg.NuWeight for event in nu: if event.EventNumber in events: print event.EventNumber if True:#fit_converged(event): count_reco +=1 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 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 for event in t: w = weight(event,s) if fit_converged(event) and goodTruthVtx(event,s): nAfterCuts[0] +=1 nAfterCuts_weighted[0] += w 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 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 print for s in ['pimu', 'mumunu', 'tiny']:#studies.keys():#['pimu', 'mumunu', 'bg']: for list in results[s]: print s,list,results[s][list] print 'Background:' print '\t\\begin{table}' print '\t\\centering\\footnotesize' 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\\footnotesize' 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}' 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]