RecoTrack/test/analyseMVA.py

0001 #!/usr/bin/env python3
0002
0003 from builtins import range
0004 import array
0005 import collections
0006 import itertools
0007
0008 import ROOT
0009
0010 from Validation.RecoTrack.plotting.ntuple import *
0011 from Validation.RecoTrack.plotting.ntuplePlotting import *
0012
0013 # This file is a sketch of a simple analysis for MVA debugging
0014
0015 def selectMVA(track):
0016     mva = track.mva()
0017     return mva > 0.35 and mva < 0.6
0018
0019 def selectTrue(track):
0020     return track.nMatchedTrackingParticles() > 0
0021
0022 def selectFake(track):
0023     return track.nMatchedTrackingParticles() == 0
0024
0025 def main():
0026     ROOT.TH1.AddDirectory(False)
0027
0028     ntuple_new = TrackingNtuple("trackingNtuple.root")
0029     ntuple_old = TrackingNtuple("trackingNtuple_oldMVA.root")
0030
0031     common = dict(
0032         #ratio=True,
0033     )
0034     common_ylog = dict(
0035         ylog=True,
0036         ymin=0.5,
0037     )
0038     common_ylog.update(common)
0039
0040     opts = dict(
0041         mva = dict(xtitle="MVA", **common_ylog),
0042         pt  = dict(xtitle="p_{T}", xlog=True, **common_ylog),
0043         eta = dict(xtitle="#eta", **common),
0044         relpterr = dict(xtitle="p_{T} error / p_{T}", **common_ylog),
0045         absdxy = dict(xtitle="|d_{xy}(BS)|", **common_ylog),
0046         absdz = dict(xtitle="|d_{z}(BS)|", **common_ylog),
0047         absdxypv = dict(xtitle="|d_{xy}(closest PV)|", **common_ylog),
0048         absdzpv = dict(xtitle="|d_{z}(closest PV)|", **common_ylog),
0049         nhits = dict(xtitle="hits", **common),
0050         nlayers = dict(xtitle="layers", **common),
0051         nlayers3D = dict(xtitle="3D layers", **common),
0052         nlayersLost = dict(xtitle="lost layers", **common),
0053         minlost = dict(xtitle="min(inner, outer) lost layers", **common_ylog),
0054         lostmidfrac = dict(xtitle="(lost hits) / (lost + valid hits)", **common),
0055         ndof = dict(xtitle="ndof", **common),
0056         chi2 = dict(xtitle="chi2/ndof", **common_ylog),
0057         chi2_1Dmod = dict(xtitle="chi2/ndof with 1D modification", **common_ylog),
0058     )
0059
0060     if True:
0061         histos_new = histos(ntuple_new)
0062         histos_old = histos(ntuple_old)
0063         drawMany("newMVA_vs_oldMVA", [histos_old, histos_new], opts=opts)
0064
0065     if True:
0066         histos_new = histos(ntuple_new, selector=selectMVA)
0067         histos_old = histos(ntuple_old, selector=selectMVA)
0068         drawMany("newMVA_vs_oldMVA_mvaselected", [histos_old, histos_new], opts=opts)
0069
0070     if True:
0071         histos_new = histos(ntuple_new, selector=selectTrue)
0072         histos_old = histos(ntuple_old, selector=selectTrue)
0073         drawMany("newMVA_vs_oldMVA_true", [histos_old, histos_new], opts=opts)
0074
0075     if True:
0076         histos_new = histos(ntuple_new, selector=lambda t: selectTrue(t) and selectMVA(t))
0077         histos_old = histos(ntuple_old, selector=lambda t: selectTrue(t) and selectMVA(t))
0078         drawMany("newMVA_vs_oldMVA_true_mvaselected", [histos_old, histos_new], opts=opts)
0079
0080     if True:
0081         histos_new = histos(ntuple_new, selector=selectFake)
0082         histos_old = histos(ntuple_old, selector=selectFake)
0083         drawMany("newMVA_vs_oldMVA_fake", [histos_old, histos_new], opts=opts)
0084
0085     if True:
0086         histos_new = histos(ntuple_new, selector=lambda t: selectFake(t) and selectMVA(t))
0087         histos_old = histos(ntuple_old, selector=lambda t: selectFake(t) and selectMVA(t))
0088         drawMany("newMVA_vs_oldMVA_fake_mvaselected", [histos_old, histos_new], opts=opts)
0089
0090     if True:
0091         (histos_old, histos_new) = histos2(ntuple_old, ntuple_new, selectMVA)
0092         drawMany("newMVA_vs_oldMVA_mvaSelectedNew", [histos_old, histos_new], opts=opts)
0093
0094     if True:
0095         (histos_old, histos_new) = histos2(ntuple_old, ntuple_new, lambda t: selectTrue(t) and selectMVA(t))
0096         drawMany("newMVA_vs_oldMVA_true_mvaSelectedNew", [histos_old, histos_new], opts=opts)
0097
0098     if True:
0099         (histos_old, histos_new) = histos2(ntuple_old, ntuple_new, lambda t: selectFake(t) and selectMVA(t))
0100         drawMany("newMVA_vs_oldMVA_fake_mvaSelectedNew", [histos_old, histos_new], opts=opts)
0101
0102
0103 def makeHistos():
0104     h = collections.OrderedDict()
0105     def addTH(name, *args, **kwargs):
0106         _h = ROOT.TH1F(name, name, *args)
0107         if kwargs.get("xlog", False):
0108             axis = _h.GetXaxis()
0109             bins = axis.GetNbins()
0110             minLog10 = math.log10(axis.GetXmin())
0111             maxLog10 = math.log10(axis.GetXmax())
0112             width = (maxLog10-minLog10)/bins
0113             new_bins = array.array("d", [0]*(bins+1))
0114             new_bins[0] = 10**minLog10
0115             mult = 10**width
0116             for i in range(1, bins+1):
0117                 new_bins[i] = new_bins[i-1]*mult
0118             axis.Set(bins, new_bins)
0119         h[name] = _h
0120
0121     addTH("mva", 80, -1, 1)
0122     addTH("pt", 40, 0.1, 1000, xlog=True)
0123     addTH("eta", 60, -3, 3)
0124     addTH("relpterr", 20, 0, 1)
0125
0126     addTH("absdxy", 50, 0, 1)
0127     addTH("absdz", 30, 0, 15)
0128     addTH("absdxypv", 50, 0., 0.5)
0129     addTH("absdzpv", 20, 0, 1)
0130
0131     addTH("nhits", 41, -0.5, 40.5)
0132     addTH("nlayers", 26, -0.5, 25.5)
0133     addTH("nlayers3D", 26, -0.5, 25.5)
0134     addTH("nlayersLost", 6, -0.5, 5.5)
0135     addTH("minlost", 6, -0.5, 5.5)
0136     addTH("lostmidfrac", 20, 0, 1)
0137
0138     addTH("ndof", 20, 0, 20)
0139     addTH("chi2", 40, 0, 20)
0140     addTH("chi2_1Dmod", 40, 0, 20)
0141
0142     return h
0143
0144 def fillHistos(h, track):
0145     h["mva"].Fill(track.mva())
0146     h["pt"].Fill(track.pt())
0147     h["eta"].Fill(track.eta())
0148     h["ndof"].Fill(track.ndof())
0149     h["nlayers"].Fill(track.nPixelLay()+track.nStripLay())
0150     h["nlayers3D"].Fill(track.n3DLay())
0151     h["nlayersLost"].Fill(track.nLostLay())
0152     h["chi2"].Fill(track.nChi2())
0153     h["chi2_1Dmod"].Fill(track.nChi2_1Dmod())
0154     h["relpterr"].Fill(track.ptErr()/track.pt())
0155     h["nhits"].Fill(track.nValid())
0156     h["minlost"].Fill(min(track.nInnerLost(), track.nOuterLost()))
0157     h["lostmidfrac"].Fill(track.nInvalid() / (track.nValid() + track.nInvalid()))
0158
0159     h["absdxy"].Fill(abs(track.dxy()))
0160     h["absdz"].Fill(abs(track.dz()))
0161     h["absdxypv"].Fill(abs(track.dxyClosestPV()))
0162     h["absdzpv"].Fill(abs(track.dzClosestPV()))
0163
0164 def histos(ntuple, selector=None):
0165     h = makeHistos()
0166
0167     for event in ntuple:
0168         for track in event.tracks():
0169             if selector is not None and not selector(track):
0170                 continue
0171             fillHistos(h, track)
0172
0173     return h
0174
0175 def histos2(ntuple1, ntuple2, selector2=None):
0176     # assume the two ntuples have the very same tracks except possibly
0177     # for their parameters
0178     h1 = makeHistos()
0179     h2 = makeHistos()
0180
0181     for (event1, event2) in itertools.izip(ntuple1, ntuple2):
0182         #print event1.eventIdStr(), event2.eventIdStr()
0183         if event1.eventId() != event2.eventId():
0184             raise Exception("Inconsistent events %s != %s" % (event1.eventIdStr(), event2.eventIdStr()))
0185
0186         for (track1, track2) in itertools.izip(event1.tracks(), event2.tracks()):
0187
0188             if selector2 is not None and not selector2(track2):
0189                 continue
0190
0191             fillHistos(h1, track1)
0192             fillHistos(h2, track2)
0193     return (h1, h2)
0194
0195
0196 if __name__ == "__main__":
0197     main()