Newer
Older
FairShipTools / distribsForHNLandBG_byEvent.py
import sys,os
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)

datafolder = '/afs/cern.ch/user/e/egrave/private/ShipData/UsedForTP'
fpimu = r.TFile(datafolder + '/pimu/ShipAna_newGen_caloFix.root', 'read')
pimu = fpimu.Get('t')
pimu_geo = tools.searchForNodes3_xyz_dict(datafolder + '/pimu/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(datafolder + '/mumunu/ShipAna_newGen_caloFix.root', 'read')#r.TFile('../DATA/MUMUNU/ShipAna.root','read')
mumunu = fmumunu.Get('t')
mumunu_geo = tools.searchForNodes3_xyz_dict(datafolder + '/mumunu/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(datafolder + '/eenu-100MeV/ShipAna_newGen_caloFix.root', 'read')
tiny = ftiny.Get('t')
tiny_geo = tools.searchForNodes3_xyz_dict(datafolder + '/eenu-100MeV/geofile_full.10.0.Pythia8-TGeant4-0.1GeV-011.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')

# Neutrinos from yandex
#fnu = r.TFile('../DATA/nuFromYandex/newGen_15340000_nuBg.root')
#fnu = r.TFile.Open('root://eoslhcb.cern.ch//eos/ship/data/neutrinoBackground/Yandex_for_TP/newGen_11560000_nuBg_total.root','read')
fnu = r.TFile.Open('root://eoslhcb.cern.ch//eos/ship/data/neutrinoBackground/20150702/outputJune2015-nu-16700000-addWeight_100mradcut.root','read')
nu = fnu.Get('t')
nu_geo = tools.searchForNodes3_xyz_dict(datafolder + '/nuBgGeo/geofile_full.10.0.Genie-TGeant4.root')#neutrino77.root')

fdp = r.TFile(datafolder + '/../HNLasDP/ShipAna_newGen.root')
dp = fdp.Get('t')
dp_geo = tools.searchForNodes3_xyz_dict(datafolder + '/../HNLasDP/geofile_full.10.0.Pythia8-TGeant4.root')

# For a correct normalisation, we subtract the events that have an HNL vertex produced outside of the 60m fiducial volume
r.gROOT.SetBatch(r.kTRUE)
ndp = dp.Draw('weight', 'BadTruthVtx==0 && weight>0')
r.gROOT.SetBatch(r.kFALSE)


studies = {
    'pimu': {'data':pimu, 'geo':pimu_geo, 'seen':[], 'ntot':20000-125},#136},
    'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'seen':[], 'ntot':20000-280},#289},
    'tiny': {'data': tiny, 'geo':tiny_geo, 'seen':[], 'ntot':20000-43},#66},
    'nu': {'data': nu, 'geo':nu_geo, 'seen':[], 'ntot':1.534e7},
    'fakedp': {'data': dp, 'geo':dp_geo, 'seen':[], 'ntot':ndp},
    #'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?

"""
flux_nu = 1.05e+11 * 2.e+20/ 5.e+13#1.61e+11 * 2.e+20/ 4.e+13
flux_anti = 7.38e+10 * 2.e+20 / 5.e+13#1.02e+11 * 2.e+20 / 4.e+13
#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 =175*10e+3#entries_anti = 200000-2
weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu)
weight_anti ="(weight * %s / %s)"%(flux_anti,entries_anti)
"""

# ATTENTION: NEW WEIGHTS SINCE 2 JUL 2015
# WEIGHT FORMULA HAS CHANGED
# BUT I DID NOT UNDERSTAND IT :)
# Neutrinos
flux_nu = (1.05e+11 * 2.e+20 / (5.e+13)) * 6.7e-39 * 2.
entries_nu = 16700000
flux_anti = (1.05e+11 * 2.e+20 / (5.e+13)) * 6.7e-39
entries_anti = 1750000


class MyParticle(object):
    """ Simulating an ntuple filled by particle from one filled by event :) """
    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.__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)




class Cuts(object):
    def __init__(self, cutConfigList, study, byEvent=True): 
        self.byEvent = byEvent
        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, 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, noB, eff)): break
        else:
            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_vtxSq_in_ellipse(self, event, study, params):
            return vtxSqInAcceptance(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 cut_has_ecal(self, event, study, params):
            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
    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 has_n_muons_in_station(event, n, station):
    if n>2:
        print "has_n_muons_in_station ERROR: we have maximum 2 tracks!!"
        sys.exit()
    if station>2:
        print "has_n_muons_in_station ERROR: we looked only at the first and second stations!!"
        sys.exit()
    if n==1 and station==1:
        return event.Has1Muon1
    if n==2 and station==1:
        return event.Has2Muon1
    if n==1 and station==2:
        return event.Has1Muon2
    if n==2 and station==2:
        return event.Has2Muon2

def has_ecalDeposit(event):
    return event.HasEcal


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 vtxSqInAcceptance(event, study):
    xSq = event.vtxxSqr
    ySq = event.vtxySqr
    #if 'cosmics' in study or 'bg' in study:
    #   print study, x, y
    return tools.inEllipse(r.TMath.Sqrt(xSq), r.TMath.Sqrt(ySq), 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


class cutsWithDraw(object):
    def __init__(self):
        self.noB = 'NoB_'#''
        self.opt = '' # '_eff' # 'Ethr'
        self.nu = False
    def setNoB(self, noB):
        # to ignore or take into account Katerina's modifications
        self.noB = noB
    def recoed(self):
        cut = " recoed==1 && Sum$(BadTruthVtx==0) "
        if self.nu:
            cut = " recoed==1 "
        string = "Event reconstructed"
        return cut, string
    def nRecoed(self, num):
        cut = " nRecoed==%s "%num
        if num == 1:
            strn = "1 HNL candidate"
        elif num>1:
            strn = "%s HNL candidates"%num
        string = strn
        return cut, string
    def doca(self, th, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(%sDOCA < %s) "%(noB, th)
        if th == int(round(th)):
            th = int(round(th))
        string = "DOCA $<$ %s~cm"%th
        return cut, string
    def ip(self, th, noB=None):
        if not noB:
            noB = self.noB
        cut =  " Sum$(%sIP0 < %s) "%(noB, th)
        string = "IP $<$ %.1f m"%(th/100.)
        return cut, string
    def z(self, zmin, zmax, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(%svtxz > %s && %svtxz < %s) "%(noB, zmin, noB, zmax)
        string = "$z_{straw\, veto} < z_{vtx} < z_{trk1}$"
        return cut, string
    def ellipse(self, a, b, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(( (%svtxx*%svtxx)/(%s*%s) + (%svtxy*%svtxy)/(%s*%s) ) < 1.) "%(noB,noB, a,a, noB,noB, b,b)
        string = "Vertex $\\in $ ellipse"
        return cut, string
    def fiducial(self, zmin, zmax, a, b, noB=None):
        cutz = self.z(zmin, zmax, noB)[0]
        cute = self.ellipseSq(a, b, noB)[0]
        string = "Vtx in fiducial vol."
        cut = cutz+"&&"+cute
        return cut, string
    def ellipseSq(self, a, b, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(( (%svtxxSqr)/(%s*%s) + (%svtxySqr)/(%s*%s) ) < 1.) "%(noB, a,a, noB, b,b)
        string = "Vertex $\\in $ ellipse"
        return cut, string
    def converged(self):
        cut = " Sum$(DaughtersFitConverged==1) "
        string = "Track fit converged"
        return cut, string
    def goodtracks(self):
        cut = " Sum$(DaughtersAlwaysIn==1) "
        string = "Tracks in fiducial vol."# $\\in$ ellipse"
        return cut, string
    def redChi2(self, th, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(%sMaxDaughtersRedChi2<%s) "%(noB, th)
        if th == int(round(th)):
            th = int(round(th))
        string = "$\\chi^2/\\text{N.d.f.} < %s $"%th
        return cut, string
    def ndf(self, th, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(%sMinDaughtersNdf>%s) "%(noB, th)
        string = "$\\text{N.d.f.} > %s$"%th
        return cut, string
    def muon1(self, num):
        cut = " Sum$(Has%sMuon1==1) "%num
        if num==1:
            strn = "1 muon"
        elif num==2:
            strn = "2 muons"
        string = strn+" in 1$^{\\text{st}}$ muon station"
        return cut, string
    def muon2(self, num):
        cut = " Sum$(Has%sMuon2==1) "%num
        if num==1:
            strn = "1 muon"
        elif num==2:
            strn = "2 muons"
        string = strn+" in 2$^{\\text{nd}}$ muon station"
        return cut, string
    def ecal(self):
        cut = " Sum$(HasEcal==1) "
        string = "150~MeV in Ecal"
        return cut, string
    def notVetoed(self, th, noB=None, opt=None):
        if not noB:
            noB = self.noB
        if not opt:
            opt = self.opt
        cut = " Sum$(%sIP0 < %s) &&"%(noB, th)+" Sum$(numberOfHitLSSegments == 0) && Sum$(strawVetoAny%s==0) && Sum$(upstreamVetoAny%s==0) && Sum$(RPCany%s==0) "%(opt, opt, opt)
        string = "Event not vetoed"
        return cut, string
    def daughtersP(self, th, noB=None):
        if not noB:
            noB = self.noB
        cut = " Sum$(%sDaughtersMinP>%s) "%(noB, th)
        if th == int(round(th)):
            th = int(round(th))
        string = "Daughters $P > %s$ GeV"%th
        return cut, string


def countWithDraw(study, answers=False):
    t = studies[study]['data']
    geo = studies[study]['geo']
    ntot = studies[study]['ntot']
    tc = cutsWithDraw()
    if study == 'nu':
        tc.setNoB('NoB_')
        tc.nu = True
    zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim']
    zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim']
    if (study == 'pimu') or (study == 'fakedp'):
        cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], 
                     tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0], tc.doca(30.)[0], tc.ip(250.)[0] ]
        cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], 
                     tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1], tc.doca(30.)[1], tc.ip(250.)[1] ]
        if answers:
            # answering to referees questions
            cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], 
                         tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0],
                         tc.daughtersP(1.5)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.doca(1.)[0], tc.ip(10.)[0],
                         tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0] ]
            cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], 
                         tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1],
                         tc.daughtersP(1.5)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.doca(1.)[1], tc.ip(10.)[1],
                         tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1] ]
        if answers == 2:
            # answering to referees questions (NO VETO)
            cuts =     [ tc.recoed()[0], 
                         tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0],
                         tc.daughtersP(1.5)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.doca(1.)[0], tc.ip(10.)[0],
                         tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0] ]
            cutnames = [ tc.recoed()[1], 
                         tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1],
                         tc.daughtersP(1.5)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.doca(1.)[1], tc.ip(10.)[1],
                         tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1] ]

    elif (study == 'nu'):
        cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], 
                     tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0], tc.doca(30.)[0], tc.ip(250.)[0] ]
        cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], 
                     tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1], tc.doca(30.)[1], tc.ip(250.)[1] ]
        if answers:
            # answering to referees questions
            cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], 
                         tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0],
                         tc.daughtersP(1.5)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.doca(1.)[0], tc.ip(10.)[0],
                         tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0] ]
            cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], 
                         tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1],
                         tc.daughtersP(1.5)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.doca(1.)[1], tc.ip(10.)[1],
                         tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1] ]
        if answers == 2:
            # answering to referees questions (NO VETO)
            cuts =     [ tc.recoed()[0], 
                         tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0],
                         tc.daughtersP(1.5)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.doca(1.)[0], tc.ip(10.)[0],
                         tc.ecal()[0], tc.muon1(1)[0], tc.muon2(1)[0] ]
            cutnames = [ tc.recoed()[1], 
                         tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1],
                         tc.daughtersP(1.5)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.doca(1.)[1], tc.ip(10.)[1],
                         tc.ecal()[1], tc.muon1(1)[1], tc.muon2(1)[1] ]
    elif (study == 'mumunu'):
        cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], 
                     tc.muon1(2)[0], tc.muon2(2)[0], tc.doca(30.)[0], tc.ip(250.)[0] ]
        cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], 
                     tc.muon1(2)[1], tc.muon2(2)[1], tc.doca(30.)[1], tc.ip(250.)[1] ]
    elif (study == 'tiny'):
        cuts =     [ tc.recoed()[0], tc.notVetoed(1000.)[0], tc.redChi2(5.)[0], tc.ndf(25)[0], tc.fiducial(zmin, zmax, 250., 500.)[0], tc.goodtracks()[0], 
                     tc.ecal()[0], tc.doca(30.)[0], tc.ip(250.)[0] ]
        cutnames = [ tc.recoed()[1], tc.notVetoed(1000.)[1], tc.redChi2(5.)[1], tc.ndf(25)[1], tc.fiducial(zmin, zmax, 250., 500.)[1], tc.goodtracks()[1], 
                     tc.ecal()[1], tc.doca(30.)[1], tc.ip(250.)[1] ]
    # if not cutnames:
    #    cutnames = [ "converged", "not vetoed", "chi2/n<5", "ndf>15", "z", "ellipse", "good tracks", "ecal", "muon1", "muon2", "doca<30", "IP<250" ]
    # cutStrings = [  "Track fit converged",
    #                 "Event not vetoed",
    #                 "$\\xi^2/ndf < 5 $",
    #                 "$ndf > 15$",
    #                 "$z_{straw veto} < z_{vtx} < z_{Trk1}$",
    #                 "Vertex $\\in $ ellipse",
    #                 "Tracks $\\in$ ellipse",
    #                 "150~MeV in Ecal",
    #                 "1$^{\\text{st}}$ muon station",
    #                 "2$^{\\text{nd}}$ muon station",
    #                 "DOCA $<$ 30~cm",
    #                 "IP $<$ 2.5~m",
    #                 ]
    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+1]]))
        #if study == 'nu':
        #    t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")")
        #elif study == 'nubar':
        #    t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_anti, entries_anti) + str("&&".join([c for c in cuts[:i+1]])) + ")")
        #else:
        #    t.Draw("1.>>hw","(weight/%s)*("%studies[study]['ntot'] + str("&&".join([c for c in cuts[:i+1]])) + ")")
        if study == 'nu':
            t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_nu, entries_nu) + str("&&".join([c for c in cuts[:i+1]])) + ")")
        elif study == 'nubar':
            t.Draw("1.>>hw","(weight*%s/%s)*("%(flux_anti, entries_anti) + str("&&".join([c for c in cuts[:i+1]])) + ")")
        else:
            t.Draw("1.>>hw","(weight/%s)*("%ntot + str("&&".join([c for c in cuts[:i+1]])) + ")")
        w = t.GetHistogram().GetSumOfWeights()#Integral()#Integral('width')#GetSum()
        #w = r.gDirectory.Get('hw').Integral()
        entries.append(n)
        weights.append(w)
        print cutnames[i], '\t\t\t', n, '\t', w, cuts[i]
    print
    r.gROOT.SetBatch(r.kFALSE)
    return entries, weights, cutnames

