Newer
Older
FairShipTools / BarbaraTP / int_plotTP.py
@Ubuntu Ubuntu on 28 May 2015 8 KB status just after the tp
from ROOT import *

gROOT.SetBatch(kTRUE)

'''f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root")
t = f.Get("t")

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
entries_nu = 1.345e+7#200000-2
entries_anti = 200000-2
weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu)
'''

what = "nu"


if what == "anti":
    f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_anti-nuBg.root")
    t = f.Get("t")
else:
    #f = TFile("/afs/cern.ch/work/a/austyuzh/public/4barbara/newGen_15340000_nuBg.root")#newGen/ShipAna_Barbara_nu_th45MeV_LSSegments45MeV_FINALTEST.root")
    f = TFile("../../DATA/nuFromYandex/newGen_15340000_nuBg.root")
    t = f.Get("t")



#b = TBrowser()

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
entries_nu = 1.534e+7#1.345e+7#200000-2
entries_anti =175*10e+4
if what=="anti":
    
    weight_nu = "(weight * %s / %s)"%(flux_anti,entries_anti)

else:
    weight_nu = "(weight * %s / %s)"%(flux_nu,entries_nu)


res = {'recoed':{}, 'not-vetoed':{}, 'selected':{}}
entries = t.GetEntries()
bbb = {}
for entry in xrange(entries):
    t.GetEntry(entry)
    intEl = t.nuIntNumSimpl
    if intEl in [12,8,11,0,9,16]:
        inter = 'nuDet'
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0
    elif intEl in [4,10,7,2]:
        inter = "vesselWalls"
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0
    elif intEl in [14,18,19,17]:
        inter = 'TrackSyst'
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0
    elif intEl in [15,]:
        inter = 'strawVeto'
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0
    elif intEl in [5,1]:
        inter = 'cave'
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0
    elif intEl in [6,3,13]:
        inter = 'vesselLids'
        #res['recoed'][inter]=0
        #res['not-vetoed'][inter]=0
        #res['selected'][inter]=0

    val = (1*t.weight*flux_nu/entries_nu)
        
    if inter in res['recoed']:
        res['recoed'][inter]+=val
    else:
        res['recoed'][inter]=val
        
    if (t.numberOfHitLSSegments == 0  and t.strawVetoAny_eff==0 and t.upstreamVetoAny_eff==0 and t.RPCany_eff==0 ):
        flag = False
        for p in t.NoB_IP0:
            if p < 1000.:
                flag = True
        if flag:
            if inter in res['not-vetoed']:
                res['not-vetoed'][inter]+=val
            else:
                res['not-vetoed'][inter]=val

    track_f = False
    cutElipse_f = False
    docaVtxIp_f = False
    doca_f = False
    ip_f = False
    vessel_f = False
    vtx_f = False
    for a in t.NoB_MaxDaughtersRedChi2:
        if a<5.:
            track_f=True
            if '1' in bbb:
                bbb['1']+=val
            else:
                bbb['1']=val
    assert(len(t.NoB_vtxxSqr)==len(t.NoB_vtxySqr))
    for (a,b) in zip(t.NoB_vtxxSqr,t.NoB_vtxySqr):
        if (a/(245.*245.) + b/(490.*490.) ) < 1.: 
            cutElipse_f = True
            #if t.event == 453 and t.jobID==2750601:
                #print t.event, t.jobID, t.event+t.jobID, a,b
                #print "all"
                #for (a,b) in zip(t.NoB_vtxxSqr,t.NoB_vtxySqr):
                #    print a, b
                #assert(False)
            if '2' in bbb:
                bbb['2']+=val
            else:
                bbb['2']=val

    for a in t.NoB_DOCA:
        if a < 30.: 
            doca_f=True
            if '3' in bbb:
                bbb['3']+=val
            else:
                bbb['3']=val
    for a in t.NoB_IP0:
        if a< 250.: 
            ip_f = True
            if '4' in bbb:
                bbb['4']+=val
            else:
                bbb['4']=val
    for a in t.NoB_vtxz:
        if a > -1958. and a<2568: 
            vessel_f=True
            if '5' in bbb:
                bbb['5']+=val
            else:
                bbb['5']=val
    if fabs(max(t.NoB_vtxz)-min(t.NoB_vtxz))<30:
        vtx_f = True
        if '6' in bbb:
            bbb['6']+=val
        else:
            bbb['6']=val
       
    if track_f and cutElipse_f and doca_f and ip_f and vessel_f and vtx_f:
        if inter in res['selected']:
            res['selected'][inter]+=val
        else:
            res['selected'][inter]=val

'''        
trackQ = "Sum$(NoB_MaxDaughtersRedChi2<5.)"
cutEllipse = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.))"
docaVtxIp = "Sum$(NoB_DOCA < 30.) && Sum$(NoB_IP0 < 250.) && Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)"
vtxz  = "fabs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30"
offline_cut = "%s && %s && %s && %s"%(docaVtxIp,vtxz,trackQ,cutEllipse)
'''

