Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-26 02:34:40

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