def tableWithDraw(answers=False):
    results = {}
    for study in ['nu', 'pimu', 'mumunu', 'tiny', 'fakedp']:
        results[study] = {}
        numbers, numbersWithWeights, cuts = countWithDraw(study, answers)
        results[study]['numbers'] = numbers
        results[study]['numbersWithWeights'] = numbersWithWeights
        results[study]['cuts'] = cuts
    return results


def plotSigBg(cSaver, name, sig_tree, sig_ntot, bg_tree, what, cuts, xaxis, legSize, log, sig_cuts=None, bg_cuts=None):
    if not sig_cuts:
        sig_cuts = cuts
    if not bg_cuts:
        bg_cuts = cuts.replace("&& Sum$(BadTruthVtx==0) ", "")
    # colors
    bgcolor = r.kRed
    sigcolor = r.kBlack
    sigfill = r.kGray
    bgstyle = 1#9
    r.gStyle.SetOptTitle(r.kFALSE)
    # draw normalized
    c1 = r.TCanvas(name, name)
    whatbg = what
    #if what in ["DOCA", "vtxz"]:
    #    whatbg = "NoB_"+what
    if "z" in name:
        sig_tree.Draw(what+'>>h_sig_%s(100,-15500.,4500.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm')
    elif "pimu-IP0-lin" in name:
        sig_tree.Draw(what+'>>h_sig_%s(100,0.,20.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm')
    elif "IP0-lin" in name:
        sig_tree.Draw(what+'>>h_sig_%s(100,0.,250.)'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm')
    else:
        sig_tree.Draw(what+'>>h_sig_%s'%name, "(weight/%s)*("%sig_ntot + sig_cuts + ")", 'norm')
    hsig = r.gDirectory.Get('h_sig_%s'%name)
    hsig.SetMarkerColor(sigcolor)
    hsig.SetLineColor(sigcolor)
    hsig.SetFillColor(sigfill)
    hsig.GetXaxis().SetTitle(xaxis)
    hsig.GetYaxis().SetTitle("a. u.")
    hsig.GetYaxis().SetDecimals()
    if (("mumunu-mass" in name) or ("mumunu-IP0" in name) or ("tiny-mass" in name) or ("tiny-IP0" in name)) and not ("lin" in name):
        hsig.GetYaxis().SetRangeUser(0.0001,5.)
    #if "z" in name:
    #    print name
    #    hsig.GetXaxis().SetLimits(-17000.,4500.)
    #    hsig.Draw()
    hbg = None
    if not "no-bg-mass" in name:
        bg_tree.Draw(whatbg+'>>h_bg_%s'%name, "(weight*%s/%s)*("%(flux_nu, entries_nu) + bg_cuts + ")", 'norm same')
        hbg = r.gDirectory.Get('h_bg_%s'%name)
        hbg.SetMarkerColor(bgcolor)
        hbg.SetLineColor(bgcolor)
        hbg.SetLineStyle(bgstyle)
    if "y" in log:
        c1.SetLogy()
    if "x" in log:
        c1.SetLogx()
    newLegSize = list(legSize)
    if ("tiny" in name) and not ("2" in name):
        newLegSize[2] += 0.02
        newLegSize[3] += 0.02
    leg = r.TLegend(*newLegSize)
    leg.SetTextFont(42)
    leg.SetLineColor(r.kWhite)
    leg.SetBorderSize(0)
    sigstring = "HNL #rightarrow #pi#mu"
    if ('eenu' in name) or ('tiny' in name):
        sigstring = sigstring.replace("#pi#mu", "ee#nu")
    if ('mumu' in name):
        sigstring = sigstring.replace("#pi#mu", "#mu#mu#nu")
    leg.AddEntry(hsig, sigstring, 'lf' )
    if not "no-bg-mass" in name:
        leg.AddEntry(hbg, "Neutrino BG", 'l' )
    leg.Draw("same")
    c1.Modified(); c1.Update()
    print "\n", name, hsig.GetMean(), hsig.GetRMS()
    cSaver[name] = (c1, leg, hsig, hbg)
    return c1, leg, hsig, hbg

cSaver = []
def accPlotsTP(study):
    if not os.path.exists('accPlotsTP'):
        os.system('mkdir accPlotsTP')
    r.gROOT.ProcessLine(".x newGen/mystyle.C")
    t = studies[study]['data']
    geo = studies[study]['geo']
    ntot = studies[study]['ntot']
    tc = cutsWithDraw()
    zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim']
    zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim']
    # BG stuff
    bg_t = studies['nu']['data']
    bg_geo = studies['nu']['geo']
    bg_ntot = studies['nu']['ntot']
    bg_tc = cutsWithDraw()
    bg_tc.setNoB('NoB_')
    bg_zmin = geo['Veto_5']['z']['pos']+geo['Veto_5']['z']['dim']
    bg_zmax = geo['Tr1_1']['z']['pos']-geo['Tr1_1']['z']['dim']
    plot_ip = plotSigBg(study+"-IP0", t, ntot, bg_t, "IP0", tc.notVetoed(1000.)[0], "Impact parameter to target [cm]")[0]
    cSaver.append(plot_ip)
    raw_input('')




# per particle
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}},
]

