test/fakeAnalysis/graphics.py

0001 #!/usr/bin/env python3
0002
0003 from __future__ import print_function
0004 from builtins import range
0005 import ROOT
0006 from array import array
0007 from copy import copy
0008 from collections import OrderedDict
0009
0010 from Validation.RecoTrack.plotting.ntuple import *
0011 import analysis
0012
0013 from math import sqrt, copysign, sin, cos, pi
0014
0015 class EventPlotter(object):
0016     '''
0017     This class plots histograms and graphical objects with ROOT.
0018     The 3D and 2D graphics objects (except histograms) are stored in member lists plots_2D and plots_3D
0019     and all of them can be drawn with Draw method.
0020     '''
0021     def __init__(self):
0022         self.plots_2D = []
0023     self.plots_3D = []
0024     c = ROOT.TColor()
0025         self.colors_G = [c.GetColor(0,255,0), c.GetColor(0,185,0), c.GetColor(50,255,50), \
0026                    c.GetColor(0,100,0), c.GetColor(50,155,0), c.GetColor(0,70,155), \
0027                c.GetColor(0,255,0), c.GetColor(0,255,0), c.GetColor(0,255,0)]
0028     self.colors_B = [c.GetColor(0,0,255), c.GetColor(0,0,155), c.GetColor(50,50,255), \
0029                    c.GetColor(0,0,80), c.GetColor(50,0,155), c.GetColor(0,70,155), \
0030                c.GetColor(0,0,255), c.GetColor(0,0,255), c.GetColor(0,0,255)]
0031
0032     def Reset(self):
0033         self.plots_2D = []
0034     self.plots_3D = []
0035
0036     ###### OLD UNNECESSARY FUNCTIONS ######
0037
0038     def PlotEvent3DHits(self, event, flag="PSG"):
0039     '''
0040     Plots the 3D hits of an event.
0041     flag is an string which defines which hit types to plot
0042     (p for pixel, s for strip, g for glued)
0043     '''
0044     if('p' in flag or 'P' in flag):
0045         pixel_hits = event.pixelHits()
0046         pix_coords = []
0047         for hit in pixel_hits:
0048         pix_coords.append(hit.z())
0049         pix_coords.append(hit.x())
0050         pix_coords.append(hit.y())
0051         pix_plot = ROOT.TPolyMarker3D(len(pix_coords)/3, array('f', pix_coords), 1)
0052         pix_plot.SetMarkerColor(4)
0053
0054     if('s' in flag or 'S' in flag):
0055         strip_hits = event.stripHits()
0056         str_coords = []
0057         for hit in strip_hits:
0058         str_coords.append(hit.z())
0059         str_coords.append(hit.x())
0060         str_coords.append(hit.y())
0061         str_plot = ROOT.TPolyMarker3D(len(str_coords)/3, array('f', str_coords), 1)
0062         str_plot.SetMarkerColor(2)
0063
0064     if('g' in flag or 'G' in flag):
0065         glued_hits = event.gluedHits()
0066         glu_coords = []
0067         for hit in glued_hits:
0068         glu_coords.append(hit.z())
0069         glu_coords.append(hit.x())
0070         glu_coords.append(hit.y())
0071         glu_plot = ROOT.TPolyMarker3D(len(glu_coords)/3, array('f', glu_coords), 1)
0072         glu_plot.SetMarkerColor(3)
0073
0074     if('p' in flag or 'P' in flag): self.plots_3D.append(pix_plot)
0075     if('s' in flag or 'S' in flag): self.plots_3D.append(str_plot)
0076     if('g' in flag or 'G' in flag): self.plots_3D.append(glu_plot)
0077
0078     def PlotXY(self, event, limits=[-1000,1000], flag="PSG"):
0079     '''
0080     Plots the hits of an event in an XY plane.
0081     flag is an string which defines which hit types to plot
0082     (p for pixel, s for strip, g for glued)
0083     '''
0084
0085     if('p' in flag or 'P' in flag):
0086         pixel_hits = event.pixelHits()
0087         pix_x = []
0088         pix_y = []
0089         for hit in pixel_hits:
0090         if(limits[0] < hit.z() < limits[1]):
0091             pix_x.append(hit.x())
0092             pix_y.append(hit.y())
0093         pix_plot = ROOT.TGraph(len(pix_x), array('f', pix_x), array('f', pix_y))
0094         pix_plot.SetMarkerColor(4)
0095
0096     if('s' in flag or 'S' in flag):
0097         strip_hits = event.stripHits()
0098         str_x = []
0099         str_y = []
0100         for hit in strip_hits:
0101         if(limits[0] < hit.z() < limits[1]):
0102             str_x.append(hit.x())
0103             str_y.append(hit.y())
0104         str_plot = ROOT.TGraph(len(str_x), array('f', str_x), array('f', str_y))
0105         str_plot.SetMarkerColor(2)
0106
0107     if('g' in flag or 'G' in flag):
0108         glued_hits = event.gluedHits()
0109         glu_x = []
0110         glu_y = []
0111         for hit in glued_hits:
0112         if(limits[0] < hit.z() < limits[1]):
0113             glu_x.append(hit.x())
0114             glu_y.append(hit.y())
0115         glu_plot = ROOT.TGraph(len(glu_x), array('f', glu_x), array('f', glu_y))
0116         glu_plot.SetMarkerColor(3)
0117
0118     plot = ROOT.TMultiGraph()
0119
0120     if('p' in flag or 'P' in flag): plot.Add(pix_plot,"P")
0121     if('s' in flag or 'S' in flag): plot.Add(str_plot,"P")
0122     if('g' in flag or 'G' in flag): plot.Add(glu_plot,"P")
0123     self.plots_2D.append(plot)
0124
0125     def PlotZY(self, event, limits=[-1000,1000], flag="PSG"):
0126     '''
0127     Plots the hits of an event in an ZR plane.
0128     flag is an string which defines which hit types to plot
0129     (p for pixel, s for strip, g for glued)
0130     '''
0131
0132     if('p' in flag or 'P' in flag):
0133         pixel_hits = event.pixelHits()
0134         pix_z = []
0135         pix_y = []
0136         for hit in pixel_hits:
0137         if(limits[0] < hit.z() < limits[1]):
0138             pix_z.append(hit.z())
0139             pix_y.append(hit.y())
0140         pix_plot = ROOT.TGraph(len(pix_z), array('f', pix_z), array('f', pix_y))
0141         pix_plot.SetMarkerColor(4)
0142
0143     if('s' in flag or 'S' in flag):
0144         strip_hits = event.stripHits()
0145         str_z = []
0146         str_y = []
0147         for hit in strip_hits:
0148         if(limits[0] < hit.z() < limits[1]):
0149             str_z.append(hit.z())
0150             str_y.append(hit.y())
0151         str_plot = ROOT.TGraph(len(str_z), array('f', str_z), array('f', str_y))
0152         str_plot.SetMarkerColor(2)
0153
0154     if('g' in flag or 'G' in flag):
0155         glued_hits = event.gluedHits()
0156         glu_z = []
0157         glu_y = []
0158         for hit in glued_hits:
0159         if(limits[0] < hit.z() < limits[1]):
0160             glu_z.append(hit.z())
0161             glu_y.append(hit.y())
0162         glu_plot = ROOT.TGraph(len(glu_z), array('f', glu_z), array('f', glu_y))
0163         glu_plot.SetMarkerColor(3)
0164
0165     plot = ROOT.TMultiGraph()
0166
0167     if('p' in flag or 'P' in flag): plot.Add(pix_plot,"P")
0168     if('s' in flag or 'S' in flag): plot.Add(str_plot,"P")
0169     if('g' in flag or 'G' in flag): plot.Add(glu_plot,"P")
0170     self.plots_2D.append(plot)
0171
0172     def PlotTracksXY(self, tracks):
0173     ''' Plots tracks like polyline graphs in the XY plane. tracks is an iterable for tracks. '''
0174     plot = ROOT.TMultiGraph()
0175
0176     for track in tracks:
0177         X = []; Y = [];
0178         for hit in track.hits():
0179         if(hit.isValidHit()):
0180             X.append(hit.x())
0181             Y.append(hit.y())
0182         plot.Add(ROOT.TGraph(len(X),array("f",X),array("f",Y)),"L")
0183
0184     self.plots_2D.append(plot)
0185
0186     def PlotTracksZY(self, tracks):
0187     ''' Plots tracks like polyline graphs in the ZY plane. tracks is an iterable for tracks. '''
0188     plot = ROOT.TMultiGraph()
0189
0190     for track in tracks:
0191         Y = []; Z = [];
0192         for hit in track.hits():
0193         if(hit.isValidHit()):
0194             Y.append(hit.y())
0195             Z.append(hit.z())
0196         plot.Add(ROOT.TGraph(len(Z),array("f",Z),array("f",Y)),"L")
0197
0198     self.plots_2D.append(plot)
0199
0200     def PlotTracks3D(self, tracks):
0201     for track in tracks:
0202         X = []; Y = []; Z = [];
0203         for hit in track.hits():
0204         if(hit.isValidHit()):
0205             X.append(hit.x())
0206             Y.append(hit.y())
0207             Z.append(hit.z())
0208         if(X):
0209         self.plots_3D.append(ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y)))
0210
0211     def PlotPixelTracks3D(self, tracks):
0212     for track in tracks:
0213         X = []; Y = []; Z = [];
0214         for hit in track.pixelHits():
0215         if(hit.isValidHit()):
0216             X.append(hit.x())
0217             Y.append(hit.y())
0218             Z.append(hit.z())
0219         if(X):
0220         self.plots_3D.append(ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y)))
0221
0222     def PlotPixelGluedTracks3D(self, tracks):
0223     for track in tracks:
0224         X = []; Y = []; Z = [];
0225         for hit in track.hits():
0226         if(hit.isValidHit() and hit.hitType != 1):
0227             X.append(hit.x())
0228             Y.append(hit.y())
0229             Z.append(hit.z())
0230         if(X):
0231         self.plots_3D.append(ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y)))
0232
0233     def PlotTrack3D(self, track, color = 1):
0234     '''Plots a single track as a polyline and prints track info'''
0235     # Not so hasardous experimental edit:
0236     #hits = sorted([hit for hit in track.hits()], key = lambda hit: hit.index())
0237     #print [hit.index() for hit in hits]
0238     X = []; Y = []; Z = [];
0239     for hit in track.hits(): #hits: #track.hits():
0240         if(hit.isValidHit()):
0241         X.append(hit.x())
0242         Y.append(hit.y())
0243         Z.append(hit.z())
0244     if(not X):
0245         print("Track has no valid points")
0246         return
0247     plot = ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y))
0248     plot.SetLineColor(color)
0249     self.plots_3D.append(plot)
0250
0251     ''' Uncomment to print track info
0252     print "Track parameters:"
0253     print "px  : " + str(track.px())
0254     print "py  : " + str(track.py())
0255     print "pz  : " + str(track.pz())
0256     print "pt  : " + str(track.pt())
0257     print "eta : " + str(track.eta())
0258     print "phi : " + str(track.phi())
0259     print "dxy : " + str(track.dxy())
0260     print "dz  : " + str(track.dz())
0261     print "q   : " + str(track.q())
0262     '''
0263
0264     ###### METHODS USED BY PLOTFAKES ######
0265
0266     def TrackHelix(self, track, color = 1, style = 0):
0267     '''
0268     Creates a THelix object which can be plotted with Draw() method.
0269
0270     NOTE: The helixes are better drawn with the vertical z-axis.
0271     Even then the helixes are not drawn exactly correct.
0272     '''
0273     if isinstance(track, TrackingParticle):
0274         phi = track.pca_phi()
0275         dxy = track.pca_dxy()
0276         dz = track.pca_dz()
0277     else:
0278         phi = track.phi()
0279         dxy = track.dxy()
0280         dz = track.dz()
0281     xyz = array("d", [-dxy*ROOT.TMath.Sin(phi), dxy*ROOT.TMath.Cos(phi), dz])# [dz, -dxy*ROOT.TMath.Sin(phi), dxy*ROOT.TMath.Cos(phi)])
0282     v = array("d", [track.py(), track.px(), track.pz()]) #[track.px(), track.py(), track.pz()]) #[track.pz(), track.px(), track.py()])
0283     w = 0.3*3.8*track.q()*0.01 #Angular frequency = 0.3*B*q*hattuvakio, close enough
0284     z_last = dz
0285     for hit in track.hits():
0286         if(hit.isValidHit()): z_last = hit.z()
0287
0288     helix = ROOT.THelix(xyz, v, w)#, array("d", [dz, z_last]))#, ROOT.kHelixX, array("d", [1,0,0]))
0289     helix.SetAxis(array("d", [-1,0,0]))
0290     helix.SetRange(z_last, dz, ROOT.kHelixX)
0291     helix.SetLineColor(color)
0292     if style == 1: helix.SetLineStyle(9)
0293     if style == 2: helix.SetLineStyle(7)
0294     return helix
0295
0296     def Plot3DHelixes(self, tracks, color = 1, style = 0):
0297     ''' Plots 3D helixes from a track iterable '''
0298     for track in tracks:
0299         if(track.hits()):
0300         self.plots_3D.append(self.TrackHelix(track, color, style))
0301
0302     def Plot3DHelix(self, track, color = 1, style = 0):
0303     ''' Plots a single track helix '''
0304     if(track.hits()):
0305         self.plots_3D.append(self.TrackHelix(track, color, style))
0306
0307     def Plot3DHits(self, track, color = 1, style = 0):
0308     '''
0309     Plots the 3D hits from a track.
0310     '''
0311     pix_coords = []
0312     for hit in track.pixelHits():
0313         if hit.isValidHit():
0314         pix_coords.append(hit.z())
0315         pix_coords.append(hit.x())
0316         pix_coords.append(hit.y())
0317         if pix_coords:
0318         pix_plot = ROOT.TPolyMarker3D(len(pix_coords)/3, array('f', pix_coords), 2)
0319         pix_plot.SetMarkerColor(color)
0320         if style == 1: pix_plot.SetMarkerStyle(5)
0321         if style == 2: pix_plot.SetMarkerStyle(4)
0322         self.plots_3D.append(pix_plot)
0323
0324     for hit in track.gluedHits():
0325         if hit.isValidHit():
0326         x = hit.x(); y = hit.y(); z = hit.z()
0327         if hit.isBarrel():
0328             X = [x, x]
0329             Y = [y, y]
0330             Z = [z - sqrt(hit.zz()), z + sqrt(hit.zz())]
0331         else:
0332             X = [x - copysign(sqrt(hit.xx()),x), x + copysign(sqrt(hit.xx()),x)]
0333             Y = [y - copysign(sqrt(hit.yy()),y), y + copysign(sqrt(hit.yy()),y)]
0334             Z = [hit.z(), hit.z()]
0335         glu_plot = ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y))
0336         #glu_plot.SetLineStyle(2)
0337         if style == 1: glu_plot.SetLineStyle(2)
0338         if style == 2: glu_plot.SetLineStyle(3)
0339         glu_plot.SetLineColor(color)
0340         self.plots_3D.append(glu_plot)
0341
0342     for hit in track.stripHits():
0343         if hit.isValidHit():
0344         x = hit.x(); y = hit.y(); z = hit.z()
0345         if hit.isBarrel():
0346             X = [x, x]
0347             Y = [y, y]
0348             Z = [z - 1.5*sqrt(hit.zz()), z + 1.5*sqrt(hit.zz())]
0349         else:
0350             X = [x - 1.5*copysign(sqrt(hit.xx()),x), x + 1.5*copysign(sqrt(hit.xx()),x)]
0351             Y = [y - 1.5*copysign(sqrt(hit.yy()),y), y + 1.5*copysign(sqrt(hit.yy()),y)]
0352             Z = [hit.z(), hit.z()]
0353         str_plot = ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y))
0354         if style == 1: str_plot.SetLineStyle(2)
0355         if style == 2: str_plot.SetLineStyle(3)
0356         str_plot.SetLineColor(color)
0357         self.plots_3D.append(str_plot)
0358
0359     def PlotVertex3D(self, vertex, color=1):
0360     ''' Plots a single vertex object '''
0361     plot = ROOT.TPolyMarker3D(1, array('f', [vertex.z(), vertex.x(), vertex.y()]),3)
0362     plot.SetMarkerColor(color)
0363         self.plots_3D.append(plot)
0364
0365     def PlotPoint3D(self, point, color=1):
0366         ''' Plots a single 3D point from a 3-element list '''
0367     plot = ROOT.TPolyMarker3D(1, array('f', [point[2], point[0], point[1]]),3)
0368     plot.SetMarkerColor(color)
0369         self.plots_3D.append(plot)
0370
0371     def PlotTrackingFail(self, match):
0372         '''
0373     Plots the end of tracking hits for a fake.
0374     This is drawn with thick red lines from last shared hit to following fake and particle hits,
0375     and the particle hit is indicated with a purple star marker.
0376     '''
0377     X = array('f', [match.fake_end_loc[0], match.last_loc[0], match.particle_end_loc[0]])
0378     Y = array('f', [match.fake_end_loc[1], match.last_loc[1], match.particle_end_loc[1]])
0379     Z = array('f', [match.fake_end_loc[2], match.last_loc[2], match.particle_end_loc[2]])
0380         plot = ROOT.TPolyLine3D(3, Z, X, Y)
0381     plot.SetLineWidth(3)
0382     plot.SetLineColor(2)
0383         self.plots_3D.append(plot)
0384     self.PlotPoint3D(match.last_loc, 2)
0385     self.PlotPoint3D(match.particle_end_loc, 6)
0386
0387     def PlotFakes(self, event, reconstructed = 1, fake_filter = None, end_filter = None, last_filter = None, particle_end_filter = None):
0388     '''
0389         This is the main function to plot fakes in 3D with related TrackingParticles and tracks.
0390     Fake track is drawn as red, TrackingParticles as green and reconstructed tracks as blue.
0391     The detector scheme is also drawn with black circles.
0392     The decay vertices (yellow) and reconstructed collision vertices (colour of the track) are also drawn.
0393     Alongside with the fake plotting, the track and particle informations are printed.
0394
0395     Set reconstructed to false for not drawing other reconstructed tracks related to a fake,
0396     use filters (lists of class indexes or detector layer strings) to filter out everything else from plotting.
0397     '''
0398     iterative = 1 # simple edits for not needing to refactorise code
0399
0400         fakes = analysis.FindFakes(event)
0401     if iterative:
0402             # Plot fakes one by one
0403         for fake in fakes:
0404         fake_info = analysis.FakeInfo(fake)
0405                 # Check for filter
0406         if fake_filter and fake_info.fake_class not in fake_filter:
0407             continue
0408             if last_filter or particle_end_filter:
0409             found_flag = False
0410             for match in fake_info.matches:
0411             if (not end_filter or match.end_class in end_filter) and (not last_filter or last_filter in match.last_str) and (not particle_end_filter or particle_end_filter in match.particle_end_str):
0412                 found_flag = True
0413                 if not found_flag:
0414             continue
0415
0416         self.Reset()
0417         self.Plot3DHelixes([fake],2)
0418         self.Plot3DHits(fake, 2)
0419         #self.PlotVertex3D(vertex,2)
0420         analysis.PrintTrackInfo(fake, None, 0, fake_info)
0421                 if fake_info.matches:
0422             for match in fake_info.matches:
0423             self.PlotTrackingFail(match)
0424             #self.PlotPoint3D(fake_info.end_loc, 2)
0425
0426                 # Find real particle tracks which include fake tracks hits
0427         icol = 0
0428         particle_inds = []
0429         particles = []
0430         reco_inds = []
0431         recos = []
0432         for hit in fake.hits():
0433             if hit.isValidHit() and hit.nSimHits() >= 0:
0434             for simHit in hit.simHits():
0435                 particle = simHit.trackingParticle()
0436                 if particle.index() not in particle_inds:
0437                 particle_inds.append(particle.index())
0438                 particles.append(particle)
0439                 if reconstructed and particle.nMatchedTracks() > 0:
0440                 for info in particle.matchedTrackInfos():
0441                     track = info.track()
0442                     if track.index() not in reco_inds:
0443                     reco_inds.append(track.index())
0444                     recos.append(track)
0445
0446
0447             # Find reconstructed tracks included in fakes hits
0448             if hit.isValidHit() and reconstructed and hit.ntracks() > 0:
0449             for track in hit.tracks():
0450                 #track = info.track()
0451                 if (track.index() not in reco_inds) and track.index() != fake.index():
0452                 reco_inds.append(track.index())
0453                 recos.append(track)
0454
0455                 # Plot the particles and reconstructed tracks
0456         icol =  0
0457         self.ParticleTest(particles, draw=False)
0458                 for particle in particles:
0459             self.Plot3DHelix(particle, self.colors_G[icol], 1)
0460             self.plots_3D[-1].SetLineStyle(5)
0461             self.Plot3DHits(particle, self.colors_G[icol], 1)
0462             self.PlotVertex3D(particle.parentVertex(),self.colors_G[icol])
0463                     # EXPERIMENTAL LINE:
0464                     #self.PlotTrack3D(particle, self.colors_G[icol])
0465
0466             for decay in particle.decayVertices():
0467             self.PlotVertex3D(decay, 5)
0468             analysis.PrintTrackInfo(particle, fake)
0469             icol += 1
0470         icol = 0
0471         for track in recos:
0472             self.Plot3DHelix(track,self.colors_B[icol],2)
0473             self.Plot3DHits(track,self.colors_B[icol],2)
0474             #self.PlotVertex3D(vertex,self.colors_B[icol])
0475             analysis.PrintTrackInfo(track, fake)
0476             icol += 1
0477
0478         if hit.isValidHit() and hit.z() >= 0: self.PlotDetectorRange("p",4)
0479         elif hit.isValidHit() and hit.z() < 0: self.PlotDetectorRange("n",4)
0480         else: self.PlotDetectorRange("b",4)
0481
0482                 print("****************************\n")
0483             self.Draw()
0484             return
0485
0486     def DrawTrackTest(self, track):
0487     self.PlotDetectorRange("b", 4)
0488     self.PlotTrack3D(track)
0489
0490     def ParticleTest(self, particles, draw=False):
0491     for particle in particles:
0492         print("\nPARTICLE " + str(particle.index()))
0493         for hit in particle.hits():
0494         tof = -1
0495         for simHit in hit.simHits():
0496             if simHit.trackingParticle().index() == particle.index():
0497             #if len(info.tof()):
0498             tof = simHit.tof()
0499         print("Index: " + str(hit.index()) + ", Layer: " + str(hit.layerStr()) + ", TOF: " + str(tof) +\
0500             "     XY distance: " + str(sqrt(hit.x()**2 + hit.y()**2)) + ", Z: " + str(hit.z()))
0501         self.DrawTrackTest(particle)
0502         if draw:
0503         self.Draw()
0504         self.Reset()
0505
0506     def PlotEndOfTrackingPoints3D(self, last, fail, fake_classes = [], color = ROOT.kBlue):
0507     alpha = 0.01
0508         for i in range(len(last)):
0509             if fail and last[i] != fail[i]:
0510         X = array('f', [last[i][0], fail[i][0]])
0511         Y = array('f', [last[i][1], fail[i][1]])
0512         Z = array('f', [last[i][2], fail[i][2]])
0513         line = ROOT.TPolyLine3D(2, Z, X, Y)
0514         line.SetLineWidth(1)
0515         line.SetLineColorAlpha(color, alpha)
0516         self.plots_3D.append(line)
0517         else:
0518         point = ROOT.TPolyMarker3D(1, array('f', [last[i][2], last[i][0], last[i][1]]))
0519         point.SetMarkerColorAlpha(color, alpha)
0520         self.plots_3D.append(point)
0521
0522     def PlotEndOfTracking3D(self, last, fake_ends, particle_ends, end_classes = [], fake_classes = [], end_mask = [], fake_mask = []):
0523     '''
0524     Plots the end of tracking hits in 3D
0525     '''
0526     alpha = 0.01
0527         for i in range(len(last)):
0528         if (not end_mask or end_classes[i] in end_mask) and (not fake_mask or fake_classes[i] in fake_mask):
0529         point = ROOT.TPolyMarker3D(1, array('f', [last[i][2], last[i][0], last[i][1]]))
0530         point.SetMarkerColorAlpha(4, alpha)
0531         self.plots_3D.append(point)
0532
0533         point = ROOT.TPolyMarker3D(1, array('f', [fake_ends[i][2], fake_ends[i][0], fake_ends[i][1]]))
0534         point.SetMarkerColorAlpha(2, alpha)
0535         self.plots_3D.append(point)
0536
0537         point = ROOT.TPolyMarker3D(1, array('f', [particle_ends[i][2], particle_ends[i][0], particle_ends[i][1]]))
0538         point.SetMarkerColorAlpha(3, alpha)
0539         self.plots_3D.append(point)
0540
0541     def PlotEndOfTrackingPoints2D(self, last, fail, classes = [], color = ROOT.kBlue):
0542     alpha = 0.5
0543     plot = ROOT.TMultiGraph()
0544
0545         for i in range(len(last)):
0546             if fail and last[i] != fail[i]:
0547         #c = ROOT.TCanvas("s","S",1000,1000)
0548         R = array('f', [sqrt(last[i][0]**2 + last[i][1]**2), sqrt(fail[i][0]**2 + fail[i][1]**2)])
0549         Z = array('f', [last[i][2], fail[i][2]])
0550         #line = ROOT.TPolyLine(2, Z, R, "L")
0551         line = ROOT.TGraph(2, Z, R)
0552         line.SetLineWidth(1)
0553         line.SetLineColorAlpha(color, alpha)
0554         plot.Add(line, "L")
0555         #self.plots_2D.append(line)
0556         else:
0557         #point = ROOT.TMarker(array("f", [last[i][2]]), array("f", [sqrt(last[i][0]**2 + last[i][1]**2)]), 1)
0558         #point = ROOT.TMarker(last[i][2], sqrt(last[i][0]**2 + last[i][1]**2), 1)
0559         point = ROOT.TGraph(1, array("f", [last[i][2]]), array("f", [sqrt(last[i][0]**2 + last[i][1]**2)]))
0560         #point.SetDrawOption("*")
0561         point.SetMarkerColorAlpha(color, alpha)
0562         plot.Add(point, "P")
0563         #self.plots_2D.append(point)
0564     self.plots_2D.append(plot)
0565
0566
0567     def PlotDetectorRange(self, area="b", plates = 6):
0568     '''
0569     Plots the detector schematic layout.
0570     Area  means which part to plot (p = positive side, n = negative side, b = both)
0571     Plates means how many detector circle schemes to plot per side.
0572     '''
0573         r = [17, 17, 57, 57, 112, 112]
0574     z = [30, 70, 70, 120, 120, 280]
0575     for i in range(plates):
0576         X = []; Y = []; Z = []; ZN = []
0577         for t in range(0,101):
0578         X.append(r[i]*cos(2*pi*t/100))
0579         Y.append(r[i]*sin(2*pi*t/100))
0580         Z.append(z[i]*1.0)
0581         ZN.append(-1.0*z[i])
0582         plot = ROOT.TPolyLine3D(len(X),array("f",Z),array("f",X),array("f",Y),"S")
0583         nplot = ROOT.TPolyLine3D(len(X),array("f",ZN),array("f",X),array("f",Y),"S")
0584         plot.SetLineStyle(3)
0585         nplot.SetLineStyle(3)
0586         if area == "p": self.plots_3D.append(plot)
0587         if area == "n": self.plots_3D.append(nplot)
0588         if area == "b":
0589         self.plots_3D.append(plot)
0590         self.plots_3D.append(nplot)
0591
0592     ###### HISTOGRAM PLOTTING METHODS ######
0593
0594     def ClassHistogram(self, data, name = "Histogram", labels = False, ready = False, logy = False, xmin = None, xmax = None, nbin = None, small = False):
0595         '''
0596         A method for drawing histograms with several types of data.
0597     Data can be in two different forms:
0598         A dictionary (used for labeled data)
0599         Raw numerical data
0600
0601     In the case of the dictionary, labels (ordered dictionary) should contain class indexes with class names.
0602     Also if ready = True, the observations for classes are already calculated.
0603     '''
0604     if ready:
0605         if labels:
0606         keys = [key for key in data]
0607         hist_labels = [labels[key] for key in keys]
0608         hist_data = [data[key] for key in keys]
0609         else:
0610         hist_labels = [key for key in data]
0611         hist_data = [data[key] for key in hist_labels]
0612
0613     if small: c1 = ROOT.TCanvas(name, name, 500, 500)
0614     else: c1 = ROOT.TCanvas(name, name, 700, 500)
0615     c1.SetGrid()
0616     #c1.SetTopMargin(0.15)
0617     if ready: hist = ROOT.TH1F(name, name, 30,0,30)#, 3,0,3)
0618     else:
0619         if xmin != None and xmax != None and nbin != None:
0620         hist = ROOT.TH1F(name, name, nbin, xmin, xmax)
0621         elif "fraction" in name or isinstance(data[0], float):
0622         hist = ROOT.TH1F(name, name, 41,0,1.025)#100,0,5)#max(data)+0.0001)#11,0,1.1)#21,0,1.05)#41*(max(data)+1)/max(data))#, max(data)+1,-0.5,max(data)+0.5) #### EXPERIMENT!!!111
0623         elif isinstance(data[0], int): hist = ROOT.TH1F(name, name, max(data)-min(data)+1, min(data)-0.5,max(data)+0.5)
0624         else: hist = ROOT.TH1F(name, name, 3,0,3)
0625         #c1.SetLogy()
0626         #hist.SetBinsLength(1)
0627     if logy: c1.SetLogy()
0628
0629         hist.SetStats(0);
0630         hist.SetFillColor(38);
0631     hist.SetCanExtend(ROOT.TH1.kAllAxes);
0632     if ready: # cumulative results already calculated in data
0633         for i in range(len(hist_labels)):
0634         #hist.Fill(hist_labels[i],hist_data[i])
0635         [hist.Fill(hist_labels[i],1) for j in range(hist_data[i])]
0636     elif not ready and labels:
0637         for entry in data: hist.Fill(labels[entry],1)
0638         hist.LabelsOption(">")
0639     else: # data is in samples, not in cumulative results
0640         for entry in data: hist.Fill(entry,1)
0641         hist.LabelsOption("a")
0642     hist.LabelsDeflate()
0643     if name == "Classification of all fakes": pass
0644     #hist.SetLabelSize(0.05)
0645     #hist.GetXaxis().SetRange(0,10)
0646     hist.Draw()
0647
0648     input("Press any key to continue")
0649
0650     def EndOfTrackingHistogram(self, end_list, hist_type="end_class", end_mask = [], fake_mask = [], normalised = False, BPix3mask = False):
0651     ''' Plots several types of histograms related to end of tracking analysis '''
0652     data = []
0653         if hist_type == "end_class": # Classification of end of trackings
0654         for end in end_list:
0655         if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0656             data.append(end.end_class)
0657         self.ClassHistogram(data, "End of tracking classes", analysis.end_class_names, False)
0658     elif hist_type == "particle_pdgId":
0659         ids = OrderedDict()
0660         for end in end_list:
0661         if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0662             data.append(end.particle_pdgId)
0663             if str(end.particle_pdgId) not in ids:
0664             ids[end.particle_pdgId] = str(end.particle_pdgId)
0665         self.ClassHistogram(data, "Particle pdgId", ids, False)
0666     elif "invalid" not in hist_type:
0667         for end in end_list:
0668         if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0669             if hist_type == "last_layer": data.append(end.last_str.split()[0]) # split()[0] for taking only the name of the layer, skipping reason for invalid hit
0670                 elif hist_type == "fake_end_layer": data.append(end.fake_end_str.split()[0])
0671             elif hist_type == "particle_end_layer": data.append(end.particle_end_str.split()[0])
0672
0673             data_dict = copy(analysis.layer_data_tmp)
0674         for entry in data:
0675         #print entry
0676         data_dict[analysis.layer_names_rev[entry]] += 1
0677
0678         if normalised:
0679         norm_cff = analysis.Get_Normalisation_Coefficients()
0680         for i, v in data_dict.items():
0681             data_dict[i] = int(round(v*norm_cff[i]))
0682
0683         name = ""
0684         if hist_type == "last_layer": name = "End of tracking layers"
0685         elif hist_type == "fake_end_layer": name = "Fake end of tracking layers"
0686         elif hist_type == "particle_end_layer": name = "Particle end of tracking layers"
0687         if normalised: name += " normalised"
0688
0689         self.ClassHistogram(data_dict, name, analysis.layer_names , True)
0690
0691     else:
0692         if hist_type == "invalid_reason_fake":
0693         data_dict = copy(analysis.invalid_types_tmp)
0694         for end in end_list:
0695             if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0696                         if "missing" in end.fake_end_str: data_dict[1] += 1; #print end.fake_end_str
0697             elif "inactive" in end.fake_end_str: data_dict[2] += 1; #print end.fake_end_str
0698             elif "bad" in end.fake_end_str: data_dict[3] += 1; #print end.fake_end_str
0699             elif "missing_inner" in end.fake_end_str: data_dict[4] += 1; #print end.fake_end_str
0700             elif "missing_outer" in end.fake_end_str: data_dict[5] += 1; #print end.fake_end_str
0701             else: data_dict[0] += 1 # valid hit
0702
0703         self.ClassHistogram(data_dict, "Invalid hit reasons in fake EOT", analysis.invalid_types, True)
0704         elif hist_type == "invalid_reason_last":
0705         data_dict = copy(analysis.invalid_types_tmp)
0706         for end in end_list:
0707             if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0708                         if "missing" in end.last_str: data_dict[1] += 1
0709             elif "inactive" in end.last_str: data_dict[2] += 1
0710             elif "bad" in end.last_str: data_dict[3] += 1
0711             elif "missing_inner" in end.last_str: data_dict[4] += 1
0712             elif "missing_outer" in end.last_str: data_dict[5] += 1
0713             else: data_dict[0] += 1 # valid hit
0714
0715         self.ClassHistogram(data_dict, "Invalid hit reasons in EOT last shared hit", analysis.invalid_types, True)
0716         elif hist_type == "invalid_reason_particle":
0717         data_dict = copy(analysis.invalid_types_tmp)
0718         for end in end_list:
0719             if (not end_mask or end.end_class in end_mask) and (not fake_mask or end.fake_class in fake_mask) and (not BPix3mask or "BPix3" in end.last_str):
0720                         if "missing" in end.particle_end_str: data_dict[1] += 1
0721             elif "inactive" in end.particle_end_str: data_dict[2] += 1
0722             elif "bad" in end.particle_end_str: data_dict[3] += 1
0723             elif "missing_inner" in end.particle_end_str: data_dict[4] += 1
0724             elif "missing_outer" in end.particle_end_str: data_dict[5] += 1
0725             else: data_dict[0] += 1 # valid hit
0726
0727         self.ClassHistogram(data_dict, "Invalid hit reasons in particle EOT", analysis.invalid_types, True)
0728
0729     def ResolutionHistograms(self, res_data, tracks = "fake", fake_mask = [], draw = True, skip = True, normalised = False):
0730     '''
0731     Returns (and possibly draws) resolution histograms for the ResolutionData object.
0732         '''
0733     if tracks == "fake":
0734         if fake_mask:
0735         res_tmp = analysis.ResolutionData()
0736         for i in range(len(res_data)):
0737             if res_data.fake_class in fake_mask:
0738             res_tmp.pt.append(res_data.pt[i])
0739             res_tmp.eta.append(res_data.eta[i])
0740             res_tmp.Lambda.append(res_data.Lambda[i])
0741             res_tmp.cotTheta.append(res_data.cotTheta[i])
0742             res_tmp.phi.append(res_data.phi[i])
0743             res_tmp.dxy.append(res_data.dxy[i])
0744             res_tmp.dz.append(res_data.dz[i])
0745             res_tmp.fake_class.append(res_data.fake_class[i])
0746
0747         res_data = res_tmp
0748         color = ROOT.kRed
0749     elif tracks == "true":
0750         if normalised: color = ROOT.kBlue
0751         else: color = 38
0752     else:
0753         print("Unspecified data type")
0754         return
0755
0756     c1 = ROOT.TCanvas("Resolution_histograms","Resolution histograms", 1000, 900)
0757     c1.Divide(3,3)
0758     c1.SetGrid()
0759     #c1.SetLogy()
0760         histograms = []
0761
0762     c1.cd(1)
0763     histograms.append( self.Histogram1D([res.pt_delta/res.pt_rec for res in res_data], "Resolution of #Delta p_{t}/p_{t,rec} in " + tracks + " tracks", 100, -0.15, 0.15, "p_{t}/p_{t,rec}", color, draw, normalised) )
0764     c1.GetPad(1).SetLogy()
0765     c1.cd(2)
0766     if not skip: histograms.append( self.Histogram1D(res_data.eta, "Resolution of #eta in " + tracks + " tracks", 100, -0.02, 0.02, "#eta", color, draw, normalised) )
0767     c1.GetPad(2).SetLogy()
0768     c1.cd(3)
0769     if not skip: histograms.append( self.Histogram1D(res_data.Lambda, "Resolution of #lambda in " + tracks + " tracks", 100, -0.015, 0.015, "#lambda", color, draw, normalised) )
0770     c1.GetPad(3).SetLogy()
0771     c1.cd(4)
0772     histograms.append( self.Histogram1D(res_data.cotTheta, "Resolution of cot(#theta) in " + tracks + " tracks", 100, -0.03, 0.03, "cot(#theta)", color, draw, normalised) )
0773     c1.GetPad(4).SetLogy()
0774     c1.cd(5)
0775     histograms.append( self.Histogram1D(res_data.phi, "Resolution of #phi in " + tracks + " tracks", 100, -0.03, 0.03, "#phi", color, draw, normalised) )
0776     c1.GetPad(5).SetLogy()
0777     c1.cd(6)
0778     histograms.append( self.Histogram1D(res_data.dxy, "Resolution of dxy in " + tracks + " tracks", 100, -0.15, 0.15, "dxy", color, draw, normalised) )
0779     c1.GetPad(6).SetLogy()
0780     c1.cd(7)
0781     histograms.append( self.Histogram1D(res_data.dz, "Resolution of dz in " + tracks + " tracks", 100, -0.15, 0.15, "dz", color, draw, normalised) )
0782     c1.GetPad(7).SetLogy()
0783     #c1.SetLogy()
0784     if draw:
0785         c1.Draw()
0786         input("Press enter to continue")
0787
0788     return histograms
0789
0790     def ResolutionHistogramsCombine(self, res_data_trues, res_data_fakes, fake_mask = [], draw = True):
0791     '''
0792         Combines resolutionhistograms for true and matched fake tracks to a single set of histograms.
0793     '''
0794     true_hists = self.ResolutionHistograms(res_data_trues, "true", [], False, True)
0795     fake_hists = self.ResolutionHistograms(res_data_fakes, "fake", fake_mask, False, True)
0796
0797     c1 = ROOT.TCanvas("Resolution_histograms","Resolution histograms", 700, 900)
0798     c1.Divide(2,3)
0799     c1.SetGrid()
0800         histograms = []
0801
0802     for i in range(len(true_hists)):
0803         c1.cd(i+1)
0804         histograms.append( self.Histogram1DStacked(true_hists[i], fake_hists[i], true_hists[i].GetName() + " stacked") )
0805         c1.GetPad(i+1).SetLogy()
0806
0807
0808     c1.cd(i+2)
0809         leg = ROOT.TLegend(0.,0.2,0.99,0.8)
0810     leg.AddEntry(true_hists[0], "True tracks", "F")
0811     leg.AddEntry(fake_hists[0], "Fake tracks", "F")
0812     leg.Draw()
0813
0814
0815     if draw:
0816         c1.Draw()
0817         input("Press enter to continue")
0818
0819     return histograms
0820
0821     def ResolutionHistogramsNormalised(self, res_data_trues, res_data_fakes, fake_mask = [], draw = True):
0822     '''
0823         Plots the resolution histograms of true tracks and matched fake tracks, normalised to have an area of 1.
0824     '''
0825     true_hists = self.ResolutionHistograms(res_data_trues, "true", [], False, True, True)
0826     fake_hists = self.ResolutionHistograms(res_data_fakes, "fake", fake_mask, False, True, True)
0827
0828     c1 = ROOT.TCanvas("Resolution_histograms","Resolution histograms", 1100, 700) #700, 900)
0829     c1.Divide(3,2)
0830     c1.SetGrid()
0831         #histograms = []
0832
0833     for i in range(len(true_hists)):
0834         c1.cd(i+1)
0835         #true_hists[i].SetTitle("#splitline{True tracks: " + true_hists[i].GetTitle() + "}{Fake tracks: " + fake_hists[i].GetTitle() + "}")
0836         true_hists[i].SetTitle("Resolution of " + fake_hists[i].GetXaxis().GetTitle())
0837         true_hists[i].GetXaxis().SetTitle("#Delta" + true_hists[i].GetXaxis().GetTitle())
0838         #true_hists[i].SetTitleOffset(1.2)
0839         #c1.SetTopMargin(0.9)
0840         true_hists[i].DrawNormalized("")
0841         fake_hists[i].DrawNormalized("Same")
0842         #histograms.append( self.Histogram1DStacked(true_hists[i], fake_hists[i], true_hists[i].GetName() + " stacked") )
0843         #c1.GetPad(i+1).SetLogy()
0844
0845
0846     c1.cd(i+2)
0847         leg = ROOT.TLegend(0.1,0.3,0.9,0.7)
0848     leg.AddEntry(true_hists[0], "True tracks", "L")
0849     leg.AddEntry(fake_hists[0], "Fake tracks", "L")
0850     leg.Draw()
0851
0852
0853     if draw:
0854         c1.Draw()
0855         input("Press enter to continue")
0856
0857     #return histograms
0858
0859     def ResolutionHistogramsNormalisedCumulative(self, res_data_trues, res_data_fakes, fake_mask = [], draw = True):
0860     '''
0861         Plots the cumulative resolution histograms of true tracks and matched fake tracks from the normalised histograms from the previous method.
0862     '''
0863     true_hists = self.ResolutionHistograms(res_data_trues, "true", [], False, True, True)
0864     fake_hists = self.ResolutionHistograms(res_data_fakes, "fake", fake_mask, False, True, True)
0865
0866     c1 = ROOT.TCanvas("Resolution_histograms","Resolution histograms", 1100, 700) #700, 900)
0867     c1.Divide(3,2)
0868     c1.SetGrid()
0869         #histograms = []
0870
0871     for i in range(len(true_hists)):
0872         c1.cd(i+1)
0873         true_hists[i].SetTitle("Cumulative resolution of " + fake_hists[i].GetXaxis().GetTitle())
0874         true_hists[i].GetXaxis().SetTitle("#Delta" + true_hists[i].GetXaxis().GetTitle())
0875         #true_hists[i].SetTitleOffset(1.2)
0876         #c1.SetTopMargin(0.9)
0877         for res in res_data_trues:
0878         if i == 0 and res.pt_delta/res.pt_rec < -0.15: true_hists[i].Fill(-0.15)
0879         if i == 1 and res.cotTheta < -0.03: true_hists[i].Fill(-0.03)
0880         if i == 2 and res.phi < -0.03: true_hists[i].Fill(-0.03)
0881         if i == 3 and res.dxy < -0.15: true_hists[i].Fill(-0.15)
0882         if i == 4 and res.dz < -0.15: true_hists[i].Fill(-0.15)
0883         for res in res_data_fakes:
0884         if i == 0 and res.pt_delta/res.pt_rec < -0.15: fake_hists[i].Fill(-0.15)
0885         if i == 1 and res.cotTheta < -0.03: fake_hists[i].Fill(-0.03)
0886         if i == 2 and res.phi < -0.03: fake_hists[i].Fill(-0.03)
0887         if i == 3 and res.dxy < -0.15: fake_hists[i].Fill(-0.15)
0888         if i == 4 and res.dz < -0.15: fake_hists[i].Fill(-0.15)
0889         true_hists[i].Scale(1.0/len(res_data_trues.pt_rec))
0890         fake_hists[i].Scale(1.0/len(res_data_fakes.pt_rec))
0891         #c1.GetPad(i+1).Clear()
0892         true_hists[i].GetCumulative().Draw("")
0893         fake_hists[i].GetCumulative().Draw("Same")
0894
0895     c1.cd(i+2)
0896         leg = ROOT.TLegend(0.1,0.3,0.9,0.7)
0897     leg.AddEntry(true_hists[0], "True tracks", "L")
0898     leg.AddEntry(fake_hists[0], "Fake tracks", "L")
0899     leg.Draw()
0900
0901
0902     if draw:
0903         c1.Draw()
0904         input("Press enter to continue")
0905
0906     def Histogram1D(self, data, name = "Histogram", nbins = 100, xmin = -1, xmax = 1, xlabel = "", color = 38, draw = True, normalised = False):
0907     ''' Creates a single histogram of resolutions for one parameter list '''
0908
0909     hist = ROOT.TH1F(name, name, nbins, xmin, xmax)
0910
0911     hist.SetStats(0);
0912     if normalised:
0913         hist.SetLineColor(color)
0914         hist.SetLineWidth(1)
0915     else:
0916         hist.SetFillColor(color);
0917     hist.GetXaxis().SetTitle(xlabel)
0918     #hist.SetCanExtend(ROOT.TH1.kAllAxes);
0919     for entry in data: hist.Fill(entry,1)
0920
0921     f = ROOT.TF1("f1", "gaus", xmin, xmax)
0922     hist.Fit(f,"0")
0923
0924     hist.SetTitle("#sigma(#Delta" + xlabel + ") = " + str(f.GetParameter("Sigma"))[0:8])
0925
0926     if draw: hist.Draw()
0927     return hist
0928
0929     #input("Press any key to continue")
0930
0931     def HistogramSigmaPt(self, data, name = "Sigma(Pt) vs Pt"): # Old function
0932
0933     c1 = ROOT.TCanvas(name, name, 700, 500)
0934     c1.SetLogx()
0935     c1.SetLogy()
0936
0937     hist = ROOT.TH2F(name, name, 10, array("d",[10**(i/10.0) for i in range(-10, 12, 2)]), 100, -0.15, 0.15)
0938
0939     for res in data:
0940         hist.Fill(res.pt_sim, res.pt_delta/res.pt_rec)
0941
0942     hist.FitSlicesY()
0943     hist_2 = ROOT.gDirectory.Get(name + "_2")
0944
0945     hist_2.GetXaxis().SetTitle("p_{t,sim}")
0946     hist_2.GetYaxis().SetTitle("#sigma(#Delta p_{t}/p_{t,res})")
0947
0948     hist_2.Draw()
0949
0950     input("Press enter to continue")
0951
0952     def HistogramSigmaPtCombined(self, data1, data2, name = "Sigma(Pt) vs Pt"): # Old function
0953
0954     c1 = ROOT.TCanvas(name, name, 700, 500)
0955     c1.SetLogx()
0956     c1.SetLogy()
0957
0958     hist = ROOT.TH2F(name, name, 10, array("d",[10**(i/10.0) for i in range(-10, 12, 2)]), 100, -0.15, 0.15)
0959
0960     for res in data1:
0961         hist.Fill(res.pt_sim, res.pt_delta/res.pt_rec)
0962     for res in data2:
0963         hist.Fill(res.pt_sim, res.pt_delta/res.pt_rec)
0964
0965     hist.FitSlicesY()
0966     hist_2 = ROOT.gDirectory.Get(name + "_2")
0967
0968     hist_2.GetXaxis().SetTitle("p_{t,sim}")
0969     hist_2.GetYaxis().SetTitle("#sigma(#Delta p_{t}/p_{t,rec})")
0970
0971     hist_2.Draw()
0972
0973     input("Press enter to continue")
0974
0975     def HistogramSigmaPtDual(self, data1, data2, parameter = "pt", pad = None, name = "Sigma(parameter) vs Pt"):
0976     '''
0977     Plots the sigma parameter from Gaussian fits to a single parameter resolution with respect to simulated particle pt.
0978
0979     data1 = True tracks, data2 = Fake track
0980     '''
0981
0982
0983     if not pad:
0984         c1 = ROOT.TCanvas(name, name, 700, 500)
0985         c1.SetLogx()
0986         c1.SetLogy()
0987     else:
0988         pad.SetLogx()
0989         pad.SetLogy()
0990
0991     hist1 = ROOT.TH2F(name, name, 10, array("d",[10**(i/10.0) for i in range(-10, 12, 2)]), 100, -0.15, 0.15)
0992     hist2 = ROOT.TH2F(name + "_h2", name + "_h2", 10, array("d",[10**(i/10.0) for i in range(-10, 12, 2)]), 100, -0.15, 0.15)
0993     hist3 = ROOT.TH2F(name + "_h3", name + "_h3", 10, array("d",[10**(i/10.0) for i in range(-10, 12, 2)]), 100, -0.15, 0.15)
0994
0995     #print min(data.pt_sim)
0996     for res in data1:
0997         if parameter == "pt": hist1.Fill(res.pt_sim, res.pt_delta/res.pt_rec)
0998         if parameter == "cotTheta": hist1.Fill(res.pt_sim, res.cotTheta)
0999         if parameter == "phi": hist1.Fill(res.pt_sim, res.phi)
1000         if parameter == "dxy": hist1.Fill(res.pt_sim, res.dxy)
1001         if parameter == "dz": hist1.Fill(res.pt_sim, res.dz)
1002     for res in data2:
1003         if parameter == "pt": hist2.Fill(res.pt_sim, res.pt_delta/res.pt_rec)
1004         if parameter == "cotTheta": hist2.Fill(res.pt_sim, res.cotTheta)
1005         if parameter == "phi": hist2.Fill(res.pt_sim, res.phi)
1006         if parameter == "dxy": hist2.Fill(res.pt_sim, res.dxy)
1007         if parameter == "dz": hist2.Fill(res.pt_sim, res.dz)
1008
1009         hist3.Add(hist1)
1010     hist3.Add(hist2)
1011
1012     hist1.FitSlicesY()
1013     hist1_2 = ROOT.gDirectory.Get(name + "_2")
1014         hist1_2.SetLineColor(ROOT.kBlue)
1015     #hist1_2.GetYaxis().SetRange(0.01, 1)
1016
1017     hist2.FitSlicesY()
1018     hist2_2 = ROOT.gDirectory.Get(name + "_h2" + "_2")
1019     hist2_2.SetLineColor(ROOT.kRed)
1020
1021     hist2_2.GetXaxis().SetTitle("p_{t,sim}")
1022     hist2_2.SetTitleSize(0.04,"xy")
1023     if parameter == "pt":
1024         hist2_2.GetYaxis().SetTitle("#sigma(#Delta p_{t}/p_{t,rec})")
1025         hist2_2.SetMaximum(1.5)
1026         hist2_2.SetMinimum(0.003)
1027     if parameter == "cotTheta":
1028         hist2_2.GetYaxis().SetTitle("#sigma(#Delta cot(#theta))")
1029         hist2_2.SetMaximum(1.0)
1030         hist2_2.SetMinimum(0.0005)
1031     if parameter == "phi":
1032         hist2_2.GetYaxis().SetTitle("#sigma(#Delta #phi)")
1033         hist2_2.SetMaximum(1.0)
1034         hist2_2.SetMinimum(0.0005)
1035     if parameter == "dxy":
1036         hist2_2.GetYaxis().SetTitle("#sigma(#Delta dxy)")
1037         hist2_2.SetMaximum(1.3)
1038         hist2_2.SetMinimum(0.0005)
1039     if parameter == "dz":
1040         hist2_2.GetYaxis().SetTitle("#sigma(#Delta dz)")
1041         hist2_2.SetMaximum(1.5)
1042         hist2_2.SetMinimum(0.0005)
1043
1044
1045     hist2_2.SetStats(0)
1046     hist2_2.SetTitle("")
1047
1048     hist3.FitSlicesY()
1049     hist3_2 = ROOT.gDirectory.Get(name + "_h3" + "_2")
1050     hist3_2.SetLineColor(ROOT.kViolet)
1051
1052     hist2_2.Draw()
1053     hist3_2.Draw("Same")
1054     hist1_2.Draw("Same")
1055
1056     '''
1057     leg = ROOT.TLegend(0.25,0.75,0.45,0.9)
1058     leg.AddEntry(hist1_2, "True tracks", "L")
1059     leg.AddEntry(hist2_2, "Fake tracks", "L")
1060     leg.AddEntry(hist3_2, "Both", "L")
1061     leg.Draw()
1062
1063     input("Press enter to continue")
1064     '''
1065     return [hist1_2, hist2_2, hist3_2]
1066
1067     def ResolutionsSigmaVsPt(self, res_data_trues, res_data_fakes, fake_mask = [], draw = True):
1068         '''
1069         Plots the sigma parameter from Gaussian fits to all parameter resolutions with respect to simulated particle pt.
1070     '''
1071
1072     c1 = ROOT.TCanvas("Resolution_sigma_vs_pt","Resolution_sigma_vs_pt", 1100, 700) #700, 900)
1073     c1.Divide(3,2)
1074     c1.SetGrid()
1075         hists_tmp = []
1076
1077     pad = c1.cd(1)
1078     hists_tmp.append(self.HistogramSigmaPtDual(res_data_trues, res_data_fakes, "pt", pad, "pt"))
1079     pad = c1.cd(2)
1080     hists_tmp.append(self.HistogramSigmaPtDual(res_data_trues, res_data_fakes, "cotTheta", pad, "cotTheta"))
1081     pad = c1.cd(3)
1082     hists_tmp.append(self.HistogramSigmaPtDual(res_data_trues, res_data_fakes, "phi", pad, "phi"))
1083     pad = c1.cd(4)
1084     hists_tmp.append(self.HistogramSigmaPtDual(res_data_trues, res_data_fakes, "dxy", pad, "dxy"))
1085     pad = c1.cd(5)
1086     hists = self.HistogramSigmaPtDual(res_data_trues, res_data_fakes, "dz", pad, "dz")
1087     hists_tmp.append(hists)
1088
1089     c1.cd(6)
1090         leg = ROOT.TLegend(0.1,0.3,0.9,0.7)
1091     leg.AddEntry(hists[0], "True tracks", "L")
1092     leg.AddEntry(hists[1], "Fake tracks", "L")
1093     leg.AddEntry(hists[2], "Both", "L")
1094     leg.Draw()
1095
1096
1097     if draw:
1098         c1.Draw()
1099         input("Press enter to continue")
1100
1101
1102
1103     def Histogram1DStacked(self, hist1, hist2, name): # Old function
1104         hist_stack = ROOT.THStack(name, name)
1105     hist_stack.Add(hist2)
1106     hist_stack.Add(hist1)
1107
1108     hist = ROOT.TH1F(hist1)
1109     hist.Add(hist2)
1110     f = ROOT.TF1("f1", "gaus", hist.GetXaxis().GetXmin(), hist.GetXaxis().GetXmax())
1111     hist.Fit(f,"0")
1112     hist_stack.SetTitle("#sigma(" + hist.GetXaxis().GetTitle() + ") = " + str(f.GetParameter("Sigma")))
1113     hist_stack.Draw()#("nostack")
1114     hist_stack.GetXaxis().SetTitle(hist.GetXaxis().GetTitle())
1115
1116     #hist_stack.Draw("nostack")
1117     return hist_stack
1118
1119     def Histogram1D_Draw(self, data, name = "Histogram", nbins = 10): # Old function
1120     c1 = ROOT.TCanvas(name, name, 700, 500)
1121     #if "pt " in name: c1.SetLogx()
1122     c1.SetGrid()
1123
1124     hist = ROOT.TH1F(name, name, nbins, min(data),max(data))
1125     c1.SetLogy()
1126
1127     hist.SetStats(0);
1128         hist.SetFillColor(38);
1129     #hist.SetCanExtend(ROOT.TH1.kAllAxes);
1130     for entry in data: hist.Fill(entry,1)
1131     hist.Draw()
1132
1133     def Draw(self):
1134         ''' Draws everything else except histograms. '''
1135     if self.plots_3D:
1136         c3 = ROOT.TCanvas("3D Plot","3D Plot", 1200, 1200)
1137         axis = ROOT.TAxis3D()
1138         axis.SetXTitle("z")
1139         axis.SetYTitle("x")
1140         axis.SetZTitle("y")
1141         #axis.GetZaxis().RotateTitle()
1142
1143         for plot in self.plots_3D:
1144         if plot: plot.Draw()
1145         #axis.SetAxisRange(-278,278,"Z")
1146         axis.Draw()
1147         axis.ToggleZoom()
1148
1149     if self.plots_2D:
1150         c2 = ROOT.TCanvas("2D Plot","2D Plot", 1200, 1200)
1151         #axis = ROOT.TAxis()
1152
1153         for plot in self.plots_2D:
1154         plot.Draw("A")
1155         #axis.Draw()
1156
1157         '''
1158         plots = ROOT.TMultiGraph()
1159         for plot in self.plots_2D:
1160         plots.Add(plot) #plot.Draw("A")
1161             plots.Draw("A")
1162         '''
1163
1164         input("Press any key to continue")
1165
1166