diff --git a/distribsForHNLandBG_byEvent.py b/distribsForHNLandBG_byEvent.py index 92991a2..8f0ea20 100755 --- a/distribsForHNLandBG_byEvent.py +++ b/distribsForHNLandBG_byEvent.py @@ -76,15 +76,15 @@ -fpimu = r.TFile('../DATA/NewPIMU/ShipAna_Barbara_wholeStat.root','read')#r.TFile('../DATA/MUPI/ShipAna.root','read') +fpimu = r.TFile('../DATA/NewPIMU/ShipAna_newGen.root','read')#r.TFile('../DATA/MUPI/ShipAna.root','read') pimu = fpimu.Get('t') 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_Barbara_wholeStat.root','read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read') +fmumunu = r.TFile('../DATA/NewMUMUNU/ShipAna_newGen.root','read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read') mumunu = fmumunu.Get('t') 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_Barbara_wholeStat.root', 'read') +ftiny = r.TFile('../DATA/SmallHNLs/ShipAna_newGen.root', 'read') tiny = ftiny.Get('t') tiny_geo = tools.searchForNodes3_xyz_dict('../DATA/SmallHNLs/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-011.root') @@ -100,6 +100,10 @@ #nu = fnu.Get('ShipAna') #nu_geo = tools.searchForNodes3_xyz_dict('../DATA/NewKLONG/geofile_full.10.0.Genie-TGeant4_nu.root')#neutrino77.root') +# Neutrinos from yandex +fnu = r.TFile('../DATA/nuFromYandex/newGen_15340000_nuBg.root') +nu = fnu.Get('t') +nu_geo = tools.searchForNodes3_xyz_dict('../DATA/nuFromYandex/geofile_full.10.0.Genie-TGeant4.root')#neutrino77.root') # For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume @@ -107,6 +111,7 @@ 'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-136},#1000}, 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-289},#7000}, 'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-66}, + 'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':1.534e7}, #'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} @@ -121,29 +126,46 @@ # - DOCA < 10 cm # - has no veto? +flux_nu = (1.09e+12 * 2.e+20)/ 5.e+13 +flux_anti = 6.98e+11 * 2.e+20 / 5.e+13 +entries_nu = 1.534e7#200000-2 +entries_anti = 200000-2 +weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu) +weight_anti ="(weight * %s / %s)"%(flux_anti,entries_anti) + class MyParticle(object): """ Simulating an ntuple filled by particle from one filled by event :) """ - def __init__(self, event, index): - self.EventNumber = event.event - self.DaughtersFitConverged = [event.DaughtersFitConverged[index], event.DaughtersFitConverged[index]] - chi2red = [event.MaxDaughtersRedChi2[index], event.MaxDaughtersRedChi2[index]] - ndf = [event.MinDaughtersNdf[index], event.MinDaughtersNdf[index]] + def __init__(self, event, index, noB='', eff=''): + self.EventNumber = event.__getattr__('event') + self.DaughtersFitConverged = [event.__getattr__('DaughtersFitConverged')[index], event.__getattr__('DaughtersFitConverged')[index]] + chi2red = [event.__getattr__('%sMaxDaughtersRedChi2'%noB)[index], event.__getattr__('%sMaxDaughtersRedChi2'%noB)[index]] + ndf = [event.__getattr__('%sMinDaughtersNdf'%noB)[index], event.__getattr__('%sMinDaughtersNdf'%noB)[index]] self.DaughtersChi2 = [chi2red[i]*ndf[i] for i in xrange(len(ndf))] self.DaughtersNPoints = ndf - self.DOCA = event.DOCA[index] - self.vtxz = event.vtxz[index] - self.vtxy = event.vtxy[index] - self.vtxx = event.vtxx[index] - self.IP0 = event.IP0[index] - self.BadTruthVtx = event.BadTruthVtx - self.NParticles = event.nRecoed - self.DaughtersAlwaysIn = event.DaughtersAlwaysIn[index] - self.Has1Muon1 = event.Has1Muon1[index] - self.Has2Muon1 = event.Has2Muon1[index] - self.Has1Muon2 = event.Has1Muon2[index] - self.Has2Muon2 = event.Has2Muon2[index] - self.HasEcal = event.HasEcal[index] - self.HNLw = event.weight + self.DOCA = event.__getattr__('%sDOCA'%noB)[index] + self.vtxz = event.__getattr__('%svtxz'%noB)[index] + try: + self.vtxy = event.__getattr__('%svtxy'%noB)[index] + self.vtxx = event.__getattr__('%svtxx'%noB)[index] + except: + self.vtxy = r.TMath.Sqrt(event.__getattr__('%svtxySqr'%noB)[index]) + self.vtxx = r.TMath.Sqrt(event.__getattr__('%svtxxSqr'%noB)[index]) + self.IP0 = event.__getattr__('%sIP0'%noB)[index] + self.BadTruthVtx = event.__getattr__('BadTruthVtx') + self.NParticles = event.__getattr__('nRecoed') + self.DaughtersAlwaysIn = event.__getattr__('DaughtersAlwaysIn')[index] + self.Has1Muon1 = event.__getattr__('Has1Muon1')[index] + self.Has2Muon1 = event.__getattr__('Has2Muon1')[index] + self.Has1Muon2 = event.__getattr__('Has1Muon2')[index] + self.Has2Muon2 = event.__getattr__('Has2Muon2')[index] + self.HasEcal = event.__getattr__('HasEcal')[index] + self.HNLw = event.__getattr__('weight') + self.NuWeight = event.__getattr__('weight')*flux_nu/entries_nu + self.numberOfHitLSSegments = event.__getattr__('numberOfHitLSSegments') + self.strawVetoAny = event.__getattr__('strawVetoAny%s'%eff) + self.upstreamVetoAny = event.__getattr__('upstreamVetoAny%s'%eff) + self.RPCany = event.__getattr__('RPCany%s'%eff) + @@ -156,12 +178,12 @@ 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): + def process(self, eventList, noB, eff=''): if self.byEvent: for event in eventList: for index in xrange(event.nRecoed): for cut in self.cutList: - if not cut.count(MyParticle(event, index)): break + if not cut.count(MyParticle(event, index, noB, eff)): break else: for event in eventList: for cut in self.cutList: @@ -242,10 +264,17 @@ return has_ecalDeposit(event) def cut_has_n_muons_in(self, event, study, params): return has_n_muons_in_station(event, params['n'], params['station']) + def cut_not_vetoed(self, event, study, params): + return was_not_vetoed(event, study) def compute_weight(self, event, study): return weight(event, study) +def was_not_vetoed(event, study): + if study in ['pimu', 'mumunu', 'tiny']: + return True + flag = (event.numberOfHitLSSegments == 0) and (event.strawVetoAny==0) and (event.upstreamVetoAny==0) and (event.RPCany==0) + return flag def fit_converged(event): flags = event.DaughtersFitConverged @@ -436,51 +465,76 @@ class cutsWithDraw(object): def __init__(self): self.noB = '' - cutEllipse = "( (vtxx*vtxx)/(250.*250.) + (vtxy*vtxy)/(500.*500.) ) < 1." + self.opt = '' # '_eff' # 'Ethr' def setNoB(self, noB): # to ignore or take into account Katerina's modifications self.noB = noB - def doca(self, th, noB=self.noB): + def doca(self, th, noB=None): + if not noB: + noB = self.noB return " Sum$(%sDOCA < %s) "%(noB, th) - def ip(self, th, noB=self.noB): + def ip(self, th, noB=None): + if not noB: + noB = self.noB return " Sum$(%sIP0 < %s) "%(noB, th) - def z(self, zmin, zmax, noB=self.noB): + def z(self, zmin, zmax, noB=None): + if not noB: + noB = self.noB return " Sum$(%svtxz > %s) && Sum$(%svtxz < %s) "%(noB, zmin, noB, zmax) - def ellipse(self, a, b): - return " Sum$(( (vtxx*vtxx)/(%s*%s) + (vtxy*vtxy)/(%s*%s) ) < 1.) "%(a,a,b,b) + def ellipse(self, a, b, noB=None): + if not noB: + noB = self.noB + return " Sum$(( (%svtxx*%svtxx)/(%s*%s) + (%svtxy*%svtxy)/(%s*%s) ) < 1.) "%(noB,noB, a,a, noB,noB, b,b) + def ellipseSq(self, a, b, noB=None): + if not noB: + noB = self.noB + return " Sum$(( (%svtxxSqr)/(%s*%s) + (%svtxySqr)/(%s*%s) ) < 1.) "%(noB, a,a, noB, b,b) def converged(self): return " Sum$(DaughtersFitConverged==1) " def goodtracks(self): return " Sum$(DaughtersAlwaysIn==1) " - def redChi2(self, th, noB=self.noB): + def redChi2(self, th, noB=None): + if not noB: + noB = self.noB return " Sum$(%sMaxDaughtersRedChi2<%s) "%(noB, th) - def ndf(self, th, noB=self.noB): + def ndf(self, th, noB=None): + if not noB: + noB = self.noB return " Sum$(%sMinDaughtersNdf>%s) "%(noB, th) def muon1(self, num): return " Sum$(Has%sMuon1==1) "%num def muon2(self, num): return " Sum$(Has%sMuon2==1) "%num + def notVetoed(self, opt=None): + if not opt: + opt = self.opt + return " Sum$(numberOfHitLSSegments == 0) && Sum$(strawVetoAny%s==0) && Sum$(upstreamVetoAny%s==0) && Sum$(RPCany%s==0) "%(opt, opt, opt) -def countWithDraw(): - t = studies['pimu']['data'] - geo = studies['pimu']['geo'] - ntot = studies['pimu']['ntot'] + +def countWithDraw(study): + t = studies[study]['data'] + geo = studies[study]['geo'] + ntot = studies[study]['ntot'] tc = cutsWithDraw() + tc.setNoB('NoB_') zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim'] - zmin = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] - cuts = [ tc.converged(), tc.redChi2(5.), tc.z(zmin, zmax), tc.ellipse(250., 500.), tc.goodtracks(), - tc.doca(30.), tc.ip(250.), tc.muon1(1), tc.muon2(1) ] - cutnames = [ "converged", "chi2/n<5", "z", "ellipse", "good tracks", "doca<30", "ip<250", "muon1", "muon2" ] + zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim'] + cuts = [ tc.notVetoed(), tc.converged(), tc.redChi2(5.), tc.z(zmin, zmax), tc.ellipseSq(250., 500.), tc.goodtracks(), + tc.muon1(1), tc.muon2(1), tc.doca(30.) ] + cuts.extend([tc.ip(float(i)) for i in range(100, 1000, 50)][::-1]) + cutnames = [ "not vetoed", "converged", "chi2/n<5", "z", "ellipse", "good tracks", "muon1", "muon2", "doca<30" ] + cutnames.extend(["IP<%s"%i for i in range(100, 1000, 50)][::-1]) entries = [] weights = [] print + r.gROOT.SetBatch(r.kTRUE) for i,cut in enumerate(cuts): - n = t.Draw("event", "&&".join([c for c in cuts[:i]])) - t.Draw("event","weight*("+"&&".join([c for c in cuts[:i]])+")") + n = t.Draw("event", "&&".join([c for c in cuts[:i+1]])) + t.Draw("event","weight*%s/%s("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")") w = t.GetHistogram().GetSum() / ntot entries.append(n) weights.append(w) - print cutname[i], '\t\t', n, '\t', w + print cutnames[i], '\t\t', n, '\t', w print @@ -516,13 +570,50 @@ {'name': 'ip', 'params':{'ipmax':250.}}, ] +cutConfigDict = { 'sig': [{'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_n_muons_in', 'params':{'n':1,'station':1}}, + {'name': 'has_n_muons_in', 'params':{'n':1,'station':2}}, + {'name': 'has_ecal', 'params':{}}, + {'name': 'doca', 'params':{'docamax':30.}}, + {'name': 'ip', 'params':{'ipmax':250.}}], + + 'bg': [{'name': 'not_vetoed', '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_n_muons_in', 'params':{'n':1,'station':1}}, + {'name': 'has_n_muons_in', 'params':{'n':1,'station':2}}, + {'name': 'has_ecal', 'params':{}}, + {'name': 'doca', 'params':{'docamax':30.}}, + {'name': 'ip', 'params':{'ipmax':250.}}] + +} + def numbers_v2(): mycuts, results = {}, {} for study in studies: print 'Reading %s...'%study - mycuts[study] = Cuts(cutConfigList, study).process(studies[study]['data']) + if study in ['tiny', 'pimu', 'mumunu']: + cutConfigList = cutConfigDict['sig'] + noB = '' + else: + cutConfigList = cutConfigDict['bg'] + noB = 'NoB_' + mycuts[study] = Cuts(cutConfigList, study).process(studies[study]['data'], noB) for cut in mycuts[study].cutList: + assert( len(cut.eventNumbers) == len(list(set(cut.eventNumbers))) ) 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:] ]