# new analysis sw - per event
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_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.}},
]

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

header = {  'pimu': ("$HNL \\to \\pi \\mu$", "tab:recpimu"),
            'mumunu': ("$HNL \\to \\mu \\mu \\nu$", "tab:recmmnu"),
            'tiny': ("$HNL \\to e e \\nu$", "tab:receenu"),
            'nu': ("neutrino-induced background", "tab:recnubkg"),
            'fakedp': ("$\\gamma' \\to \\mu\\mu$", "tab:darkphoton"),
         }

def makeTex(results):
    for study in results.keys():
        weighted_string = "Acceptance"
        if study == 'nu':
            weighted_string = "Events / 5 years"
        results[study]['efficiencies'] = [ results[study]['numbersWithWeights'][i]/results[study]['numbersWithWeights'][i-1] if results[study]['numbersWithWeights'][i-1] != 0 else 0. for i in xrange(1, len(results[study]['numbersWithWeights']))]
        results[study]['effStrings'] = ['-'] + [ "{:.1f} \\%".format(float(e) * 100) for e in results[study]['efficiencies'] ]
        #for s in results[study]['effStrings']:
        #    print s
        print '\n\n\n'
        print '\t\\begin{table}[!!!!h]'
        print '\t\\begin{center}'
        print '\t\t\\caption{ Effect of the offline selection on %s \\label{%s}}'%(header[study][0],header[study][1])
        print '\t\t\\begin{tabular}{|l|c|c|c|}'
        print '\t\t\t\\hline'
        print '\t\t\t\\textbf{Selection} & \\textbf{Entries} & \\textbf{%s} & \\textbf{Selection efficiency} \\\\'%weighted_string
        print '\t\t\t\\hline'
        assert(len(results[study]['cuts']) == len(results[study]['numbers']) == len(results[study]['numbersWithWeights']))
        for i in xrange(len(results[study]['cuts'])):
            print '\t\t\t', results[study]['cuts'][i], ' & ', '%s'%int(results[study]['numbers'][i]), ' & ', '%s'%latex_float(results[study]['numbersWithWeights'][i]), ' & ', results[study]['effStrings'][i], '\\\\'
            print '\t\t\t\\hline'
        print '\t\t\\end{tabular}'
        print '\t\\end{center}'
        print '\t\\end{table}'
        print '\n\n'


def makeTexOld(results, cuts=None):
    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'
    if not cuts:
        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]