c1 = "Sum$(NoB_MaxDaughtersRedChi2<5.)"
c2 = "Sum$(( (NoB_vtxxSqr)/(245.*245.) + (NoB_vtxySqr)/(490.*490.) ) < 1.)"
c3 = "Sum$(NoB_DOCA < 30.)"
c4 = "Sum$(NoB_IP0 < 250.)" 
c5 = "Sum$(NoB_vtxz > -1958. && NoB_vtxz < 2568.)"
c6  = "TMath::Abs(Max$(NoB_vtxz)-Min$(NoB_vtxz))<30"
offline_cut = "%s && %s && %s && %s && %s && %s"%(c1,c2,c3,c4,c5,c6)

t.Draw("event","(%s)*%s"%(c1,weight_nu))
print bbb['1'],t.GetHistogram().GetSum()

t.Draw("event","(%s)*%s"%(c2,weight_nu))
print bbb['2'],t.GetHistogram().GetSum()

t.Draw("event","(%s)*%s"%(c3,weight_nu))
print bbb['3'],t.GetHistogram().GetSum() 

t.Draw("event","(%s)*%s"%(c4,weight_nu))
print bbb['4'],t.GetHistogram().GetSum()   

t.Draw("event","(%s)*%s"%(c5,weight_nu))
print bbb['5'],t.GetHistogram().GetSum() 

t.Draw("event","(%s)*%s"%(c6,weight_nu))
print bbb['6'],t.GetHistogram().GetSum() 


print res['selected']
    
a = 0
tot = {}
for tt in res:
    a = 0
    for k in res[tt]:
        a+=res[tt][k]
    tot[tt]=a

t.Draw("event","(%s)*%s"%(offline_cut,weight_nu))
print tot['selected'],t.GetHistogram().GetSum()  
    
for k in res['recoed']:
    if not k in res['not-vetoed']:
        res['not-vetoed'][k] = 0
    if not k in res['selected']:
        res['selected'][k] = 0

    
print "\\begin{table}"
print "   \\centering"
print "   \\begin{tabular}{lrrr}"
print "       \\hline"
print "       detector & reconstructed (\\%) & not vetoed (\\%) & selected (\\%) \\\\"
print "       \\hline"
print "       $\\nu$ detector & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['nuDet'],res['recoed']['nuDet']/tot['recoed']*100.,res['not-vetoed']['nuDet'],res['not-vetoed']['nuDet']/tot['not-vetoed']*100.,res['selected']['nuDet'],res['selected']['nuDet']/tot['selected']*100.)
print "       vessel lids & %.0f (%.1f)& %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['vesselLids'],res['recoed']['vesselLids']/tot['recoed']*100.,res['not-vetoed']['vesselLids'],res['not-vetoed']['vesselLids']/tot['not-vetoed']*100.,res['selected']['vesselLids'],res['selected']['vesselLids']/tot['selected']*100.)
print "       vessel walls & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']["vesselWalls"],res['recoed']["vesselWalls"]/tot['recoed']*100.,res['not-vetoed']["vesselWalls"],res['not-vetoed']["vesselWalls"]/tot['not-vetoed']*100.,res['selected']["vesselWalls"],res['selected']["vesselWalls"]/tot['selected']*100.)
#print "       straw veto & %.0f (%.0f) & %.0f (%.0f)& %.0f (%.0f)\\\\"%(res['recoed']['strawVeto'],res['recoed']['strawVeto']/tot['recoed'],res['not-vetoed']['strawVeto'],res['not-vetoed']['strawVeto']/tot['not-vetoed'],res['selected']['strawVeto'],res['selected']['strawVeto']/tot['selected'])
print "       tracking system & %.0f (%.1f) & %.1f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['TrackSyst'],res['recoed']['TrackSyst']/tot['recoed']*100.,res['not-vetoed']['TrackSyst'],res['not-vetoed']['TrackSyst']/tot['not-vetoed']*100.,res['selected']['TrackSyst'],res['selected']['vesselLids']/tot['selected']*100.)
print "       cave & %.0f (%.1f) & %.0f (%.1f)& %.0f (%.1f)\\\\"%(res['recoed']['cave'],res['recoed']['cave']/tot['recoed']*100.,res['not-vetoed']['cave'],res['not-vetoed']['cave']/tot['not-vetoed']*100.,res['selected']['cave'],res['selected']['cave']/tot['selected']*100.)
print "       \\hline"
print "       \\textbf{Total} & \\textbf{%.0f}  & \\textbf{%.0f} & \\textbf{%.0f} \\\\"%(tot['recoed'],tot['not-vetoed'],tot['selected'])
print "   \\end{tabular}"
print "\\end{table}"